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.;
42 m_RigidTransformMatrix.SetIdentity();
43 m_EdgePaddingValue=itk::NumericTraits<OutputPixelType>::Zero;//density images
46 m_Transform=TransformType::New();
47 m_Resampler=ResampleImageFilterType::New();
48 m_Interpolator= InterpolatorType::New();
49 m_ExtractImageFilter=ExtractImageFilterType::New();
50 m_FlipImageFilter=FlipImageFilterType::New();
52 // Parameters for output
53 this->m_OutputSpacing.Fill(0.8);
54 this->m_OutputOrigin.Fill(0.0);
55 this->m_OutputSize.Fill( 512 );
60 //-----------------------------------------------------------------------
61 // Set output info from image
62 //-----------------------------------------------------------------------
63 template <class InputImageType, class OutputImageType>
65 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
66 ::SetOutputParametersFromImage ( OutputImagePointer image )
68 this->SetOutputOrigin( image->GetOrigin() );
69 this->SetOutputSpacing( image->GetSpacing() );
70 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
74 //-----------------------------------------------------------------------
75 // Set output info from image
76 //-----------------------------------------------------------------------
77 template <class InputImageType, class OutputImageType>
79 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
80 ::SetOutputParametersFromImage ( OutputImageConstPointer image )
82 this->SetOutputOrigin( image->GetOrigin() );
83 this->SetOutputSpacing( image->GetSpacing() );
84 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
88 //=========================================================================================================================
90 template <class InputImageType, class OutputImageType>
92 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
93 ::Initialize(void) throw (itk::ExceptionObject)
96 //====================================================================
97 // Define the transform: composition of rigid transfo around origin and a centered rotation
98 // The center of rotation should be placed at the isocenter!
100 if (m_Verbose) std::cout<<"The isocenter is at "<< m_IsoCenter <<"..." <<std::endl;
102 // Get the rotation parameter array
103 itk::Array<double> rotationParameters(3);
104 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
105 if (m_Verbose)std::cout<<"The projection angle is "<< m_ProjectionAngle <<"° (0° being lateral left)..."<< std::endl;
106 rotationParameters[0]= 0.;
107 rotationParameters[1]= 0.;
108 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
110 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
111 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
112 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
113 typename InputImageType::PointType transformedCenter = rigidRotation * m_IsoCenter + rigidTranslation;
115 // Calculate the centered rotation matrix
116 itk::Matrix<double,4,4> centeredRotationMatrix = GetCenteredRotationMatrix3D(rotationParameters,transformedCenter);
118 // Compose this rotation with the rigid transform matrix
119 itk::Matrix<double,4,4> finalTransform = m_RigidTransformMatrix * centeredRotationMatrix ;
122 itk::Matrix<double,3,3> finalRotation = GetRotationalPartMatrix3D(finalTransform);
123 m_Transform->SetMatrix(finalRotation);
124 if (m_Verbose)std::cout<<"The final rotation is "<<finalRotation<<"..."<<std::endl;
126 // Set the translation
127 itk::Vector<double,3> finalTranslation = GetTranslationPartMatrix3D(finalTransform);
128 m_Transform->SetTranslation(finalTranslation);
129 if (m_Verbose)std::cout<<"The final translation is "<<finalTranslation<<"..."<<std::endl;
132 //====================================================================
133 //Define a resample filter for the projection
134 if (m_NumberOfThreadsIsGiven) m_Resampler->SetNumberOfThreads(m_NumberOfThreads);
135 m_Resampler->SetDefaultPixelValue(m_EdgePaddingValue);
136 if (m_Verbose)std::cout<<"The edge padding value is "<<m_EdgePaddingValue<<"..."<<std::endl;
138 //JV Original raycast does not take into account origin, but does implicit shift to center of volume
139 //JV patched in RayCastInterpolateImageFunctionWithOrigin
140 m_Interpolator->SetTransform(m_Transform);
141 m_Interpolator->SetThreshold(0.);
143 // 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.
144 typename InterpolatorType::InputPointType focalpoint;
145 focalpoint[0]= m_IsoCenter[0] + m_SourceToAxis;
146 focalpoint[1]= m_IsoCenter[1];
147 focalpoint[2]= m_IsoCenter[2];
148 m_Interpolator->SetFocalPoint(focalpoint);
149 if (m_Verbose)std::cout<<"The focalpoint is at "<< focalpoint <<"..."<< std::endl;
151 // Connect the interpolator and the transform with the m_Resampler
152 m_Resampler->SetInterpolator( m_Interpolator );
153 m_Resampler->SetTransform( m_Transform );
155 // Describe the output projection image
156 typename InputImageType::SizeType sizeOuput;
157 sizeOuput[0] = 1; // number of pixels along X of the 2D DRR image
158 sizeOuput[1] = m_OutputSize[0]; // number of pixels along Y of the 2D DRR image
159 sizeOuput[2] = m_OutputSize[1]; // number of pixels along Z of the 2D DRR image
160 m_Resampler->SetSize( sizeOuput );
161 if (m_Verbose)std::cout<<"The output size is "<< m_OutputSize <<"..."<< std::endl;
163 typename InputImageType::SpacingType spacingOutput;
164 spacingOutput[0] = 1; // pixel spacing along X of the 2D DRR image [mm]
165 spacingOutput[1] = m_OutputSpacing[0]; // pixel spacing along Y of the 2D DRR image [mm]
166 spacingOutput[2] = m_OutputSpacing[1]; // pixel spacing along Y of the 2D DRR image [mm]
167 m_Resampler->SetOutputSpacing( spacingOutput );
168 if (m_Verbose)std::cout<<"The output spacing is "<< m_OutputSpacing <<"..."<< std::endl;
170 // The position of the DRR is specified, we presume that for an angle of 0° the flatpanel is located at the negative x-axis
171 // JV -1 seems to correspond better with shearwarp of Simon Rit
172 typename InterpolatorType::InputPointType originOutput;
173 originOutput[0] = m_IsoCenter[0]- (m_SourceToScreen - m_SourceToAxis);
175 originOutput[1] = m_IsoCenter[1]-static_cast<double>(sizeOuput[1]-1)*spacingOutput[1]/2.0 - m_PanelShift[0];
176 originOutput[2] = m_IsoCenter[2]-static_cast<double>(sizeOuput[2]-1)*spacingOutput[2]/2.0 - m_PanelShift[1];
177 m_Resampler->SetOutputOrigin( originOutput );
178 if (m_Verbose)std::cout<<"The origin of the flat panel is at "<< originOutput <<",..."<< std::endl;
180 // We define the region to be extracted. Projection was in the YZ plane, X should be set to zero
181 typename InternalImageType::RegionType::SizeType sizeTemp=sizeOuput;
183 typename InternalImageType::IndexType startTemp; //=temp->GetLargestPossibleRegion().GetIndex();
186 typename InternalImageType::RegionType desiredRegion;
187 desiredRegion.SetSize( sizeTemp );
188 desiredRegion.SetIndex( startTemp );
189 m_ExtractImageFilter->SetExtractionRegion( desiredRegion );
191 // Prepare the Flip Image filter
192 typename FlipImageFilterType::FlipAxesArrayType flipArray;
195 m_FlipImageFilter->SetFlipAxes( flipArray );
197 // Initialization complete
198 m_IsInitialized=true;
202 //=========================================================================================================================
204 template <class InputImageType, class OutputImageType>
206 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
210 //==================================================================
211 // If geometry changed reinitialize
212 if (! m_IsInitialized) this->Initialize();
215 m_Resampler->SetInput( m_Input );
217 //==================================================================
218 // Execute the filter
219 if (m_Verbose)std::cout<<"Starting the projection..."<<std::endl;
221 m_Resampler->Update();
223 catch( itk::ExceptionObject & err )
225 std::cerr << "ExceptionObject caught! Projection failed!" << std::endl;
226 std::cerr << err << std::endl;
229 //==================================================================
230 // Make a 2D image out of it
231 typename InternalImageType::Pointer temp=m_Resampler->GetOutput();
232 m_ExtractImageFilter->SetInput(temp);
234 // We should flip the Y axis
235 m_FlipImageFilter->SetInput(m_ExtractImageFilter->GetOutput());
236 m_FlipImageFilter->Update();
239 m_Output= m_FlipImageFilter->GetOutput();
241 m_Output->SetOrigin(m_OutputOrigin);
244 //=========================================================================================================================
246 template <class InputImageType, class OutputImageType>
247 typename OutputImageType::Pointer
248 ConeBeamProjectImageFilter<InputImageType, OutputImageType>
251 //JV it is not safe to repeatedly call getoutput/ always call Update first