[Insight-users] MPI Fast Marching

Fri Oct 23 11:54:50 EDT 2009

On Tue, 2009-10-20 at 12:10 -0400, Luis Ibanez wrote:
> Hi Kevin,
> 
> This is great !
> 
> I hope you are planning on submitting it to the Insight Journal.

Is MPI available in the IJ environment?

It is my desire to submit it to the insight journal.

Realizing my desire during work hours however is up to my boss and
whether or not the lab gets any credit from funding agencies or Ohio
University come grant time.

Sounds like Open Access Week Day #3 policy 1 to me.

> 
> 
> Regarding your point about the float comparison, I agree with
> you in that an exact comparison is unrealistic.  You may want
> to consider the use of the following "epsilon" values:
> 
>               vcl_numeric_limits<double>::epsilon()
>               vcl_numeric_limits<float>::epsilon()
> 
> as tolerances for the comparison.
> 
> 
> Another option is to use the value:
> 
>             itk::NumericTraits< float >::min()
> 

Thank you for these comments.

These numbers are tiny.

I am thinking about the re-setting of the trial and alive points between
iterations.

I can easily see the level plus or minus a tolerance falling completely
between the values at each point.

A 1d ASCII art speed image :

.5 .5 1 1 1 1.5 1.5 2 2

split into two padded images :

.5 .5 1 1 1
          1 1.5 1.5 2 2

fast marching seeded at both ends with value 0 :

0  2  3 4 5
          2.8 1.8 1.2 .5 0

Now both pieces need to recompute their trial and alive points based on
the minimum value of their boundaries after MPI communication. For the
left piece 2.8 plus or minus a tiny value will not equal any value in
the left piece, so I made a neighborhood iterator walk the image and set
the points as; alive if its value and all its neighbor's values are less
than 2.8, trial if its value is less than 2.8 and any of its neighbors
values are greater than 2.8, and nothing was set if the value was above
2.8.

A  T  N N T
          A   A   A   A  A

The next round of fast marching would produce something like :

0  2  3 3.8 2.8
            2.8 1.8 1.2 .5 0

Does this seem like the right way to go about it?

I've attached my most recent code with huge chunks of wrong stuff
commented out.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: itkImageRegionNonCubeSplitter.h
Type: text/x-chdr
Size: 3886 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091023/0710c62f/attachment-0001.h>
-------------- next part --------------

#ifndef _itkImageRegionNonCubeSplitter_txx
#define _itkImageRegionNonCubeSplitter_txx
#include "itkImageRegionNonCubeSplitter.h"
#include <math.h>

