program main

!*****************************************************************************80
!
!! MAIN is the main program for REGION_PRB.
!
!  Discussion:
!
!    REGION_PRB reads voxel data, processes it, and writes it back out.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 September 2004
!
!  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 ( kind = 8 ) 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 ( kind = 8 ) dvoxel(nx,ny,nz)
  integer sum
  integer thresh

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'REGION_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  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.75D+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)' ) ' '

  if ( .false. ) then
    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 ( *, '(4i6)' ) i, j, k, l
          end if
        end do
      end do
    end do
  end if

  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

    do

      iregion = iregion - 1

      if ( nlist <= 0 ) then
        exit
      end if

      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 ( *, '(3i6)' ) i, j, k

      end do

    end do

  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 )
!
!  Build a new array, DVOXEL.
!
!  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, dvoxel, c, center, ivoxel, max_region )
!
!  Copy the real ARRAY into the integer IVOXEL.
!
  call dvoxel_to_ivoxel ( nx, ny, nz, dvoxel, ivoxel )
!
!  Do a 0/nonzero plot of the data.
!
  if ( .false. ) then
    call ivoxel_plot3 ( nx, ny, nz, ivoxel )
  end if
!
!  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.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end