**************************  part-proc.f  *******************************
c this version of the program looks for any particles that are big
c enough to pass icutoff, but get rejected for either being a possible
c multi-particle (CM is not on a particle voxel) or else the 
c volumes don't agree within 3 or 5 %. 
c For these particles only, the particle has a different numbering
c system and ONLY a voxel VRML image is output.

c takes gray scale (8 bit) images from XCT and does automatic
c thresholding, 3 phase, with Otsu algorithm and also applies
c a limited watershed split, LWS=5. Automatically adjust if it is
c an interior scan (iscan=1) or and exterior scan (iscan=0).

c also identifies pores in individual particles, whether SH or non-SH,
c and saves them in files as individual particles, with numbering
c tied to the particle they come from (new as of April 30 2018). 
c spnum is number of pores found in a given particle
c savepnum(**) is an arry containing the voxel size of the pores found in
c a given particle, max label in savepnum is spnum
c as pores are found, their voxels are sequentially written into a dummyfile
c when all the pores are found in a given particle,
c then this dummyfile is rewound and the individual pores are all
c written out into separate files containing the header of the 
c particle number and the number of the pore in that particle
c pore files are similar to non-SH particle voxel files
c Everything else is the same as part-proc-3d-all.f

c  This program works on a 3-D file where all the particles have been labelled
c  such that particle 1 has label 1, particle 2 has label 2, etc.
c  Program identifies particles, and checks whether or not the particle is a
c  valid one or not (if center of mass is not inside particle, then particle
c  is thought to be a possible multi-particle and is not analyzed).
c  Once a valid particle is identified, then r(theta,phi) is interpolated
c  and the anm coefficients are computed and output.  Then a subroutine is called
c  to output a VRML file image of the particle.
c  Each file has the number of the particle, as well as an overall identifier
c  for the kind of particles examined (their origin).

c This version (7/16/08) segments, using numbers seen in ImagePro, and 
c also does the 3-D limited watershed split, that was originally written
c in the program LWS.f

c  (USER) Dimensions of main arrays
      real*8 aaa1,aaa0,aaa2,aaa3
c add 2 to each dimension of particle image
      integer*2 pix(1024,1024,1024),old(100000000,3),new(100000000,3)
c note that pixsmall is dimensioned 0.3 of pix in each direction. This
c allows for a particle that is 0.3 the size of the system, which is unlikely.
      integer*2 pixsmall(900,900,900)
      integer*2 burned,phase
      real p(0:200,0:500)
      real xg(1000),wg(1000),phi,pi
      real theta,lenx,leny,lenz
      complex y(0:200,-250:250),a(0:200,-250:250)
      complex yphi(0:200,-250:250)
      complex ytheta(0:200,-250:250)
      integer*8 ns
      integer*4 partvol(1000000,2),spnum,savepnum(1000000)
      integer*4 porevol(10000000)
      integer*4 partvolR(1000000,3)
      integer*2 save(100000000,4)
      integer*2 in(6),jn(6),kn(6),nijk
      integer*4 isave,inew,iold
      character*12 filestat
      character*36 fnonSH
      character*30 fSH
      character*16 micfile
      character*13 namef
      character*36 namef13
      character*24 namefff
      character*16 nameff

c  (USER) Unit = 9 is the input file, unit 7 is the output file
c  unit 8 is the particle data file (extent, etc.)
      open(unit=8,file='ASTM-AMPM-2-particles-data.dat')
      open(unit=88,file='ASTM-AMPM-2-part-volume.dat')
      open(unit=89,file='ASTM-AMPM-2-pore-volume.dat')
      open(unit=19,file='ASTM-AMPM-2-nonSH-internal-porosity.dat')
      open(unit=49,file='ASTM-AMPM-2-SH-internal-porosity.dat')
      open(unit=11,file='gauss120.dat')
      open(unit=27,file='ASTM-AMPM-2-sysconfig.dat')

      ng=120
      write(8,*) ' files opened'
      call flush(8)
c read in Gaussian points and weights
      do 925 i=1,ng
      read(11,*) xg(i),wg(i)
925   continue
      write(8,*) ' Gaussian file read'
      call flush(8)
      pi=4.0*atan(1.0)

c nxs, etc. are approximately 0.3 the size of the system
      nxs=900
      nys=900
      nzs=900
c save and old, new array dimensions
      nnsave=100000000
c note that save(m,4)= 0 is not a surface voxel, =1 is a surface voxel
      nnarray=100000000
      nsplit=500000000
c      nsplit=50
      iii=4

c  Label of burned pixel
      burned=-10000
      write(8,*) ' burned label = ',burned
      call flush(8)
c particles with volumes smaller than icutoff will not be processed
c If icutoff is chosen smaller than that value used with part-2d.f,
c then all particles stored in the particle image will be processed.
c      icutoff=512
      icutoff=5000
      jcutoff=6000000
c      icutoff=0
      
c USER
c determine number that filenames start with (numstar+1)
      numstar=0
c numbering pores, sequentially for all
      ipore=0
      numtot=numstar
c use to number non-SH particles
      ibad=0
      
c loop in 9 ASTM-AMPM-2 microstructures
      do 9900 isys=1,9
c      do 9900 isys=1,1
c (USER)sys size ns = nx ny nz, root name for files, and voxel scale
      read(27,*) filestat,nx,ny,nz,sx,iscan,nphase
      sy=sx
      sz=sx
      if(isys.eq.29) then
c start Apollo 14 systems, so start renumbering from zero
       numtot=0
       ibad=0
       ipore=0
      end if
c build particle structure filename
      micfile(1:12)=filestat
      micfile(13:16)='.mic'
      open(unit=9,file=micfile)
      write(8,*) ' microstructure file is  ',micfile
      write(8,*) ' actual image size ',nx,ny,nz
      write(8,*) ' sx sy sz = ',sx,sy,sz
      if(iscan.eq.1) write(8,*) ' interior scan'
      if(iscan.eq.0) write(8,*) ' exterior scan'
      call flush(8)
      ns=nx*ny*nz

c USER
c  Read in microstructure file (system of matrix and particles)
      write(8,*) ' before read'

      call flush(8)
c initialize pix
      do 340 k=1,1024
      do 340 j=1,1024
      do 340 i=1,1024
      if(iscan.eq.1) pix(i,j,k)=3
      if(iscan.eq.0) pix(i,j,k)=0
340   continue
      do 330 k=1,nz
      do 330 j=ny,1,-1
      do 330 i=1,nx
      read(9,*) pix(i,j,k)
330   continue
      close(9)
      write(8,*) ' after initialization and read'
      call flush(8)

      if(isys.eq.1) then
c output one original slice just for check
      k=500
      write(namef(5:9),'(I5)'),10000+k
      namef(1:5)='OriA-'
      namef(10:13)='.pgm'
      open(unit=14,file=namef)
      write(14,133) 
      write(14,134) 
      write(14,157) 
133   format('P2')
134   format('980 1010')
157   format('255')

      do 741 j=ny,1,-1
      do 741 i=1,nx
      write(14,22) pix(i,j,k)
22    format(i3)
741   continue
      close(14)

c convert pgm file to tiff

      namef13(1:8)='convert '
      write(namef13(13:17),'(I5)'),10000+k
      namef13(9:13)='OriA-'
      namef13(18:22)='.pgm '

      write(namef13(27:31),'(I5)'),10000+k
      namef13(23:27)='OriA-'
      namef13(32:36)='.tiff'

      call system(namef13)
c move tiff files to tiff/
      namefff(1:3)='mv '
      write(namefff(8:12),'(I5)'),10000+k
      namefff(4:8)='OriA-'
      namefff(13:24)='.tiff tiff/.'
      call system(namefff)

c delete pgm files
      nameff(1:3)='rm '
      write(nameff(8:12),'(I5)'),10000+k
      nameff(4:8)='OriA-'
      nameff(13:16)='.pgm'
      call system(nameff)
      end if

c segment the image using automatic method

c  Direction labels to check for burning path (nearest neighbor information
c  in 3-D digital system).
      in(1)=-1
      in(2)=1
      in(3)=0
      in(4)=0
      in(5)=0
      in(6)=0
      jn(1)=0
      jn(2)=0
      jn(3)=-1
      jn(4)=1
      jn(5)=0
      jn(6)=0
      kn(1)=0
      kn(2)=0
      kn(3)=0
      kn(4)=0
      kn(5)=-1
      kn(6)=1

      call pixbin(pix,nx,ny,nz,iscan,in,jn,kn,nphase)

      if(isys.eq.1) then
c output one slice after pixbin before LWS
      k=500
      write(namef(5:9),'(I5)'),10000+k
      namef(1:5)='PixA-'
      namef(10:13)='.pgm'
      open(unit=14,file=namef)
      write(14,133) 
      write(14,134) 
      write(14,137) 
137   format('1')

      do 942 j=ny,1,-1
      do 942 i=1,nx
      write(14,22) pix(i,j,k)
942   continue
      close(14)

c convert pgm file to tiff

      namef13(1:8)='convert '
      write(namef13(13:17),'(I5)'),10000+k
      namef13(9:13)='PixA-'
      namef13(18:22)='.pgm '

      write(namef13(27:31),'(I5)'),10000+k
      namef13(23:27)='PixA-'
      namef13(32:36)='.tiff'

      call system(namef13)
c move tiff files to tiff/
      namefff(1:3)='mv '
      write(namefff(8:12),'(I5)'),10000+k
      namefff(4:8)='PixA-'
      namefff(13:24)='.tiff tiff/.'
      call system(namefff)

c delete pgm files
      nameff(1:3)='rm '
      write(nameff(8:12),'(I5)'),10000+k
      nameff(4:8)='PixA-'
      nameff(13:16)='.pgm'
      call system(nameff)
      end if

c now do limited watershed split. Comment out subroutine 
c call if not desired. The variable lwss is the 
c number of LWS steps. 

      lwss=3
      call LWS(pix,nx,ny,nz,lwss,nnarray,iii,old,new)

      if(isys.eq.1) then
c output one slice after pixbin before LWS
      k=500
      write(namef(5:9),'(I5)'),10000+k
      namef(1:5)='LWSA-'
      namef(10:13)='.pgm'
      open(unit=14,file=namef)
      write(14,133) 
      write(14,134) 
      write(14,137) 

      do 943 j=ny,1,-1
      do 943 i=1,nx
      write(14,22) pix(i,j,k)
943   continue
      close(14)

c convert pgm file to tiff

      namef13(1:8)='convert '
      write(namef13(13:17),'(I5)'),10000+k
      namef13(9:13)='LWSA-'
      namef13(18:22)='.pgm '

      write(namef13(27:31),'(I5)'),10000+k
      namef13(23:27)='LWSA-'
      namef13(32:36)='.tiff'

      call system(namef13)
c move tiff files to tiff/
      namefff(1:3)='mv '
      write(namefff(8:12),'(I5)'),10000+k
      namefff(4:8)='LWSA-'
      namefff(13:24)='.tiff tiff/.'
      call system(namefff)

c delete pgm files
      nameff(1:3)='rm '
      write(nameff(8:12),'(I5)'),10000+k
      nameff(4:8)='LWSA-'
      nameff(13:16)='.pgm'
      call system(nameff)
      end if

c This next part of the program does the actual burning = particle
c identification. Potentially good particles have their small intermal
c pores cleaned up, then are mathematically analyzed.  Particle whose
c interpolated volumes are near enough to their digital volumes then have
c their anm coefficients computed and stored, and then optionally a VRML
c file is generated that contains the spherical harmonic image and the
c raw voxel image.

      countd=0.0
c ipore is the counter for the pores found in particles
c start on 2,nx; 2,ny; 2,nz with ntouch = 0. That way, if i1 or j1 or k1 touch edge,
c ntouch will register and change to 1 and invalidate that particle.
      do 1000 k=2,nz-1
      do 1000 j=2,ny-1
      do 1000 i=2,nx-1

