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 __clitkConeBeamProjectImageFilter_txx
19 #define __clitkConeBeamProjectImageFilter_txx
24 //=========================================================================================================================
26 template <class InputImageType, class OutputImageType>
27 ConeBeamProjectImageFilter<InputImageType, OutputImageType>::
28 ConeBeamProjectImageFilter()
30 // Initialize with the default values (Elekta Synergie)
32 m_IsInitialized=false;
33 m_NumberOfThreadsIsGiven=false;
36 m_IsoCenter.Fill(0.0);
37 m_SourceToScreen=1536.;
40 m_RigidTransformMatrix.SetIdentity();
41 m_EdgePaddingValue=itk::NumericTraits<OutputPixelType>::Zero;//density images
44 m_Transform=TransformType::New();
45 m_Resampler=ResampleImageFilterType::New();
46 m_Interpolator= InterpolatorType::New();
47 m_ExtractImageFilter=ExtractImageFilterType::New();
48 m_FlipImageFilter=FlipImageFilterType::New();
50 // Parameters for output
51 this->m_OutputSpacing.Fill(0.8);
52 this->m_OutputOrigin.Fill(0.0);
53 this->m_OutputSize.Fill( 512 );
58 //-----------------------------------------------------------------------
59 // Set output info from image
60 //-----------------------------------------------------------------------
61 template <class InputImageType, class OutputImageType>
63 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
64 ::SetOutputParametersFromImage ( OutputImagePointer image )
66 this->SetOutputOrigin( image->GetOrigin() );
67 this->SetOutputSpacing( image->GetSpacing() );
68 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
72 //-----------------------------------------------------------------------
73 // Set output info from image
74 //-----------------------------------------------------------------------
75 template <class InputImageType, class OutputImageType>
77 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
78 ::SetOutputParametersFromImage ( OutputImageConstPointer image )
80 this->SetOutputOrigin( image->GetOrigin() );
81 this->SetOutputSpacing( image->GetSpacing() );
82 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
86 //=========================================================================================================================
88 template <class InputImageType, class OutputImageType>
90 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
91 ::Initialize(void) throw (itk::ExceptionObject)
94 //====================================================================
95 // Define the transform: composition of rigid transfo around origin and a centered rotation
96 // The center of rotation should be placed at the isocenter!
98 if (m_Verbose) std::cout<<"The isocenter is at "<< m_IsoCenter <<"..." <<std::endl;
100 // Get the rotation parameter array
101 itk::Array<double> rotationParameters(3);
102 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
103 if (m_Verbose)std::cout<<"The projection angle is "<< m_ProjectionAngle <<"° (0° being lateral left)..."<< std::endl;
104 rotationParameters[0]= 0.;
105 rotationParameters[1]= 0.;
106 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
108 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
109 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
110 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
111 typename InputImageType::PointType transformedCenter = rigidRotation * m_IsoCenter + rigidTranslation;
113 // Calculate the centered rotation matrix
114 itk::Matrix<double,4,4> centeredRotationMatrix = GetCenteredRotationMatrix3D(rotationParameters,transformedCenter);
116 // Compose this rotation with the rigid transform matrix
117 itk::Matrix<double,4,4> finalTransform = m_RigidTransformMatrix * centeredRotationMatrix ;
120 itk::Matrix<double,3,3> finalRotation = GetRotationalPartMatrix3D(finalTransform);
121 m_Transform->SetMatrix(finalRotation);
122 if (m_Verbose)std::cout<<"The final rotation is "<<finalRotation<<"..."<<std::endl;
124 // Set the translation
125 itk::Vector<double,3> finalTranslation = GetTranslationPartMatrix3D(finalTransform);
126 m_Transform->SetTranslation(finalTranslation);
127 if (m_Verbose)std::cout<<"The final translation is "<<finalTranslation<<"..."<<std::endl;
130 //====================================================================
131 //Define a resample filter for the projection
132 if (m_NumberOfThreadsIsGiven) m_Resampler->SetNumberOfThreads(m_NumberOfThreads);
133 m_Resampler->SetDefaultPixelValue(m_EdgePaddingValue);
134 if (m_Verbose)std::cout<<"The edge padding value is "<<m_EdgePaddingValue<<"..."<<std::endl;
136 //JV Original raycast does not take into account origin, but does implicit shift to center of volume
137 //JV patched in RayCastInterpolateImageFunctionWithOrigin
138 m_Interpolator->SetTransform(m_Transform);
139 m_Interpolator->SetThreshold(0.);
141 // Focalpoint: we presume that for an angle of 0° the kV is lateral left (for the patient on his back), or on the positive X-axis.
142 typename InterpolatorType::InputPointType focalpoint;
143 focalpoint[0]= m_IsoCenter[0] + m_SourceToAxis;
144 focalpoint[1]= m_IsoCenter[1];
145 focalpoint[2]= m_IsoCenter[2];
146 m_Interpolator->SetFocalPoint(focalpoint);
147 if (m_Verbose)std::cout<<"The focalpoint is at "<< focalpoint <<"..."<< std::endl;
149 // Connect the interpolator and the transform with the m_Resampler
150 m_Resampler->SetInterpolator( m_Interpolator );
151 m_Resampler->SetTransform( m_Transform );
153 // Describe the output projection image
154 typename InputImageType::SizeType sizeOuput;
155 sizeOuput[0] = 1; // number of pixels along X of the 2D DRR image
156 sizeOuput[1] = m_OutputSize[0]; // number of pixels along Y of the 2D DRR image
157 sizeOuput[2] = m_OutputSize[1]; // number of pixels along Z of the 2D DRR image
158 m_Resampler->SetSize( sizeOuput );
159 if (m_Verbose)std::cout<<"The output size is "<< m_OutputSize <<"..."<< std::endl;
161 typename InputImageType::SpacingType spacingOutput;
162 spacingOutput[0] = 1; // pixel spacing along X of the 2D DRR image [mm]
163 spacingOutput[1] = m_OutputSpacing[0]; // pixel spacing along Y of the 2D DRR image [mm]
164 spacingOutput[2] = m_OutputSpacing[1]; // pixel spacing along Y of the 2D DRR image [mm]
165 m_Resampler->SetOutputSpacing( spacingOutput );
166 if (m_Verbose)std::cout<<"The output size is "<< m_OutputSpacing <<"..."<< std::endl;
168 // The position of the DRR is specified, we presume that for an angle of 0° the flatpanel is located at the negative x-axis
169 // JV -1 seems to correspond better with shearwarp of Simon Rit
170 typename InterpolatorType::InputPointType originOutput;
171 originOutput[0] = m_IsoCenter[0]- (m_SourceToScreen - m_SourceToAxis);
173 originOutput[1] = m_IsoCenter[1]-static_cast<double>(sizeOuput[1]-1)*spacingOutput[1]/2.0 - m_PanelShift;
174 originOutput[2] = m_IsoCenter[2]-static_cast<double>(sizeOuput[2]-1)*spacingOutput[2]/2.0;
175 m_Resampler->SetOutputOrigin( originOutput );
176 if (m_Verbose)std::cout<<"The origin of the flat panel is at "<< originOutput <<",..."<< std::endl;
178 // We define the region to be extracted. Projection was in the YZ plane, X should be set to zero
179 typename InternalImageType::RegionType::SizeType sizeTemp=sizeOuput;
181 typename InternalImageType::IndexType startTemp; //=temp->GetLargestPossibleRegion().GetIndex();
184 typename InternalImageType::RegionType desiredRegion;
185 desiredRegion.SetSize( sizeTemp );
186 desiredRegion.SetIndex( startTemp );
187 m_ExtractImageFilter->SetExtractionRegion( desiredRegion );
189 // Prepare the Flip Image filter
190 typename FlipImageFilterType::FlipAxesArrayType flipArray;
193 m_FlipImageFilter->SetFlipAxes( flipArray );
195 // Initialization complete
196 m_IsInitialized=true;
200 //=========================================================================================================================
202 template <class InputImageType, class OutputImageType>
204 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
208 //==================================================================
209 // If geometry changed reinitialize
210 if (! m_IsInitialized) this->Initialize();
213 m_Resampler->SetInput( m_Input );
215 //==================================================================
216 // Execute the filter
217 if (m_Verbose)std::cout<<"Starting the projection..."<<std::endl;
219 m_Resampler->Update();
221 catch( itk::ExceptionObject & err )
223 std::cerr << "ExceptionObject caught! Projection failed!" << std::endl;
224 std::cerr << err << std::endl;
227 //==================================================================
228 // Make a 2D image out of it
229 typename InternalImageType::Pointer temp=m_Resampler->GetOutput();
230 m_ExtractImageFilter->SetInput(temp);
232 // We should flip the Y axis
233 m_FlipImageFilter->SetInput(m_ExtractImageFilter->GetOutput());
234 m_FlipImageFilter->Update();
237 m_Output= m_FlipImageFilter->GetOutput();
239 m_Output->SetOrigin(m_OutputOrigin);
242 //=========================================================================================================================
244 template <class InputImageType, class OutputImageType>
245 typename OutputImageType::Pointer
246 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
249 //JV it is not safe to repeatedly call getoutput/ always call Update first