[vtkusers] I need a way to quickly but semi-randomly access image data

[Mark Roden]

Hi all,

I'm looking to speed up the vtkDijkstraImageGeodesicPath::BuildAdjacency
method.  Right now, on a 512x512 image, that takes many seconds to compute.

It looks like the biggest cost (apart from a sqrt and a division at each
pixel) is the line


  double cost = this->ImageWeight*(
    image->GetScalarComponentAsDouble( ijk1[0], ijk1[1], ijk1[2], 0 ) +
    image->GetScalarComponentAsDouble( ijk1[0], ijk1[1], ijk1[2], 0 ) );

Because GetScalarComponentAsDouble is a huge time sink when done over all
pixels at once.  I'm hoping to unroll the loop there.

I know that all of the pixels I'm going to access are safe.  I want to get a
double or floating point value, but I don't know the type of the array in
the image (or do I?).

So far, I have:

  vtkDataArray* theImageData = image->GetPointData()->GetScalars();
  vtkIdType increments[3];
  image->GetArrayIncrements(theImageData, increments);
  void* imageArrayPointer = image->GetArrayPointer(theImageData, theOrigin);

  double result1, result2;
  const int theSize = static_cast<int>(this->Internals->Adjacency.size());
  for( int u = 0; u < theSize; ++u ) {

      int v = (*it).first;
      double p1[3];
      image->GetPoint( u, p1 );
      double p2[3];
      image->GetPoint( v, p2 );

      int ijk1[3];
      int ijk2[3];

      for (i=0; i<3; i++)
        {
          d = p1[i] - origin[i];
          doubleLoc = d / spacing[i];
          // Floor for negative indexes.
          ijk1[i] = static_cast<int>(floor(doubleLoc));
          d = p2[i] - origin[i];
          doubleLoc = d / spacing[i];
          // Floor for negative indexes.
          ijk2[i] = static_cast<int>(floor(doubleLoc));
        }
      switch (image->GetScalarType())
        {
        vtkTemplateMacro(vtkImageDataConvertScalar(imageArrayPointer[ijk1[0]
+ ijk1[1]*increments[0] + ijk1[2]*increments[1]],&result1));
        vtkTemplateMacro(vtkImageDataConvertScalar(imageArrayPointer[ijk2[0]
+ ijk2[1]*increments[0] + ijk2[2]*increments[1]],&result1));
        default:
        {
          vtkErrorMacro("Unknown Scalar type " << this->GetScalarType());
        }
      }

which is cobbled together from the various methods that are used in the
computation of static edge costs in that file.

Of course, that's not right.  What are the magic words to make this go?  Do
I need to template this function to the possible types of vtkImageData
types? I'm thinking something like:

template <type>
type GetValueAtPoint(int loc){
  return ArrayOfType[loc];
}

I mean, even better if I could get an iterator or something that I could
just increment to the appropriate location.

Thanks,
Mark
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20101124/8e447ed5/attachment.htm>