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