subroutine face_print ( i, j, k, nx, ny, nz, face, iunit )
!
!*******************************************************************************
!
!! FACE_PRINT prints the nodes that define one face of a voxel.
!
!
!  Discussion:
!
!    FACE_PRINT prints a face using the format appropriate for an OBJ
!    file, and is, in fact, a utility routine for IVOXEL_TO_OBJ.
!
!  Modified:
!
!    23 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, K, the indices of a voxel.
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, character ( len = 2 ) FACE, specifies the voxel face.
!    '-I' is the lower I=constant side, '+I' is the upper I=constant side;
!    '-J', '+J', '-K' and '+K' are also valid.
!
!    Input, integer IUNIT, the output unit number.
!
  implicit none
!
  character ( len = 2 ) face
  integer i
  integer iunit
  integer j
  integer k
  integer n000
  integer n001
  integer n010
  integer n011
  integer n100
  integer n101
  integer n110
  integer n111
  integer nx
  integer ny
  integer nz
!
  call voxel_nodes ( i, j, k, nx, ny, nz, n000, n001, n010, n011, n100, &
    n101, n110, n111 )

  if ( face == '+I' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n100, n110, n111, n101
  else if ( face == '-I' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n000, n001, n011, n010
  else if ( face == '+J' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n010, n011, n111, n110
  else if ( face == '-J' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n000, n100, n101, n001
  else if ( face == '+K' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n001, n101, n111, n011
  else if ( face == '-K' ) then
    write ( iunit, '(a,2x,i7,2x,i7,2x,i7,2x,i7)' ) 'f', n000, n010, n110, n100
  end if

  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine ivoxel_bound_print ( nx, ny, nz, ivoxel, iunit, num_face )
!
!*******************************************************************************
!
!! IVOXEL_BOUND_PRINT writes bounding faces to an OBJ file.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, integer IUNIT, the output unit number.
!
!    Output, integer NUM_FACE, the number of faces printed.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer iunit
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  integer num_face
!
  num_face = 0
!
!  I-planes
!
  do k = 1, nz
    do j = 1, ny
      do i = 0, nx

        if ( i == 0 ) then

          if ( ivoxel(i+1,j,k) /= 0 ) then
            call face_print ( i+1, j, k, nx, ny, nz, '-I', iunit )
            num_face = num_face + 1
          end if

        else if ( i < nx ) then

          if ( ivoxel(i,j,k) == 0 .and. ivoxel(i+1,j,k) /= 0 ) then
            call face_print ( i+1, j, k, nx, ny, nz, '-I', iunit )
            num_face = num_face + 1
          else if ( ivoxel(i,j,k) /= 0 .and. ivoxel(i+1,j,k) == 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+I', iunit )
            num_face = num_face + 1
          end if

        else if ( i == nx ) then

          if ( ivoxel(i,j,k) /= 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+I', iunit )
            num_face = num_face + 1
          end if

        end if

      end do
    end do
  end do
!
!  J-planes
!
  do k = 1, nz
    do i = 1, nx
      do j = 0, ny

        if ( j == 0 ) then

          if ( ivoxel(i,j+1,k) /= 0 ) then
            call face_print ( i, j+1, k, nx, ny, nz, '-J', iunit )
            num_face = num_face + 1
          end if

        else if ( j < ny ) then

          if ( ivoxel(i,j,k) == 0 .and. ivoxel(i,j+1,k) /= 0 ) then
            call face_print ( i, j+1, k, nx, ny, nz, '-J', iunit )
            num_face = num_face + 1
          else if ( ivoxel(i,j,k) /= 0 .and. ivoxel(i,j+1,k) == 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+J', iunit )
            num_face = num_face + 1
          end if

        else if ( j == ny ) then

          if ( ivoxel(i,j,k) /= 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+J', iunit )
            num_face = num_face + 1
          end if

        end if

      end do
    end do
  end do
!
!  K-planes
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        if ( k == 0 ) then

          if ( ivoxel(i,j,k+1) /= 0 ) then
            call face_print ( i, j, k+1, nx, ny, nz, '-K', iunit )
            num_face = num_face + 1
          end if

        else if ( k < nz ) then

          if ( ivoxel(i,j,k) == 0 .and. ivoxel(i,j,k+1) /= 0 ) then
            call face_print ( i, j, k+1, nx, ny, nz, '-K', iunit )
            num_face = num_face + 1
          else if ( ivoxel(i,j,k) /= 0 .and. ivoxel(i,j,k+1) == 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+K', iunit )
            num_face = num_face + 1
          end if

        else if ( k == nz ) then

          if ( ivoxel(i,j,k) /= 0 ) then
            call face_print ( i, j, k, nx, ny, nz, '+K', iunit )
            num_face = num_face + 1
          end if

        end if

      end do
    end do
  end do

  return
end
subroutine ivoxel_count_positive ( nx, ny, nz, ivoxel, num_pos )
!
!*******************************************************************************
!
!! IVOXEL_COUNT_POSITIVE counts the positive entries in a voxel array.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Output, integer NUM_POS, the number of positive entries in the voxel array.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer j
  integer k
  integer num_pos
  integer ivoxel(nx,ny,nz)
!
  num_pos = 0

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

        if ( ivoxel(i,j,k) > 0 ) then
          num_pos = num_pos + 1
        end if

      end do
    end do
  end do

  return
end
subroutine ivoxel_plot ( nx, ny, nz, ivoxel )
!
!*******************************************************************************
!
!! IVOXEL_PLOT prints out a typewriter plot of the Z slices.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  character ( len = 64 ) string
!
  do k = 1, nz

    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'Z plane ', k

    do j = 1, ny

      do i = 1, nx

        if ( ivoxel(i,j,k) == 0 ) then
          string(i:i)=' '
        else if ( ivoxel(i,j,k) == 158 ) then
          string(i:i) = '.'
        else
          string(i:i) = '*'
        end if

      end do

      write ( *, '(a)' ) trim ( string )

    end do

  end do

  return
end
subroutine ivoxel_plot2 ( nx, ny, nz, ivoxel )
!
!*******************************************************************************
!
!! IVOXEL_PLOT2 prints out a typewriter plot of the regions.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  character ( len = 64 ) string
!
  do k = 1, nz

    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'Z plane ', k

    do j = 1, ny

      string = ' '

      do i = 1, nx

        if ( ivoxel(i,j,k) /= 0 ) then
          write ( string(i:i), '(i1)' ) ivoxel(i,j,k)
        end if

      end do

      write ( *, '(a)' ) trim ( string )

    end do

  end do

  return
end
subroutine ivoxel_plot3 ( nx, ny, nz, ivoxel )
!
!*******************************************************************************
!
!! IVOXEL_PLOT3 prints out a typewriter plot of data/100.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  character ( len = 64 ) string
!
  do k = 1, nz

    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'Z plane ', k

    do j = 1, ny

      string = ' '

      do i = 1, nx

        if ( ivoxel(i,j,k) /= 0 ) then
          write ( string(i:i), '(i1)' ) ivoxel(i,j,k) / 100
        end if

      end do

      write ( *, '(a)' ) trim ( string )

    end do

  end do

  return
end
subroutine ivoxel_read ( nx, ny, nz, ivoxel, filename )
!
!*******************************************************************************
!
!! IVOXEL_READ reads MRI data from an ASCII file, one item per line.
!
!
!  Modified:
!
!    03 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Output, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, character ( len = * ) FILENAME, the name of the file to be read.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  character ( len = 80 ) filename
  integer i
  integer ios
  integer iunit
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  integer nrec
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = filename, status = 'old', iostat = ios, &
    form = 'formatted', access = 'sequential' )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IVOXEL_READ - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file.'
    stop
  end if
!
!  Read the data.
!
  nrec = 0

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

        read ( iunit, *, iostat = ios ) ivoxel(i,j,k)

        if ( ios /= 0 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'IVOXEL_READ - Fatal error!'
          write ( *, '(a,i6)' ) '  END or ERR reading record ', nrec+1
          stop
        end if

        nrec = nrec + 1

      end do
    end do
  end do
!
!  Close the file.
!
  close ( unit = iunit )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6,a)' ) 'IVOXEL_READ read ', nrec, ' records.'

  return
end
subroutine ivoxel_sum ( nx, ny, nz, ivoxel, vsum )
!
!*******************************************************************************
!
!! IVOXEL_SUM sums the entries in a voxel array.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Output, integer VSUM, the sum of the entries in the voxel array.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer ivoxel(nx,ny,nz)
  integer vsum
!
  vsum = sum ( ivoxel(1:nx,1:ny,1:nz) )

  return
end
subroutine ivoxel_thicken ( nx, ny, nz, ivoxel )
!
!*******************************************************************************
!
!! IVOXEL_THICKEN "thickens" the voxels.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer i2
  integer ivoxel(nx,ny,nz)
  integer ihi
  integer ilo
  integer j
  integer j2
  integer jhi
  integer jlo
  integer k
  integer k2
  integer khi
  integer klo
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        if ( ivoxel(i,j,k) > 0 ) then

          ilo = max ( i-1, 1 )
          ihi = min ( i+1, nx )
          jlo = max ( j-1, 1 )
          jhi = min ( j+1, nx )
          klo = max ( k-1, 1 )
          khi = min ( k+1, nx )

          do i2 = ilo, ihi
            do j2 = jlo, jhi
              do k2 = klo, khi
                if ( ivoxel(i2,j2,k2) == 0 ) then
                  ivoxel(i2,j2,k2) = - ivoxel(i,j,k)
                end if
              end do
            end do
          end do

        end if

      end do
    end do
  end do
!
!  The new voxels were marked with negative values, so that we
!  didn't allow yet more voxels to grow from THEM.
!  Now that we're done, we can reset the values in these voxels
!  to positive values, so they join their region.
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        if ( ivoxel(i,j,k) < 0 ) then
          ivoxel(i,j,k) = - ivoxel(i,j,k)
        end if

      end do
    end do
  end do

  return
end
subroutine ivoxel_thresh ( nx, ny, nz, ivoxel, thresh )
!
!*******************************************************************************
!
!! IVOXEL_THRESH zeroes out array entries below a given threshhold.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input/output, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, integer THRESH, the threshhold value.  Any voxel with a value
!    less than THRESH is reset to 0.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer i
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  integer nkeep
  integer nthresh
  integer thresh
!
  nkeep = 0
  nthresh = 0

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

        if ( 0 < ivoxel(i,j,k) ) then

          if ( ivoxel(i,j,k) < thresh ) then

            nthresh = nthresh + 1
            ivoxel(i,j,k) = 0

          else

            nkeep = nkeep + 1

          end if

        end if

      end do
    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'IVOXEL_THRESH' 
  write ( *, '(a,i6)' ) '  Elements kept             ', nkeep
  write ( *, '(a,i6)' ) '  Elements zeroed out       ', nthresh
  write ( *, '(a,i6)' ) '  which were below THRESH = ', thresh

  return
end
subroutine ivoxel_to_obj ( nx, ny, nz, ivoxel, file_name )
!
!*******************************************************************************
!
!! IVOXEL_TO_OBJ writes out an OBJ file from a voxel array.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to be created.
!    An extension of '.obj' is recommended.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  character ( len = * ) file_name
  integer ios
  integer iunit
  integer ivoxel(nx,ny,nz)
  integer num_face
  integer num_node
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'replace', &
    form = 'formatted', iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IVOXEL_TO_OBJ - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file.'
    stop
  end if
!
!  Print the header.
!
  write ( iunit, '(a)' ) '# shape.obj created by john'
!
!  Print the node position data.
!
  write ( iunit, '(a)' ) ' '
  call node_print ( nx, ny, nz, iunit )
!
!  Print the face data.
!
  write ( iunit, '(a)' ) ' '

  call ivoxel_bound_print ( nx, ny, nz, ivoxel, iunit, num_face )
!
!  Close the file.
!
  close ( unit = iunit )
!
!  Report
!
  num_node = ( nx + 1 ) * ( ny + 1 ) * ( nz + 1 )
     
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'IVOXEL_TO_OBJ'
  write ( *, '(a)' ) '  OBJ file created.'
  write ( *, '(a,i6)' ) '  Number of nodes: ', num_node
  write ( *, '(a,i6)' ) '  Number of faces: ', num_face

  return
end
subroutine ivoxel_to_region ( nx, ny, nz, ivoxel, list, maxlist, nlist, &
  nregion )
!
!*******************************************************************************
!
!! IVOXEL_TO_REGION arranges a set of voxels into contiguous regions.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the number of voxels in the X, Y and
!    Z directions.
!
!    Input/output, integer IVOXEL(NX,NY,NZ).
!
!    On input, IVOXEL(I,J,K) has the values:
!      0, if the voxel is OFF;
!      anything else, if the voxel is ON.
!
!    On output, IVOXEL(I,J,K) has the values:
!      0, if the voxel is off,
!      N, if the voxel is ON, and is part of region N.
!
!    Output, integer LIST(MAXLIST), contains, in stack form, a list
!    of the indices of the elements in each region.
!
!    The number of elements in NREGION is NELEM = LIST(NLIST).  The 
!    (I,J,K) indices of the last element in this region are in 
!    LIST(NLIST-3) through LIST(NLIST-1), and the first element is 
!    listed in LIST(NLIST-3*NELEM), LIST(NLIST-3*NELEM+1), 
!    LIST(NLIST-3*NELEM+2).
!
!    The number of elements in NREGION-1 is listed in 
!    LIST(NLIST-3*NELEM-1), and so on.
!
!    Input, integer MAXLIST, the maximum length of the array used to
!    list the elements of the regions.
!
!    Output, integer NLIST, the number of entries of LIST that were used.
!    However, if NLIST > MAXLIST, then there was not enough space in
!    LIST to store the data properly, and LIST should not be used,
!    although the data in IVOXEL should be correct.
!
!    Output, integer NREGION, the number of regions discovered.
!
  implicit none
!
  integer, parameter :: maxstack = 5000
  integer maxlist
  integer nx
  integer ny
  integer nz
!
  integer i
  integer i2
  integer ibase
  integer ihi
  integer ilo
  integer ivoxel(nx,ny,nz)
  integer j
  integer j2
  integer jbase
  integer jhi
  integer jlo
  integer k
  integer k2
  integer kbase
  integer khi
  integer klo
  integer list(maxlist)
  integer nabes
  integer ncan
  integer nelements
  integer nlist
  integer nregion
  integer nstack
  integer stack(maxstack)
!
!  Reset all nonzero entries of IVOXEL to -1.
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        if ( ivoxel(i,j,k) /= 0 ) then
          ivoxel(i,j,k) = -1
        end if

      end do
    end do
  end do
!
!  Start the number of items in the region list at 0.
!
  nlist = 0
!
!  Start the number of regions at 0.
!
  nregion = 0
!
!  The stack begins empty.
!
  nstack = 0
!
!  Search for an unused "ON" voxel from which we can "grow" a new region.
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz
!
!  We found a voxel that is "ON", and does not belong to any region.
!
        if ( ivoxel(i,j,k) == -1 ) then
!
!  Increase the number of regions.
!
          nregion = nregion + 1
!
!  Add this voxel to the region.
!
          ivoxel(i,j,k) = nregion
!
!  Add this voxel to the stack.
!
          if ( nstack + 4 > maxstack ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'IVOXEL_TO_REGION - Fatal error!'
            write ( *, '(a)' ) '  The internal stack overflowed.'
            write ( *, '(a)' ) '  The algorithm has failed.'
            stop
          end if

          stack(nstack+1) = i
          stack(nstack+2) = j
          stack(nstack+3) = k

          stack(nstack+4) = 1

          nstack = nstack + 4
!
!  Add this voxel to the description of the region.
!
          nelements = 1

          if ( nlist + 3 <= maxlist ) then
            list(nlist+1) = i
            list(nlist+2) = j
            list(nlist+3) = k
          end if

          nlist = nlist + 3

10            continue
!
!  Find all neighbors of BASE that are "ON" but unused.
!  Mark them as belonging to this region, and stack their indices.
!
          ibase = stack(nstack-3)
          jbase = stack(nstack-2)
          kbase = stack(nstack-1)

          ilo = max ( ibase-1, 1 )
          ihi = min ( ibase+1, nx )
          jlo = max ( jbase-1, 1 )
          jhi = min ( jbase+1, ny )
          klo = max ( kbase-1, 1 )
          khi = min ( kbase+1, nz )

          nabes = 0

          do i2 = ilo, ihi
            do j2 = jlo, jhi
              do k2 = klo, khi
!
!  We found a neighbor to our current search point, which is "ON" and unused.
!
                if ( ivoxel(i2,j2,k2) == -1 ) then
!
!  Increase the number of neighbors.
!
                  nabes = nabes + 1
!
!  Mark the neighbor as belonging to the region.
!
                  ivoxel(i2,j2,k2) = nregion
!
!  Add the neighbor to the stack.
!
                  if ( nstack+3 > maxstack ) then
                    write ( *, '(a)' ) ' '
                    write ( *, '(a)' ) 'IVOXEL_TO_REGION - Fatal error!'
                    write ( *, '(a)' ) '  The internal stack overflowed.'
                    write ( *, '(a)' ) '  The algorithm has failed.'
                    stop
                  end if

                  stack(nstack+1) = i2
                  stack(nstack+2) = j2
                  stack(nstack+3) = k2

                  nstack = nstack+3
!
!  Add the neighbor to the description of the region.
!
                  nelements = nelements + 1

                  if ( nlist+3 <= maxlist ) then
                    list(nlist+1) = i2
                    list(nlist+2) = j2
                    list(nlist+3) = k2
                  end if

                  nlist = nlist + 3

                end if

              end do
            end do
          end do
!
!  If any new neighbors were found, take the last one as the basis
!  for a deeper search.
!
          if ( nabes > 0 ) then

            if ( nstack+1 > maxstack ) then
              write ( *, '(a)' ) ' '
              write ( *, '(a)' ) 'IVOXEL_TO_REGION - Fatal error!'
              write ( *, '(a)' ) '  The internal stack overflowed.'
              write ( *, '(a)' ) '  The algorithm has failed.'
              stop
            end if

            stack(nstack+1) = nabes
            nstack = nstack + 1
            go to 10

          end if
!
!  If the current search point had no new neighbors, drop it from the stack.
!
          ncan = stack(nstack) - 1
          nstack = nstack - 3
          stack(nstack) = ncan 
!
!  If there are still any unused candidates at this level, take the
!  last one as the basis for a deeper search.
!
          if ( stack(nstack) > 0 ) then
            go to 10
          end if
!
!  If there are no more unused candidates at this level, then we need
!  to back up a level in the stack.  If there are any candidates at
!  that earlier level, then we can still do more searching.
!
          nstack = nstack - 1

          if ( nstack > 0 ) then
            go to 10
          end if
!
!  If we have exhausted the stack, we have completed this region.
!  Tag the number of elements to the end of the region description list.
!
          nlist = nlist + 1
          if ( nlist <= maxlist ) then
            list(nlist) = nelements
          end if

        end if

      end do
    end do
  end do
!
!  Print some warnings.
!
  if ( nlist > maxlist ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IVOXEL_TO_REGION - Warning!'
    write ( *, '(a)' ) '  MAXLIST was too small to list the regions.'
    write ( *, '(a)' ) '  Do not try to use the LIST array!'
    write ( *, '(a)' ) '  The IVOXEL data is OK, however.'
  end if

  return
end
subroutine ivoxel_write ( nx, ny, nz, ivoxel, filename )
!
!*******************************************************************************
!
!! IVOXEL_WRITE writes MRI data to an ASCII file, one item per line.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, character ( len = * ) FILENAME, the name of the file to be created.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  character ( len = * ) filename
  integer i
  integer ios
  integer iunit
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  integer nrec
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = filename, status = 'replace', iostat = ios, &
    form = 'formatted', access = 'sequential' )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IVOXEL_WRITE - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file.'
    stop
  end if

  nrec = 0

  do k = 1, nz
    do j = 1, ny
      do i = 1, nx
        write ( iunit, * ) ivoxel(i,j,k)
        nrec = nrec + 1
      end do
    end do
  end do

  close ( unit = iunit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'IVOXEL_WRITE'
  write ( *, '(a,i6)' ) '  Number of records written: ', nrec

  return
end
subroutine node_print ( nx, ny, nz, iunit )
!
!*******************************************************************************
!
!! NODE_PRINT prints the nodes that define one face of a voxel.
!
!
!  Discussion:
!
!    NODE_PRINT prints the node information using the format appropriate 
!    for an OBJ file, and is, in fact, a utility routine for IVOXEL_TO_OBJ.
!
!  Modified:
!
!    23 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IUNIT, the output unit number.
!
  implicit none
!
  integer i
  integer iunit
  integer j
  integer k
  integer nx
  integer ny
  integer nz
  real x
  real xmax
  real y
  real ymax
  real z
  real zmax
!
  xmax = real ( nx )
  ymax = real ( ny )
  zmax = real ( nz )

  do k = 0, nz
    z = zmax * real ( k  ) / real ( nz )
    do j = 0, ny
      y = ymax * real ( j ) / real ( ny )
      do i = 0, nx
        x = xmax * real ( i ) / real ( nx )
        write ( iunit, '(a,2x,f8.4,2x,f8.4,2x,f8.4)' ) 'v', x, y, z
      end do
    end do
  end do

  return
end
subroutine region_blank ( nx, ny, nz, ivoxel, center, iregion, max_region, &
  nregion )
!
!*******************************************************************************
!
!! REGION_BLANK zeroes out voxels in a particular numbered region.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer MAX_REGION, the maximum number of regions for which
!    storage has been allocated.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer max_region
  integer nx
  integer ny
  integer nz
!
  integer center(4,max_region)
  integer i
  integer ivoxel(nx,ny,nz)
  integer iregion
  integer j
  integer k
  integer nregion
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        if ( ivoxel(i,j,k) == iregion ) then
          ivoxel(i,j,k) = 0
        else if ( ivoxel(i,j,k) > iregion ) then
          ivoxel(i,j,k) = ivoxel(i,j,k) - 1
        end if

      end do
    end do
  end do
!
!  Shift the CENTER data down one index.
!
  do j = iregion+1, nregion
    do i = 1, 4
      center(i,j-1) = center(i,j)
    end do
  end do
!
!  Update the number of regions.
!
  nregion = nregion - 1

  return
end
subroutine region_center ( nx, ny, nz, ivoxel, max_region, nregion, center )
!
!*******************************************************************************
!
!! REGION_CENTER computes the centers of mass of the regions.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
!    Input, integer MAX_REGION, the maximum number of regions for which
!    storage has been allocated.
!
!    Input, integer NREGION, the number of regions.
!
!    Output, real CENTER(4,MAX_REGION), contains the coordinates of the
!    center of mass of region I in CENTER(1,I), CENTER(2,I), CENTER(3,I),
!    and the number of voxels in region I in CENTER(4,I).
!
  implicit none
!
  integer max_region
  integer nx
  integer ny
  integer nz
!
  real a
  real b
  integer center(4,max_region)
  integer i
  integer ivoxel(nx,ny,nz)
  integer iregion
  integer j
  integer k
  integer nregion
!
  center(1:4,1:max_region) = 0
!
!  Add each (I,J,K) to its region's total.
!
  do i = 1, nx
    do j = 1, ny
      do k = 1, nz

        iregion = ivoxel(i,j,k)

        if ( iregion /= 0 ) then

          if ( 0 < iregion .and. iregion < max_region ) then
            center(1,iregion) = center(1,iregion) + i
            center(2,iregion) = center(2,iregion) + j
            center(3,iregion) = center(3,iregion) + k
            center(4,iregion) = center(4,iregion) + 1
          end if

        end if

      end do
    end do
  end do
!
!  Now normalize each center of mass by the number of voxels in the
!  region.  We round to the nearest integer.
!
  do iregion = 1, nregion
    if ( center(4,iregion) /= 0 ) then
      do i = 1, 3
        a = real ( center(i,iregion) ) 
        b = real ( center(4,iregion) )
        center(i,iregion) = nint ( a / b )
      end do
    end if
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CMASS:'
  write ( *, '(a)' ) '  Region, Center(I,J,K), Voxels:'
  write ( *, '(a)' ) ' '

  do iregion = 1, nregion
    write ( *, '(i6,4i6)' ) iregion, center(1:4,iregion)
  end do

  return
end
subroutine rvoxel_to_ivoxel ( nx, ny, nz, rvoxel, ivoxel )
!
!*******************************************************************************
!
!! RVOXEL_TO_IVOXEL copies real voxel data into an integer voxel array.
!
!
!  Modified:
!
!    21 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, real RVOXEL(NX,NY,NZ), an array of real voxel data.
!
!    Output, integer IVOXEL(NX,NY,NZ), a copy of the real data, rounded
!    using the nearest-integer function NINT.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer ivoxel(nx,ny,nz)
  real rvoxel(nx,ny,nz)
!
  ivoxel(1:nx,1:ny,1:nz) = nint ( rvoxel(1:nx,1:ny,1:nz) )

  return
end
subroutine transport ( nx, ny, nz, rvoxel, c, center, ivoxel, max_region )
!
!*******************************************************************************
!
!! TRANSPORT transports voxels to the boundary, counts intermediate hits.
!
!
!  Discussion:
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Input, integer MAX_REGION, the maximum number of regions for which
!    storage has been allocated.
!
!    Input, integer IVOXEL(NX,NY,NZ), an array of voxel data.
!
  implicit none
!
  integer max_region
  integer nx
  integer ny
  integer nz
!
  integer a(3)
  integer b(3)
  integer c(3)
  integer center(4,max_region)
  integer i
  integer ivoxel(nx,ny,nz)
  integer iregion
  integer j
  integer k
  real percent
  real rvoxel(nx,ny,nz)
!
  rvoxel(1:nx,1:ny,1:nz) = 0.0E+00

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

        if ( ivoxel(i,j,k) /= 0 ) then

          iregion = ivoxel(i,j,k)
          a(1) = i
          a(2) = j
          a(3) = k
          b(1) = center(1,iregion)
          b(2) = center(2,iregion)
          b(3) = center(3,iregion)
          percent = 100.0E+00 / real ( center(4,iregion) )

          call transvox ( nx, ny, nz, a, rvoxel, b, c, iregion, percent )

        end if

      end do
    end do
  end do

  return
end
subroutine transvox ( nx, ny, nz, a, rvoxel, b, c, iregion, percent )
!
!*******************************************************************************
!
!! TRANSVOX transports one voxel to the boundary.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
  implicit none
!
  integer nx
  integer ny
  integer nz
!
  integer a(3)
  integer b(3)
  integer c(3)
  integer i
  integer i2
  integer inc
  integer iregion
  integer j
  integer j2
  integer jnc
  integer k
  integer k2
  integer knc
  real percent
  real rvoxel(nx,ny,nz)
!
  i = a(1)
  j = a(2)
  k = a(3)
!
!  Mark the image of the voxel.
!
  if ( rvoxel(i,j,k) == 0.0E+00 ) then
    rvoxel(i,j,k) = 100.0E+00 * real ( iregion )
  end if

  rvoxel(i,j,k) = rvoxel(i,j,k) + percent

  inc = b(1) - c(1)
  jnc = b(2) - c(2)
  knc = b(3) - c(3)
!
!  Starting at (I,J,K), transport the voxel out to the boundary.
!
  i2 = i
  j2 = j
  k2 = k

  do

    call voxel_index_step ( i, j, k, i2, j2, k2, inc, jnc, knc )

    if ( i2 < 1 .or. nx < i2 .or. &
         j2 < 1 .or. ny < j2 .or. &
         k2 < 1 .or. nz < k2 ) then
      exit
    end if

    if ( rvoxel(i2,j2,k2) == 0.0E+00 ) then
      rvoxel(i2,j2,k2) = 100.0E+00 * real ( iregion )
    end if

    rvoxel(i2,j2,k2) = rvoxel(i2,j2,k2) + percent

  end do

  return
end
subroutine voxel_index_step ( i1, j1, k1, i2, j2, k2, inc, jnc, knc )
!
!*******************************************************************************
!
!! VOXEL_INDEX_STEP computes indices of voxels along a line from a given point.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I1, J1, K1, the coordinates of the base voxel from
!    which the line begins.
!
!    Input/output, integer I2, J2, K2.
!
!    On input, these are the coordinates of the current voxel on 
!    the line.  For the first call, these might be I1, J1 and K1.
!
!    On output, these are the coordinates of the next voxel along 
!    the line.
!
!    Input, integer INC, JNC, KNC, the increments to the voxels.
!    These values define the direction along which the line proceeds.
!    However, the voxels on the line will typically be incremented
!    by a fractional value of the vector (INC,JNC,KNC), and the
!    result is essentially rounded.
!
!    If you input INC = JNC = KNC, then no movement is possible,
!    and none is made.
!
  implicit none
!
  real alpha
  real alphai
  real alphaj
  real alphak
  real, parameter :: big = 100000.0E+00
  integer i1
  integer i2
  integer inc
  integer j1
  integer j2
  integer jnc
  integer k1
  integer k2
  integer knc
!
!  Assuming for the moment that (I,J,K) can take on real values,
!  points on the line have the form:
!
!    I = I1 + alpha * inc
!    J = J1 + alpha * jnc
!    K = K1 + alpha * knc
!
  if ( inc == 0 .and. jnc == 0 .and. knc == 0 ) then
    return
  end if

  alpha = 0.0E+00
!
!  Compute the smallest ALPHA that will change I2, J2 or K2 by +-0.5.
!
  if ( inc > 0 ) then
    alphai = real ( i2 - i1 + 0.5E+00 ) / real ( inc )
  else if ( inc < 0 ) then
    alphai = real ( i2 - i1 - 0.5E+00 ) / real ( inc )
  else
    alphai = big
  end if

  if ( jnc > 0 ) then
    alphaj = real ( j2 - j1 + 0.5E+00 ) / real ( jnc )
  else if ( jnc < 0 ) then
    alphaj = real ( j2 - j1 - 0.5E+00 ) / real ( jnc )
  else
    alphaj = big
  end if

  if ( knc > 0 ) then
    alphak = real ( k2 - k1 + 0.5E+00 ) / real ( knc )
  else if ( knc < 0 ) then
    alphak = real ( k2 - k1 - 0.5E+00 ) / real ( knc )
  else
    alphaj = big
  end if
!
!  The ALPHA of smallest positive magnitude represents the closest 
!  next voxel.
!
  alpha = big

  if ( alphai > 0.0E+00 ) then
    alpha = min ( alpha, alphai )
  end if

  if ( alphaj > 0.0E+00 ) then
    alpha = min ( alpha, alphaj )
  end if

  if ( alphak > 0.0E+00 ) then
    alpha = min ( alpha, alphak )
  end if
!
!  Move to the new voxel.  Whichever index just made the half
!  step must be forced to take a whole step.
!
  if ( alpha == alphai ) then
    i2 = i2 + sign ( 1, inc )
    j2 = j1 + nint ( alpha * jnc )
    k2 = k1 + nint ( alpha * knc )
  else if ( alpha == alphaj ) then
    i2 = i1 + nint ( alpha * inc )
    j2 = j2 + sign ( 1, jnc )
    k2 = k1 + nint ( alpha * knc )
  else if ( alpha == alphak ) then
    i2 = i1 + nint ( alpha * inc )
    j2 = j1 + nint ( alpha * jnc )
    k2 = k2 + sign ( 1, knc )
  end if

  return
end
subroutine voxel_nodes ( i, j, k, nx, ny, nz, n000, n001, n010, &
  n011, n100, n101, n110, n111 )
!
!*******************************************************************************
!
!! VOXEL_NODES returns the indices of the nodes of a voxel.
!
!
!  Diagram:
!
!                           n011-----n111
!                            /|      /|
!                           / |     / |
!   ^                     n001----n101|
!   K                      |  |    |  |
!   |                      |  |    |  |
!   T                      |  |    |  |
!   h                      |  |    |  |
!   i                      |  |    |  |
!   r                      |n010---|-n110
!   d                      | /     | /
!   |                      |/      |/
!   |  I Numbered first->-n000----n100
!   |                     /
!   |              J Numbered
!                   second
!
!  Modified:
!
!    23 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, K, the indices of a voxel.
!
!    Input, integer NX, NY, NZ, the X, Y and Z dimensions of the voxel data.
!
!    Output, integer N000, N001, N010, N011, N100, N101, N110, N111, the
!    indices of the 8 nodes associated with a voxel.
!
  implicit none
!
  integer i
  integer j
  integer k
  integer n000
  integer n001
  integer n010
  integer n011
  integer n100
  integer n101
  integer n110
  integer n111
  integer nx
  integer ny
  integer nz
!
  n000 = ( k - 1 ) * ( ny + 1 ) * ( nx + 1 ) + ( j - 1 ) * ( nx + 1 ) + i

  n100 = n000 + 1
  n010 = n000 + ( nx + 1 )
  n110 = n010 + 1

  n001 = n000 + ( ny + 1 ) * ( nx + 1 )
  n101 = n001 + 1
  n011 = n001 + ( nx + 1 )
  n111 = n011 + 1

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end