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
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: SampleToHistogramProjectionFilter.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
// Software Guide : BeginLatex
//
// \index{Statistics!Projecting measurement vectors to 1-D histogram}
// \index{itk::Statistics::Sample\-To\-Histogram\-Projection\-Filter}
//
// The \subdoxygen{Statistics}{SampleToHistogramProjectionFilter} projects
// measurement vectors of a sample onto a vector and fills up a 1-D
// \subdoxygen{Statistics}{Histogram}. The histogram will be formed around the
// mean value set by the \code{SetMean()} method. The histogram's measurement
// values are the distance between the mean and the projected measurement
// vectors normalized by the standard deviation set by the
// \code{SetStandardDeviation()} method. Such histogram can be used to
// analyze the multi-dimensional distribution or examine the
// \emph{goodness-of-fit} of a projected distribution (histogram) with its
// expected distribution.
//
// We will use the ListSample as the input sample.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkListSample.h"
#include "itkSampleToHistogramProjectionFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We need another header for measurement vectors. We are going to use
// the \doxygen{Vector} class which is a subclass of the \doxygen{FixedArray}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkVector.h"
// Software Guide : EndCodeSnippet
int main()
{
// Software Guide : BeginLatex
//
// The following code snippet will create a ListSample object
// with two-component int measurement vectors and put the measurement
// vectors: [1,1] - 1 time, [2,2] - 2 times, [3,3] - 3 times, [4,4] -
// 4 times, [5,5] - 5 times into the \code{listSample}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int MeasurementVectorLength = 2;
typedef int MeasurementType;
typedef itk::Vector< MeasurementType , MeasurementVectorLength > MeasurementVectorType;
typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;
SampleType::Pointer sample = SampleType::New();
sample->SetMeasurementVectorSize( MeasurementVectorLength );
MeasurementVectorType mv;
for ( unsigned int i = 1 ; i < 6 ; i++ )
{
for ( unsigned int j = 0 ; j < 2 ; j++ )
{
mv[j] = ( MeasurementType ) i;
}
for ( unsigned int j = 0 ; j < i ; j++ )
{
sample->PushBack(mv);
}
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We create a histogram that has six bins. The histogram's range is
// [-2, 2). Since the \code{sample} has measurement vectors between
// [1, 1] and [5,5], The histogram does not seem to cover the whole
// range. However, the SampleToHistogramProjectionFilter
// normalizes the measurement vectors with the given mean and the
// standard deviation. Therefore, the projected value is approximately
// the distance between the measurement vector and the mean divided by
// the standard deviation.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::Statistics::Histogram< float, 1 > HistogramType;
HistogramType::Pointer histogram = HistogramType::New();
HistogramType::SizeType size;
size.Fill(6);
HistogramType::MeasurementVectorType lowerBound;
HistogramType::MeasurementVectorType upperBound;
lowerBound[0] = -2;
upperBound[0] = 2;
histogram->Initialize( size, lowerBound, upperBound );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
// We use the \code{SetInputSample(sample*)} and the
// \code{SetHistogram(histogram*)} methods to set the input
// sample and the output histogram that have been created.
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::Statistics::SampleToHistogramProjectionFilter<SampleType, float>
ProjectorType;
ProjectorType::Pointer projector = ProjectorType::New();
projector->SetInputSample( sample );
projector->SetHistogram( histogram );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// As mentioned above, this class projects measurement vectors onto the
// projection axis with normalization using the mean and standard
// deviation.
// \begin{equation}
// y = \frac{\sum^{d}_{i=0} (x_{i} - \mu_{i})\alpha_{i}}{\sigma}
// \end{equation}
// where, $y$ is the projected value, $x$ is the $i$th component of the
// measurement vector, $\mu_{i}$ is the $i$th component of the mean vector,
// $\alpha_{i}$ is the $i$th component of the projection axis (a
// vector), and $\sigma$ is the standard deviation.
//
// If the bin overlap value is set by the \code{SetHistogramBinOverlap()}
// method and it is greater than 0.001, the frequency will be weighted based
// on its closeness of the bin boundaries. In other words, even if a
// measurement vector falls into a bin, depending on its closeness to the
// adjacent bins, the frequencies of the adjacent bins will be also updated
// with weights. If we do not want to use the bin overlapping function, we do
// not call the \code{SetHistogramBinOverlap(double)} method. The default
// value for the histogram bin overlap is zero, so without calling the
// method, the filter will not use bin overlapping \cite{Aylward1997a}
// \cite{Aylward1997b}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ProjectorType::MeanType mean( MeasurementVectorLength );
mean[0] = 3.66667;
mean[1] = 3.66667;
double standardDeviation = 3;
ProjectorType::ArrayType projectionAxis( MeasurementVectorLength );
projectionAxis[0] = 1;
projectionAxis[1] = 1;
projector->SetMean( &mean );
projector->SetStandardDeviation( &standardDeviation );
projector->SetProjectionAxis( &projectionAxis );
projector->SetHistogramBinOverlap( 0.25 );
projector->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We print out the updated histogram after the projection.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
float fSum = 0.0;
HistogramType::Iterator iter = histogram->Begin();
while ( iter != histogram->End() )
{
std::cout << "instance identifier = " << iter.GetInstanceIdentifier()
<< "\t measurement vector = "
<< iter.GetMeasurementVector()
<< "\t frequency = "
<< iter.GetFrequency() << std::endl;
fSum += iter.GetFrequency();
++iter;
}
std::cout << " sum of frequency = " << fSum << std::endl;
// Software Guide : EndCodeSnippet
return 0;
}
|
