1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: itkWarpMeshFilterTest.cxx
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkWarpMeshFilter.h"
#include "itkRegularSphereMeshSource.h"
#include "itkMesh.h"
#include "itkImage.h"
int itkWarpMeshFilterTest(int, char* [] )
{
const unsigned int Dimension = 3;
// Declare the type of the input and output mesh
typedef itk::DefaultStaticMeshTraits<
double, Dimension, Dimension, double, double, double> MeshTraits;
typedef itk::Mesh<double,Dimension,MeshTraits> MeshType;
// declare triangle mesh source
typedef itk::RegularSphereMeshSource< MeshType > SphereMeshSourceType;
typedef SphereMeshSourceType::PointType PointType;
typedef SphereMeshSourceType::VectorType VectorType;
SphereMeshSourceType::Pointer sphereMeshSource = SphereMeshSourceType::New();
PointType center;
center[0] = 25.0;
center[1] = 25.0;
center[2] = 25.0;
VectorType scale;
scale[0] = 20.0;
scale[1] = 20.0;
scale[2] = 20.0;
sphereMeshSource->SetCenter(center);
sphereMeshSource->SetScale(scale);
sphereMeshSource->SetResolution(2);
// Declare the deformation field
typedef itk::Vector< double, Dimension> VectorType;
typedef itk::Image< VectorType, Dimension > DeformationFieldType;
DeformationFieldType::Pointer deformationField = DeformationFieldType::New();
DeformationFieldType::IndexType start;
start[0] = 0;
start[1] = 0;
start[2] = 0;
DeformationFieldType::SizeType size;
size[0] = 25;
size[1] = 25;
size[2] = 25;
DeformationFieldType::RegionType region;
region.SetSize( size );
region.SetIndex( start );
deformationField->SetRegions( region );
DeformationFieldType::PointType origin;
origin[0] = 0.0;
origin[1] = 0.0;
origin[2] = 0.0;
deformationField->SetOrigin( origin );
DeformationFieldType::SpacingType spacing;
spacing[0] = 2.0;
spacing[1] = 2.0;
spacing[2] = 2.0;
deformationField->SetSpacing( spacing );
deformationField->Allocate();
VectorType simpleVector; // set up an uniform vector field
simpleVector[0] = 3.0;
simpleVector[1] = 5.0;
simpleVector[2] = 7.0;
deformationField->FillBuffer( simpleVector );
// Declare the Warping filter
typedef itk::WarpMeshFilter<
MeshType, MeshType, DeformationFieldType >
WarpFilterType;
WarpFilterType::Pointer warpFilter = WarpFilterType::New();
warpFilter->SetInput( sphereMeshSource->GetOutput() );
warpFilter->SetDeformationField( deformationField );
try
{
warpFilter->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception: " << excp << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
MeshType::Pointer outputMesh = warpFilter->GetOutput();
MeshType::ConstPointer inputMesh = warpFilter->GetInput();
MeshType::PointsContainerPointer outPoints = outputMesh->GetPoints();
MeshType::PointsContainerConstPointer inPoints = inputMesh->GetPoints();
MeshType::PointsContainer::ConstIterator inputPoint = inPoints->Begin();
MeshType::PointsContainer::ConstIterator outputPoint = outPoints->Begin();
MeshType::PointsContainer::ConstIterator lastInputPoint = inPoints->End();
MeshType::PointsContainer::ConstIterator lastOutputPoint = outPoints->End();
const double tolerance = 1e-8;
while( inputPoint != lastInputPoint && outputPoint != lastOutputPoint )
{
for(unsigned int i=0; i<Dimension; i++)
{
const double distance = outputPoint.Value()[i] - inputPoint.Value()[i];
if( vcl_fabs( distance - simpleVector[i] ) > tolerance )
{
std::cerr << "Filter failed" << std::endl;
std::cerr << "Expected displacement = " << simpleVector[i] << std::endl;
std::cerr << "but got = " << distance << std::endl;
return EXIT_FAILURE;
}
}
++inputPoint;
++outputPoint;
}
return EXIT_SUCCESS;
}
|
