[Insight-developers] Neighborhood iterator refactoring request for comment

[Joshua Cates]

Hi Folks,

I am refactoring the neighborhood iterators to optimize them with respect
to current use and expected future use.  I expect to change their
implementation slightly (mostly transparent to the user), but also I want
to add some new functionality.

Below is a summary of the functionality that both exists in the iterators
and that I think would be nice to have based on how I've used the
iterators and how I've seen others use them.

Items which are new or represent changes are marked with +.


itk neighborhood iterators

(1)+ A single iterator class handles forward, backward, and random access,
and automatically enables/disables bounds checking at construction time 
by analyzing its iteration region.  This is different from original 
design where increasing functionality was added through subclassing.

(2) The API supports access of neighbor pixel values through
    (a) A scalar index into a linear array.  it.GetPixel(14)
   +(b) An itk::Offset from the center index.  
        it.GetPixel(itk::Offset<3>(1,0,0))
   +(c) A distance and direction along a principle axis. it.GetNext(0, 1)
    (d) Returning/setting entire neighborhood at once. it.GetNeighborhood()

(3) The API supports querying the following information about a neighbor 
pixel.
   +(a) Its image index (itk::Index) location.  
        it.GetIndex(14), it.GetIndex(itk::Offset<3>(1,0,0), etc.
   +(b) Its offset (itk::Offset) from the center of the neighborhood. 
        it.GetOffset(14), etc.
   +(c) What else do people need to know? Perhaps spatial distance from the
        center pixel?

(4) The API supports querying the following information about the 
neighborhood
    (a) Stride lengths along each neighborhood axis.  GetStride(1).
    (b) Extent of the neighborhood (radius)
    (c) Slices (std::slice) along axes.  GetSlice(2)
    (d) Size of the underlying neighborhood container, i.e. how many 
        pixels are the neighborhood.
    (e) Image location (itk::Index) of the center of the neighborhood.

(5)+ A new iterator class extends (1) to allow the creation of arbitrarily
     shaped iterators. Shaped neighborhood iterators have the following
     properties.
   (a) Iterators are optimized for linear array access, versus indexed
       lookups, though both should be possible.  This means I expect more use
       for things like neighborhood searching, summation of values, etc.
   (b) Respond to all access methods in (2).  Attempted access 
       of undefined pixel locations throws an exception.
   (c) Respond to all queries in (3). Querying on an undefined neighbor
       throws an exception.
   (d) Responds to a restricted subset of (4). GetSlice and GetStride  don't
       make sense in this context, for example.
   (e) Active or inactive neighbor locations are specified by offsets from  the
       center of the neighborhood.
   (f) Neighbor locations can be activated or deactivated after construction
       to support reshaping the neighborhood on-the-fly?
   (g) What else?

Attached below are some use case scenarios I put together which exercise
some of the existing and proposed API. Please send me your additional use
cases.

I welcome comments from anyone.  I plan to start working on these changes
next week.

Thanks,

Josh.
 

--------------------------------------------------------------------------------

typedef itk::Image<float, 3> ImageType;
ImageType::Pointer image;
  
  .
  .
  .

//
// In a 3x3 neighborhood mask, create a shaped neighborhood
//
ShapedNeighborhoodIterator<ImageType>::RadiusType radius = {1, 1, 1};
ShapedNeighborhoodIterator<ImageType> it(radius, image, 
image->GetRequestedRegion);

//
// activate the center cell
//
it->ActivateIndex(itk::Offset<3>(0,0,0));

//
// activate a list of cells which represent the 6-neighbors
//
std::list<itk::Offset<3> > active_cells;
active_cells.push_back(itk::Offset<3>(0,-1,0));
active_cells.push_back(itk::Offset<3>(0,1,0));
active_cells.push_back(itk::Offset<3>(-1,0,0));
active_cells.push_back(itk::Offset<3>(1,0,0));
active_cells.push_back(itk::Offset<3>(0,0,1));
active_cells.push_back(itk::Offset<3>(0,0,-1));

it->ActivateIndex(active_cells);


//
// Calculate all z partial derivatives using a non-shaped neighborhood
// 
itk::ConstNeighborhoodIterator<ImageType> it(region, image, 
image->GetRequestedRegion());
itk::NeighborhoodInnerProduct<ImageType> IP;
itk::DerivativeOperator<3> d_op;
d_op->SetOrder(1);
d_op->SetDirection(2);
d_op->CreateDirectional();
d_op->FlipAxes();

std::slice z_slice = it.GetSlice(2);
for (; !it.IsAtEnd(); ++it) { float df_dz = IP(z_slice, it, d_op); }


// 
// Calculate "half" forward derivatives in y direction
//
for (it.GoToBegin(); !it.IsAtEnd(); ++it) { float dy_forward = 
it.GetNext(2) - it.GetCenter(); }


//
// Sum pixel values in neighborhoods
//
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
  float sum = 0.0f;
  for (unsigned i = 0; i < it.Size(); ++i) { sum += it.GetPixel(i); }
}

// 
// Search the neighborhood for a value and and find the itk::Offset of that 
// neighbor from the center of the neighborhood.  here used to implement a 
// "flood fill" type iterator.  Can use shaped or non-shaped iterator. 
// 
std::list<itk::Index<3> > indicies_to_visit;
indicies_to_visit.push_back( itk::Index<3>(x_seed,y_seed,z_seed) ); 

while (! indicies_to_visit.empty() ) 
{
  itk::Index idx = indicies_to_visit.front()
  indicies_to_visit.pop_front();
  this->do_some_stuff(idx);
  
  it.SetLocation(idx); // or perhaps: it+=(idx - it.GetIndex());
  for (unsigned i =0; i < it.Size(); ++i)   
    {
      if ( it.GetPixel(i) == it.GetCenterPixel() )
        {
          if (this->already_visited( it.GetIndex(i) ) == false)
            indicies_to_visit.push_back( it.GetIndex(i) );
        }
    }
}



______________________________
 Josh Cates			
 School of Computer Science	
 University of Utah
 Email: cates@sci.utah.edu
 Phone: (801) 587-7697
 URL:   www.cs.utk.edu/~cates