![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iProjectorIRIS 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "iProjectorIRIS.hh" 00019 #include "sOutputManager.hh" 00020 00021 // ===================================================================== 00022 // --------------------------------------------------------------------- 00023 // --------------------------------------------------------------------- 00024 // ===================================================================== 00025 00026 iProjectorIRIS::iProjectorIRIS() : vProjector() 00027 { 00028 mp_pathToIDRFFiles = NULL; 00029 mp_IDRF = NULL; 00030 m2p_IDRF_CDFs = NULL; 00031 m_nbLinesPerLOR = -1; 00032 m_nVoxDepthIDRF = -1; 00033 m_nVoxTransaxialIDRF = -1; 00034 m_nVoxAxialIDRF = -1; 00035 m_nVoxXYZIDRF = -1; 00036 m_nBetaAnglesIDRF = -1; 00037 m_nAlphaAnglesIDRF = -1; 00038 m_sizeVoxDepthIDRF = -1.; 00039 m_sizeVoxTransaxialIDRF = -1.; 00040 m_sizeVoxAxialIDRF = -1.; 00041 m_stepBetaAnglesIDRF = -1.; 00042 m_stepAlphaAnglesIDRF = -1.; 00043 // This projector is not compatible with SPECT attenuation correction because 00044 // the voxels contributing to the line are not strictly ordered with respect to 00045 // their distance to point 2 (due to the use of multiple lines that are 00046 // stack one after the other) 00047 m_compatibleWithSPECTAttenuationCorrection = false; 00048 } 00049 00050 // ===================================================================== 00051 // --------------------------------------------------------------------- 00052 // --------------------------------------------------------------------- 00053 // ===================================================================== 00054 00055 iProjectorIRIS::~iProjectorIRIS() 00056 { 00057 if (mp_pathToIDRFFiles) delete[] mp_pathToIDRFFiles; 00058 00059 if (m_nAlphaAnglesIDRF>0 && m_nBetaAnglesIDRF>0 && m2p_IDRF_CDFs) 00060 for (int a=0 ; a<m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF ; a++) 00061 if (m2p_IDRF_CDFs[a]) delete[] m2p_IDRF_CDFs[a]; 00062 00063 if (m2p_IDRF_CDFs) delete[] m2p_IDRF_CDFs; 00064 00065 if (mp_IDRF) delete[] mp_IDRF; 00066 } 00067 00068 // ===================================================================== 00069 // --------------------------------------------------------------------- 00070 // --------------------------------------------------------------------- 00071 // ===================================================================== 00072 00073 int iProjectorIRIS::ReadConfigurationFile(const string& a_configurationFile) 00074 { 00075 if (ReadDataASCIIFile(a_configurationFile, "number lines per LOR", &m_nbLinesPerLOR, 1, 1) ) 00076 { 00077 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the number of lines per LOR from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00078 return 1; 00079 } 00080 00081 if (ReadDataASCIIFile(a_configurationFile, "number voxels depth", &m_nVoxDepthIDRF, 1, 1) || 00082 ReadDataASCIIFile(a_configurationFile, "number voxels transaxial", &m_nVoxTransaxialIDRF, 1, 1) || 00083 ReadDataASCIIFile(a_configurationFile, "number voxels axial", &m_nVoxAxialIDRF, 1, 1) ) 00084 { 00085 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the IDRF volume dimensions (nb voxels) from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00086 return 1; 00087 } 00088 00089 // Total number of voxels in the IDRF 00090 m_nVoxXYZIDRF = m_nVoxDepthIDRF*m_nVoxTransaxialIDRF*m_nVoxAxialIDRF; 00091 00092 if (ReadDataASCIIFile(a_configurationFile, "size voxels depth", &m_sizeVoxDepthIDRF, 1, 1) || 00093 ReadDataASCIIFile(a_configurationFile, "size voxels transaxial", &m_sizeVoxTransaxialIDRF, 1, 1) || 00094 ReadDataASCIIFile(a_configurationFile, "size voxels axial", &m_sizeVoxAxialIDRF, 1, 1) ) 00095 { 00096 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the voxel sizes of the IDRF volumes from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00097 return 1; 00098 } 00099 00100 if (ReadDataASCIIFile(a_configurationFile, "number beta angles", &m_nBetaAnglesIDRF, 1, 1) || 00101 ReadDataASCIIFile(a_configurationFile, "number alpha angles", &m_nAlphaAnglesIDRF, 1, 1) ) 00102 { 00103 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the number of alpha/beta angles from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00104 return 1; 00105 } 00106 00107 if (ReadDataASCIIFile(a_configurationFile, "step beta angles", &m_stepBetaAnglesIDRF, 1, 1) || 00108 ReadDataASCIIFile(a_configurationFile, "step alpha angles", &m_stepAlphaAnglesIDRF, 1, 1) ) 00109 { 00110 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the alpha/beta angles steps from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00111 return 1; 00112 } 00113 00114 // Convert angle steps in radians 00115 m_stepAlphaAnglesIDRF *= M_PI/180.; 00116 m_stepBetaAnglesIDRF *= M_PI/180.; 00117 00118 mp_pathToIDRFFiles = new string[m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF]; 00119 00120 mp_IDRF = new float[m_nVoxXYZIDRF]; 00121 m2p_IDRF_CDFs = new float*[m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF]; 00122 00123 for(int a=0 ; a<m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF ; a++) 00124 m2p_IDRF_CDFs[a] = new float[m_nVoxXYZIDRF]; 00125 00126 if (ReadDataASCIIFile(a_configurationFile, "IDRFs path", mp_pathToIDRFFiles, 1, m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF, 1) ) 00127 { 00128 Cerr("***** iProjectorIRIS::ReadConfigurationFile() -> An error occurred while trying to read the path to the IDRF files, from the IRIS configuration file at :" << a_configurationFile << " !" << endl); 00129 Cerr("***** ("<< m_nAlphaAnglesIDRF*m_nBetaAnglesIDRF << " lines expected)" << endl); 00130 return 1; 00131 } 00132 00133 // Normal end 00134 return 0; 00135 } 00136 00137 // ===================================================================== 00138 // --------------------------------------------------------------------- 00139 // --------------------------------------------------------------------- 00140 // ===================================================================== 00141 00142 int iProjectorIRIS::ReadOptionsList(const string& a_optionsList) 00143 { 00144 // Configuration is possible only through a configuration file 00145 Cerr("***** iProjectorIRIS::ReadOptionsList() -> This projector should be initialized with a configuration file, and not parameters (-proj IRIS:path_to_configuration_file) !" << endl); 00146 return 1; 00147 } 00148 00149 // ===================================================================== 00150 // --------------------------------------------------------------------- 00151 // --------------------------------------------------------------------- 00152 // ===================================================================== 00153 00154 void iProjectorIRIS::ShowHelpSpecific() 00155 { 00156 cout << "Multiray projector requiring Intrinsic Detector Response Functions (IDRFs) models" << endl; 00157 cout << "The projector will generate an user-defined number of lines using from pairs of random points generated using the IDRFs models" << endl; 00158 cout << "The lines will be rendered using the Incremental Siddon projector (incrementalSiddon)" << endl; 00159 cout << "This projector requires an initialization file and Monte-Carlo generated IDRFs. For more informations, please read [article]" << endl; 00160 cout << "INITIALIZATION FILE FIELDS :" << endl; 00161 cout << "number lines per LOR : number of lines generated to estimate the CDRF (int32)" << endl; 00162 cout << "number voxels depth : number of voxels in the depth direction in the IDRF volumes (int32)" << endl; 00163 cout << "number voxels transaxial : number of voxels in the transaxial direction in the IDRF volumes (int32)" << endl; 00164 cout << "number voxels axial : number of voxels in the axial direction in the IDRF volumes (int32)" << endl; 00165 cout << "size voxels depth : size of voxels in the depth direction of the IDRF volumes (float32)" << endl; 00166 cout << "size voxels transaxial : size of voxels in the transaxial direction of the IDRF volume (float32)" << endl; 00167 cout << "size voxels axial : size of voxels in the axial direction of the IDRF volume (float32)" << endl; 00168 cout << "number beta angles : number of axial angles in the IDRF volumes (int32)" << endl; 00169 cout << "number alpha angles : number of transaxial angles in the IDRF volumes (int32)" << endl; 00170 cout << "step beta angles : axial angular steps in the IDRF volumes (float32)" << endl; 00171 cout << "step alpha angles : transaxial angular steps in the IDRF volumes (float32)" << endl; 00172 cout << "IDRFs path : path to the IDRF .vol files (string)" << endl; 00173 cout << " : each path should be entered on a separated line" << endl; 00174 cout << " " << endl; 00175 cout << "Each .vol files should contain a header of 6 float32 fields :" << endl; 00176 cout << "number voxels axial ; number voxels transaxial ; number voxels depth ; size voxels axial ; size voxels transaxial ; size voxels depth" << endl; 00177 cout << "... followed by number voxels depth*number voxels transaxial*number voxels axial float32 elements corresponding to the IDRF coefficients" << endl; 00178 } 00179 00180 // ===================================================================== 00181 // --------------------------------------------------------------------- 00182 // --------------------------------------------------------------------- 00183 // ===================================================================== 00184 00185 int iProjectorIRIS::CheckSpecificParameters() 00186 { 00187 if(mp_pathToIDRFFiles==NULL) 00188 { 00189 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> path to IDRF files not initialized !" << endl); 00190 return 1; 00191 } 00192 if(mp_IDRF==NULL || m2p_IDRF_CDFs==NULL) 00193 { 00194 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> IDRF coefficients not initialized !" << endl); 00195 return 1; 00196 } 00197 if(m_nbLinesPerLOR<0) 00198 { 00199 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> number of lines per LOR not initialized !" << endl); 00200 return 1; 00201 } 00202 if(m_nVoxDepthIDRF<0 || m_nVoxTransaxialIDRF<0 || m_nVoxAxialIDRF<0 || m_nVoxXYZIDRF<0) 00203 { 00204 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> number of voxels for the IDRF volumes not initialized !" << endl); 00205 return 1; 00206 } 00207 if(m_sizeVoxDepthIDRF<0 || m_sizeVoxTransaxialIDRF<0 || m_sizeVoxAxialIDRF<0) 00208 { 00209 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> voxel sizes of the IDRF volumes not initialized !" << endl); 00210 return 1; 00211 } 00212 if(m_nBetaAnglesIDRF<0 || m_nAlphaAnglesIDRF<0) 00213 { 00214 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> number of angles for the IDRF volumes not initialized !" << endl); 00215 return 1; 00216 } 00217 if(m_stepBetaAnglesIDRF<0 || m_stepAlphaAnglesIDRF<0) 00218 { 00219 Cerr("***** iProjectorIRIS::CheckSpecificParameters() -> step of the IDRF angles not initialized !" << endl); 00220 return 1; 00221 } 00222 // Normal end 00223 return 0; 00224 } 00225 00226 // ===================================================================== 00227 // --------------------------------------------------------------------- 00228 // --------------------------------------------------------------------- 00229 // ===================================================================== 00230 00231 int iProjectorIRIS::InitializeSpecific() 00232 { 00233 // Verbose 00234 if (m_verbose>=2) Cout("iProjectorIRIS::InitializeSpecific() -> Use IRIS projector" << endl); 00235 00236 // Read IDRFs coefficients from the raw data files, and write them on a single vector 00237 for(int a=0 ; a<m_nAlphaAnglesIDRF ; a++) // loop on alpha angles 00238 for(int b=0 ; b<m_nBetaAnglesIDRF ; b++) // loop on beta angles 00239 { 00240 // Open & check file existence 00241 ifstream IDRF_file(mp_pathToIDRFFiles[a*m_nBetaAnglesIDRF+b].c_str(), ios::binary| ios::in); 00242 if (!IDRF_file.is_open()) 00243 { 00244 Cerr("***** iProjectorIRIS::InitializeSpecific() -> Error reading the IRDF .vol file at the path: '" << mp_pathToIDRFFiles[a*m_nBetaAnglesIDRF+b] << "'" << endl); 00245 return 1; 00246 } 00247 // Check file size consistency 00248 IDRF_file.seekg(0, ios::end); 00249 size_t size_in_bytes = IDRF_file.tellg(); 00250 00251 if(((size_t)((6+m_nVoxXYZIDRF)*sizeof(float))) != size_in_bytes) 00252 { 00253 Cerr("***** iProjectorIRIS::InitializeSpecific() -> Error : Size of the .vol " << mp_pathToIDRFFiles[a*m_nBetaAnglesIDRF+b] << " is not consistent with the information provided by the user configuration file!" << endl); 00254 Cerr("***** iProjectorIRIS::InitializeSpecific() -> Expected size : "<< (6+m_nVoxXYZIDRF)*sizeof(float) << endl); 00255 Cerr("***** iProjectorIRIS::InitializeSpecific() -> Actual size : "<< size_in_bytes << endl << endl); 00256 return 1; 00257 } 00258 // Check header of the volume 00259 float vol_file_header[6]; 00260 // Go back to first character 00261 IDRF_file.seekg(0); 00262 for(int i=0 ; i<6 ; i++) 00263 IDRF_file.read((char*)&vol_file_header[i], sizeof(float)); 00264 // TODO : Take out header from .vol files. Just use the other checks and delete this condition 00265 if ((int)vol_file_header[0]!=m_nVoxAxialIDRF || 00266 (int)vol_file_header[1]!=m_nVoxTransaxialIDRF || 00267 (int)vol_file_header[2]!=m_nVoxDepthIDRF || 00268 fabs(vol_file_header[3]-m_sizeVoxAxialIDRF) > 0.0001 || 00269 fabs(vol_file_header[4]-m_sizeVoxTransaxialIDRF) > 0.0001 || 00270 fabs(vol_file_header[5]-m_sizeVoxDepthIDRF) > 0.0001 ) 00271 { 00272 Cerr("***** iProjectorIRIS::InitializeSpecific() -> Inconsistency between initialized and read header file in the .vol: '" << mp_pathToIDRFFiles[a*m_nBetaAnglesIDRF+b] << "'" << endl); 00273 return 1; 00274 } 00275 // Read raw data 00276 for(int i=0 ; i<m_nVoxXYZIDRF ; i++) 00277 { 00278 IDRF_file.read((char*)&vol_file_header[0], sizeof(float)); 00279 mp_IDRF[i] = vol_file_header[0]; 00280 } 00281 if(ComputeIDRF_CDF(a*m_nBetaAnglesIDRF + b) ) 00282 { 00283 Cerr("***** iProjectorIRIS::InitializeSpecific() -> An error occurs while computing the coefficient of the following IDRF volume: alpha " << a*m_stepAlphaAnglesIDRF << " beta " << b*m_stepBetaAnglesIDRF << endl); 00284 return 1; 00285 } 00286 IDRF_file.close(); 00287 } 00288 00289 // Normal end 00290 return 0; 00291 } 00292 00293 // ===================================================================== 00294 // --------------------------------------------------------------------- 00295 // --------------------------------------------------------------------- 00296 // ===================================================================== 00297 00298 INTNB iProjectorIRIS::EstimateMaxNumberOfVoxelsPerLine() 00299 { 00300 // Find the maximum number of voxels along a given dimension 00301 INTNB max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxX(); 00302 if (mp_ImageDimensionsAndQuantification->GetNbVoxY()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxY(); 00303 if (mp_ImageDimensionsAndQuantification->GetNbVoxZ()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxZ(); 00304 // We should have at most 4 voxels in a given plane, so multiply by 4 00305 // (note: this is not true however it ensures no overflow and is already quite optimized for RAM usage !) 00306 max_nb_voxels_in_dimension *= 4; 00307 // Finally multiply by the number of lines 00308 max_nb_voxels_in_dimension *= ((INTNB)m_nbLinesPerLOR); 00309 // Return the value 00310 return max_nb_voxels_in_dimension; 00311 } 00312 00313 // ===================================================================== 00314 // --------------------------------------------------------------------- 00315 // --------------------------------------------------------------------- 00316 // ===================================================================== 00317 00318 int iProjectorIRIS::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine ) 00319 { 00320 #ifdef CASTOR_DEBUG 00321 if (!m_initialized) 00322 { 00323 Cerr("***** iProjectorIRIS::ProjectWithoutTOF() -> Called while not initialized !" << endl); 00324 Exit(EXIT_DEBUG); 00325 } 00326 #endif 00327 00328 #ifdef CASTOR_VERBOSE 00329 if (m_verbose>=10) 00330 { 00331 string direction = ""; 00332 if (a_direction==FORWARD) direction = "forward"; 00333 else direction = "backward"; 00334 Cout("iProjectorIRIS::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl); 00335 } 00336 #endif 00337 00338 // Get central positions of the crystals 00339 FLTNB* event1 = ap_ProjectionLine->GetPosition1(); 00340 FLTNB* event2 = ap_ProjectionLine->GetPosition2(); 00341 00342 // Local variables to recover crystal central positions and orientations 00343 float x1, y1, z1, crystal_alpha1, x2, y2, z2, crystal_alpha2; 00344 00345 x1 = event1[0]; 00346 y1 = event1[1]; 00347 z1 = event1[2]; 00348 x2 = event2[0]; 00349 y2 = event2[1]; 00350 z2 = event2[2]; 00351 00352 // Get transaxial orientations of the two crystals // CHECK HERE 00353 FLTNB* orientation1 = ap_ProjectionLine->GetOrientation1(); 00354 FLTNB* orientation2 = ap_ProjectionLine->GetOrientation2(); 00355 //crystal_alpha1 = acosf(orientation1[0]); //LUT_alpha[id1]; 00356 //crystal_alpha2 = acosf(orientation2[0]); //LUT_alpha[id2]; 00357 crystal_alpha1 = atan2f(orientation1[0], orientation1[1]); //LUT_alpha[id1]; //HERE 00358 crystal_alpha2 = atan2f(orientation2[0], orientation2[1]); //LUT_alpha[id2]; //HERE 00359 00360 // Compute transaxial incident angles (alpha) between the LOR and the crystals 00361 float alpha1, alpha2, alpha_lor; 00362 00363 // Transaxial incident angle of the LOR 00364 alpha_lor = atan2f(x1-x2, y1-y2); // alpha_lor = atan2f( y1-y2 , x1-x2 ); //CHECK HERE 00365 00366 //if (ap_ProjectionLine->GetIndex1() == 0 && ap_ProjectionLine->GetIndex2() <= 615) 00367 // cout << crystal_alpha1 << " " << crystal_alpha2 << " rsector1: " << ap_ProjectionLine->GetIndex1()/22 << " rsector2: " << ap_ProjectionLine->GetIndex2()/22 << " angle_LOR: " << alpha_lor << endl; 00368 00369 00370 alpha1 = remainderf(alpha_lor-M_PI-crystal_alpha1, M_PI); 00371 if(alpha1>(M_PI/2.)) // Check if this operation and other surrounding should require double instead of float 00372 alpha1 -= M_PI; 00373 00374 alpha2 = remainderf(alpha_lor-M_PI-crystal_alpha2, M_PI); 00375 if(alpha2>(M_PI/2.)) 00376 alpha2 -= M_PI; 00377 00378 //alpha1 = remainderf(alpha_lor-Pi-crystal_alpha1, M_PI_2); 00379 //alpha2 = remainderf(alpha_lor-crystal_alpha2, M_PI_2); 00380 00381 // Compute axial incident angles (beta) between the LOR and the crystals 00382 float beta1, beta2, H; 00383 00384 H = sqrtf((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); 00385 00386 beta1 = atan2f( (z1-z2), H ); //CHECK HERE 00387 00388 beta1 = remainderf(beta1, M_PI); 00389 if(beta1>(M_PI/2.)) 00390 beta1 -= M_PI; 00391 00392 beta2 = -beta1; 00393 beta2 = remainderf(beta2, M_PI); 00394 if(beta2>(M_PI/2.)) 00395 beta2 -= M_PI; 00396 00397 00399 for (int i_line=0; i_line<m_nbLinesPerLOR; ++i_line) 00400 { 00401 // Generate two random points using the IDRFs. 00402 float X1, Y1, Z1; 00403 float X2, Y2, Z2; 00404 float pos1[3], pos2[3]; 00405 if (GenerateIRISRdmPos(pos1, alpha1, beta1) ) 00406 { 00407 Cerr("***** iProjectorIRIS::ProjectWithoutTOF() -> An error occurred while trying to generate random position with the IRIS projector" << endl); 00408 return 1; 00409 } 00410 00411 X1 = pos1[0]; 00412 Y1 = pos1[1]; 00413 Z1 = pos1[2]; 00414 00415 //X1 = event1[0]; 00416 //Y1 = event1[1]; 00417 //Z1 = event1[2]; 00418 00419 if (GenerateIRISRdmPos(pos2, alpha2, beta2) ) 00420 { 00421 Cerr("***** iProjectorIRIS::ProjectWithoutTOF() -> An error occurred while trying to generate random position with the IRIS projector" << endl); 00422 return 1; 00423 } 00424 00425 X2 = pos2[0]; 00426 Y2 = pos2[1]; 00427 Z2 = pos2[2]; 00428 00429 //X2 = event2[0]; 00430 //Y2 = event2[1]; 00431 //Z2 = event2[2]; 00432 00433 // Move the random points in the frame of reference of the scanner. 00434 float tmp_X, tmp_Y; 00435 tmp_X = X1; tmp_Y = Y1; 00436 X1 = tmp_X*orientation1[1]-tmp_Y*orientation1[0]+x1; //HERE 00437 Y1 = tmp_X*orientation1[0]+tmp_Y*orientation1[1]+y1; //HERE 00438 Z1 += z1; 00439 00440 tmp_X = X2; tmp_Y = Y2; 00441 X2 = tmp_X*orientation2[1]-tmp_Y*orientation2[0]+x2; //HERE 00442 Y2 = tmp_X*orientation2[0]+tmp_Y*orientation2[1]+y2; //HERE 00443 Z2 += z2; 00444 00445 00446 // SIDDON ALGORITHM 00447 00448 // ************************************** 00449 // STEP 1: LOR length calculation 00450 // ************************************** 00451 double length_LOR = ap_ProjectionLine->GetLength(); 00452 00453 // ************************************** 00454 // STEP 2: Compute entrance voxel indexes 00455 // ************************************** 00456 double alphaFirst[] = { 0.0, 0.0, 0.0 }; 00457 double alphaLast[] = { 0.0, 0.0, 0.0 }; 00458 00459 double alphaMin = 0.0f, alphaMax = 1.0f; 00460 double delta_pos[] = { 00461 X2 - X1, // event2[ 0 ] - event1[ 0 ] 00462 Y2 - Y1, // event2[ 1 ] - event1[ 1 ] 00463 Z2 - Z1 // event2[ 2 ] - event1[ 2 ] 00464 }; 00465 00466 // Computation of alphaMin et alphaMax values (entrance and exit point of the ray) 00467 for( int i = 0; i < 3; ++i ) 00468 { 00469 if( delta_pos[ i ] != 0.0 ) 00470 { 00471 alphaFirst[i] = (-mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00472 alphaLast[i] = ( mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00473 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i])); 00474 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i])); 00475 } 00476 } 00477 00478 // if alphaMax is less than or equal to alphaMin no intersection 00479 // and return an empty buffer 00480 if( alphaMax <= alphaMin ) return 0; 00481 00482 // Now we have to find the indices of the particular plane 00483 // (iMin,iMax), (jMin,jMax), (kMin,kMax) 00484 int iMin = 0, iMax = 0; 00485 int jMin = 0, jMax = 0; 00486 int kMin = 0, kMax = 0; 00487 00488 // For the X-axis 00489 if( delta_pos[ 0 ] > 0.0f ) 00490 { 00491 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMin * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00492 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00493 } 00494 else if( delta_pos[ 0 ] < 0.0f ) 00495 { 00496 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMax * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00497 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00498 } 00499 if( delta_pos[ 0 ] == 0 ) 00500 { 00501 iMin = 1, iMax = 0; 00502 } 00503 00504 // For the Y-axis 00505 if( delta_pos[ 1 ] > 0 ) 00506 { 00507 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMin * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00508 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00509 } 00510 else if( delta_pos[ 1 ] < 0 ) 00511 { 00512 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMax * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00513 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00514 } 00515 else if( delta_pos[ 1 ] == 0 ) 00516 { 00517 jMin = 1, jMax = 0; 00518 } 00519 00520 // For the Z-axis 00521 if( delta_pos[ 2 ] > 0 ) 00522 { 00523 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMin * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00524 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00525 } 00526 else if( delta_pos[ 2 ] < 0 ) 00527 { 00528 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMax * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00529 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00530 } 00531 else if( delta_pos[ 2 ] == 0 ) 00532 { 00533 kMin = 1, kMax = 0; 00534 } 00535 00536 // Computing the last term n number of intersection 00537 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 ) 00538 + ( kMax - kMin + 1 ); 00539 00540 // Array storing the first alpha in X, Y and Z 00541 double alpha_XYZ[ 3 ] = { 1.0f, 1.0f, 1.0f }; 00542 00543 // Computing the first alpha in X 00544 if( delta_pos[ 0 ] > 0 ) 00545 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMin - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00546 else if( delta_pos[ 0 ] < 0 ) 00547 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMax - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00548 00549 // Computing the first alpha in Y 00550 if( delta_pos[ 1 ] > 0 ) 00551 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMin - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00552 else if( delta_pos[ 1 ] < 0 ) 00553 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMax - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00554 00555 // Computing the first alpha in Z 00556 if( delta_pos[ 2 ] > 0 ) 00557 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMin - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00558 else if( delta_pos[ 2 ] < 0 ) 00559 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMax - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00560 00561 // Computing the alpha updating 00562 double alpha_u[ 3 ] = { 00563 mp_sizeVox[0] / std::fabs( delta_pos[ 0 ] ), 00564 mp_sizeVox[1] / std::fabs( delta_pos[ 1 ] ), 00565 mp_sizeVox[2] / std::fabs( delta_pos[ 2 ] ) 00566 }; 00567 00568 // Computing the index updating 00569 INTNB index_u[ 3 ] = { 00570 (delta_pos[ 0 ] < 0) ? -1 : 1, 00571 (delta_pos[ 1 ] < 0) ? -1 : 1, 00572 (delta_pos[ 2 ] < 0) ? -1 : 1 00573 }; 00574 00575 // Check which alpha is the min/max and increment 00576 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00577 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00578 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00579 00580 // Computing the minimum value in the alpha_XYZ buffer 00581 double const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ], 00582 (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) ); 00583 00584 // Computing the first index of intersection 00585 double const alphaMid = ( min_alpha_XYZ + alphaMin ) / 2.0f; 00586 INTNB index_ijk[ 3 ] = { 00587 1 + (int)( ( event1[ 0 ] + alphaMid * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ), 00588 1 + (int)( ( event1[ 1 ] + alphaMid * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ), 00589 1 + (int)( ( event1[ 2 ] + alphaMid * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ) 00590 }; 00591 00592 INTNB const w = mp_nbVox[ 0 ]; 00593 INTNB const wh = w * mp_nbVox[ 1 ]; 00594 00595 // Loop over the number of plane to cross 00596 double alpha_c = alphaMin; 00597 FLTNB coeff = 0.0f; 00598 INTNB numVox = 0; 00599 for( int nP = 0; nP < n - 1; ++nP ) 00600 { 00601 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] ) 00602 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) ) 00603 { 00604 // Storing values 00605 if( ( alpha_XYZ[ 0 ] >= alphaMin ) 00606 && ( index_ijk[ 0 ] - 1 >= 0 ) 00607 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00608 && ( index_ijk[ 1 ] - 1 >= 0 ) 00609 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00610 && ( index_ijk[ 2 ] - 1 >= 0 ) 00611 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00612 { 00613 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR; 00614 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00615 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/m_nbLinesPerLOR); // TODO Normalizations according to the number of lines are performed in these lines. Perhaps we should do this at a upper level (vProjector) 00616 } 00617 00618 // Increment values 00619 alpha_c = alpha_XYZ[ 0 ]; 00620 alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00621 index_ijk[ 0 ] += index_u[ 0 ]; 00622 } 00623 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] ) 00624 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) ) 00625 { 00626 // Storing values 00627 if( ( alpha_XYZ[ 1 ] >= alphaMin ) 00628 && ( index_ijk[ 0 ] - 1 >= 0 ) 00629 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00630 && ( index_ijk[ 1 ] - 1 >= 0 ) 00631 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00632 && ( index_ijk[ 2 ] - 1 >= 0 ) 00633 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00634 { 00635 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR; 00636 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00637 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/m_nbLinesPerLOR); 00638 } 00639 00640 // Increment values 00641 alpha_c = alpha_XYZ[ 1 ]; 00642 alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00643 index_ijk[ 1 ] += index_u[ 1 ]; 00644 } 00645 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] ) 00646 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) ) 00647 { 00648 // Storing values 00649 if( ( alpha_XYZ[ 2 ] >= alphaMin ) 00650 && ( index_ijk[ 0 ] - 1 >= 0 ) 00651 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00652 && ( index_ijk[ 1 ] - 1 >= 0 ) 00653 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00654 && ( index_ijk[ 2 ] - 1 >= 0 ) 00655 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00656 { 00657 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR; 00658 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00659 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/m_nbLinesPerLOR); 00660 } 00661 00662 // Increment values 00663 alpha_c = alpha_XYZ[ 2 ]; 00664 alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00665 index_ijk[ 2 ] += index_u[ 2 ]; 00666 } 00667 } 00668 00669 } 00670 00671 return 0; 00672 } 00673 00674 // ===================================================================== 00675 // --------------------------------------------------------------------- 00676 // --------------------------------------------------------------------- 00677 // ===================================================================== 00678 00679 int iProjectorIRIS::ProjectWithTOFPos(int a_Projector, oProjectionLine* ap_ProjectionLine) 00680 { 00681 Cerr("***** iProjectorIRIS::ProjectWithTOFPos() -> Not yet implemented !" << endl); 00682 return 1; 00683 } 00684 00685 // ===================================================================== 00686 // --------------------------------------------------------------------- 00687 // --------------------------------------------------------------------- 00688 // ===================================================================== 00689 00690 int iProjectorIRIS::ProjectWithTOFBin(int a_Projector, oProjectionLine* ap_ProjectionLine) 00691 { 00692 Cerr("***** iProjectorIRIS::ProjectWithTOFBin() -> Not yet implemented !" << endl); 00693 return 1; 00694 } 00695 00696 // ===================================================================== 00697 // --------------------------------------------------------------------- 00698 // --------------------------------------------------------------------- 00699 // ===================================================================== 00700 00701 int iProjectorIRIS::ComputeIDRF_CDF(int a_angleId) 00702 { 00703 // Compute cumulated sum 00704 m2p_IDRF_CDFs[a_angleId][0] = mp_IDRF[0]; 00705 00706 for(int i = 1; i<m_nVoxXYZIDRF; ++i) 00707 m2p_IDRF_CDFs[a_angleId][i] = m2p_IDRF_CDFs[a_angleId][i-1] + mp_IDRF[i]; 00708 00709 // Normalization 00710 if(m2p_IDRF_CDFs[a_angleId][m_nVoxXYZIDRF-1] > 0) // check if not null 00711 for(int i = 0; i<m_nVoxXYZIDRF; ++i) 00712 m2p_IDRF_CDFs[a_angleId][i] /= m2p_IDRF_CDFs[a_angleId][m_nVoxXYZIDRF-1]; 00713 00714 return 0; 00715 } 00716 00717 // ===================================================================== 00718 // --------------------------------------------------------------------- 00719 // --------------------------------------------------------------------- 00720 // ===================================================================== 00721 00722 int iProjectorIRIS::GenerateIRISRdmPos(float ap_generatedPos[3], float a_alpha, float a_beta) 00723 { 00724 sRNG* p_RNG = sRNG::GetInstance(); 00725 00726 float rdm_pos[3]; 00727 00728 int alpha_id = (int) (fabsf(a_alpha)/m_stepAlphaAnglesIDRF); 00729 int beta_id = (int) (fabsf(a_beta) /m_stepBetaAnglesIDRF); 00730 00731 alpha_id = min(alpha_id, m_nAlphaAnglesIDRF-1); 00732 beta_id = min(beta_id , m_nBetaAnglesIDRF -1); 00733 00734 int id = FindGreaterValue(m2p_IDRF_CDFs[alpha_id*m_nBetaAnglesIDRF+beta_id], p_RNG->GenerateRdmNber(), m_nVoxXYZIDRF); 00735 00736 // IDRFs indices sorting are axial/transaxial/DOI, DOI_id being the smallest index. 00737 int axial_id = id / (m_nVoxDepthIDRF*m_nVoxTransaxialIDRF); 00738 int trans_id = (id-axial_id*m_nVoxTransaxialIDRF*m_nVoxDepthIDRF) / m_nVoxDepthIDRF; 00739 int DOI_id = id-axial_id*m_nVoxTransaxialIDRF*m_nVoxDepthIDRF - trans_id*m_nVoxDepthIDRF; 00740 00741 // Use random number on each axis to generate a random position on the IDRF voxel 00742 rdm_pos[0] = (DOI_id + p_RNG->GenerateRdmNber() - m_nVoxDepthIDRF/2.0f )*m_sizeVoxDepthIDRF; // CHECK : should choose surface central position by default when using IRIS, or adress the difference here 00743 rdm_pos[1] = (trans_id + p_RNG->GenerateRdmNber() - m_nVoxTransaxialIDRF/2.0f )*m_sizeVoxTransaxialIDRF; 00744 rdm_pos[2] = (axial_id + p_RNG->GenerateRdmNber() - m_nVoxAxialIDRF/2.0f )*m_sizeVoxAxialIDRF; 00745 00746 if(a_alpha<0) rdm_pos[1] = -rdm_pos[1]; 00747 if(a_beta<0) rdm_pos[2] = -rdm_pos[2]; 00748 00749 ap_generatedPos[0] = rdm_pos[0]; 00750 ap_generatedPos[1] = rdm_pos[1]; 00751 ap_generatedPos[2] = rdm_pos[2]; 00752 00753 return 0; 00754 } 00755 00756 // ===================================================================== 00757 // --------------------------------------------------------------------- 00758 // --------------------------------------------------------------------- 00759 // ===================================================================== 00760 00761 int iProjectorIRIS::FindGreaterValue(float *ap_val, float a_key, int a_maxValue, int a_minStart, int a_maxStart) 00762 { 00763 int min = a_minStart, max, mid; 00764 00765 max = a_maxStart ? a_maxStart : a_maxValue; 00766 00767 int cpt=0; 00768 00769 while (min < max) 00770 { 00771 mid = (min + max) >> 1; //(min+max)/2 00772 if (a_key > ap_val[mid]) 00773 { 00774 min = mid + 1; 00775 } 00776 else 00777 { 00778 max = mid; 00779 } 00780 00781 cpt++; 00782 } 00783 00784 return min; 00785 }