1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18 #ifndef _itkBSplineDecompositionImageFilterWithOBD_txx
19 #define _itkBSplineDecompositionImageFilterWithOBD_txx
20 #include "itkBSplineDecompositionImageFilterWithOBD.h"
21 #include "itkImageRegionConstIteratorWithIndex.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkProgressReporter.h"
24 #include "itkVector.h"
32 template <class TInputImage, class TOutputImage>
33 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
34 ::BSplineDecompositionImageFilterWithOBD()
37 for(unsigned int i=0; i<ImageDimension; i++)
40 m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
41 m_IteratorDirection = 0;
42 this->SetSplineOrder(SplineOrder);
47 * Standard "PrintSelf" method
49 template <class TInputImage, class TOutputImage>
51 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
56 Superclass::PrintSelf( os, indent );
57 os << indent << "Spline Order: " << m_SplineOrder << std::endl;
62 template <class TInputImage, class TOutputImage>
64 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
65 ::DataToCoefficients1D()
68 // See Unser, 1993, Part II, Equation 2.5,
69 // or Unser, 1999, Box 2. for an explaination.
73 if (m_DataLength[m_IteratorDirection] == 1) { //Required by mirror boundaries
77 // Compute overall gain
78 for (int k = 0; k < m_NumberOfPoles; k++) {
79 // Note for cubic splines lambda = 6
80 c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
84 for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) {
88 // loop over all poles
89 for (int k = 0; k < m_NumberOfPoles; k++) {
90 // causal initialization
91 this->SetInitialCausalCoefficient(m_SplinePoles[k]);
93 for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) {
94 m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
97 // anticausal initialization
98 this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
99 // anticausal recursion
100 for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) {
101 m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
109 template <class TInputImage, class TOutputImage>
111 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
112 ::SetSplineOrder(unsigned int SplineOrder)
114 if (SplineOrder == m_SplineOrder) {
117 m_SplineOrder = SplineOrder;
124 template <class TInputImage, class TOutputImage>
126 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
127 ::SetSplineOrders(SizeType SplineOrders)
129 if (SplineOrders == m_SplineOrders)
133 m_SplineOrders=SplineOrders;
137 template <class TInputImage, class TOutputImage>
139 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
142 /* See Unser, 1997. Part II, Table I for Pole values */
143 // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
145 switch (m_SplineOrder) {
149 m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
159 m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
163 m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
164 m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
168 m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
170 m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
174 // SplineOrder not implemented yet.
175 ExceptionObject err(__FILE__, __LINE__);
176 err.SetLocation( ITK_LOCATION);
177 err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
184 template <class TInputImage, class TOutputImage>
186 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
187 ::SetInitialCausalCoefficient(double z)
189 /* begining InitialCausalCoefficient */
190 /* See Unser, 1999, Box 2 for explaination */
192 double sum, zn, z2n, iz;
193 unsigned long horizon;
195 /* this initialization corresponds to mirror boundaries */
196 horizon = m_DataLength[m_IteratorDirection];
198 if (m_Tolerance > 0.0) {
199 horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
201 if (horizon < m_DataLength[m_IteratorDirection]) {
202 /* accelerated loop */
203 sum = m_Scratch[0]; // verify this
204 for (unsigned int n = 1; n < horizon; n++) {
205 sum += zn * m_Scratch[n];
212 z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
213 sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
215 for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) {
216 sum += (zn + z2n) * m_Scratch[n];
220 m_Scratch[0] = sum / (1.0 - zn * zn);
225 template <class TInputImage, class TOutputImage>
227 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
228 ::SetInitialAntiCausalCoefficient(double z)
230 // this initialization corresponds to mirror boundaries
231 /* See Unser, 1999, Box 2 for explaination */
232 // Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
233 m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
234 (z / (z * z - 1.0)) *
235 (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
239 template <class TInputImage, class TOutputImage>
241 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
242 ::DataToCoefficientsND()
244 OutputImagePointer output = this->GetOutput();
246 Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
248 unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
250 ProgressReporter progress(this, 0, count, 10);
252 // Initialize coeffient array
253 this->CopyImageToImage(); // Coefficients are initialized to the input data
255 for (unsigned int n=0; n < ImageDimension; n++) {
256 m_IteratorDirection = n;
257 // Loop through each dimension
259 //JV Set the correct order by dimension!
260 SetSplineOrder(m_SplineOrders[n]);
262 // Initialize iterators
263 OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
264 CIterator.SetDirection( m_IteratorDirection );
265 // For each data vector
266 while ( !CIterator.IsAtEnd() ) {
267 // Copy coefficients to scratch
268 this->CopyCoefficientsToScratch( CIterator );
270 // Perform 1D BSpline calculations
271 this->DataToCoefficients1D();
273 // Copy scratch back to coefficients.
274 // Brings us back to the end of the line we were working on.
275 CIterator.GoToBeginOfLine();
276 this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
277 CIterator.NextLine();
278 progress.CompletedPixel();
285 * Copy the input image into the output image
287 template <class TInputImage, class TOutputImage>
289 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
293 typedef ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
294 typedef ImageRegionIterator< TOutputImage > OutputIterator;
295 typedef typename TOutputImage::PixelType OutputPixelType;
297 InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
298 OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
301 outIt = outIt.Begin();
303 while ( !outIt.IsAtEnd() ) {
304 outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
313 * Copy the scratch to one line of the output image
315 template <class TInputImage, class TOutputImage>
317 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
318 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
320 typedef typename TOutputImage::PixelType OutputPixelType;
322 while ( !Iter.IsAtEndOfLine() ) {
323 Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
332 * Copy one line of the output image to the scratch
334 template <class TInputImage, class TOutputImage>
336 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
337 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
340 while ( !Iter.IsAtEndOfLine() ) {
341 m_Scratch[j] = static_cast<double>( Iter.Get() ) ;
349 * GenerateInputRequestedRegion method.
351 template <class TInputImage, class TOutputImage>
353 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
354 ::GenerateInputRequestedRegion()
356 // this filter requires the all of the input image to be in
358 InputImagePointer inputPtr = const_cast< TInputImage * > ( this->GetInput() );
360 inputPtr->SetRequestedRegionToLargestPossibleRegion();
366 * EnlargeOutputRequestedRegion method.
368 template <class TInputImage, class TOutputImage>
370 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
371 ::EnlargeOutputRequestedRegion(
375 // this filter requires the all of the output image to be in
377 TOutputImage *imgData;
378 imgData = dynamic_cast<TOutputImage*>( output );
380 imgData->SetRequestedRegionToLargestPossibleRegion();
388 template <class TInputImage, class TOutputImage>
390 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
394 // Allocate scratch memory
395 InputImageConstPointer inputPtr = this->GetInput();
396 m_DataLength = inputPtr->GetBufferedRegion().GetSize();
398 unsigned long maxLength = 0;
399 for ( unsigned int n = 0; n < ImageDimension; n++ ) {
400 if ( m_DataLength[n] > maxLength ) {
401 maxLength = m_DataLength[n];
404 m_Scratch.resize( maxLength );
406 // Allocate memory for output image
407 OutputImagePointer outputPtr = this->GetOutput();
408 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
409 outputPtr->Allocate();
411 // Calculate actual output
412 this->DataToCoefficientsND();