Hi All,<br><br>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 <b>bold</b> font. The rest is just what you can find in vtkImageGradientMagnitude). <br>
<br>It sounds like a threading issue because the range of the nested iterating loops over pixels are not the same as the input image?!<br><br>I would be so grateful for any clue.<br><br>Regards,<br><br>Mehdi<br><br><br>template <class T><br>
void vtkImageAndImageGradientMagnitudeExecute(vtkImageAndImageGradientMagnitude *self,<br> vtkImageData *inData, T *inPtr,<br> vtkImageData *outData, T *outPtr,<br>
int outExt[6], int id)<br>{<br> int idxC, idxX, idxY, idxZ;<br> int maxC, maxX, maxY, maxZ;<br> vtkIdType inIncX, inIncY, inIncZ;<br> vtkIdType outIncX, outIncY, outIncZ;<br> unsigned long count = 0;<br>
unsigned long target;<br> int axesNum;<br> int *wholeExtent;<br> vtkIdType *inIncs;<br> double r[3], d, sum;<br> int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;<br> int *inExt = inData->GetExtent();<br>
<br> // find the region to loop over<br> maxC = outData->GetNumberOfScalarComponents();<br> maxX = outExt[1] - outExt[0];<br> maxY = outExt[3] - outExt[2];<br> maxZ = outExt[5] - outExt[4];<br> target = static_cast<unsigned long>((maxZ+1)*(maxY+1)/50.0);<br>
target++;<br><br> // Get the dimensionality of the gradient.<br> axesNum = self->GetDimensionality();<br> <br> // Get increments to march through data <br> inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);<br>
outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);<br><br> // The data spacing is important for computing the gradient.<br> inData->GetSpacing(r);<br> r[0] = 0.5 / r[0];<br> r[1] = 0.5 / r[1];<br>
r[2] = 0.5 / r[2];<br><br> // get some other info we need<br> inIncs = inData->GetIncrements(); <br> wholeExtent = inData->GetExtent(); <br><br> // Move the starting pointer to the correct location.<br> inPtr += (outExt[0]-inExt[0])*inIncs[0] +<br>
(outExt[2]-inExt[2])*inIncs[1] +<br> (outExt[4]-inExt[4])*inIncs[2];<br> <br> // Loop through ouput pixels<br> for (idxZ = 0; idxZ <= maxZ; idxZ++)<br> {<br> useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];<br>
useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];<br> for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)<br> {<br> if (!id) <br> {<br> if (!(count%target))<br>
{<br> self->UpdateProgress(count/(50.0*target));<br> }<br> count++;<br> }<br> useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];<br> useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];<br>
for (idxX = 0; idxX <= maxX; idxX++)<br> {<br> useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];<br> useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];<br>
for (idxC = 0; idxC < maxC; idxC++)<br> {<br> // do X axis<br> d = static_cast<double>(inPtr[useXMin]);<br> d -= static_cast<double>(inPtr[useXMax]);<br> d *= r[0]; // multiply by the data spacing<br>
sum = d * d;<br> // do y axis<br> d = static_cast<double>(inPtr[useYMin]);<br> d -= static_cast<double>(inPtr[useYMax]);<br> d *= r[1]; // multiply by the data spacing<br>
sum += (d * d);<br> if (axesNum == 3)<br> {<br> // do z axis<br> d = static_cast<double>(inPtr[useZMin]);<br> d -= static_cast<double>(inPtr[useZMax]);<br>
d *= r[2]; // multiply by the data spacing<br> sum += (d * d);<br> }<br> //*outPtr = static_cast<T>(sqrt(sum));<br><br> // cuting out the gradient circule from the image<br>
int currentPoint[3] = {idxX, idxY, idxZ};<br> int currentExtent[3] = {maxX, maxY, maxZ};<br><br> <br> <br><b> if ( self->PointIsInCircle( currentPoint, currentExtent ) )<br> {<br>
*outPtr = static_cast<T>( sqrt( sum ) );<br> }<br> else<br> {<br> *outPtr = *inPtr;<br> }<br> <br> outPtr++;<br> inPtr++;<br> }<br>
}<br><br></b> outPtr += outIncY;<br> inPtr += inIncY;<br> }<br> outPtr += outIncZ;<br> inPtr += inIncZ;<br> }<br> std::cout << "I am running without crash!" << std::endl;<br>
}<br><br><br><b>//----------------------------------------------------------------------------<br>bool vtkImageAndImageGradientMagnitude::PointIsInCircle(int currentPoint[3], int currentExtent[3])<br>{<br> int tmp = (currentExtent[0] < currentExtent[1])? currentExtent[0]: currentExtent[1];<br>
<br> if ( GetDimensionality() > 2)<br> {<br> radius = (tmp < currentExtent[2])? tmp: currentExtent[2];<br> }<br> else<br> {<br> radius = tmp;<br> }<br> <br> radius *= 0.60;<br>
<br> int centX = currentExtent[0];<br> int centY = currentExtent[1];<br> //int centY = 0;<br><br> int length = pow((centX - currentPoint[0])*1.0, 2) + pow((centY - currentPoint[1]*1.0), 2);<br><br> if ( GetDimensionality() > 2)<br>
{<br> int centZ = currentExtent[2] / 2;<br> length += pow((centZ - currentPoint[2])*1.0, 2);<br> }<br><br> return (length < pow(radius*1.0, 2));<br>}<br></b><br>