1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkImageNonMaximumSuppression.cxx
5 
6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7   All rights reserved.
8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10      This software is distributed WITHOUT ANY WARRANTY; without even
11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12      PURPOSE.  See the above copyright notice for more information.
13 
14 =========================================================================*/
15 #include "vtkImageNonMaximumSuppression.h"
16 
17 #include "vtkDataArray.h"
18 #include "vtkImageData.h"
19 #include "vtkInformation.h"
20 #include "vtkInformationVector.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkStreamingDemandDrivenPipeline.h"
23 #include "vtkPointData.h"
24 
25 #include <cmath>
26 
27 vtkStandardNewMacro(vtkImageNonMaximumSuppression);
28 
29 //----------------------------------------------------------------------------
30 // Construct an instance of vtkImageNonMaximumSuppression filter.
vtkImageNonMaximumSuppression()31 vtkImageNonMaximumSuppression::vtkImageNonMaximumSuppression()
32 {
33   this->Dimensionality= 2;
34   this->HandleBoundaries = 1;
35   this->SetNumberOfInputPorts(2);
36 }
37 
38 //----------------------------------------------------------------------------
39 // This method is passed a region that holds the image extent of this filters
40 // input, and changes the region to hold the image extent of this filters
41 // output.
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)42 int vtkImageNonMaximumSuppression::RequestInformation (
43   vtkInformation * vtkNotUsed(request),
44   vtkInformationVector **inputVector,
45   vtkInformationVector *outputVector)
46 {
47   // get the info objects
48   vtkInformation* outInfo = outputVector->GetInformationObject(0);
49   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
50 
51   int extent[6];
52   int idx;
53 
54   inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent);
55   if ( ! this->HandleBoundaries)
56   {
57     // shrink output image extent.
58     for (idx = 0; idx < this->Dimensionality; ++idx)
59     {
60       extent[idx*2] += 1;
61       extent[idx*2+1] -= 1;
62     }
63   }
64 
65   outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
66 
67   return 1;
68 }
69 
70 //----------------------------------------------------------------------------
71 // This method computes the input extent necessary to generate the output.
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)72 int vtkImageNonMaximumSuppression::RequestUpdateExtent (
73   vtkInformation * vtkNotUsed(request),
74   vtkInformationVector **inputVector,
75   vtkInformationVector *outputVector)
76 {
77   // get the info objects
78   vtkInformation* outInfo = outputVector->GetInformationObject(0);
79   vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
80   vtkInformation* inInfo2 = inputVector[1]->GetInformationObject(0);
81 
82   int *wholeExtent;
83   int idx;
84 
85   // get the whole image for input 2
86   int inExt[6];
87   outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),inExt);
88   wholeExtent = inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
89   inInfo2->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),inExt,6);
90 
91   // grow input image extent for input 0
92   for (idx = 0; idx < this->Dimensionality; ++idx)
93   {
94     inExt[idx*2] -= 1;
95     inExt[idx*2+1] += 1;
96     if (this->HandleBoundaries)
97     {
98       // we must clip extent with whole extent if we hanlde boundaries.
99       if (inExt[idx*2] < wholeExtent[idx*2])
100       {
101         inExt[idx*2] = wholeExtent[idx*2];
102       }
103       if (inExt[idx*2 + 1] > wholeExtent[idx*2 + 1])
104       {
105         inExt[idx*2 + 1] = wholeExtent[idx*2 + 1];
106       }
107     }
108   }
109   inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),inExt,6);
110 
111   return 1;
112 }
113 
114 //----------------------------------------------------------------------------
115 // This templated function executes the filter for any type of data.
116 // Handles the two input operations
117 template <class T>
vtkImageNonMaximumSuppressionExecute(vtkImageNonMaximumSuppression * self,vtkImageData * in1Data,T * in1Ptr,vtkImageData * in2Data,T * in2Ptr,vtkImageData * outData,T * outPtr,int outExt[6],int id)118 void vtkImageNonMaximumSuppressionExecute(vtkImageNonMaximumSuppression *self,
119                                           vtkImageData *in1Data, T *in1Ptr,
120                                           vtkImageData *in2Data,
121                                           T *in2Ptr,
122                                           vtkImageData *outData,
123                                           T *outPtr,
124                                           int outExt[6], int id)
125 {
126   int idxC, idxX, idxY, idxZ;
127   int maxC, maxX, maxY, maxZ;
128   vtkIdType inIncX, inIncY, inIncZ;
129   vtkIdType in2IncX, in2IncY, in2IncZ;
130   vtkIdType outIncX, outIncY, outIncZ;
131   unsigned long count = 0;
132   unsigned long target;
133   int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;
134   double d, normalizeFactor, vector[3], *ratio;
135   int neighborA, neighborB;
136   int *wholeExtent;
137   vtkIdType inIncs[3];
138   int axesNum;
139 
140   vector[0] = 0.0;
141   vector[1] = 0.0;
142   vector[2] = 0.0;
143 
144   // find the region to loop over
145   maxC = outData->GetNumberOfScalarComponents();
146   maxX = outExt[1] - outExt[0];
147   maxY = outExt[3] - outExt[2];
148   maxZ = outExt[5] - outExt[4];
149   target = static_cast<unsigned long>((maxZ+1)*(maxY+1)/50.0);
150   target++;
151 
152   // Get the dimensionality of the gradient.
153   axesNum = self->GetDimensionality();
154   // get some other info we need
155   in1Data->GetIncrements(inIncs);
156   wholeExtent = in1Data->GetExtent();
157 
158   // Get increments to march through data
159   in1Data->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
160   in2Data->GetContinuousIncrements(outExt, in2IncX, in2IncY, in2IncZ);
161   outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);
162 
163 
164   // Gradient is computed with data spacing (world coordinates)
165   ratio = in2Data->GetSpacing();
166 
167   // Loop through output pixels
168   for (idxZ = 0; idxZ <= maxZ; idxZ++)
169   {
170     useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];
171     useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];
172     for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)
173     {
174       useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];
175       useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];
176       if (!id)
177       {
178         if (!(count%target))
179         {
180           self->UpdateProgress(count/(50.0*target));
181         }
182         count++;
183       }
184       for (idxX = 0; idxX <= maxX; idxX++)
185       {
186         useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];
187         useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];
188 
189         // calculate the neighbors
190         d = vector[0] = static_cast<double>(*in2Ptr) * ratio[0];
191         normalizeFactor = (d * d);
192         d = vector[1] = static_cast<double>(in2Ptr[1]) * ratio[1];
193         normalizeFactor += (d * d);
194         if (axesNum == 3)
195         {
196           d = vector[2] = static_cast<double>(in2Ptr[2]) * ratio[2];
197           normalizeFactor += (d * d);
198         }
199         if (normalizeFactor != 0.0)
200         {
201           normalizeFactor = 1.0 / sqrt(normalizeFactor);
202         }
203         // Vector points positive along this idx?
204         // (can point along multiple axes)
205         d = vector[0] * normalizeFactor;
206 
207         if (d > 0.5)
208         {
209           neighborA = useXMax;
210           neighborB = useXMin;
211         }
212         else if (d < -0.5)
213         {
214           neighborB = useXMax;
215           neighborA = useXMin;
216         }
217         else
218         {
219           neighborA = 0;
220           neighborB = 0;
221         }
222         d = vector[1] * normalizeFactor;
223         if (d > 0.5)
224         {
225           neighborA += useYMax;
226           neighborB += useYMin;
227         }
228         else if (d < -0.5)
229         {
230           neighborB += useYMax;
231           neighborA += useYMin;
232         }
233         if (axesNum == 3)
234         {
235           d = vector[2] * normalizeFactor;
236           if (d > 0.5)
237           {
238             neighborA += useZMax;
239             neighborB += useZMin;
240           }
241           else if (d < -0.5)
242           {
243             neighborB += useZMax;
244             neighborA += useZMin;
245           }
246         }
247 
248         // now process the components
249         for (idxC = 0; idxC < maxC; idxC++)
250         {
251           // Pixel operation
252           // Set Output Magnitude
253           if (in1Ptr[neighborA] > *in1Ptr || in1Ptr[neighborB] > *in1Ptr)
254           {
255             *outPtr = 0;
256           }
257           else
258           {
259             *outPtr = *in1Ptr;
260             // also check for them being equal is neighbor with larger ptr
261             if ((neighborA > neighborB)&&(in1Ptr[neighborA] == *in1Ptr))
262             {
263               *outPtr = 0;
264             }
265             else if ((neighborB > neighborA)&&(in1Ptr[neighborB] == *in1Ptr))
266             {
267               *outPtr = 0;
268             }
269           }
270           outPtr++;
271           in1Ptr++;
272         }
273         in2Ptr += axesNum;
274       }
275       outPtr += outIncY;
276       in1Ptr += inIncY;
277       in2Ptr += in2IncY;
278     }
279     outPtr += outIncZ;
280     in1Ptr += inIncZ;
281     in2Ptr += in2IncZ;
282   }
283 }
284 
285 //----------------------------------------------------------------------------
286 // This method is passed a input and output regions, and executes the filter
287 // algorithm to fill the output from the inputs.
288 // It just executes a switch statement to call the correct function for
289 // the regions data types.
ThreadedRequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * vtkNotUsed (outputVector),vtkImageData *** inData,vtkImageData ** outData,int outExt[6],int id)290 void vtkImageNonMaximumSuppression::ThreadedRequestData(
291   vtkInformation * vtkNotUsed( request ),
292   vtkInformationVector ** vtkNotUsed( inputVector ),
293   vtkInformationVector * vtkNotUsed( outputVector ),
294   vtkImageData ***inData,
295   vtkImageData **outData,
296   int outExt[6], int id)
297 {
298   void *in1Ptr;
299   void *in2Ptr;
300   void *outPtr;
301 
302   if (id == 0)
303   {
304     if (outData[0]->GetPointData()->GetScalars())
305     {
306       outData[0]->GetPointData()->GetScalars()->SetName("SuppressedMaximum");
307     }
308   }
309 
310   in1Ptr = inData[0][0]->GetScalarPointerForExtent(outExt);
311   in2Ptr = inData[1][0]->GetScalarPointerForExtent(outExt);
312   outPtr = outData[0]->GetScalarPointerForExtent(outExt);
313 
314   // this filter expects that input is the same type as output.
315   if (inData[0][0]->GetScalarType() != outData[0]->GetScalarType() ||
316       inData[1][0]->GetScalarType() != outData[0]->GetScalarType())
317   {
318     vtkErrorMacro(<< "Execute: input ScalarType, " <<
319     inData[0][0]->GetScalarType()
320     << ", must match out ScalarType " << outData[0]->GetScalarType());
321     return;
322   }
323 
324   switch (inData[0][0]->GetScalarType())
325   {
326     vtkTemplateMacro(
327       vtkImageNonMaximumSuppressionExecute(this, inData[0][0],
328                                            static_cast<VTK_TT *>(in1Ptr),
329                                            inData[1][0],
330                                            static_cast<VTK_TT *>(in2Ptr),
331                                            outData[0],
332                                            static_cast<VTK_TT *>(outPtr),
333                                            outExt, id));
334     default:
335       vtkErrorMacro(<< "Execute: Unknown ScalarType");
336       return;
337   }
338 }
339 
PrintSelf(ostream & os,vtkIndent indent)340 void vtkImageNonMaximumSuppression::PrintSelf(ostream& os, vtkIndent indent)
341 {
342   this->Superclass::PrintSelf(os,indent);
343 
344   os << indent << "Dimensionality: " << this->Dimensionality << "\n";
345 
346   os << indent << "HandleBoundaries: " << (this->HandleBoundaries ? "On\n" : "Off\n");
347 }
348