c don't choose any pore pixel or previously burned and identified particle
c pixel, which is still equal to the value of burned. or background pixel
      if(pix(i,j,k).le.0.or.pix(i,j,k).eq.3) goto 1000
      ipart=pix(i,j,k)
      write(8,*) 
      write(8,*) ' found particle '
      write(8,*) ' at i j k = ',i,j,k,' ipart = ',ipart
      call flush(8)

      iold=0
      isave=0
c  find a first particle pixel
      if(pix(i,j,k).eq.ipart) then
      pix(i,j,k)=burned
      iold=iold+1
      isave=isave+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
c save particle voxel coordinates
      save(iold,1)=i
      save(iold,2)=j
      save(iold,3)=k
      save(iold,4)=0
      do 995 iin=1,6
      ii2=i+in(iin)
      jj2=j+jn(iin)
      kk2=k+kn(iin)
      if(pix(ii2,jj2,kk2).eq.0) save(iold,4)=1
995   continue
      ntouch=0
      end if

c  Now start building up new burned pixels from old set of burned pixels,
c  thus propagating the fire and identifying new particle pixels connnected
c  to the ones already burned.

60    inew=0
      do 100 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
      kk=old(ijk,3)
c  check all six nearest neighbors of previously burned pixel
      do 90 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

c  Hard boundary conditions (can't pass beyond boundary)

c don't analyze particles that touch the edges - probably artificially
c cut, so shape is not right
      if(k1.eq.1.or.k1.eq.nz.or.i1.eq.1.or.i1.eq.nx.or.j1.eq.1.
     &      or.j1.eq.ny) then
        if(pix(i1,j1,k1).eq.ipart) ntouch=1
      end if
c if particle touches the outer circular border, then could be cut off by virtual core, delete
      if(pix(i1,j1,k1).eq.3) ntouch=1

      if(k1.lt.1.or.k1.gt.nz) then
      goto 90
      end if
      if(j1.lt.1.or.j1.gt.ny) then
      goto 90
      end if
      if(i1.lt.1.or.i1.gt.nx) then
      goto 90
      end if

c  Store (i,j,k) labels of newly burned pixels in array new().
      if(pix(i1,j1,k1).eq.ipart) then
      pix(i1,j1,k1)=burned
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      new(inew,3)=k1
      end if
90    continue
100   continue

c  If new pixels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 150 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
      old(ijk,3)=new(ijk,3)
      save(ijk+isave,1)=new(ijk,1)
      save(ijk+isave,2)=new(ijk,2)
      save(ijk+isave,3)=new(ijk,3)
      save(ijk+isave,4)=0
      do 895 iin=1,6
      ii2=new(ijk,1)+in(iin)
      jj2=new(ijk,2)+jn(iin)
      kk2=new(ijk,3)+kn(iin)
      if(pix(ii2,jj2,kk2).eq.0) save(ijk+isave,4)=1
895   continue

150   continue
      isave=isave+inew
      goto 60
      end if

c  CHECK 1 - is particle too small?

c  If particle is not big enough, go to restore, wipe it out,
c  and try again.  Do not store the particle. icutoff is the minimum
c  number of pixels that a particle must have to be stored.
c      if(isave.le.10) goto 444
      if(isave.lt.icutoff.or.isave.gt.jcutoff) then
       write(8,*) ' found volume = ',isave
       if(isave.lt.icutoff) write(8,*) ' rejected because 
     &     smaller than icutoff = ',icutoff
       if(isave.gt.jcutoff) write(8,*) ' rejected because 
     &     larger than jcutoff = ',jcutoff
      call flush(8)
       goto 444
      end if

c  CHECK 2 - particle was touching image or virtual core boundaries, so probably
c  wrong shape from being cut off. Do not store the particle.
      if(ntouch.eq.1) then
       write(8,*) ' rejected-touching boundary or virtual core boundary'
       write(8,*) ' set particle to zero = background'
      call flush(8)
       goto 444
      end if

c  CHECK 2 - is particle a probable multi-particle? If so, fill in
c internal pores before declaring to be non-SH.
c ***********************************************************************
c at this point, I want to save images of the non-SH particles. They
c could be multi-particles, or they could have a large, central
c pore, so that the CM is not in the particle.
c In this case, I want to run small first, to clean out that
c central pore, since the particle may be SH without that 
c central pore. Without doing this, I then am making a difference
c between central and non-central pores, which I should not.
c This will serve for ALL particles, so I don't have to run small again.
c************************************************************************

c Find voxel that contains the true center of mass. Call small() and
c fill in small internal pores. If this voxel is then
c inside the particle, then go ahead and calculate anm, see if
c the particle is an SH particle or not, depending on match of 
c SH volume to voxel volume. If, after small, the center of mass
c voxel is still outside the particle, then particle must be 
c a multi-particle and should be classified as non-SH.

      write(8,*) ' particle size = ',isave
      xcm=0.0
      ycm=0.0
      zcm=0.0
      iright=0
      jright=0
      kright=0
      ileft=100000
      jleft=100000
      kleft=100000
      do 500 ii=1,isave
      xcm=xcm+save(ii,1)
      ycm=ycm+save(ii,2)
      zcm=zcm+save(ii,3)
      if(save(ii,1).lt.ileft) ileft=save(ii,1)
      if(save(ii,2).lt.jleft) jleft=save(ii,2)
      if(save(ii,3).lt.kleft) kleft=save(ii,3)
      if(save(ii,1).gt.iright) iright=save(ii,1)
      if(save(ii,2).gt.jright) jright=save(ii,2)
      if(save(ii,3).gt.kright) kright=save(ii,3)
500   continue
c this is the true center of mass, in voxel units
      xcm=xcm/float(isave)
      ycm=ycm/float(isave)
      zcm=zcm/float(isave)
c  now identify which voxel contains this center of mass
      icm=xcm+0.5
      jcm=ycm+0.5
      kcm=zcm+0.5
      write(8,*) ' real cm = ',xcm,ycm,zcm
      write(8,*) ' voxel cm = ',icm,jcm,kcm

      if(ileft.eq.1.or.iright.eq.nx.or.jleft.eq.1.or.jright.eq.ny
     &   .or.kleft.eq.1.or.kright.eq.nz) then
c this was taken care of before, but this is a double check
      write(8,*) ' particle touches edge, artificial flattening'
      write(8,*) ' do not extract it, but just erase it'
      write(8,*) 'Actual center of particle = ',xcm,ycm,zcm
      write(8,*) 'Pixel center of particle = ',icm,jcm,kcm
      write(8,*) 'Value at center of particle = ',pix(icm,jcm,kcm)
      write(8,*) ' Min pixel labels (ijk) ',ileft,jleft,kleft
      write(8,*) ' Max pixel labels (ijk) ',iright,jright,kright
      call flush(8)
      goto 444
      end if

c check number of burned in defined box
      check=0.0
      do 1270 k3=kleft,kright
      do 1270 j3=jleft,jright
      do 1270 i3=ileft,iright
      if(pix(i3,j3,k3).eq.burned) check=check+1.0
1270  continue
      write(8,*) ' check burned before call small = ',check
      write(8,*) ' before call small '
      write(8,*) ' Min pixel labels (ijk) ',ileft,jleft,kleft
      write(8,*) ' Max pixel labels (ijk) ',iright,jright,kright
      write(8,*) ' pix(ileft,jleft,kleft)= ',pix(ileft,jleft,kleft)
      
c find any small internal pores and wipe them out,
c recording them in a file before doing any other analysis
c this file wil eventually be given a name that contains
c the actual particle name, whether SH or nonSH

      write(8,*) 'going into small volume = ',isave
      call flush(8)
      call small(pix,pixsmall,nx,ny,nz,nxs,nys,nzs,ileft,jleft,
     &            kleft,iright,jright,kright,burned,nnarray,
     &            old,new,iii,isave,save,nnsave,bcount,ipore,
     &            sx,sy,sz,porevol,filestat,savepnum,spnum)

      write(8,*) 'out of  small'
      call flush(8)
c note - all internal pores come out of small with value
c of burned, isave and save have been updated

c this is the case where, even after filling in 
c the small internal pores, the voxel containing the true
c CM is still not in the particle
c this must mean that particle is made up of joined particles
c that still have space between them in which the CM is located.
c this kind of particle is not SH and so must be saved voxel by voxel.

      if(pix(icm,jcm,kcm).ne.burned) then
      countd=countd+1.0
      write(8,*) ' voxel volume = ',isave
      write(8,*) ' CM NOT IN PARTICLE '
      write(8,*) ' pix(icm,jcm,kcm) = ',pix(icm,jcm,kcm)
      write(8,*) ' probable multi-particle'
      write(8,*) ' save as voxel collection'
c use ibad to number the file
      ibad=ibad+1
      iibadd=10
      partvolR(ibad,1)=ipart
      partvolR(ibad,2)=isave
      partvolR(ibad,3)=iibadd
      write(8,*) ' type of non-SH particle ',iibadd
      write(8,*) ' non-SH particle number ',ibad
      write(8,*) 'Actual center of particle = ',xcm,ycm,zcm
      write(8,*) 'Pixel center of particle = ',icm,jcm,kcm
      write(8,*) ' Min pixel labels (ijk) ',ileft,jleft,kleft
      write(8,*) ' Max pixel labels (ijk) ',iright,jright,kright
      write(8,*) ' going into voxpixR'
      call flush(8)
      call voxpixR(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,numtot,lenx,kright,sx,sy,sz,ibad,iibadd,filestat)
      write(8,*) ' coming out of voxpixR'
      call flush(8)
      call nonSHout(isave,save,ibad,nnsave,iii,filestat,sx)
c create pore filename using filestat and ibad and iibadd
      fnonSH(1:12)=filestat
      write(fnonSH(19:24),'(I6)') 100000+ibad
      write(fnonSH(30:32),'(I3)') 100+iibadd
      fnonSH(13:19)='-nonSH-'
      fnonSH(25:30)='-pore-'
      fnonSH(33:36)='.dat'
c record porosity, even if zero
       write(19,819) fnonSH,isave,sx,
     &     100.*bcount/float(isave)
819   format(a36,i10,f12.8,' %porosity=',f12.8)
       call flush(19)
c now output volumes of pores found in particle only if bcount > 0
      if(bcount.gt.0.0) then
      open(unit=77,file=fnonSH)
      write(77,771) spnum
      do 5322 i44=1,spnum
      write(77,771) i44,savepnum(i44)
771   format(i6,i10)
5322  continue
      close(77)
      end if
      goto 444
      end if
      
c  possible SH particle, go ahead, record and analyze
c last check is if interpolated volume is close to
c pixel volume - criterion is 3 % difference
c Less than 3 % difference between voxel volume and SH volume
c means that particle is determined to be SH. If more than 3 % 
c difference, then particle is determined to be nonSH and is stored
c voxel by voxel.

c Small black internal pores that particle
c might contain have already been cleaned out. 
      
      write(8,*) 'Now do anm process'
      write(8,*) 'Actual center of particle = ',xcm,ycm,zcm
      write(8,*) 'Pixel center of particle = ',icm,jcm,kcm
      write(8,*) ' Min pixel labels (ijk) ',ileft,jleft,kleft
      write(8,*) ' Max pixel labels (ijk) ',iright,jright,kright
      call flush(8)

      xr=iright-icm
      xl=icm-ileft
      yr=jright-jcm
      yl=jcm-jleft
      zr=kright-kcm
      zl=kcm-kleft
      lenx=2*max(xl,xr)
      leny=2*max(yl,yr)
      lenz=2*max(zl,zr)
      write(8,*) 'x range = ',lenx
      write(8,*) 'y range = ',leny
      write(8,*) 'z range = ',lenz
      call flush(8)

c  Interpolate r(theta,phi) and compute anm coefficients, output
c resolution used (mm/pixel length)

      write(8,197) sx,sy,sz
197   format(' Scale factors (xyz) = ',3f10.6)
      volpix=isave*sx*sy*sz
      write(8,*) 'Particle pixel vol*sx*sy*sz = ',volpix
      write(8,*) 'Voxel volume = ',isave
      if(isave.gt.100000) write(8,*) 'LARGE volume'
      call flush(8)

c  nnn is the number of anm coefficients to use
c probable good particle,since it passed all the previous tests
c go ahead and record anm's if its interpolated volume is within
c 3 % of the pixel volume
c Particle label is decremented by one if it turns out to be a bad particle 
      nnn=26
      numtot=numtot+1
      write(8,*) ' going into anm'
      call flush(8)
      call anm(nnn,p,y,a,yphi,ytheta,pix,xg,wg,ng,sx,sy,sz,lenx,leny,
     &lenz,nx,ny,nz,xcm,ycm,zcm,burned,volpix,volume,numtot,filestat)
      write(8,*) ' coming out of anm'
      call flush(8)
      vvv=100.*abs(volume-volpix)/volpix

c if particle does not pass the interpolated volume test, then drop the particle
c count by 1 and do not make the VRML image
c however, do count it as a non-SH particle and number it with ibad
c and output a voxpix image, named according to this number
      if(vvv.gt.3.) then
c use ibad to number the file
      ibad=ibad+1
      iibadd=30
      write(8,*) ' type of non-SH particle ',iibadd
      partvolR(ibad,1)=ipart
      partvolR(ibad,2)=isave
      partvolR(ibad,3)=iibadd
      write(8,*) ' non-SH particle number ',ibad
      write(8,*) ' voxel volume not equal to SH volume = ',isave
      write(8,*) ' going into voxpixR'
      call flush(8)
      call voxpixR(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,numtot,lenx,kright,sx,sy,sz,ibad,iibadd,filestat)
      write(8,*) ' coming out of voxpix'
      call flush(8)
      call nonSHout(isave,save,ibad,nnsave,iii,filestat,sx)
      numtot=numtot-1
c create filename using filestat and ibad and iibadd
      fnonSH(1:12)=filestat
      write(fnonSH(19:24),'(I6)') 100000+ibad
      write(fnonSH(30:32),'(I3)') 100+iibadd
      fnonSH(13:19)='-nonSH-'
      fnonSH(25:30)='-pore-'
      fnonSH(33:36)='.dat'
c record total porosity of particle, even if zero
       write(19,819) fnonSH,isave,sx,
     &     100.*bcount/float(isave)
       call flush(19)
c now output volumes of pores found in particle if bcount > 0
      if(bcount.gt.0.0) then
      open(unit=77,file=fnonSH)
      write(77,771) spnum
      do 5323 i44=1,spnum
      write(77,771) i44,savepnum(i44)
5323  continue
      close(77)
      end if
      goto 444
      end if

c since vvv.le.3, then particle is an SH particle, go ahead and record
c appropriately. anm's recorded, so make VRML images
      partvol(numtot,1)=ipart
      partvol(numtot,2)=isave
c now output volumes of pores found in particle if bcount > 0
c create filename using filestat and numtot
      fSH(1:12)=filestat
      write(fSH(16:21),'(I6)') 100000+numtot
      fSH(13:16)='-SH-'
      fSH(22:26)='-pore'
      fSH(27:30)='.dat'
c record porosity of particle, even if zero
      write(49,817) fSH,isave,sx,
     &     100.*bcount/float(isave)
817   format(a30,i10,f12.8,' %porosity=',f12.8)
       call flush(49)
      if(bcount.gt.0.0) then
      open(unit=77,file=fSH)
      write(77,771) spnum
      do 5324 i44=1,spnum
      write(77,771) i44,savepnum(i44)
5324  continue
      close(77)
      end if

c  Compute VRML file, output
c  First set resolution of VRML image
      ntheta=101
      nphi=101
      write(8,*) ' going into vrml'
      call flush(8)
      call vrmlout(p,xg,wg,ng,y,a,nnn,ntheta,nphi,numtot,filestat)
      write(8,*) ' coming out of vrml'
      call flush(8)

c Particle is a good one, anm's recorded, so make a VRML image of raw voxels
c particle is labelled in array save(isave,3)
c name of file looks like spherical harmonic VRML file, but with "vox-" in front

      write(8,*) ' going into voxpix'
      call flush(8)

      call voxpix(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,numtot,lenx,kright,sx,sy,sz,filestat)

      write(8,*) ' coming out of voxpix'
      call flush(8)


c  Wipe out burned pixels, so particle is not found again, set to
c  background (0 in this case)
      write(8,*) ' going into restore'
      call flush(8)
444   call restore(pix,0,burned,nx,ny,nz,save,isave,nnsave,iii)
      write(8,*) ' coming out of restore'
      call flush(8)

1000  continue

9900  continue

1111  write(8,*)
      icc=countd
      write(8,*) ' Value of vol cutoffs = ',icutoff,jcutoff
      write(8,*) ' Number of particles analyzed = ',numtot-numstar
c      write(8,*) ' Number of possible multi-particles = ',icc
      write(8,*) 
      write(8,*) 'Pixel labels,volumes of all SH particles analyzed'
      write(8,*) 
      numstar=numstar+1
      do 2345 ii=numstar,numtot
      write(8,1919) ii,partvol(ii,2)
      write(88,1919) ii,partvol(ii,2)
1919  format(i6,i10,i5)
2345  continue
      write(8,*) 'Pixel labels,volumes of all non-SH particles'
      write(8,*) ' non-SH particles found = ',ibad
      write(8,*) 
      call flush(8)
      write(8,*) ' 10 = CM 20 = rcomp>1 30 = anm volume'
      write(8,*) ' particle number - voxel volume - type of non-SH'
      write(8,*) 
      do 2445 ii=1,ibad
      write(8,1919) ii,partvolR(ii,2),partvolR(ii,3)
2445  continue
      call flush(8)

6666  continue
      end

      subroutine anm(nnn,p,y,a,yphi,ytheta,pix,xg,wg,ng,
     &         sx,sy,sz,lenx,leny,lenz,nx,ny,nz,xcm,ycm,zcm,burned,
     &          volpix,volume,numtot,filestat)

c  Interpolate surface function (rr(theta,phi)) of shape 
c  (0 = matrix, ipart = particle). Take regular increments of theta, phi
c  according to Gaussian quadrature points,
c  get rr(theta,phi) corresponding to these points.
c  Go along (theta,phi) vector (r unit vector), approximately 
c  find surface points by taking pairs inside and outside object.
c  Use this numerical characterization of r(theta,phi) to compute
c  the spherical harmonic coefficients (anm)
c  output into new file with labelling scheme

         real rr(1000,1000),len,p(0:200,0:500)
         real xg(1000),wg(1000),phi,pi
         real theta,lenx,leny,lenz
         complex y(0:200,-250:250),a(0:200,-250:250)
         complex yphi(0:200,-250:250)
         complex ytheta(0:200,-250:250)
         integer*2 pix(1024,1024,1024)
         integer*2 burned
         character*26 filename
         character*12 filestat

c sx, sy, sz are the scale factors

         pi=4.0*atan(1.0)
         write(8,*) ' filestat in anm = ',filestat
         call flush(8)

c  Interpolate r(theta,phi) for a regular arrangement of (theta,phi)
c  cm coords have been computed in subroutine fire
         xcm=xcm*sx
         ycm=ycm*sy
         zcm=zcm*sz
         write(8,*) ' xcm ycm zcm scaled in anm =',xcm,ycm,zcm

         len=min(lenx*sx,leny*sy,lenz*sz)
         if(len.eq.0.0) then
         len=0.5*max(lenx*sx,leny*sy,lenz*sz)
         write(8,*) ' len = 0 '
         end if
         write(8,*) ' len = ',len
         cutoff=len/100000.0
         write(8,*) ' cutoff in anm = ',cutoff
         write(8,*) ' start interpolation numtot = ',numtot
         call flush(8)

         do 210 i=1,ng
         do 200 j=1,ng
         theta=0.5*pi*xg(i)+0.5*pi
         phi=pi*xg(j)+pi
         phi=pi-phi
         dd=len/10.
         rstart=0.0

45       continue
         do 40 ii=1,10
         rr1=rstart+(ii-1)*dd
         rr2=rr1+dd
c  position with respect to center of mass
         x1=rr1*sin(theta)*cos(phi)
         y1=rr1*sin(theta)*sin(phi)
         z1=rr1*cos(theta)
         x2=rr2*sin(theta)*cos(phi)
         y2=rr2*sin(theta)*sin(phi)
         z2=rr2*cos(theta)

c  actual coordinate
         x1=x1+xcm
         x2=x2+xcm
         y1=y1+ycm
         y2=y2+ycm
         z1=z1+zcm
         z2=z2+zcm

c  nearest pixel center
         i1=int(x1/sx+0.5)
         i2=int(x2/sx+0.5)
         j1=int(y1/sy+0.5)
         j2=int(y2/sy+0.5)
         k1=int(z1/sz+0.5)
         k2=int(z2/sz+0.5)

c  check if 2nd point goes beyond image box, on the first try
           if(i1.gt.nx.or.i1.lt.1.or.j1.gt.ny.or.j1.lt.1.or.
     &      k1.gt.nz.or.k1.lt.1) then
             write(8,*) ' i,j,k beyond image bd for 1',i1,j1,k1
c             dd=dd/2.0
c             goto 45
           end if
           if(i2.gt.nx.or.i2.lt.1.or.j2.gt.ny.or.j2.lt.1.or.
     &      k2.gt.nz.or.k2.lt.1) then
             dd=dd/2.0
             write(8,*) ' i,j,k beyond image bd for 2 ',i2,j2,k2
             write(8,*) ' rr1 rr2 ',rr1,rr2
             goto 45
           end if

c  if pair of points are both within particle, go to next pair
      if(pix(i1,j1,k1).eq.burned.and.pix(i2,j2,k2).eq.burned) then
           if(ii.lt.10) goto 40
             if(ii.eq.10) then
             rstart=rr2
             goto 45
             end if
           end if

c  if pair crossed boundary, reduce step size, rename rstart at first point
c  of pair, and do over.  If step size is small enough, then end.
         if(pix(i1,j1,k1).eq.burned.and.pix(i2,j2,k2).ne.burned) then
          if(abs(dd).lt.cutoff) then
          goto 49
          end if
          rstart=rr1
          dd=dd/10.0
          goto 45
         end if

40       continue

49       rr(i,j)=0.5*(rr1+rr2)

200      continue
210      continue

c  compute spherical harmonic expansion

c  nnn is how many y's to use in series (n limit)
         a(0,0)=cmplx(0.,0.)
         do 390 n=1,nnn
         do 390 m=-n,n,1
         a(n,m)=cmplx(0.0,0.0)
390      continue

         factor=pi*pi/2.
         volume=0.0
         do 400 i=1,ng
         do 400 j=1,ng
         theta=0.5*pi*xg(i)+0.5*pi
         phi=pi*xg(j)+pi
         call harm(theta,phi,nnn,p,y)
         a(0,0)=a(0,0)+rr(i,j)*wg(i)*wg(j)*sin(theta)*
     &        cmplx(real(y(0,0)),-aimag(y(0,0))) 
         volume=volume+wg(i)*wg(j)*sin(theta)*
     &    (rr(i,j)**3)/3.0

         do 450 n=1,nnn
         do 450 m=n,-n,-1
         a(n,m)=a(n,m)+rr(i,j)*wg(i)*wg(j)*sin(theta)*
     &        cmplx(real(y(n,m)),-aimag(y(n,m))) 
450      continue
         
400      continue

         a(0,0)=a(0,0)*factor
         volume=volume*factor
         write(8,*) ' particle ',numtot,' inter. vol. = ',volume
c if volumes do not agree within 3 %, then bad interpolation
         vvv=100.*abs(volume-volpix)/volpix
         if(vvv.gt.3.) then
         write(8,*) ' volumes do not agree within 5 %, must be'
         write(8,*) ' bad interpolation, do not store anm coeffs.'
         return
         end if

      filename(1:12)=filestat
      write(filename(17:22),'(I6)') 100000+numtot
      filename(13:17)='-anm-'
      filename(23:26)='.dat'

c 15 is the unit for the anm data file to be created
         open(unit=15,file=filename)
         write(15,12) 0,0,a(0,0)
         do 460 n=1,nnn
         do 460 m=n,-n,-1
         a(n,m)=a(n,m)*factor
         write(15,12) n,m,a(n,m)
12       format(2i4,2f20.10)
460      continue
         close(15)
        
         return

         end

c  this is the factorial function
       function fac(j)
       integer j

       if(j.eq.0) then
       fac=1.
       return
       end if
       if(j.eq.1) then
       fac=1.
       return
       end if

       fac=1.
       do 10 i=1,j
       fac=fac*float(i)
10     continue
       return
       end

       subroutine harm(theta,phi,nnn,p,y)
c  computes spherical harmonics (complex)
c  for a value of x = cos(theta), phi=angle phi
c  so -1 < x < 1, P(n,m), -n < m < n, 0 < n .
c  uses two recursion relations, plus exact formulas
c  for the associated Legendre functions up to n=8

       real p(0:200,0:500),x,s,fac
       complex y(0:200,-250:250)
       real phi,theta,pi,xn,xm

       pi=4.0*atan(1.0)
       x=cos(theta)
       s=sqrt(1.-x*x)

       do 50 n=0,nnn
       do 50 m=0,n
       p(n,m)=0.0
50     continue

       p(0,0)=1.0
       p(1,0)=x
       p(1,1)=s
       p(2,0)=0.5*(3.*x*x-1.)
       p(2,1)=3.*x*s
       p(2,2)=3.*(1.-x*x)
       p(3,0)=0.5*x*(5.*x*x-3.)
       p(3,1)=1.5*(5.*x*x-1.)*s
       p(3,2)=15.*x*(1.-x*x)
       p(3,3)=15.*(s**3)
       p(4,0)=0.125*(35.*(x**4)-30.*x*x+3.)
       p(4,1)=2.5*(7.*x*x*x-3.*x)*s
       p(4,2)=7.5*(7.*x*x-1.)*(1.-x*x)
       p(4,3)=105.*x*(s**3)
       p(4,4)=105.*((1.-x*x)**2)
       p(5,0)=0.125*x*(63.*(x**4)-70.*x*x+15.)
       p(5,1)=0.125*15.*s*(21.*(x**4)-14.*x*x+1.)
       p(5,2)=0.5*105.*x*(1.-x*x)*(3.*x*x-1.) 
       p(5,3)=0.5*105.*(s**3)*(9.*x*x-1.)
       p(5,4)=945.*x*((1.-x*x)**2)
       p(5,5)=945.*(s**5)
       p(6,0)=0.0625*(231.*(x**6)-315.*(x**4)+105.*x*x-5.)
       p(6,1)=0.125*21.*x*(33.*(x**4)-30.*x*x+5.)*s
       p(6,2)=0.125*105.*(1.-x*x)*(33.*(x**4)-18.*x*x+1.)
       p(6,3)=0.5*315.*(11.*x*x-3.)*x*(s**3)
       p(6,4)=0.5*945.*(1.-x*x)*(1.-x*x)*(11.*x*x-1.)
       p(6,5)=10395.*x*(s**5)
       p(6,6)=10395.*((1.-x*x)**3)
       p(7,0)=0.0625*x*(429.*(x**6)-693.*(x**4)+315.*x*x-35.)
       p(7,1)=0.0625*7.*s*(429.*(x**6)-495.*(x**4)+135.*x*x-5.)
       p(7,2)=0.125*63.*x*(1.-x*x)*(143.*(x**4)-110.*x*x+15.)
       p(7,3)=0.125*315.*(s**3)*(143.*(x**4)-66.*x*x+3.)
       p(7,4)=0.5*3465.*x*(1.-x*x)*(1.-x*x)*(13.*x*x-3.)
       p(7,5)=0.5*10395.*(s**5)*(13.*x*x-1.)
       p(7,6)=135135.*x*(1.-x*x)*(1.-x*x)*(1.-x*x)
       p(7,7)=135135.*(s**7)
       p(8,0)=(1./128.)*(6435.*(x**8)-12012.*(x**6)+6930.*(x**4)-
     &   1260.*x*x+35.)
       p(8,1)=0.0625*9.*x*s*(715.*(x**6)-1001.*(x**4)+385.*x*x-35.)
       p(8,2)=0.0625*315.*(1.-x*x)*(143.*(x**6)-143.*(x**4)+33.*x*x-1.)
       p(8,3)=0.125*3465.*x*(s**3)*(39.*(x**4)-26.*x*x+3.)
       p(8,4)=0.125*10395.*(1.-x*x)*(1.-x*x)*(65.*(x**4)-26.*x*x+1.)
       p(8,5)=0.5*135135.*x*(s**5)*(5.*x*x-1.)
       p(8,6)=0.5*135135.*((1.-x*x)**3)*(15.*x*x-1.)
       p(8,7)=2027025.*x*(s**7)
       p(8,8)=2027025.*((1.-x*x)**4)

c  now generate spherical harmonics for n=0,8 (follows Arfken)
       do 300 n=0,8

c does n=0 separately
       if(n.eq.0) then
       y(0,0)=1./sqrt(4.*pi)
       goto 300
       end if

       do 200 m=n,-n,-1
       if(m.ge.0) then
       y(n,m)=((-1.)**m)*sqrt( ((2*n+1)/4./pi)*fac(n-m)/fac(n+m) )
     &    *p(n,m)*cmplx(cos(m*phi),sin(m*phi))
       goto 200
       end if
       if(m.lt.0) then
       mm=-m
       y(n,m)=((-1.)**mm)*cmplx(real(y(n,mm)),-aimag(y(n,mm)))
       goto 200
       end if
200    continue
300    continue
     
c  use recursion relations for n >= 9
c Do recursion on spherical harmonics, as they are better behaved numerically
c than the associated Legendre functions

       do 100 n=9,nnn
       do 99 m=0,n-2
       xn=float(n-1)
       xm=float(m)
       y(n,m)=(2.*xn+1.)*x*y(n-1,m)-sqrt( (2.*xn+1.)*(xn**2-xm**2)
     &        /(2.*xn-1.) )*y(n-2,m)
       y(n,m)=y(n,m)/sqrt( (2.*xn+1.)*((xn+1.)**2-xm**2)/(2.*xn+3.) )
99     continue
c  now do top m=n exactly
       nn=2*n-1
       p(n,n)=s**n
       do 97 i=1,nn,2
       p(n,n)=p(n,n)*float(i)
97     continue
       y(n,n)=((-1.)**n)*sqrt( ((2*n+1)/4./pi)*fac(n-n)/fac(n+n) )
     &    *p(n,n)*cmplx(cos(n*phi),sin(n*phi))

c  now do second to the top m=n-1 using the exact m=n, and the 
c  recursive m=n-2 found previously
       xm=float(n-1)
       xn=float(n)
       y(n,n-1)=-cmplx(cos(phi),sin(phi))*y(n,n-2)*
     &          ( xn*(xn+1.)-xm*(xm-1.) )
     &         /sqrt((xn+xm)*(xn-xm+1.))
       y(n,n-1)=y(n,n-1)-cmplx(cos(phi),-sin(phi))*y(n,n)
     &           *sqrt((xn-xm)*(xn+xm+1.)) 
       y(n,n-1)=y(n,n-1)*s/2./xm/x
100    continue
c now fill in -m terms
       do 220 n=0,nnn
       do 220 m=-1,-n,-1
       mm=-m
       y(n,m)=((-1.)**mm)*cmplx(real(y(n,mm)),-aimag(y(n,mm)))
220    continue

2222   continue

       return
       end

       subroutine surharm(theta,phi,nnn,p,y,yphi,ytheta,yphiphi,
     &                     ythth,yphith)

c  computes spherical harmonics (complex)
c  for a value of x = cos(theta), phi=angle phi
c  so -1 < x < 1, P(n,m), -n < m < n, 0 < n .
c  uses the spherical harmonics calculated in harm 

       real p(0:200,0:500),x,s,phi
       complex y(0:200,-250:250),yphith(0:200,-250:250)
       complex yphi(0:200,-250:250),yphiphi(0:200,-250:250)
       complex ytheta(0:200,-250:250),ythth(0:200,-250:250)
       real pqa(100),theta
       integer ipqa(100)

       pi=4.0*atan(1.0)
       x=cos(theta)
       s=sqrt(1.-x*x)

c  Define various partial derivatives of the y(n,m)'s using the y(n,m)
c  functions themselves (following Arfken)

2222   continue

       do 300 n=0,nnn

c does n=0 separately, derivative of n=0 constant term is identically zero
       if(n.eq.0) then
       yphi(0,0)=cmplx(0.0,0.0)
       yphiphi(0,0)=cmplx(0.0,0.0)
       yphith(0,0)=cmplx(0.0,0.0)
       ytheta(0,0)=cmplx(0.0,0.0)
       ythth(0,0)=cmplx(0.0,0.0)
       goto 300
       end if

c  do m=n separately and exactly, below works for all other values of m
       do 200 m=n-1,-n,-1
       xn=float(n)
       xm=float(m)
       aaa=sqrt((2.*xn+1.)*(xn**2-xm**2)/(2.*xn-1.))
       yphi(n,m)=cmplx(0.,1.)*m*y(n,m)
       yphiphi(n,m)=-(m**2)*y(n,m)
       ytheta(n,m)=xn*x*y(n,m)/s-aaa*y(n-1,m)/s
       ythth(n,m)=((xm*xm-xn*x*x)/s/s-xn*(xn+1.))*y(n,m)+
     &             aaa*x*y(n-1,m)/s/s
       yphith(n,m)=ytheta(n,m)*cmplx(0.,1.)*m
200    continue
c  now do m=n exactly
       m=n
       xn=float(n)
       xm=float(m)
       yphi(n,m)=cmplx(0.,1.)*m*y(n,m)
       yphiphi(n,m)=-(m**2)*y(n,m)
       ytheta(n,m)=xn*x*y(n,m)/s
       ythth(n,m)=xn*(xn*x*x-1.)*y(n,m)/s/s
       yphith(n,m)=ytheta(n,m)*cmplx(0.,1.)*m
300    continue
     
       return
       end
       subroutine surfarea(pix,nx,ny,nz,surface)
       integer*2 pix(1024,1024,1024)
       
c  compute surface area
      surface=0.0
      do 300 k=1,nz
      do 300 j=1,ny
      do 300 i=1,nx 
      i2=i-1
      j2=j-1
      k2=k-1
      if(i2.lt.1) i2=i2+nx
      if(j2.lt.1) j2=j2+ny
      if(k2.lt.1) k2=k2+nz
      m1=nxy*(k-1)+nx*(j-1)+i2
      m2=nxy*(k-1)+nx*(j2-1)+i
      m3=nxy*(k2-1)+nx*(j-1)+i
      if(pix(i,j,k).eq.1.and.pix(i2,j,k).eq.2) then
      surface=surface+1.0
      end if
      if(pix(i,j,k).eq.2.and.pix(i2,j,k).eq.1) then
      surface=surface+1.0
      end if
      if(pix(i,j,k).eq.1.and.pix(i,j2,k).eq.2) then
      surface=surface+1.0
      end if
      if(pix(i,j,k).eq.2.and.pix(i,j2,k).eq.1) then
      surface=surface+1.0
      end if
      if(pix(i,j,k).eq.1.and.pix(i,j,k2).eq.2) then
      surface=surface+1.0
      end if
      if(pix(i,j,k).eq.2.and.pix(i,j,k2).eq.1) then
      surface=surface+1.0
      end if
300   continue

      return
      end
      subroutine vrmlout(p,xg,wg,ng,y,a,nnn,ntheta,nphi,ipart,filestat)
c  Reads in anm coefficients, outputs VRML file for visualization

         real rr(1000,1000),len,p(0:200,0:500)
         real xg(1000),wg(1000),phi,pi,rsave(200,200)
         real theta,I11,I22,I33,I12,I13,I23
         complex y(0:200,-250:250),a(0:200,-250:250)
         complex r
         character*22 filename
         character*12 filestat

         pi=4.0*atan(1.0)

c  note - in this subroutine, I use ipart but the value of numtot is 
c  passed to the subroutine
      filename(1:12)=filestat
      write(filename(13:18),'(I6)') 100000+ipart
      filename(13:13)='-'
      filename(19:22)='.wrl'

c 17 is the unit for the VRML file to be created
         open(unit=17,file=filename)

c write out header of VRML file

         write(17,11)
11       format('#VRML V2.0 utf8')
         write(17,12)
12       format('NavigationInfo {')
         write(17,13)
13       format('type ["EXAMINE","WALK","FLY"]')
         write(17,14)
14       format('}')
         write(17,15)
15       format('Group {')
         write(17,16)
16       format('children [')
         write(17,17)
17       format('Shape {')
         write(17,18)
18       format('geometry IndexedFaceSet {') 
         write(17,19)
19       format('solid TRUE') 
         write(17,20) 
20       format('ccw FALSE') 
         write(17,21)
21       format('coord Coordinate{') 
         write(17,22)
22       format('point [') 

         zmax=0.0
         do 420 i=1,ntheta
         do 420 j=1,nphi
         theta=(float(i)-1.0)*pi/float(ntheta)
         if(i.eq.1) theta=0.001*pi
         if(i.eq.ntheta) theta=0.999*pi
         phi=2.*pi*(float(j)-1.0)/float(nphi)
         call harm(theta,phi,nnn+2,p,y)

         r=a(0,0)*y(0,0)
         do 82 n=1,nnn
         do 82 m=n,-n,-1
         r=r+a(n,m)*y(n,m)
82       continue
         xx=real(r)*sin(theta)*cos(phi)
         yy=real(r)*sin(theta)*sin(phi)
         zz=real(r)*cos(theta)
         if(zmax.lt.zz) zmax=zz

         if(i.eq.1.and.j.eq.1) then
         write(17,97) 0.,0.,zz
         end if
         write(17,97) xx,yy,zz
97       format(3f20.10)
         if(i.eq.ntheta.and.j.eq.nphi) write(17,97) 0.,0.,zz

420      continue
         write(17,23)
23       format(']')
         write(17,24)
24       format('}')
         write(17,26)
26       format('coordIndex [')

c first do top end cap in triangles
         do 1900 j=1,nphi-1
         write(17,*) 0,j+1,j,-1 
1900     continue
         write(17,*) 0,1,ntheta,-1
c  done with top end cap

         do 1800 i=1,ntheta-1
         do 1700 j=1,nphi
         number=ntheta*(i-1)+j
         if(j.eq.nphi) then
         num=number+1-nphi
         write(17,*) number,num,num+ntheta,number+ntheta,-1
         goto 1700
         end if
         write(17,*) number,number+1,number+1+ntheta,number+ntheta,-1
1700     continue
1800     continue

c finally do bottom end cap in triangles
         number=ntheta*nphi+1
         do 1910 j=1,nphi-1
         num=(ntheta-1)*nphi+j
         write(17,*) number,num,num+1,-1
1910     continue
         write(17,*) number,number-1,(ntheta-1)*nphi+1,-1
c  done with bottom end cap
         write(17,27)
27       format(']')
         write(17,47)
47       format('creaseAngle 0.8')
c         write(17,28)
c28       format('normal Normal {')
c         write(17,29)
c29       format('vector [')
c**** x y z of normal vectors********
c         write(17,30) 
c30       format(']')
c         write(17,31)
c31       format('}')
         write(17,32)
32       format('}')
         write(17,33)
33       format('appearance Appearance {')
         write(17,34)
34       format('material Material {')
         write(17,35)
35       format('diffuseColor 0.923106 0.923106 0.923106')
         write(17,36)
36       format('}')
         write(17,36)
         write(17,36)
         write(17,37)
37       format('Viewpoint {')
         write(17,38) 0.0,0.0,7.0*zmax
38       format('position ',2x,f6.4,2x,f6.4,2x,f9.4)
         write(17,36)
         write(17,39)
39       format(']')
         write(17,36)

c         close(17)
         return
         end
c This subroutine restores the burned pixels of a particular particle
c (= burned) back to their original, unburned
c value (phase). This is only done for a particular particle to avoid 
c doing it for the entire system, for each particle.
c Subroutine either wipes out a particle that was too small, or
c sets the value of a good particle to its new negative value.

      subroutine restore(pix,phase,burned,nx,ny,nz,save,
     &                   isave,nnsave,iii)
      integer*2 pix(1024,1024,1024),save(nnsave,iii)
      integer*2 burned,phase

      do 10 m=1,isave
      i=save(m,1)
      j=save(m,2)
      k=save(m,3)
      if(pix(i,j,k).eq.burned) pix(i,j,k)=phase
10    continue

      return
      end
      subroutine small(pix,pixsmall,nx,ny,nz,nxs,nys,nzs,ileft,jleft,
     &            kleft,iright,jright,kright,burned,nnarray,
     &            old,new,iii,isave,save,nnsave,bcount,ipore,
     &            sx,sy,sz,porevol,filestat,savepnum,spnum)

      integer*2 pix(1024,1024,1024),pixsmall(nxs,nys,nzs)
      integer*2 old(nnarray,iii),new(nnarray,iii)
      integer*2 phase,burned,irestore,tempburn,poreburn
      integer*2 in(6),jn(6),kn(6)
c      integer*2 save(nnsave,iii),savep(nnsave,iii)
      integer*2 save(nnsave,iii),savep(1000000,4)
      integer*4 inew,iold
      integer*4 porevol(10000000),savepnum(1000000),spnum
      character*12 filestat

c principle of subroutine:  Start with a box that is two bigger in each
c direction than ileft, iright, etc. Turn all outer planes into 0, then 
c build a "particle" or fire inward from each side.
c When "particle" or fire is stopped,
c then just need to turn inner black pores into burned particle phase,
c and compare pixsmall to pix. If pix=0 and pixsmall=burned, then pix=burned.
c Otherwise, no change in pix.

      pi=4.0*atan(1.0)
c  Direction labels to check for burning path (nearest neighbor information
c  in 3-D digital system).
      in(1)=-1
      in(2)=1
      in(3)=0
      in(4)=0
      in(5)=0
      in(6)=0
      jn(1)=0
      jn(2)=0
      jn(3)=-1
      jn(4)=1
      jn(5)=0
      jn(6)=0
      kn(1)=0
      kn(2)=0
      kn(3)=0
      kn(4)=0
      kn(5)=-1
      kn(6)=1

c initialize pixsmall to pore

      do 777 k=1,nzs
      do 777 j=1,nys
      do 777 i=1,nxs
      pixsmall(i,j,k)=0
777   continue

c Map appropriate volume of pix into pixsmall
      kr=kright+1-(kleft-1)+1
      jr=jright+1-(jleft-1)+1
      ir=iright+1-(ileft-1)+1
      write(8,*) ' in subroutine small'
      write(8,*) ' ir jr kr = ',ir,jr,kr
      write(8,*) ' ileft iright = ',ileft,iright
      write(8,*) ' jleft jright = ',jleft,jright
      write(8,*) ' kleft kright = ',kleft,kright
      write(8,*) ' pix(ileft,jleft,kleft)= ',pix(ileft,jleft,kleft)

      check=0.0
      checkpix=0.0
      do 333 k=kleft-1,kright+1
      do 333 j=jleft-1,jright+1
      do 333 i=ileft-1,iright+1
      i1=i-(ileft-1)+1
      j1=j-(jleft-1)+1
      k1=k-(kleft-1)+1
      if(pix(i,j,k).eq.burned) checkpix=checkpix+1.0
      if(pix(i,j,k).eq.burned.or.pix(i,j,k).eq.0) then
      pixsmall(i1,j1,k1)=pix(i,j,k)
      end if
      if(pixsmall(i1,j1,k1).eq.burned) check=check+1.0
c      if(pix(i,j,k).eq.1) then
c      write(8,*) i,j,k,i1,j1,k1,pix(i,j,k)
c      end if
333   continue
      write(8,*) ' burned vol of pix = ',checkpix
      write(8,*) ' after pix-pixsmall mapping, burned vol = ',check

      tempburn=-1500
      poreburn=-2500
      write(8,*) ' burned = ',burned
      write(8,*) ' tempburn = ',tempburn
      write(8,*) ' poreburn = ',poreburn
      write(8,*) 
c  find first matrix pixels around faces of particle system
      do 990 k=1,kr
      do 990 j=1,jr
      i=1
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      i=ir
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
990   continue
      do 991 k=1,kr
      do 991 i=1,ir
      j=1
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      j=jr
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
991   continue
      do 992 j=1,jr
      do 992 i=1,ir
      k=1
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      k=kr
      pixsmall(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
992   continue

c  Now start building up new pixels from old set of pixels,
c  thus propagating the fire inward and identifying new particle 
c  pixels connnected to the ones already burned.

      icount=1
60    inew=0
      do 100 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
      kk=old(ijk,3)
c  check all six nearest neighbors of previously burned pixel
      do 90 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

c  Hard boundary conditions (can't pass beyond boundary)

      if(k1.lt.1.or.k1.gt.kr) goto 90
      if(j1.lt.1.or.j1.gt.jr) goto 90
      if(i1.lt.1.or.i1.gt.ir) goto 90

c  Store (i,j,k) labels of newly burned pixels in array new().
      if(pixsmall(i1,j1,k1).ne.burned.and.
     &      pixsmall(i1,j1,k1).ne.tempburn) then
      pixsmall(i1,j1,k1)=tempburn
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      new(inew,3)=k1
      end if
90    continue
100   continue

c  If new pixels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 150 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
      old(ijk,3)=new(ijk,3)
150   continue
      icount=icount+1
      goto 60
      end if

6060    continue

c do quick count of pores, if none, go ahead and do the 
c finishing steps. If some exist, then need to identify the individual particles,
c change them to something else (poreburn) and record pores as
c particles.

c *****
c following step has not been done at present, numbering of 
c pore as particles is independent of the particle that they came from.
c Have to figure out how to number pores tied to the number of the 
c particle they came out of. Shouldn't be too many in each particle, so
c allow for up to 100.
c ****

c do quick count of total pore volume in this particle
c set everything that is not burned or tempburn to 0=pore
c this is to prevent something funny, like a small solid particle
c inside the pore, to cause trouble
       bcount=0.0
       acount=0.0
       dcount=0.0
       do 225 k=1,kr
       do 225 j=1,jr
       do 225 i=1,ir
      if(pixsmall(i,j,k).ne.burned.and.pixsmall(i,j,k).ne.tempburn) then
       pixsmall(i,j,k)=0
       end if
       if(pixsmall(i,j,k).eq.0) then
       bcount=bcount+1.0
       end if
       if(pixsmall(i,j,k).eq.burned) then
       acount=acount+1.0
       end if
       if(pixsmall(i,j,k).eq.tempburn) then
       dcount=dcount+1.0
       end if
225    continue
       write(8,*) ' after tempburn bcount (pores) = ',bcount
       write(8,*) ' after tempburn tempburn = ',dcount
       write(8,*) ' after tempburn burned = ',acount

c if no pores in this particle, then jump to 5500 and finish subroutine small
       if(bcount.eq.0) then
       write(8,*) ' no internal pores found '
       goto 5500
       end if
c if there are pores, then need to count and record each one
c don't need to do SH analysis, just record voxels

c look for an internal particle pore, set pix=0 to poreburn as found

       spnum=0
       do 250 k=1,kr
       do 250 j=1,jr
       do 250 i=1,ir
       if(pixsmall(i,j,k).ne.0) goto 250
       if(pixsmall(i,j,k).eq.0) then
       ipart=pixsmall(i,j,k)
        write(8,*) 
        write(8,*) ' found pore in particle '
        write(8,*) ' at i j k = ',i,j,k,' ipart = ',ipart
        call flush(8)
        ipore=ipore+1
       end if

      iold=0
      isavep=0
c  find a first particle pixel
      if(pixsmall(i,j,k).eq.ipart) then
      pixsmall(i,j,k)=poreburn
      iold=iold+1
      isavep=isavep+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
c save particle voxel coordinates
      savep(iold,1)=i
      savep(iold,2)=j
      savep(iold,3)=k
      savep(iold,4)=0
      end if

c  Now start building up new pore pixels from old set of pixels,
c  thus propagating the fire inward and identifying new pore particle 
c  pixels connnected to the ones already burned.

      icount=1
61    inew=0
      do 101 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
      kk=old(ijk,3)
c  check all six nearest neighbors of previously burned pixel
      do 91 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

c  Hard boundary conditions (can't pass beyond boundary)

      if(k1.lt.1.or.k1.gt.kr) goto 91
      if(j1.lt.1.or.j1.gt.jr) goto 91
      if(i1.lt.1.or.i1.gt.ir) goto 91

c  Store (i,j,k) labels of newly burned pixels in array new().
      if(pixsmall(i1,j1,k1).eq.0) then
      pixsmall(i1,j1,k1)=poreburn
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      new(inew,3)=k1
      end if
91    continue
101   continue

c  If new pixels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 151 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
      old(ijk,3)=new(ijk,3)
      savep(ijk+isavep,1)=new(ijk,1)
      savep(ijk+isavep,2)=new(ijk,2)
      savep(ijk+isavep,3)=new(ijk,3)
      savep(ijk+isavep,4)=0
151   continue
      isavep=isavep+inew
      goto 61
      end if
c record size of pore in savepnum, increment spnum by 1
      if(isavep.lt.8)then
      ipore=ipore-1
      goto 555
      end if
      savepnum(ipore)=isavep
      spnum=spnum+1
      porevol(ipore)=isavep
      write(89,*) ipore,isavep,(6.*isavep/pi)**(1./3.)
      call flush(89)

555   continue
c now a pore as particle has been recorded. For vrmlpore, need
c to calculate icm, jcm, kcm and kright=largest value of k in
c pore as particle

c      xcm=0.0
c      ycm=0.0
c      zcm=0.0
c      kright=0
c      do 7300 m=1,isavep
c      xcm=xcm+savep(m,1)
c      ycm=ycm+savep(m,2)
c      zcm=zcm+savep(m,3)
c      if(kright.lt.savep(m,3)) kright=savep(m,3)
c7300  continue
c      xcm=0.5+xcm/float(isavep)
c      ycm=0.5+ycm/float(isavep)
c      zcm=0.5+zcm/float(isavep)
cc integer version of cm of pore as particle
c      icm=xcm
c      jcm=ycm
c      kcm=zcm

c poreout outputs the voxel locations of the pore that was found
c      call poreout(isavep,savep,ipore,nnsave,iii,filestat)
c vrmlpore outputs the VRML voxel image of the pore that was found
c      call vrmlpore(pixsmall,ir,jr,kr,icm,jcm,kcm,poreburn,savep,
c     &    isavep,nnsave,iii,kright,sx,sy,sz,ipore,filestat)

250   continue

c note - any small black pores must be internal, since outside
c was filled in by tempburn from the outer walls inward.
c therefore, any new particle voxels added cannot be on the surface,
c so that save(n,4) does not need to be updated (whether it is a 
c surface voxel or not)

c  Found all of matrix outside of particle, now further process it by
c  changing all black pore particles (0) to particle, 
c  since they can only be inside the actual particle
c  (fills in small black holes in  main particle).
c Change the appropriate values of pix at the same time.

c regular way to finish after any pores have been counted
5500   bcount=0.0
       iisave=0
       do 200 k=1,kr
       do 200 j=1,jr
       do 200 i=1,ir
       ii=i+ileft-2
       jj=j+jleft-2
       kk=k+kleft-2
       if(kk.lt.1.or.kk.gt.nz.or.jj.lt.1.or.jj.gt.ny.or.ii.lt.1.or
     &    .ii.gt.nx) write(8,*) ii,jj,kk,' ii jj kk'
       if(pixsmall(i,j,k).eq.0.or.pixsmall(i,j,k).eq.poreburn) then
       pix(ii,jj,kk)=burned
       bcount=bcount+1.0
       isave=isave+1
       save(isave,1)=ii
       save(isave,2)=jj
       save(isave,3)=kk
       save(isave,4)=0
       end if
       if(pix(ii,jj,kk).eq.burned) iisave=iisave+1
200    continue

      write(8,*) bcount,' small black pore pixels removed'
      write(8,*) ' isave = ',isave,' iisave = ',iisave
      return
      end
      subroutine voxpix(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,numtot,lenx,kright,sx,sy,sz,filestat)

      integer*4 face(500000000,4)
      real coord(500000000,3)
      real lenx
      integer*2 pix(1024,1024,1024)
      integer*2 burned
      integer*2 save(nnsave,iii)
      character*50 filename
      character*12 filestat

c scan over all voxels in particle

c mp is number of points, nf is number of faces
c note: I will have many duplicates of coords, so what, as long
c as there is unique labelling
      mp=0
      nf=0
      icm=icm-1.7*lenx
      do 1000 m=1,isave
      i=save(m,1)
      j=save(m,2)
      k=save(m,3)

c +i surfaces
      if(pix(i+1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -i surfaces
      if(pix(i-1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c +j surfaces
      if(pix(i,j+1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c -j surfaces
      if(pix(i,j-1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c +k surfaces
      if(pix(i,j,k+1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -k surfaces
      if(pix(i,j,k-1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if
1000  continue

c now output the VRML file (it is simple)

c write out header of VRML file

         write(17,11)
11       format('#VRML V2.0 utf8')
c         write(17,12)
c12       format('NavigationInfo {')
c         write(17,13)
c13       format('type ["EXAMINE","WALK","FLY"]')
c         write(17,14)
c14       format('}')
         write(17,17)
17       format('Shape {')
         write(17,18)
18       format('geometry IndexedFaceSet {') 
         write(17,19)
19       format('coord Coordinate{point[')
         do 100 m=1,mp
      if(m.lt.mp) write(17,20) -sx*coord(m,1),sy*coord(m,2),
     &        sz*coord(m,3)
      if(m.eq.mp) write(17,21) -sx*coord(m,1),sx*coord(m,2),
     &        sz*coord(m,3)
20       format(3f12.6,',')
21       format(3f12.6)
100      continue
         write(17,22)
22       format(']}')
         write(17,23)
23       format('coordIndex[')
         do 200 n=1,nf
      if(n.lt.nf) write(17,24) face(n,1),face(n,2),face(n,3),face(n,4)
      if(n.eq.nf) write(17,25) face(n,1),face(n,2),face(n,3),face(n,4)
24       format(4i8,' -1,')
25       format(4i8,']}')
200      continue
         write(17,26)
26       format('appearance Appearance{')
         write(17,27)
27       format('material Material{')
         write(17,28)
28       format('diffuseColor 0.923 0.923 0.96 }}}')
         zzk=4.*sz*(kright+1-kcm)
         write(17,29) 0.,0.,zzk
29       format('Viewpoint{position ',3f8.2,'}')        
         close(17)

         return
         end
      subroutine LWS(pix,nx,ny,nz,lwss,nnarray,iii,old,new)

c**************************  LWS.f  ********************************
c  This program works on a 3-D file and applies a limited watershed
c  split in 3-D.
c*******************************************************************

c  (USER) Dimensions of main arrays
      integer*2 pix(1024,1024,1024),old(nnarray,iii),new(nnarray,iii)
      integer*8 ns,mstep(7),LL(6)
      integer*2 in(6),jn(6),kn(6)
      integer*4 isave,inew,iold
      integer*2 array(7,100000000,4)

      ns=nx*ny*nz
c  Direction labels to check for burning path (nearest neighbor information
c  in 3-D digital system).
      in(1)=-1
      in(2)=1
      in(3)=0
      in(4)=0
      in(5)=0
      in(6)=0
      jn(1)=0
      jn(2)=0
      jn(3)=-1
      jn(4)=1
      jn(5)=0
      jn(6)=0
      kn(1)=0
      kn(2)=0
      kn(3)=0
      kn(4)=0
      kn(5)=-1
      kn(6)=1

c This part erodes and stores eroded voxels
c Assumes that system is segmented

      isum=0
      do 4500 mmm=1,lwss
      mstep(mmm)=0
      do 4000 k=1,nz
      do 4000 j=1,ny
      do 4000 i=1,nx
      if(pix(i,j,k).ne.1) goto 4000 
      do 490 n=1,6
      i1=i+in(n)
      j1=j+jn(n)
      k1=k+kn(n)
      if(k1.lt.1.or.k1.gt.nz) then
      goto 490
      end if
      if(j1.lt.1.or.j1.gt.ny) then
      goto 490
      end if
      if(i1.lt.1.or.i1.gt.nx) then
      goto 490
      end if

      if(pix(i1,j1,k1).eq.0) then
      mstep(mmm)=mstep(mmm)+1
      array(mmm,mstep(mmm),1)=i
      array(mmm,mstep(mmm),2)=j
      array(mmm,mstep(mmm),3)=k
      goto 4000
      end if
490   continue

4000  continue

c now remove eroded voxels
      mn1=1
      mn2=mstep(mmm)
      do 4100 mn=mn1,mn2
      ii=array(mmm,mn,1)
      jj=array(mmm,mn,2)
      kk=array(mmm,mn,3)
      pix(ii,jj,kk)=0
4100  continue

      write(8,*) ' mmm = ',mmm,' no. voxels removed = ',mstep(mmm)
      isum=isum+mstep(mmm)
4500  continue
      write(8,*) ' total no. voxels removed = ',isum
      
c This next part of the program does the actual burning = particle
c identification and labelling. 
c Need to identify each particle and give each a unique label

      label=0
      isum=0

      do 1000 k=1,nz
      do 1000 j=1,ny
      do 1000 i=1,nx

      if(pix(i,j,k).le.0) goto 1000
      label=label-1
      ipart=pix(i,j,k)

      iold=0
      isave=0
c  find a first particle voxel
      if(pix(i,j,k).eq.ipart) then
      pix(i,j,k)=label
      iold=iold+1
      isave=isave+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      ntouch=0
      end if

c  Now start building up new burned voxels from old set of burned voxels,
c  thus propagating the fire and identifying new particle voxels connnected
c  to the ones already burned.

60    inew=0
      do 100 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
      kk=old(ijk,3)
c  check all six nearest neighbors of previously burned voxel
      do 90 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

c  Hard boundary conditions (can't pass beyond boundary)

      if(k1.lt.1.or.k1.gt.nz) then
      goto 90
      end if
      if(j1.lt.1.or.j1.gt.ny) then
      goto 90
      end if
      if(i1.lt.1.or.i1.gt.nx) then
      goto 90
      end if

c  Store (i,j,k) labels of newly burned voxels in array new().
      if(pix(i1,j1,k1).eq.ipart) then
      pix(i1,j1,k1)=label
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      new(inew,3)=k1
      end if
90    continue
100   continue

c  If new voxels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 150 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
      old(ijk,3)=new(ijk,3)
150   continue
      goto 60
      end if

1000  continue

c all particles have been labelled now, with unique negative numbers
c pore space is still 0

c now put back eroded voxels, in reverse order, but don't allow any new connections
c between labelled particles. Put back voxels one at a time, so that each
c new voxel checks existing ones plus the ones that have been put
c back previously.
c what if all the neighbors of the replaced voxel are pore? Then don't
c put it back at all, since it is just noise, since it does not connect
c to any existing particle
c it is also possible that voxel only connects to one kind of particle
c so should be replaced
 
      do 8000 mmm=lwss,1,-1
      mn1=1
      mn2=mstep(mmm)
      do 7900 mn=mn1,mn2
      ii=array(mmm,mn,1)
      jj=array(mmm,mn,2)
      kk=array(mmm,mn,3)

c  check all six nearest neighbors of replaced voxel
      icheck=0
      mylabel=0
      do 690 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

      LL(n)=1000000

c  Hard boundary conditions (can't pass beyond boundary)

      if(k1.lt.1.or.k1.gt.nz) then
      goto 690
      end if
      if(j1.lt.1.or.j1.gt.ny) then
      goto 690
      end if
      if(i1.lt.1.or.i1.gt.nx) then
      goto 690
      end if

      LL(n)=pix(i1,j1,k1)
      if(LL(n).lt.0) mylabel=LL(n)
690   continue
c if all the neighbors were pores, then mylabel still equals 0

c if two kinds of particles are neighbors, then don't put this voxel back
c icheck = 0 means put it back, icheck = 1 means don't put it back

      icheck=0
      do 691 n=1,6
      if(LL(n).lt.0) then

       do 692 nn=1,6
       if(LL(nn).ge.0) goto 692
       if(LL(n).ne.LL(nn)) then
       icheck=1
       goto 693
       end if
692    continue

      end if

691   continue

c if all neighbors were 0, then icheck and mylabel remain 0
c LL(n)=1000000 if all neighbors were 0
693   if(icheck.eq.0) then
        pix(ii,jj,kk)=mylabel
        if(mylabel.ne.0) isum=isum+1
      end if

7900  continue
8000  continue

      do 9000 k=1,nz
      do 9000 j=1,ny
      do 9000 i=1,nx
      if(pix(i,j,k).lt.0) pix(i,j,k)=1
9000  continue
      write(8,*) ' total no. voxels replaced = ',isum

      return
      end
      subroutine voxpixR(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,numtot,lenx,kright,sx,sy,sz,ibad,iibadd,filestat)
      integer*4 face(500000000,4),ibad
      real coord(500000000,3)
      real lenx
      integer*2 pix(1024,1024,1024)
      integer*2 burned
      integer*2 save(nnsave,iii)
      character*31 filename
      character*12 filestat

c scan over all voxels in particle

c mp is number of points, nf is number of faces
c note: I will have many duplicates of coords, so what, as long
c as there is unique labelling
      mp=0
      nf=0
c      icm=icm-1.7*lenx
      do 1000 m=1,isave
      i=save(m,1)
      j=save(m,2)
      k=save(m,3)

c +i surfaces
      if(pix(i+1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -i surfaces
      if(pix(i-1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c +j surfaces
      if(pix(i,j+1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c -j surfaces
      if(pix(i,j-1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c +k surfaces
      if(pix(i,j,k+1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -k surfaces
      if(pix(i,j,k-1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if
1000  continue

c now output the VRML file (it is simple)

      filename(1:12)=filestat
      write(filename(19:21),'(I3)') 100+iibadd
      filename(13:19)='-nonSH-'
      write(filename(22:27),'(I6)') 100000+ibad
      filename(22:22)='-'
      filename(28:31)='.wrl'

c  I want file names involving system type, non-SH particle number
c (ibad), and type of non-SH (iibadd = 10, 20, or 30)
c 17 is the unit for the VRML file to be created
         open(unit=17,file=filename)

c write out header of VRML file

         write(17,11)
11       format('#VRML V2.0 utf8')
c         write(17,12)
c12       format('NavigationInfo {')
c         write(17,13)
c13       format('type ["EXAMINE","WALK","FLY"]')
c         write(17,14)
c14       format('}')
         write(17,17)
17       format('Shape {')
         write(17,18)
18       format('geometry IndexedFaceSet {') 
         write(17,19)
19       format('coord Coordinate{point[')
         do 100 m=1,mp
      if(m.lt.mp) write(17,20) -sx*coord(m,1),sy*coord(m,2),
     &        sz*coord(m,3)
      if(m.eq.mp) write(17,21) -sx*coord(m,1),sx*coord(m,2),
     &        sz*coord(m,3)
20       format(3f12.6,',')
21       format(3f12.6)
100      continue
         write(17,22)
22       format(']}')
         write(17,23)
23       format('coordIndex[')
         do 200 n=1,nf
      if(n.lt.nf) write(17,24) face(n,1),face(n,2),face(n,3),face(n,4)
      if(n.eq.nf) write(17,25) face(n,1),face(n,2),face(n,3),face(n,4)
24       format(4i8,' -1,')
25       format(4i8,']}')
200      continue
         write(17,26)
26       format('appearance Appearance{')
         write(17,27)
27       format('material Material{')
         write(17,28)
28       format('diffuseColor 0.923 0.923 0.96 }}}')
         zzk=4.*sz*(kright+1-kcm)
         write(17,29) 0.,0.,zzk
29       format('Viewpoint{position ',3f8.2,'}')        
         close(17)

         return
         end
      subroutine vrmlpore(pix,nx,ny,nz,icm,jcm,kcm,burned,save,isave,
     &    nnsave,iii,kright,sx,sy,sz,ibad,filestat)
      integer*4 face(50000000,4),ibad
      real coord(50000000,3)
      integer*2 pix(1024,1024,1024)
      integer*2 burned
      integer*2 save(nnsave,iii)
      character*27 filename
      character*12 filestat

c scan over all voxels in particle

c mp is number of points, nf is number of faces
c note: I will have many duplicates of coords, so what, as long
c as there is unique labelling
      mp=0
      nf=0
      do 1000 m=1,isave
      i=save(m,1)
      j=save(m,2)
      k=save(m,3)

c +i surfaces
      if(pix(i+1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -i surfaces
      if(pix(i-1,j,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c +j surfaces
      if(pix(i,j+1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if

c -j surfaces
      if(pix(i,j-1,k).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c +k surfaces
      if(pix(i,j,k+1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k+1-kcm
      face(nf,4)=mp-1
      end if

c -k surfaces
      if(pix(i,j,k-1).ne.burned) then
      nf=nf+1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,1)=mp-1
      mp=mp+1
      coord(mp,1)=i-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,2)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j+1-jcm 
      coord(mp,3)=k-kcm
      face(nf,3)=mp-1
      mp=mp+1
      coord(mp,1)=i+1-icm 
      coord(mp,2)=j-jcm 
      coord(mp,3)=k-kcm
      face(nf,4)=mp-1
      end if
1000  continue

c now output the VRML file (it is simple)

      filename(1:12)=filestat
      write(filename(18:23),'(I6)') 100000+ibad
      filename(13:18)='-pore-'
      filename(24:27)='.wrl'

c  I want file names involving pore as particle number
c 17 is the unit for the VRML file to be created
         open(unit=17,file=filename)

c write out header of VRML file

         write(17,11)
11       format('#VRML V2.0 utf8')
c         write(17,12)
c12       format('NavigationInfo {')
c         write(17,13)
c13       format('type ["EXAMINE","WALK","FLY"]')
c         write(17,14)
c14       format('}')
         write(17,17)
17       format('Shape {')
         write(17,18)
18       format('geometry IndexedFaceSet {') 
         write(17,19)
19       format('coord Coordinate{point[')
         do 100 m=1,mp
      if(m.lt.mp) write(17,20) -sx*coord(m,1),sy*coord(m,2),
     &        sz*coord(m,3)
      if(m.eq.mp) write(17,21) -sx*coord(m,1),sx*coord(m,2),
     &        sz*coord(m,3)
20       format(3f12.6,',')
21       format(3f12.6)
100      continue
         write(17,22)
22       format(']}')
         write(17,23)
23       format('coordIndex[')
         do 200 n=1,nf
      if(n.lt.nf) write(17,24) face(n,1),face(n,2),face(n,3),face(n,4)
      if(n.eq.nf) write(17,25) face(n,1),face(n,2),face(n,3),face(n,4)
24       format(4i8,' -1,')
25       format(4i8,']}')
200      continue
         write(17,26)
26       format('appearance Appearance{')
         write(17,27)
27       format('material Material{')
         write(17,28)
28       format('diffuseColor 0.923 0.923 0.96 }}}')
         zzk=4.*sz*(kright+1-kcm)
         write(17,29) 0.,0.,zzk
29       format('Viewpoint{position ',3f8.2,'}')        
         close(17)

         return
         end
         subroutine nonSHout(isave,save,ibad,nnsave,iii,filestat,sx)
         integer*2 save(nnsave,iii)
         character*33 filename
         character*12 filestat

      filename(1:12)=filestat
      write(filename(24:29),'(I6)') 100000+ibad
      filename(13:24)='-nonSH-part-'
      filename(30:33)='.dat'

c 41 is the unit for the save non-SH particle to be created
         open(unit=41,file=filename)
         write(41,1) isave,sx
1        format(i10,f9.6)
         do 100 n=1,isave
         write(41,2) save(n,1),save(n,2),save(n,3),save(n,4)
2        format(4i10)
100      continue
         close(41)
         return
         end
         subroutine poreout(isave,save,ipore,nnsave,iii,filestat)
         integer*2 save(nnsave,iii)
         integer*4 ipore
         character*27 filename
         character*12 filestat

      filename(1:12)=filestat
      write(filename(18:23),'(I6)') 100000+ipore
      filename(13:18)='-pore-'
      filename(24:27)='.dat'

c 41 is the unit for the pore to be created
         open(unit=41,file=filename)
         write(41,1) isave
1        format(i10)
         do 100 n=1,isave
         write(41,2) save(n,1),save(n,2),save(n,3),save(n,4)
2        format(4i10)
100      continue
         close(41)
         return
         end
      subroutine pixbin(pix,nx,ny,nz,iscan,in,jn,kn,nphase)

c this version is designed to be incorporated into my
c particle analysis code, the most current version is pp.f
c input is a gray scale microstructure, read into pix(i,j,k)
c program changes this into a binary version, also
c loaded into pix(i,j,k)

c uses Otsu to segment a powder/epoxy image, 4 phases
c A limited watershed split is done by pp.f in 3D, with LWS=5

      real*8 aaa1,aaa0,aaa2,xnx,yny,znz
c f(i) = number fraction of voxels that have gray scale value= i
      real*8 pi,gray(1024),f(0:255)
      integer*2 da1(1024,1024,1024),pix(1024,1024,1024)
      integer*4 nx,ny,nz
      integer*2 in(6),jn(6),kn(6)
      integer*2 old(30000000,4),new(30000000,4)
      integer*2 tempburn
      integer*4 inew,iold

c USER

      pi=4.0*atan(1.0)
      xnx=nx
      yny=ny
      znz=nz

c iscan = 1 for interior scan, which requires removal of particles
c that touch the virtual cylindrical surface, and iscan = 0 for 
c exterior scans of the whole sample that do not require this removal
c iscan is the last number on each line of the *sysconfig.dat file

c initialize da1(i,j,k)
      do 224 k=1,1024
      do 224 j=1,1024
      do 224 i=1,1024
c      if(iscan.eq.1) da1(i,j,k)=3
c      if(iscan.eq.0) da1(i,j,k)=0
      da1(i,j,k)=0
224   continue

c initialize da1 to pix
      do 324 k=1,nz
      do 324 j=1,ny
      do 324 i=1,nx
      da1(i,j,k)=pix(i,j,k)
324   continue

c now do Otsu thresholding on da1
c this version does 3-phase, since there is black outside the circle
c (interior scan on Versa), gray epoxy, and white particles
c even on exterior scans, there are still air bubbles that are
c darker than the epoxy, and there is air outside the sample.

c First make each slice have the same average gray scale
c this corrects for distortions in the analysis
c seems to make a difference in the Otsu analysis

      do 7550 k=1,nz
      gray(k)=0.0
      do 7551 j=1,ny
      do 7551 i=1,nx
      gray(k)=gray(k)+da1(i,j,k)
7551  continue
      gray(k)=gray(k)/float(nx)/float(ny)
      write(8,*) 'average gray scale of ',k,gray(k)
7550  continue
      gray500=gray(nz/2)
      call flush(8)

      do 7560 k=1,nz
      ratio=gray500/gray(k)
      do 7561 j=1,ny
      do 7561 i=1,nx
      da1(i,j,k)=int(ratio*da1(i,j,k)+0.5)
      if(da1(i,j,k).gt.255) da1(i,j,k)=255
7561  continue
7560  continue

c initialize f(i)
      do 4549 m=0,255
      f(m)=0.0
4549  continue

c generate f(i) for gray scale image
      frac=1.0/float(nx*ny*nz)
      do 4550 k=1,nz
      do 4550 j=1,ny
      do 4550 i=1,nx
      ijk=da1(i,j,k)
      f(ijk)=f(ijk)+frac
4550  continue

      write(8,*) ' computed gray scale histogram'
      volfrac=0.0
      do 4555 m=0,255
      volfrac=volfrac+f(m)
      write(8,*) m,f(m),volfrac
4555  continue
      write(8,*) 
      call flush(8)
      
c go over all gray scale values, determine the threshold that maximizes
c the interphase variance = w1*w2*(u1**2-u2**2)
c use 4 phases = outside sample+ pore, straw wall, epoxy, sand
c if interior scan, then use three phases - pore, epoxy, sand
c nphase=3 or 4 is the key

      if(nphase.eq.3) then
      sigmax=0.0
c mt is the trial threshold, m just sums over the gray scales
      do 4542 mt=2,255
      do 4542 ms=2,255
      w1=0.0
      w2=0.0
      w3=0.0
      u1=0.0
      u2=0.0
      u3=0.0
      do 4562 m=0,255
      if(m.le.ms) then
      w1=w1+f(m)
      u1=u1+m*f(m)
      end if
      if(m.gt.ms.and.m.le.mt) then
      w2=w2+f(m)
      u2=u2+m*f(m)
      end if
      if(m.gt.mt) then
      w3=w3+f(m)
      u3=u3+m*f(m)
      end if
4562  continue
      if(w1.eq.0.0) goto 4542
      u1=u1/w1
      if(w2.eq.0.0) goto 4542
      u2=u2/w2
      if(w3.eq.0.0) goto 4542
      u3=u3/w3
      sig=w1*w2*((u1-u2)**2)
      sig=sig+w1*w3*((u1-u3)**2)
      sig=sig+w2*w3*((u2-u3)**2)
c      write(8,*) ms,mt,mu,sig
      if(sigmax.lt.sig) then
      sigmax=sig
      mthmax=mt
      mshmax=ms
      end if
4542  continue
      write(8,*) 
      write(8,*) ' max value of sigma = ',sigmax
      write(8,*) ' Otsu thresholds = ',mshmax,mthmax

c now segment image using mshmax, mthmax, and muhmax
c just do inside boundary like Bernsen does,since I want to compare
c to Bernsen
      aaa0=0.0
      aaa1=0.0
      aaa2=0.0
      do 4602 k=1,nz
      do 4602 j=1,ny
      do 4602 i=1,nx
      if(da1(i,j,k).le.mshmax) then
      da1(i,j,k)=0
      goto 2234
      end if
      if(da1(i,j,k).gt.mshmax.and.da1(i,j,k).le.mthmax) then
      da1(i,j,k)=1
      goto 2234
      end if
      if(da1(i,j,k).gt.mthmax) then
      da1(i,j,k)=4
      goto 2234
      end if
2234  if(da1(i,j,k).eq.0) aaa0=aaa0+1.0 
      if(da1(i,j,k).eq.1) aaa1=aaa1+1.0 
      if(da1(i,j,k).eq.4) aaa3=aaa3+1.0 
4602  continue

      aaa0=aaa0/xnx/yny/znz
      aaa1=aaa1/xnx/yny/znz
      aaa3=aaa3/xnx/yny/znz
      write(8,*) ' outside aaa0 = ',aaa0
      write(8,*) ' epoxy aaa1 = ',aaa1
      write(8,*) ' sand aaa3 = ',aaa3
      call flush(8)

      end if

      if(nphase.eq.4) then
      sigmax=0.0
c mt is the trial threshold, m just sums over the gray scales
      do 4540 mu=2,255
      do 4540 mt=2,255
      do 4540 ms=2,255
      w1=0.0
      w2=0.0
      w3=0.0
      w4=0.0
      u1=0.0
      u2=0.0
      u3=0.0
      u4=0.0
      do 4560 m=0,255
      if(m.le.ms) then
      w1=w1+f(m)
      u1=u1+m*f(m)
      end if
      if(m.gt.ms.and.m.le.mt) then
      w2=w2+f(m)
      u2=u2+m*f(m)
      end if
      if(m.gt.mt.and.m.le.mu) then
      w3=w3+f(m)
      u3=u3+m*f(m)
      end if
      if(m.gt.mu) then
      w4=w4+f(m)
      u4=u4+m*f(m)
      end if
4560  continue
      if(w1.eq.0.0) goto 4540
      u1=u1/w1
      if(w2.eq.0.0) goto 4540
      u2=u2/w2
      if(w3.eq.0.0) goto 4540
      u3=u3/w3
      if(w4.eq.0.0) goto 4540
      u4=u4/w4
      sig=w1*w2*((u1-u2)**2)
      sig=sig+w1*w3*((u1-u3)**2)
      sig=sig+w1*w4*((u1-u4)**2)
      sig=sig+w2*w3*((u2-u3)**2)
      sig=sig+w2*w4*((u2-u4)**2)
      sig=sig+w3*w4*((u3-u4)**2)
c      write(8,*) ms,mt,mu,sig
      if(sigmax.lt.sig) then
      sigmax=sig
      mthmax=mt
      mshmax=ms
      muhmax=mu
      end if
4540  continue
      write(8,*) 
      write(8,*) ' max value of sigma = ',sigmax
      write(8,*) ' Otsu thresholds = ',mshmax,mthmax,muhmax

c now segment image using mshmax, mthmax, and muhmax
c just do inside boundary like Bernsen does,since I want to compare
c to Bernsen
      aaa0=0.0
      aaa1=0.0
      aaa2=0.0
      aaa3=0.0
      do 4600 k=1,nz
      do 4600 j=1,ny
      do 4600 i=1,nx
      if(da1(i,j,k).le.mshmax) then
      da1(i,j,k)=0
      goto 2233
      end if
      if(da1(i,j,k).gt.mshmax.and.da1(i,j,k).le.mthmax) then
      da1(i,j,k)=1
      goto 2233
      end if
      if(da1(i,j,k).gt.mthmax.and.da1(i,j,k).le.muhmax) then
      da1(i,j,k)=2
      goto 2233
      end if
      if(da1(i,j,k).gt.muhmax) then
      da1(i,j,k)=4
      goto 2233
      end if
2233  if(da1(i,j,k).eq.0) aaa0=aaa0+1.0 
      if(da1(i,j,k).eq.1) aaa1=aaa1+1.0 
      if(da1(i,j,k).eq.2) aaa2=aaa2+1.0 
      if(da1(i,j,k).eq.4) aaa3=aaa3+1.0 
4600  continue

      aaa0=aaa0/xnx/yny/znz
      aaa1=aaa1/xnx/yny/znz
      aaa2=aaa2/xnx/yny/znz
      aaa3=aaa3/xnx/yny/znz
      write(8,*) ' outside aaa0 = ',aaa0
      write(8,*) ' straw aaa1 = ',aaa1
      write(8,*) ' epoxy aaa2 = ',aaa2
      write(8,*) ' sand aaa3 = ',aaa3
      call flush(8)
      end if

c for interior scans, do the following
c for exterior scans, just skip
      if(iscan.eq.0) goto 5544
c now get rid of da1(i,j,k)=4 (white) particles that touch da1(i,j,k)=0
c background (black)
c use a version of subroutine small to do so

c principle of subroutine:  
c Turn all outer planes into 0, then 
c build a "particle" or fire inward from each side. When a real particle
c is touched, burn this as well, but don't burn through epoxy matrix.
c When "particle" or fire is stopped,
c then just need to turn all burned voxels into black pores

      tempburn=-1500
c  find first matrix voxels around faces of particle system
      do 990 k=1,nz
      do 990 j=1,ny
      i=1
      da1(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      i=nx
      da1(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
990   continue
      do 991 k=1,nz
      do 991 i=1,nx
      j=1
      da1(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
      j=ny
      da1(i,j,k)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      old(iold,3)=k
991   continue

c  Now start building up new voxels from old set of voxels,
c  thus propagating the fire inward and identifying new particle 
c  voxels connnected to the ones already burned.

60    inew=0
      do 100 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
      kk=old(ijk,3)
c  check all six nearest neighbors of previously burned voxel
      do 90 n=1,6
      i1=ii+in(n)
      j1=jj+jn(n)
      k1=kk+kn(n)

c  Hard boundary conditions (can't pass beyond boundary)

      if(k1.lt.1.or.k1.gt.nz) goto 90
      if(j1.lt.1.or.j1.gt.ny) goto 90
      if(i1.lt.1.or.i1.gt.nx) goto 90

c  Store (i,j,k) labels of newly burned voxels in array new().
c if two kinds of particles, both brighter than epoxy,
c then use 4 phases and treat 4 and 2 as particles
      if(da1(i1,j1,k1).eq.4.or.da1(i1,j1,k1).eq.2.or.
     &      da1(i1,j1,k1).eq.0) then
c      if(da1(i1,j1,k1).eq.4.or.
c     &      da1(i1,j1,k1).eq.0) then
      da1(i1,j1,k1)=tempburn
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      new(inew,3)=k1
      end if
90    continue
100   continue

c  If new voxels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 150 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
      old(ijk,3)=new(ijk,3)
150   continue
      goto 60
      end if

6060    continue

c done burning
c now set all burned voxels to zero, which will include the particles
c touching the outside of the virtual cylinder
       do 225 k=1,nz
       do 225 j=1,ny
       do 225 i=1,nx
       if(da1(i,j,k).eq.tempburn) then
       da1(i,j,k)=0
       end if
225    continue

c exterior scans skip to here
5544   continue

c reset all 1s to 0s, 2s to 0, 3s to 0, and all
c 4s to 1s to make a binary image
c of white (1) particles and black (0) background
c if two kinds of particles, 2 and 4 after segmention,
c then set both to 1 = particle

      do 2200 k=1,1024
      do 2200 j=1,1024
      do 2200 i=1,1024
      if(da1(i,j,k).eq.1) da1(i,j,k)=0
c      if(da1(i,j,k).eq.2) da1(i,j,k)=0
      if(da1(i,j,k).eq.2) da1(i,j,k)=1
      if(da1(i,j,k).eq.3) da1(i,j,k)=0
      if(da1(i,j,k).eq.4) da1(i,j,k)=1
2200  continue

c transfer new binary microstructure back into pix
      do 2300 k=1,nz
      do 2300 j=ny,1,-1
      do 2300 i=1,nx
      pix(i,j,k)=da1(i,j,k)
2300  continue

      return
      end