CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class vProjector 00004 00005 - separators: 00006 - doxygen: 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "vProjector.hh" 00019 #include "vScanner.hh" 00020 #include "vDataFile.hh" 00021 00022 // ===================================================================== 00023 // --------------------------------------------------------------------- 00024 // --------------------------------------------------------------------- 00025 // ===================================================================== 00026 00027 vProjector::vProjector() 00028 { 00029 // Affect default values 00030 mp_Scanner = NULL; 00031 mp_ImageDimensionsAndQuantification = NULL; 00032 mp_sizeVox[0] = -1.; 00033 mp_sizeVox[1] = -1.; 00034 mp_sizeVox[2] = -1.; 00035 mp_nbVox[0] = -1; 00036 mp_nbVox[1] = -1; 00037 mp_nbVox[2] = -1; 00038 m_nbVoxXY = -1; 00039 mp_halfFOV[0] = -1.; 00040 mp_halfFOV[1] = -1.; 00041 mp_halfFOV[2] = -1.; 00042 m_applyTOF = -1; 00043 m_applyPOI = false; 00044 m_compatibleWithSPECTAttenuationCorrection = false; 00045 m_verbose = 0; 00046 m_checked = false; 00047 m_initialized = false; 00048 } 00049 00050 // ===================================================================== 00051 // --------------------------------------------------------------------- 00052 // --------------------------------------------------------------------- 00053 // ===================================================================== 00054 00055 vProjector::~vProjector() 00056 { 00057 } 00058 00059 // ===================================================================== 00060 // --------------------------------------------------------------------- 00061 // --------------------------------------------------------------------- 00062 // ===================================================================== 00063 00064 int vProjector::SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification* ap_ImageDimensionsAndQuantification) 00065 { 00066 // Check that the parameter is not NULL 00067 if (ap_ImageDimensionsAndQuantification==NULL) 00068 { 00069 Cerr("***** vProjector::SetImageDimensionsAndQuantification() -> Input image dimensions object is null !" << endl); 00070 return 1; 00071 } 00072 // Affect image dimensions 00073 mp_sizeVox[0] = ap_ImageDimensionsAndQuantification->GetVoxSizeX(); 00074 mp_sizeVox[1] = ap_ImageDimensionsAndQuantification->GetVoxSizeY(); 00075 mp_sizeVox[2] = ap_ImageDimensionsAndQuantification->GetVoxSizeZ(); 00076 mp_nbVox[0] = ap_ImageDimensionsAndQuantification->GetNbVoxX(); 00077 mp_nbVox[1] = ap_ImageDimensionsAndQuantification->GetNbVoxY(); 00078 mp_nbVox[2] = ap_ImageDimensionsAndQuantification->GetNbVoxZ(); 00079 m_nbVoxXY = mp_nbVox[0] * mp_nbVox[1]; 00080 mp_halfFOV[0] = mp_sizeVox[0] * ((FLTNB)mp_nbVox[0]) / 2.; 00081 mp_halfFOV[1] = mp_sizeVox[1] * ((FLTNB)mp_nbVox[1]) / 2.; 00082 mp_halfFOV[2] = mp_sizeVox[2] * ((FLTNB)mp_nbVox[2]) / 2.; 00083 mp_ImageDimensionsAndQuantification = ap_ImageDimensionsAndQuantification; 00084 // Normal end 00085 return 0; 00086 } 00087 00088 // ===================================================================== 00089 // --------------------------------------------------------------------- 00090 // --------------------------------------------------------------------- 00091 // ===================================================================== 00092 00093 int vProjector::SetTOFAndPOIOptions(int a_dataType, bool a_applyTOF, bool a_applyPOI) 00094 { 00095 // Case 1: list-mode data type 00096 if (a_dataType==MODE_LIST) 00097 { 00098 // Affect TOF option 00099 if (a_applyTOF) m_applyTOF = USE_TOFPOS; 00100 else m_applyTOF = USE_NOTOF; 00101 // Affect POI option 00102 m_applyPOI = a_applyPOI; 00103 } 00104 // Case 2: histogram mode 00105 else 00106 { 00107 // Cannot use POI with histogram mode 00108 if (a_applyPOI) 00109 { 00110 Cerr("***** vProjector::SetTOFAndPOIOptions() -> POI (position of interaction) has no sense with histogram mode !" << endl); 00111 return 1; 00112 } 00113 // Affect TOF option 00114 if (a_applyTOF) m_applyTOF = USE_TOFBIN; 00115 else m_applyTOF = USE_NOTOF; 00116 } 00117 // Normal end 00118 return 0; 00119 } 00120 00121 // ===================================================================== 00122 // --------------------------------------------------------------------- 00123 // --------------------------------------------------------------------- 00124 // ===================================================================== 00125 00126 void vProjector::ShowCommonHelp() 00127 { 00128 // Return when using MPI and mpi_rank is not 0 00129 #ifdef CASTOR_MPI 00130 int mpi_rank = 0; 00131 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00132 if (mpi_rank!=0) return; 00133 #endif 00134 // Show help 00135 cout << "------------------------------------------------------------------" << endl; 00136 cout << "----- Common options for all projectors" << endl; 00137 cout << "------------------------------------------------------------------" << endl; 00138 cout << "Currently no options implemented." << endl; 00139 } 00140 00141 // ===================================================================== 00142 // --------------------------------------------------------------------- 00143 // --------------------------------------------------------------------- 00144 // ===================================================================== 00145 00146 void vProjector::ShowHelp() 00147 { 00148 // Call the specific help function from the children 00149 ShowHelpSpecific(); 00150 // Then, say if the child projector in use is compatible with SPECT attenuation correction or not 00151 if (m_compatibleWithSPECTAttenuationCorrection) cout << "This projector is compatible with SPECT attenuation correction." << endl; 00152 else cout << "This projector is NOT compatible with SPECT attenuation correction." << endl; 00153 } 00154 00155 // ===================================================================== 00156 // --------------------------------------------------------------------- 00157 // --------------------------------------------------------------------- 00158 // ===================================================================== 00159 00160 int vProjector::ReadCommonOptionsList(const string& a_optionsList) 00161 { 00162 // Currently no common options implemented 00163 return 0; 00164 } 00165 00166 // ===================================================================== 00167 // --------------------------------------------------------------------- 00168 // --------------------------------------------------------------------- 00169 // ===================================================================== 00170 00171 int vProjector::CheckParameters() 00172 { 00173 // Check scanner 00174 if (mp_Scanner==NULL) 00175 { 00176 Cerr("***** vProjector::CheckParameters() -> vScanner is null !" << endl); 00177 return 1; 00178 } 00179 // Check image dimensions 00180 if (mp_ImageDimensionsAndQuantification==NULL) 00181 { 00182 Cerr("***** vProjector::CheckParameters() -> oImageDimensionsAndQuantification is null !" << endl); 00183 return 1; 00184 } 00185 if (mp_sizeVox[0]<=0. || mp_sizeVox[1]<=0. || mp_sizeVox[2]<=0.) 00186 { 00187 Cerr("***** vProjector::CheckParameters() -> One or more voxel sizes is negative or null !" << endl); 00188 return 1; 00189 } 00190 if (mp_nbVox[0]<=0 || mp_nbVox[1]<=0 || mp_nbVox[2]<=0) 00191 { 00192 Cerr("***** vProjector::CheckParameters() -> One or more number of voxels is negative or null !" << endl); 00193 return 1; 00194 } 00195 // Check TOF 00196 if ( m_applyTOF!=USE_TOFPOS && m_applyTOF!=USE_TOFBIN && m_applyTOF!=USE_NOTOF ) 00197 { 00198 Cerr("***** vProjector::CheckParameters() -> TOF flag is incorrect or not set !" << endl); 00199 return 1; 00200 } 00201 // Check verbose level 00202 if (m_verbose<0) 00203 { 00204 Cerr("***** vProjector::CheckParameters() -> Verbose level is negative !" << endl); 00205 return 1; 00206 } 00207 // Check parameters of the child class 00208 if (CheckSpecificParameters()) 00209 { 00210 Cerr("***** vProjector::CheckParameters() -> An error occurred while checking parameters of the child projector class !" << endl); 00211 return 1; 00212 } 00213 // Normal end 00214 m_checked = true; 00215 return 0; 00216 } 00217 00218 // ===================================================================== 00219 // --------------------------------------------------------------------- 00220 // --------------------------------------------------------------------- 00221 // ===================================================================== 00222 00223 int vProjector::Initialize() 00224 { 00225 // First check that the parameters has been checked ! 00226 if (!m_checked) 00227 { 00228 Cerr("***** vProjector::Initialize() -> Must call the CheckParameters() function before initialization !" << endl); 00229 return 1; 00230 } 00231 // Call the intialize function specific to the children 00232 if (InitializeSpecific()) 00233 { 00234 Cerr("***** vProjector::Initialize() -> A problem occured while calling the specific initialization of the child projector !" << endl); 00235 return 1; 00236 } 00237 // Normal end 00238 m_initialized = true; 00239 return 0; 00240 } 00241 00242 // ===================================================================== 00243 // --------------------------------------------------------------------- 00244 // --------------------------------------------------------------------- 00245 // ===================================================================== 00246 00247 INTNB vProjector::EstimateMaxNumberOfVoxelsPerLine() 00248 { 00249 return mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); 00250 } 00251 00252 // ===================================================================== 00253 // --------------------------------------------------------------------- 00254 // --------------------------------------------------------------------- 00255 // ===================================================================== 00256 00257 int vProjector::Project(int a_direction, oProjectionLine* ap_ProjectionLine, uint32_t* ap_index1, uint32_t* ap_index2, int a_nbIndices) 00258 { 00259 #ifdef CASTOR_DEBUG 00260 if (!m_initialized) 00261 { 00262 Cerr("***** vProjector::Project() -> Called while not initialized !" << endl); 00263 Exit(EXIT_DEBUG); 00264 } 00265 #endif 00266 00267 #ifdef CASTOR_VERBOSE 00268 if (m_verbose>=10) 00269 { 00270 string direction = ""; 00271 if (a_direction==FORWARD) direction = "forward"; 00272 else direction = "backward"; 00273 Cout("vProjector::Project() -> Project line '" << ap_ProjectionLine << "' for " << direction << " projector" << endl); 00274 } 00275 #endif 00276 00277 // --------------------------------------------------------------------------------------- 00278 // First: Get cartesian coordinates from the scanner and average positions if compression. 00279 // Here, we differentiate the case with no compression (i.e. a_nbIndices==1) from the 00280 // compression case, because it will avoid to perform useless computation without 00281 // compression. However, it produces some duplication of code parts as a compromise. 00282 // --------------------------------------------------------------------------------------- 00283 00284 // ______________________________________________________________________________________ 00285 // Case1: no compression (i.e. a_nbIndices==1) 00286 if (a_nbIndices==1) 00287 { 00288 // Set indices of the line 00289 ap_ProjectionLine->SetIndex1(((int)(ap_index1[0]))); 00290 ap_ProjectionLine->SetIndex2(((int)(ap_index2[0]))); 00291 // Get positions and orientations from the scanner (mean depth of interaction intrinsicaly taken into account), taking POI into account if any 00292 if (mp_Scanner->GetPositionsAndOrientations( ((int)(ap_index1[0])), ((int)(ap_index2[0])), 00293 ap_ProjectionLine->GetPosition1(), ap_ProjectionLine->GetPosition2(), 00294 ap_ProjectionLine->GetOrientation1(), ap_ProjectionLine->GetOrientation2(), 00295 ap_ProjectionLine->GetPOI1(), ap_ProjectionLine->GetPOI2() )) 00296 { 00297 Cerr("***** vProjector::Project() -> A problem occured while getting positions and orientations from scanner !" << endl); 00298 return 1; 00299 } 00300 } 00301 // ______________________________________________________________________________________ 00302 // Case2: compression (i.e. a_nbIndices>1) 00303 else 00304 { 00305 // Set default indices of the line to -1 when compression 00306 ap_ProjectionLine->SetIndex1(-1); 00307 ap_ProjectionLine->SetIndex2(-1); 00308 // Buffer pointers for positions and orientations handled by the projection line 00309 FLTNB* position1 = ap_ProjectionLine->GetPosition1(); 00310 FLTNB* position2 = ap_ProjectionLine->GetPosition2(); 00311 FLTNB* orientation1 = ap_ProjectionLine->GetOrientation1(); 00312 FLTNB* orientation2 = ap_ProjectionLine->GetOrientation2(); 00313 FLTNB* buffer_position1 = ap_ProjectionLine->GetBufferPosition1(); 00314 FLTNB* buffer_position2 = ap_ProjectionLine->GetBufferPosition2(); 00315 FLTNB* buffer_orientation1 = ap_ProjectionLine->GetBufferOrientation1(); 00316 FLTNB* buffer_orientation2 = ap_ProjectionLine->GetBufferOrientation2(); 00317 // Zero the position and orientation 00318 for (int i=0; i<3; i++) 00319 { 00320 position1[i] = 0.; 00321 position2[i] = 0.; 00322 orientation1[i] = 0.; 00323 orientation2[i] = 0.; 00324 } 00325 // Loop on provided indices 00326 for (int l=0; l<a_nbIndices; l++) 00327 { 00328 // Get positions and orientations from scanner (mean depth of interaction intrinsicaly taken into account), taking POI into account if any 00329 if (mp_Scanner->GetPositionsAndOrientations( ((int)(ap_index1[l])), ((int)(ap_index2[l])), 00330 buffer_position1, buffer_position2, 00331 buffer_orientation1, buffer_orientation2, 00332 ap_ProjectionLine->GetPOI1(), ap_ProjectionLine->GetPOI2() )) 00333 { 00334 Cerr("***** vProjector::Project() -> A problem occured while getting positions and orientations from scanner !" << endl); 00335 return 1; 00336 } 00337 // Add those contributions to the mean position and orientation 00338 for (int i=0; i<3; i++) 00339 { 00340 position1[i] += buffer_position1[i]; 00341 position2[i] += buffer_position2[i]; 00342 orientation1[i] += buffer_orientation1[i]; 00343 orientation2[i] += buffer_orientation2[i]; 00344 } 00345 } 00346 // Average positions and orientations 00347 for (int i=0; i<3; i++) 00348 { 00349 position1[i] /= ((FLTNB)a_nbIndices); 00350 position2[i] /= ((FLTNB)a_nbIndices); 00351 orientation1[i] /= ((FLTNB)a_nbIndices); 00352 orientation2[i] /= ((FLTNB)a_nbIndices); 00353 } 00354 } 00355 00356 // -------------------------------------------------------------- 00357 // Second: Modify the end points coordinates from common options, 00358 // random, offset, and LOR displacements. 00359 // ----------------------------------------------------------- 00360 00361 // Apply common options TODO 00362 00363 // Apply LORs displacement TODO 00364 00365 // Apply offset 00366 ap_ProjectionLine->ApplyOffset(); 00367 00368 // ----------------------------------------------------------- 00369 // Third: project the line 00370 // ----------------------------------------------------------- 00371 00372 // Compute LOR length 00373 ap_ProjectionLine->ComputeLineLength(); 00374 00375 // Switch on different TOF options 00376 switch (m_applyTOF) 00377 { 00378 case USE_NOTOF: 00379 if (ProjectWithoutTOF( a_direction, ap_ProjectionLine )) 00380 { 00381 Cerr("***** vProjector::Project() -> A problem occured while projecting a line without time-of-flight !" << endl); 00382 return 1; 00383 } 00384 break; 00385 case USE_TOFPOS: 00386 if (ProjectWithTOFPos( a_direction, ap_ProjectionLine )) 00387 { 00388 Cerr("***** vProjector::Project() -> A problem occured while projecting a line with time-of-flight position !" << endl); 00389 return 1; 00390 } 00391 break; 00392 case USE_TOFBIN: 00393 if (ProjectWithTOFBin( a_direction, ap_ProjectionLine )) 00394 { 00395 Cerr("***** vProjector::Project() -> A problem occured while projecting a line with binned time-of-flight !" << endl); 00396 return 1; 00397 } 00398 break; 00399 // No default 00400 } 00401 00402 #ifdef CASTOR_VERBOSE 00403 if (m_verbose>=10) 00404 { 00405 Cout("vProjector::Project() -> Exit function" << endl); 00406 } 00407 #endif 00408 00409 // End 00410 return 0; 00411 } 00412 00413 // ===================================================================== 00414 // --------------------------------------------------------------------- 00415 // --------------------------------------------------------------------- 00416 // =====================================================================