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