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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "itkPolyLineParametricPath.h"
#include "itkChainCodePath.h"
#include "itkFourierSeriesPath.h"
#include "itkPathToChainCodePathFilter.h"
#include "itkChainCodeToFourierSeriesPathFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkExtractOrthogonalSwath2DImageFilter.h"
#include "itkDerivativeImageFilter.h"
#include "itkOrthogonalSwath2DPathFilter.h"
#include "itkPathToImageFilter.h"
int itkOrthogonalSwath2DPathFilterTest(int, char*[])
{
typedef itk::Image<unsigned char, 2> ImageType;
typedef itk::Image<double, 2> FloatImageType;
typedef itk::PolyLineParametricPath<2> InPathType;
typedef itk::ChainCodePath<2> ChainPathType;
typedef itk::FourierSeriesPath<2> FSPathType;
typedef InPathType::VertexType VertexType;
typedef ImageType::IndexType IndexType;
typedef itk::OrthogonallyCorrected2DParametricPath OutputPathType;
// pre-process the image
typedef itk::RescaleIntensityImageFilter<ImageType,FloatImageType> CastFilterType;
typedef itk::DiscreteGaussianImageFilter<FloatImageType,FloatImageType> SmoothFilterType;
typedef itk::ExtractOrthogonalSwath2DImageFilter<FloatImageType> SwathFilterType;
typedef itk::DerivativeImageFilter<FloatImageType,FloatImageType> MeritFilterType;
// pre-process the path
typedef itk::PathToChainCodePathFilter<InPathType,ChainPathType> PathFilter1Type;
typedef itk::ChainCodeToFourierSeriesPathFilter<ChainPathType,FSPathType> PathFilter2Type;
// test the filter
typedef itk::OrthogonalSwath2DPathFilter<FSPathType,FloatImageType> TestFilterType;
// save the results
typedef itk::PathToImageFilter< FSPathType, ImageType > Output1FilterType;
typedef itk::PathToImageFilter< OutputPathType, ImageType > Output2FilterType;
typedef itk::RescaleIntensityImageFilter<FloatImageType,ImageType> Output3FilterType;
// Setup the path
std::cout << "Making a square Path with v0 at (24,24) -> (24,104) -> (104,104) -> (104,24)" << std::endl;
InPathType::Pointer inPath;
VertexType v;
inPath = InPathType::New();
v.Fill(24);
inPath->AddVertex(v);
v[0]=24;
v[1]=104;
inPath->AddVertex(v);
v.Fill(104);
inPath->AddVertex(v);
v[0]=104;
v[1]=24;
inPath->AddVertex(v);
v.Fill(24);
inPath->AddVertex(v);
// Setup the first path filter
PathFilter1Type::Pointer pathFilter1 = PathFilter1Type::New();
pathFilter1->SetInput(inPath);
// Setup the second path filter
PathFilter2Type::Pointer pathFilter2 = PathFilter2Type::New();
pathFilter2->SetInput(pathFilter1->GetOutput());
pathFilter2->SetNumberOfHarmonics(7); // make a nice, round, path for the swath
// Setup the image
std::cout << "Making a 64x64 black square centered in a 128x128 white image"<<std::endl;
ImageType::Pointer inImage = ImageType::New();
IndexType start;
start[0]=0;
start[1]=0;
ImageType::SizeType size;
size[0]=128;
size[1]=128;
ImageType::RegionType region;
region.SetSize(size);
region.SetIndex(start);
inImage->SetRegions(region);
double spacing[ ImageType::ImageDimension ];
spacing[0]=1.0;
spacing[1]=1.0;
inImage->SetSpacing(spacing);
inImage->Allocate();
typedef itk::ImageRegionIterator<ImageType> ImageItType;
ImageItType it( inImage, inImage->GetRequestedRegion() );
it.GoToBegin();
IndexType pixelIndex;
while( !it.IsAtEnd() )
{
pixelIndex = it.GetIndex();
if( pixelIndex[0] >= int(size[0]/4) && pixelIndex[0] < int(size[0]*3/4) &&
pixelIndex[1] >= int(size[1]/4) && pixelIndex[1] < int(size[1]*3/4) )
it.Set(0);
else
it.Set(255);
++it;
}
// Cast the input image into a double image
CastFilterType::Pointer castFilter = CastFilterType::New();
castFilter->SetInput( inImage );
castFilter->SetOutputMinimum(0);
castFilter->SetOutputMaximum(1.0);
// Smooth the (double pixel type) input image
SmoothFilterType::Pointer smoothFilter = SmoothFilterType::New();
smoothFilter->SetInput( castFilter->GetOutput() );
double gaussianVariance = 1.0;
// We want a fast 3x3 kernel. Gausian Operator will not truncate its kernel
// width to any less than a 5x5 kernel (kernel width of 3 for 1 center pixel +
// 2 edge pixels). However, GausianOperator always uses at least a 3x3
// kernel, and so setting the maximum error to 1.0 (no limit) will make it
// stop growing the kernel at the desired 3x3 size.
double maxError = 0.9;
smoothFilter->SetUseImageSpacingOff();
smoothFilter->SetVariance( gaussianVariance );
smoothFilter->SetMaximumError( maxError );
// Extract the swath image
SwathFilterType::Pointer swathFilter = SwathFilterType::New();
swathFilter->SetImageInput( smoothFilter->GetOutput() );
// swathFilter->SetImageInput( castFilter->GetOutput() );
swathFilter->SetPathInput( pathFilter2->GetOutput() );
size[0]=512;
size[1]=16*2+1; // the top 1 and bottom 1 rows are dropped when smoothing
swathFilter->SetSize(size);
// Find the vertical gradient of the swath image
MeritFilterType::Pointer meritFilter = MeritFilterType::New();
meritFilter->SetInput( swathFilter->GetOutput() );
meritFilter->SetOrder( 1 ); // first partial derivative
meritFilter->SetDirection( 1 ); // d/dy
// Setup the test filter
std::cerr << "Creating the test filter" << std::endl;
TestFilterType::Pointer testFilter = TestFilterType::New();
std::cerr << "Setting up the test filter" << std::endl;
testFilter->SetPathInput(pathFilter2->GetOutput());
testFilter->SetImageInput(meritFilter->GetOutput());
OutputPathType::Pointer outPath = testFilter->GetOutput();
// Setup the input & output path images
std::cerr << "Setting up the output path image" << std::endl;
Output1FilterType::Pointer output1Filter = Output1FilterType::New();
Output2FilterType::Pointer output2Filter = Output2FilterType::New();
size[0]=128;
size[1]=128;
output1Filter->SetSize(size); // same size as the input image
output2Filter->SetSize(size); // same size as the input image
output1Filter->SetPathValue(255);
output2Filter->SetPathValue(255);
output1Filter->SetInput( testFilter->GetPathInput() );
output2Filter->SetInput( testFilter->GetOutput() );
ImageType::Pointer inPathImage = output1Filter->GetOutput();
ImageType::Pointer outImage = output2Filter->GetOutput();
// Setup the swath merit output image
Output3FilterType::Pointer output3Filter = Output3FilterType::New();
output3Filter->SetInput( meritFilter->GetOutput() );
output3Filter->SetOutputMinimum(0);
output3Filter->SetOutputMaximum(255);
ImageType::Pointer swathMeritImage = output3Filter->GetOutput();
// Testing PrintSelf
std::cerr << testFilter << std::endl;
// Update the pipeline
std::cerr << "Running the Pipeline: ";
outImage->Update();
std::cerr << "[DONE]" << std::endl;
/*
// Save the output images
itk::ImageFileWriter<ImageType>::Pointer writer
= itk::ImageFileWriter<ImageType>::New();
writer->SetInput( output1Filter->GetOutput() );
writer->SetFileName( "OrthSwathFilterPathIn.png" );
writer->Write();
writer->SetInput( output2Filter->GetOutput() );
writer->SetFileName( "OrthSwathFilterPathOut.png" );
writer->Write();
writer->SetInput( output3Filter->GetOutput() );
writer->SetFileName( "OrthSwathFilterSwathMerit.png" );
writer->Write();
*/
return EXIT_SUCCESS;
}
|
