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