forked from InsightSoftwareConsortium/ITK
-
Notifications
You must be signed in to change notification settings - Fork 0
/
NeighborhoodIterators5.cxx
217 lines (187 loc) · 7.77 KB
/
NeighborhoodIterators5.cxx
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
216
217
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: NeighborhoodIterators5.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 "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkGaussianOperator.h"
#include "itkNeighborhoodInnerProduct.h"
// Software Guide : BeginLatex
//
// This example introduces slice-based neighborhood processing. A slice, in
// this context, is a 1D path through an ND neighborhood. Slices are defined
// for generic arrays by the \code{std::slice} class as a start index, a step
// size, and an end index. Slices simplify the implementation of certain
// neighborhood calculations. They also provide a mechanism for taking inner
// products with subregions of neighborhoods.
//
// Suppose, for example, that we want to take partial derivatives in the $y$
// direction of a neighborhood, but offset those derivatives by one pixel
// position along the positive $x$ direction. For a $3\times3$, 2D
// neighborhood iterator, we can construct an \code{std::slice}, \code{(start =
// 2, stride = 3, end = 8)}, that represents the neighborhood offsets $(1,
// -1)$, $(1, 0)$, $(1, 1)$ (see Figure~\ref{fig:NeighborhoodIteratorFig2}). If we
// pass this slice as an extra argument to the
// \doxygen{NeighborhoodInnerProduct} function, then the inner product is taken
// only along that slice. This ``sliced'' inner product with a 1D
// \doxygen{DerivativeOperator} gives the desired derivative.
//
// The previous separable Gaussian filtering example can be rewritten using
// slices and slice-based inner products. In general, slice-based processing
// is most useful when doing many different calculations on the same
// neighborhood, where defining multiple iterators as in
// Section~\ref{sec:NeighborhoodExample4} becomes impractical or inefficient.
// Good examples of slice-based neighborhood processing can be found in any of
// the ND anisotropic diffusion function objects, such as
// \doxygen{CurvatureNDAnisotropicDiffusionFunction}.
//
// Software Guide : EndLatex
int main( int argc, char ** argv )
{
if ( argc < 4 )
{
std::cerr << "Missing parameters. " << std::endl;
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0]
<< " inputImageFile outputImageFile sigma"
<< std::endl;
return -1;
}
typedef float PixelType;
typedef itk::Image< PixelType, 2 > ImageType;
typedef itk::ImageFileReader< ImageType > ReaderType;
typedef itk::ConstNeighborhoodIterator< ImageType > NeighborhoodIteratorType;
typedef itk::ImageRegionIterator< ImageType> IteratorType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( argv[1] );
try
{
reader->Update();
}
catch ( itk::ExceptionObject &err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
ImageType::Pointer output = ImageType::New();
output->SetRegions(reader->GetOutput()->GetRequestedRegion());
output->Allocate();
itk::NeighborhoodInnerProduct<ImageType> innerProduct;
typedef itk::NeighborhoodAlgorithm
::ImageBoundaryFacesCalculator< ImageType > FaceCalculatorType;
FaceCalculatorType faceCalculator;
FaceCalculatorType::FaceListType faceList;
FaceCalculatorType::FaceListType::iterator fit;
IteratorType out;
NeighborhoodIteratorType it;
// Software Guide: BeginLatex
//
// The first difference between this example and the previous example is that
// the Gaussian operator is only initialized once. Its direction is not
// important because it is only a 1D array of coefficients.
//
// Software Guide: EndLatex
// Software Guide : BeginCodeSnippet
itk::GaussianOperator< PixelType, 2 > gaussianOperator;
gaussianOperator.SetDirection(0);
gaussianOperator.SetVariance( ::atof(argv[3]) * ::atof(argv[3]) );
gaussianOperator.CreateDirectional();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Next we need to define a radius for the iterator. The radius in all
// directions matches that of the single extent of the Gaussian operator,
// defining a square neighborhood.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
NeighborhoodIteratorType::RadiusType radius;
radius.Fill( gaussianOperator.GetRadius()[0] );
// Software Guide EndCodeSnippet
// Software Guide : BeginLatex
//
// The inner product and face calculator are defined for the main processing
// loop as before, but now the iterator is reinitialized each iteration with
// the square \code{radius} instead of the radius of the operator. The
// inner product is taken using a slice along the axial direction corresponding
// to the current iteration. Note the use of \code{GetSlice()} to return the
// proper slice from the iterator itself. \code{GetSlice()} can only be used
// to return the slice along the complete extent of the axial direction of a
// neighborhood.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ImageType::Pointer input = reader->GetOutput();
faceList = faceCalculator(input, output->GetRequestedRegion(), radius);
for (unsigned int i = 0; i < ImageType::ImageDimension; ++i)
{
for ( fit=faceList.begin(); fit != faceList.end(); ++fit )
{
it = NeighborhoodIteratorType( radius, input, *fit );
out = IteratorType( output, *fit );
for (it.GoToBegin(), out.GoToBegin(); ! it.IsAtEnd(); ++it, ++out)
{
out.Set( innerProduct(it.GetSlice(i), it, gaussianOperator) );
}
}
// Swap the input and output buffers
if (i != ImageType::ImageDimension - 1)
{
ImageType::Pointer tmp = input;
input = output;
output = tmp;
}
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// This technique produces exactly the same results as the previous example. A
// little experimentation, however, will reveal that it is less efficient since
// the neighborhood iterator is keeping track of extra, unused pixel locations
// for each iteration, while the previous example only references those pixels
// that it needs. In cases, however, where an algorithm takes multiple
// derivatives or convolution products over the same neighborhood, slice-based
// processing can increase efficiency and simplify the implementation.
//
// Software Guide : EndLatex
typedef unsigned char WritePixelType;
typedef itk::Image< WritePixelType, 2 > WriteImageType;
typedef itk::ImageFileWriter< WriteImageType > WriterType;
typedef itk::RescaleIntensityImageFilter< ImageType,
WriteImageType > RescaleFilterType;
RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
rescaler->SetOutputMinimum( 0 );
rescaler->SetOutputMaximum( 255 );
rescaler->SetInput(output);
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( argv[2] );
writer->SetInput( rescaler->GetOutput() );
try
{
writer->Update();
}
catch ( itk::ExceptionObject &err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
return 0;
}