![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 /* 00002 Implementation of class file gOptions.hh 00003 00004 - separators: X 00005 - doxygen: X 00006 - default initialization: none 00007 - CASTOR_DEBUG: none 00008 - CASTOR_VERBOSE: none (no verbose member variable) 00009 */ 00010 00011 00020 #include "gVariables.hh" 00021 #include "gDataConversionUtilities.hh" 00022 #include <iomanip> 00023 00024 /* 00025 Miscelleanous functionalities 00026 */ 00027 00028 // ===================================================================== 00029 // --------------------------------------------------------------------- 00030 // --------------------------------------------------------------------- 00031 // ===================================================================== 00032 /* 00033 \fn CheckGATECommand 00034 \param a_key: string containing a GATE command to recover 00035 \param a_line : string containing the line to check 00036 \brief Check if the line contains the provided GATE command. In this case, \n 00037 parse the line and returns the values in a vector of strings 00038 \details Values are converted in mm if 'cm' is found in the line 00039 \return the vector of strings containing the elements of the line. 00040 */ 00041 vector<string> CheckGATECommand(const string& a_key, const string& a_line) 00042 { 00043 vector<string> values; 00044 00045 // cut any part after comment symbol 00046 string line = a_line; 00047 00048 line = (line.find("#") != string::npos) ? 00049 line.substr(0, line.find_first_of("#")) : 00050 line; 00051 00052 size_t foundAdress = line.find(a_key); 00053 00054 if (foundAdress != string::npos) 00055 { 00056 values = Split(a_line); 00057 values.erase (values.begin()); 00058 } 00059 00060 ConvertValuesTomm(values); 00061 return values; 00062 } 00063 00064 00065 00066 // ===================================================================== 00067 // --------------------------------------------------------------------- 00068 // --------------------------------------------------------------------- 00069 // ===================================================================== 00070 /* 00071 \fn split 00072 \param a_line : string to split 00073 \brief Split the line provided in parameter into a vector of strings (separator is blankspace) 00074 \return vector of string containing the splitted elements of the line 00075 */ 00076 vector<string> Split(string a_line) 00077 { 00078 string buf; 00079 stringstream ss(a_line); 00080 00081 vector<string> tokens; 00082 00083 while (ss >> buf) 00084 tokens.push_back(buf); 00085 00086 return tokens; 00087 } 00088 00089 00090 00091 // ===================================================================== 00092 // --------------------------------------------------------------------- 00093 // --------------------------------------------------------------------- 00094 // ===================================================================== 00095 /* 00096 \fn ConvertValuesTomm 00097 \param ap_v : vector of strings 00098 \brief Check if the vector of strings passed in parameter contains the 'cm' unit \n 00099 In this case, convert all its values in mm 00100 */ 00101 void ConvertValuesTomm(vector<string>& ap_v) 00102 { 00103 // loop on values 00104 for(uint16_t i=0; i < ap_v.size(); i++) 00105 // check if values were provided in cm 00106 if (ap_v[i] == "cm") 00107 // Convert all values to mm (skip command key) 00108 for (int j=0; j<i; j++) 00109 { 00110 float val; 00111 ConvertFromString(ap_v[j], &val); 00112 ap_v[j] = toString(10 * val); 00113 } 00114 } 00115 00116 00117 // ===================================================================== 00118 // --------------------------------------------------------------------- 00119 // --------------------------------------------------------------------- 00120 // ===================================================================== 00121 /* 00122 \fn WriteVector 00123 \param file : the output file 00124 \param a_key : key to write 00125 \param a_vals : vector containing the key values 00126 \brief Write the key and its values in the file provided in parameter 00127 \return 0 if success, positive value otherwise 00128 */ 00129 template <class T> 00130 int WriteVector(ofstream& file, const string& a_key, vector <T> a_vals) 00131 { 00132 int n = a_vals.size(); 00133 00134 if (n > 0) 00135 { 00136 file << a_key ; 00137 for (int i=0; i < n; i++) 00138 { 00139 stringstream ss; 00140 ss << a_vals[i]; 00141 if (i == n-1) 00142 file << ss.str() << endl; 00143 else 00144 file << ss.str() << ","; 00145 } 00146 } 00147 else 00148 file << a_key << "0" <<endl; 00149 00150 return 0; 00151 } 00152 00153 template int WriteVector(ofstream& file, const string& a_key, vector <double> a_vals); 00154 00155 00156 00157 // ===================================================================== 00158 // --------------------------------------------------------------------- 00159 // --------------------------------------------------------------------- 00160 // ===================================================================== 00161 /* 00162 \fn WriteVector 00163 \param file : the output file 00164 \param a_key : key to write 00165 \param a_vals : vector containing the key values 00166 \brief Write the key and its values in the file provided in parameter 00167 \return 0 if success, positive value otherwise 00168 */ 00169 int WriteVector(ofstream& file, const string& a_key, vector <string> a_vals) 00170 { 00171 int n = a_vals.size(); 00172 00173 if (n > 0) 00174 { 00175 file << a_key ; 00176 for (int i=0; i < n; i++) 00177 { 00178 if (i == n-1) 00179 file << a_vals[i] << endl; 00180 else 00181 file << a_vals[i] << ","; 00182 } 00183 } 00184 else 00185 file << a_key << "0" <<endl; 00186 00187 return 0; 00188 } 00189 00190 00191 00192 // ===================================================================== 00193 // --------------------------------------------------------------------- 00194 // --------------------------------------------------------------------- 00195 // ===================================================================== 00196 /* 00197 \fn WriteVector 00198 \param file : the output file 00199 \param a_key : key to write 00200 \param a_vals : vector containing the key values 00201 \brief Write the key and its values contained in a 2 level vector of strings in the file provided in parameter 00202 \return 0 if success, positive value otherwise 00203 */ 00204 int WriteVector(ofstream& file, const string& a_key, vector <vector<string> > a_vals) 00205 { 00206 int n = a_vals.size(); 00207 file << a_key; 00208 for (int i=0; i < n; i++) 00209 for (int j=0; j < 3; j++) 00210 if (i == n-1 && j == 2) 00211 file << a_vals[i][j] << endl; 00212 else 00213 file << a_vals[i][j] << ","; 00214 00215 return 0; 00216 } 00217 00218 00219 00220 00221 // ===================================================================== 00222 // --------------------------------------------------------------------- 00223 // --------------------------------------------------------------------- 00224 // ===================================================================== 00225 /* 00226 \fn GetGATEMacFiles(const string& a_pathMac, vector<string> ap_pathToMacFiles) 00227 \param a_pathMac : path to a GATE main macro file 00228 \param ap_pathToMacFiles : array containing the paths of macro files 00229 \brief Extract the paths to each macro file contained in the main macro file. 00230 \return 0 if success, positive value otherwise 00231 */ 00232 int GetGATEMacFiles(const string& a_pathMac, vector<string> &ap_pathToMacFiles) 00233 { 00234 ifstream mac_file(a_pathMac, ios::in); 00235 00236 if(mac_file) 00237 { 00238 string quickLine; 00239 00240 while(getline(mac_file, quickLine)) 00241 { 00242 vector <string> values; 00243 00244 values = CheckGATECommand("/control/execute", quickLine); 00245 00246 if (values.size()>0) 00247 ap_pathToMacFiles.push_back(GetPathOfFile(a_pathMac)+values[0]); 00248 } 00249 } 00250 else 00251 { 00252 Cerr("***** GetGATEMacFiles() -> Couldn't open mac file "<< a_pathMac << " !" << endl); 00253 return -1; 00254 } 00255 00256 mac_file.close(); 00257 return 0; 00258 } 00259 00260 00261 00262 00263 // ===================================================================== 00264 // --------------------------------------------------------------------- 00265 // --------------------------------------------------------------------- 00266 // ===================================================================== 00267 /* 00268 \fn GetGATESystemType 00269 \param a_pathMac : path to a GATE macro file 00270 \brief Read a GATE macro file and identify the system type from the 'gate/systems/' command lines 00271 \return system type, as described by the macro GATE_SYS_TYPE 00272 */ 00273 int GetGATESystemType(const string& a_pathMac) 00274 { 00275 int system_type = GATE_SYS_UNKNOWN; 00276 00277 // Recover path to all macro files from the main mac file 00278 vector<string> path_mac_files; 00279 path_mac_files.push_back(a_pathMac); 00280 00281 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 00282 { 00283 Cerr("***** GetGATESystemType() -> Error occured when trying to recover paths to GATE macro files !" << endl); 00284 return 1; 00285 } 00286 00287 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 00288 cout << f << " : " << path_mac_files[f] << endl; 00289 00290 // Loop on all macro files 00291 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 00292 { 00293 ifstream mac_file(path_mac_files[f], ios::in); 00294 00295 if(mac_file) 00296 { 00297 string line; 00298 while(getline(mac_file, line)) 00299 { 00300 // Check a 'gate/systems' command line 00301 if(line.find("/gate/systems/") != string::npos ) 00302 { 00303 if(line.find("/gate/systems/ecat") != string::npos ) 00304 system_type = GATE_SYS_ECAT; 00305 00306 if(line.find("/gate/systems/cylindricalPET") != string::npos ) 00307 system_type = GATE_SYS_CYLINDRICAL; 00308 00309 if(line.find("/gate/systems/SPECThead") != string::npos ) 00310 system_type = GATE_SYS_SPECT; 00311 00312 if(line.find("/gate/systems/OPET") != string::npos || 00313 line.find("/gate/systems/CTSCANNER") != string::npos || 00314 line.find("/gate/systems/CPET") != string::npos || 00315 line.find("/gate/systems/ecatAccel") != string::npos || 00316 line.find("/gate/systems/OpticalSystem") != string::npos ) 00317 { 00318 Cerr("unsupported system detected (line = " << line <<") ! "<< endl); 00319 Cerr("supported systems for this script are cylindricalPET, SPECThead, and ecat" << endl); 00320 } 00321 } 00322 } 00323 } 00324 else 00325 { 00326 Cerr("***** GetGATESystemType() -> Error : Couldn't open mac file "<< path_mac_files[f] << " !" << endl); 00327 return -1; 00328 } 00329 00330 mac_file.close(); 00331 } 00332 00333 return system_type; 00334 } 00335 00336 00337 00338 00339 // ===================================================================== 00340 // --------------------------------------------------------------------- 00341 // --------------------------------------------------------------------- 00342 // ===================================================================== 00343 /* 00344 \fn GetGATEAliasesCylindrical(vector<string> path_mac_files 00345 string& rsector_name, 00346 string& module_name, 00347 string& submodule_name, 00348 string& crystal_name, 00349 vector<string>& layers_name, 00350 int vb ) 00351 \param path_mac_files 00352 \param rsector_name 00353 \param module_name 00354 \param submodule_name 00355 \param crystal_name 00356 \param layers_name 00357 \param vb 00358 \brief Loop over a list of path to GATE macro files passed in parameter 00359 to recover aliases of the different parts of a cylindricalPET 00360 \return 0 if success, positive value otherwise 00361 */ 00362 int GetGATEAliasesCylindrical(vector<string> path_mac_files, 00363 string& rsector_name, 00364 string& module_name, 00365 string& submodule_name, 00366 string& crystal_name, 00367 vector<string>& layers_name, 00368 int vb ) 00369 { 00370 int depth = 0; 00371 00372 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 00373 { 00374 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 00375 00376 if(mac_file) 00377 { 00378 string quickLine; 00379 while(getline(mac_file, quickLine)) 00380 { 00381 vector <string> values; 00382 00383 values = CheckGATECommand("/gate/systems/cylindricalPET/rsector/attach", quickLine); 00384 00385 if (values.size()>0) 00386 { 00387 rsector_name = values[0]; 00388 depth++; 00389 } 00390 00391 values = CheckGATECommand("/gate/systems/cylindricalPET/module/attach", quickLine); 00392 if (values.size()>0) 00393 { 00394 module_name = values[0]; 00395 depth++; 00396 } 00397 00398 values = CheckGATECommand("/gate/systems/cylindricalPET/submodule/attach", quickLine); 00399 if (values.size()>0) 00400 { 00401 submodule_name = values[0]; 00402 depth++; 00403 } 00404 00405 values = CheckGATECommand("/gate/systems/cylindricalPET/crystal/attach", quickLine); 00406 if (values.size()>0) 00407 { 00408 crystal_name = values[0]; 00409 depth++; 00410 } 00411 00412 // layers 00413 for(int l=0 ; l<GATE_NB_MAX_LAYERS ; l++) 00414 { 00415 stringstream ss; 00416 ss << l; 00417 values = CheckGATECommand("/gate/systems/cylindricalPET/layer"+ss.str()+"/attach", quickLine); 00418 if (values.size()>0) 00419 { 00420 layers_name.push_back(values[0]); 00421 depth++; 00422 } 00423 } 00424 00425 } 00426 } 00427 else 00428 { 00429 Cerr("***** GetGATEAliasesCylindrical()->Couldn't open mac file "<< path_mac_files[f] << " !" << endl); 00430 return 1; 00431 } 00432 00433 mac_file.close(); 00434 } 00435 00436 00437 // Check we have the required number of elements for the system 00438 if (depth < 2) 00439 { 00440 Cout("***** GetGATEAliasesCylindrical() :: Error : Missing elements in the system architecture" << endl << 00441 " At least two of the following lines are required :" << endl << 00442 " - /gate/systems/cylindricalPET/rsector/attach" << endl << 00443 " - /gate/systems/cylindricalPET/module/attach" << endl << 00444 " - /gate/systems/cylindricalPET/submodule/attach" << endl << 00445 " - /gate/systems/cylindricalPET/crystal/attach" << endl << 00446 " - /gate/systems/cylindricalPET/layeri[i=0..3]/attach" << endl); 00447 return 1; 00448 } 00449 else 00450 { 00451 if(vb >= 2) 00452 { 00453 if(!rsector_name.empty()) Cout("Detected rsector container's name : " << rsector_name << endl); 00454 if(!module_name.empty()) Cout("Detected module container's name : " << module_name << endl); 00455 if(!submodule_name.empty()) Cout("Detected submodule container's name : " << submodule_name << endl); 00456 if(!crystal_name.empty()) Cout("Detected crystal container's name : " << crystal_name << endl); 00457 for(size_t l=0 ; l<layers_name.size() ; l++) 00458 if(!layers_name[l].empty()) Cout("Detected layer #"<< l << " container's name : " << layers_name[l] << endl); 00459 } 00460 } 00461 00462 return 0; 00463 } 00464 00465 00466 00467 00468 // ===================================================================== 00469 // --------------------------------------------------------------------- 00470 // --------------------------------------------------------------------- 00471 // ===================================================================== 00472 /* 00473 \fn GetGATEAliasesEcat(vector<string> path_mac_files, 00474 string& block_name, 00475 string& crystal_name, 00476 int vb ) 00477 \param path_mac_files 00478 \param block_name 00479 \param crystal_name 00480 \param vb 00481 \brief Loop over a list of path to GATE macro files passed in parameter 00482 to recover aliases of the different parts of an ecat system 00483 \return 0 if success, positive value otherwise 00484 */ 00485 int GetGATEAliasesEcat(vector<string> path_mac_files, 00486 string& block_name, 00487 string& crystal_name, 00488 int vb ) 00489 { 00490 int depth = 0; 00491 00492 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 00493 { 00494 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 00495 00496 if(mac_file) 00497 { 00498 // Reading the .mac file line by line and finding the names of the different parts of the architecture 00499 string quickLine; 00500 while(getline(mac_file, quickLine)) 00501 { 00502 vector <string> values; 00503 00504 values = CheckGATECommand("/gate/systems/ecat/block/attach", quickLine); 00505 if (values.size()>0) 00506 { 00507 block_name = values[0]; 00508 depth++; 00509 } 00510 00511 values = CheckGATECommand("/gate/systems/ecat/crystal/attach", quickLine); 00512 if (values.size()>0) 00513 { 00514 crystal_name = values[0]; 00515 depth++; 00516 } 00517 } 00518 } 00519 else 00520 { 00521 Cerr("***** GetGATEAliasesEcat()-> Couldn't open mac file "<< path_mac_files[f].c_str()<< " !" << endl); 00522 return 1; 00523 } 00524 } 00525 00526 00527 // Check we have the required number of elements for the system 00528 if (depth < 2) 00529 { 00530 Cerr("***** GetGATEAliasesEcat() :: Error : Missing elements in the system architecture" << endl 00531 << " The following lines are required :" << endl 00532 << " - /gate/systems/ecat/block/attach" << endl 00533 << " - /gate/systems/ecat/crystal/attach" << endl); 00534 return 1; 00535 } 00536 else 00537 { 00538 if(vb >= 2) 00539 Cout("First container's name (usually block) is : " << block_name << endl 00540 << "Second container's name (usually crystal) is : " << crystal_name << endl); 00541 } 00542 00543 return 0; 00544 } 00545 00546 00547 00548 00549 00550 // ===================================================================== 00551 // --------------------------------------------------------------------- 00552 // --------------------------------------------------------------------- 00553 // ===================================================================== 00554 /* 00555 \fn GetGATEAliasesSPECT(vector<string> path_mac_files, 00556 string& base_name, 00557 string& crystal_name, 00558 string& pixel_name, 00559 int vb ) 00560 \param path_mac_files 00561 \param base_name 00562 \param crystal_name 00563 \param pixel_name 00564 \param vb 00565 \brief Loop over a list of path to GATE macro files passed in parameter 00566 to recover aliases of the different parts of an ecat system 00567 \return 0 if success, positive value otherwise 00568 */ 00569 int GetGATEAliasesSPECT(vector<string> path_mac_files, 00570 string& base_name, 00571 string& crystal_name, 00572 string& pixel_name, 00573 int vb ) 00574 { 00575 int depth = 0; 00576 00577 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 00578 { 00579 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 00580 00581 if(mac_file) 00582 { 00583 // Reading the .mac file line by line and finding the names of the different parts of the architecture 00584 string quickLine; 00585 00586 while(getline(mac_file, quickLine)) 00587 { 00588 vector <string> values; 00589 00590 values = CheckGATECommand("/gate/systems/SPECThead/base/attach", quickLine); 00591 if (values.size()>0) 00592 { 00593 base_name = values[0]; 00594 depth++; 00595 } 00596 00597 values = CheckGATECommand("/gate/systems/SPECThead/crystal/attach", quickLine); 00598 if (values.size()>0) 00599 { 00600 crystal_name = values[0]; 00601 depth++; 00602 } 00603 00604 values = CheckGATECommand("/gate/systems/SPECThead/pixel/attach", quickLine); 00605 if (values.size()>0) 00606 { 00607 pixel_name = values[0]; 00608 depth++; 00609 } 00610 } 00611 00612 } 00613 else 00614 { 00615 Cerr("***** GetGATEAliasesSPECT()-> Couldn't open mac file "<< path_mac_files[f].c_str()<< " !" << endl); 00616 return 1; 00617 } 00618 } 00619 00620 // Check we have the required number of elements for the system 00621 if (depth < 1) 00622 { 00623 Cerr("***** GetGATEAliasesSPECT() :: Error : Missing elements in the system architecture" << endl 00624 << " The following line is required :" << endl 00625 << " - /gate/systems/SPECThead/crystal/attach" << endl); 00626 return 1; 00627 } 00628 else 00629 { 00630 if(vb >= 2) 00631 Cout("Crystal container's name is : " << crystal_name << endl); 00632 } 00633 00634 return 0; 00635 } 00636 00637 00638 00639 00640 // ===================================================================== 00641 // --------------------------------------------------------------------- 00642 // --------------------------------------------------------------------- 00643 // ===================================================================== 00644 /* 00645 \fn ConvertIDecat 00646 \param nBlocksPerRing 00647 \param nBlocksLine 00648 \param nCrystalsTransaxial 00649 \param nCrystalsAxial 00650 \param crystalID 00651 \param blockID 00652 \brief Compute a CASToR crystal index of a GATE ecat system from its indexes (block/crystal) and the system layout 00653 \return CASToR crystal index 00654 */ 00655 uint32_t ConvertIDecat( int32_t nBlocksPerRing, 00656 int32_t nBlocksLine, 00657 int32_t nCrystalsTransaxial, 00658 int32_t nCrystalsAxial, 00659 int32_t crystalID, 00660 int32_t blockID) 00661 { 00662 int32_t nCrystalsPerRing = nBlocksPerRing * nCrystalsTransaxial; 00663 00664 int32_t ringID = (int32_t)( blockID/nBlocksPerRing ) * nCrystalsAxial 00665 + (int32_t)( crystalID/nCrystalsTransaxial ); 00666 00667 int32_t castorID = nCrystalsPerRing * ringID 00668 + nCrystalsTransaxial*( blockID % nBlocksPerRing ) 00669 + crystalID % nCrystalsTransaxial; 00670 00671 return castorID; 00672 } 00673 00674 00675 00676 // ===================================================================== 00677 // --------------------------------------------------------------------- 00678 // --------------------------------------------------------------------- 00679 // ===================================================================== 00680 /* 00681 \fn ConvertIDSPECTRoot1 00682 \param a_headID : head index as provided in the root file 00683 \param a_rotAngle : rotation angle (deg) as provided in the root file 00684 \param a_angStep : angular step (deg) computed from the macro file 00685 \param a_nProjectionsByHead : total number of projections for each head 00686 \brief Compute a CASToR projection index of a GATE SPECThead system 00687 \return CASToR crystal index 00688 */ 00689 uint32_t ConvertIDSPECTRoot1( int32_t a_headID, 00690 float_t a_rotAngle, 00691 float_t a_angStep, 00692 uint32_t a_nProjectionsByHead) 00693 { 00694 // Compute angular index from the angle position of the head and the angular step 00695 int32_t angID = round(a_rotAngle/a_angStep); 00696 00697 // Get final index for this head 00698 int32_t castorID = a_headID*a_nProjectionsByHead + angID; 00699 00700 return castorID; 00701 } 00702 00703 00704 00705 00706 // ===================================================================== 00707 // --------------------------------------------------------------------- 00708 // --------------------------------------------------------------------- 00709 // ===================================================================== 00710 /* 00711 \fn ConvertIDSPECTRoot2 00712 \param a_nbSimulatedPixels 00713 \param a_nPixTrs 00714 \param a_nPixAxl 00715 \param a_headID 00716 \param a_crystalID 00717 \param a_pixelID 00718 \param a_rotAngle 00719 \param a_headAngPitch 00720 \param a_crystalSizeAxl 00721 \param a_crystalSizeTrs 00722 \param a_gPosX 00723 \param a_gPosY 00724 \param a_gPosZ 00725 \brief Compute a CASToR crystal index of a GATE SPECThead system 00726 \return CASToR crystal index. Return a higher number than the max index if error 00727 */ 00728 uint32_t ConvertIDSPECTRoot2( uint32_t a_nbSimulatedPixels, 00729 uint32_t a_nPixTrs, 00730 uint32_t a_nPixAxl, 00731 int32_t a_headID, 00732 int32_t a_crystalID, 00733 int32_t a_pixelID, 00734 float_t a_rotAngle, 00735 float_t a_headAngPitch, 00736 float_t a_crystalSizeAxl, 00737 float_t a_crystalSizeTrs, 00738 float_t a_gPosX, 00739 float_t a_gPosY, 00740 float_t a_gPosZ) 00741 { 00742 int32_t castorID = 0; 00743 00744 // we have a pixel matrix, just return the pixelID 00745 if (a_nbSimulatedPixels > 1) 00746 { 00747 castorID = a_pixelID; 00748 } 00749 // Compute the pixelID according to the global XYZ positions in the crystal 00750 // Cheap implementation. Does not take the DOI into account (would depend on collimator) 00751 else 00752 { 00753 // Compute pixel sizes 00754 FLTNB sizePixTrs = a_crystalSizeTrs/a_nPixTrs; 00755 FLTNB sizePixAxl = a_crystalSizeAxl/a_nPixAxl; 00756 00757 // Compute axial ID 00758 uint32_t axialID = (int32_t)(( a_gPosZ + a_nPixAxl/2*sizePixAxl) / sizePixAxl); 00759 00760 // Compute transaxial ID 00761 // Compute head position angle ( in GATE referential) 00762 float_t ang = a_headID*a_headAngPitch + a_rotAngle; 00763 00764 // Get transaxial position 00765 float_t sin_a = sin(-ang*M_PI/180); 00766 float_t cos_a = cos(-ang*M_PI/180); 00767 float_t trs_pos = a_gPosX*sin_a + a_gPosY*cos_a ; 00768 00769 // Compute transaxial ID 00770 uint32_t transID = (int32_t)(( trs_pos + a_nPixTrs/2*sizePixTrs) / sizePixTrs); 00771 00772 if(axialID >= 0 && axialID < a_nPixAxl && 00773 transID >= 0 && transID < a_nPixTrs ) 00774 { 00775 castorID = axialID*a_nPixAxl + transID; 00776 } 00777 else // return more than the max crystal ID to check for errors 00778 { 00779 castorID = a_nPixAxl*a_nPixTrs; 00780 } 00781 } 00782 00783 return castorID; 00784 } 00785 00786 00787 00788 00789 // ===================================================================== 00790 // --------------------------------------------------------------------- 00791 // --------------------------------------------------------------------- 00792 // ===================================================================== 00793 /* 00794 \fn ConvertIDcylindrical 00795 \param nRsectorsPerRing 00796 \param nModulesTransaxial 00797 \param nModulesAxial 00798 \param nSubmodulesTransaxial 00799 \param nSubmodulesAxial 00800 \param nCrystalsTransaxial 00801 \param nCrystalsAxial 00802 \param nLayers 00803 \param nb_crystal_per_layer 00804 \param nLayersRptTransaxial 00805 \param nLayersRptAxial 00806 \param layerID 00807 \param crystalID 00808 \param submoduleID 00809 \param moduleID 00810 \param rsectorID 00811 \brief Compute a CASToR crystal index of a GATE cylindricalPET system from its indexes (rsector/module/submodule/crystal) and the system layout 00812 \return CASToR crystal index 00813 */ 00814 uint32_t ConvertIDcylindrical(uint32_t nRsectorsPerRing, 00815 uint32_t nModulesTransaxial, 00816 uint32_t nModulesAxial, 00817 uint32_t nSubmodulesTransaxial, 00818 uint32_t nSubmodulesAxial, 00819 uint32_t nCrystalsTransaxial, 00820 uint32_t nCrystalsAxial, 00821 uint8_t nLayers, 00822 uint32_t* nCrystalPerLayer, 00823 vector<uint32_t> nLayersRptTransaxial, 00824 vector<uint32_t> nLayersRptAxial, 00825 int32_t layerID, 00826 int32_t crystalID, 00827 int32_t submoduleID, 00828 int32_t moduleID, 00829 int32_t rsectorID) 00830 { 00831 // Castor ID definition 00832 uint32_t castorID = 0; 00833 uint8_t layer = 0; 00834 00835 if(nLayersRptTransaxial.size()==0 && nLayersRptAxial.size()==0) 00836 { 00837 // layerID represents the actual layer level 00838 00839 // add the number of crystals contained in previous layers as 00840 // CASToR indexes all crystals of a layer ring before the next layer 00841 layer = layerID; 00842 00843 for(int l=0 ; l<layer; l++) 00844 castorID += nCrystalPerLayer[l]; 00845 00846 int32_t nTrsCrystalsPerSubmodule = nCrystalsTransaxial; 00847 int32_t nTrsCrystalsPerModule = nTrsCrystalsPerSubmodule * nSubmodulesTransaxial; 00848 int32_t nTrsCrystalsPerRsector = nTrsCrystalsPerModule * nModulesTransaxial; 00849 int32_t nCrystalsPerRing = nTrsCrystalsPerRsector * nRsectorsPerRing; 00850 00851 int32_t ringID = (int32_t)(moduleID/nModulesTransaxial) * nSubmodulesAxial * nCrystalsAxial 00852 + (int32_t)(submoduleID/nSubmodulesTransaxial) * nCrystalsAxial 00853 + (int32_t)(crystalID/nCrystalsTransaxial); 00854 00855 castorID += nCrystalsPerRing * ringID 00856 + nTrsCrystalsPerRsector * rsectorID 00857 + nTrsCrystalsPerModule * ( moduleID % nModulesTransaxial ) 00858 + nTrsCrystalsPerSubmodule * ( submoduleID % nSubmodulesTransaxial ) 00859 + crystalID % nCrystalsTransaxial; 00860 } 00861 00862 else 00863 { 00864 // layerID represents a crystal layer element 00865 00866 // Get the total number of crystals in the first layer 00867 uint32_t sum_crystals = nLayersRptTransaxial[layer] 00868 * nLayersRptAxial[layer]; 00869 00870 // Get the layer which the crystal belongs to 00871 while (layerID >= (int32_t)sum_crystals) 00872 { 00873 layer++; 00874 sum_crystals += nLayersRptTransaxial[layer] 00875 * nLayersRptAxial[layer]; 00876 } 00877 00878 // add the number of crystals contained in previous layers as 00879 // CASToR indexes all crystals of a layer ring before the next layer 00880 for(int l=0 ; l<layer ; l++) 00881 castorID += nCrystalPerLayer[l]; 00882 00883 int32_t nTrsCrystalsPerSubmodule = nCrystalsTransaxial * nLayersRptTransaxial[layer]; 00884 int32_t nTrsCrystalsPerModule = nTrsCrystalsPerSubmodule * nSubmodulesTransaxial; 00885 int32_t nTrsCrystalsPerRsector = nTrsCrystalsPerModule * nModulesTransaxial; 00886 int32_t nCrystalsPerRing = nTrsCrystalsPerRsector * nRsectorsPerRing; 00887 00888 int32_t ringID = (int32_t)(moduleID/nModulesTransaxial) * nSubmodulesAxial * nCrystalsAxial * nLayersRptAxial[layer] 00889 + (int32_t)(submoduleID/nSubmodulesTransaxial) * nCrystalsAxial * nLayersRptAxial[layer] 00890 + (int32_t)(crystalID/nCrystalsTransaxial) * nLayersRptAxial[layer]; 00891 00892 if(nLayersRptTransaxial.empty() || nLayersRptAxial.empty() ) 00893 ringID += (int32_t)(layerID/nLayersRptTransaxial[layer]); 00894 00895 castorID += nCrystalsPerRing * ringID 00896 + nTrsCrystalsPerRsector * rsectorID 00897 + nTrsCrystalsPerModule * ( moduleID % nModulesTransaxial ) 00898 + nTrsCrystalsPerSubmodule * ( submoduleID % nSubmodulesTransaxial ) 00899 + crystalID % nCrystalsTransaxial 00900 + layerID % nLayersRptTransaxial[layer]; 00901 00902 } 00903 00904 00905 return castorID; 00906 } 00907 00908 00909 00910 // ===================================================================== 00911 // --------------------------------------------------------------------- 00912 // --------------------------------------------------------------------- 00913 // ===================================================================== 00914 /* 00915 \fn ComputeKindGATEEvent 00916 \param eventID1 00917 \param eventID2 00918 \param comptonPhantom1 00919 \param comptonPhantom2 00920 \param rayleighPhantom1 00921 \param rayleighPhantom2 00922 \brief Determine kind of a given coincidence event, from its attributes 00923 \return kind of the coincidence : unknown (=0, default), true(=1), single scat(=2), multiple scat(=3), random(=4)) (default value =0) 00924 */ 00925 int ComputeKindGATEEvent(uint32_t eventID1, uint32_t eventID2, 00926 int comptonPhantom1, int comptonPhantom2, 00927 int rayleighPhantom1, int rayleighPhantom2) 00928 { 00929 if (eventID1 != eventID2) 00930 //random 00931 return KIND_RDM; 00932 else 00933 { 00934 if (comptonPhantom1 == 0 && comptonPhantom2 == 0 && 00935 rayleighPhantom1 == 0 && rayleighPhantom2 == 0) 00936 //true 00937 return KIND_TRUE; 00938 else 00939 { 00940 if (comptonPhantom1 == 1 || comptonPhantom2 == 1 || 00941 rayleighPhantom1 == 1 || rayleighPhantom2 == 1) 00942 //single scat 00943 return KIND_SCAT; 00944 if (comptonPhantom1 > 1 || comptonPhantom2 > 1 || 00945 rayleighPhantom1 > 1 || rayleighPhantom2 > 1) 00946 //multiple scat 00947 return KIND_MSCAT; 00948 } 00949 } 00950 // unknown 00951 return KIND_UNKNOWN; 00952 } 00953 00954 00955 00956 // ===================================================================== 00957 // --------------------------------------------------------------------- 00958 // --------------------------------------------------------------------- 00959 // ===================================================================== 00960 /* 00961 \fn ReadMacCylindrical 00962 \param a_pathMac : path to a macro file 00963 \param nLayers : nb of crystal layers 00964 \param nCrystalsAxial : nb of axial crystals (in a submodule) 00965 \param nCrystalsTransaxial : nb of transaxial crystals (in a submodule) 00966 \param nLayersRptAxial : Array containing the number of axial crystals in each layer 00967 \param nLayersRptTransaxial : Array containing the number of transaxial crystals in each layer 00968 \param nSubmodulesAxial : nb of axial submodules (in a module) 00969 \param nSubmodulesTransaxial : nb of transaxial submodules (in a module) 00970 \param nModulesAxial : nb of axial modules (in a rsector) 00971 \param nModulesTransaxial : nb of transaxial modules (in a rsector) 00972 \param nRsectorsPerRing : nb of rsectors per ring 00973 \param start_time_ms : acquisition start time converted in ms 00974 \param duration_ms : acquisition duration converted in ms 00975 \param vb : verbosity 00976 \brief Recover informations about the scanner element of a cylindricalPET system and acquisition duration, from a GATE macro file 00977 \return 0 if success, positive value otherwise 00978 */ 00979 int ReadMacCylindrical(string a_pathMac, 00980 uint8_t &nLayers, 00981 uint32_t *nb_crystal_per_layer, 00982 uint32_t &nCrystalsTot, 00983 uint32_t &nCrystalsAxial, 00984 uint32_t &nCrystalsTransaxial, 00985 vector<uint32_t> &nLayersRptAxial, 00986 vector<uint32_t> &nLayersRptTransaxial, 00987 uint32_t &nSubmodulesAxial, 00988 uint32_t &nSubmodulesTransaxial, 00989 uint32_t &nModulesAxial, 00990 uint32_t &nModulesTransaxial, 00991 uint32_t &nRsectorsPerRing, 00992 uint32_t &start_time_ms, 00993 uint32_t &duration_ms, 00994 int vb) 00995 { 00996 vector<string> path_mac_files; 00997 path_mac_files.push_back(a_pathMac); 00998 00999 // Recover path to all macro files from the main mac file 01000 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 01001 { 01002 Cerr("***** GetGATESystemType() -> Error occured when trying to recover paths to GATE macro files !" << endl); 01003 return 1; 01004 } 01005 01006 string rsector_name = ""; 01007 string module_name = ""; 01008 string submodule_name = ""; 01009 string crystal_name = ""; 01010 vector <string> layers_name; 01011 bool is_rsector_Y_axis = false; 01012 01013 // Recover aliases of the different parts of the architecture 01014 if(GetGATEAliasesCylindrical(path_mac_files, rsector_name, module_name, submodule_name, crystal_name, layers_name, vb) ) 01015 { 01016 Cerr("***** GetGATESystemType() -> Error occured when trying to recover aliases for the elements of the cylindricalPET !" << endl); 01017 return 1; 01018 } 01019 01020 // Recover nb of detected layers 01021 nLayers = layers_name.size(); 01022 01023 // Loop to recover all other informations 01024 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 01025 { 01026 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 01027 01028 string line; 01029 double time_start =-1., 01030 time_stop =-1., 01031 time_slices =-1.; 01032 01033 while(getline(mac_file, line)) 01034 { 01035 vector <string> values; 01036 string kword =""; 01037 01038 kword = "/gate/"+rsector_name+"/placement/setTranslation"; 01039 values = CheckGATECommand(kword, line); 01040 01041 // Check where the first rsector has been created 01042 if (values.size()>0) 01043 { 01044 FLTNB rsector_pos_X =0., 01045 rsector_pos_Y =0.; 01046 01047 if(ConvertFromString(values[0], &rsector_pos_X) || 01048 ConvertFromString(values[1], &rsector_pos_Y) ) 01049 { 01050 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 01051 return 1; 01052 } 01053 01054 // Get the axis on which the first rsector is positionned 01055 if(rsector_pos_Y!=0) is_rsector_Y_axis = true; 01056 } 01057 01058 01059 kword = "/gate/"+rsector_name+"/ring/setRepeatNumber"; 01060 values = CheckGATECommand(kword, line); 01061 if (values.size()>0) 01062 { 01063 if(ConvertFromString( values[0] , &nRsectorsPerRing) ) 01064 { 01065 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01066 return 1; 01067 } 01068 } 01069 01070 kword = is_rsector_Y_axis ? 01071 "/gate/"+module_name+"/cubicArray/setRepeatNumberX": 01072 "/gate/"+module_name+"/cubicArray/setRepeatNumberY"; 01073 01074 values = CheckGATECommand(kword, line); 01075 if (values.size()>0) 01076 { 01077 if(ConvertFromString( values[0] , &nModulesTransaxial) ) 01078 { 01079 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01080 return 1; 01081 } 01082 } 01083 01084 kword = "/gate/"+module_name+"/cubicArray/setRepeatNumberZ"; 01085 values = CheckGATECommand(kword, line); 01086 if (values.size()>0) 01087 { 01088 if(ConvertFromString( values[0] , &nModulesAxial) ) 01089 { 01090 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01091 return 1; 01092 } 01093 } 01094 01095 kword = is_rsector_Y_axis ? 01096 "/gate/"+submodule_name+"/cubicArray/setRepeatNumberX": 01097 "/gate/"+submodule_name+"/cubicArray/setRepeatNumberY"; 01098 01099 values = CheckGATECommand(kword, line); 01100 if (values.size()>0) 01101 { 01102 if(ConvertFromString( values[0] , &nSubmodulesTransaxial) ) 01103 { 01104 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01105 return 1; 01106 } 01107 } 01108 01109 kword = "/gate/"+submodule_name+"/cubicArray/setRepeatNumberZ"; 01110 values = CheckGATECommand(kword, line); 01111 if (values.size()>0) 01112 { 01113 if(ConvertFromString( values[0] , &nSubmodulesAxial) ) 01114 { 01115 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01116 return 1; 01117 } 01118 } 01119 01120 01121 kword = is_rsector_Y_axis ? 01122 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberX": 01123 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberY"; 01124 01125 kword = "/gate/"+crystal_name+"/cubicArray/setRepeatNumberY"; 01126 values = CheckGATECommand(kword, line); 01127 if (values.size()>0) 01128 { 01129 if(ConvertFromString( values[0] , &nCrystalsTransaxial) ) 01130 { 01131 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01132 return 1; 01133 } 01134 } 01135 01136 kword = "/gate/"+crystal_name+"/cubicArray/setRepeatNumberZ"; 01137 values = CheckGATECommand(kword, line); 01138 if (values.size()>0) 01139 { 01140 if(ConvertFromString( values[0] , &nCrystalsAxial) ) 01141 { 01142 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01143 return 1; 01144 } 01145 } 01146 01147 // Check if there are any repeaters on crystal layers 01148 for(int l=0 ; l<nLayers ; l++) 01149 { 01150 kword = is_rsector_Y_axis ? 01151 "/gate/"+layers_name[l]+"/cubicArray/setRepeatNumberX": 01152 "/gate/"+layers_name[l]+"/cubicArray/setRepeatNumberY"; 01153 01154 values = CheckGATECommand(kword, line); 01155 if (values.size()>0) 01156 { 01157 int32_t val; 01158 if(ConvertFromString( values[0] , &val) ) 01159 { 01160 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01161 return 1; 01162 } 01163 nLayersRptTransaxial.push_back(val); 01164 } 01165 01166 kword = "/gate/"+layers_name[l]+"/cubicArray/setRepeatNumberZ"; 01167 values = CheckGATECommand(kword, line); 01168 if (values.size()>0) 01169 { 01170 int32_t val; 01171 if(ConvertFromString( values[0] , &val) ) 01172 { 01173 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01174 return 1; 01175 } 01176 nLayersRptAxial.push_back(val); 01177 } 01178 } 01179 01180 01181 kword = "/gate/application/setTimeStart"; 01182 values = CheckGATECommand(kword, line); 01183 if (values.size()>0) 01184 { 01185 if(ConvertFromString( values[0] , &time_start) ) 01186 { 01187 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01188 return 1; 01189 } 01190 01191 // Convert in ms 01192 if (values.size()>1) 01193 { 01194 if(values[1] == "s") time_start *= 1000; 01195 } 01196 else 01197 Cerr("***** dataConversionUtilities::readMacCylindrical()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01198 01199 start_time_ms = time_start; 01200 } 01201 01202 kword = "/gate/application/setTimeStop"; 01203 values = CheckGATECommand(kword, line); 01204 if (values.size()>0) 01205 { 01206 if(ConvertFromString( values[0] , &time_stop) ) 01207 { 01208 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01209 return 1; 01210 } 01211 01212 // Convert in ms 01213 if (values.size()>1) 01214 { 01215 if(values[1] == "s") time_stop *= 1000; 01216 } 01217 else 01218 Cerr("***** dataConversionUtilities::readMacCylindrical()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01219 } 01220 01221 kword = "/gate/application/addSlice"; 01222 values = CheckGATECommand(kword, line); 01223 if (values.size()>0) 01224 { 01225 double time_slice_tmp=0; 01226 01227 if(ConvertFromString( values[0] , &time_slice_tmp) ) 01228 { 01229 Cerr("***** dataConversionUtilities::readMacCylindrical()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01230 return 1; 01231 } 01232 01233 // Convert in ms 01234 if (values.size()>1) 01235 { 01236 if(values[1] == "s") time_slice_tmp *= 1000; 01237 } 01238 else 01239 Cerr("***** dataConversionUtilities::readMacCylindrical()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01240 01241 time_slices += time_slice_tmp; 01242 } 01243 } 01244 01245 // Compute duration in ms 01246 duration_ms = (time_slices>0) ? 01247 (uint32_t)(time_slices-time_start) : 01248 duration_ms; 01249 01250 duration_ms = (time_start>=0 && time_stop>=0) ? 01251 (uint32_t)(time_stop-time_start) : 01252 duration_ms ; 01253 01254 mac_file.close(); 01255 } 01256 01257 // Computing the total number of crystals in the scanner 01258 if(nLayers == 0) 01259 { 01260 nCrystalsTot = nRsectorsPerRing 01261 * nModulesTransaxial * nModulesAxial 01262 * nSubmodulesTransaxial * nSubmodulesAxial 01263 * nCrystalsTransaxial * nCrystalsAxial; 01264 } 01265 else 01266 for(int l=0 ; l<nLayers ; l++) 01267 { 01268 uint32_t nb_crystals_layer = nRsectorsPerRing 01269 * nModulesTransaxial * nModulesAxial 01270 * nSubmodulesTransaxial * nSubmodulesAxial 01271 * nCrystalsTransaxial * nCrystalsAxial; 01272 01273 // Add layer elements if repeaters have been used on layers 01274 if(nLayersRptTransaxial.size()>0 || nLayersRptAxial.size()>0 ) 01275 nb_crystals_layer *= nLayersRptTransaxial[l] * nLayersRptAxial[l]; 01276 01277 nb_crystal_per_layer[l] = nb_crystals_layer; 01278 01279 nCrystalsTot += nb_crystals_layer; 01280 } 01281 01282 01283 return 0; 01284 } 01285 01286 01287 01288 // ===================================================================== 01289 // --------------------------------------------------------------------- 01290 // --------------------------------------------------------------------- 01291 // ===================================================================== 01292 /* 01293 \fn ReadMacECAT 01294 \param a_pathMac : path to a macro file 01295 \param nCrystalsAxial : nb of axial crystals 01296 \param nCrystalsTransaxial : nb of transaxial crystals 01297 \param nBlocksLine : nb of axial blocks 01298 \param nBlocksPerRing : nb of blocks per ring 01299 \param start_time_ms : acquisition start time converted in ms 01300 \param duration_ms : acquisition duration converted in ms 01301 \param vb : verbosity 01302 \brief Recover informations about the scanner element of an ECAT system, and acquisition duration, from a GATE macro file 01303 \return 0 if success, positive value otherwise 01304 */ 01305 int ReadMacECAT( string a_pathMac, 01306 uint32_t &nCrystalsTot, 01307 uint32_t &nCrystalsAxial, 01308 uint32_t &nCrystalsTransaxial, 01309 uint32_t &nBlocksLine, 01310 uint32_t &nBlocksPerRing, 01311 uint32_t &start_time_ms, 01312 uint32_t &duration_ms, 01313 int vb) 01314 { 01315 // Recover path to all macro files from the main mac file 01316 vector<string> path_mac_files; 01317 path_mac_files.push_back(a_pathMac); 01318 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 01319 { 01320 Cerr("***** GetGATESystemType() -> Error occured when trying to recover paths to GATE macro files !" << endl); 01321 return 1; 01322 } 01323 01324 string block_name = "block"; 01325 string crystal_name = "crystal"; 01326 bool is_block_Y_axis = false; 01327 01328 // Recover aliases of the different parts of the architecture 01329 if(GetGATEAliasesEcat(path_mac_files, block_name, crystal_name, vb) ) 01330 { 01331 Cerr("***** GetGATESystemType() -> Error occured when trying to recover aliases for the elements of the ecat !" << endl); 01332 return 1; 01333 } 01334 01335 01336 // 2nd loop to recover all other informations 01337 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 01338 { 01339 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 01340 01341 string line; 01342 double time_start=-1., 01343 time_stop=-1., 01344 time_slices=-1.; 01345 01346 while(getline(mac_file, line)) 01347 { 01348 vector <string> values; 01349 string kword =""; 01350 01351 kword = "/gate/"+block_name+"/placement/setTranslation"; 01352 values = CheckGATECommand(kword, line); 01353 if (values.size()>0) 01354 { 01355 FLTNB block_pos_X=0., 01356 block_pos_Y=0.; 01357 01358 if(ConvertFromString(values[0], &block_pos_X) || 01359 ConvertFromString(values[1], &block_pos_Y) ) 01360 { 01361 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 01362 return 1; 01363 } 01364 01365 // Get the axis on which the first rsector is positionned 01366 if(block_pos_Y!=0) is_block_Y_axis = true; 01367 } 01368 01369 kword = "/gate/"+block_name+"/ring/setRepeatNumber"; 01370 values = CheckGATECommand(kword, line); 01371 if (values.size()>0) 01372 { 01373 if(ConvertFromString( values[0] , &nBlocksPerRing) ) 01374 { 01375 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01376 return 1; 01377 } 01378 } 01379 01380 kword = "/gate/"+block_name+"/linear/setRepeatNumber"; 01381 values = CheckGATECommand(kword, line); 01382 if (values.size()>0) 01383 { 01384 if(ConvertFromString( values[0] , &nBlocksLine) ) 01385 { 01386 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01387 return 1; 01388 } 01389 } 01390 01391 kword = is_block_Y_axis ? 01392 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberX": 01393 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberY"; 01394 01395 values = CheckGATECommand(kword, line); 01396 if (values.size()>0) 01397 { 01398 if(ConvertFromString( values[0] , &nCrystalsTransaxial) ) 01399 { 01400 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01401 return 1; 01402 } 01403 } 01404 01405 kword = "/gate/"+crystal_name+"/cubicArray/setRepeatNumberZ"; 01406 values = CheckGATECommand(kword, line); 01407 if (values.size()>0) 01408 { 01409 if(ConvertFromString( values[0] , &nCrystalsAxial) ) 01410 { 01411 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01412 return 1; 01413 } 01414 } 01415 01416 01417 kword = "/gate/application/setTimeStart"; 01418 values = CheckGATECommand(kword, line); 01419 if (values.size()>0) 01420 { 01421 if(ConvertFromString( values[0] , &time_start) ) 01422 { 01423 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01424 return 1; 01425 } 01426 01427 // Convert in ms 01428 if (values.size()>1) 01429 { 01430 if(values[1] == "s") time_start *= 1000; 01431 } 01432 else 01433 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01434 01435 start_time_ms = time_start; 01436 } 01437 01438 kword = "/gate/application/setTimeStop"; 01439 values = CheckGATECommand(kword, line); 01440 if (values.size()>0) 01441 { 01442 if(ConvertFromString( values[0] , &time_stop) ) 01443 { 01444 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01445 return 1; 01446 } 01447 01448 // Convert in ms 01449 if (values.size()>1) 01450 { 01451 if(values[1] == "s") time_stop *= 1000; 01452 } 01453 else 01454 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01455 } 01456 01457 kword = "/gate/application/addSlice"; 01458 values = CheckGATECommand(kword, line); 01459 if (values.size()>0) 01460 { 01461 double time_slice_tmp=0; 01462 01463 if(ConvertFromString( values[0] , &time_slice_tmp) ) 01464 { 01465 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01466 return 1; 01467 } 01468 01469 // Convert in ms 01470 if (values.size()>1) 01471 { 01472 if(values[1] == "s") time_slice_tmp *= 1000; 01473 } 01474 else 01475 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01476 01477 time_slices = time_slice_tmp; 01478 } 01479 } 01480 01481 // Compute duration in ms 01482 duration_ms = (time_slices>0) ? 01483 (uint32_t)(time_slices-time_start) : 01484 duration_ms; 01485 01486 duration_ms = (time_start>=0 && time_stop>=0) ? 01487 (uint32_t)(time_stop-time_start) : 01488 duration_ms ; 01489 01490 mac_file.close(); 01491 } 01492 01493 // Computing the total number of crystals in the scanner 01494 nCrystalsTot = nCrystalsTransaxial * nCrystalsAxial 01495 * nBlocksLine * nBlocksPerRing; 01496 01497 return 0; 01498 } 01499 01500 01501 01502 // ===================================================================== 01503 // --------------------------------------------------------------------- 01504 // --------------------------------------------------------------------- 01505 // ===================================================================== 01506 /* 01507 \fn ReadMacSPECT 01508 \param a_pathMac : path to a macro file 01509 \param distToDetector : distance between center of rotation and detector surface 01510 \param nHeads : nb of cameras 01511 \param nPixAxl : nb of transaxial pixels 01512 \param nPixTrs : nb of axial pixels 01513 \param crystalSizeAxl : crystal axial dimension 01514 \param crystalSizeTrs : crystal transaxial dimension 01515 \param head1stAngle : head first angle 01516 \param headAngPitch : angular pitch between heads 01517 \param headRotSpeed : angle between projection 01518 \param headRotDirection : rotation direction of the head (CW or CCW) 01519 \param start_time_ms : acquisition start time converted in ms 01520 \param duration_ms : acquisition duration converted in ms 01521 \param vb : verbosity 01522 \brief Recover informations about the scanner element of an ECAT system, and acquisition duration, from a GATE macro file 01523 \return 0 if success, positive value otherwise 01524 */ 01525 int ReadMacSPECT( string a_pathMac, 01526 float_t &a_distToDetector, 01527 uint32_t &a_nHeads, 01528 uint32_t &a_nPixAxl, 01529 uint32_t &a_nPixTrs, 01530 float_t &a_crystalSizeAxl, 01531 float_t &a_crystalSizeTrs, 01532 uint32_t &a_nProjectionsTot, 01533 uint32_t &a_nProjectionsByHead, 01534 float_t &a_head1stAngle, 01535 float_t &a_headAngPitch, 01536 float_t &a_headAngStepDeg, 01537 int &a_headRotDirection, 01538 uint32_t &a_start_time_ms, 01539 uint32_t &a_duration_ms, 01540 int vb) 01541 { 01542 // Recover path to all macro files from the main mac file 01543 vector<string> path_mac_files; 01544 path_mac_files.push_back(a_pathMac); 01545 01546 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 01547 { 01548 Cerr("***** GetGATESystemType() -> Error occured when trying to recover paths to GATE macro files !" << endl); 01549 return 1; 01550 } 01551 01552 string head_name = "SPECThead"; 01553 string crystal_name = "crystal"; 01554 string pixel_name = "pixel"; 01555 string head_orbit_name = ""; 01556 bool is_head_Y_axis = false; 01557 01558 // Recover aliases of the different parts of the architecture 01559 if(GetGATEAliasesSPECT(path_mac_files, head_name, crystal_name, pixel_name, vb) ) 01560 { 01561 Cerr("***** GetGATESystemType() -> Error occured when trying to recover aliases for the elements of the ecat !" << endl); 01562 return 1; 01563 } 01564 01565 // 2nd loop to recover all other informations 01566 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 01567 { 01568 ifstream mac_file(path_mac_files[f].c_str(), ios::in); 01569 01570 string line; 01571 double time_start=-1., 01572 time_stop=-1., 01573 time_slice=-1, 01574 time_slices=-1., 01575 time_slice_ms=-1., 01576 head_rot_speed =-1.; 01577 01578 01579 while(getline(mac_file, line)) 01580 { 01581 vector <string> values; 01582 string kword =""; 01583 01584 kword = "/gate/"+head_name+"/placement/setTranslation"; 01585 values = CheckGATECommand(kword, line); 01586 if (values.size()>0) 01587 { 01588 FLTNB head_pos_X=0., 01589 head_pos_Y=0.; 01590 01591 if(ConvertFromString(values[0], &head_pos_X) || 01592 ConvertFromString(values[1], &head_pos_Y) ) 01593 { 01594 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01595 return 1; 01596 } 01597 01598 // Get the axis on which the first rsector is positionned 01599 if(head_pos_Y!=0) is_head_Y_axis = true; 01600 01601 a_distToDetector = is_head_Y_axis ? abs(head_pos_Y) : abs(head_pos_X) ; 01602 } 01603 01604 01605 01606 kword = "/gate/"+head_name+"/ring/setRepeatNumber"; 01607 values = CheckGATECommand(kword, line); 01608 if (values.size()>0) 01609 { 01610 if(ConvertFromString(values[0], &a_nHeads)) 01611 { 01612 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01613 return 1; 01614 } 01615 } 01616 01617 kword = "/gate/"+head_name+"/ring/setFirstAngle"; 01618 values = CheckGATECommand(kword, line); 01619 if (values.size()>0) 01620 { 01621 if(ConvertFromString(values[0], &a_head1stAngle)) 01622 { 01623 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01624 return 1; 01625 } 01626 } 01627 01628 kword = "/gate/"+head_name+"/ring/setAngularPitch"; 01629 values = CheckGATECommand(kword, line); 01630 if (values.size()>0) 01631 { 01632 if(ConvertFromString(values[0], &a_headAngPitch)) 01633 { 01634 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01635 return 1; 01636 } 01637 } 01638 01639 kword = "/gate/"+head_name+"/moves/insert"; 01640 values = CheckGATECommand(kword, line); 01641 if (values.size()>0) 01642 { 01643 if(ConvertFromString(values[0], &head_orbit_name)) 01644 { 01645 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01646 return 1; 01647 } 01648 } 01649 01650 kword = "/gate/"+head_name+"/"+head_orbit_name+"/setSpeed"; 01651 values = CheckGATECommand(kword, line); 01652 if (values.size()>0) 01653 { 01654 if(ConvertFromString(values[0], &head_rot_speed)) 01655 { 01656 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01657 return 1; 01658 } 01659 } 01660 01661 kword = "/gate/"+head_name+"/"+head_orbit_name+"/setPoint2"; 01662 values = CheckGATECommand(kword, line); 01663 if (values.size()>0) 01664 { 01665 int rot_direction; 01666 if(ConvertFromString(values[2], &rot_direction)) 01667 { 01668 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 01669 return 1; 01670 } 01671 01672 a_headRotDirection = rot_direction>0 ? GEO_ROT_CW : GEO_ROT_CCW; 01673 } 01674 01675 // --- Crystals --- 01676 kword = is_head_Y_axis ? 01677 "/gate/"+crystal_name+"/geometry/setXLength": 01678 "/gate/"+crystal_name+"/geometry/setYLength"; 01679 01680 values = CheckGATECommand(kword, line); 01681 if (values.size()>0) 01682 { 01683 if(ConvertFromString(values[0], &a_crystalSizeTrs) ) 01684 { 01685 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 01686 return 1; 01687 } 01688 01689 } 01690 01691 kword = "/gate/"+crystal_name+"/geometry/setZLength"; 01692 values = CheckGATECommand(kword, line); 01693 if (values.size()>0) 01694 { 01695 if(ConvertFromString(values[0], &a_crystalSizeAxl) ) 01696 { 01697 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 01698 return 1; 01699 } 01700 } 01701 01702 01703 // --- Pixels --- 01704 kword = is_head_Y_axis ? 01705 "/gate/"+pixel_name+"/cubicArray/setRepeatNumberX": 01706 "/gate/"+pixel_name+"/cubicArray/setRepeatNumberY"; 01707 01708 01709 values = CheckGATECommand(kword, line); 01710 if (values.size()>0) 01711 { 01712 if(ConvertFromString(values[0], &a_nPixTrs) ) 01713 { 01714 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 01715 return 1; 01716 } 01717 } 01718 01719 kword = "/gate/"+pixel_name+"/cubicArray/setRepeatNumberZ"; 01720 values = CheckGATECommand(kword, line); 01721 if (values.size()>0) 01722 { 01723 if(ConvertFromString(values[0], &a_nPixAxl) ) 01724 { 01725 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 01726 return 1; 01727 } 01728 } 01729 01730 kword = "/gate/application/setTimeStart"; 01731 values = CheckGATECommand(kword, line); 01732 if (values.size()>0) 01733 { 01734 if(ConvertFromString( values[0] , &time_start) ) 01735 { 01736 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01737 return 1; 01738 } 01739 01740 // Convert in ms 01741 if (values.size()>1) 01742 { 01743 if(values[1] == "s") time_start *= 1000; 01744 } 01745 else 01746 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01747 01748 a_start_time_ms = time_start; 01749 } 01750 01751 kword = "/gate/application/setTimeSlice"; 01752 values = CheckGATECommand(kword, line); 01753 if (values.size()>0) 01754 { 01755 if(ConvertFromString( values[0] , &time_slice) ) 01756 { 01757 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01758 return 1; 01759 } 01760 01761 // Convert in ms 01762 if (values.size()>1) 01763 { 01764 if(values[1] == "s") time_slice *= 1000; 01765 } 01766 else 01767 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01768 01769 time_slice_ms = time_slice; 01770 } 01771 01772 kword = "/gate/application/setTimeStop"; 01773 values = CheckGATECommand(kword, line); 01774 if (values.size()>0) 01775 { 01776 if(ConvertFromString( values[0] , &time_stop) ) 01777 { 01778 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01779 return 1; 01780 } 01781 01782 // Convert in ms 01783 if (values.size()>1) 01784 { 01785 if(values[1] == "s") time_stop *= 1000; 01786 } 01787 else 01788 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01789 } 01790 01791 kword = "/gate/application/addSlice"; 01792 values = CheckGATECommand(kword, line); 01793 if (values.size()>0) 01794 { 01795 double time_slice_tmp=0; 01796 01797 if(ConvertFromString( values[0] , &time_slice_tmp) ) 01798 { 01799 Cerr("***** dataConversionUtilities::readMacECAT()-> Error occured while trying to get value from entry '"<< kword << "' in file " << path_mac_files[f] << "!"); 01800 return 1; 01801 } 01802 01803 // Convert in ms 01804 if (values.size()>1) 01805 { 01806 if(values[1] == "s") time_slice_tmp *= 1000; 01807 } 01808 else 01809 Cerr("***** dataConversionUtilities::readMacECAT()-> WARNING : can't read unit of '"<< kword << ". Assuming time in seconds"); 01810 01811 time_slices = time_slice_tmp; 01812 } 01813 } 01814 01815 // Compute duration in ms 01816 a_duration_ms = (time_slices>0) ? 01817 (uint32_t)(time_slices-time_start) : 01818 a_duration_ms; 01819 01820 a_duration_ms = (time_start>=0 && time_stop>=0) ? 01821 (uint32_t)(time_stop-time_start) : 01822 a_duration_ms ; 01823 01824 01825 // Compute the nb of projections and total number of crystals in the system 01826 a_nProjectionsByHead = a_duration_ms / time_slice_ms; 01827 a_nProjectionsTot = a_nHeads*a_nProjectionsByHead; 01828 01829 // Compute angular pitch if not provided in the mac files (=-1) 01830 a_headAngPitch = (a_headAngPitch<0) ? 01831 360./a_nHeads : 01832 a_headAngPitch ; 01833 01834 a_headAngStepDeg = head_rot_speed*time_slice_ms/1000.; 01835 01836 if(head_rot_speed<0 || time_slice_ms<0) 01837 { 01838 Cerr("***** GetGATESystemType() -> Error couldn't find lines '/gate/"+head_name+"/"+head_orbit_name+"/setSpeed' and '/gate/application/setTimeSlice' !" << endl); 01839 Cerr(" These informations are mandatory to compute the projection angle step." << endl); 01840 return 1; 01841 } 01842 01843 mac_file.close(); 01844 } 01845 return 0; 01846 } 01847 01848 01849 01850 01851 // ===================================================================== 01852 // --------------------------------------------------------------------- 01853 // --------------------------------------------------------------------- 01854 // ===================================================================== 01855 /* 01856 \fn ReadIntfSPECT 01857 \param a_pathIntf : path to the interfile header 01858 \param distToDetector : distance between center of rotation and detector surface 01859 \param nHeads : nb of cameras 01860 \param nPixAxl : nb of transaxial pixels 01861 \param nPixTrs : nb of axial pixels 01862 \param crystalSizeAxl : crystal axial dimension 01863 \param crystalSizeTrs : crystal transaxial dimension 01864 \param head1stAngle : head first angle 01865 \param headAngPitch : angular pitch between heads 01866 \param headRotSpeed : angle between projection 01867 \param headRotDirection : rotation direction of the head (CW or CCW) 01868 \param start_time_ms : acquisition start time converted in ms 01869 \param duration_ms : acquisition duration converted in ms 01870 \param vb : verbosity 01871 \brief Recover informations about the scanner element of an ECAT system, and acquisition duration, from a GATE macro file 01872 \return 0 if success, positive value otherwise 01873 */ 01874 int ReadIntfSPECT(string a_pathIntf, 01875 float_t &a_distToDetector, 01876 uint32_t &a_nHeads, 01877 uint32_t &a_nPixAxl, 01878 uint32_t &a_nPixTrs, 01879 float_t &a_crystalSizeAxl, 01880 float_t &a_crystalSizeTrs, 01881 uint32_t &a_nProjectionsTot, 01882 uint32_t &a_nProjectionsByHead, 01883 float_t &a_head1stAngle, 01884 float_t &a_headAngPitchDeg, 01885 float_t &a_headAngStepDeg, 01886 int &a_headRotDirection, 01887 uint32_t &a_start_time_ms, 01888 uint32_t &a_duration_ms, 01889 int vb) 01890 { 01891 // Read all required information in Interfile header 01892 01893 string key; 01894 01895 key = "matrix size [1]"; 01896 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_nPixTrs, 1, KEYWORD_MANDATORY)) 01897 { 01898 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01899 Cerr(" Either key not found or conversion error" << endl); 01900 return 1; 01901 } 01902 01903 key = "matrix size [2]"; 01904 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_nPixAxl, 1, KEYWORD_MANDATORY)) 01905 { 01906 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01907 Cerr(" Either key not found or conversion error" << endl); 01908 return 1; 01909 } 01910 01911 FLTNB size_pix_trs = 1., 01912 size_pix_axl = 1.; 01913 01914 key = "scaling factor (mm/pixel) [1]"; 01915 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &size_pix_trs, 1, KEYWORD_MANDATORY)) 01916 { 01917 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01918 Cerr(" Either key not found or conversion error" << endl); 01919 return 1; 01920 } 01921 01922 key = "scaling factor (mm/pixel) [2]"; 01923 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &size_pix_axl, 1, KEYWORD_MANDATORY)) 01924 { 01925 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01926 Cerr(" Either key not found or conversion error" << endl); 01927 return 1; 01928 } 01929 01930 key = "total number of images"; 01931 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_nProjectionsTot, 1, KEYWORD_MANDATORY)) 01932 { 01933 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01934 Cerr(" Either key not found or conversion error" << endl); 01935 return 1; 01936 } 01937 01938 key = "number of projections"; 01939 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_nProjectionsByHead, 1, KEYWORD_MANDATORY)) 01940 { 01941 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01942 Cerr(" Either key not found or conversion error" << endl); 01943 return 1; 01944 } 01945 01946 key = "number of detector heads"; 01947 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_nHeads, 1, KEYWORD_MANDATORY)) 01948 { 01949 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01950 Cerr(" Either key not found or conversion error" << endl); 01951 return 1; 01952 } 01953 01954 key = "study duration (sec)"; 01955 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_duration_ms, 1, KEYWORD_MANDATORY)) 01956 { 01957 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01958 Cerr(" Either key not found or conversion error" << endl); 01959 return 1; 01960 } 01961 // convert duration in ms 01962 a_duration_ms *= 1000; 01963 01964 01965 FLTNB size_crystal_X =0., 01966 size_crystal_Y =0.; 01967 01968 key = "crystal x dimension (cm)"; 01969 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &size_crystal_X, 1, KEYWORD_OPTIONAL) == 1) 01970 { 01971 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01972 Cerr(" Either key not found or conversion error" << endl); 01973 return 1; 01974 } 01975 01976 key = "crystal y dimension (cm)"; 01977 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &size_crystal_Y, 1, KEYWORD_OPTIONAL) == 1) 01978 { 01979 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01980 Cerr(" Either key not found or conversion error" << endl); 01981 return 1; 01982 } 01983 01984 key = "crystal z dimension (cm)"; 01985 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_crystalSizeAxl, 1, KEYWORD_OPTIONAL) == 1) 01986 { 01987 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 01988 Cerr(" Either key not found or conversion error" << endl); 01989 return 1; 01990 } 01991 01992 // Just pick the larger value to get the transaxial crystal size 01993 // convert in mm (values in cm in Interfile) 01994 if(size_crystal_X>0 && size_crystal_Y>0) 01995 a_crystalSizeTrs = (size_crystal_X>size_crystal_Y) ? size_crystal_X*10. : size_crystal_Y*10.; 01996 a_crystalSizeAxl = (a_crystalSizeAxl>0) ? a_crystalSizeAxl*10. : a_crystalSizeAxl ; 01997 01998 01999 FLTNB head_pos_X =-1., 02000 head_pos_Y =-1., 02001 head_pos_Z =-1.; 02002 02003 key = "head x translation (cm)"; 02004 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &head_pos_X, 1, KEYWORD_OPTIONAL) == 1) 02005 { 02006 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02007 Cerr(" Either key not found or conversion error" << endl); 02008 return 1; 02009 } 02010 02011 key = "head y translation (cm)"; 02012 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &head_pos_Y, 1, KEYWORD_OPTIONAL) == 1) 02013 { 02014 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02015 Cerr(" Either key not found or conversion error" << endl); 02016 return 1; 02017 } 02018 02019 key = "head z translation (cm)"; 02020 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &head_pos_Z, 1, KEYWORD_OPTIONAL) == 1) 02021 { 02022 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02023 Cerr(" Either key not found or conversion error" << endl); 02024 return 1; 02025 } 02026 02027 // Pick the largest value to initialize distance betweeen COR and detector, converted in mm (cm by default) 02028 if(head_pos_X > head_pos_Y) 02029 a_distToDetector = head_pos_X > head_pos_Z ? head_pos_X*10 : head_pos_Z*10; 02030 else 02031 a_distToDetector = head_pos_Y > head_pos_Z ? head_pos_Y*10 : head_pos_Z*10; 02032 02033 02034 key = "direction of rotation"; 02035 string rot_dir; 02036 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &rot_dir, 1, KEYWORD_MANDATORY)) 02037 { 02038 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02039 Cerr(" Either key not found or conversion error" << endl); 02040 return 1; 02041 } 02042 02043 a_headRotDirection = (rot_dir == "CCW") ? GEO_ROT_CCW : GEO_ROT_CW ; 02044 02045 // Could have several keys with this name (each for each head) 02046 // It will recover the value corresponding to the first match 02047 key = "start angle"; 02048 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &a_head1stAngle, 1, KEYWORD_MANDATORY)) 02049 { 02050 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02051 Cerr(" Either key not found or conversion error" << endl); 02052 return 1; 02053 } 02054 02055 key = "extent of rotation"; 02056 uint32_t extent_rotation =360; 02057 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &extent_rotation, 1, KEYWORD_MANDATORY)) 02058 { 02059 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02060 Cerr(" Either key not found or conversion error" << endl); 02061 return 1; 02062 } 02063 02064 a_headAngStepDeg = (FLTNB)extent_rotation / a_nProjectionsByHead; 02065 02066 float_t first_angle = 0; 02067 float_t second_angle = extent_rotation; 02068 02069 key = "start angle"; 02070 if(IntfKeyGetValueFromFile(a_pathIntf.c_str(), key.c_str(), &first_angle, 1, KEYWORD_MANDATORY)) 02071 { 02072 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02073 Cerr(" Either key not found or conversion error" << endl); 02074 return 1; 02075 } 02076 02077 key = "start angle"; 02078 if(IntfKeyGetRecurringValueFromFile(a_pathIntf.c_str(), key.c_str(), &second_angle, 1, KEYWORD_MANDATORY, 1)) 02079 { 02080 Cerr("***** castor-GATERootToCastor :: Error when trying to read key: '"<< key << "' in interfile : " << a_pathIntf << "!" << endl); 02081 Cerr(" Either key not found or conversion error" << endl); 02082 return 1; 02083 } 02084 02085 a_headAngPitchDeg = second_angle - first_angle; 02086 02087 return 0; 02088 } 02089 02090 02091 02092 02093 // ===================================================================== 02094 // --------------------------------------------------------------------- 02095 // --------------------------------------------------------------------- 02096 // ===================================================================== 02097 /* 02098 \fn CreateGeomWithCylindrical() 02099 \param a_pathMac : string containing the path to a GATE macro file 02100 \param a_pathGeom : string containing the path to a CASToR output geom file 02101 \brief Read a GATE macro file containing the description of a cylindricalPET system, and convert it to a geom file 02102 \return 0 if success, positive value otherwise 02103 */ 02104 int CreateGeomWithCylindrical(string a_pathMac, string a_pathGeom) 02105 { 02106 /* Declaring variables */ 02107 02108 // Scanner 02109 string modality; 02110 string scanner_name = GetFileFromPath(a_pathGeom.substr(0,a_pathGeom.find_first_of(".geom"))); 02111 double scanner_radius = -1.; 02112 uint32_t number_of_elements = 0; 02113 string description = "PET system extracted from GATE macro: " + a_pathMac; 02114 02115 // Rsectors 02116 string rsector_name = "rsector"; 02117 uint32_t number_of_rsectors = 1; 02118 double rsector_first_angle = 0.; 02119 double rsector_angular_span = 360.; 02120 double rsector_pos_X = 0.; 02121 double rsector_pos_Y = 0.; 02122 uint32_t rsector_nb_zshifts = 0; 02123 vector <double> vec_rsector_Z_Shift; 02124 bool is_rsector_Y_axis = false; 02125 02126 // Modules 02127 string module_name = "module"; 02128 uint32_t number_of_modules_trans = 1; 02129 uint32_t number_of_modules_axial = 1; 02130 double module_step_trans = 0.; 02131 double module_step_axial = 0.; 02132 double module_size_Y; 02133 double module_size_Z; 02134 02135 // Submodules 02136 string submodule_name = "submodule"; 02137 uint32_t number_of_submodules_trans = 1; 02138 uint32_t number_of_submodules_axial = 1; 02139 double submodule_step_trans = 0.; 02140 double submodule_step_axial = 0.; 02141 double submodule_size_Y; 02142 double submodule_size_Z; 02143 02144 // Crystals 02145 string crystal_name = "crystal"; 02146 uint32_t number_of_crystals_trans = 1; 02147 uint32_t number_of_crystals_axial = 1; 02148 double crystal_step_trans = 0.; 02149 double crystal_step_axial = 0.; 02150 double crystal_size_depth = 0.; 02151 double crystal_size_trans = 0.; 02152 double crystal_size_axial = 0.; 02153 02154 02155 // Layers 02156 int number_of_layers = 0; 02157 string n_layers = "1"; //default value for number_of_layers 02158 vector <uint32_t> number_of_lyr_elts_trans; 02159 vector <uint32_t> number_of_lyr_elts_axial; 02160 02161 vector <string> layers_names; 02162 vector <vector <double> > layers_positions; 02163 vector <double> layers_size_depth; 02164 vector <double> layers_size_trans; 02165 vector <double> layers_size_axial; 02166 02167 vector <double> layers_step_trans; 02168 vector <double> layers_step_axial; 02169 02170 // Optional 02171 uint32_t voxels_number_trans; 02172 uint32_t voxels_number_axial; 02173 double fov_trans; 02174 double fov_axial; 02175 double mean_depth_of_interaction = -1.; 02176 double min_angle_diff = 0.; 02177 02178 02179 02180 02181 02182 // If we have multiple layers, we need to enter multiple values in the .geom. 02183 vector <string> vec_scanner_radius; 02184 02185 // Rsector vectors 02186 vector <string> vec_number_of_rsectors; 02187 vector <string> vec_rsector_first_angle; 02188 02189 02190 // Module vectors 02191 vector <string> vec_number_of_modules_trans; 02192 vector <string> vec_number_of_modules_axial; 02193 vector <string> vec_module_gap_trans; 02194 vector <string> vec_module_gap_axial; 02195 02196 // Submodule vectors 02197 vector <string> vec_number_of_submodules_trans; 02198 vector <string> vec_number_of_submodules_axial; 02199 vector <string> vec_submodule_gap_trans; 02200 vector <string> vec_submodule_gap_axial; 02201 02202 // Crystal vectors 02203 vector <string> vec_number_of_crystals_trans; 02204 vector <string> vec_number_of_crystals_axial; 02205 vector <string> vec_crystal_gap_trans; 02206 vector <string> vec_crystal_gap_axial; 02207 02208 // Optionnal 02209 vector <string> vec_mean_depth_of_interaction; 02210 vector <string> vec_min_angle_diff; 02211 02212 vector<string> path_mac_files; 02213 path_mac_files.push_back(a_pathMac); 02214 02215 // Recover path to all macro files from the main mac file 02216 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 02217 { 02218 Cerr("***** GetGATESystemType() -> Error occured when trying to recover paths to GATE macro files !" << endl); 02219 return 1; 02220 } 02221 02222 // Recover aliases of the different parts of the architecture 02223 if(GetGATEAliasesCylindrical(path_mac_files, rsector_name, module_name, submodule_name, crystal_name, layers_names, 2) ) 02224 { 02225 Cerr("***** GetGATESystemType() -> Error occured when trying to recover aliases for the elements of the cylindricalPET !" << endl); 02226 return 1; 02227 } 02228 02229 // Recover nb of detected layers 02230 n_layers = layers_names.size(); 02231 02232 02233 // Loop to recover all other informations 02234 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 02235 { 02236 ifstream systemMac(path_mac_files[f].c_str(), ios::in); 02237 02238 // Reading the .mac file line by line and update the variables if a matching adress is found 02239 string line; 02240 while(getline(systemMac, line)) 02241 { 02242 vector <string> values; 02243 02244 // Scanner 02245 modality = "PET"; 02246 02247 string entry = ""; 02248 02249 entry = "/gate/"+rsector_name+"/placement/setTranslation"; 02250 values = CheckGATECommand(entry, line); 02251 02252 // Assumes that first rsector located on the X or Y axis 02253 if (values.size()>0) 02254 { 02255 if(ConvertFromString(values[0], &rsector_pos_X) || 02256 ConvertFromString(values[1], &rsector_pos_Y) ) 02257 { 02258 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02259 return 1; 02260 } 02261 02262 // Check first rsector cartesian coordinates is 0 on X or Y axis 02263 // Throw error otherwise 02264 if(rsector_pos_X!=0 && rsector_pos_Y !=0) 02265 { 02266 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02267 Cerr(" Rsector cartesian coordinates on either the X or Y axis expected to be equal to 0 "<< endl); 02268 return 1; 02269 } 02270 02271 // Get the axis on which the first rsector is positionned 02272 if(rsector_pos_Y!=0) is_rsector_Y_axis = true; 02273 02274 scanner_radius = is_rsector_Y_axis ? abs(rsector_pos_Y) : abs(rsector_pos_X) ; 02275 } 02276 02277 // axial FOV computed from rsector z-length 02278 entry = "/gate/"+rsector_name+"/geometry/setZLength"; 02279 values = CheckGATECommand(entry, line); 02280 if (values.size()>0) 02281 { 02282 if(ConvertFromString(values[0], &fov_axial) ) 02283 { 02284 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02285 return 1; 02286 } 02287 02288 // 4mm voxels by default 02289 voxels_number_axial = fov_axial/4 + 1; 02290 } 02291 02292 // --- Rsectors --- 02293 02294 entry = "/gate/"+rsector_name+"/ring/setModuloNumber"; 02295 values = CheckGATECommand(entry, line); 02296 if (values.size()>0) 02297 { 02298 if(ConvertFromString(values[0], &rsector_nb_zshifts) ) 02299 { 02300 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02301 return 1; 02302 } 02303 } 02304 02305 entry = "/gate/"+rsector_name+"/ring/setRepeatNumber"; 02306 values = CheckGATECommand(entry, line); 02307 if (values.size()>0) 02308 { 02309 if(ConvertFromString(values[0], &number_of_rsectors) ) 02310 { 02311 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02312 return 1; 02313 } 02314 } 02315 02316 entry = "/gate/"+rsector_name+"/ring/setFirstAngle"; 02317 values = CheckGATECommand(entry, line); 02318 if (values.size()>0) 02319 { 02320 if(ConvertFromString(values[0], &rsector_first_angle) ) 02321 { 02322 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02323 return 1; 02324 } 02325 } 02326 02327 entry = "/gate/"+rsector_name+"/ring/setAngularSpan"; 02328 values = CheckGATECommand(entry, line); 02329 if (values.size()>0) 02330 { 02331 if(ConvertFromString(values[0], &rsector_angular_span) ) 02332 { 02333 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02334 return 1; 02335 } 02336 } 02337 02338 // Up to 8 zshifts in Gate 02339 for (int i=1; i <8; i++) 02340 { 02341 entry = "/gate/"+rsector_name+"/ring/setZShift"+toString(i); 02342 values = CheckGATECommand(entry, line); 02343 if (values.size()>0) 02344 { 02345 double zshift = 0.; 02346 if(ConvertFromString(values[0], &zshift) ) 02347 { 02348 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02349 return 1; 02350 } 02351 02352 vec_rsector_Z_Shift.push_back(zshift); 02353 } 02354 } 02355 02356 02357 02358 // --- Modules --- 02359 02360 entry = is_rsector_Y_axis ? 02361 "/gate/"+module_name+"/cubicArray/setRepeatNumberX": 02362 "/gate/"+module_name+"/cubicArray/setRepeatNumberY"; 02363 02364 values = CheckGATECommand(entry, line); 02365 if (values.size()>0) 02366 { 02367 if(ConvertFromString(values[0], &number_of_modules_trans) ) 02368 { 02369 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02370 return 1; 02371 } 02372 } 02373 02374 entry = "/gate/"+module_name+"/cubicArray/setRepeatNumberZ"; 02375 values = CheckGATECommand(entry, line); 02376 if (values.size()>0) 02377 { 02378 if(ConvertFromString(values[0], &number_of_modules_axial) ) 02379 { 02380 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02381 return 1; 02382 } 02383 } 02384 02385 entry = is_rsector_Y_axis ? 02386 "/gate/"+module_name+"/geometry/setXLength": 02387 "/gate/"+module_name+"/geometry/setYLength"; 02388 02389 values = CheckGATECommand(entry, line); 02390 if (values.size()>0) 02391 { 02392 if(ConvertFromString(values[0], &module_size_Y) ) 02393 { 02394 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02395 return 1; 02396 } 02397 } 02398 02399 02400 entry = "/gate/"+module_name+"/geometry/setZLength"; 02401 values = CheckGATECommand(entry, line); 02402 if (values.size()>0) 02403 { 02404 if(ConvertFromString(values[0], &module_size_Z) ) 02405 { 02406 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02407 return 1; 02408 } 02409 } 02410 02411 entry = "/gate/"+module_name+"/cubicArray/setRepeatVector"; 02412 values = CheckGATECommand(entry, line); 02413 if (values.size()>0) 02414 { 02415 string trs_step = is_rsector_Y_axis ? 02416 values[0] : 02417 values[1] ; 02418 02419 if(ConvertFromString(trs_step, &module_step_trans) || 02420 ConvertFromString(values[2], &module_step_axial)) 02421 { 02422 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02423 return 1; 02424 } 02425 } 02426 02427 02428 // --- Submodules --- 02429 entry = is_rsector_Y_axis ? 02430 "/gate/"+submodule_name+"/cubicArray/setRepeatNumberX": 02431 "/gate/"+submodule_name+"/cubicArray/setRepeatNumberY"; 02432 02433 02434 values = CheckGATECommand(entry, line); 02435 if (values.size()>0) 02436 { 02437 if(ConvertFromString(values[0], &number_of_submodules_trans) ) 02438 { 02439 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02440 return 1; 02441 } 02442 } 02443 02444 entry = "/gate/"+submodule_name+"/cubicArray/setRepeatNumberZ"; 02445 values = CheckGATECommand(entry, line); 02446 if (values.size()>0) 02447 { 02448 if(ConvertFromString(values[0], &number_of_submodules_axial) ) 02449 { 02450 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02451 return 1; 02452 } 02453 } 02454 02455 entry = is_rsector_Y_axis ? 02456 "/gate/"+submodule_name+"/geometry/setXLength": 02457 "/gate/"+submodule_name+"/geometry/setYLength"; 02458 02459 values = CheckGATECommand(entry, line); 02460 if (values.size()>0) 02461 { 02462 if(ConvertFromString(values[0], &submodule_size_Y) ) 02463 { 02464 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02465 return 1; 02466 } 02467 } 02468 02469 entry = "/gate/"+submodule_name+"/geometry/setZLength"; 02470 values = CheckGATECommand(entry, line); 02471 if (values.size()>0) 02472 { 02473 if(ConvertFromString(values[0], &submodule_size_Z) ) 02474 { 02475 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02476 return 1; 02477 } 02478 } 02479 02480 entry = "/gate/"+submodule_name+"/cubicArray/setRepeatVector"; 02481 values = CheckGATECommand(entry, line); 02482 if (values.size()>0) 02483 { 02484 string trs_step = is_rsector_Y_axis ? 02485 values[0] : 02486 values[1] ; 02487 02488 if(ConvertFromString(trs_step, &submodule_step_trans) || 02489 ConvertFromString(values[2], &submodule_step_axial) ) 02490 { 02491 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02492 return 1; 02493 } 02494 } 02495 02496 02497 02498 02499 // --- Crystals --- 02500 entry = is_rsector_Y_axis ? 02501 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberX": 02502 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberY"; 02503 02504 02505 values = CheckGATECommand(entry, line); 02506 if (values.size()>0) 02507 { 02508 if(ConvertFromString(values[0], &number_of_crystals_trans) ) 02509 { 02510 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02511 return 1; 02512 } 02513 } 02514 02515 entry = "/gate/"+crystal_name+"/cubicArray/setRepeatNumberZ"; 02516 values = CheckGATECommand(entry, line); 02517 if (values.size()>0) 02518 { 02519 if(ConvertFromString(values[0], &number_of_crystals_axial) ) 02520 { 02521 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02522 return 1; 02523 } 02524 } 02525 02526 entry = "/gate/"+crystal_name+"/geometry/setXLength"; 02527 values = CheckGATECommand(entry, line); 02528 if (values.size()>0) 02529 { 02530 double x_length; 02531 if(ConvertFromString(values[0], &x_length) ) 02532 { 02533 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02534 return 1; 02535 } 02536 02537 if (is_rsector_Y_axis) 02538 crystal_size_trans = x_length; 02539 else 02540 crystal_size_depth = x_length; 02541 } 02542 02543 entry = "/gate/"+crystal_name+"/geometry/setYLength"; 02544 values = CheckGATECommand(entry, line); 02545 if (values.size()>0) 02546 { 02547 double y_length; 02548 if(ConvertFromString(values[0], &y_length) ) 02549 { 02550 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02551 return 1; 02552 } 02553 02554 if (is_rsector_Y_axis) 02555 crystal_size_depth = y_length; 02556 else 02557 crystal_size_trans = y_length; 02558 02559 } 02560 02561 entry = "/gate/"+crystal_name+"/geometry/setZLength"; 02562 values = CheckGATECommand(entry, line); 02563 if (values.size()>0) 02564 { 02565 if(ConvertFromString(values[0], &crystal_size_axial) ) 02566 { 02567 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02568 return 1; 02569 } 02570 } 02571 02572 entry = "/gate/"+crystal_name+"/cubicArray/setRepeatVector"; 02573 values = CheckGATECommand(entry, line); 02574 if (values.size()>0) 02575 { 02576 string trs_step = is_rsector_Y_axis ? 02577 values[0] : 02578 values[1] ; 02579 02580 if(ConvertFromString(trs_step, &crystal_step_trans) || 02581 ConvertFromString(values[2], &crystal_step_axial) ) 02582 { 02583 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02584 return 1; 02585 } 02586 } 02587 02588 02589 02590 02591 // --- Layers --- 02592 entry = "/gate/"+crystal_name+"/daughters/name"; 02593 values = CheckGATECommand(entry, line); 02594 if (values.size()>0) 02595 { 02596 layers_names.push_back(values[0]); 02597 number_of_layers++; 02598 } 02599 02600 // loop on the layers to get their specific information as they are integrated in layers_names. 02601 for (int i=0; i < number_of_layers; i++) 02602 { 02603 // assumes that first block is positionned on the x-axis 02604 entry = "/gate/"+layers_names[i]+"/placement/setTranslation"; 02605 values = CheckGATECommand(entry, line); 02606 02607 if (values.size()>0) 02608 { 02609 vector<double> layer_pos; 02610 for(int d=0 ; d<3 ; d++) 02611 { 02612 double pos; 02613 if(ConvertFromString(values[d], &pos) ) 02614 { 02615 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02616 return 1; 02617 } 02618 layer_pos.push_back(pos); 02619 } 02620 02621 if(is_rsector_Y_axis) 02622 { 02623 double pos = layer_pos[1]; 02624 layer_pos[1] = layer_pos[0]; 02625 layer_pos[0] = pos; 02626 } 02627 02628 layers_positions.push_back(layer_pos); 02629 } 02630 02631 02632 // assumes that first block is positionned on the x-axis 02633 entry = "/gate/"+layers_names[i]+"/geometry/setXLength"; 02634 values = CheckGATECommand(entry, line); 02635 02636 if (values.size()>0) 02637 { 02638 double xlength; 02639 if(ConvertFromString(values[0], &xlength) ) 02640 { 02641 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02642 return 1; 02643 } 02644 02645 if (is_rsector_Y_axis) 02646 layers_size_trans.push_back(xlength); 02647 else 02648 layers_size_depth.push_back(xlength); 02649 } 02650 02651 // assumes that first block is positionned on the x-axis 02652 entry = "/gate/"+layers_names[i]+"/geometry/setYLength"; 02653 values = CheckGATECommand(entry, line); 02654 if (values.size()>0) 02655 { 02656 double ylength; 02657 if(ConvertFromString(values[0], &ylength) ) 02658 { 02659 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02660 return 1; 02661 } 02662 if (is_rsector_Y_axis) 02663 layers_size_depth.push_back(ylength); 02664 else 02665 layers_size_trans.push_back(ylength); 02666 } 02667 02668 entry = "/gate/"+layers_names[i]+"/geometry/setZLength"; 02669 values = CheckGATECommand(entry, line); 02670 if (values.size()>0) 02671 { 02672 double zlength; 02673 if(ConvertFromString(values[0], &zlength) ) 02674 { 02675 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02676 return 1; 02677 } 02678 layers_size_axial.push_back(zlength); 02679 } 02680 02681 // Check if a repeater if applied to the layers elements 02682 // In this case, the total number of crystals for a layer[i] 02683 // will be crystal elts*layers elts[i] 02684 entry = is_rsector_Y_axis ? 02685 "/gate/"+layers_names[i]+"/cubicArray/setRepeatNumberX": 02686 "/gate/"+layers_names[i]+"/cubicArray/setRepeatNumberY"; 02687 02688 values = CheckGATECommand(entry, line); 02689 if (values.size()>0) 02690 { 02691 uint32_t step_trs; 02692 if(ConvertFromString(values[0], &step_trs)) 02693 { 02694 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02695 return 1; 02696 } 02697 number_of_lyr_elts_trans.push_back(step_trs); 02698 } 02699 02700 entry = "/gate/"+layers_names[i]+"/cubicArray/setRepeatNumberZ"; 02701 values = CheckGATECommand(entry, line); 02702 if (values.size()>0) 02703 { 02704 uint32_t step_axl; 02705 if(ConvertFromString(values[0], &step_axl)) 02706 { 02707 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02708 return 1; 02709 } 02710 number_of_lyr_elts_axial.push_back(step_axl); 02711 } 02712 02713 02714 entry = "/gate/"+layers_names[i]+"/cubicArray/setRepeatVector"; 02715 values = CheckGATECommand(entry, line); 02716 if (values.size()>0) 02717 { 02718 string trs_step = is_rsector_Y_axis ? 02719 values[0] : 02720 values[1] ; 02721 02722 double step_trs, step_axl; 02723 if(ConvertFromString(trs_step, &step_trs) || 02724 ConvertFromString(values[2], &step_axl) ) 02725 { 02726 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02727 return 1; 02728 } 02729 02730 02731 layers_step_trans.push_back(step_trs); 02732 layers_step_axial.push_back(step_axl); 02733 } 02734 02735 } 02736 02737 02738 entry = "/gate/digitizer/Coincidences/minSectorDifference"; 02739 values = CheckGATECommand(entry, line); 02740 if (values.size()>0) 02741 { 02742 FLTNB min_sector_diff= 0.; 02743 02744 if(ConvertFromString(values[0], &min_sector_diff)) 02745 { 02746 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 02747 return 1; 02748 } 02749 02750 min_angle_diff = 360./number_of_rsectors*min_sector_diff; 02751 } 02752 02753 } 02754 systemMac.close(); 02755 02756 } 02757 02758 // Compute total number of crystal elements 02759 if(number_of_layers == 0) 02760 { 02761 number_of_elements = number_of_rsectors 02762 * number_of_modules_trans 02763 * number_of_modules_axial 02764 * number_of_submodules_trans 02765 * number_of_submodules_axial 02766 * number_of_crystals_trans 02767 * number_of_crystals_axial; 02768 } 02769 else 02770 for(int l=0 ; l<number_of_layers ; l++) 02771 { 02772 int32_t nb_crystals_layer = number_of_rsectors 02773 * number_of_modules_trans 02774 * number_of_modules_axial 02775 * number_of_submodules_trans 02776 * number_of_submodules_axial 02777 * number_of_crystals_trans 02778 * number_of_crystals_axial; 02779 02780 // Add layer elements if repeaters have been used on layers 02781 if(number_of_lyr_elts_trans.size()>0 || number_of_lyr_elts_axial.size()>0 ) 02782 nb_crystals_layer *= number_of_lyr_elts_trans[l] * number_of_lyr_elts_axial[l]; 02783 02784 number_of_elements += nb_crystals_layer; 02785 } 02786 02787 02788 // Transaxial FOV defined as scanner radius / 2 02789 fov_trans = scanner_radius/2; 02790 // 4mm voxels by default 02791 voxels_number_trans = fov_trans/2 + 1; 02792 02793 02794 // Compute gaps 02795 if (module_step_axial - module_size_Z >= 0) 02796 module_step_axial = module_step_axial - module_size_Z; 02797 02798 if (module_step_trans - module_size_Y >= 0) 02799 module_step_trans =module_step_trans - module_size_Y; 02800 02801 if (submodule_step_axial - submodule_size_Z >= 0) 02802 submodule_step_axial = submodule_step_axial - submodule_size_Z; 02803 02804 if (submodule_step_trans - submodule_size_Y >= 0) 02805 submodule_step_trans = submodule_step_trans - submodule_size_Y; 02806 02807 if (crystal_step_axial - crystal_size_axial >= 0) 02808 crystal_step_axial = crystal_step_axial - crystal_size_axial; 02809 02810 if (crystal_step_trans - crystal_size_trans >= 0) 02811 crystal_step_trans = crystal_step_trans - crystal_size_trans; 02812 02813 02814 // Compute first angle 02815 // CASToR x and y axis for rotation inverted in comparison with GATE 02816 rsector_first_angle = round(atan2f(rsector_pos_X , rsector_pos_Y) * 180. / M_PI) 02817 - rsector_first_angle; 02818 02819 02820 // Convert angular span to the CASToR convention 02821 if(rsector_angular_span != 360.) 02822 rsector_angular_span *= (number_of_rsectors+1)/number_of_rsectors; 02823 02824 // Computing the other parts of variables, if their is more than one layer 02825 if (number_of_layers > 0) 02826 { 02827 for (int l=0; l < number_of_layers ; l++) 02828 { 02829 // Scanner vectors 02830 vec_scanner_radius.push_back( toString(scanner_radius+layers_positions[l][0]) ); 02831 02832 // Rsector vectors 02833 vec_number_of_rsectors.push_back( toString(number_of_rsectors) ); 02834 vec_rsector_first_angle.push_back( toString(rsector_first_angle) ); 02835 02836 // Module vectors 02837 vec_number_of_modules_trans.push_back( toString(number_of_modules_trans) ); 02838 vec_number_of_modules_axial.push_back( toString(number_of_modules_axial) ); 02839 vec_module_gap_trans.push_back( toString(module_step_trans) ); 02840 vec_module_gap_axial.push_back( toString(module_step_axial) ); 02841 02842 // Submodule vectors 02843 vec_number_of_submodules_trans.push_back( toString(number_of_submodules_trans) ); 02844 vec_number_of_submodules_axial.push_back( toString(number_of_submodules_axial) ); 02845 vec_submodule_gap_trans.push_back( toString(submodule_step_trans) ); 02846 vec_submodule_gap_axial.push_back( toString(submodule_step_axial) ); 02847 02848 // Crystal vectors 02849 uint32_t nb_tot_trs_cry = (number_of_lyr_elts_trans.size()>0) ? 02850 number_of_lyr_elts_trans[l]*number_of_crystals_trans : 02851 number_of_crystals_trans ; 02852 02853 uint32_t nb_tot_axl_cry = (number_of_lyr_elts_axial.size()>0) ? 02854 number_of_lyr_elts_axial[l]*number_of_crystals_axial : 02855 number_of_crystals_axial ; 02856 02857 vec_number_of_crystals_trans.push_back( toString(nb_tot_trs_cry) ); 02858 vec_number_of_crystals_axial.push_back( toString(nb_tot_axl_cry) ); 02859 02860 02861 02862 if(layers_step_trans.size()>0 && 02863 layers_step_trans.size()>0) 02864 { 02865 if (layers_step_trans[l] - layers_size_trans[l] >= 0) 02866 crystal_step_trans = layers_step_trans[l] - layers_size_trans[l]; 02867 02868 if (layers_step_axial[l] - layers_size_axial[l] >= 0) 02869 crystal_step_axial = layers_step_axial[l] - layers_size_axial[l]; 02870 } 02871 else 02872 { 02873 crystal_step_trans = crystal_step_trans ; 02874 crystal_step_axial = crystal_step_axial ; 02875 } 02876 02877 vec_crystal_gap_trans.push_back( toString(crystal_step_trans) ); 02878 vec_crystal_gap_axial.push_back( toString(crystal_step_axial) ); 02879 02880 // Optionnal 02881 vec_mean_depth_of_interaction.push_back( toString(mean_depth_of_interaction) ); 02882 vec_min_angle_diff.push_back( toString(min_angle_diff) ); 02883 } 02884 } 02885 else 02886 { 02887 // Layer dimensions = crystal dimensions 02888 layers_size_depth.push_back( crystal_size_depth ); 02889 layers_size_trans.push_back( crystal_size_trans ); 02890 layers_size_axial.push_back( crystal_size_axial ); 02891 02892 02893 // Scanner vectors 02894 vec_scanner_radius.push_back( toString(scanner_radius) ); 02895 02896 // Rsector vectors 02897 vec_number_of_rsectors.push_back( toString(number_of_rsectors) ); 02898 vec_rsector_first_angle.push_back( toString(rsector_first_angle) ); 02899 02900 // Module vectors 02901 vec_number_of_modules_trans.push_back( toString(number_of_modules_trans) ); 02902 vec_number_of_modules_axial.push_back( toString(number_of_modules_axial) ); 02903 vec_module_gap_trans.push_back( toString(module_step_trans) ); 02904 vec_module_gap_axial.push_back( toString(module_step_axial) ); 02905 02906 // Submodule vectors 02907 vec_number_of_submodules_trans.push_back( toString(number_of_submodules_trans) ); 02908 vec_number_of_submodules_axial.push_back( toString(number_of_submodules_axial) ); 02909 vec_submodule_gap_trans.push_back( toString(submodule_step_trans) ); 02910 vec_submodule_gap_axial.push_back( toString(submodule_step_axial) ); 02911 02912 // Crystal vectors 02913 vec_number_of_crystals_trans.push_back( toString(number_of_crystals_trans) ); 02914 vec_number_of_crystals_axial.push_back( toString(number_of_crystals_axial) ); 02915 vec_crystal_gap_trans.push_back( toString(crystal_step_trans) ); 02916 vec_crystal_gap_axial.push_back( toString(crystal_step_axial) ); 02917 02918 // Optionnal 02919 vec_mean_depth_of_interaction.push_back( toString(mean_depth_of_interaction) ); 02920 vec_min_angle_diff.push_back( toString(min_angle_diff) ); 02921 02922 // Update the real number of layers 02923 number_of_layers = 1; 02924 } 02925 02926 02927 // Write the .geom file 02928 ofstream fileGeom(a_pathGeom.c_str(), ios::out | ios::trunc); 02929 if(fileGeom) 02930 { 02931 fileGeom <<"# comments" << endl; 02932 fileGeom <<"# Y _________ "<< endl; 02933 fileGeom <<"# | / _ \\ \\ "<< endl; 02934 fileGeom <<"# | | / \\ | |"<< endl; 02935 fileGeom <<"# Z ____| | | | | |"<< endl; 02936 fileGeom <<"# \\ | | | | |" << endl; 02937 fileGeom <<"# \\ | \\_/ | |" << endl; 02938 fileGeom <<"# X \\___/_____/" << endl; 02939 fileGeom <<"# positions in millimeters"<< endl; 02940 fileGeom <<"# scanner axis is z" << endl; 02941 fileGeom <<"# Use comma without space as separator in the tables." << endl; 02942 02943 fileGeom << ""<< endl; 02944 02945 // Mandatory fields 02946 fileGeom << "# MANDATORY FIELDS"<< endl; 02947 fileGeom << "modality : " << modality << endl; 02948 fileGeom << "scanner name : " << scanner_name << endl; 02949 fileGeom << "number of elements : " << number_of_elements << endl; 02950 fileGeom << "number of layers : " << number_of_layers << endl; 02951 fileGeom << "" << endl; 02952 fileGeom << "voxels number transaxial : " << voxels_number_trans << endl; 02953 fileGeom << "voxels number axial : " << voxels_number_axial << endl; 02954 02955 fileGeom << "field of view transaxial : " << fov_trans << endl; 02956 fileGeom << "field of view axial : " << fov_axial << endl << endl; 02957 fileGeom << "description : " << description << endl; 02958 fileGeom << "" << endl; 02959 WriteVector(fileGeom,"scanner radius : ",vec_scanner_radius); 02960 WriteVector(fileGeom,"number of rsectors : ",vec_number_of_rsectors); 02961 WriteVector(fileGeom,"number of crystals transaxial : ",vec_number_of_crystals_trans); 02962 WriteVector(fileGeom, "number of crystals axial : ",vec_number_of_crystals_axial); 02963 fileGeom << ""<< endl; 02964 WriteVector(fileGeom, "layers size depth : ", layers_size_depth); 02965 WriteVector(fileGeom, "layers size transaxial : ", layers_size_trans); 02966 WriteVector(fileGeom, "layers size axial : ", layers_size_axial); 02967 fileGeom << ""<< endl; 02968 fileGeom << ""<< endl; 02969 02970 // Optional fields 02971 fileGeom << "# OPTIONAL FIELDS"<< endl; 02972 WriteVector(fileGeom,"rsectors first angle : ",vec_rsector_first_angle); 02973 WriteVector(fileGeom,"number of modules transaxial : ",vec_number_of_modules_trans); 02974 WriteVector(fileGeom, "number of modules axial : ",vec_number_of_modules_axial); 02975 WriteVector(fileGeom,"module gap transaxial : ",vec_module_gap_trans); 02976 WriteVector(fileGeom,"module gap axial : ",vec_module_gap_axial); 02977 WriteVector(fileGeom,"number of submodules transaxial : ",vec_number_of_submodules_trans); 02978 WriteVector(fileGeom, "number of submodules axial : ",vec_number_of_submodules_axial); 02979 WriteVector(fileGeom,"submodule gap transaxial : ",vec_submodule_gap_trans); 02980 WriteVector(fileGeom,"submodule gap axial : ",vec_submodule_gap_axial); 02981 WriteVector(fileGeom,"crystal gap transaxial : ",vec_crystal_gap_trans); 02982 WriteVector(fileGeom,"crystal gap axial : ",vec_crystal_gap_axial); 02983 WriteVector(fileGeom, "mean depth of interaction : ", vec_mean_depth_of_interaction); 02984 fileGeom << ""<< endl; 02985 02986 // Write only if different from default values 02987 if(min_angle_diff > 0.) fileGeom << "min angle difference : " << min_angle_diff << endl; 02988 if(rsector_angular_span != 360.) fileGeom << "rsectors angular span : " << rsector_angular_span << endl; 02989 if(rsector_nb_zshifts > 0) fileGeom << "rsectors nbZShift :" << rsector_nb_zshifts << endl; 02990 if(!vec_rsector_Z_Shift.empty()) WriteVector(fileGeom, "rsectors ZShift : ", vec_rsector_Z_Shift); 02991 02992 fileGeom.close(); 02993 02994 cout << "Output geom file written at :" << a_pathGeom << endl; 02995 } 02996 else 02997 { 02998 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Couldn't open geom file for writing "<< a_pathGeom << " !" << endl); 02999 return 1; 03000 } 03001 03002 03003 03004 return 0; 03005 } 03006 03007 03008 03009 // ===================================================================== 03010 // --------------------------------------------------------------------- 03011 // --------------------------------------------------------------------- 03012 // ===================================================================== 03013 /* 03014 \fn CreateGeomWithECAT() 03015 \param a_pathMac : string containing the path to a GATE macro file 03016 \param a_pathGeom : string containing the path to a CASToR output geom file 03017 \brief Read a GATE macro file containing the description of an ecat system, and convert it to a geom file 03018 \return 0 if success, positive value otherwise 03019 */ 03020 int CreateGeomWithECAT(string a_pathMac, string a_pathGeom) 03021 { 03022 /* Declaring variables */ 03023 03024 // Scanner 03025 string modality; 03026 string scanner_name = GetFileFromPath(a_pathGeom.substr(0,a_pathGeom.find_first_of(".geom"))); 03027 double scanner_radius = -1.; 03028 uint32_t number_of_elements = 0; 03029 string description = "ECAT system extracted from GATE macro: " + a_pathMac; 03030 03031 // blocks 03032 string block_name = "block"; 03033 uint32_t number_of_blocks = 1; 03034 uint32_t number_of_blocks_trans = 1; 03035 uint32_t number_of_blocks_axial = 1; 03036 double block_step_trans = 0.; 03037 double block_step_axial = 0.; 03038 double block_size_Y; 03039 double block_size_Z; 03040 double block_pos_X = 0.; 03041 double block_pos_Y = 0.; 03042 double block_first_angle = 0.; 03043 double block_angular_span = 360.; 03044 uint32_t block_nb_zshifts = 0; 03045 vector <double> vec_block_Z_Shift; 03046 bool is_block_Y_axis = false; 03047 03048 03049 03050 // Crystals 03051 string crystal_name = "crystal"; 03052 uint32_t number_of_crystals_trans = 1; 03053 uint32_t number_of_crystals_axial = 1; 03054 double crystal_step_trans = 0.; 03055 double crystal_step_axial = 0.; 03056 double crystal_size_depth = 0.; 03057 double crystal_size_trans = 0.; 03058 double crystal_size_axial = 0.; 03059 03060 // Optional 03061 uint32_t voxels_number_trans; 03062 uint32_t voxels_number_axial; 03063 double fov_trans; 03064 double fov_axial; 03065 double min_angle_diff = 0.; 03066 03067 vector<string> path_mac_files; 03068 path_mac_files.push_back(a_pathMac); 03069 03070 // Recover path to all macro files from the main mac file 03071 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 03072 { 03073 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Error occured when trying to recover paths to GATE macro files !" << endl); 03074 return 1; 03075 } 03076 03077 03078 // Recover aliases of the different parts of the architecture 03079 if(GetGATEAliasesEcat(path_mac_files, block_name, crystal_name, 2) ) 03080 { 03081 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Error occured when trying to recover aliases for the elements of the ecat !" << endl); 03082 return 1; 03083 } 03084 03085 03086 // Loop to recover all other informations 03087 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 03088 { 03089 ifstream systemMac(path_mac_files[f].c_str(), ios::in); 03090 03091 string line; 03092 while(getline(systemMac, line)) 03093 { 03094 // Scanner 03095 modality = "PET"; 03096 vector <string> values; 03097 string entry = ""; 03098 03099 03100 entry = "/gate/"+block_name+"/placement/setTranslation"; 03101 values = CheckGATECommand(entry, line); 03102 if (values.size()>0) 03103 { 03104 if(ConvertFromString(values[0], &block_pos_X) || 03105 ConvertFromString(values[1], &block_pos_Y) ) 03106 { 03107 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03108 return 1; 03109 } 03110 03111 // Check first block cartesian coordinates is 0 on X or Y axis 03112 // Throw error otherwise 03113 if(block_pos_X!=0 && block_pos_Y !=0) 03114 { 03115 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03116 Cerr(" Block cartesian coordinates on either the X or Y axis expected to be equal to 0 "<< endl); 03117 return 1; 03118 } 03119 03120 // Get the axis on which the first rsector is positionned 03121 if(block_pos_Y!=0) is_block_Y_axis = true; 03122 03123 scanner_radius = is_block_Y_axis ? abs(block_pos_Y) : abs(block_pos_X) ; 03124 } 03125 03126 entry = "/gate/"+block_name+"/geometry/setZLength"; 03127 values = CheckGATECommand(entry, line); 03128 if (values.size()>0) 03129 { 03130 if(ConvertFromString(values[0], &fov_axial)) 03131 { 03132 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03133 return 1; 03134 } 03135 // 4mm voxels by default 03136 voxels_number_axial = fov_axial/4 + 1; 03137 } 03138 03139 03140 // blocks 03141 entry = "/gate/"+block_name+"/ring/setRepeatNumber"; 03142 values = CheckGATECommand(entry, line); 03143 if (values.size()>0) 03144 { 03145 if(ConvertFromString(values[0], &number_of_blocks)) 03146 { 03147 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03148 return 1; 03149 } 03150 } 03151 03152 03153 entry = "/gate/"+block_name+"/ring/setFirstAngle"; 03154 values = CheckGATECommand(entry, line); 03155 if (values.size()>0) 03156 { 03157 if(ConvertFromString(values[0], &block_first_angle)) 03158 { 03159 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03160 return 1; 03161 } 03162 } 03163 03164 03165 entry = "/gate/"+block_name+"/ring/setAngularSpan"; 03166 values = CheckGATECommand(entry, line); 03167 if (values.size()>0) 03168 { 03169 if(ConvertFromString(values[0], &block_angular_span)) 03170 { 03171 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03172 return 1; 03173 } 03174 } 03175 03176 for (int i=1; i<8; i++) 03177 { 03178 entry = "/gate/"+block_name+"/ring/setZShift"+toString(i); 03179 values = CheckGATECommand(entry, line); 03180 if (values.size()>0) 03181 { 03182 double zshift; 03183 if(ConvertFromString(values[0], &zshift)) 03184 { 03185 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03186 return 1; 03187 } 03188 vec_block_Z_Shift.push_back(zshift); 03189 } 03190 } 03191 03192 03193 03194 03195 // Modules 03196 entry = "/gate/"+block_name+"/linear/setRepeatNumber"; 03197 values = CheckGATECommand(entry, line); 03198 if (values.size()>0) 03199 { 03200 if(ConvertFromString(values[0], &number_of_blocks_axial)) 03201 { 03202 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03203 return 1; 03204 } 03205 } 03206 03207 entry = is_block_Y_axis ? 03208 "/gate/"+block_name+"/geometry/setXLength": 03209 "/gate/"+block_name+"/geometry/setYLength"; 03210 03211 values = CheckGATECommand(entry, line); 03212 if (values.size()>0) 03213 { 03214 if(ConvertFromString(values[0], &block_size_Y)) 03215 { 03216 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03217 return 1; 03218 } 03219 } 03220 03221 entry = "/gate/"+block_name+"/geometry/setZLength"; 03222 values = CheckGATECommand(entry, line); 03223 if (values.size()>0) 03224 { 03225 if(ConvertFromString(values[0], &block_size_Z)) 03226 { 03227 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03228 return 1; 03229 } 03230 } 03231 03232 entry = "/gate/"+block_name+"/linear/setRepeatVector"; 03233 values = CheckGATECommand(entry, line); 03234 if (values.size()>0) 03235 { 03236 if(ConvertFromString(values[2], &block_step_axial)) 03237 { 03238 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03239 return 1; 03240 } 03241 } 03242 03243 03244 03245 // Crystals 03246 03247 // assumes that first block is positionned on the x-axis 03248 entry = is_block_Y_axis ? 03249 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberX": 03250 "/gate/"+crystal_name+"/cubicArray/setRepeatNumberY"; 03251 03252 values = CheckGATECommand(entry, line); 03253 if (values.size()>0) 03254 { 03255 if(ConvertFromString(values[0], &number_of_crystals_trans)) 03256 { 03257 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03258 return 1; 03259 } 03260 } 03261 03262 entry = "/gate/"+crystal_name+"/cubicArray/setRepeatNumberZ"; 03263 values = CheckGATECommand(entry, line); 03264 if (values.size()>0) 03265 { 03266 if(ConvertFromString(values[0], &number_of_crystals_axial)) 03267 { 03268 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03269 return 1; 03270 } 03271 } 03272 03273 03274 entry = "/gate/"+crystal_name+"/cubicArray/setRepeatVector"; 03275 values = CheckGATECommand(entry, line); 03276 if (values.size()>0) 03277 { 03278 string trs_step = is_block_Y_axis ? 03279 values[0] : 03280 values[1] ; 03281 03282 if(ConvertFromString(values[1], &crystal_step_trans) || 03283 ConvertFromString(values[2], &crystal_step_axial) ) 03284 { 03285 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03286 return 1; 03287 } 03288 } 03289 03290 03291 // assumes that first block is positionned on the x-axis 03292 entry = "/gate/"+crystal_name+"/geometry/setXLength"; 03293 values = CheckGATECommand(entry, line); 03294 if (values.size()>0) 03295 { 03296 double x_length; 03297 if(ConvertFromString(values[0], &x_length)) 03298 { 03299 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03300 return 1; 03301 } 03302 03303 if (is_block_Y_axis) 03304 crystal_size_trans = x_length; 03305 else 03306 crystal_size_depth = x_length; 03307 } 03308 03309 03310 // assumes that first block is positionned on the x-axis 03311 entry = "/gate/"+crystal_name+"/geometry/setYLength"; 03312 values = CheckGATECommand(entry, line); 03313 if (values.size()>0) 03314 { 03315 double y_length; 03316 if(ConvertFromString(values[0], &y_length)) 03317 { 03318 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03319 return 1; 03320 } 03321 03322 if (is_block_Y_axis) 03323 crystal_size_depth = y_length; 03324 else 03325 crystal_size_trans = y_length; 03326 } 03327 03328 entry = "/gate/"+crystal_name+"/geometry/setZLength"; 03329 values = CheckGATECommand(entry, line); 03330 if (values.size()>0) 03331 { 03332 if(ConvertFromString(values[0], &crystal_size_axial)) 03333 { 03334 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03335 return 1; 03336 } 03337 } 03338 03339 03340 03341 03342 entry = "/gate/digitizer/Coincidences/minSectorDifference"; 03343 values = CheckGATECommand(entry, line); 03344 if (values.size()>0) 03345 { 03346 double min_sector_diff; 03347 if(ConvertFromString(values[0], &min_sector_diff)) 03348 { 03349 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Conversion error occured while trying to parse line: " << line<< endl); 03350 return 1; 03351 } 03352 min_angle_diff = 360/number_of_blocks*min_sector_diff; 03353 } 03354 03355 } 03356 systemMac.close(); 03357 } 03358 03359 number_of_elements = number_of_blocks * 03360 number_of_blocks_axial * 03361 number_of_crystals_trans * 03362 number_of_crystals_axial; 03363 03364 // Transaxial FOV defined as scanner radius / 2 03365 fov_trans = scanner_radius/2; 03366 // 4mm voxels by default 03367 voxels_number_trans = fov_trans/2 + 1; 03368 03369 if (crystal_step_axial - crystal_size_axial >= 0) 03370 crystal_step_axial = crystal_step_axial - crystal_size_axial; 03371 03372 if (crystal_step_trans - crystal_size_trans >= 0) 03373 crystal_step_trans = crystal_step_trans - crystal_size_trans; 03374 03375 if (block_step_axial - block_size_Z >= 0) 03376 block_step_axial = block_step_axial - block_size_Z; 03377 03378 if (block_step_trans - block_size_Y >= 0) 03379 block_step_trans = block_step_trans - block_size_Y; 03380 03381 // CASToR x and y axis for rotation inverted in comparison with GATE 03382 block_first_angle = round(atan2f(block_pos_X , block_pos_Y) * 180. / M_PI) 03383 - block_first_angle; 03384 03385 03386 // Writing the .geom file 03387 ofstream fileGeom(a_pathGeom.c_str(), ios::out | ios::trunc); 03388 03389 if(fileGeom) 03390 { 03391 fileGeom <<"# comments" << endl; 03392 fileGeom <<"# Y _________ "<< endl; 03393 fileGeom <<"# | / _ \\ \\ "<< endl; 03394 fileGeom <<"# | | / \\ | |"<< endl; 03395 fileGeom <<"# Z ____| | | | | |"<< endl; 03396 fileGeom <<"# \\ | | | | |" << endl; 03397 fileGeom <<"# \\ | \\_/ | |" << endl; 03398 fileGeom <<"# X \\___/_____/" << endl; 03399 fileGeom <<"# positions in millimeters"<< endl; 03400 fileGeom <<"# scanner axis is z" << endl; 03401 fileGeom <<"# Use comma without space as separator in the tables." << endl; 03402 03403 03404 fileGeom << "modality : " << modality << endl; 03405 fileGeom << "scanner name : " << scanner_name << endl; 03406 fileGeom << "number of elements : " << number_of_elements << endl; 03407 fileGeom << "number of layers : " << "1" << endl; 03408 fileGeom << "" << endl; 03409 fileGeom << "# default reconstruction parameters" << endl; 03410 fileGeom << "voxels number transaxial : " << voxels_number_trans << endl; 03411 fileGeom << "voxels number axial : " << voxels_number_axial << endl; 03412 03413 fileGeom << "field of view transaxial : " << fov_trans << endl; 03414 fileGeom << "field of view axial : " << fov_axial << endl; 03415 fileGeom << "" << endl; 03416 fileGeom << "description : " << description << endl; 03417 fileGeom << "" << endl; 03418 fileGeom << "scanner radius : " << scanner_radius << endl; 03419 fileGeom << "number of rsectors : " << number_of_blocks << endl; 03420 fileGeom << "number of crystals transaxial : " << number_of_crystals_trans << endl; 03421 fileGeom << "number of crystals axial : " << number_of_crystals_axial << endl; 03422 03423 fileGeom << ""<< endl; 03424 fileGeom << "# The 4 following parameters could be defined in arrays (SizeLayer1,SizeLayer2,SizeLayer3,etc..) if their is more than one layer"<< endl; 03425 fileGeom << "crystals size depth : " << crystal_size_depth << endl; 03426 fileGeom << "crystals size transaxial : " << crystal_size_trans << endl; 03427 fileGeom << "crystals size axial : " << crystal_size_axial << endl; 03428 fileGeom << ""<< endl; 03429 03430 // Optional fields 03431 03432 fileGeom << "rsectors first angle : " << block_first_angle << endl; 03433 fileGeom << "number of modules transaxial : " << number_of_blocks_trans << endl; 03434 fileGeom << "number of modules axial : " << number_of_blocks_axial << endl; 03435 fileGeom << "module gap transaxial : " << block_step_trans << endl; 03436 fileGeom << "module gap axial : " << block_step_axial << endl; 03437 fileGeom << "crystal gap transaxial : " << crystal_step_trans << endl; 03438 fileGeom << "crystal gap axial : " << crystal_step_axial << endl; 03439 fileGeom << ""<< endl; 03440 03441 // Write only if different from default values 03442 if(min_angle_diff > 0.) fileGeom << "min angle difference : " << min_angle_diff << endl; 03443 if(block_angular_span != 360.) fileGeom << "rsectors angular span : " << block_angular_span << endl; 03444 if(block_nb_zshifts > 0) fileGeom << "rsectors nbZShift :" << block_nb_zshifts << endl; 03445 if(!vec_block_Z_Shift.empty()) WriteVector(fileGeom, "rsectors ZShift : ", vec_block_Z_Shift); 03446 03447 fileGeom.close(); 03448 03449 Cout("Output geom file written at :" << a_pathGeom << endl); 03450 } 03451 else 03452 { 03453 Cerr("***** dataConversionUtilities::CreateGeomWithECAT()-> Couldn't open geom file for writing "<< a_pathGeom << " !" << endl); 03454 return 1; 03455 } 03456 03457 return 0; 03458 } 03459 03460 03461 03462 03463 // ===================================================================== 03464 // --------------------------------------------------------------------- 03465 // --------------------------------------------------------------------- 03466 // ===================================================================== 03467 /* 03468 \fn CreateGeomWithSPECT() 03469 \param a_pathMac : string containing the path to a GATE macro file 03470 \param a_pathGeom : string containing the path to a CASToR output geom file 03471 \brief Read a GATE macro file containing the description of a SPECThead system, and convert it to a geom file 03472 \return 0 if success, positive value otherwise 03473 */ 03474 int CreateGeomWithSPECT(string a_pathMac, string a_pathGeom) 03475 { 03476 /* Declaring variables */ 03477 string modality; 03478 string scanner_name = GetFileFromPath(a_pathGeom.substr(0,a_pathGeom.find_first_of(".geom"))); 03479 FLTNB scanner_radius = -1.; 03480 string description = "SPECT camera extracted from GATE macro: " + a_pathMac; 03481 03482 // head(s) 03483 string head_name = "SPECThead"; // not used. Correspond to SPECThead 03484 uint32_t number_of_heads = 0; 03485 FLTNB head_pos_X = 0.; 03486 FLTNB head_pos_Y = 0.; 03487 FLTNB head_first_angle = 0.; 03488 FLTNB head_angular_pitch = -1.; 03489 string head_orbit_name = ""; 03490 FLTNB head_rotation_speed = 0.; 03491 bool is_head_Y_axis = false; 03492 03493 // crystal 03494 string crystal_name = "crystal"; 03495 FLTNB crystal_size_trans = 0.; 03496 FLTNB crystal_size_axial = 0.; 03497 FLTNB crystal_depth = 0; 03498 03499 // pixel 03500 string pixel_name = "pixel"; 03501 uint32_t number_of_pixels_trans = 1; 03502 uint32_t number_of_pixels_axial = 1; 03503 FLTNB pix_size_trans = 0.; 03504 FLTNB pix_size_axial = 0.; 03505 FLTNB pix_step_trans = 0.; 03506 FLTNB pix_step_axial = 0.; 03507 03508 // collimator parameters 03509 string focal_model_trs = "constant"; 03510 uint16_t nb_coeff_model_trs = 1; 03511 FLTNB coeff_model_trs = 0.; 03512 string focal_model_axl = "constant"; 03513 uint16_t nb_coeff_model_axl = 1; 03514 FLTNB coeff_model_axl = 0.; 03515 03516 // others 03517 uint32_t voxels_number_trans; 03518 uint32_t voxels_number_axial; 03519 double fov_trans; 03520 double fov_axial; 03521 03522 vector<string> path_mac_files; 03523 path_mac_files.push_back(a_pathMac); 03524 03525 // Recover path to all macro files from the main mac file 03526 if(GetGATEMacFiles(a_pathMac , path_mac_files)) 03527 { 03528 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Error occured when trying to recover paths to GATE macro files !" << endl); 03529 return 1; 03530 } 03531 03532 03533 // Recover aliases of the different parts of the architecture 03534 if(GetGATEAliasesSPECT(path_mac_files, head_name, crystal_name, pixel_name, 2) ) 03535 { 03536 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Error occured when trying to recover aliases for the elements of the SPECThead !" << endl); 03537 return 1; 03538 } 03539 03540 // Loop to recover all other informations 03541 for(uint16_t f=0 ; f<path_mac_files.size() ; f++) 03542 { 03543 ifstream systemMac(path_mac_files[f].c_str(), ios::in); 03544 03545 string line; 03546 while(getline(systemMac, line)) 03547 { 03548 // Scanner 03549 modality = "SPECT_CONVERGENT"; 03550 vector <string> values; 03551 string entry = ""; 03552 03553 // assumes that first block is positionned on the x-axis 03554 entry = "/gate/"+head_name+"/placement/setTranslation"; 03555 values = CheckGATECommand(entry, line); 03556 if (values.size()>0) 03557 { 03558 if(ConvertFromString(values[0], &head_pos_X) || 03559 ConvertFromString(values[1], &head_pos_Y) ) 03560 { 03561 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03562 return 1; 03563 } 03564 03565 // Check first block cartesian coordinates is 0 on X or Y axis 03566 // Throw error otherwise 03567 if(head_pos_X!=0 && head_pos_Y !=0) 03568 { 03569 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03570 Cerr(" Block cartesian coordinates on either the X or Y axis expected to be equal to 0 "<< endl); 03571 return 1; 03572 } 03573 03574 // Get the axis on which the first rsector is positionned 03575 if(head_pos_Y!=0) is_head_Y_axis = true; 03576 03577 scanner_radius = is_head_Y_axis ? abs(head_pos_Y) : abs(head_pos_X) ; 03578 } 03579 03580 03581 entry = "/gate/"+head_name+"/geometry/setZLength"; 03582 values = CheckGATECommand(entry, line); 03583 if (values.size()>0) 03584 { 03585 if(ConvertFromString(values[0], &fov_axial)) 03586 { 03587 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03588 return 1; 03589 } 03590 // 4mm voxels by default 03591 voxels_number_axial = fov_axial/4 + 1; 03592 } 03593 03594 03595 entry = "/gate/"+head_name+"/ring/setRepeatNumber"; 03596 values = CheckGATECommand(entry, line); 03597 if (values.size()>0) 03598 { 03599 if(ConvertFromString(values[0], &number_of_heads)) 03600 { 03601 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03602 return 1; 03603 } 03604 } 03605 03606 03607 entry = "/gate/"+head_name+"/ring/setFirstAngle"; 03608 values = CheckGATECommand(entry, line); 03609 if (values.size()>0) 03610 { 03611 if(ConvertFromString(values[0], &head_first_angle)) 03612 { 03613 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03614 return 1; 03615 } 03616 } 03617 03618 03619 entry = "/gate/"+head_name+"/ring/setAngularPitch"; 03620 values = CheckGATECommand(entry, line); 03621 if (values.size()>0) 03622 { 03623 if(ConvertFromString(values[0], &head_angular_pitch)) 03624 { 03625 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03626 return 1; 03627 } 03628 } 03629 03630 entry = "/gate/"+head_name+"/moves/insert"; 03631 values = CheckGATECommand(entry, line); 03632 if (values.size()>0) 03633 { 03634 if(ConvertFromString(values[0], &head_orbit_name)) 03635 { 03636 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03637 return 1; 03638 } 03639 } 03640 03641 entry = "/gate/"+head_orbit_name+"/setSpeed"; 03642 values = CheckGATECommand(entry, line); 03643 if (values.size()>0) 03644 { 03645 if(ConvertFromString(values[0], &head_rotation_speed)) 03646 { 03647 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Conversion error occured while trying to parse line: " << line<< endl); 03648 return 1; 03649 } 03650 } 03651 03652 03653 03654 // --- Crystals --- 03655 entry = "/gate/"+crystal_name+"/geometry/setXLength"; 03656 values = CheckGATECommand(entry, line); 03657 if (values.size()>0) 03658 { 03659 double x_length; 03660 if(ConvertFromString(values[0], &x_length) ) 03661 { 03662 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03663 return 1; 03664 } 03665 03666 if (is_head_Y_axis) 03667 crystal_size_trans = x_length; 03668 else 03669 crystal_depth = x_length; 03670 } 03671 03672 entry = "/gate/"+crystal_name+"/geometry/setYLength"; 03673 values = CheckGATECommand(entry, line); 03674 if (values.size()>0) 03675 { 03676 double y_length; 03677 if(ConvertFromString(values[0], &y_length) ) 03678 { 03679 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03680 return 1; 03681 } 03682 03683 if (is_head_Y_axis) 03684 crystal_depth = y_length; 03685 else 03686 crystal_size_trans = y_length; 03687 03688 } 03689 03690 entry = "/gate/"+crystal_name+"/geometry/setZLength"; 03691 values = CheckGATECommand(entry, line); 03692 if (values.size()>0) 03693 { 03694 if(ConvertFromString(values[0], &crystal_size_axial) ) 03695 { 03696 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03697 return 1; 03698 } 03699 } 03700 03701 03702 03703 03704 03705 // --- Pixels --- 03706 entry = is_head_Y_axis ? 03707 "/gate/"+pixel_name+"/cubicArray/setRepeatNumberX": 03708 "/gate/"+pixel_name+"/cubicArray/setRepeatNumberY"; 03709 03710 03711 values = CheckGATECommand(entry, line); 03712 if (values.size()>0) 03713 { 03714 if(ConvertFromString(values[0], &number_of_pixels_trans) ) 03715 { 03716 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03717 return 1; 03718 } 03719 } 03720 03721 entry = "/gate/"+pixel_name+"/cubicArray/setRepeatNumberZ"; 03722 values = CheckGATECommand(entry, line); 03723 if (values.size()>0) 03724 { 03725 if(ConvertFromString(values[0], &number_of_pixels_axial) ) 03726 { 03727 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03728 return 1; 03729 } 03730 } 03731 03732 03733 entry = "/gate/"+pixel_name+"/geometry/setXLength"; 03734 values = CheckGATECommand(entry, line); 03735 if (values.size()>0) 03736 { 03737 double x_length; 03738 if(ConvertFromString(values[0], &x_length) ) 03739 { 03740 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03741 return 1; 03742 } 03743 03744 if (is_head_Y_axis) 03745 pix_size_trans = x_length; 03746 } 03747 03748 entry = "/gate/"+pixel_name+"/geometry/setYLength"; 03749 values = CheckGATECommand(entry, line); 03750 if (values.size()>0) 03751 { 03752 double y_length; 03753 if(ConvertFromString(values[0], &y_length) ) 03754 { 03755 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03756 return 1; 03757 } 03758 03759 if (!is_head_Y_axis) 03760 pix_size_trans = y_length; 03761 03762 } 03763 03764 entry = "/gate/"+pixel_name+"/geometry/setZLength"; 03765 values = CheckGATECommand(entry, line); 03766 if (values.size()>0) 03767 { 03768 if(ConvertFromString(values[0], &pix_size_axial) ) 03769 { 03770 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03771 return 1; 03772 } 03773 } 03774 03775 03776 entry = "/gate/"+pixel_name+"/cubicArray/setRepeatVector"; 03777 values = CheckGATECommand(entry, line); 03778 if (values.size()>0) 03779 { 03780 string trs_step = is_head_Y_axis ? 03781 values[0] : 03782 values[1] ; 03783 03784 if(ConvertFromString(trs_step, &pix_step_trans) || 03785 ConvertFromString(values[2], &pix_step_axial) ) 03786 { 03787 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03788 return 1; 03789 } 03790 } 03791 03792 03793 // collimator parameter 03794 entry = is_head_Y_axis ? 03795 "/gate/fanbeam/geometry/setFocalDistanceY": 03796 "/gate/fanbeam/geometry/setFocalDistanceX"; 03797 03798 entry = "/gate/fanbeam/geometry/setFocalDistanceX"; 03799 values = CheckGATECommand(entry, line); 03800 if (values.size()>0) 03801 { 03802 if(ConvertFromString(values[0], &coeff_model_trs) ) 03803 { 03804 Cerr("***** dataConversionUtilities::CreateGeomWithCylindrical()-> Conversion error occured while trying to parse line: " << line<< endl); 03805 return 1; 03806 } 03807 03808 focal_model_trs = "polynomial"; 03809 } 03810 03811 } 03812 systemMac.close(); 03813 } 03814 03815 // Transaxial FOV defined as scanner radius / 2 03816 fov_trans = scanner_radius/2; 03817 // 4mm voxels by default 03818 voxels_number_trans = fov_trans/2 + 1; 03819 03820 uint32_t nb_pixels = number_of_pixels_axial * number_of_pixels_trans; 03821 03822 pix_size_axial = nb_pixels>1 ? pix_size_axial : crystal_size_axial; 03823 pix_size_trans = nb_pixels>1 ? pix_size_trans : crystal_size_trans; 03824 03825 // Compute gaps 03826 if (pix_step_axial - pix_size_axial >= 0) 03827 pix_step_axial = pix_step_axial - pix_size_axial; 03828 03829 if (pix_step_trans - pix_size_trans >= 0) 03830 pix_step_trans = pix_step_trans - pix_size_trans; 03831 03832 03833 // CASToR x and y axis for rotation inverted in comparison with GATE 03834 head_first_angle = round(atan2f(head_pos_X , head_pos_Y) * 180. / M_PI) 03835 - head_first_angle; 03836 03837 03838 // Writing the .geom file 03839 ofstream fileGeom(a_pathGeom.c_str(), ios::out | ios::trunc); 03840 03841 if(fileGeom) 03842 { 03843 fileGeom << "modality : " << modality << endl; 03844 fileGeom << "scanner name : " << scanner_name << endl; 03845 fileGeom << "number of detector heads : " << number_of_heads << endl; 03846 fileGeom << "trans number of pixels : " << number_of_pixels_trans << endl; 03847 fileGeom << "trans pixel size : " << pix_size_trans << endl; 03848 fileGeom << "trans gap size : " << pix_step_trans << endl; 03849 fileGeom << "axial number of pixels : " << number_of_pixels_axial << endl; 03850 fileGeom << "axial pixel size : " << pix_size_axial << endl; 03851 fileGeom << "axial gap size : " << pix_step_axial << endl; 03852 03853 fileGeom << "detector depth : " << crystal_depth << endl; 03854 03855 fileGeom << "scanner radius : " << scanner_radius; 03856 for(size_t h=1 ; h<number_of_heads ; h++) 03857 fileGeom << "," << scanner_radius; 03858 fileGeom << endl; 03859 03860 fileGeom << "# Collimator configuration : "<< endl << endl; 03861 for(size_t h=0 ; h<number_of_heads ; h++) 03862 { 03863 fileGeom << "head" << h+1 << ":" << endl; 03864 fileGeom << "trans focal model: " << focal_model_trs << endl; 03865 fileGeom << "trans number of coef model: " << nb_coeff_model_trs << endl; 03866 fileGeom << "trans parameters: " << coeff_model_trs << endl; 03867 fileGeom << "axial focal model: " << focal_model_axl << endl; 03868 fileGeom << "axial number of coef model: " << nb_coeff_model_axl << endl; 03869 fileGeom << "axial parameters: " << coeff_model_axl << endl; 03870 fileGeom << endl; 03871 } 03872 03873 fileGeom << "" << endl; 03874 fileGeom << "# default reconstruction parameters" << endl; 03875 fileGeom << "voxels number transaxial : " << voxels_number_trans << endl; 03876 fileGeom << "voxels number axial : " << voxels_number_axial << endl; 03877 03878 fileGeom << "field of view transaxial : " << fov_trans << endl; 03879 fileGeom << "field of view axial : " << fov_axial << endl << endl ; 03880 fileGeom << ""<< endl; 03881 03882 fileGeom << "# description" << endl; 03883 fileGeom << "description : " << description << endl; 03884 03885 fileGeom.close(); 03886 03887 Cout("Output geom file written at :" << a_pathGeom << endl); 03888 } 03889 else 03890 { 03891 Cerr("***** dataConversionUtilities::CreateGeomWithSPECT()-> Couldn't open geom file for writing "<< a_pathGeom << " !" << endl); 03892 return 1; 03893 } 03894 03895 return 0; 03896 }