1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageGradient.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 "vtkImageGradient.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 <math.h>
26 #include <vtksys/ios/sstream>
27
28 vtkStandardNewMacro(vtkImageGradient);
29
30 //----------------------------------------------------------------------------
31 // Construct an instance of vtkImageGradient fitler.
vtkImageGradient()32 vtkImageGradient::vtkImageGradient()
33 {
34 this->HandleBoundaries = 1;
35 this->Dimensionality = 2;
36
37 // by default process active point scalars
38 this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
39 vtkDataSetAttributes::SCALARS);
40 }
41
42 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)43 void vtkImageGradient::PrintSelf(ostream& os, vtkIndent indent)
44 {
45 this->Superclass::PrintSelf(os, indent);
46 os << indent << "HandleBoundaries: " << this->HandleBoundaries << "\n";
47 os << indent << "Dimensionality: " << this->Dimensionality << "\n";
48 }
49
50 //----------------------------------------------------------------------------
RequestInformation(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)51 int vtkImageGradient::RequestInformation(vtkInformation*,
52 vtkInformationVector** inputVector,
53 vtkInformationVector* outputVector)
54 {
55 // Get input and output pipeline information.
56 vtkInformation* outInfo = outputVector->GetInformationObject(0);
57 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
58
59 // Get the input whole extent.
60 int extent[6];
61 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
62
63 // Shrink output image extent by one pixel if not handling boundaries.
64 if(!this->HandleBoundaries)
65 {
66 for(int idx = 0; idx < this->Dimensionality; ++idx)
67 {
68 extent[idx*2] += 1;
69 extent[idx*2 + 1] -= 1;
70 }
71 }
72
73 // Store the new whole extent for the output.
74 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);
75
76 // Set the number of point data componets to the number of
77 // components in the gradient vector.
78 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_DOUBLE,
79 this->Dimensionality);
80
81 return 1;
82 }
83
84 //----------------------------------------------------------------------------
85 // This method computes the input extent necessary to generate the output.
RequestUpdateExtent(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)86 int vtkImageGradient::RequestUpdateExtent(vtkInformation*,
87 vtkInformationVector** inputVector,
88 vtkInformationVector* outputVector)
89 {
90 // Get input and output pipeline information.
91 vtkInformation* outInfo = outputVector->GetInformationObject(0);
92 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
93
94 // We need one extra ghost level
95 int ugl = outInfo->Get(
96 vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
97 inInfo->Set(
98 vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), ugl + 1);
99
100 // Get the input whole extent.
101 int wholeExtent[6];
102 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
103
104 // Get the requested update extent from the output.
105 int inUExt[6];
106 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt);
107
108 // In order to do central differencing we need one more layer of
109 // input pixels than we are producing output pixels.
110 for(int idx = 0; idx < this->Dimensionality; ++idx)
111 {
112 inUExt[idx*2] -= 1;
113 inUExt[idx*2+1] += 1;
114
115 // If handling boundaries instead of shrinking the image then we
116 // must clip the needed extent within the whole extent of the
117 // input.
118 if (this->HandleBoundaries)
119 {
120 if (inUExt[idx*2] < wholeExtent[idx*2])
121 {
122 inUExt[idx*2] = wholeExtent[idx*2];
123 }
124 if (inUExt[idx*2 + 1] > wholeExtent[idx*2 + 1])
125 {
126 inUExt[idx*2 + 1] = wholeExtent[idx*2 + 1];
127 }
128 }
129 }
130
131 // Store the update extent needed from the intput.
132 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt, 6);
133
134 return 1;
135 }
136
137 //----------------------------------------------------------------------------
138 // This execute method handles boundaries.
139 // it handles boundaries. Pixels are just replicated to get values
140 // out of extent.
141 template <class T>
vtkImageGradientExecute(vtkImageGradient * self,vtkImageData * inData,T * inPtr,vtkImageData * outData,double * outPtr,int outExt[6],int id)142 void vtkImageGradientExecute(vtkImageGradient *self,
143 vtkImageData *inData, T *inPtr,
144 vtkImageData *outData, double *outPtr,
145 int outExt[6], int id)
146 {
147 int idxX, idxY, idxZ;
148 int maxX, maxY, maxZ;
149 vtkIdType inIncX, inIncY, inIncZ;
150 vtkIdType outIncX, outIncY, outIncZ;
151 unsigned long count = 0;
152 unsigned long target;
153 int axesNum;
154 int *inExt = inData->GetExtent();
155 int *wholeExtent;
156 vtkIdType *inIncs;
157 double r[3], d;
158 int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;
159
160 // find the region to loop over
161 maxX = outExt[1] - outExt[0];
162 maxY = outExt[3] - outExt[2];
163 maxZ = outExt[5] - outExt[4];
164 target = static_cast<unsigned long>((maxZ+1)*(maxY+1)/50.0);
165 target++;
166
167 // Get the dimensionality of the gradient.
168 axesNum = self->GetDimensionality();
169
170 // Get increments to march through data
171 inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
172 outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);
173
174 // The data spacing is important for computing the gradient.
175 // central differences (2 * ratio).
176 // Negative because below we have (min - max) for dx ...
177 inData->GetSpacing(r);
178 r[0] = -0.5 / r[0];
179 r[1] = -0.5 / r[1];
180 r[2] = -0.5 / r[2];
181
182 // get some other info we need
183 inIncs = inData->GetIncrements();
184 wholeExtent = inData->GetExtent();
185
186 // Move the pointer to the correct starting position.
187 inPtr += (outExt[0]-inExt[0])*inIncs[0] +
188 (outExt[2]-inExt[2])*inIncs[1] +
189 (outExt[4]-inExt[4])*inIncs[2];
190
191 // Loop through output pixels
192 for (idxZ = 0; idxZ <= maxZ; idxZ++)
193 {
194 useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];
195 useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];
196 for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)
197 {
198 if (!id)
199 {
200 if (!(count%target))
201 {
202 self->UpdateProgress(count/(50.0*target));
203 }
204 count++;
205 }
206 useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];
207 useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];
208 for (idxX = 0; idxX <= maxX; idxX++)
209 {
210 useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];
211 useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];
212
213 // do X axis
214 d = static_cast<double>(inPtr[useXMin]);
215 d -= static_cast<double>(inPtr[useXMax]);
216 d *= r[0]; // multiply by the data spacing
217 *outPtr = d;
218 outPtr++;
219
220 // do y axis
221 d = static_cast<double>(inPtr[useYMin]);
222 d -= static_cast<double>(inPtr[useYMax]);
223 d *= r[1]; // multiply by the data spacing
224 *outPtr = d;
225 outPtr++;
226 if (axesNum == 3)
227 {
228 // do z axis
229 d = static_cast<double>(inPtr[useZMin]);
230 d -= static_cast<double>(inPtr[useZMax]);
231 d *= r[2]; // multiply by the data spacing
232 *outPtr = d;
233 outPtr++;
234 }
235 inPtr++;
236 }
237 outPtr += outIncY;
238 inPtr += inIncY;
239 }
240 outPtr += outIncZ;
241 inPtr += inIncZ;
242 }
243 }
244
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)245 int vtkImageGradient::RequestData(
246 vtkInformation* request,
247 vtkInformationVector** inputVector,
248 vtkInformationVector* outputVector)
249 {
250 // Shrink the update extent to the input extent. Input extent
251 // can be smaller than update extent when there is an piece
252 // request (UPDATE_NUMBER_OF_PIECES() > 1).
253 // Since the superclass and this class make decisions based
254 // on UPDATE_EXTENT(), this is the quickest way of making this
255 // filter in distributed parallel mode.
256 // In the future, this logic should move up the hierarchy so
257 // other imaging classes can benefit from it.
258 vtkImageData* input = vtkImageData::GetData(inputVector[0]);
259 vtkInformation* outInfo = outputVector->GetInformationObject(0);
260 int ue[6];
261 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ue);
262 int ue2[6];
263 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ue2);
264 int* ie = input->GetExtent();
265 for (int i=0; i<3; i++)
266 {
267 if (ue[2*i] < ie[2*i])
268 {
269 ue2[2*i] = ie[2*i];
270 }
271 if (ue[2*i+1] > ie[2*i+1])
272 {
273 ue2[2*i+1] = ie[2*i+1];
274 }
275 }
276 outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ue2, 6);
277
278 if (!this->Superclass::RequestData(request, inputVector, outputVector))
279 {
280 return 0;
281 }
282 vtkImageData* output = vtkImageData::GetData(outputVector);
283 vtkDataArray* outArray = output->GetPointData()->GetScalars();
284 vtksys_ios::ostringstream newname;
285 newname << (outArray->GetName()?outArray->GetName():"")
286 << "Gradient";
287 outArray->SetName(newname.str().c_str());
288 // Why not pass the original array?
289 if (this->GetInputArrayToProcess(0, inputVector))
290 {
291 output->GetPointData()->AddArray(
292 this->GetInputArrayToProcess(0, inputVector));
293 }
294 // Restore the previous update extent. See code above for details.
295 outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ue, 6);
296 return 1;
297 }
298
299 //----------------------------------------------------------------------------
300 // This method contains a switch statement that calls the correct
301 // templated function for the input data type. This method does handle
302 // boundary conditions.
ThreadedRequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector *,vtkImageData *** inData,vtkImageData ** outData,int outExt[6],int threadId)303 void vtkImageGradient::ThreadedRequestData(vtkInformation*,
304 vtkInformationVector** inputVector,
305 vtkInformationVector*,
306 vtkImageData*** inData,
307 vtkImageData** outData,
308 int outExt[6],
309 int threadId)
310 {
311 // Get the input and output data objects.
312 vtkImageData* input = inData[0][0];
313 vtkImageData* output = outData[0];
314
315 // The ouptut scalar type must be double to store proper gradients.
316 if(output->GetScalarType() != VTK_DOUBLE)
317 {
318 vtkErrorMacro("Execute: output ScalarType is "
319 << output->GetScalarType() << "but must be double.");
320 return;
321 }
322
323 vtkDataArray* inputArray = this->GetInputArrayToProcess(0, inputVector);
324 if (!inputArray)
325 {
326 vtkErrorMacro("No input array was found. Cannot execute");
327 return;
328 }
329
330 // Gradient makes sense only with one input component. This is not
331 // a Jacobian filter.
332 if(inputArray->GetNumberOfComponents() != 1)
333 {
334 vtkErrorMacro(
335 "Execute: input has more than one component. "
336 "The input to gradient should be a single component image. "
337 "Think about it. If you insist on using a color image then "
338 "run it though RGBToHSV then ExtractComponents to get the V "
339 "components. That's probably what you want anyhow.");
340 return;
341 }
342
343 void* inPtr = inputArray->GetVoidPointer(0);
344 double* outPtr = static_cast<double *>(
345 output->GetScalarPointerForExtent(outExt));
346 switch(inputArray->GetDataType())
347 {
348 vtkTemplateMacro(
349 vtkImageGradientExecute(this, input, static_cast<VTK_TT*>(inPtr),
350 output, outPtr, outExt, threadId)
351 );
352 default:
353 vtkErrorMacro("Execute: Unknown ScalarType " << input->GetScalarType());
354 return;
355 }
356 }
357