[vtkusers] xxxxxSPAMxxxxx vtkImageGradientMagnitude

[mehdi esteghamatian]

Hi All,

I have modified vtkImageGradientMagnitude so that it calculates edge
magnitude for only one circle that is placed at the center of an input
image. However, for the rest of the image the values of the image pixel is
calculated only. As I ran the algorithm over an input image it shows two
circles which is strange. The implementation is like below (modified parts
are in *bold* font. The rest is just what you can find in
vtkImageGradientMagnitude).

It sounds like a threading issue because the range of the nested iterating
loops over pixels are not the same as the input image?!

I would be so grateful for any clue.

Regards,

Mehdi


template <class T>
void
vtkImageAndImageGradientMagnitudeExecute(vtkImageAndImageGradientMagnitude
*self,
                                      vtkImageData *inData, T *inPtr,
                                      vtkImageData *outData, T *outPtr,
                                      int outExt[6], int id)
{
  int idxC, idxX, idxY, idxZ;
  int maxC, maxX, maxY, maxZ;
  vtkIdType inIncX, inIncY, inIncZ;
  vtkIdType outIncX, outIncY, outIncZ;
  unsigned long count = 0;
  unsigned long target;
  int axesNum;
  int *wholeExtent;
  vtkIdType *inIncs;
  double r[3], d, sum;
  int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;
  int *inExt = inData->GetExtent();

  // find the region to loop over
  maxC = outData->GetNumberOfScalarComponents();
  maxX = outExt[1] - outExt[0];
  maxY = outExt[3] - outExt[2];
  maxZ = outExt[5] - outExt[4];
  target = static_cast<unsigned long>((maxZ+1)*(maxY+1)/50.0);
  target++;

  // Get the dimensionality of the gradient.
  axesNum = self->GetDimensionality();

  // Get increments to march through data
  inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
  outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);

  // The data spacing is important for computing the gradient.
  inData->GetSpacing(r);
  r[0] = 0.5 / r[0];
  r[1] = 0.5 / r[1];
  r[2] = 0.5 / r[2];

  // get some other info we need
  inIncs = inData->GetIncrements();
  wholeExtent = inData->GetExtent();

  // Move the starting pointer to the correct location.
  inPtr += (outExt[0]-inExt[0])*inIncs[0] +
           (outExt[2]-inExt[2])*inIncs[1] +
           (outExt[4]-inExt[4])*inIncs[2];

  // Loop through ouput pixels
  for (idxZ = 0; idxZ <= maxZ; idxZ++)
    {
    useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];
    useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];
    for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)
      {
      if (!id)
        {
        if (!(count%target))
          {
          self->UpdateProgress(count/(50.0*target));
          }
        count++;
        }
      useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];
      useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];
      for (idxX = 0; idxX <= maxX; idxX++)
        {
        useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];
        useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];
        for (idxC = 0; idxC < maxC; idxC++)
          {
          // do X axis
          d = static_cast<double>(inPtr[useXMin]);
          d -= static_cast<double>(inPtr[useXMax]);
          d *= r[0]; // multiply by the data spacing
          sum = d * d;
          // do y axis
          d = static_cast<double>(inPtr[useYMin]);
          d -= static_cast<double>(inPtr[useYMax]);
          d *= r[1]; // multiply by the data spacing
          sum += (d * d);
          if (axesNum == 3)
            {
            // do z axis
            d = static_cast<double>(inPtr[useZMin]);
            d -= static_cast<double>(inPtr[useZMax]);
            d *= r[2]; // multiply by the data spacing
            sum += (d * d);
            }
          //*outPtr = static_cast<T>(sqrt(sum));

          // cuting out the gradient circule from the image
          int currentPoint[3] = {idxX, idxY, idxZ};
          int currentExtent[3] = {maxX, maxY, maxZ};



*         if ( self->PointIsInCircle( currentPoint, currentExtent ) )
         {
            *outPtr = static_cast<T>( sqrt( sum ) );
         }
         else
         {
            *outPtr = *inPtr;
         }

          outPtr++;
          inPtr++;
          }
        }

*      outPtr += outIncY;
      inPtr += inIncY;
      }
    outPtr += outIncZ;
    inPtr += inIncZ;
    }
  std::cout << "I am running without crash!" << std::endl;
}


*
//----------------------------------------------------------------------------
bool vtkImageAndImageGradientMagnitude::PointIsInCircle(int currentPoint[3],
int currentExtent[3])
{
    int tmp = (currentExtent[0] < currentExtent[1])? currentExtent[0]:
currentExtent[1];

    if ( GetDimensionality() > 2)
    {
        radius = (tmp < currentExtent[2])? tmp: currentExtent[2];
    }
    else
    {
        radius = tmp;
    }

    radius *= 0.60;

    int centX = currentExtent[0];
    int centY = currentExtent[1];
    //int centY = 0;

    int length = pow((centX - currentPoint[0])*1.0, 2) + pow((centY -
currentPoint[1]*1.0), 2);

    if ( GetDimensionality() > 2)
    {
        int centZ = currentExtent[2] / 2;
        length += pow((centZ - currentPoint[2])*1.0, 2);
    }

    return (length < pow(radius*1.0, 2));
}
*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20100112/0885c01d/attachment.htm>