Orfeo Toolbox  4.0
itkCumulativeGaussianOptimizer.cxx
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef _itkCumulativeGaussianOptimizer_cxx
19 #define _itkCumulativeGaussianOptimizer_cxx
20 
22 #include "itkMath.h"
23 
24 namespace itk
25 {
27 {
28  // Set some initial values for the variables.
29  m_ComputedMean = 0;
33  m_UpperAsymptote = 0;
34  m_LowerAsymptote = 0;
35  m_OffsetForMean = 0;
37  m_Verbose = 0;
38  m_FitError = 0;
41  m_StopConditionDescription << this->GetNameOfClass() << ": Constructed";
42 }
43 
45 {
46  delete m_FinalSampledArray;
47 }
48 
51 ::ExtendGaussian(MeasureType *originalArray, MeasureType *extendedArray, int startingPointForInsertion)
52 {
53  // Use the parameters from originalArray to construct a Gaussian in
54  // extendedArray
55  // shifting the mean to the right by startingPointForInsertion.
56  double mean = startingPointForInsertion + m_ComputedMean;
57  double sd = m_ComputedStandardDeviation;
58  double amplitude = m_ComputedAmplitude;
59 
60  m_OffsetForMean = startingPointForInsertion;
61 
62  for ( int i = 0; i < (int)( extendedArray->GetNumberOfElements() ); i++ )
63  {
64  extendedArray->put( i, amplitude * vcl_exp( -( vcl_pow( ( i - mean ), 2 ) / ( 2 * vcl_pow(sd, 2) ) ) ) );
65  }
66  // Then insert the originalArray over the middle section of extendedArray.
67  for ( int i = 0; i < (int)( originalArray->GetNumberOfElements() ); i++ )
68  {
69  extendedArray->put( i + startingPointForInsertion, originalArray->get(i) );
70  }
71  return extendedArray;
72 }
73 
74 double
77 {
78  // Given two arrays array1 and array2 of equal length, calculate the average
79  // sum of squared
80  // differences between them.
81  int size = array1->GetNumberOfElements();
82  double sum = 0;
83 
84  for ( int i = 0; i < size; i++ )
85  {
86  sum = sum + ( array1->get(i) - array2->get(i) ) * ( array1->get(i) - array2->get(i) );
87  }
88  return ( sum / size );
89 }
90 
91 void
94 {
95  // Measure the parameters of the sampled Gaussian curve and use these
96  // parameters to
97  // construct an extended curve. Then measure the parameters of the extended
98  // curve, which
99  // should be closer to the original Gaussian parameters and recalculate the
100  // extended portion
101  // of the curve. Iterate these last two steps until the average sum of squared
102  // differences
103  // between 2 iterations converge within differenceTolerance.
104  MeasureGaussianParameters(sampledGaussianArray);
105 
106  if ( m_Verbose )
107  {
108  PrintComputedParameterHeader();
109  PrintComputedParameters();
110  }
111 
112  int sampledGaussianArraySize = sampledGaussianArray->GetNumberOfElements();
113  int extendedArraySize = 3 * sampledGaussianArraySize;
114  MeasureType *extendedArray = new MeasureType();
115  extendedArray->SetSize(extendedArraySize);
116  MeasureType *extendedArrayCopy = new MeasureType();
117  extendedArrayCopy->SetSize(extendedArraySize);
118 
119  double averageSumOfSquaredDifferences = m_DifferenceTolerance;
120 
121  extendedArray = ExtendGaussian(sampledGaussianArray, extendedArray, sampledGaussianArraySize);
122 
123  MeasureGaussianParameters(extendedArray);
124  bool smallChangeBetweenIterations = false;
125  while ( averageSumOfSquaredDifferences >= m_DifferenceTolerance )
126  {
127  for ( int j = 0; j < extendedArraySize; j++ )
128  {
129  extendedArrayCopy->put( j, extendedArray->get(j) );
130  }
131  extendedArray = RecalculateExtendedArrayFromGaussianParameters(sampledGaussianArray,
132  extendedArray,
133  sampledGaussianArraySize);
134 
135  MeasureGaussianParameters(extendedArray);
136  if ( m_Verbose )
137  {
138  PrintComputedParameters();
139  }
140  double temp = averageSumOfSquaredDifferences;
141  averageSumOfSquaredDifferences = FindAverageSumOfSquaredDifferences(extendedArray, extendedArrayCopy);
142 
143  // Stop if there is a very very very small change between iterations.
144  if ( vcl_fabs(temp - averageSumOfSquaredDifferences) <= m_DifferenceTolerance )
145  {
146  m_StopConditionDescription.str("");
147  m_StopConditionDescription << this->GetNameOfClass() << ": "
148  << "Change between iterations ("
149  << vcl_fabs(temp - averageSumOfSquaredDifferences)
150  << ") is less than DifferenceTolerance ("
151  << m_DifferenceTolerance
152  << ").";
153  break;
154  }
155  }
156  if ( !smallChangeBetweenIterations )
157  {
158  m_StopConditionDescription.str("");
159  m_StopConditionDescription << this->GetNameOfClass() << ": "
160  << "Average sum of squared differences ("
161  << averageSumOfSquaredDifferences
162  << ") is less than DifferenceTolerance ("
163  << m_DifferenceTolerance
164  << ").";
165  }
166 
167  // Update the mean calculation.
168  m_ComputedMean = m_ComputedMean - m_OffsetForMean;
169 
170  delete extendedArray;
171  delete extendedArrayCopy;
172 }
173 
176 {
177  // Assuming the input array is Gaussian, compute the mean, SD, amplitude, and
178  // change in intensity.
179  m_ComputedMean = 0;
180  m_ComputedStandardDeviation = 0;
181  m_ComputedAmplitude = 0;
182  m_ComputedTransitionHeight = 0;
183 
184  double sum = 0;
185 
186  // Calculate the mean.
187  for ( int i = 0; i < (int)( array->GetNumberOfElements() ); i++ )
188  {
189  m_ComputedMean += i * array->get(i);
190  sum += array->get(i);
191  }
192  // Assertion fails if number of samples <=2 or UpperAsymptote==LowerAsymptote
193  // improper behavior if number of samples == 3.
194  itkAssertInDebugAndIgnoreInReleaseMacro(sum != 0);
195  m_ComputedMean /= sum;
196 
197  // Calculate the standard deviation
198  for ( int i = 0; i < (int)( array->GetNumberOfElements() ); i++ )
199  {
200  m_ComputedStandardDeviation += array->get(i) * vcl_pow( ( i - m_ComputedMean ), 2 );
201  }
202  m_ComputedStandardDeviation = vcl_sqrt(m_ComputedStandardDeviation / sum);
203 
204  // For the ERF, sum is the difference between the lower and upper intensities.
205  m_ComputedTransitionHeight = sum;
206 
207  // Calculate the amplitude.
208  m_ComputedAmplitude = sum / ( m_ComputedStandardDeviation * vcl_sqrt(2 * vnl_math::pi) );
209 }
210 
211 void
214 {
215  std::cerr << "Mean\t" << "SD\t" << "Amp\t" << "Transition" << std::endl;
216 }
217 
220 {
221  std::cerr << m_ComputedMean - m_OffsetForMean << "\t" // Printed mean is
222  // shifted.
223  << m_ComputedStandardDeviation << "\t"
224  << m_ComputedAmplitude << "\t"
225  << m_ComputedTransitionHeight << std::endl;
226 }
227 
231  MeasureType *extendedArray,
232  int startingPointForInsertion)
233 {
234  // From the Gaussian parameters stored with the extendedArray,
235  // recalculate the extended portion of the extendedArray,
236  // leaving the inserted original array unchaged.
237  double mean = m_ComputedMean;
238  double sd = m_ComputedStandardDeviation;
239  double amplitude = m_ComputedAmplitude;
240 
241  for ( int i = 0; i < (int)( extendedArray->GetNumberOfElements() ); i++ )
242  {
243  // Leave the original inserted array unchanged.
244  if ( i < startingPointForInsertion
245  || i >= startingPointForInsertion + (int)( originalArray->GetNumberOfElements() ) )
246  {
247  extendedArray->put( i, amplitude * vcl_exp( -( vcl_pow( ( i - mean ), 2 ) / ( 2 * vcl_pow(sd, 2) ) ) ) );
248  }
249  }
250  return extendedArray;
251 }
252 
253 void
255 ::SetDataArray(MeasureType *cumGaussianArray)
256 {
257  m_CumulativeGaussianArray = cumGaussianArray;
258 }
259 
260 void
263 {
264  this->InvokeEvent( StartEvent() );
265  m_StopConditionDescription.str("");
266  m_StopConditionDescription << this->GetNameOfClass() << ": Running";
267 
268  // Declare arrays.
269  int cumGaussianArraySize = m_CumulativeGaussianArray->GetNumberOfElements();
270  int sampledGaussianArraySize = cumGaussianArraySize;
271  // int cumGaussianArrayCopySize = cumGaussianArraySize;
272 
273  MeasureType *sampledGaussianArray = new MeasureType();
274  sampledGaussianArray->SetSize(sampledGaussianArraySize);
275 
276  MeasureType *cumGaussianArrayCopy = new MeasureType();
277  cumGaussianArrayCopy->SetSize(cumGaussianArraySize);
278 
279  // Make a copy of the Cumulative Gaussian sampled data array.
280  for ( int j = 0; j < cumGaussianArraySize; j++ )
281  {
282  cumGaussianArrayCopy->put( j, m_CumulativeGaussianArray->get(j) );
283  }
284  // Take the derivative of the data array resulting in a Gaussian array.
285  MeasureType *derivative = new MeasureType();
286  derivative->SetSize(cumGaussianArraySize - 1);
287 
288  for ( int i = 1; i < (int)( derivative->GetNumberOfElements() + 1 ); i++ )
289  {
290  derivative->put( i - 1, m_CumulativeGaussianArray->get(i) - m_CumulativeGaussianArray->get(i - 1) );
291  }
292  m_CumulativeGaussianArray = derivative;
293 
294  // Iteratively recalculate and resample the Gaussian array.
295  FindParametersOfGaussian(m_CumulativeGaussianArray);
296 
297  // Generate new Gaussian array with final parameters.
298  for ( int i = 0; i < sampledGaussianArraySize; i++ )
299  {
300  sampledGaussianArray->put( i, m_ComputedAmplitude
301  * vcl_exp( -( vcl_pow( ( i - m_ComputedMean ),
302  2 ) / ( 2 * vcl_pow(m_ComputedStandardDeviation, 2) ) ) ) );
303  }
304  // Add 0.5 to the mean of the sampled Gaussian curve to make up for the 0.5
305  // shift during derivation, then take the integral of the Gaussian sample
306  // to produce a Cumulative Gaussian.
307 
308  for ( int i = sampledGaussianArraySize - 1; i > 0; i-- )
309  {
310  sampledGaussianArray->put( i - 1, sampledGaussianArray->get(i) - sampledGaussianArray->get(i - 1) );
311  }
312  m_ComputedMean += 0.5;
313 
314  // Find the best vertical shift that minimizes the least square error.
315  double c = VerticalBestShift(cumGaussianArrayCopy, sampledGaussianArray);
316 
317  // Add constant c to array.
318  for ( int i = 0; i < (int)( sampledGaussianArray->GetNumberOfElements() ); i++ )
319  {
320  sampledGaussianArray->put(i, sampledGaussianArray->get(i) + c);
321  }
322  // Calculate the mean, standard deviation, lower and upper asymptotes of the
323  // sampled Cumulative Gaussian.
324  int floorOfMean = (int)( m_ComputedMean );
325  double yFloorOfMean = sampledGaussianArray->get(floorOfMean);
326  double yCeilingOfMean = sampledGaussianArray->get(floorOfMean + 1);
327  double y = ( m_ComputedMean - floorOfMean ) * ( yCeilingOfMean - yFloorOfMean ) + yFloorOfMean;
328  m_UpperAsymptote = y + m_ComputedTransitionHeight / 2;
329  m_LowerAsymptote = y - m_ComputedTransitionHeight / 2;
330 
331  m_FinalSampledArray = new MeasureType();
332  m_FinalSampledArray->SetSize( sampledGaussianArray->GetNumberOfElements() );
333  for ( int i = 0; i < (int)( m_FinalSampledArray->GetNumberOfElements() ); i++ )
334  {
335  m_FinalSampledArray->put( i, sampledGaussianArray->get(i) );
336  }
337  // Calculate the least square error as a measure of goodness of fit.
338  m_FitError = static_cast< CostFunctionType * >( m_CostFunction.GetPointer() )->CalculateFitError(sampledGaussianArray);
339 
340  delete sampledGaussianArray;
341  delete cumGaussianArrayCopy;
342  delete derivative;
343 }
344 
346 {
347  for ( int i = 0; i < (int)( array->GetNumberOfElements() ); i++ )
348  {
349  std::cerr << i << " " << array->get(i) << std::endl;
350  }
351 }
352 
353 double
355 ::VerticalBestShift(MeasureType *originalArray, MeasureType *newArray)
356 {
357  // Find the constant to minimize the sum of squares of the difference between
358  // original Array and newArray+c
359  // Proof of algorithm:
360  // Let A = the original array.
361  // Let B = the new array.
362  // Let n = the number of elements in each array (note they must be the
363  // same).
364  // We want to mimimize sum(((Bi+c) - (Ai))^2).
365  // So we take the derivative with respect to c and equate this derivative
366  // to 0.
367  // d/dc sum(((Bi+c) - (Ai))^2) dc = 0
368  // => sum (2(Bi+c - Ai)) = 0
369  // => sum (Bi + c - Ai) = 0
370  // => (sum(Bi)) + (sum(c)) - (sum(Ai)) = 0
371  // => nC = sum(Ai) - sum(Bi)
372  // => C = (sum(Ai) - sum(Bi)) / n
373 
374  double c = 0;
375  int size = originalArray->GetNumberOfElements();
376 
377  for ( int i = 0; i < size; i++ )
378  {
379  c += originalArray->get(i);
380  }
381  for ( int i = 0; i < size; i++ )
382  {
383  c -= newArray->get(i);
384  }
385  return ( c / size );
386 }
387 
388 const std::string
391 {
392  return m_StopConditionDescription.str();
393 }
394 
395 void
397 ::PrintSelf(std::ostream & os, Indent indent) const
398 {
399  Superclass::PrintSelf(os, indent);
400 
401  os << indent << "Difference Tolerance = " << m_DifferenceTolerance << std::endl;
402  os << indent << "Computed Mean = " << m_ComputedMean << std::endl;
403  os << indent << "Computed Standard Deviation = " << m_ComputedStandardDeviation << std::endl;
404  os << indent << "Computed Amplitude = " << m_ComputedAmplitude << std::endl;
405  os << indent << "Computed Transition Height = " << m_ComputedTransitionHeight << std::endl;
406 
407  os << indent << "Upper Asymptote = " << m_UpperAsymptote << std::endl;
408  os << indent << "Lower Asymptote = " << m_LowerAsymptote << std::endl;
409  os << indent << "Offset For Mean = " << m_OffsetForMean << std::endl;
410  os << indent << "Verbose = " << m_Verbose << std::endl;
411  os << indent << "Fit Error = " << m_FitError << std::endl;
412 
413  os << indent << "StopConditionDescription: " << m_StopConditionDescription.str() << std::endl;
414  if ( m_FinalSampledArray )
415  {
416  os << indent << "Final Sampled Array = " << m_FinalSampledArray << std::endl;
417  }
418  else
419  {
420  os << indent << "Final Sampled Array = [not defined] " << std::endl;
421  }
422 }
423 } // end namespace itk
424 
425 #endif

Generated at Sat Mar 8 2014 14:31:02 for Orfeo Toolbox with doxygen 1.8.3.1