namespace itk
{

/**
 *
 */
template <unsigned int VImageDimension>
unsigned int 
ImageRegionNonCubeSplitter<VImageDimension>
::GetNumberOfSplits(const RegionType &region, unsigned int requestedNumber)
{
  const SizeType &regionSize = region.GetSize();

  // if a given region has fewer pixels than requestedNumber, then
  // only split number of pixels times
  unsigned int i, numPixels;
  numPixels = 1;
  for (i=0; i < VImageDimension; ++i)
    numPixels *= regionSize[i];
  if (numPixels > requestedNumber)
    return requestedNumber;

  return numPixels;
}

/**
 *
 */
template <unsigned int VImageDimension>
ImageRegion<VImageDimension>
ImageRegionNonCubeSplitter<VImageDimension>
::GetSplit(unsigned int i, unsigned int numberOfPieces,
           const RegionType &region)
{
  RegionType splitRegion;
  IndexType splitIndex;
  SizeType splitSize, regionSize;

  // Initialize the splitRegion to the requested region
  splitRegion = region;
  splitIndex = splitRegion.GetIndex();
  splitSize = splitRegion.GetSize();

  regionSize = region.GetSize();

  // Keep the number of pieces from going over the 
  // number of pixels.
  // The user will probably be surprised by this.
  unsigned long numPixels;
  numPixels = 1;
  for ( unsigned long dim = 0;
    dim < VImageDimension; 
    ++dim )
    numPixels *= regionSize[dim];
  if ( numberOfPieces > numPixels )
    numberOfPieces = numPixels;

  // I'm not as smart as whoever wrote the VTK splitter.
  // I'm going to make a big array of sizes and indecies 
  // and search through it a lot. I'll justify it by 
  // saying that hopefully I'll need all of the splits 
  // at once.

  // For some reason I can't make an array of RegionType, IndexType, 
  // or even SizeType::SizeValueType.
  typedef itk::Array< unsigned long > ArrayOfSizeType;
  ArrayOfSizeType sizes( numberOfPieces * VImageDimension );

  typedef itk::Array< unsigned long > ArrayOfIndexType;
  ArrayOfIndexType indecies( numberOfPieces * VImageDimension );

  unsigned long num_pieces_so_far = 0;

  // Set the first piece.
  for (unsigned long dim = 0; dim < VImageDimension; ++dim)
    {
    sizes[ num_pieces_so_far * VImageDimension + dim ] 
      = region.GetSize()[ dim ];

    indecies[ num_pieces_so_far * VImageDimension + dim ] 
      = region.GetIndex()[ dim ];
    }
  ++num_pieces_so_far;

  // Split untill we've made the right number.
  while ( num_pieces_so_far < numberOfPieces )
    {
    // Find the longest dim.
    unsigned long long_dim = 0;
    unsigned long long_dim_piece = 0;
    unsigned long long_dim_dim = 0;
    // Loop over image dim
    for ( unsigned long dim = 0;
      dim < VImageDimension;
      ++dim )
      {
      // Loop over piece
      for ( unsigned long piece = 0;
        piece < num_pieces_so_far;
        ++piece )
        {
        // Check the dim size
        if ( sizes[piece * VImageDimension + dim] > long_dim )
          {
          long_dim = sizes[piece * VImageDimension + dim];
          long_dim_piece = piece;
          long_dim_dim = dim;
          } // End if long_dim
        } // End for piece
      } // End for dim
    // Copy the longest piece
    for ( unsigned long dim = 0;
      dim < VImageDimension;
      ++dim )
      {
      sizes[num_pieces_so_far * VImageDimension + dim] 
        = sizes[long_dim_piece * VImageDimension + dim];
      indecies[num_pieces_so_far * VImageDimension + dim] 
        = indecies[long_dim_piece * VImageDimension + dim];
      }
    // Split the longest dim of the copy
    sizes[num_pieces_so_far * VImageDimension + long_dim_dim] 
      -= sizes[long_dim_piece * VImageDimension + long_dim_dim] / 2;
    indecies[num_pieces_so_far * VImageDimension + long_dim_dim] 
      += sizes[long_dim_piece * VImageDimension + long_dim_dim] / 2;
    // Split the origional piece
    sizes[long_dim_piece * VImageDimension + long_dim_dim] /= 2;
    ++num_pieces_so_far;
    } // End while num_pieces_so_far

  // Copy out the piece
  for ( unsigned long dim = 0; dim < VImageDimension; ++dim )
    {
    splitIndex[dim] = indecies[i * VImageDimension + dim];
    splitSize[dim] = sizes[i * VImageDimension + dim];
    }
  // set the split region ivars
  splitRegion.SetIndex( splitIndex );
  splitRegion.SetSize( splitSize );

  itkDebugMacro("  Split Piece: " << std::endl << splitRegion );

  return splitRegion;
}

/**
 *
 */
template <unsigned int VImageDimension>
void 
ImageRegionNonCubeSplitter<VImageDimension>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
}

} // end namespace itk

#endif
-------------- next part --------------
A non-text attachment was scrubbed...
Name: MPIFastMarching.cxx
Type: text/x-c++src
Size: 15107 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091023/0710c62f/attachment-0001.cxx>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: CMakeLists.txt
Type: text/x-cmake
Size: 445 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091023/0710c62f/attachment-0001.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091023/0710c62f/attachment-0001.pgp>