program region_prb
!
!*******************************************************************************
!
!! REGION_PRB reads voxel data, processes it, and writes it back out.
!
!
!  Modified:
!
!    23 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: maxlist = 8000
  integer, parameter :: max_region = 20
  integer, parameter :: nx = 64
  integer, parameter :: ny = 64
  integer, parameter :: nz = 26
!
  real ave_pos
  integer c(3)
  integer center(4,max_region)
  character ( len = 80 ) filename
  integer i
  integer iregion
  integer ivoxel(nx,ny,nz)
  integer j
  integer k
  integer l
  integer list(maxlist)
  integer nelements
  integer nlist
  integer nregion
  integer num_pos
  real rvoxel(nx,ny,nz)
  integer sum
  integer thresh
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'REGION_PRB'
  write ( *, '(a)' ) '  Set up regions in MRI data.'
!
!  Read the MRI data.
!
  filename = 'roi.ascii'

  call ivoxel_read ( nx, ny, nz, ivoxel, filename )
!
!  Make an OBJ file of the data.
!
  filename = 'roi.obj'

  call ivoxel_to_obj ( nx, ny, nz, ivoxel, filename )
!
!  Find the average nonzero value.
!
  call ivoxel_sum ( nx, ny, nz, ivoxel, sum )

  call ivoxel_count_positive ( nx, ny, nz, ivoxel, num_pos )

  ave_pos = real ( sum ) / real ( num_pos )

! call ivoxel_plot ( nx, ny, nz, ivoxel )
!
!  Zero out all voxels below a given threshhold value.
!
  thresh = nint ( 0.75E+00 * ave_pos )

  call ivoxel_thresh ( nx, ny, nz, ivoxel, thresh )
!
!  Thicken the voxels.
!
  call ivoxel_thicken ( nx, ny, nz, ivoxel )
!
!  Now find the regions.
!
  call ivoxel_to_region ( nx, ny, nz, ivoxel, list, maxlist, nlist, nregion )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6,a)' ) 'The voxels are grouped into ', nregion, ' regions.'
  write ( *, '(a)' ) ' '
!     write ( *, '(a)' ) 'The nonzero array elements are:'
!     write ( *, '(a)' ) ' '

!     do i = 1, nx
!       do j = 1, ny
!         do k = 1, nz
!           l = ivoxel(i,j,k)
!           if ( l /= 0 ) then
!             write ( *, * ) i, j, k, l
!           end if
!         end do
!       end do
!     end do

  if ( nlist > maxlist ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'The stack-based list of regions is unusable.'

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'The stack-based list of regions is:'
    write ( *, '(a)' ) ' '

    iregion = nregion + 1

10      continue

    iregion = iregion - 1

    if ( nlist > 0 ) then

      nelements = list(nlist)
      nlist = nlist - 1

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Region ', iregion, ' includes ', nelements, ' voxels:'

      do l = 1, nelements

        k = list(nlist)
        nlist = nlist - 1
        j = list(nlist)
        nlist = nlist - 1
        i = list(nlist)
        nlist = nlist - 1
        write ( *, * ) i, j, k

      end do

      go to 10

    end if

  end if
!
!  Make a typewriter plot of the regions in Z slices.
!
!     call ivoxel_plot2 ( nx, ny, nz, ivoxel )
!
!  Compute center of mass of each region.
!
  call region_center ( nx, ny, nz, ivoxel, max_region, nregion, center )
!
!  Save the center of mass of the central region.
!
  iregion = 4
  do i = 1, 3
    c(i) = center(i,iregion)
  end do
!
!  Zero out region #4, slide other region numbers down 1.
!
  iregion = 4
  call region_blank ( nx, ny, nz, ivoxel, center, iregion, max_region, nregion )
!
!  Write out an ASCII MRI file containing the voxels, marked by region only.
!
  filename = 'roi_region.ascii'

  call ivoxel_write ( nx, ny, nz, ivoxel, filename )
!
!  Now build a new array, RVOXEL.
!
!  For each region IREGION, construct the line from the center of
!  the central region through its own center.
!
!  For each voxel in region IREGION, move outward along this line
!  until you reach the boundary.  Add PERCENT to each voxel in ARRAY
!  through which the voxel passes, where PERCENT is the percentage of
!  the total region volume represented by one voxel.  And if this
!  is the first time any voxel has passed through this voxel,
!  add 100 * IREGION to it.
!
  call transport ( nx, ny, nz, rvoxel, c, center, ivoxel, max_region )
!
!  Now copy the real ARRAY into the integer IVOXEL.
!
  call rvoxel_to_ivoxel ( nx, ny, nz, rvoxel, ivoxel )
!
!  Do a 0/nonzero plot of the data.
!
!     call ivoxel_plot3 ( nx, ny, nz, ivoxel )
!
!  Write data to a file.
!
  filename = 'roi_colors.ascii'

  call ivoxel_write ( nx, ny, nz, ivoxel, filename )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'REGION_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end