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 _clitkVectorBSplineDecompositionImageFilter_txx
19 #define _clitkVectorBSplineDecompositionImageFilter_txx
20 #include "clitkVectorBSplineDecompositionImageFilter.h"
22 #include "itkImageRegionConstIteratorWithIndex.h"
23 #include "itkImageRegionIterator.h"
24 #include "itkProgressReporter.h"
25 #include "itkVector.h"
33 template <class TInputImage, class TOutputImage>
34 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
35 ::VectorBSplineDecompositionImageFilter()
39 m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
40 m_IteratorDirection = 0;
41 this->SetSplineOrder(SplineOrder);
46 * Standard "PrintSelf" method
48 template <class TInputImage, class TOutputImage>
50 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
53 itk::Indent indent) const
55 Superclass::PrintSelf( os, indent );
56 os << indent << "Spline Order: " << m_SplineOrder << std::endl;
61 template <class TInputImage, class TOutputImage>
63 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
64 ::DataToCoefficients1D()
67 // See Unser, 1993, Part II, Equation 2.5,
68 // or Unser, 1999, Box 2. for an explaination.
72 if (m_DataLength[m_IteratorDirection] == 1) { //Required by mirror boundaries
76 // Compute overall gain
77 for (int k = 0; k < m_NumberOfPoles; k++) {
78 // Note for cubic splines lambda = 6
79 c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
83 for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) {
87 // loop over all poles
88 for (int k = 0; k < m_NumberOfPoles; k++) {
89 // causal initialization
90 this->SetInitialCausalCoefficient(m_SplinePoles[k]);
92 for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) {
93 m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
96 // anticausal initialization
97 this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
98 // anticausal recursion
99 for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) {
100 m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
108 template <class TInputImage, class TOutputImage>
110 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
111 ::SetSplineOrder(unsigned int SplineOrder)
113 if (SplineOrder == m_SplineOrder) {
116 m_SplineOrder = SplineOrder;
123 template <class TInputImage, class TOutputImage>
125 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
128 /* See Unser, 1997. Part II, Table I for Pole values */
129 // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
131 switch (m_SplineOrder) {
134 m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
144 m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
148 m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
149 m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
153 m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
155 m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
159 // SplineOrder not implemented yet.
160 itk::ExceptionObject err(__FILE__, __LINE__);
161 err.SetLocation( ITK_LOCATION);
162 err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
169 template <class TInputImage, class TOutputImage>
171 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
172 ::SetInitialCausalCoefficient(double z)
174 /* begining InitialCausalCoefficient */
175 /* See Unser, 1999, Box 2 for explaination */
177 itk::Vector<double, VectorDimension> sum;
178 double zn, z2n, iz; //sum
179 unsigned long horizon;
181 /* this initialization corresponds to mirror boundaries */
182 horizon = m_DataLength[m_IteratorDirection];
184 if (m_Tolerance > 0.0) {
185 horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
187 if (horizon < m_DataLength[m_IteratorDirection]) {
188 /* accelerated loop */
189 sum = m_Scratch[0]; // verify this
190 for (unsigned int n = 1; n < horizon; n++) {
191 sum += zn * m_Scratch[n];
198 z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
199 sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
201 for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) {
202 sum += (zn + z2n) * m_Scratch[n];
206 m_Scratch[0] = sum / (1.0 - zn * zn);
211 template <class TInputImage, class TOutputImage>
213 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
214 ::SetInitialAntiCausalCoefficient(double z)
216 // this initialization corresponds to mirror boundaries
217 /* See Unser, 1999, Box 2 for explaination */
218 // Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
219 m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
220 (z / (z * z - 1.0)) *
221 (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
225 template <class TInputImage, class TOutputImage>
227 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
228 ::DataToCoefficientsND()
230 OutputImagePointer output = this->GetOutput();
232 itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
234 unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
236 itk::ProgressReporter progress(this, 0, count, 10);
238 // Initialize coeffient array
239 this->CopyImageToImage(); // Coefficients are initialized to the input data
241 for (unsigned int n=0; n < ImageDimension; n++) {
242 m_IteratorDirection = n;
243 // Loop through each dimension
245 // Initialize iterators
246 OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
247 CIterator.SetDirection( m_IteratorDirection );
248 // For each data vector
249 while ( !CIterator.IsAtEnd() ) {
250 // Copy coefficients to scratch
251 this->CopyCoefficientsToScratch( CIterator );
254 // Perform 1D BSpline calculations
255 this->DataToCoefficients1D();
257 // Copy scratch back to coefficients.
258 // Brings us back to the end of the line we were working on.
259 CIterator.GoToBeginOfLine();
260 this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
261 CIterator.NextLine();
262 progress.CompletedPixel();
269 * Copy the input image into the output image
271 template <class TInputImage, class TOutputImage>
273 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
277 typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
278 typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
279 typedef typename TOutputImage::PixelType OutputPixelType;
281 InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
282 OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
285 outIt = outIt.Begin();
287 while ( !outIt.IsAtEnd() ) {
288 for (unsigned int i=0; i< VectorDimension; i++) {
289 v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
300 * Copy the scratch to one line of the output image
302 template <class TInputImage, class TOutputImage>
304 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
305 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
307 typedef typename TOutputImage::PixelType OutputPixelType;
310 while ( !Iter.IsAtEndOfLine() ) {
311 for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
321 * Copy one line of the output image to the scratch
323 template <class TInputImage, class TOutputImage>
325 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
326 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
329 itk::Vector<double, VectorDimension> v;
330 while ( !Iter.IsAtEndOfLine() ) {
331 for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
340 * GenerateInputRequestedRegion method.
342 template <class TInputImage, class TOutputImage>
344 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
345 ::GenerateInputRequestedRegion()
347 // this filter requires the all of the input image to be in
349 InputImagePointer inputPtr = const_cast< TInputImage * > ( this->GetInput() );
351 inputPtr->SetRequestedRegionToLargestPossibleRegion();
357 * EnlargeOutputRequestedRegion method.
359 template <class TInputImage, class TOutputImage>
361 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
362 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
365 // this filter requires the all of the output image to be in
367 TOutputImage *imgData;
368 imgData = dynamic_cast<TOutputImage*>( output );
370 imgData->SetRequestedRegionToLargestPossibleRegion();
378 template <class TInputImage, class TOutputImage>
380 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
383 DD("VectorBSplineDecompositionImageFilter GenerateData()");
384 // Allocate scratch memory
385 InputImageConstPointer inputPtr = this->GetInput();
386 m_DataLength = inputPtr->GetBufferedRegion().GetSize();
388 unsigned long maxLength = 0;
389 for ( unsigned int n = 0; n < ImageDimension; n++ ) {
390 if ( m_DataLength[n] > maxLength ) {
391 maxLength = m_DataLength[n];
394 m_Scratch.resize( maxLength );
396 // Allocate memory for output image
397 OutputImagePointer outputPtr = this->GetOutput();
398 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
399 outputPtr->Allocate();
401 // Calculate actual output
402 this->DataToCoefficientsND();