1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkShepardMethod.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 "vtkShepardMethod.h"
16
17 #include "vtkFloatArray.h"
18 #include "vtkImageData.h"
19 #include "vtkMath.h"
20 #include "vtkInformation.h"
21 #include "vtkInformationVector.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkStreamingDemandDrivenPipeline.h"
24 #include "vtkPointData.h"
25
26 vtkStandardNewMacro(vtkShepardMethod);
27
28 // Construct with sample dimensions=(50,50,50) and so that model bounds are
29 // automatically computed from input. Null value for each unvisited output
30 // point is 0.0. Maximum distance is 0.25.
vtkShepardMethod()31 vtkShepardMethod::vtkShepardMethod()
32 {
33 this->MaximumDistance = 0.25;
34
35 this->ModelBounds[0] = 0.0;
36 this->ModelBounds[1] = 0.0;
37 this->ModelBounds[2] = 0.0;
38 this->ModelBounds[3] = 0.0;
39 this->ModelBounds[4] = 0.0;
40 this->ModelBounds[5] = 0.0;
41
42 this->SampleDimensions[0] = 50;
43 this->SampleDimensions[1] = 50;
44 this->SampleDimensions[2] = 50;
45
46 this->NullValue = 0.0;
47 }
48
49 // Compute ModelBounds from input geometry.
ComputeModelBounds(double origin[3],double spacing[3])50 double vtkShepardMethod::ComputeModelBounds(double origin[3],
51 double spacing[3])
52 {
53 double *bounds, maxDist;
54 int i, adjustBounds=0;
55
56 // compute model bounds if not set previously
57 if ( this->ModelBounds[0] >= this->ModelBounds[1] ||
58 this->ModelBounds[2] >= this->ModelBounds[3] ||
59 this->ModelBounds[4] >= this->ModelBounds[5] )
60 {
61 adjustBounds = 1;
62 vtkDataSet *ds = vtkDataSet::SafeDownCast(this->GetInput());
63 // ds better be non null otherwise something is very wrong here
64 bounds = ds->GetBounds();
65 }
66 else
67 {
68 bounds = this->ModelBounds;
69 }
70
71 for (maxDist=0.0, i=0; i<3; i++)
72 {
73 if ( (bounds[2*i+1] - bounds[2*i]) > maxDist )
74 {
75 maxDist = bounds[2*i+1] - bounds[2*i];
76 }
77 }
78 maxDist *= this->MaximumDistance;
79
80 // adjust bounds so model fits strictly inside (only if not set previously)
81 if ( adjustBounds )
82 {
83 for (i=0; i<3; i++)
84 {
85 this->ModelBounds[2*i] = bounds[2*i] - maxDist;
86 this->ModelBounds[2*i+1] = bounds[2*i+1] + maxDist;
87 }
88 }
89
90 // Set volume origin and data spacing
91 for (i=0; i<3; i++)
92 {
93 origin[i] = this->ModelBounds[2*i];
94 spacing[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i])
95 / (this->SampleDimensions[i] - 1);
96 }
97
98 return maxDist;
99 }
100
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)101 int vtkShepardMethod::RequestInformation (
102 vtkInformation * vtkNotUsed(request),
103 vtkInformationVector ** vtkNotUsed( inputVector ),
104 vtkInformationVector *outputVector)
105 {
106 // get the info objects
107 vtkInformation* outInfo = outputVector->GetInformationObject(0);
108
109 int i;
110 double ar[3], origin[3];
111
112 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
113 0, this->SampleDimensions[0]-1,
114 0, this->SampleDimensions[1]-1,
115 0, this->SampleDimensions[2]-1);
116
117 for (i=0; i < 3; i++)
118 {
119 origin[i] = this->ModelBounds[2*i];
120 if ( this->SampleDimensions[i] <= 1 )
121 {
122 ar[i] = 1;
123 }
124 else
125 {
126 ar[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i])
127 / (this->SampleDimensions[i] - 1);
128 }
129 }
130 outInfo->Set(vtkDataObject::ORIGIN(),origin,3);
131 outInfo->Set(vtkDataObject::SPACING(),ar,3);
132
133 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_FLOAT, 1);
134 return 1;
135 }
136
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)137 int vtkShepardMethod::RequestData(
138 vtkInformation* vtkNotUsed( request ),
139 vtkInformationVector** inputVector,
140 vtkInformationVector* outputVector)
141 {
142 // get the input
143 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
144 vtkDataSet *input = vtkDataSet::SafeDownCast(
145 inInfo->Get(vtkDataObject::DATA_OBJECT()));
146
147 // get the output
148 vtkInformation *outInfo = outputVector->GetInformationObject(0);
149 vtkImageData *output = vtkImageData::SafeDownCast(
150 outInfo->Get(vtkDataObject::DATA_OBJECT()));
151
152 // We need to allocate our own scalars since we are overriding
153 // the superclasses "Execute()" method.
154 output->SetExtent(
155 outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
156 output->AllocateScalars(outInfo);
157
158 vtkIdType ptId, i;
159 int j, k;
160 double *px, x[3], s, *sum, spacing[3], origin[3];
161
162 double maxDistance, distance2, inScalar;
163 vtkDataArray *inScalars;
164 vtkIdType numPts, numNewPts, idx;
165 int min[3], max[3];
166 int jkFactor;
167 vtkFloatArray *newScalars =
168 vtkFloatArray::SafeDownCast(output->GetPointData()->GetScalars());
169
170 vtkDebugMacro(<< "Executing Shepard method");
171
172 // Check input
173 //
174 if ( (numPts=input->GetNumberOfPoints()) < 1 )
175 {
176 vtkErrorMacro(<<"Points must be defined!");
177 return 1;
178 }
179
180 if ( (inScalars = input->GetPointData()->GetScalars()) == NULL )
181 {
182 vtkErrorMacro(<<"Scalars must be defined!");
183 return 1;
184 }
185
186 newScalars->SetName(inScalars->GetName());
187
188 // Allocate
189 //
190 numNewPts = this->SampleDimensions[0] * this->SampleDimensions[1]
191 * this->SampleDimensions[2];
192
193 sum = new double[numNewPts];
194 for (i=0; i<numNewPts; i++)
195 {
196 newScalars->SetComponent(i,0,0.0);
197 sum[i] = 0.0;
198 }
199
200 maxDistance = this->ComputeModelBounds(origin,spacing);
201 outInfo->Set(vtkDataObject::ORIGIN(),origin,3);
202 outInfo->Set(vtkDataObject::SPACING(),spacing,3);
203
204
205 // Traverse all input points.
206 // Each input point affects voxels within maxDistance.
207 //
208 for (ptId=0; ptId < numPts; ptId++)
209 {
210 if ( ! (ptId % 1000) )
211 {
212 vtkDebugMacro(<<"Inserting point #" << ptId);
213 this->UpdateProgress (ptId/numPts);
214 if (this->GetAbortExecute())
215 {
216 break;
217 }
218 }
219
220 px = input->GetPoint(ptId);
221 inScalar = inScalars->GetComponent(ptId,0);
222
223 for (i=0; i<3; i++) //compute dimensional bounds in data set
224 {
225 double amin = static_cast<double>(
226 (px[i] - maxDistance) - origin[i]) / spacing[i];
227 double amax = static_cast<double>(
228 (px[i] + maxDistance) - origin[i]) / spacing[i];
229 min[i] = static_cast<int>(amin);
230 max[i] = static_cast<int>(amax);
231
232 if (min[i] < amin)
233 {
234 min[i]++; // round upward to nearest integer to get min[i]
235 }
236 if (max[i] > amax)
237 {
238 max[i]--; // round downward to nearest integer to get max[i]
239 }
240
241 if (min[i] < 0)
242 {
243 min[i] = 0; // valid range check
244 }
245 if (max[i] >= this->SampleDimensions[i])
246 {
247 max[i] = this->SampleDimensions[i] - 1;
248 }
249 }
250
251 for (i=0; i<3; i++) //compute dimensional bounds in data set
252 {
253 min[i] = static_cast<int>(
254 static_cast<double>((px[i] - maxDistance) - origin[i]) / spacing[i]);
255 max[i] = static_cast<int>(
256 static_cast<double>((px[i] + maxDistance) - origin[i]) / spacing[i]);
257 if (min[i] < 0)
258 {
259 min[i] = 0;
260 }
261 if (max[i] >= this->SampleDimensions[i])
262 {
263 max[i] = this->SampleDimensions[i] - 1;
264 }
265 }
266
267 jkFactor = this->SampleDimensions[0]*this->SampleDimensions[1];
268 for (k = min[2]; k <= max[2]; k++)
269 {
270 x[2] = spacing[2] * k + origin[2];
271 for (j = min[1]; j <= max[1]; j++)
272 {
273 x[1] = spacing[1] * j + origin[1];
274 for (i = min[0]; i <= max[0]; i++)
275 {
276 x[0] = spacing[0] * i + origin[0];
277 idx = jkFactor*k + this->SampleDimensions[0]*j + i;
278
279 distance2 = vtkMath::Distance2BetweenPoints(x,px);
280
281 if ( distance2 == 0.0 )
282 {
283 sum[idx] = VTK_DOUBLE_MAX;
284 newScalars->SetComponent(idx,0,VTK_FLOAT_MAX);
285 }
286 else
287 {
288 s = newScalars->GetComponent(idx,0);
289 sum[idx] += 1.0 / distance2;
290 newScalars->SetComponent(idx,0,s+(inScalar/distance2));
291 }
292 }
293 }
294 }
295 }
296
297 // Run through scalars and compute final values
298 //
299 for (ptId=0; ptId<numNewPts; ptId++)
300 {
301 s = newScalars->GetComponent(ptId,0);
302 if ( sum[ptId] != 0.0 )
303 {
304 newScalars->SetComponent(ptId,0,s/sum[ptId]);
305 }
306 else
307 {
308 newScalars->SetComponent(ptId,0,this->NullValue);
309 }
310 }
311
312 // Update self
313 //
314 delete [] sum;
315
316 return 1;
317 }
318
319 // Set the i-j-k dimensions on which to sample the distance function.
SetSampleDimensions(int i,int j,int k)320 void vtkShepardMethod::SetSampleDimensions(int i, int j, int k)
321 {
322 int dim[3];
323
324 dim[0] = i;
325 dim[1] = j;
326 dim[2] = k;
327
328 this->SetSampleDimensions(dim);
329 }
330
331 // Set the i-j-k dimensions on which to sample the distance function.
SetSampleDimensions(int dim[3])332 void vtkShepardMethod::SetSampleDimensions(int dim[3])
333 {
334 int dataDim, i;
335
336 vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << ","
337 << dim[1] << "," << dim[2] << ")");
338
339 if ( dim[0] != this->SampleDimensions[0] ||
340 dim[1] != this->SampleDimensions[1] ||
341 dim[2] != this->SampleDimensions[2] )
342 {
343 if ( dim[0]<1 || dim[1]<1 || dim[2]<1 )
344 {
345 vtkErrorMacro (<< "Bad Sample Dimensions, retaining previous values");
346 return;
347 }
348
349 for (dataDim=0, i=0; i<3 ; i++)
350 {
351 if (dim[i] > 1)
352 {
353 dataDim++;
354 }
355 }
356
357 if ( dataDim < 3 )
358 {
359 vtkErrorMacro(<<"Sample dimensions must define a volume!");
360 return;
361 }
362
363 for ( i=0; i<3; i++)
364 {
365 this->SampleDimensions[i] = dim[i];
366 }
367
368 this->Modified();
369 }
370 }
371
FillInputPortInformation(int vtkNotUsed (port),vtkInformation * info)372 int vtkShepardMethod::FillInputPortInformation(
373 int vtkNotUsed(port), vtkInformation* info)
374 {
375 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
376 return 1;
377 }
378
PrintSelf(ostream & os,vtkIndent indent)379 void vtkShepardMethod::PrintSelf(ostream& os, vtkIndent indent)
380 {
381 this->Superclass::PrintSelf(os,indent);
382
383 os << indent << "Maximum Distance: " << this->MaximumDistance << "\n";
384
385 os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", "
386 << this->SampleDimensions[1] << ", "
387 << this->SampleDimensions[2] << ")\n";
388
389 os << indent << "ModelBounds: \n";
390 os << indent << " Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n";
391 os << indent << " Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n";
392 os << indent << " Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n";
393
394 os << indent << "Null Value: " << this->NullValue << "\n";
395
396 }
397