1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageSobel3D.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 "vtkImageSobel3D.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 <math.h>
24
25 vtkStandardNewMacro(vtkImageSobel3D);
26
27 //----------------------------------------------------------------------------
28 // Construct an instance of vtkImageSobel3D fitler.
vtkImageSobel3D()29 vtkImageSobel3D::vtkImageSobel3D()
30 {
31 this->KernelSize[0] = 3;
32 this->KernelSize[1] = 3;
33 this->KernelSize[2] = 3;
34 this->KernelMiddle[0] = 1;
35 this->KernelMiddle[1] = 1;
36 this->KernelMiddle[2] = 1;
37 this->HandleBoundaries = 1;
38 }
39
40 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)41 void vtkImageSobel3D::PrintSelf(ostream& os, vtkIndent indent)
42 {
43 this->Superclass::PrintSelf(os, indent);
44 }
45
46 //----------------------------------------------------------------------------
RequestInformation(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)47 int vtkImageSobel3D::RequestInformation (vtkInformation *request,
48 vtkInformationVector **inputVector,
49 vtkInformationVector *outputVector)
50 {
51 int retval =
52 this->Superclass::RequestInformation(request, inputVector, outputVector);
53 vtkInformation *outInfo = outputVector->GetInformationObject(0);
54 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_DOUBLE, 3);
55 return retval;
56 }
57
58 //----------------------------------------------------------------------------
59 // This execute method handles boundaries.
60 // it handles boundaries. Pixels are just replicated to get values
61 // out of extent.
62 template <class T>
vtkImageSobel3DExecute(vtkImageSobel3D * self,vtkImageData * inData,T * inPtr,vtkImageData * outData,int * outExt,double * outPtr,int id,vtkInformation * inInfo)63 void vtkImageSobel3DExecute(vtkImageSobel3D *self,
64 vtkImageData *inData, T *inPtr,
65 vtkImageData *outData, int *outExt,
66 double *outPtr, int id, vtkInformation *inInfo)
67 {
68 double r0, r1, r2, *r;
69 // For looping though output (and input) pixels.
70 int min0, max0, min1, max1, min2, max2;
71 int outIdx0, outIdx1, outIdx2;
72 vtkIdType outInc0, outInc1, outInc2;
73 double *outPtr0, *outPtr1, *outPtr2, *outPtrV;
74 vtkIdType inInc0, inInc1, inInc2;
75 T *inPtr0, *inPtr1, *inPtr2;
76 // For sobel function convolution (Left Right incs for each axis)
77 vtkIdType inInc0L, inInc0R, inInc1L, inInc1R, inInc2L, inInc2R;
78 T *inPtrL, *inPtrR;
79 double sum;
80 // Boundary of input image
81 int inWholeMin0, inWholeMax0, inWholeMin1, inWholeMax1;
82 int inWholeMin2, inWholeMax2;
83 int inWholeExt[6];
84 unsigned long count = 0;
85 unsigned long target;
86
87 // Get boundary information
88 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExt);
89 inWholeMin0 = inWholeExt[0];
90 inWholeMax0 = inWholeExt[1];
91 inWholeMin1 = inWholeExt[2];
92 inWholeMax1 = inWholeExt[3];
93 inWholeMin2 = inWholeExt[4];
94 inWholeMax2 = inWholeExt[5];
95
96 // Get information to march through data (skip component)
97 inData->GetIncrements(inInc0, inInc1, inInc2);
98 outData->GetIncrements(outInc0, outInc1, outInc2);
99 min0 = outExt[0]; max0 = outExt[1];
100 min1 = outExt[2]; max1 = outExt[3];
101 min2 = outExt[4]; max2 = outExt[5];
102
103 // We want the input pixel to correspond to output
104 inPtr = static_cast<T *>(inData->GetScalarPointer(min0,min1,min2));
105
106 // The data spacing is important for computing the gradient.
107 // Scale so it has the same range as gradient.
108 r = inData->GetSpacing();
109 r0 = 0.060445 / r[0];
110 r1 = 0.060445 / r[1];
111 r2 = 0.060445 / r[2];
112
113 target = static_cast<unsigned long>((max2-min2+1)*(max1-min1+1)/50.0);
114 target++;
115
116 // loop through pixels of output
117 outPtr2 = outPtr;
118 inPtr2 = inPtr;
119 for (outIdx2 = min2; outIdx2 <= max2; ++outIdx2)
120 {
121 inInc2L = (outIdx2 == inWholeMin2) ? 0 : -inInc2;
122 inInc2R = (outIdx2 == inWholeMax2) ? 0 : inInc2;
123
124 outPtr1 = outPtr2;
125 inPtr1 = inPtr2;
126 for (outIdx1 = min1; !self->AbortExecute && outIdx1 <= max1; ++outIdx1)
127 {
128 if (!id)
129 {
130 if (!(count%target))
131 {
132 self->UpdateProgress(count/(50.0*target));
133 }
134 count++;
135 }
136 inInc1L = (outIdx1 == inWholeMin1) ? 0 : -inInc1;
137 inInc1R = (outIdx1 == inWholeMax1) ? 0 : inInc1;
138
139 outPtr0 = outPtr1;
140 inPtr0 = inPtr1;
141 for (outIdx0 = min0; outIdx0 <= max0; ++outIdx0)
142 {
143 inInc0L = (outIdx0 == inWholeMin0) ? 0 : -inInc0;
144 inInc0R = (outIdx0 == inWholeMax0) ? 0 : inInc0;
145
146 // compute vector.
147 outPtrV = outPtr0;
148 // 12 Plane
149 inPtrL = inPtr0 + inInc0L;
150 inPtrR = inPtr0 + inInc0R;
151 sum = 2.0 * (*inPtrR - *inPtrL);
152 sum += static_cast<double>(inPtrR[inInc1L] + inPtrR[inInc1R]
153 + inPtrR[inInc2L] + inPtrR[inInc2R]);
154 sum += static_cast<double>(0.586 * (inPtrR[inInc1L+inInc2L] + inPtrR[inInc1L+inInc2R]
155 + inPtrR[inInc1R+inInc2L] + inPtrR[inInc1R+inInc2R]));
156 sum -= static_cast<double>(inPtrL[inInc1L] + inPtrL[inInc1R]
157 + inPtrL[inInc2L] + inPtrL[inInc2R]);
158 sum -= static_cast<double>(0.586 * (inPtrL[inInc1L+inInc2L] + inPtrL[inInc1L+inInc2R]
159 + inPtrL[inInc1R+inInc2L] + inPtrL[inInc1R+inInc2R]));
160 *outPtrV = sum * r0;
161 ++outPtrV;
162 // 02 Plane
163 inPtrL = inPtr0 + inInc1L;
164 inPtrR = inPtr0 + inInc1R;
165 sum = 2.0 * (*inPtrR - *inPtrL);
166 sum += static_cast<double>(inPtrR[inInc0L] + inPtrR[inInc0R]
167 + inPtrR[inInc2L] + inPtrR[inInc2R]);
168 sum += static_cast<double>(0.586 * (inPtrR[inInc0L+inInc2L] + inPtrR[inInc0L+inInc2R]
169 + inPtrR[inInc0R+inInc2L] + inPtrR[inInc0R+inInc2R]));
170 sum -= static_cast<double>(inPtrL[inInc0L] + inPtrL[inInc0R]
171 + inPtrL[inInc2L] + inPtrL[inInc2R]);
172 sum -= static_cast<double>(0.586 * (inPtrL[inInc0L+inInc2L] + inPtrL[inInc0L+inInc2R]
173 + inPtrL[inInc0R+inInc2L] + inPtrL[inInc0R+inInc2R]));
174 *outPtrV = sum * r1;
175 ++outPtrV;
176 // 01 Plane
177 inPtrL = inPtr0 + inInc2L;
178 inPtrR = inPtr0 + inInc2R;
179 sum = 2.0 * (*inPtrR - *inPtrL);
180 sum += static_cast<double>(inPtrR[inInc0L] + inPtrR[inInc0R]
181 + inPtrR[inInc1L] + inPtrR[inInc1R]);
182 sum += static_cast<double>(0.586 * (inPtrR[inInc0L+inInc1L] + inPtrR[inInc0L+inInc1R]
183 + inPtrR[inInc0R+inInc1L] + inPtrR[inInc0R+inInc1R]));
184 sum -= static_cast<double>(inPtrL[inInc0L] + inPtrL[inInc0R]
185 + inPtrL[inInc1L] + inPtrL[inInc1R]);
186 sum -= static_cast<double>(0.586 * (inPtrL[inInc0L+inInc1L] + inPtrL[inInc0L+inInc1R]
187 + inPtrL[inInc0R+inInc1L] + inPtrL[inInc0R+inInc1R]));
188 *outPtrV = static_cast<double>(sum * r2);
189 ++outPtrV;
190
191 outPtr0 += outInc0;
192 inPtr0 += inInc0;
193 }
194 outPtr1 += outInc1;
195 inPtr1 += inInc1;
196 }
197 outPtr2 += outInc2;
198 inPtr2 += inInc2;
199 }
200 }
201
202 //----------------------------------------------------------------------------
203 // This method contains a switch statement that calls the correct
204 // templated function for the input Data type. The output Data
205 // must be of type double. This method does handle boundary conditions.
206 // The third axis is the component axis for the output.
ThreadedRequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * vtkNotUsed (outputVector),vtkImageData *** inData,vtkImageData ** outData,int outExt[6],int id)207 void vtkImageSobel3D::ThreadedRequestData(
208 vtkInformation *vtkNotUsed(request),
209 vtkInformationVector **inputVector,
210 vtkInformationVector *vtkNotUsed(outputVector),
211 vtkImageData ***inData,
212 vtkImageData **outData,
213 int outExt[6], int id)
214 {
215 void *inPtr, *outPtr;
216 int inExt[6], wholeExt[6];
217
218 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
219 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExt);
220 this->InternalRequestUpdateExtent(inExt, outExt, wholeExt);
221
222 inPtr = inData[0][0]->GetScalarPointerForExtent(inExt);
223 outPtr = outData[0]->GetScalarPointerForExtent(outExt);
224
225 // this filter cannot handle multi component input.
226 if (inData[0][0]->GetNumberOfScalarComponents() != 1)
227 {
228 vtkWarningMacro("Expecting input with only one compenent.\n");
229 }
230
231 // this filter expects that output is type double.
232 if (outData[0]->GetScalarType() != VTK_DOUBLE)
233 {
234 vtkErrorMacro(<< "Execute: output ScalarType, "
235 << vtkImageScalarTypeNameMacro(outData[0]->GetScalarType())
236 << ", must be double");
237 return;
238 }
239
240 switch (inData[0][0]->GetScalarType())
241 {
242 vtkTemplateMacro(
243 vtkImageSobel3DExecute(this, inData[0][0],
244 static_cast<VTK_TT *>(inPtr), outData[0], outExt,
245 static_cast<double *>(outPtr),id, inInfo));
246 default:
247 vtkErrorMacro(<< "Execute: Unknown ScalarType");
248 return;
249 }
250 }
251