1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageGradientMagnitude.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 "vtkImageGradientMagnitude.h"
16
17 #include "vtkDataArray.h"
18 #include "vtkImageData.h"
19 #include "vtkInformation.h"
20 #include "vtkInformationVector.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkPointData.h"
23 #include "vtkStreamingDemandDrivenPipeline.h"
24
25 #include <cmath>
26
27 vtkStandardNewMacro(vtkImageGradientMagnitude);
28
29 //------------------------------------------------------------------------------
30 // Construct an instance of vtkImageGradientMagnitude filter.
vtkImageGradientMagnitude()31 vtkImageGradientMagnitude::vtkImageGradientMagnitude()
32 {
33 this->SetNumberOfInputPorts(1);
34 this->SetNumberOfOutputPorts(1);
35 this->Dimensionality = 2;
36 this->HandleBoundaries = 1;
37 }
38
39 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)40 void vtkImageGradientMagnitude::PrintSelf(ostream& os, vtkIndent indent)
41 {
42 this->Superclass::PrintSelf(os, indent);
43 os << indent << "HandleBoundaries: " << this->HandleBoundaries << "\n";
44 os << indent << "Dimensionality: " << this->Dimensionality << "\n";
45 }
46
47 //------------------------------------------------------------------------------
48 // This method is passed a region that holds the image extent of this filters
49 // input, and changes the region to hold the image extent of this filters
50 // output.
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)51 int vtkImageGradientMagnitude::RequestInformation(vtkInformation* vtkNotUsed(request),
52 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
53 {
54 int extent[6];
55 int idx;
56
57 // get the info objects
58 vtkInformation* outInfo = outputVector->GetInformationObject(0);
59 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
60
61 // invalid setting, it has not been set, so default to whole Extent
62 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
63
64 if (!this->HandleBoundaries)
65 {
66 // shrink output image extent.
67 for (idx = 0; idx < this->Dimensionality; ++idx)
68 {
69 extent[idx * 2] += 1;
70 extent[idx * 2 + 1] -= 1;
71 }
72 }
73
74 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);
75
76 return 1;
77 }
78
79 //------------------------------------------------------------------------------
80 // This method computes the input extent necessary to generate the output.
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)81 int vtkImageGradientMagnitude::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
82 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
83 {
84 int wholeExtent[6];
85 int idx;
86
87 // get the info objects
88 vtkInformation* outInfo = outputVector->GetInformationObject(0);
89 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
90
91 // invalid setting, it has not been set, so default to whole Extent
92 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
93 int inUExt[6];
94 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt);
95
96 // grow input whole extent.
97 for (idx = 0; idx < this->Dimensionality; ++idx)
98 {
99 inUExt[idx * 2] -= 1;
100 inUExt[idx * 2 + 1] += 1;
101 if (this->HandleBoundaries)
102 {
103 // we must clip extent with whole extent is we handle boundaries.
104 if (inUExt[idx * 2] < wholeExtent[idx * 2])
105 {
106 inUExt[idx * 2] = wholeExtent[idx * 2];
107 }
108 if (inUExt[idx * 2 + 1] > wholeExtent[idx * 2 + 1])
109 {
110 inUExt[idx * 2 + 1] = wholeExtent[idx * 2 + 1];
111 }
112 }
113 }
114 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt, 6);
115
116 return 1;
117 }
118
119 //------------------------------------------------------------------------------
120 // This execute method handles boundaries.
121 // it handles boundaries. Pixels are just replicated to get values
122 // out of extent.
123 template <class T>
vtkImageGradientMagnitudeExecute(vtkImageGradientMagnitude * self,vtkImageData * inData,T * inPtr,vtkImageData * outData,T * outPtr,int outExt[6],int id)124 void vtkImageGradientMagnitudeExecute(vtkImageGradientMagnitude* self, vtkImageData* inData,
125 T* inPtr, vtkImageData* outData, T* outPtr, int outExt[6], int id)
126 {
127 int idxC, idxX, idxY, idxZ;
128 int maxC, maxX, maxY, maxZ;
129 vtkIdType inIncX, inIncY, inIncZ;
130 vtkIdType outIncX, outIncY, outIncZ;
131 unsigned long count = 0;
132 unsigned long target;
133 int axesNum;
134 int* wholeExtent;
135 vtkIdType inIncs[3];
136 double r[3], d, sum;
137 int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;
138 int* inExt = inData->GetExtent();
139
140 // find the region to loop over
141 maxC = outData->GetNumberOfScalarComponents();
142 maxX = outExt[1] - outExt[0];
143 maxY = outExt[3] - outExt[2];
144 maxZ = outExt[5] - outExt[4];
145 target = static_cast<unsigned long>((maxZ + 1) * (maxY + 1) / 50.0);
146 target++;
147
148 // Get the dimensionality of the gradient.
149 axesNum = self->GetDimensionality();
150
151 // Get increments to march through data
152 inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
153 outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);
154
155 // The data spacing is important for computing the gradient.
156 inData->GetSpacing(r);
157 r[0] = 0.5 / r[0];
158 r[1] = 0.5 / r[1];
159 r[2] = 0.5 / r[2];
160
161 // get some other info we need
162 inData->GetIncrements(inIncs);
163 wholeExtent = inData->GetExtent();
164
165 // Move the starting pointer to the correct location.
166 inPtr += (outExt[0] - inExt[0]) * inIncs[0] + (outExt[2] - inExt[2]) * inIncs[1] +
167 (outExt[4] - inExt[4]) * inIncs[2];
168
169 // Loop through output pixels
170 for (idxZ = 0; idxZ <= maxZ; idxZ++)
171 {
172 useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];
173 useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];
174 for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)
175 {
176 if (!id)
177 {
178 if (!(count % target))
179 {
180 self->UpdateProgress(count / (50.0 * target));
181 }
182 count++;
183 }
184 useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];
185 useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];
186 for (idxX = 0; idxX <= maxX; idxX++)
187 {
188 useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];
189 useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];
190 for (idxC = 0; idxC < maxC; idxC++)
191 {
192 // do X axis
193 d = static_cast<double>(inPtr[useXMin]);
194 d -= static_cast<double>(inPtr[useXMax]);
195 d *= r[0]; // multiply by the data spacing
196 sum = d * d;
197 // do y axis
198 d = static_cast<double>(inPtr[useYMin]);
199 d -= static_cast<double>(inPtr[useYMax]);
200 d *= r[1]; // multiply by the data spacing
201 sum += (d * d);
202 if (axesNum == 3)
203 {
204 // do z axis
205 d = static_cast<double>(inPtr[useZMin]);
206 d -= static_cast<double>(inPtr[useZMax]);
207 d *= r[2]; // multiply by the data spacing
208 sum += (d * d);
209 }
210 *outPtr = static_cast<T>(sqrt(sum));
211 outPtr++;
212 inPtr++;
213 }
214 }
215 outPtr += outIncY;
216 inPtr += inIncY;
217 }
218 outPtr += outIncZ;
219 inPtr += inIncZ;
220 }
221 }
222
223 //------------------------------------------------------------------------------
224 // This method contains a switch statement that calls the correct
225 // templated function for the input data type. The output data
226 // must match input type. This method does handle boundary conditions.
ThreadedExecute(vtkImageData * inData,vtkImageData * outData,int outExt[6],int id)227 void vtkImageGradientMagnitude::ThreadedExecute(
228 vtkImageData* inData, vtkImageData* outData, int outExt[6], int id)
229 {
230 void* inPtr;
231 void* outPtr = outData->GetScalarPointerForExtent(outExt);
232 inPtr = inData->GetScalarPointer();
233
234 // this filter expects that input is the same type as output.
235 if (inData->GetScalarType() != outData->GetScalarType())
236 {
237 vtkErrorMacro(<< "Execute: input data type, " << inData->GetScalarType()
238 << ", must match out ScalarType " << outData->GetScalarType());
239 return;
240 }
241
242 switch (inData->GetScalarType())
243 {
244 vtkTemplateMacro(vtkImageGradientMagnitudeExecute(this, inData, static_cast<VTK_TT*>(inPtr),
245 outData, static_cast<VTK_TT*>(outPtr), outExt, id));
246 default:
247 vtkErrorMacro(<< "Execute: Unknown ScalarType");
248 return;
249 }
250 }
251