1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: TestImageDataToStructuredGridFilter.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 "vtkUniformGrid.h"
16 #include "vtkStructuredGrid.h"
17 #include "vtkCell.h"
18 #include "vtkMath.h"
19 #include "vtkPointData.h"
20 #include "vtkCellData.h"
21 #include "vtkDataArray.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkImageToStructuredGrid.h"
24
25 #include <cmath>
26 #include <limits>
27 #include <string>
28 #include <cstring>
29
30 // Description:
31 // Performs safe division a/b which also checks for underflow & overflow
SafeDivision(const double a,const double b)32 double SafeDivision( const double a, const double b )
33 {
34 // Catch overflow
35 if( (b < 1) && (a > ( b*std::numeric_limits<double>::max() ) ) )
36 {
37 return std::numeric_limits<double>::max();
38 }
39
40 // Catch underflow
41 if( (a == static_cast<double>(0.0)) ||
42 ( (b > 1) && (a < b*std::numeric_limits<double>::max() ) ) )
43 {
44 return( static_cast<double>( 0.0 ) );
45 }
46
47 return( a/b );
48 }
49
50 // Description:
51 // Checks if two given floating numbers are equivalent.
52 // This algorithm is based on Knuth, The of Computer Programming (vol II).
FloatNumberEquals(double a,double b,double TOL)53 bool FloatNumberEquals( double a, double b, double TOL )
54 {
55 double adiff = std::abs( a-b );
56 double d1 = SafeDivision( adiff,std::abs(a) );
57 double d2 = SafeDivision( adiff,std::abs(b) );
58 if( (d1 <= TOL) || (d2 <= TOL) )
59 {
60 return true;
61 }
62 return false;
63 }
64
65 // Description:
66 // Checks if the given two points a & b are equal
ArePointsEqual(double a[3],double b[3],double TOL=1e-9)67 bool ArePointsEqual( double a[3], double b[3], double TOL=1e-9 )
68 {
69 for( int i=0; i < 3; ++i )
70 {
71 if( !FloatNumberEquals( a[i], b[i], TOL ) )
72 {
73 return false;
74 }
75 }
76 return true;
77 }
78
79 // Description:
80 // Constructs a uniform grid instance with the given spacing & dimensions
81 // at the user-supplied origin.
GetGrid(double * origin,double * spacing,int * ndim)82 vtkUniformGrid* GetGrid( double *origin, double *spacing, int *ndim )
83 {
84 vtkUniformGrid *grd = vtkUniformGrid::New();
85 grd->Initialize();
86 grd->SetOrigin( origin );
87 grd->SetSpacing( spacing );
88 grd->SetDimensions( ndim );
89
90 vtkDoubleArray* pntData = vtkDoubleArray::New();
91 pntData->SetName( "XYZ-NODE" );
92 pntData->SetNumberOfComponents( 1 );
93 pntData->SetNumberOfTuples( grd->GetNumberOfPoints() );
94
95 double node[3];
96 for( int pntIdx=0; pntIdx < grd->GetNumberOfPoints(); ++pntIdx )
97 {
98 grd->GetPoint( pntIdx, node );
99 pntData->SetValue(pntIdx, (node[0]+node[1]+node[2]) );
100 } // END for all points
101 grd->GetPointData()->AddArray( pntData );
102 pntData->Delete();
103
104
105 vtkDoubleArray* xyz = vtkDoubleArray::New( );
106 xyz->SetName( "XYZ-CELL" );
107 xyz->SetNumberOfComponents( 1 );
108 xyz->SetNumberOfTuples( grd->GetNumberOfCells() );
109
110 for( int cellIdx=0; cellIdx < grd->GetNumberOfCells(); ++cellIdx )
111 {
112 vtkCell* myCell = grd->GetCell( cellIdx);
113
114 vtkPoints *cellPoints = myCell->GetPoints();
115
116 double xyzCenter[3];
117 xyzCenter[0] = 0.0;
118 xyzCenter[1] = 0.0;
119 xyzCenter[2] = 0.0;
120
121 for( int cp=0; cp < cellPoints->GetNumberOfPoints(); ++cp )
122 {
123 double pnt[3];
124 cellPoints->GetPoint( cp, pnt );
125 xyzCenter[0] += pnt[0];
126 xyzCenter[1] += pnt[1];
127 xyzCenter[2] += pnt[2];
128 } // END for all cell points
129
130 xyzCenter[0] = xyzCenter[0] / (cellPoints->GetNumberOfPoints());
131 xyzCenter[1] = xyzCenter[1] / (cellPoints->GetNumberOfPoints());
132 xyzCenter[2] = xyzCenter[2] / (cellPoints->GetNumberOfPoints());
133
134 double f = xyzCenter[0]*xyzCenter[0] +
135 xyzCenter[1]*xyzCenter[1] + xyzCenter[2]*xyzCenter[2];
136 xyz->SetTuple1(cellIdx,f);
137 } // END for all cells
138
139 grd->GetCellData()->AddArray(xyz);
140 xyz->Delete();
141 return grd;
142 }
143
144 // Description
145 // Checks if the given image data-set is equivalent to the structured grid
146 // data-set.
DataSetsEqual(vtkImageData * img,vtkStructuredGrid * sg)147 bool DataSetsEqual( vtkImageData* img, vtkStructuredGrid* sg )
148 {
149 // 0. Check dimensions
150 int imgdim[3]; int sgdim[3];
151 img->GetDimensions( imgdim );
152 sg->GetDimensions( sgdim );
153 for( int i=0; i < 3; ++i )
154 {
155 if( imgdim[i] != sgdim[i] )
156 {
157 return false;
158 }
159 }
160
161 // 1. Check Number of elements
162 if( img->GetNumberOfCells() != sg->GetNumberOfCells() )
163 {
164 return false;
165 }
166
167 // 2. Check Number of points
168 if( img->GetNumberOfPoints() != sg->GetNumberOfPoints() )
169 {
170 return false;
171 }
172
173 // 3. Check Point equality
174 double pnt1[3]; double pnt2[3];
175 for( int pntIdx=0; pntIdx < img->GetNumberOfPoints(); ++pntIdx )
176 {
177 img->GetPoint( pntIdx, pnt1 );
178 sg->GetPoint( pntIdx,pnt2 );
179 if( !ArePointsEqual( pnt1, pnt2 ) )
180 {
181 return false;
182 }
183 }
184
185 // 4. Check Point data equality
186 if( img->GetPointData()->GetNumberOfArrays() !=
187 sg->GetPointData()->GetNumberOfArrays() )
188 {
189 return false;
190 }
191 int dataArrayIdx = 0;
192 for(;dataArrayIdx<img->GetPointData()->GetNumberOfArrays();++dataArrayIdx )
193 {
194 vtkDataArray* array1 = img->GetPointData()->GetArray( dataArrayIdx );
195 vtkDataArray* array2 = sg->GetPointData()->GetArray( dataArrayIdx );
196 if( array1->GetNumberOfComponents() != array2->GetNumberOfComponents() )
197 {
198 return false;
199 }
200 if( array1->GetNumberOfTuples() != array2->GetNumberOfTuples() )
201 {
202 return false;
203 }
204 if( std::strcmp( array1->GetName(), array2->GetName() ) != 0 )
205 {
206 return false;
207 }
208 } // END for all point arrays
209
210 // 5. Check Cell data equality
211 if( img->GetCellData()->GetNumberOfArrays() !=
212 sg->GetCellData()->GetNumberOfArrays() )
213 {
214 return false;
215 }
216
217 dataArrayIdx = 0;
218 for(;dataArrayIdx<img->GetCellData()->GetNumberOfArrays();++dataArrayIdx)
219 {
220 vtkDataArray* array1 = img->GetCellData()->GetArray( dataArrayIdx );
221 vtkDataArray* array2 = sg->GetCellData()->GetArray( dataArrayIdx );
222 if( array1->GetNumberOfComponents() != array2->GetNumberOfComponents() )
223 {
224 return false;
225 }
226 if( array1->GetNumberOfTuples() != array2->GetNumberOfTuples() )
227 {
228 return false;
229 }
230 if( std::strcmp( array1->GetName(), array2->GetName() ) != 0 )
231 {
232 return false;
233 }
234 } // END for all cell arrays
235
236 return true;
237 }
238
TestImageDataToStructuredGrid(int,char * [])239 int TestImageDataToStructuredGrid(int,char *[])
240 {
241 int rval = 0;
242
243 double origin[3] = {0.0,0.0,0.0};
244 double spacing[3] = {0.5,0.2,0.0};
245 int ndim[3] = {10, 10, 1};
246 vtkUniformGrid* img1 = GetGrid( origin, spacing, ndim );
247
248 vtkImageToStructuredGrid* myFilter = vtkImageToStructuredGrid::New( );
249 myFilter->SetInputData( img1 );
250 img1->Delete();
251 myFilter->Update();
252 vtkStructuredGrid* sg1 = myFilter->GetOutput();
253
254 if( !DataSetsEqual( img1, sg1 ) )
255 {
256 rval = 1;
257 }
258
259 myFilter->Delete();
260 return( rval );
261 }
262