![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oInterfileIO 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: none 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00012 00019 #include "oInterfileIO.hh" 00020 #include <iomanip> 00021 00022 int User_Endianness = -1; 00025 // ===================================================================== 00026 // --------------------------------------------------------------------- 00027 // --------------------------------------------------------------------- 00028 // ===================================================================== 00029 /* 00030 \fn IntfKeyGetValueFromFile 00031 \param a_pathToHeader : path to the interfile header 00032 \param a_key : the key to recover 00033 \param T* ap_return : template array in which the data will be recovered 00034 \param int a_nbElts : number of elements to recover 00035 \param int a_mandatoryFlag : flag indicating if the data to recover if mandatory (true) or optionnal (false) 00036 \brief Look for "a_nbElts" elts in the "a_pathToHeader" interfile header matching the "a_keyword" key 00037 passed as parameter and return the corresponding value(s) in the "ap_return" templated array. 00038 \details If more than one elements are to be recovered, the function first check that 00039 the key has a correct Interfile key layout (brackets and commas : {,,}) 00040 Depending on the mandatoryFlag, the function will return an error (flag > 0) 00041 or a warning (flag = 0) if the key is not found 00042 \return 0 if success, and positive value otherwise (1 if error, 2 if key not found). 00043 */ 00044 template<class T> 00045 int IntfKeyGetValueFromFile(const string& a_pathToHeader, 00046 const string& a_key, 00047 T* ap_return, 00048 int a_nbElts, 00049 int a_mandatoryFlag) 00050 { 00051 ifstream input_file(a_pathToHeader.c_str(), ios::in); 00052 string line; 00053 00054 // Check file 00055 if (input_file) 00056 { 00057 while(!input_file.eof()) 00058 { 00059 getline(input_file, line); 00060 00061 if(line.empty()) 00062 continue; 00063 00064 Intf_key Key; 00065 00066 // Process the Key at this line 00067 if (IntfRecoverKey(&Key, line) ) 00068 { 00069 Cerr("***** IntfKeyGetValueFromFile()-> Error : couldn't correctly read interfile key '"<< line << "' !" << endl); 00070 return 1; 00071 } 00072 //if (line.find(a_keyword) != string::npos) 00073 if(IntfCheckKeyMatch(Key, a_key)) 00074 { 00075 00076 if(a_nbElts == 1) // 1 elt is required, just return the result 00077 { 00078 if (ConvertFromString(Key.kvalue, &ap_return[0])) 00079 { 00080 Cerr("***** IntfKeyGetValueFromFile()-> Exception when trying to read tag '" << a_key << "' in file '" << a_pathToHeader << "'." << endl); 00081 return 1; 00082 } 00083 00084 return 0; 00085 } 00086 else // Several elements required 00087 { 00088 // First check we have an array 00089 if (!IntfKeyIsArray(Key)) 00090 { 00091 Cerr("***** IntfKeyGetValueFromFile() -> " << a_nbElts << " requested for interfile key " << a_key << " , but this key is not an array !" << endl); 00092 return 1; 00093 } 00094 else 00095 { 00096 // Check the number of elements in the key. 00097 if(IntfKeyGetArrayNbElts(Key) != a_nbElts) 00098 { 00099 Cerr("***** IntfKeyGetValueFromFile() -> Nb of elements to recover (=" << a_nbElts << ") does not correspond to the number of elements found in the key '" 00100 << a_key << "' (" << IntfKeyGetArrayNbElts(Key) << ") !" << endl); 00101 return 1; 00102 } 00103 00104 // Read array key 00105 if (IntfKeyGetArrayElts(Key, ap_return) ) 00106 { 00107 Cerr("***** IntfKeyGetValueFromFile() -> " << a_nbElts << " requested for interfile key " << a_key << " , but this key is not an array !" << endl); 00108 return 1; 00109 } 00110 00111 return 0; 00112 } 00113 } 00114 } 00115 } 00116 00117 // Tag not found, throw an error message if the tag is mandatory 00118 if (a_mandatoryFlag > 0) 00119 { 00120 Cerr("***** IntfKeyGetValueFromFile()-> Error when reading Interfile '" << a_pathToHeader << "'. Key '" << a_key << "' was not found." << endl); 00121 return KEYWORD_MANDATORY_NOT_FOUND; 00122 } 00123 else 00124 { 00125 return KEYWORD_OPTIONAL_NOT_FOUND; 00126 } 00127 00128 } 00129 else 00130 { 00131 Cerr("***** IntfKeyGetValueFromFile()-> Couldn't find or read data-file '"<< a_pathToHeader << "' !" << endl); 00132 return 1; 00133 } 00134 } 00135 00136 // Templated functions definitions 00137 template int IntfKeyGetValueFromFile<string>(const string& a_pathToHeader, const string& a_key, string* ap_return, int a_nbElts, int a_mandatoryFlag); 00138 template int IntfKeyGetValueFromFile<int>(const string& a_pathToHeader, const string& a_key, int* ap_return, int a_nbElts, int a_mandatoryFlag); 00139 template int IntfKeyGetValueFromFile<int64_t>(const string& a_pathToHeader, const string& a_key, int64_t* ap_return, int a_nbElts, int a_mandatoryFlag); 00140 template int IntfKeyGetValueFromFile<float>(const string& a_pathToHeader, const string& a_key, float* ap_return, int a_nbElts, int a_mandatoryFlag); 00141 template int IntfKeyGetValueFromFile<double>(const string& a_pathToHeader, const string& a_key, double* ap_return, int a_nbElts, int a_mandatoryFlag); 00142 template int IntfKeyGetValueFromFile<long double>(const string& a_pathToHeader, const string& a_key, long double* ap_return, int a_nbElts, int a_mandatoryFlag); 00143 template int IntfKeyGetValueFromFile<uint8_t>(const string& a_pathToHeader, const string& a_key, uint8_t* ap_return, int a_nbElts, int a_mandatoryFlag); 00144 template int IntfKeyGetValueFromFile<uint16_t>(const string& a_pathToHeader, const string& a_key, uint16_t* ap_return, int a_nbElts, int a_mandatoryFlag); 00145 template int IntfKeyGetValueFromFile<uint32_t>(const string& a_pathToHeader, const string& a_key, uint32_t* ap_return, int a_nbElts, int a_mandatoryFlag); 00146 template int IntfKeyGetValueFromFile<bool>(const string& a_pathToHeader, const string& a_key, bool* ap_return, int a_nbElts, int a_mandatoryFlag); 00147 00148 00149 /* 00150 \fn IntfKeyGetValueFromFile(const string& a_pathToHeader, const string& a_key, T* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences) 00151 \param a_pathToHeader : path to the interfile header 00152 \param a_key : the key to recover 00153 \param T* ap_return : template array in which the data will be recovered 00154 \param int a_nbElts : number of elements to recover 00155 \param int a_mandatoryFlag : flag indicating if the data to recover if mandatory (true) or optionnal (false) 00156 \param int a_nbOccurences : number of occurences of the field before recovering the value 00157 \brief Look for "a_nbElts" elts in the "a_pathToHeader" interfile header matching the "a_keyword" key passed as parameter and return the corresponding value(s) in the "ap_return" templated array.\n 00158 Parameter "a_nbOccurences" can be set to recover a specific occurrence of a recurring value 00159 \details If more than one elements are to be recovered, the function first check the key has a correct Interfile kay layout (brackets and commas : {,,})\n 00160 Depending on the mandatoryFlag, the function will return an error (flag > 0) or a warning (flag = 0) if the key is not found 00161 \return 0 if success, and positive value otherwise (1 if error, 2 if key not found). 00162 */ 00163 template <typename T> int IntfKeyGetRecurringValueFromFile(const string& a_pathToHeader, 00164 const string& a_key, 00165 T* ap_return, 00166 int a_nbElts, 00167 int a_mandatoryFlag, 00168 uint16_t a_nbOccurrences) 00169 { 00170 ifstream input_file(a_pathToHeader.c_str(), ios::in); 00171 string line; 00172 uint16_t nb_occurences_cur =0; 00173 00174 // Check file 00175 if (input_file) 00176 { 00177 while(!input_file.eof()) 00178 { 00179 getline(input_file, line); 00180 00181 if(line.empty()) 00182 continue; 00183 00184 Intf_key Key; 00185 00186 // Process the Key at this line 00187 if (IntfRecoverKey(&Key, line) ) 00188 { 00189 Cerr("***** IntfKeyGetValueFromFile()-> Error : couldn't correctly read interfile key '"<< line << "' !" << endl); 00190 return 1; 00191 } 00192 00193 00194 00195 //if (line.find(a_keyword) != string::npos) 00196 if(IntfCheckKeyMatch(Key, a_key)) 00197 { 00198 // Check if we reached the correct number of occurence of the key 00199 // Skip otherwise 00200 if(nb_occurences_cur < a_nbOccurrences) 00201 { 00202 nb_occurences_cur++; 00203 continue; 00204 } 00205 00206 if(a_nbElts == 1) // 1 elt is required, just return the result 00207 { 00208 if (ConvertFromString(Key.kvalue, &ap_return[0])) 00209 { 00210 Cerr("***** IntfKeyGetValueFromFile()-> Exception when trying to read tag '" << a_key << "' in file '" << a_pathToHeader << "'." << endl); 00211 return 1; 00212 } 00213 00214 return 0; 00215 } 00216 else // Several elements required 00217 { 00218 // First check we have an array 00219 if (!IntfKeyIsArray(Key)) 00220 { 00221 Cerr("***** IntfKeyGetValueFromFile() -> " << a_nbElts << " requested for interfile key " << a_key << " , but this key is not an array !" << endl); 00222 return 1; 00223 } 00224 else 00225 { 00226 // Check the number of elements in the key. 00227 if(IntfKeyGetArrayNbElts(Key) != a_nbElts) 00228 { 00229 Cerr("***** IntfKeyGetValueFromFile() -> Nb of elements to recover (=" << a_nbElts << ") does not correspond to the number of elements found in the key '" 00230 << a_key << "' (" << IntfKeyGetArrayNbElts(Key) << ") !" << endl); 00231 return 1; 00232 } 00233 00234 // Read array key 00235 if (IntfKeyGetArrayElts(Key, ap_return) ) 00236 { 00237 Cerr("***** IntfKeyGetValueFromFile() -> " << a_nbElts << " requested for interfile key " << a_key << " , but this key is not an array !" << endl); 00238 return 1; 00239 } 00240 00241 return 0; 00242 } 00243 } 00244 } 00245 } 00246 00247 // Tag not found, throw an error message if the tag is mandatory 00248 if (a_mandatoryFlag > 0) 00249 { 00250 Cerr("***** IntfKeyGetValueFromFile()-> Error when reading Interfile '" << a_pathToHeader << "'. Key '" << a_key << "' was not found." << endl); 00251 return KEYWORD_MANDATORY_NOT_FOUND; 00252 } 00253 else 00254 { 00255 return KEYWORD_OPTIONAL_NOT_FOUND; 00256 } 00257 00258 } 00259 else 00260 { 00261 Cerr("***** IntfKeyGetValueFromFile()-> Couldn't find or read data-file '"<< a_pathToHeader << "' !" << endl); 00262 return 1; 00263 } 00264 } 00265 00266 // Templated functions definitions 00267 template int IntfKeyGetRecurringValueFromFile<string>(const string& a_pathToHeader, const string& a_key, string* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00268 template int IntfKeyGetRecurringValueFromFile<int>(const string& a_pathToHeader, const string& a_key, int* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00269 template int IntfKeyGetRecurringValueFromFile<int64_t>(const string& a_pathToHeader, const string& a_key, int64_t* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00270 template int IntfKeyGetRecurringValueFromFile<float>(const string& a_pathToHeader, const string& a_key, float* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00271 template int IntfKeyGetRecurringValueFromFile<double>(const string& a_pathToHeader, const string& a_key, double* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00272 template int IntfKeyGetRecurringValueFromFile<long double>(const string& a_pathToHeader, const string& a_key, long double* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00273 template int IntfKeyGetRecurringValueFromFile<uint8_t>(const string& a_pathToHeader, const string& a_key, uint8_t* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00274 template int IntfKeyGetRecurringValueFromFile<uint16_t>(const string& a_pathToHeader, const string& a_key, uint16_t* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00275 template int IntfKeyGetRecurringValueFromFile<uint32_t>(const string& a_pathToHeader, const string& a_key, uint32_t* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00276 template int IntfKeyGetRecurringValueFromFile<bool>(const string& a_pathToHeader, const string& a_key, bool* ap_return, int a_nbElts, int a_mandatoryFlag, uint16_t a_nbOccurrences); 00277 00278 00279 00280 00281 // ===================================================================== 00282 // --------------------------------------------------------------------- 00283 // --------------------------------------------------------------------- 00284 // ===================================================================== 00285 /* 00286 \fn IntfReadImage(const string& a_pathToHeaderFile, FLTNB* ap_ImgMatrix, Intf_fields* ap_IF, int vb, bool a_lerpFlag) 00287 \param a_pathToHeaderFile : path to the header file 00288 \param ap_ImgMatrix : 1 dimensional image matrix which will recover the image. 00289 \param ap_IF : Intf_fields structure containing image metadata 00290 \param vb : Verbosity level 00291 \param a_lerpFlag : if true, enable linear interpolation of the image if img dimensions differ from the reconstruction dimensions 00292 \brief Main function dedicated to Interfile 3D image loading \n 00293 Get image information from a provided Intf_fields structure. 00294 \details Call the main functions dedicated to Interfile reading : IntfReadHeader(), IntfCheckConsistency(), then IntfReadImage() 00295 \return 0 if success, positive value otherwise. 00296 */ 00297 int IntfReadProjectionImage( const string& a_pathToHeaderFile, 00298 FLTNB* ap_ImgMatrix, 00299 Intf_fields* ap_IF, 00300 int vb, 00301 bool a_lerpFlag) 00302 { 00303 if(vb >= 3) Cout("IntfReadProjectionImage()-> Read Interfile header : "<< a_pathToHeaderFile << endl); 00304 00305 // Init Interfile Key structure 00306 IntfKeyInitFields(ap_IF); 00307 00308 // Recover image infos from the header file 00309 if(IntfReadHeader(a_pathToHeaderFile, ap_IF, vb) ) 00310 { 00311 Cerr("***** IntfReadProjectionImage()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00312 return 1; 00313 } 00314 00315 // Specific changes for projection data. //TODO : deal with that in IntfReadHeader, specification regarding the nature of the data 00316 ap_IF->mtx_size[2] = ap_IF->nb_total_imgs; 00317 00318 int nb_tot_pixels = ap_IF->mtx_size[0] 00319 *ap_IF->mtx_size[1] 00320 *ap_IF->mtx_size[2]; 00321 00322 // Read image data 00323 ifstream img_file(ap_IF->path_to_image.c_str(), ios::binary | ios::in); 00324 if(img_file) 00325 { 00326 // Get the position of the beginning of the image data 00327 uint32_t offset = ap_IF->data_offset; 00328 00329 // Call the right IntfReadData() function according to the pixel type, and read data 00330 IntfGetPixelTypeAndReadData(*ap_IF, &img_file, ap_ImgMatrix, NULL, &offset, nb_tot_pixels, vb); 00331 } 00332 else 00333 { 00334 Cerr("***** IntfReadProjectionImage()-> Error occurred while trying to read the image file at the path: '"<< ap_IF->path_to_image << "' !" << endl); 00335 return 1; 00336 } 00337 00338 return 0; 00339 } 00340 00341 00342 00343 00344 00345 00346 00347 00348 00349 00350 // ===================================================================== 00351 // --------------------------------------------------------------------- 00352 // --------------------------------------------------------------------- 00353 // ===================================================================== 00354 /* 00355 \fn IntfReadImage 00356 \param a_pathToHeaderFile : path to the header file 00357 \param ap_ImgMatrix : 1 dimensional image matrix which will recover the image. 00358 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 00359 \param vb : Verbosity level 00360 \param a_lerpFlag : if true, enable linear interpolation of the image if img dimensions differ from the reconstruction dimensions 00361 \brief Main function dedicated to Interfile 3D image loading 00362 \details Call the main functions dedicated to Interfile reading : IntfReadHeader(), IntfCheckConsistency(), then IntfReadImage() 00363 \return 0 if success, positive value otherwise. 00364 */ 00365 int IntfReadImage( const string& a_pathToHeaderFile, 00366 FLTNB* ap_ImgMatrix, 00367 oImageDimensionsAndQuantification* ap_ID, 00368 int vb, 00369 bool a_lerpFlag) 00370 { 00371 if(vb >= 3) Cout("IntfReadImage()-> Read Interfile header : "<< a_pathToHeaderFile << endl); 00372 00373 // Init Interfile Key structure 00374 Intf_fields Img_fields; 00375 IntfKeyInitFields(&Img_fields); 00376 00377 // Recover image infos from the header file 00378 if(IntfReadHeader(a_pathToHeaderFile, &Img_fields, vb) ) 00379 { 00380 Cerr("***** IntfReadImage()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00381 return 1; 00382 } 00383 00384 // Check consistency between image interfile data and the algorithm data 00385 if(IntfCheckConsistency(&Img_fields, ap_ID, vb, a_lerpFlag) ) 00386 { 00387 Cerr("***** IntfReadImage()-> Error : while checking consistencies between the reconstruction parameters and the interfile keys in the header '"<< a_pathToHeaderFile << "' !" << endl); 00388 return 1; 00389 } 00390 00391 // If interpolation required, allocate matrix with original image dimensions. 00392 // Todo : find a way to free memory if error occurs in following functions 00393 FLTNB* pimg_erp = NULL; 00394 IntfAllocInterpImg(&pimg_erp, Img_fields); 00395 00396 // Read image data 00397 ifstream img_file(Img_fields.path_to_image.c_str(), ios::binary | ios::in); 00398 if(img_file) 00399 { 00400 // Get the position of the beginning of the image data 00401 uint32_t offset = Img_fields.data_offset; 00402 00403 // Call the right IntfReadData() function according to the pixel type, and read data 00404 IntfGetPixelTypeAndReadData(Img_fields, &img_file, ap_ImgMatrix, pimg_erp, &offset, ap_ID->GetNbVoxXYZ(), vb); 00405 } 00406 else 00407 { 00408 Cerr("***** IntfReadImage()-> Error occurred while trying to read the image file at the path: '"<< Img_fields.path_to_image << "' !" << endl); 00409 return 1; 00410 } 00411 00412 // If interpolation was required, deallocate image matrix 00413 if(pimg_erp) delete[] pimg_erp; 00414 00415 return 0; 00416 } 00417 00418 00419 00420 00421 // ===================================================================== 00422 // --------------------------------------------------------------------- 00423 // --------------------------------------------------------------------- 00424 // ===================================================================== 00425 /* 00426 \fn IntfReadImage 00427 \param a_pathToHeaderFile : path to the main header file 00428 \param a4p_ImgMatrix : 4 dimensional image matrix which will recover the image. 00429 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 00430 \param vb : Verbosity level 00431 \param a_lerpFlag : if true, enable linear interpolation of the image if img dimensions differ from the reconstruction dimensions 00432 \brief Main function dedicated to Interfile 5D (1D+1D+1D time + 3D) image loading 00433 \details Check is the main header file is a metaheader associated to several image files, or a unique interfile header 00434 Depending on the type of file input (metaheader or unique file), read the group of image files or the unique provided image file 00435 \todo : Check if image 3D dimensions are different from an image to another ? 00436 (very unlikely, but it would cause segfault if interpolation is enabled) 00437 \return 0 if success, positive value otherwise. 00438 */ 00439 int IntfReadImage( const string& a_pathToHeaderFile, 00440 FLTNB**** a4p_ImgMatrix, 00441 oImageDimensionsAndQuantification* ap_ID, 00442 int vb, 00443 bool a_lerpFlag) 00444 { 00445 if(vb >= 5) Cout("IntfReadImage()****" << endl); 00446 00447 // Init Interfile Key structure 00448 Intf_fields Img_fields; 00449 IntfKeyInitFields(&Img_fields); 00450 00451 // List containing single/multiple interfile headers of the images to read 00452 vector<string> lpath_to_headers; 00453 00454 // First check if the provided file is a metaheader with multiple files 00455 // and recover the list of header file 00456 if(IntfIsMHD(a_pathToHeaderFile, lpath_to_headers) < 0 ) // error 00457 { 00458 Cerr("***** IntfReadImage()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00459 return 1; 00460 } 00461 00462 // --- Case 1 : We have one single data file containing all the data --- // 00463 if(lpath_to_headers.size() == 1) 00464 { 00465 // Recover image infos from either metaheader or single header file 00466 if(IntfReadHeader(a_pathToHeaderFile, &Img_fields, vb) ) 00467 { 00468 Cerr("***** IntfReadImage()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00469 return 1; 00470 } 00471 00472 // Check consistency between image interfile data and the algorithm data 00473 if(IntfCheckConsistency(&Img_fields, ap_ID, vb, a_lerpFlag) ) 00474 { 00475 Cerr("***** IntfReadImage()-> Error : while checking consistencies between the reconstruction parameters and the interfile keys in the header '"<< a_pathToHeaderFile << "' !" << endl); 00476 return 1; 00477 } 00478 00479 00480 00481 // If interpolation required, instanciate matrix with original image dimensions. 00482 // Todo : find a way to free memory if error occurs in following functions 00483 FLTNB* pimg_erp = NULL; 00484 //pimg_erp = new FLTNB[ Img_fields.mtx_size[0]*Img_fields.mtx_size[1]*Img_fields.mtx_size[2] ]; 00485 IntfAllocInterpImg(&pimg_erp, Img_fields); 00486 00487 ifstream img_file(Img_fields.path_to_image.c_str(), ios::binary | ios::in); 00488 00489 if(img_file) 00490 { 00491 // Get the position of the beginning of the image data 00492 uint32_t offset = Img_fields.data_offset; 00493 00494 // Call the right IntfReadData() function according to the pixel type, and read data. 00495 // offset is updated with the called functions within the loops 00496 00497 for(int d1=0 ; d1<ap_ID->GetNbTimeFrames() ; d1++) 00498 for(int d2=0 ; d2<ap_ID->GetNbRespGates() ; d2++) 00499 for(int d3=0 ; d3<ap_ID->GetNbCardGates() ; d3++) 00500 IntfGetPixelTypeAndReadData(Img_fields, &img_file, a4p_ImgMatrix[d1][d2][d3], pimg_erp, &offset, ap_ID->GetNbVoxXYZ(), vb); 00501 } 00502 else 00503 { 00504 Cerr("***** IntfReadImage()-> Error occurred while trying to read the image file at the path: '"<< Img_fields.path_to_image << "' !" << endl); 00505 return 1; 00506 } 00507 00508 // If interpolation was required, deallocate image matrix with original image dimensions 00509 if(pimg_erp) delete[] pimg_erp; 00510 } 00511 00512 00513 00514 00515 // --- Case 2 : We have a number of data file equal to the number of image to load --- // 00516 else 00517 { 00518 // Recover image infos from the metaheader 00519 if(IntfReadHeader(a_pathToHeaderFile, &Img_fields, vb) ) 00520 { 00521 Cerr("***** IntfReadImage()-> Error : while trying to read the interfile metaheader '"<< a_pathToHeaderFile << "' !" << endl); 00522 return 1; 00523 } 00524 00525 // Get dimensions from ImageDimensions object 00526 int dims[3]; 00527 dims[0] = ap_ID->GetNbTimeFrames(); 00528 dims[1] = ap_ID->GetNbRespGates(); 00529 dims[2] = ap_ID->GetNbCardGates(); 00530 00531 if(lpath_to_headers.size() != (uint32_t)(dims[0]*dims[1]*dims[2])) // Check we have a number of file corresponding to the number of images to load 00532 { 00533 Cerr("***** IntfReadImage()-> Error : number of interfile images ("<< lpath_to_headers.size() << 00534 ") not consistent with the number of images to load (" << dims[0]*dims[1]*dims[2]<< ") !" << endl); 00535 return 1; 00536 } 00537 00538 // If interpolation required, instanciate matrix with original image dimensions. 00539 // Todo : find a way to free memory if error occurs in following functions 00540 FLTNB* pimg_erp = NULL; 00541 IntfAllocInterpImg(&pimg_erp, Img_fields); 00542 00543 // Loop on dynamic image files 00544 for(int d1=0 ; d1<dims[0] ; d1++) 00545 for(int d2=0 ; d2<dims[1] ; d2++) 00546 for(int d3=0 ; d3<dims[2] ; d3++) 00547 { 00548 int idx_img = d1*dims[1]*dims[2] + d2*dims[2] + d3; 00549 00550 // Recover image infos from the specific interfile header 00551 // Todo : Check if image 3D dimensions are different from an image to another ? 00552 // (very unlikely, but it would cause segfault if interpolation is enabled) 00553 if(IntfReadHeader(lpath_to_headers[idx_img], &Img_fields, vb) ) 00554 { 00555 Cerr("***** IntfReadImage()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00556 return 1; 00557 } 00558 00559 // Check consistency between image interfile data and the algorithm data 00560 if(IntfCheckConsistency(&Img_fields, ap_ID, vb, a_lerpFlag) ) 00561 { 00562 Cerr("***** IntfReadImage()-> Error : while checking consistencies between the reconstruction parameters and the interfile keys in the header '" 00563 << a_pathToHeaderFile << "' !" << endl); 00564 return 1; 00565 } 00566 00567 00568 ifstream img_file(Img_fields.path_to_image.c_str(), ios::binary | ios::in); 00569 if(img_file) 00570 { 00571 // Get the position of the beginning of the image data (assuming that offset is the same for each file (if any) ) 00572 uint32_t offset = Img_fields.data_offset; 00573 00574 // Call the right IntfReadData() function according to the pixel type, and read data (offset not updated) 00575 IntfGetPixelTypeAndReadData(Img_fields, &img_file, a4p_ImgMatrix[d1][d2][d3], pimg_erp, &offset, ap_ID->GetNbVoxXYZ(), vb); 00576 } 00577 else 00578 { 00579 Cerr("***** IntfReadImage()-> Error occurred while trying to read the image file at the path: '"<< Img_fields.path_to_image << "' !" << endl); 00580 return 1; 00581 } 00582 } 00583 00584 // If interpolation was required, deallocate image matrix with original image dimensions 00585 if(pimg_erp) delete[] pimg_erp; 00586 } 00587 00588 /* If gating is enabled with separate 3D gated images, the image_duration key may be in each header 00589 The following check will be broken in this case. 00590 // Check some dynamic recovered infos 00591 if(Img_fields.nb_time_frames > 1 && 00592 Img_fields.nb_time_frames != Img_fields.image_duration.size()) 00593 { 00594 Cerr("***** IntfReadImage()-> Error : the number of provided image duration: '"<< Img_fields.image_duration.size() 00595 << "' does not match the number of frames '"<< Img_fields.nb_time_frames <<"' !" << endl); 00596 return 1; 00597 } 00598 */ 00599 return 0; 00600 } 00601 00602 00603 00604 00605 // ===================================================================== 00606 // --------------------------------------------------------------------- 00607 // --------------------------------------------------------------------- 00608 // ===================================================================== 00609 /* 00610 \fn IntfReadImgDynCoeffFile 00611 \param a_pathToHeaderFile : path to the header file 00612 \param a2p_ImgMatrix : 2 dimensional image matrix which will recover the image.. 00613 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 00614 \param a_nbFbases : Number of basis functions 00615 \param vb : verbosity 00616 \param a_lerpFlag : if true, enable linear interpolation of the image if img dimensions differ from the reconstruction dimensions 00617 \brief Function dedicated to Interfile image reading for dynamic coefficients images 00618 \details The total number of basis functions should be provided in parameters 00619 Check is the main header file is a metaheader associated to several image files, or a unique interfile header 00620 Depending on the type of file input (metaheader or unique file), read the group of image files or the unique provided image file 00621 \return 0 if success, positive value otherwise. 00622 */ 00623 int IntfReadImgDynCoeffFile(const string& a_pathToHeaderFile, 00624 FLTNB** a2p_ImgMatrix, 00625 oImageDimensionsAndQuantification* ap_ID, 00626 int a_nbFbases, 00627 int vb, 00628 bool a_lerpFlag) 00629 { 00630 if(vb >= 5) Cout("IntfReadImgDynCoeffFile" << endl); 00631 00632 // Init Interfile Key structure 00633 Intf_fields Img_fields; 00634 IntfKeyInitFields(&Img_fields); 00635 00636 // List containing single/multiple interfile headers of the images to read 00637 vector<string> lpath_to_headers; 00638 00639 // First check if the provided file is a metaheader with multiple files 00640 // and recover the list of image header file 00641 if(IntfIsMHD(a_pathToHeaderFile, lpath_to_headers) <0) // error 00642 { 00643 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00644 return 1; 00645 } 00646 00647 00648 00649 // --- Case 1 :We have one single data file containing all the data --- // 00650 if(lpath_to_headers.size() == 1) 00651 { 00652 // Recover image infos from either metaheader or single header file 00653 if(IntfReadHeader(a_pathToHeaderFile, &Img_fields, vb) ) 00654 { 00655 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00656 return 1; 00657 } 00658 00659 // Check consistency between image interfile data and the algorithm data 00660 if(IntfCheckConsistency(&Img_fields, ap_ID, vb, a_lerpFlag) ) 00661 { 00662 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while checking consistencies between the reconstruction parameters and the interfile keys in the header '"<< a_pathToHeaderFile << "' !" << endl); 00663 return 1; 00664 } 00665 00666 // If interpolation required, instanciate matrix with original image dimensions. 00667 // Todo : find a way to free memory if error occurs in following functions 00668 FLTNB* pimg_erp = NULL; 00669 IntfAllocInterpImg(&pimg_erp, Img_fields); 00670 00671 ifstream img_file(Img_fields.path_to_image.c_str(), ios::binary | ios::in); 00672 00673 if(img_file) 00674 { 00675 // Get the position of the beginning of the image data 00676 uint32_t offset = Img_fields.data_offset; 00677 00678 // Call the right IntfReadData() function according to the pixel type, and read data 00679 // offset is updated with the called functions within the loops 00680 for(int bf=0 ; bf<a_nbFbases ; bf++) 00681 IntfGetPixelTypeAndReadData(Img_fields, &img_file, a2p_ImgMatrix[bf], pimg_erp, &offset, ap_ID->GetNbVoxXYZ(), vb); 00682 } 00683 else 00684 { 00685 Cerr("***** IntfReadImgDynCoeffFile()-> Error occurred while trying to read the image file at the path: '"<< Img_fields.path_to_image << "' !" << endl); 00686 return 1; 00687 } 00688 00689 // If interpolation was required, deallocate image matrix with original image dimensions 00690 if(pimg_erp) delete[] pimg_erp; 00691 } 00692 00693 00694 00695 00696 // --- Case 2 : We have a number of data file equal to the number of image to load --- // 00697 else 00698 { 00699 // Recover image infos from the metaheader 00700 if(IntfReadHeader(a_pathToHeaderFile, &Img_fields, vb) ) 00701 { 00702 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while trying to read the interfile metaheader '"<< a_pathToHeaderFile << "' !" << endl); 00703 return 1; 00704 } 00705 00706 00707 if(lpath_to_headers.size() != (uint32_t)a_nbFbases) // Check we have a number of file corresponding to the number of images to load 00708 { 00709 Cerr("***** IntfReadImgDynCoeffFile()-> Error : number of interfile images ("<< lpath_to_headers.size() << 00710 ") not consistent with the number of parametric images (" << a_nbFbases<< ") !" << endl); 00711 return 1; 00712 } 00713 00714 // If interpolation required, instanciate matrix with original image dimensions. 00715 // Todo : find a way to free memory if error occurs in following functions 00716 FLTNB* pimg_erp = NULL; 00717 IntfAllocInterpImg(&pimg_erp, Img_fields); 00718 00719 for(int bf=0 ; bf<a_nbFbases ; bf++) 00720 { 00721 // Recover image infos from the specific interfile header 00722 if(IntfReadHeader(lpath_to_headers[bf], &Img_fields, vb) ) 00723 { 00724 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while trying to read the interfile header '"<< a_pathToHeaderFile << "' !" << endl); 00725 return 1; 00726 } 00727 00728 // Check consistency between image interfile data and the algorithm data 00729 if(IntfCheckConsistency(&Img_fields, ap_ID, vb, a_lerpFlag) ) 00730 { 00731 Cerr("***** IntfReadImgDynCoeffFile()-> Error : while checking consistencies between the reconstruction parameters and the interfile keys in the header '"<< a_pathToHeaderFile << "' !" << endl); 00732 return 1; 00733 } 00734 00735 ifstream img_file(Img_fields.path_to_image.c_str(), ios::binary | ios::in); 00736 if(img_file) 00737 { 00738 // Get the position of the beginning of the image data (assuming that offset is the same for each file (if any) ) 00739 uint32_t offset = Img_fields.data_offset; 00740 00741 // Call the right IntfReadData() function according to the pixel type, and read data (offset not updated) 00742 IntfGetPixelTypeAndReadData(Img_fields, &img_file, a2p_ImgMatrix[bf], pimg_erp, &offset, ap_ID->GetNbVoxXYZ(), vb); 00743 } 00744 else 00745 { 00746 Cerr("***** IntfReadImgDynCoeffFile()-> Error occurred while trying to read the image file at the path: '"<< Img_fields.path_to_image << "' !" << endl); 00747 return 1; 00748 } 00749 } 00750 00751 // If interpolation was required, deallocate image matrix with original image dimensions 00752 if(pimg_erp) delete[] pimg_erp; 00753 } 00754 00755 return 0; 00756 } 00757 00758 00759 00760 00761 // ===================================================================== 00762 // --------------------------------------------------------------------- 00763 // --------------------------------------------------------------------- 00764 // ===================================================================== 00765 /* 00766 \fn IntfAllocInterpImg 00767 \param a2p_img : Pointer to 1 dimensional image matrix to recover the image to interpolate 00768 \param ap_IF : Structure containing the interfile fields read in a interfile header 00769 \brief Allocate memory for an image matrix to recover an image to interpolate 00770 */ 00771 void IntfAllocInterpImg(FLTNB **a2p_img, Intf_fields a_IF) 00772 { 00773 if(a_IF.is_mtx_size_different == true) 00774 { 00775 uint32_t nb_vox = a_IF.mtx_size[0] * 00776 a_IF.mtx_size[1] * 00777 a_IF.mtx_size[2]; 00778 00779 // Allocate image and 00780 *a2p_img = new FLTNB[ nb_vox ]; 00781 ::memset(*a2p_img, 0, sizeof(FLTNB) * nb_vox); 00782 } 00783 } 00784 00785 00786 00787 00788 // ===================================================================== 00789 // --------------------------------------------------------------------- 00790 // --------------------------------------------------------------------- 00791 // ===================================================================== 00792 /* 00793 \fn IntfCheckConsistency 00794 \param ap_IF : Structure containing the interfile fields read in a interfile header 00795 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 00796 \param vb : Verbosity level 00797 \param a_lerpFlag : if true, enable linear interpolation of the image if img dimensions differ from the reconstruction dimensions 00798 \brief Check if the mandatory fields have been initialize in the ap_IF structure, and check consistencies with the reconstruction parameters 00799 \details This function also checks if the matrix size of the original image is different to the reconstruction matrix size. 00800 In this case a boolean is set up to enable interpolation during image reading 00801 \todo Add check for all mandatory parameters, and temporal image dimensions 00802 \todo Float comparison ? 00803 \return 0 if success, positive value otherwise. 00804 */ 00805 int IntfCheckConsistency(Intf_fields* ap_IF, oImageDimensionsAndQuantification* ap_ID, int vb, int a_lerpFlag) 00806 { 00807 if(vb >= 5) Cout("IntfCheckConsistency()" << endl); 00808 00809 // Check if main keys have been initialized 00810 if(ap_IF->path_to_image.size()<1 || 00811 ap_IF->nb_format == "" || 00812 ap_IF->mtx_size[0]==0 || ap_IF->mtx_size[1]==0 || ap_IF->mtx_size[2]==0 ) 00813 { 00814 Cerr("***** IntfCheckConsistency()-> Error : some mandatory keys not initialized. Cannot read the interfile image !" << endl); 00815 // Print the related key 00816 if(ap_IF->path_to_image.size()<1) 00817 Cerr(" Error when trying to read path to image data" << endl); 00818 if(ap_IF->nb_format == "") 00819 Cerr(" Error when trying to read data voxel type " << endl); 00820 if(ap_IF->mtx_size[0]==0 || ap_IF->mtx_size[1]==0 || ap_IF->mtx_size[2]==0 ) 00821 Cerr(" Error when trying to read matrix size (image dimensions) : x= " << ap_IF->mtx_size[0] << ", y= " << ap_IF->mtx_size[1] << ", z= " << ap_IF->mtx_size[2]<< endl); 00822 return 1; 00823 } 00824 00825 // If voxel size not found, throw warning 00826 if( ap_IF->vox_size[0]<0 || ap_IF->vox_size[1]<0 || ap_IF->vox_size[2]<0) 00827 { 00828 if(vb == 5) 00829 Cerr("***** IntfCheckConsistency()-> WARNING : No information found about voxel size ('scaling factor (mm/pixel)' tags.). Missing voxel sizes will be set to 1mm !" << endl); 00830 00831 if(ap_IF->vox_size[0]<0) ap_IF->vox_size[0] = 1.; 00832 if(ap_IF->vox_size[1]<0) ap_IF->vox_size[1] = 1.; 00833 if(ap_IF->vox_size[2]<0) ap_IF->vox_size[2] = 1.; 00834 if(vb == 5) 00835 Cerr(" Voxel sizes : x= " << ap_IF->vox_size[0] << ", y= " << ap_IF->vox_size[1] << ", z= " << ap_IF->vox_size[2]<< endl); 00836 } 00837 00838 00839 // If original dimensions don't match reconstructions dimensions/voxel sizes, 00840 // recover this data in Intf_fields object (if a_lerpFlag==true) or throw error (if a_lerpFlag==false) 00841 if( ((INTNB)(ap_IF->mtx_size[0])) != ap_ID->GetNbVoxX() || 00842 ((INTNB)(ap_IF->mtx_size[1])) != ap_ID->GetNbVoxY() || 00843 ((INTNB)(ap_IF->mtx_size[2])) != ap_ID->GetNbVoxZ() || 00844 !FLTNBIsEqual(ap_IF->vox_size[0], ap_ID->GetVoxSizeX(), pow(0.1, std::numeric_limits<FLTNB>::digits10+1) ) || 00845 !FLTNBIsEqual(ap_IF->vox_size[1], ap_ID->GetVoxSizeY(), pow(0.1, std::numeric_limits<FLTNB>::digits10+1) ) || 00846 !FLTNBIsEqual(ap_IF->vox_size[2], ap_ID->GetVoxSizeZ(), pow(0.1, std::numeric_limits<FLTNB>::digits10+1) ) ) 00847 { 00848 if(a_lerpFlag) 00849 { 00850 ap_IF->cmtx_size[0] = (uint32_t)(ap_ID->GetNbVoxX()); 00851 ap_IF->cmtx_size[1] = (uint32_t)(ap_ID->GetNbVoxY()); 00852 ap_IF->cmtx_size[2] = (uint32_t)(ap_ID->GetNbVoxZ()); 00853 ap_IF->cvox_size[0] = ap_ID->GetVoxSizeX(); 00854 ap_IF->cvox_size[1] = ap_ID->GetVoxSizeY(); 00855 ap_IF->cvox_size[2] = ap_ID->GetVoxSizeZ(); 00856 ap_IF->is_mtx_size_different = true; // Set this boolean to true to indicate that an interpolation step will be required 00857 } 00858 else 00859 { 00860 Cerr("***** IntfCheckConsistency()-> Error : Image dimensions don't match reconstructions dimensions/voxel sizes"); 00861 Cerr(" and linear interpolation is disabled (a_lerpFlag is false) !" << endl); 00862 Cerr("***** Recovered image dimensions (x;y;z): "<< ap_IF->mtx_size[0] <<" ; "<< ap_IF->mtx_size[1] << " ; " << ap_IF->mtx_size[2] << endl); 00863 Cerr("***** Reconstruction dimensions (x;y;z) : "<< ap_ID->GetNbVoxX() <<" ; "<< ap_ID->GetNbVoxY() << " ; " << ap_ID->GetNbVoxZ() << endl); 00864 Cerr("***** Image voxel sizes (x;y;z) : "<< ap_IF->vox_size[0] <<" ; "<< ap_IF->vox_size[1] << " ; " << ap_IF->vox_size[2] << endl); 00865 Cerr("***** Reconstruction voxel sizes (x;y;z): "<< ap_ID->GetVoxSizeX() <<" ; "<< ap_ID->GetVoxSizeY() << " ; " << ap_ID->GetVoxSizeZ() << endl); 00866 return 1; 00867 } 00868 } 00869 00870 return 0; 00871 } 00872 00873 00874 00875 00876 // ===================================================================== 00877 // --------------------------------------------------------------------- 00878 // --------------------------------------------------------------------- 00879 // ===================================================================== 00880 /* 00881 \fn IntfGetPixelTypeAndReadData 00882 \param a_IF : Interfile fields recovered from the header 00883 \param ap_iFile : Ifstream pointing to an image file 00884 \param ap_outImgMtx : 3D image matrix with reconstruction dimensions/voxel size 00885 \param ap_inImgMtx : 3D image matrix with original dimensions/voxel size 00886 \param ap_offset : Offset indicating the beginning of the data to read in the image file 00887 \param a_nbVox : A number of voxels in the 3D image matrix with reconstruction dimensions/voxel size 00888 \param vb : Verbosity level 00889 \brief The purpose of this function is to call the templated ReadData() function with the data type corresponding to the interfile image 00890 \details It uses "number format" and "number of bytes per pixel" fields to identify the correct type 00891 ASCII and bit images NOT supported 00892 \return 0 if success, positive value otherwise. 00893 */ 00894 int IntfGetPixelTypeAndReadData(Intf_fields a_IF, 00895 ifstream* ap_iFile, 00896 FLTNB* ap_outImgMatrix, 00897 FLTNB* ap_inImgMatrix, 00898 uint32_t* ap_offset, 00899 int a_nbVox, 00900 int vb) 00901 { 00902 if(vb >= 5) Cout("IntfGetPixelTypeAndReadData() " << endl); 00903 00904 // To check if an error occurred in one of the called functions 00905 int error = 0; 00906 00907 // Check the image data voxel size according to the Interfile fields 'number_format' and 'number of bytes per pixel' 00908 if (a_IF.nb_format == FLT32_str) // Float data 00909 { 00910 float flt; 00911 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &flt); 00912 } 00913 else if (a_IF.nb_format == FLT64_str) // Double data 00914 { 00915 double db; 00916 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &db); 00917 } 00918 else if (a_IF.nb_format == INT32_str) // Signed integer data 00919 { 00920 if(a_IF.nb_bytes_pixel == 1) 00921 { 00922 int8_t s_int; 00923 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &s_int); 00924 } 00925 else if(a_IF.nb_bytes_pixel == 2) 00926 { 00927 int16_t s_int; 00928 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &s_int); 00929 } 00930 else if(a_IF.nb_bytes_pixel == 8) 00931 { 00932 int64_t s_int; 00933 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &s_int); 00934 } 00935 else // default : 4 bytes 00936 { 00937 int32_t s_int; 00938 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &s_int); 00939 } 00940 } 00941 // ASCII and bit images TODO ? 00942 00943 //else if (a_IF.nb_format == ASCII_str) 00944 //{ 00945 // // Data in ASCII 00946 // uint32_t ascii; 00947 // IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &ascii); 00948 //} 00949 //else if(a_IF.nb_format == BIT_str) 00950 //{ 00951 // // Data in bits 00952 // uint32_t bit; 00953 // IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &bit); 00954 //} 00955 else // Unsigned integer data (default) 00956 { 00957 if(a_IF.nb_bytes_pixel == 1) 00958 { 00959 uint8_t u_int; 00960 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &u_int); 00961 } 00962 else if(a_IF.nb_bytes_pixel == 2) 00963 { 00964 uint16_t u_int; 00965 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &u_int); 00966 } 00967 else if(a_IF.nb_bytes_pixel == 8) 00968 { 00969 uint64_t u_int; 00970 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &u_int); 00971 } 00972 else // default : 4 bytes 00973 { 00974 uint32_t u_int; 00975 error = IntfReadData(a_IF, ap_iFile, ap_outImgMatrix, ap_inImgMatrix, ap_offset, a_nbVox, vb, &u_int); 00976 } 00977 } 00978 00979 if(error == 1) 00980 { 00981 Cerr("***** IntfGetPixelTypeAndReadData-> Error occurred when trying to read an interfile image !" << endl); 00982 return 1; 00983 } 00984 00985 return 0; 00986 } 00987 00988 00989 00990 00991 // ===================================================================== 00992 // --------------------------------------------------------------------- 00993 // --------------------------------------------------------------------- 00994 // ===================================================================== 00995 /* 00996 \fn IntfReadData 00997 \param a_IF : Interfile fields recovered from the header 00998 \param ap_iFile : Ifstream pointing to an image file 00999 \param ap_outImgMtx : 3D image matrix with reconstruction dimensions/voxel size 01000 \param ap_inImgMtx : 3D image matrix with original dimensions/voxel size 01001 \param ap_offset : Offset indicating the beginning of the data to read in the image file 01002 \param a_nbVox : A number of voxels in the 3D image matrix with reconstruction dimensions/voxel size 01003 \param vb : Verbosity level 01004 \param T* bytes : Buffer of templated size, to recover any voxel value 01005 \brief Templated function which read an image voxel by voxel and store it in the ap_outImgMtx image matrix 01006 \details Call an interpolation function if the original image dimensions/voxel sizes are different to the reconstruction dimensions/voxel sizes 01007 Manage endianness and the optionnal calibration with the rescale slope/intercept interfile keys 01008 \todo Re-orient image data from the original orientation to the default orientation (transaxial/supine/headin) 01009 \todo Perhaps check the conversion from original image type to FLTNB type has been correctly performed (maybe time-consuming) 01010 \return 0 if success, positive value otherwise. 01011 */ 01012 template <class T> 01013 int IntfReadData(Intf_fields a_IF, 01014 ifstream* ap_iFile, 01015 FLTNB* ap_outImgMatrix, 01016 FLTNB* ap_inImgMatrix, 01017 uint32_t* a_offset, 01018 int a_nbVox, 01019 int vb, 01020 T* bytes) 01021 { 01022 if(vb >= 5) Cout("IntfReadData() " << endl); 01023 01024 int nb_vox = 0; // number of voxels in the matrix which will recover the image data 01025 FLTNB* pimg = NULL; // pointer to the matrix which will recover the image data 01026 01027 // Check if interpolation is required 01028 if(a_IF.is_mtx_size_different) 01029 { 01030 // Set the number of voxels in the image to the original image dimensions as they differ from reconstruction dimensions 01031 nb_vox = a_IF.mtx_size[0]*a_IF.mtx_size[1]*a_IF.mtx_size[2]; 01032 pimg = ap_inImgMatrix; // recover the image into the image matrix initialized with interfile dimensions 01033 } 01034 else 01035 { 01036 nb_vox = a_nbVox; 01037 pimg = ap_outImgMatrix; // directly recover the image into the CASToR image matrix 01038 } 01039 01040 /* ORIGINAL CODE 01041 // Loop on image voxels 01042 for(int v=0 ; v<nb_vox ; v++) 01043 { 01044 // Todo Re-orient image data from the original orientation to the default orientation (transaxial/supine/headin) 01045 //int voxel_idx = IntfGetVoxIdxSHTOrientation(a_IF, v); 01046 int voxel_idx = v; // For now no orientations 01047 01048 // Set the steam position to the data to read and store it in the buffer 01049 ap_iFile->seekg(*a_offset); 01050 ap_iFile->read((char*)bytes, sizeof(T)); 01051 01052 // Swap data if endianness is not the same as user processor 01053 if(a_IF.endianness != User_Endianness) SwapBytes(bytes); 01054 01055 // Cast to CASToR image format 01056 FLTNB buffer = (FLTNB)*bytes; // TODO Check the casting is ok ? May slow down image reading 01057 01058 // Calibrate data using rescale slope and intercept if needed. 01059 // Then write in the image matrix 01060 pimg[voxel_idx] = buffer*a_IF.rescale_slope + a_IF.rescale_intercept; 01061 01062 // set the offset to the next position 01063 *a_offset += sizeof(T); 01064 } 01065 * */ 01066 01067 // BEGINNING OF NEW CODE (we read in a block locally allocated and freed at the end, this is XXXXXX times faster) 01068 01069 // Allocate temporary reading buffer 01070 bytes = (T*)malloc(nb_vox*sizeof(T)); 01071 // Seek to the provided position 01072 ap_iFile->seekg(*a_offset); 01073 // Read the data 01074 ap_iFile->read((char*)bytes, nb_vox*sizeof(T)); 01075 01076 // Loop over all voxels 01077 for (int v=0; v<nb_vox; v++) 01078 { 01079 // Todo Re-orient image data from the original orientation to the default orientation (transaxial/supine/headin) 01080 //int voxel_idx = IntfGetVoxIdxSHTOrientation(a_IF, v); 01081 int voxel_idx = v; // For now no orientations 01082 01083 // Swap data if endianness is not the same as user processor 01084 if (a_IF.endianness != User_Endianness) SwapBytes(&(bytes[v])); 01085 // Cast to CASToR image format 01086 FLTNB buffer = (FLTNB)(bytes[v]); 01087 // Calibrate data using rescale slope and intercept if needed. 01088 // Then write in the image matrix 01089 pimg[voxel_idx] = buffer*a_IF.rescale_slope + a_IF.rescale_intercept; 01090 } 01091 01092 free(bytes); 01093 01094 // END OF NEW CODE 01095 01096 // Call interpolation function if required 01097 if(a_IF.is_mtx_size_different) 01098 { 01099 if(ImageInterpolation(ap_inImgMatrix, ap_outImgMatrix, 01100 a_IF.mtx_size, a_IF.cmtx_size, 01101 a_IF.vox_size, a_IF.cvox_size) ) 01102 { 01103 Cerr("***** IntfReadData-> Error occurred when interpolating the interfile image to the reconstruction dimensions !" << endl); 01104 return 1; 01105 } 01106 } 01107 01108 return 0; 01109 } 01110 01111 01112 01113 01114 // ===================================================================== 01115 // --------------------------------------------------------------------- 01116 // --------------------------------------------------------------------- 01117 // ===================================================================== 01118 /* 01119 \fn ImageInterpolation 01120 \param ap_iImg : Image matrix to interpolate, with dimensions of the original interfile 01121 \param ap_oImg : Image matrix to recover the interpolated image to the reconstruction dimensions 01122 \param ap_iDimVox[3] : X,Y,Z dimensions of the image to interpolate 01123 \param ap_oDimVox[3] : X,Y,Z dimensions for the reconstruction 01124 \param ap_iSizeVox[3] : X,Y,Z voxel size of the image to interpolate 01125 \param ap_oSizeVox[3] : X,Y,Z voxel size for the reconstruction 01126 \brief Trilinear interpolation 01127 \todo Manage image position (offsets) 01128 \return 0 if success, positive value otherwise. 01129 */ 01130 int ImageInterpolation(FLTNB *ap_iImg, FLTNB *ap_oImg, 01131 const uint32_t ap_iDimVox[3], const uint32_t ap_oDimVox[3], 01132 const FLTNB ap_iSizeVox[3], const FLTNB ap_oSizeVox[3]) 01133 { 01134 // Assuming the two images are registered 01135 // todo : deal with any potential offset if the input data contains following keys 01136 // "origin (mm) [x], offset [x], first pixel offset (mm) [x] " 01137 01138 /* ORIGINAL FUNCTIONS IMPLEMENTED ONLY WITH FLOATS 01139 float const posOldImage[] = {-(float)ap_iDimVox[0]*ap_iSizeVox[0]/2 , 01140 -(float)ap_iDimVox[1]*ap_iSizeVox[1]/2 , 01141 -(float)ap_iDimVox[2]*ap_iSizeVox[2]/2 }; 01142 01143 float const posNewImage[] = {-(float)ap_oDimVox[0]*ap_oSizeVox[0]/2 , 01144 -(float)ap_oDimVox[1]*ap_oSizeVox[1]/2 , 01145 -(float)ap_oDimVox[2]*ap_oSizeVox[2]/2 }; 01146 01147 // Padding the old image in each direction 01148 uint32_t const iDimVoxPad[] 01149 { 01150 ap_iDimVox[ 0 ] + 2, 01151 ap_iDimVox[ 1 ] + 2, 01152 ap_iDimVox[ 2 ] + 2 01153 }; 01154 01155 uint32_t const nElts = iDimVoxPad[ 0 ] * 01156 iDimVoxPad[ 1 ] * 01157 iDimVoxPad[ 2 ]; 01158 01159 // Allocating a new buffer storing the padded image 01160 float *pPadIData = new float[ nElts ]; 01161 ::memset( pPadIData, 0, sizeof( float ) * nElts ); 01162 01163 for( uint32_t k = 0; k < ap_iDimVox[ 2 ]; ++k ) 01164 for( uint32_t j = 0; j < ap_iDimVox[ 1 ]; ++j ) 01165 for( uint32_t i = 0; i < ap_iDimVox[ 0 ]; ++i ) 01166 { 01167 pPadIData[ ( i + 1 ) + ( j + 1 ) * iDimVoxPad[ 0 ] 01168 + ( k + 1 ) * iDimVoxPad[ 0 ] * iDimVoxPad[ 1 ] ] = 01169 ap_iImg[ i + j * ap_iDimVox[ 0 ] 01170 + k * ap_iDimVox[ 0 ] * ap_iDimVox[ 1 ] ]; 01171 } 01172 01173 01174 // Computing the bounds in each direction depending on the pad 01175 float const boundMin[] 01176 { 01177 posOldImage[ 0 ] - ap_iSizeVox[ 0 ] / 2.0f, 01178 posOldImage[ 1 ] - ap_iSizeVox[ 1 ] / 2.0f, 01179 posOldImage[ 2 ] - ap_iSizeVox[ 2 ] / 2.0f 01180 }; 01181 01182 float const boundMax[] 01183 { 01184 posOldImage[ 0 ] + (float)ap_iDimVox[ 0 ] * ap_iSizeVox[ 0 ] 01185 + ap_iSizeVox[ 0 ] / 2.0f, 01186 posOldImage[ 1 ] + (float)ap_iDimVox[ 1 ] * ap_iSizeVox[ 1 ] 01187 + ap_iSizeVox[ 1 ] / 2.0f, 01188 posOldImage[ 2 ] + (float)ap_iDimVox[ 2 ] * ap_iSizeVox[ 2 ] 01189 + ap_iSizeVox[ 2 ] / 2.0f 01190 }; 01191 01192 // Computing and storing the position of the center of the pixel in each 01193 // direction for the initial image 01194 float **pOldCoordCenter = new float*[ 3 ]; 01195 01196 for( uint32_t i = 0; i < 3; ++i ) 01197 { 01198 pOldCoordCenter[ i ] = new float[ iDimVoxPad[ i ] ]; 01199 // Set the values 01200 for( uint32_t j = 0; j < iDimVoxPad[ i ]; ++j ) 01201 { 01202 pOldCoordCenter[ i ][ j ] = posOldImage[ i ] - ap_iSizeVox[ i ] / 2.0f 01203 + j * ap_iSizeVox[ i ]; 01204 } 01205 } 01206 01207 // Computing and storing the position of the center of the pixel in each 01208 // direction of the resampled image 01209 float **pNewCoordCenter = new float*[ 3 ]; 01210 for( uint32_t i = 0; i < 3; ++i ) 01211 { 01212 pNewCoordCenter[ i ] = new float[ ap_oDimVox[ i ] ]; 01213 // Set the values 01214 for( uint32_t j = 0; j < ap_oDimVox[ i ]; ++j ) 01215 { 01216 pNewCoordCenter[ i ][ j ] = posNewImage[ i ] + ap_oSizeVox[ i ] / 2.0f 01217 + j * ap_oSizeVox[ i ]; 01218 } 01219 } 01220 01221 // Defining parameters 01222 float const invSizeX = 1.0f / ap_iSizeVox[ 0 ]; 01223 float const invSizeY = 1.0f / ap_iSizeVox[ 1 ]; 01224 float const invSizeZ = 1.0f / ap_iSizeVox[ 2 ]; 01225 01226 // Allocating memory for the 8 pixels kernel 01227 float *pKernelData = new float[ 8 ]; 01228 01229 // Loop over the elements of the new images 01230 for( uint32_t k = 0; k < ap_oDimVox[ 2 ]; ++k ) 01231 { 01232 // Get the coordinate in Z 01233 float const z = pNewCoordCenter[ 2 ][ k ]; 01234 if( z < boundMin[ 2 ] || z > boundMax[ 2 ] ) continue; 01235 01236 // Find the bin in the old image 01237 int32_t const zBin = ( z - boundMin[ 2 ] ) * invSizeZ; 01238 01239 // Computing the 2 z composants 01240 float const zComposantI0 = invSizeZ * ( pOldCoordCenter[ 2 ][ zBin + 1 ] - z ); 01241 float const zComposantI1 = invSizeZ * ( z - pOldCoordCenter[ 2 ][ zBin ] ); 01242 01243 for( uint32_t j = 0; j < ap_oDimVox[ 1 ]; ++j ) 01244 { 01245 // Get the coordinate in Y 01246 float const y = pNewCoordCenter[ 1 ][ j ]; 01247 if( y < boundMin[ 1 ] || y > boundMax[ 1 ] ) continue; 01248 01249 // Find the bin in the old image 01250 int32_t const yBin = ( y - boundMin[ 1 ] ) * invSizeY; 01251 01252 // Computing the 2 y composants 01253 float const yComposantI0 = invSizeY * ( pOldCoordCenter[ 1 ][ yBin + 1 ] 01254 - y ); 01255 float const yComposantI1 = invSizeY * ( y 01256 - pOldCoordCenter[ 1 ][ yBin ] ); 01257 01258 for( uint32_t i = 0; i < ap_oDimVox[ 0 ]; ++i ) 01259 { 01260 // Get the coordinate in X 01261 float const x = pNewCoordCenter[ 0 ][ i ]; 01262 if( x < boundMin[ 0 ] || x > boundMax[ 0 ] ) continue; 01263 01264 // Find the bin in the old image 01265 int32_t const xBin = ( x - boundMin[ 0 ] ) * invSizeX; 01266 01267 // Computing the 2 x composants 01268 float const xComposantI0 = invSizeX * ( 01269 pOldCoordCenter[ 0 ][ xBin + 1 ] - x ); 01270 float const xComposantI1 = invSizeX * ( x 01271 - pOldCoordCenter[ 0 ][ xBin ] ); 01272 01273 // Fill the buffer storing the data 01274 for( uint32_t kk = 0; kk < 2; ++kk ) 01275 { 01276 for( uint32_t jj = 0; jj < 2; ++jj ) 01277 { 01278 for( uint32_t ii = 0; ii < 2; ++ii ) 01279 { 01280 pKernelData[ ii + jj * 2 + kk * 2 * 2 ] = 01281 pPadIData[ 01282 ( xBin + ii ) + 01283 ( yBin + jj ) * iDimVoxPad[ 0 ] + 01284 ( zBin + kk ) * iDimVoxPad[ 0 ] * iDimVoxPad[ 1 ] 01285 ]; 01286 } 01287 } 01288 } 01289 01290 // Computing the interpolation 01291 // In X 01292 float const xInterpVal0 = pKernelData[ 0 ] * xComposantI0 + 01293 pKernelData[ 1 ] * xComposantI1; 01294 01295 float const xInterpVal1 = pKernelData[ 2 ] * xComposantI0 + 01296 pKernelData[ 3 ] * xComposantI1; 01297 01298 float const xInterpVal2 = pKernelData[ 4 ] * xComposantI0 + 01299 pKernelData[ 5 ] * xComposantI1; 01300 01301 float const xInterpVal3 = pKernelData[ 6 ] * xComposantI0 + 01302 pKernelData[ 7 ] * xComposantI1; 01303 01304 // In Y 01305 float const yInterpVal0 = xInterpVal0 * yComposantI0 + 01306 xInterpVal1 * yComposantI1; 01307 01308 float const yInterpVal1 = xInterpVal2 * yComposantI0 + 01309 xInterpVal3 * yComposantI1; 01310 01311 // Final interpolation 01312 float const interpValTot = yInterpVal0 * zComposantI0 + 01313 yInterpVal1 * zComposantI1; 01314 01315 ap_oImg[ i + j * ap_oDimVox[ 0 ] 01316 + k * ap_oDimVox[ 0 ] * ap_oDimVox[ 1 ] ] = interpValTot; 01317 } 01318 } 01319 } 01320 01321 // Freeing the memory 01322 for( uint32_t i = 0; i < 3; ++i ) 01323 { 01324 delete[] pOldCoordCenter[ i ]; 01325 delete[] pNewCoordCenter[ i ]; 01326 } 01327 delete[] pOldCoordCenter; 01328 delete[] pNewCoordCenter; 01329 delete[] pPadIData; 01330 delete[] pKernelData; 01331 01332 return 0; 01333 */ 01334 01335 01336 FLTNB const posOldImage[] = {-((FLTNB)(ap_iDimVox[0]))*ap_iSizeVox[0]*((FLTNB)0.5) , 01337 -((FLTNB)(ap_iDimVox[1]))*ap_iSizeVox[1]*((FLTNB)0.5) , 01338 -((FLTNB)(ap_iDimVox[2]))*ap_iSizeVox[2]*((FLTNB)0.5) }; 01339 01340 FLTNB const posNewImage[] = {-((FLTNB)(ap_oDimVox[0]))*ap_oSizeVox[0]*((FLTNB)0.5) , 01341 -((FLTNB)(ap_oDimVox[1]))*ap_oSizeVox[1]*((FLTNB)0.5) , 01342 -((FLTNB)(ap_oDimVox[2]))*ap_oSizeVox[2]*((FLTNB)0.5) }; 01343 01344 // Padding the old image in each direction 01345 uint32_t const iDimVoxPad[] 01346 { 01347 ap_iDimVox[ 0 ] + 2, 01348 ap_iDimVox[ 1 ] + 2, 01349 ap_iDimVox[ 2 ] + 2 01350 }; 01351 01352 uint32_t const nElts = iDimVoxPad[ 0 ] * 01353 iDimVoxPad[ 1 ] * 01354 iDimVoxPad[ 2 ]; 01355 01356 // Allocating a new buffer storing the padded image 01357 FLTNB *pPadIData = new FLTNB[ nElts ]; 01358 ::memset( pPadIData, 0, sizeof( FLTNB ) * nElts ); 01359 01360 for( uint32_t k = 0; k < ap_iDimVox[ 2 ]; ++k ) 01361 for( uint32_t j = 0; j < ap_iDimVox[ 1 ]; ++j ) 01362 for( uint32_t i = 0; i < ap_iDimVox[ 0 ]; ++i ) 01363 { 01364 pPadIData[ ( i + 1 ) + ( j + 1 ) * iDimVoxPad[ 0 ] 01365 + ( k + 1 ) * iDimVoxPad[ 0 ] * iDimVoxPad[ 1 ] ] = 01366 ap_iImg[ i + j * ap_iDimVox[ 0 ] 01367 + k * ap_iDimVox[ 0 ] * ap_iDimVox[ 1 ] ]; 01368 } 01369 01370 01371 // Computing the bounds in each direction depending on the pad 01372 FLTNB const boundMin[] 01373 { 01374 posOldImage[ 0 ] - ap_iSizeVox[ 0 ] * ((FLTNB)0.5), 01375 posOldImage[ 1 ] - ap_iSizeVox[ 1 ] * ((FLTNB)0.5), 01376 posOldImage[ 2 ] - ap_iSizeVox[ 2 ] * ((FLTNB)0.5) 01377 }; 01378 01379 FLTNB const boundMax[] 01380 { 01381 posOldImage[ 0 ] + ((FLTNB)ap_iDimVox[ 0 ]) * ap_iSizeVox[ 0 ] 01382 + ap_iSizeVox[ 0 ] * ((FLTNB)0.5), 01383 posOldImage[ 1 ] + ((FLTNB)ap_iDimVox[ 1 ]) * ap_iSizeVox[ 1 ] 01384 + ap_iSizeVox[ 1 ] * ((FLTNB)0.5), 01385 posOldImage[ 2 ] + ((FLTNB)ap_iDimVox[ 2 ]) * ap_iSizeVox[ 2 ] 01386 + ap_iSizeVox[ 2 ] * ((FLTNB)0.5) 01387 }; 01388 01389 // Computing and storing the position of the center of the pixel in each 01390 // direction for the initial image 01391 FLTNB **pOldCoordCenter = new FLTNB*[ 3 ]; 01392 01393 for( uint32_t i = 0; i < 3; ++i ) 01394 { 01395 pOldCoordCenter[ i ] = new FLTNB[ iDimVoxPad[ i ] ]; 01396 // Set the values 01397 for( uint32_t j = 0; j < iDimVoxPad[ i ]; ++j ) 01398 { 01399 pOldCoordCenter[ i ][ j ] = posOldImage[ i ] - ap_iSizeVox[ i ] / 2.0 01400 + j * ap_iSizeVox[ i ]; 01401 } 01402 } 01403 01404 // Computing and storing the position of the center of the pixel in each 01405 // direction of the resampled image 01406 FLTNB **pNewCoordCenter = new FLTNB*[ 3 ]; 01407 for( uint32_t i = 0; i < 3; ++i ) 01408 { 01409 pNewCoordCenter[ i ] = new FLTNB[ ap_oDimVox[ i ] ]; 01410 // Set the values 01411 for( uint32_t j = 0; j < ap_oDimVox[ i ]; ++j ) 01412 { 01413 pNewCoordCenter[ i ][ j ] = posNewImage[ i ] + ap_oSizeVox[ i ] / 2.0 01414 + j * ap_oSizeVox[ i ]; 01415 } 01416 } 01417 01418 // Defining parameters 01419 FLTNB const invSizeX = 1.0 / ap_iSizeVox[ 0 ]; 01420 FLTNB const invSizeY = 1.0 / ap_iSizeVox[ 1 ]; 01421 FLTNB const invSizeZ = 1.0 / ap_iSizeVox[ 2 ]; 01422 01423 // Allocating memory for the 8 pixels kernel 01424 FLTNB *pKernelData = new FLTNB[ 8 ]; 01425 01426 // Loop over the elements of the new images 01427 for( uint32_t k = 0; k < ap_oDimVox[ 2 ]; ++k ) 01428 { 01429 // Get the coordinate in Z 01430 FLTNB const z = pNewCoordCenter[ 2 ][ k ]; 01431 if( z < boundMin[ 2 ] || z > boundMax[ 2 ] ) continue; 01432 01433 // Find the bin in the old image 01434 int32_t const zBin = ( z - boundMin[ 2 ] ) * invSizeZ; 01435 01436 // Computing the 2 z composants 01437 FLTNB const zComposantI0 = invSizeZ * ( pOldCoordCenter[ 2 ][ zBin + 1 ] - z ); 01438 FLTNB const zComposantI1 = invSizeZ * ( z - pOldCoordCenter[ 2 ][ zBin ] ); 01439 01440 for( uint32_t j = 0; j < ap_oDimVox[ 1 ]; ++j ) 01441 { 01442 // Get the coordinate in Y 01443 FLTNB const y = pNewCoordCenter[ 1 ][ j ]; 01444 if( y < boundMin[ 1 ] || y > boundMax[ 1 ] ) continue; 01445 01446 // Find the bin in the old image 01447 int32_t const yBin = ( y - boundMin[ 1 ] ) * invSizeY; 01448 01449 // Computing the 2 y composants 01450 FLTNB const yComposantI0 = invSizeY * ( pOldCoordCenter[ 1 ][ yBin + 1 ] 01451 - y ); 01452 FLTNB const yComposantI1 = invSizeY * ( y 01453 - pOldCoordCenter[ 1 ][ yBin ] ); 01454 01455 for( uint32_t i = 0; i < ap_oDimVox[ 0 ]; ++i ) 01456 { 01457 // Get the coordinate in X 01458 FLTNB const x = pNewCoordCenter[ 0 ][ i ]; 01459 if( x < boundMin[ 0 ] || x > boundMax[ 0 ] ) continue; 01460 01461 // Find the bin in the old image 01462 int32_t const xBin = ( x - boundMin[ 0 ] ) * invSizeX; 01463 01464 // Computing the 2 x composants 01465 FLTNB const xComposantI0 = invSizeX * ( 01466 pOldCoordCenter[ 0 ][ xBin + 1 ] - x ); 01467 FLTNB const xComposantI1 = invSizeX * ( x 01468 - pOldCoordCenter[ 0 ][ xBin ] ); 01469 01470 // Fill the buffer storing the data 01471 for( uint32_t kk = 0; kk < 2; ++kk ) 01472 { 01473 for( uint32_t jj = 0; jj < 2; ++jj ) 01474 { 01475 for( uint32_t ii = 0; ii < 2; ++ii ) 01476 { 01477 pKernelData[ ii + jj * 2 + kk * 2 * 2 ] = 01478 pPadIData[ 01479 ( xBin + ii ) + 01480 ( yBin + jj ) * iDimVoxPad[ 0 ] + 01481 ( zBin + kk ) * iDimVoxPad[ 0 ] * iDimVoxPad[ 1 ] 01482 ]; 01483 } 01484 } 01485 } 01486 01487 // Computing the interpolation 01488 // In X 01489 FLTNB const xInterpVal0 = pKernelData[ 0 ] * xComposantI0 + 01490 pKernelData[ 1 ] * xComposantI1; 01491 01492 FLTNB const xInterpVal1 = pKernelData[ 2 ] * xComposantI0 + 01493 pKernelData[ 3 ] * xComposantI1; 01494 01495 FLTNB const xInterpVal2 = pKernelData[ 4 ] * xComposantI0 + 01496 pKernelData[ 5 ] * xComposantI1; 01497 01498 FLTNB const xInterpVal3 = pKernelData[ 6 ] * xComposantI0 + 01499 pKernelData[ 7 ] * xComposantI1; 01500 01501 // In Y 01502 FLTNB const yInterpVal0 = xInterpVal0 * yComposantI0 + 01503 xInterpVal1 * yComposantI1; 01504 01505 FLTNB const yInterpVal1 = xInterpVal2 * yComposantI0 + 01506 xInterpVal3 * yComposantI1; 01507 01508 // Final interpolation 01509 FLTNB const interpValTot = yInterpVal0 * zComposantI0 + 01510 yInterpVal1 * zComposantI1; 01511 01512 ap_oImg[ i + j * ap_oDimVox[ 0 ] 01513 + k * ap_oDimVox[ 0 ] * ap_oDimVox[ 1 ] ] = interpValTot; 01514 } 01515 } 01516 } 01517 01518 // Freeing the memory 01519 for( uint32_t i = 0; i < 3; ++i ) 01520 { 01521 delete[] pOldCoordCenter[ i ]; 01522 delete[] pNewCoordCenter[ i ]; 01523 } 01524 delete[] pOldCoordCenter; 01525 delete[] pNewCoordCenter; 01526 delete[] pPadIData; 01527 delete[] pKernelData; 01528 01529 return 0; 01530 01531 01532 } 01533 01534 01535 01536 01537 // ===================================================================== 01538 // --------------------------------------------------------------------- 01539 // --------------------------------------------------------------------- 01540 // ===================================================================== 01541 /* 01542 \fn IntfWriteImgFile 01543 \param a_pathToImg : path to image basename 01544 \param ap_ImgMatrix : 1 dimensional image matrix which contains the image to write 01545 \param ap_IntfF : Intf_fields structure containing image metadata 01546 \param vb : verbosity 01547 \brief Main function dedicated to Interfile 3D image writing. \n 01548 Recover image information from a provided Intf_fields structure. 01549 \details Call the main functions dedicated to Interfile header and data writing : 01550 IntfWriteHeaderMainData() and then IntfWriteImage() 01551 \return 0 if success, positive value otherwise. 01552 */ 01553 int IntfWriteImgFile(const string& a_pathToImg, FLTNB* ap_ImgMatrix, const Intf_fields& ap_IntfF, int vb) 01554 { 01555 if (vb>=3) Cout("IntfWriteImgFile (with Intf_fields)" << endl); 01556 01557 // Write Interfile header 01558 if(IntfWriteHeaderMainData(a_pathToImg, ap_IntfF, vb) ) 01559 { 01560 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header for'"<< a_pathToImg << "' !" << endl); 01561 return 1; 01562 } 01563 01564 string path_to_image = a_pathToImg; 01565 path_to_image.append(".img"); 01566 01567 // Read Interfile image 01568 if(IntfWriteImage(path_to_image, ap_ImgMatrix, ap_IntfF.mtx_size[0] * ap_IntfF.mtx_size[1] * ap_IntfF.mtx_size[2], vb) ) 01569 { 01570 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile image '"<< path_to_image << "' !" << endl); 01571 return 1; 01572 } 01573 01574 return 0; 01575 } 01576 01577 01578 01579 01580 // ===================================================================== 01581 // --------------------------------------------------------------------- 01582 // --------------------------------------------------------------------- 01583 // ===================================================================== 01584 /* 01585 \fn IntfWriteImgFile 01586 \param a_pathToImg : path to image basename 01587 \param ap_ImgMatrix : 1 dimensional image matrix which contains the image to write 01588 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 01589 \param vb : verbosity 01590 \brief Main function dedicated to Interfile 3D image writing 01591 \details Call the main functions dedicated to Interfile header and data writing : 01592 IntfWriteHeaderMainData() and then IntfWriteImage() 01593 \todo Get metadata from a Intf_fields object ? 01594 (useful to transfer keys from read images to written images) 01595 \return 0 if success, positive value otherwise. 01596 */ 01597 int IntfWriteImgFile(const string& a_pathToImg, FLTNB* ap_ImgMatrix, oImageDimensionsAndQuantification* ap_ID, int vb) 01598 { 01599 if (vb>=3) Cout("IntfWriteImgFile (3D image)" << endl); 01600 01601 Intf_fields Img_fields; 01602 IntfKeySetFieldsOutput(&Img_fields, ap_ID); 01603 01604 // Write Interfile header 01605 if(IntfWriteHeaderMainData(a_pathToImg, Img_fields, vb) ) 01606 { 01607 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header for'"<< a_pathToImg << "' !" << endl); 01608 return 1; 01609 } 01610 01611 string path_to_image = a_pathToImg; 01612 path_to_image.append(".img"); 01613 01614 // Read Interfile image 01615 if(IntfWriteImage(path_to_image, ap_ImgMatrix, ap_ID->GetNbVoxXYZ(), vb) ) 01616 { 01617 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile image '"<< path_to_image << "' !" << endl); 01618 return 1; 01619 } 01620 01621 return 0; 01622 } 01623 01624 01625 01626 01627 // ===================================================================== 01628 // --------------------------------------------------------------------- 01629 // --------------------------------------------------------------------- 01630 // ===================================================================== 01631 /* 01632 \fn IntfWriteProjFile 01633 \param a_pathToImg : string containing the path to the image basename 01634 \param a2p_ImgMatrix : 2 dimensional image matrix which contains the image to write. 01635 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 01636 \param a_Imgfields: Structure containing information about the projected data 01637 \param vb : verbosity 01638 \brief Function dedicated to Interfile image writing for projected data 01639 \details Call the main functions dedicated to Interfile header and data writing 01640 Currently work for SPECT projected data 01641 The total number of projections should be provided in parameters 01642 Depending on the output writing mode (stored in sOutputManager), 01643 \todo Get metadata from a Intf_fields object ? 01644 (useful to transfer keys from read images to written images) 01645 \return 0 if success, positive value otherwise. 01646 */ 01647 int IntfWriteProjFile(const string& a_pathToImg, FLTNB** a2p_ImgMatrix, oImageDimensionsAndQuantification* ap_ID, Intf_fields a_Imgfields, int vb) 01648 { 01649 if (vb>=3) Cout("IntfWriteProjFile ..." << endl); 01650 01651 //Set image names 01652 vector<string> lpath_to_images; 01653 lpath_to_images.push_back(a_pathToImg); 01654 01655 if(IntfWriteHeaderMainData(a_pathToImg, a_Imgfields, vb) ) 01656 { 01657 Cerr("***** IntfWriteProjFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[0] << "' !" << endl); 01658 return 1; 01659 } 01660 01661 // Write interfile image 01662 01663 // Set dimensions 01664 uint32_t dims[2]; 01665 dims[0] = a_Imgfields.nb_projections; 01666 dims[1] = a_Imgfields.mtx_size[0]*a_Imgfields.mtx_size[1]; 01667 01668 if(IntfWriteImage(lpath_to_images, a2p_ImgMatrix, dims, vb) ) 01669 { 01670 Cerr("***** IntfWriteFile()-> Error : while trying to write the interfile image '"<< a_pathToImg << "' !" << endl); 01671 return 1; 01672 } 01673 01674 return 0; 01675 } 01676 01677 01678 01679 01680 // ===================================================================== 01681 // --------------------------------------------------------------------- 01682 // --------------------------------------------------------------------- 01683 // ===================================================================== 01684 /* 01685 \fn IntfWriteImgDynCoeffFile 01686 \param a_pathToImg : string containing the path to the image basename 01687 \param a2p_ImgMatrix : 2 dimensional image matrix which contains the image to write. 01688 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 01689 \param a_nbFbases : Number of basis functions 01690 \param vb : verbosity 01691 \brief Function dedicated to Interfile image writing for dynamic coefficients images 01692 \details Call the main functions dedicated to Interfile header and data writing 01693 The total number of basis functions should be provided in parameters 01694 Depending on the output writing mode (stored in sOutputManager), 01695 One image with one file and one will be created for the whole dynamic image 01696 or a metaheader and associated multiple 3D image raw file/headers will be generated 01697 \todo Get metadata from a Intf_fields object ? 01698 (useful to transfer keys from read images to written images) 01699 \return 0 if success, positive value otherwise. 01700 */ 01701 int IntfWriteImgDynCoeffFile(const string& a_pathToImg, FLTNB** a2p_ImgMatrix, oImageDimensionsAndQuantification* ap_ID, int a_nbFbases, int vb) 01702 { 01703 if (vb>=3) Cout("IntfWriteImgDynCoeffFile ..." << endl); 01704 01705 Intf_fields Img_fields; //TODO keep refs to Intf_fields in case we need it later (transfer keys from read images to written images) 01706 IntfKeySetFieldsOutput(&Img_fields, ap_ID); 01707 01708 // Get input from user about how dynamic files should be written 01709 // (one global image file, or separate image file for each frame/gates). 01710 sOutputManager* p_outputMgr = sOutputManager::GetInstance(); 01711 01712 //Set image names 01713 vector<string> lpath_to_images; 01714 01715 if(p_outputMgr->MergeDynImages() == true) // Case one file 01716 { 01717 lpath_to_images.push_back(a_pathToImg); 01718 01719 if(IntfWriteHeaderMainData(lpath_to_images[0], Img_fields, vb) ) 01720 { 01721 Cerr("***** IntfWriteImgDynCoeffFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[0] << "' !" << endl); 01722 return 1; 01723 } 01724 } 01725 else // Case several files 01726 { 01727 int idx_file = 0; 01728 01729 for(int bf=0 ; bf<a_nbFbases ; bf++) 01730 { 01731 //lpath_to_images[idx_file] = a_pathToImg; 01732 lpath_to_images.push_back(a_pathToImg); 01733 01734 // Add a suffix for basis functions 01735 stringstream ss; ss << bf + 1; 01736 lpath_to_images[idx_file].append("_bsf").append(ss.str()); 01737 01738 // Write header specific to this image 01739 string path_to_hdr = lpath_to_images[idx_file]; 01740 string path_to_img = lpath_to_images[idx_file]; 01741 01742 // As the content will be appended, make sure there is no existing file with such name 01743 // remove it otherwise 01744 ifstream fcheck(path_to_hdr.append(".hdr").c_str()); 01745 if(fcheck.good()) 01746 { 01747 fcheck.close(); 01748 #ifdef _WIN32 01749 string dos_instruction = "del " + path_to_hdr; 01750 system(dos_instruction.c_str()); 01751 #else 01752 remove(path_to_hdr.c_str()); 01753 #endif 01754 } 01755 01756 ofstream ofile(path_to_hdr.c_str(), ios::out | ios::app); 01757 ofile << setprecision(std::numeric_limits<FLTNB>::digits10 +1); 01758 ofile << "!INTERFILE := " << endl; 01759 ofile << "!name of data file := " << GetFileFromPath(path_to_img).append(".img") << endl; 01760 ofile << "imagedata byte order := " << IntfKeyGetEndianStr(User_Endianness) << endl; 01761 ofile << "!data offset in bytes := " << 0 << endl; 01762 01763 if(IntfWriteHeaderImgData(ofile, Img_fields, vb) ) 01764 { 01765 Cerr("***** IntfWriteImgDynCoeffFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[idx_file] << "' !" << endl); 01766 return 1; 01767 } 01768 01769 ofile.close(); 01770 idx_file++; 01771 } 01772 01773 // Write meta header 01774 if(IntfWriteMHD(a_pathToImg, lpath_to_images, Img_fields, ap_ID, vb) ) 01775 { 01776 Cerr("***** IntfWriteImgDynCoeffFile()-> Error : while trying to write the interfile header '"<< a_pathToImg << "' !" << endl); 01777 return 1; 01778 } 01779 } 01780 01781 // Write interfile image 01782 01783 // Set dimensions 01784 uint32_t dims[2]; 01785 dims[0] = a_nbFbases; 01786 dims[1] = ap_ID->GetNbVoxXYZ(); 01787 01788 if(IntfWriteImage(lpath_to_images, a2p_ImgMatrix, dims, vb) ) 01789 { 01790 Cerr("***** IntfWriteImgDynCoeffFile()-> Error : while trying to write the interfile image '"<< a_pathToImg << "' !" << endl); 01791 return 1; 01792 } 01793 01794 return 0; 01795 } 01796 01797 01798 01799 01800 // ===================================================================== 01801 // --------------------------------------------------------------------- 01802 // --------------------------------------------------------------------- 01803 // ===================================================================== 01804 /* 01805 \fn IntfWriteImgFile 01806 \param a_pathToImg : string containing the path to the image basename 01807 \param a4p_ImgMatrix : 4 dimensional image matrix which contains the image to write. 01808 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 01809 \param vb : verbosity 01810 \brief Main function dedicated to Interfile 6D (dynamic dims + 3D ) image writing 01811 \details Call the main functions dedicated to Interfile header and data writing 01812 Depending on the output writing mode (stored in sOutputManager), 01813 One image with one file and one will be created for the whole dynamic image 01814 or a metaheader and associated multiple 3D image raw file/headers will be generated 01815 \todo Get metadata from a Intf_fields object ? 01816 (useful to transfer keys from read images to written images) 01817 \return 0 if success, positive value otherwise. 01818 */ 01819 int IntfWriteImgFile(const string& a_pathToImg, FLTNB**** a4p_ImgMatrix, oImageDimensionsAndQuantification* ap_ID, int vb) 01820 { 01821 if (vb>=3) Cout("IntfWriteImgFile (dynamic/static image) ..." << endl); 01822 01823 Intf_fields Img_fields; //TODO keep refs to Intf_fields in case we need it later (transfer keys from read images to written images) 01824 IntfKeySetFieldsOutput(&Img_fields, ap_ID); 01825 01826 // Get input from user about how dynamic files should be written 01827 // (one global image file, or separate image file for each frame/gates). 01828 sOutputManager* p_outputMgr = sOutputManager::GetInstance(); 01829 01830 //Set image names 01831 vector<string> lpath_to_images; 01832 01833 if(Img_fields.nb_dims <=3 || p_outputMgr->MergeDynImages() == true) // One file to write 01834 { 01835 lpath_to_images.push_back(a_pathToImg); 01836 01837 if(IntfWriteHeaderMainData(a_pathToImg, Img_fields, vb) ) 01838 { 01839 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[0] << "' !" << endl); 01840 return 1; 01841 } 01842 } 01843 else // Case several files 01844 { 01845 int idx_file = 0; 01846 01847 for(int fr=0 ; fr<ap_ID->GetNbTimeFrames() ; fr++) 01848 for(int rg=0 ; rg<ap_ID->GetNbRespGates() ; rg++) 01849 for(int cg=0 ; cg<ap_ID->GetNbCardGates() ; cg++) 01850 { 01851 //lpath_to_images[idx_file] = a_pathToImg; 01852 lpath_to_images.push_back(a_pathToImg); 01853 // Add a suffix for time frames if more than 1 01854 if (ap_ID->GetNbTimeFrames() > 1) 01855 { 01856 stringstream ss; ss << fr + 1; 01857 lpath_to_images[idx_file].append("_fr").append(ss.str()); 01858 } 01859 // Add a suffix for respiratory gates if more than 1 01860 if (ap_ID->GetNbRespGates() > 1) 01861 { 01862 stringstream ss; ss << rg + 1; 01863 lpath_to_images[idx_file].append("_rg").append(ss.str()); 01864 } 01865 // Add a suffix for cardiac gates if more than 1 01866 if (ap_ID->GetNbCardGates() > 1) 01867 { 01868 stringstream ss; ss << cg + 1; 01869 lpath_to_images[idx_file].append("_cg").append(ss.str()); 01870 } 01871 01872 // Write header specific to this image 01873 string path_to_hdr = lpath_to_images[idx_file]; 01874 string path_to_img = lpath_to_images[idx_file]; 01875 01876 // As the content will be appended, make sure there is no existing file with such name 01877 // remove it otherwise 01878 ifstream fcheck(path_to_hdr.append(".hdr").c_str()); 01879 if(fcheck.good()) 01880 { 01881 fcheck.close(); 01882 #ifdef _WIN32 01883 string dos_instruction = "del " + path_to_hdr; 01884 system(dos_instruction.c_str()); 01885 #else 01886 remove(path_to_hdr.c_str()); 01887 #endif 01888 } 01889 01890 ofstream ofile(path_to_hdr.c_str(), ios::out | ios::app); 01891 ofile << setprecision(std::numeric_limits<FLTNB>::digits10 +1); 01892 ofile << "!INTERFILE := " << endl; 01893 ofile << "!name of data file := " << GetFileFromPath(path_to_img).append(".img") << endl; 01894 ofile << "imagedata byte order := " << IntfKeyGetEndianStr(User_Endianness) << endl; 01895 ofile << "!data offset in bytes := " << 0 << endl; 01896 01897 if(IntfWriteHeaderImgData(ofile, Img_fields, vb) ) 01898 { 01899 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[idx_file] << "' !" << endl); 01900 return 1; 01901 } 01902 01903 // Write image duration related data (this information will be replicated on each gated image) 01904 ofile << "image duration (sec) := " << (ap_ID->GetFrameTimeStopInSec(0, fr) - ap_ID->GetFrameTimeStartInSec(0, fr) )<< endl; 01905 01906 // same start time for all gates 01907 ofile << "image start time (sec) := " << ap_ID->GetFrameTimeStartInSec(0, fr) << endl; 01908 01909 ofile.close(); 01910 01911 idx_file++; 01912 } 01913 01914 // Write meta header 01915 if(IntfWriteMHD(a_pathToImg, lpath_to_images, Img_fields, ap_ID, vb) ) 01916 { 01917 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header '"<< a_pathToImg << "' !" << endl); 01918 return 1; 01919 } 01920 } 01921 01922 // Write interfile image 01923 01924 // Set dimensions 01925 uint32_t dims[4]; 01926 dims[0] = ap_ID->GetNbTimeFrames(); 01927 dims[1] = ap_ID->GetNbRespGates(); 01928 dims[2] = ap_ID->GetNbCardGates(); 01929 dims[3] = ap_ID->GetNbVoxXYZ(); 01930 01931 if(IntfWriteImage(lpath_to_images, a4p_ImgMatrix, dims, vb) ) 01932 { 01933 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile image '"<< a_pathToImg << "' !" << endl); 01934 return 1; 01935 } 01936 01937 return 0; 01938 } 01939 01940 01941 01942 01943 01944 // ===================================================================== 01945 // --------------------------------------------------------------------- 01946 // --------------------------------------------------------------------- 01947 // ===================================================================== 01948 /* 01949 \fn IntfWriteImgFile 01950 \param a_pathToImg : string containing the path to the image basename 01951 \param a4p_ImgMatrix : 4 dimensional image matrix which contains the image to write. 01952 \param ap_ID : Provide the Image dimensions object containing reconstruction dimensions 01953 \param vb : verbosity 01954 \brief Main function dedicated to Interfile 6D (dynamic dims + 3D ) image writing 01955 \details Call the main functions dedicated to Interfile header and data writing 01956 Depending on the output writing mode (stored in sOutputManager), 01957 One image with one file and one will be created for the whole dynamic image 01958 or a metaheader and associated multiple 3D image raw file/headers will be generated 01959 \todo Get metadata from a Intf_fields object ? 01960 (useful to transfer keys from read images to written images) 01961 \return 0 if success, positive value otherwise. 01962 01963 int IntfWriteSensitivityImg(const string& a_pathToImg, FLTNB**** a4p_ImgMatrix, oImageDimensionsAndQuantification* ap_ID, int vb) 01964 { 01965 if(vb >= 2) Cout("IntfWriteSensitivityImg ..." << endl); 01966 01967 Intf_fields Img_fields; //TODO keep refs to Intf_fields in case we need it later (transfer keys from read images to written images) 01968 IntfKeySetFieldsOutput(&Img_fields, ap_ID); 01969 01970 // Get input from user about how dynamic files should be written 01971 // (one global image file, or separate image file for each frame/gates). 01972 sOutputManager* p_outputMgr = sOutputManager::GetInstance(); 01973 01974 //Set image names 01975 vector<string> lpath_to_images; 01976 01977 // By default, sensitivity image are written in one file One file to write 01978 { 01979 lpath_to_images.push_back(a_pathToImg); 01980 01981 if(IntfWriteHeaderMainData(a_pathToImg, Img_fields, ap_ID, vb) ) 01982 { 01983 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile header for'"<< lpath_to_images[0] << "' !" << endl); 01984 return 1; 01985 } 01986 } 01987 01988 // Write interfile image 01989 01990 // Set dimensions 01991 uint32_t dims[4]; 01992 dims[0] = ap_ID->GetNbTimeFrames(); 01993 dims[1] = ap_ID->GetNbRespGates(); 01994 dims[2] = ap_ID->GetNbCardGates(); 01995 dims[3] = ap_ID->GetNbVoxXYZ(); 01996 01997 if(IntfWriteImage(lpath_to_images, a4p_ImgMatrix, dims, vb) ) 01998 { 01999 Cerr("***** IntfWriteImgFile()-> Error : while trying to write the interfile image '"<< a_pathToImg << "' !" << endl); 02000 return 1; 02001 } 02002 02003 return 0; 02004 } 02005 */ 02006 02007 02008 // ===================================================================== 02009 // --------------------------------------------------------------------- 02010 // --------------------------------------------------------------------- 02011 // ===================================================================== 02012 /* 02013 \fn IntfIsMHD 02014 \param a_pathToFile : String containing path to an Interfile header 02015 \param ap_lPathToImgs : pointer to a list of strings containing the path to the different images 02016 \brief Check if the string in argument contains the path to a Interfile metaheader 02017 \details Check for '!total number of data sets' key, and the associated image file names 02018 If the file is metaheader, the names of image file header will be returned in ap_lPathToImgs list 02019 It not, the file is a single header, that we add to the list as solo file 02020 \return 1 if we have a metaheader, 0 if not, negative value if error. 02021 */ 02022 int IntfIsMHD(string a_pathToFile, vector<string> &ap_lPathToImgs) 02023 { 02024 // Create file object and check existence 02025 ifstream hfile(a_pathToFile.c_str(), ios::in); 02026 string line; 02027 02028 if(hfile) 02029 { 02030 // Read until the end of file and check if we have a '!total number of data sets' key 02031 while(!hfile.eof()) 02032 { 02033 getline(hfile, line); 02034 02035 // Check if empty line. 02036 if(line.empty()) continue; 02037 02038 // Initiate a Key which will recover the potential interfile data for each line 02039 Intf_key Key; 02040 02041 // Process the Key at this line. Recover the Key, its counterpart in lower case, and the value in the structure Intf_key 02042 if (IntfRecoverKey(&Key, line) ) 02043 { 02044 Cerr("***** IntfIsMHD()-> Error : couldn't correctly read interfile key '"<< line << "' !" << endl); 02045 return -1; 02046 } 02047 02048 // --- Check if a total number of data set if provided 02049 if(IntfCheckKeyMatch(Key, "total number of data set") || 02050 IntfCheckKeyMatch(Key, "total number of data sets") || 02051 IntfCheckKeyMatch(Key, "!total number of data set") || 02052 IntfCheckKeyMatch(Key, "!total number of data sets")) 02053 { 02054 // File is MHD 02055 02056 uint16_t nb_datasets = 0; 02057 if(ConvertFromString(Key.kvalue , &nb_datasets) ) 02058 { 02059 Cerr("***** IntfIsMHD()-> Error when trying to read data from 'total number of data set' key. Recovered value = " << Key.kvalue << endl); 02060 return -1; 02061 } 02062 02063 // Get the image file path from the following lines 02064 int file_idx = 0; 02065 stringstream ss; 02066 ss << (file_idx+1); 02067 string file_key = ""; 02068 string file_keyp = "%"; 02069 string file_key_space = ""; 02070 string file_keyp_space = "%"; 02071 file_key_space.append("data set [").append(ss.str()).append("]"); 02072 file_keyp_space.append("data set [").append(ss.str()).append("]"); 02073 file_key.append("data set[").append(ss.str()).append("]"); 02074 file_keyp.append("data set[").append(ss.str()).append("]"); 02075 02076 string path_to_mhd_dir = GetPathOfFile(a_pathToFile); 02077 02078 while(!hfile.eof()) 02079 { 02080 // Read the following lines... 02081 getline(hfile, line); 02082 02083 // Check if empty line. 02084 if(line.empty()) continue; 02085 02086 // Process the Key at this line. Recover the Key, its counterpart in lower case, and the value in the structure Intf_key 02087 if (IntfRecoverKey(&Key, line) ) 02088 { 02089 Cerr("***** IntfIsMHD()-> Error : couldn't correctly read interfile key '"<< line << "' !" << endl); 02090 return -1; 02091 } 02092 02093 if(IntfCheckKeyMatch(Key, file_key) || 02094 IntfCheckKeyMatch(Key, file_keyp) || 02095 IntfCheckKeyMatch(Key, file_key_space) || 02096 IntfCheckKeyMatch(Key, file_keyp_space)) 02097 { 02098 ap_lPathToImgs.push_back(GetPathOfFile(a_pathToFile)); 02099 if(IntfKeyIsArray(Key)) 02100 { 02101 string elts_str[3]; 02102 02103 // First, check we have the correct number of elts 02104 int nb_elts = IntfKeyGetArrayNbElts(Key); 02105 if(nb_elts<0) 02106 { 02107 Cerr("***** IntfIsMHD()-> Error : when trying to read Interfile array key: '"<< line << "' !" << endl); 02108 return -1; 02109 } 02110 else if(nb_elts != 3) 02111 { 02112 Cerr("***** IntfIsMHD()-> Error : when trying to read Interfile array key: '"<< line << "' !" << endl); 02113 Cerr(" 3 elements are expected following the format {xxx, path_to_the_img_file, xxx} (xxx = ignored data)" << endl); 02114 return -1; 02115 } 02116 02117 if(IntfKeyGetArrayElts(Key, elts_str) ) 02118 { 02119 Cerr("***** IntfIsMHD()-> Error : when trying to read Interfile array key: '"<< line << "' !" << endl); 02120 return -1; 02121 } 02122 ap_lPathToImgs.at(file_idx).append(elts_str[1]); 02123 } 02124 else 02125 ap_lPathToImgs.at(file_idx).append(Key.kvalue); 02126 02127 // Compute next key to recover 02128 file_idx++; 02129 stringstream ss; 02130 ss << (file_idx+1); 02131 file_key = ""; 02132 file_keyp = "%"; 02133 file_key.append("data set [").append(ss.str()).append("]"); 02134 file_keyp.append("data set [").append(ss.str()).append("]"); 02135 02136 } 02137 } 02138 02139 // Check we have a correct nb of images 02140 if(nb_datasets != ap_lPathToImgs.size()) 02141 { 02142 // error if not found 02143 Cerr("***** IntfIsMHD()-> The number of recovered file in the metaheader ("<<ap_lPathToImgs.size()<<")"); 02144 Cerr("does not correspond to the expected number of datasets ("<<nb_datasets<<")!"<< endl); 02145 return -1; 02146 } 02147 else 02148 // Normal end, file is metaheader 02149 return 1; 02150 } 02151 } 02152 } 02153 else 02154 { 02155 Cerr("***** IntfIsMHD()-> Error : couldn't read header file '"<< a_pathToFile << "' !" << endl); 02156 return -1; 02157 } 02158 02159 // If we reach this, the file is NOT considered as metaheader the 'total number of data set' key was not found 02160 // Add the to the list as unique header file 02161 ap_lPathToImgs.push_back(a_pathToFile); 02162 return 0; 02163 } 02164 02165 02166 02167 02168 // ===================================================================== 02169 // --------------------------------------------------------------------- 02170 // --------------------------------------------------------------------- 02171 // ===================================================================== 02172 /* 02173 \fn IntfWriteMHD 02174 \param a_path : path to the Meta interfile header to write 02175 \param ap_lPathToImgs : pointer to a list of strings containing the path to the different images 02176 \param a_IntfF : Structure containing interfile keys of the image to write 02177 \param ap_ID : oImageDimensionsAndQuantification object containing additional infos about reconstruction 02178 \param vb : verbosity 02179 \brief Write an Interfile meta header at the path provided in parameter, using the field stack provided in parameter. 02180 \todo keys for multi bed 02181 \return 0 if success, positive value otherwise. 02182 */ 02183 int IntfWriteMHD( const string& a_pathToMhd, 02184 const vector<string>& ap_lPathToImgs, 02185 Intf_fields a_IntfF, 02186 oImageDimensionsAndQuantification* ap_ID, 02187 int vb) 02188 { 02189 // Get file and check existence 02190 string path_to_mhd_file = a_pathToMhd; 02191 path_to_mhd_file.append(".mhd"); 02192 ofstream ofile(path_to_mhd_file.c_str(), ios::out); 02193 ofile << setprecision(std::numeric_limits<FLTNB>::digits10 +1); 02194 02195 if(ofile) 02196 { 02197 sScannerManager* p_scanMgr = sScannerManager::GetInstance(); 02198 02199 // todo fields for multi-beds ? 02200 02201 ofile << "!INTERFILE := " << endl; 02202 //ofile << "% MetaHeader written by CASToRv1.0 := " << end; 02203 ofile << endl << "!GENERAL DATA := " << endl; 02204 ofile << "data description:=image" << endl; 02205 ofile << "!originating system := " << p_scanMgr->GetScannerName() << endl; 02206 02207 ofile << endl << "%DATA MATRIX DESCRIPTION:=" << endl; 02208 ofile << "number of time frames := " << a_IntfF.nb_time_frames << endl; 02209 ofile << "number of time windows := " << a_IntfF.nb_resp_gates * 02210 a_IntfF.nb_card_gates << endl; 02211 ofile << "number of respiratory gates := " << a_IntfF.nb_resp_gates << endl; 02212 ofile << "number of cardiac gates := " << a_IntfF.nb_card_gates << endl; 02213 02214 ofile << endl << "%DATA SET DESCRIPTION:="<< endl; 02215 02216 uint16_t nb_imgs = a_IntfF.nb_time_frames * a_IntfF.nb_resp_gates * a_IntfF.nb_card_gates ; 02217 ofile << "!total number of data sets:=" << nb_imgs << endl; 02218 02219 // Check we have the right number of strings 02220 if(ap_lPathToImgs.size() != nb_imgs) 02221 { 02222 Cerr("***** IntfWriteMHD()-> Error : nb of provided string inconsistent with the expected number of dynamic images: '"<< nb_imgs << "' !" << endl); 02223 ofile.close(); 02224 return 1; 02225 } 02226 02227 for(int ds=0 ; ds<nb_imgs ; ds++) 02228 { 02229 string path_to_header = ap_lPathToImgs.at(ds); 02230 ofile << "%data set ["<<ds+1<<"]:={0," << GetFileFromPath(path_to_header).append(".hdr") << ",UNKNOWN}"<< endl; 02231 } 02232 } 02233 else 02234 { 02235 Cerr("***** IntfWriteMHD()-> Error : couldn't find output Metaheader interfile '"<< path_to_mhd_file << "' !" << endl); 02236 return 1; 02237 } 02238 02239 ofile.close(); 02240 02241 return 0; 02242 } 02243 02244 02245 02246 02247 // ===================================================================== 02248 // --------------------------------------------------------------------- 02249 // --------------------------------------------------------------------- 02250 // ===================================================================== 02251 /* 02252 \fn IntfReadHeader 02253 \param const string& a_pathToHeaderFile 02254 \param Intf_fields* ap_IF 02255 \param int vb : Verbosity level 02256 \brief Read an Interfile header 02257 \details Initialize all mandatory fields from the Intf_fields structure passed in parameter with the related interfile key from the image interfile header 02258 \todo Check everything work for Interfile 3.3 version, extended interfile, etc. 02259 \todo Several keys from Intf_fields still needs some management 02260 \todo Deal with sinogram in interfile format (matrix size keys) 02261 \todo orientation related keys 02262 \todo time windows, energy windows keys 02263 \todo check start angle key 02264 \todo Initializations & checks at the end (ok for slice_thickness, check slice_spacing, etc..) 02265 \todo Specification of this function for projection data (currently we consider the interfile is a 3D image) 02266 \return 0 if success, positive value otherwise. 02267 */ 02268 int IntfReadHeader(const string& a_pathToHeaderFile, Intf_fields* ap_IF, int vb) 02269 { 02270 if (vb >= 4) 02271 { 02272 Cout("--------------------------------------------------------------- " << endl); 02273 Cout("IntfReadHeader()-> Start reading header interfile " << a_pathToHeaderFile << endl); 02274 Cout("--------------------------------------------------------------- " << endl); 02275 } 02276 02277 // Get file and check existence 02278 ifstream input_file(a_pathToHeaderFile.c_str(), ios::in); 02279 string line; 02280 02281 if(input_file) 02282 { 02283 // Read until the end of file 02284 while(!input_file.eof()) 02285 { 02286 getline(input_file, line); 02287 02288 // Check if empty line. 02289 if(line.empty()) 02290 continue; 02291 02292 // Initiate a Key which will recover the potential interfile data for each line 02293 Intf_key Key; 02294 02295 // Process the Key at this line. Recover the Key, its counterpart in lower case, and the value in the structure Intf_key 02296 if (IntfRecoverKey(&Key, line) ) 02297 { 02298 Cerr("***** ReadIntfHeader()-> Error : couldn't correctly read interfile key '"<< line << "' !" << endl); 02299 return 1; 02300 } 02301 02302 if (vb >=10) Cout("ReadIntfHeader()-> Key " << Key.kcase << endl); 02303 02304 02305 // --- From there, find key Id and process result --- 02306 02307 // Check version. Warning message if not 3.3 / Castor keys ? 02308 if(IntfCheckKeyMatch(Key, "version of keys")) 02309 { 02310 //if(Key.kvalue != 3.3 || 02311 // Key.kvalue != CASToRv1.0) 02312 // Cout("***** ReadIntfHeader()-> Warning : Unexpected version of Interfile keys found: " << Key.kvalue << " (supported version = 3.3)" << endl); 02313 } 02314 02315 02316 // --- path to image file 02317 // "name of data file" / Intf_fields.path_to_image 02318 else if(IntfCheckKeyMatch(Key, "name of data file")) 02319 { 02320 // Check key isn't empty, throw error otherwise 02321 if ( Key.kvalue.size()>0 ) 02322 { 02323 // Look for path separator 02324 if ( Key.kvalue.find(OS_SEP) != string::npos ) 02325 { 02326 // Assuming we have an absolute path 02327 ap_IF->path_to_image = Key.kvalue; 02328 } 02329 else 02330 { 02331 // Assuming we just have the image name 02332 // -> use the relative path of the header file, image should be in the same dir 02333 string header_file_directory = GetPathOfFile(a_pathToHeaderFile); 02334 ap_IF->path_to_image = header_file_directory.append(Key.kvalue); 02335 } 02336 } 02337 else 02338 { 02339 Cerr("***** ReadIntfHeader()-> Error : path to interfile image file is empty !" << endl); 02340 return 1; 02341 } 02342 } 02343 02344 02345 // --- Little or big endian 02346 // "imagedata byte order" / Intf_fields.endianness 02347 else if(IntfCheckKeyMatch(Key, "imagedata byte order")) 02348 { 02349 string endianness; 02350 IntfToLowerCase(&Key.kvalue); // whole string in low caps 02351 02352 endianness = Key.kvalue; 02353 02354 if (endianness == "littleendian") 02355 ap_IF->endianness = 1; 02356 else // Big endian by default 02357 ap_IF->endianness = INTF_BIG_ENDIAN; 02358 } 02359 02360 02361 // --- Data offset using "data starting block" 02362 // "data starting block" & "data offset in bytes" / Intf_fields.data_offset 02363 else if(IntfCheckKeyMatch(Key, "data starting block")) 02364 { 02365 // Throw warning indicating that the data offset has already been set 02366 if(ap_IF->data_offset>0 ) 02367 { 02368 Cerr("***** ReadIntfHeader()-> Warning : data_offset value was already set to " << ap_IF->data_offset << endl); 02369 Cerr("***** The header may contain both 'data offset in bytes' and 'data starting block' fields " << endl); 02370 } 02371 02372 if(ConvertFromString(Key.kvalue , &ap_IF->data_offset) ) 02373 { 02374 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'data starting block' key. Recovered value = " << Key.kvalue << endl); 02375 return 1; 02376 } 02377 02378 ap_IF->data_offset*= 2048L; 02379 } 02380 02381 02382 // --- Data offset using "data offset in bytes" 02383 // "data starting block" & "data offset in bytes" / Intf_fields.data_offset 02384 else if(IntfCheckKeyMatch(Key, "data offset in bytes")) 02385 { 02386 // Throw warning indicating that the data offset has already been set 02387 if(ap_IF->data_offset>0 ) 02388 { 02389 Cerr("***** ReadIntfHeader()-> Warning : data_offset value was already set to " << ap_IF->data_offset << endl); 02390 Cerr("***** The header may contain both 'data offset in bytes' and 'data starting block' fields " << endl); 02391 } 02392 02393 if(ConvertFromString(Key.kvalue , &ap_IF->data_offset) ) 02394 { 02395 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'data offset in bytes' key. Recovered value = " << Key.kvalue << endl); 02396 return 1; 02397 } 02398 } 02399 02400 02401 // --- Type of each pixel/voxel data 02402 // "number format" / Intf_fields.nb_format 02403 else if(IntfCheckKeyMatch(Key, "number format")) 02404 { 02405 if(ConvertFromString(Key.kvalue , &ap_IF->nb_format) ) 02406 { 02407 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number format' key. Recovered value = " << Key.kvalue << endl); 02408 return 1; 02409 } 02410 } 02411 02412 02413 // --- Width 02414 // "matrix size [1]" / Intf_fields.mtx_size[0] 02415 else if(IntfCheckKeyMatch(Key, "matrix size [1]")) 02416 { 02417 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[0]) ) 02418 { 02419 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [1]' key. Recovered value = " << Key.kvalue << endl); 02420 return 1; 02421 } 02422 } 02423 02424 // --- Height 02425 // "matrix size [2]" / Intf_fields.mtx_size[1] 02426 else if(IntfCheckKeyMatch(Key, "matrix size [2]")) 02427 { 02428 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[1]) ) 02429 { 02430 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [2]' key. Recovered value = " << Key.kvalue << endl); 02431 return 1; 02432 } 02433 } 02434 02435 // --- Slices : Tomographic image 02436 // "matrix size [3]" / Intf_fields.mtx_size[2] 02437 else if(IntfCheckKeyMatch(Key, "matrix size [3]")) 02438 { 02439 // Check if we have an array key 02440 if (IntfKeyIsArray(Key)) 02441 { 02442 //TODO : perhaps more convenient to use vector to deal with data whose matrix size [3] has more than 1 element (eg: sinograms) 02443 // It will depend how we decide to handle sinogram reading 02444 // For now, throw error by default as there is no implementation yet 02445 02446 // Throw error by default 02447 Cerr("***** ReadIntfHeader()-> Array reading for 'matrix size [3]' key not implemented yet. Single value expected." << endl); 02448 return 1; 02449 } 02450 else 02451 { 02452 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[2]) ) 02453 { 02454 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [3]' key. Recovered value = " << Key.kvalue << endl); 02455 return 1; 02456 } 02457 } 02458 } 02459 02460 // --- Read nb frames or other. 02461 // "matrix size [4]" / Intf_fields.mtx_size[3] 02462 // Consistency with "number of time frames" : Priority to number of time frames if this has already been initialized 02463 // TODO : this may cause issue for implementation of sinogram management 02464 else if(IntfCheckKeyMatch(Key, "matrix size [4]")) 02465 { 02466 // Could be sinogram or dynamic related (nb time frames). 02467 // Dynamic reconstructions : consistency with "number of frames" 02468 02469 if (IntfKeyIsArray(Key)) 02470 { 02471 //TODO : perhaps more convenient to use vector to deal with data whose matrix size [4] has more than 1 element (eg: sinograms) 02472 // It will depend how we decide to handle sinogram reading 02473 // For now, throw error by default as there is no implementation yet 02474 02475 // Throw error by default 02476 Cerr("***** ReadIntfHeader()-> Array reading for 'matrix size [4]' key not implemented yet. Single value expected." << endl); 02477 return 1; 02478 } 02479 else 02480 { 02481 // Check if number time frames has already been initialized 02482 if(ap_IF->nb_time_frames>1) 02483 { 02484 ap_IF->mtx_size[3] = ap_IF->nb_time_frames; 02485 Cerr("***** ReadIntfHeader()-> WARNING : both 'number of time frames' and 'matrix size [4]' keys have been provided"); 02486 Cerr(" 'number of time frames' selected by default" << endl); 02487 } 02488 else if( ConvertFromString(Key.kvalue , &ap_IF->mtx_size[3]) ) 02489 { 02490 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [4]' key. Recovered value = " << Key.kvalue << endl); 02491 return 1; 02492 } 02493 } 02494 } 02495 02496 02497 // --- Related to sinogram (not implemented yet) 02498 // "matrix size [5]" / Intf_fields.mtx_size[4] 02499 else if(IntfCheckKeyMatch(Key, "matrix size [5]")) 02500 { 02501 // Could be sinogram or dynamic related 02502 if (IntfKeyIsArray(Key)) 02503 { 02504 //TODO : perhaps more convenient to use vector to deal with data whose matrix size [5] has more than 1 element (eg: sinograms) 02505 // It will depend how we decide to handle sinogram reading 02506 // For now, throw error by default as there is no implementation yet 02507 02508 // Throw error by default 02509 Cerr("***** ReadIntfHeader()-> Array reading for 'matrix size [5]' key not implemented yet. Single value expected." << endl); 02510 return 1; 02511 } 02512 else 02513 { 02514 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[4]) ) 02515 { 02516 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [5]' key. Recovered value = " << Key.kvalue << endl); 02517 return 1; 02518 } 02519 } 02520 } 02521 02522 // --- Related to sinogram (not implemented yet) 02523 // "matrix size [6]" / Intf_fields.mtx_size[5] 02524 else if(IntfCheckKeyMatch(Key, "matrix size [6]")) 02525 { 02526 // Could be sinogram or dynamic related 02527 if (IntfKeyIsArray(Key)) 02528 { 02529 //TODO : perhaps more convenient to use vector to deal with data whose matrix size [6] has more than 1 element (eg: sinograms) 02530 // It will depend how we decide to handle sinogram reading 02531 // For now, throw error by default as there is no implementation yet 02532 02533 // Throw error by default 02534 Cerr("***** ReadIntfHeader()-> Array reading for 'matrix size [6]' key not implemented yet. Single value expected." << endl); 02535 return 1; 02536 } 02537 else 02538 { 02539 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[5]) ) 02540 { 02541 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [6]' key. Recovered value = " << Key.kvalue << endl); 02542 return 1; 02543 } 02544 } 02545 } 02546 02547 // --- Related to sinogram (not implemented yet) 02548 // "matrix size [7]" / Intf_fields.mtx_size[6] 02549 else if(IntfCheckKeyMatch(Key, "matrix size [7]")) 02550 { 02551 // Could be sinogram related 02552 if (IntfKeyIsArray(Key)) 02553 { 02554 //TODO : perhaps more convenient to use vector to deal with data whose matrix size [7] has more than 1 element (eg: sinograms) 02555 // It will depend how we decide to handle sinogram reading 02556 // For now, throw error by default as there is no implementation yet 02557 02558 // Throw error by default 02559 Cerr("***** ReadIntfHeader()-> Array reading for 'matrix size [7]' key not implemented yet. Single value expected." << endl); 02560 return 1; 02561 } 02562 else 02563 { 02564 if( ConvertFromString(Key.kvalue , &ap_IF->mtx_size[6]) ) 02565 { 02566 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'matrix size [7]' key. Recovered value = " << Key.kvalue << endl); 02567 return 1; 02568 } 02569 } 02570 02571 } 02572 02573 // --- Size voxel width 02574 // "scaling factor (mm/pixel) [1]" / Intf_fields.vox_size[0] 02575 else if(IntfCheckKeyMatch(Key, "scaling factor (mm/pixel) [1]")) 02576 { 02577 if(ConvertFromString(Key.kvalue , &ap_IF->vox_size[0]) ) 02578 { 02579 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'scaling factor (mm/pixel) [1]' key. Recovered value = " << Key.kvalue << endl); 02580 return 1; 02581 } 02582 } 02583 02584 // --- Size voxel height 02585 // "scaling factor (mm/pixel) [2]" / Intf_fields.vox_size[1] 02586 else if(IntfCheckKeyMatch(Key, "scaling factor (mm/pixel) [2]")) 02587 { 02588 if(ConvertFromString(Key.kvalue , &ap_IF->vox_size[1]) ) 02589 { 02590 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'scaling factor (mm/pixel) [2]' key. Recovered value = " << Key.kvalue << endl); 02591 return 1; 02592 } 02593 } 02594 02595 // --- Size voxel axial 02596 // "scaling factor (mm/pixel) [3]" / Intf_fields.vox_size[2] 02597 // Prefered over slice thickness by default 02598 else if(IntfCheckKeyMatch(Key, "scaling factor (mm/pixel) [3]")) 02599 { 02600 if(ConvertFromString(Key.kvalue , &ap_IF->vox_size[2]) ) 02601 { 02602 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'scaling factor (mm/pixel) [3]' key. Recovered value = " << Key.kvalue << endl); 02603 return 1; 02604 } 02605 } 02606 02607 // --- Size thickness mm 02608 // "slice thickness (pixels)" / Intf_fields.slice_thickness_mm 02609 // We assuming the slice thickness is given in mm regardless of the name 02610 else if(IntfCheckKeyMatch(Key, "slice thickness (pixels)")) 02611 { 02612 if(ConvertFromString(Key.kvalue , &ap_IF->slice_thickness_mm) ) 02613 { 02614 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'slice thickness (pixels)' key. Recovered value = " << Key.kvalue << endl); 02615 return 1; 02616 } 02617 } 02618 02619 // --- centre-centre slice separation (pixels). (used to compute slice spacing) 02620 // "process status" / Intf_fields.ctr_to_ctr_separation 02621 else if(IntfCheckKeyMatch(Key, "centre-centre slice separation (pixels)") || 02622 IntfCheckKeyMatch(Key, "center-center slice separation (pixels)") ) 02623 { 02624 if(ConvertFromString(Key.kvalue ,&ap_IF->ctr_to_ctr_separation) ) 02625 { 02626 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'centre-centre slice separation (pixels)' key."); 02627 Cerr(" Recovered value = " << Key.kvalue << endl); 02628 return 1; 02629 } 02630 Cerr("***** ReadIntfHeader()-> WARNING : 'centre-centre slice separation (pixels)' has no use in the current implementation !"<< endl); 02631 } 02632 02633 // --- Number of time frames 02634 // "number of time frames" & "number of frame groups" / Intf_fields.nb_time_frames 02635 // Consistency with matrix size 4 : priority to nb_time_frames 02636 else if(IntfCheckKeyMatch(Key, "number of time frames") || 02637 IntfCheckKeyMatch(Key, "number of frame groups")) 02638 { 02639 if( ConvertFromString(Key.kvalue , &ap_IF->nb_time_frames) ) 02640 { 02641 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of time frames' or 'number of frame groups' keys."); 02642 Cerr("Recovered value = " << Key.kvalue << endl); 02643 return 1; 02644 } 02645 } 02646 02647 02648 // --- Number of respiratory gates 02649 // "number of respiratory gates" / Intf_fields.nb_resp_gates 02650 else if(IntfCheckKeyMatch(Key, "number of respiratory gates")) 02651 { 02652 if( ConvertFromString(Key.kvalue , &ap_IF->nb_resp_gates) ) 02653 { 02654 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of respiratory gates' key. Recovered value = " << Key.kvalue << endl); 02655 return 1; 02656 } 02657 } 02658 02659 02660 // --- Number of cardiac gates 02661 // "number of cardiac gates" / Intf_fields.nb_card_gates 02662 else if(IntfCheckKeyMatch(Key, "number of cardiac gates")) 02663 { 02664 if( ConvertFromString(Key.kvalue , &ap_IF->nb_card_gates) ) 02665 { 02666 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of cardiac gates' key. Recovered value = " << Key.kvalue << endl); 02667 return 1; 02668 } 02669 } 02670 02671 // --- Total number of images 02672 // "total number of images" / Intf_fields.nb_total_imgs 02673 else if(IntfCheckKeyMatch(Key, "total number of images") ) 02674 { 02675 if(ConvertFromString(Key.kvalue , &ap_IF->nb_total_imgs) ) 02676 { 02677 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'total number of images' key. Recovered value = " << Key.kvalue << endl); 02678 return 1; 02679 } 02680 } 02681 02682 // --- Number of bytes for each pixel/voxel of the image data 02683 // "number of bytes per pixel" / Intf_fields.nb_bytes_pixel 02684 else if(IntfCheckKeyMatch(Key, "number of bytes per pixel")) 02685 { 02686 if(ConvertFromString(Key.kvalue , &ap_IF->nb_bytes_pixel) ) 02687 { 02688 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of bytes per pixel' key. Recovered value = " << Key.kvalue << endl); 02689 return 1; 02690 } 02691 } 02692 02693 // --- Slice orientations 02694 // "slice orientation" / Intf_fields.slice_orientation 02695 // TODO : This information is recovered, but not used for the current implementation 02696 else if(IntfCheckKeyMatch(Key, "slice orientation")) 02697 { 02698 string slice_orientation; 02699 if( ConvertFromString(Key.kvalue , &slice_orientation) ) 02700 { 02701 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'slice orientation' key. Recovered value = " << Key.kvalue << endl); 02702 return 1; 02703 } 02704 02705 if (slice_orientation == "transverse") 02706 ap_IF->slice_orientation = INTF_TRANSVERSE; 02707 else if (slice_orientation == "sagittal") 02708 ap_IF->slice_orientation = INTF_SAGITTAL; 02709 else if (slice_orientation == "coronal") 02710 ap_IF->slice_orientation = INTF_CORONAL; 02711 } 02712 02713 // --- Patient rotation 02714 // "patient rotation" / Intf_fields.pat_rotation 02715 // TODO : This information is recovered, but not used for the current implementation 02716 else if(IntfCheckKeyMatch(Key, "patient rotation")) 02717 { 02718 string pat_rotation; 02719 if( ConvertFromString(Key.kvalue , &pat_rotation) ) 02720 { 02721 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'patient rotation' key. Recovered value = " << Key.kvalue << endl); 02722 return 1; 02723 } 02724 02725 if (pat_rotation == "supine") 02726 ap_IF->pat_rotation = INTF_SUPINE; 02727 else if (pat_rotation == "prone") 02728 ap_IF->pat_rotation = INTF_PRONE; 02729 } 02730 02731 02732 // --- Patient orientation 02733 // "patient orientation" / Intf_fields.pat_orientation 02734 // TODO : This information is recovered, but not used for the current implementation 02735 else if(IntfCheckKeyMatch(Key, "patient orientation")) 02736 { 02737 string pat_orientation; 02738 if( ConvertFromString(Key.kvalue , &pat_orientation) ) 02739 { 02740 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'patient orientation' key. Recovered value = " << Key.kvalue << endl); 02741 return 1; 02742 } 02743 02744 if (pat_orientation == "head_in") 02745 ap_IF->pat_orientation = INTF_HEADIN; 02746 else if (pat_orientation == "feet_in") 02747 ap_IF->pat_orientation = INTF_FEETIN; 02748 } 02749 02750 02751 // --- Multiplicative calibration value 02752 // "rescale slope" / Intf_fields.rescale_slope 02753 else if(IntfCheckKeyMatch(Key, "rescale slope")) 02754 { 02755 double rs= 1.; 02756 if( ConvertFromString(Key.kvalue , &rs) ) 02757 { 02758 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'rescale slope' key. Recovered value = " << Key.kvalue << endl); 02759 return 1; 02760 } 02761 ap_IF->rescale_slope *= rs; // factorization as this variable can also be modified by quantification units 02762 02763 if (ap_IF->rescale_slope == 0.) 02764 { 02765 Cerr("***** ReadIntfHeader()-> Error : field 'resclale slope' units should be >0!" << endl); 02766 return 1; 02767 } 02768 } 02769 02770 // Several global scale factors. 02771 // --- Multiplicative calibration value 02772 // "quantification units" / Intf_fields.rescale_slope 02773 // Note : throw warning if bad conversion, as this key is sometimes 02774 // used with a string gathering the data unit (eg : kbq/cc) 02775 else if(IntfCheckKeyMatch(Key, "quantification units")) 02776 { 02777 double qu= 1.; 02778 if(ConvertFromString(Key.kvalue , &qu) ) 02779 { 02780 Cerr("***** ReadIntfHeader()-> WARNING : Error when trying to read numeric value from 'quantification units' key. Actual value = " << Key.kvalue << endl); 02781 Cerr("***** This key will be ignored" << endl); 02782 qu= 1.; 02783 } 02784 ap_IF->rescale_slope *= qu; // factorization as this variable can also be modified by rescale slope 02785 02786 if (ap_IF->rescale_slope == 0.) 02787 { 02788 Cerr("***** ReadIntfHeader()-> Error : field 'quantification units' should be >0!" << endl); 02789 return 1; 02790 } 02791 } 02792 02793 // --- Additive calibration value 02794 // "rescale intercept" / Intf_fields.rescale_intercept 02795 else if(IntfCheckKeyMatch(Key, "rescale intercept")) 02796 { 02797 if( ConvertFromString(Key.kvalue , &ap_IF->rescale_intercept) ) 02798 { 02799 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'rescale intercept' key. Recovered value = " << Key.kvalue << endl); 02800 return 1; 02801 } 02802 02803 if (ap_IF->rescale_intercept == 0.) 02804 { 02805 Cerr("***** ReadIntfHeader()-> Error : field 'resclale intercept' units should be >0!" << endl); 02806 return 1; 02807 } 02808 } 02809 02810 02811 // --- Type of projeted/reconstructed dataset (Static|Dynamic|Gated|Tomographic|Curve|ROI|GSPECT|Other) 02812 // Curve|ROI|Other are considered as Static 02813 // "type of data" / Intf_fields.data_type 02814 else if(IntfCheckKeyMatch(Key, "type of data")) 02815 { 02816 string data_str; 02817 if(ConvertFromString(Key.kvalue ,&data_str) ) 02818 { 02819 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'type of data' key. Recovered value = " << Key.kvalue << endl); 02820 return 1; 02821 } 02822 02823 ap_IF->data_type = IntfKeyGetInputImgDataType(data_str); 02824 } 02825 02826 02827 // --- Acquisition duration 02828 // "study duration (sec)" / Intf_fields.study_duration 02829 else if(IntfCheckKeyMatch(Key, "study duration (sec)")) 02830 { 02831 if(ConvertFromString(Key.kvalue ,&ap_IF->study_duration) ) 02832 { 02833 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'study duration (sec)' key. Recovered value = " << Key.kvalue << endl); 02834 return 1; 02835 } 02836 } 02837 02838 02839 // --- Image duration duration 02840 // "image_duration (sec)" / Intf_fields.image_duration 02841 // Note : check after reading all headers that the number of image durations match the actual number of frames 02842 else if(IntfCheckKeyMatch(Key, "image duration (sec)")) 02843 { 02844 FLTNB duration; 02845 if(ConvertFromString(Key.kvalue ,&duration) ) 02846 { 02847 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'study duration (sec)' key. Recovered value = " << Key.kvalue << endl); 02848 return 1; 02849 } 02850 ap_IF->image_duration.push_back(duration); 02851 } 02852 02853 // --- Image duration duration 02854 // "image_duration (sec)" / Intf_fields.image_duration 02855 // Note : check after reading all headers that the number of image durations match the actual number of frames 02856 else if(IntfCheckKeyMatch(Key, "pause between frame groups (sec)")) 02857 { 02858 FLTNB pause_duration; 02859 if(ConvertFromString(Key.kvalue ,&pause_duration) ) 02860 { 02861 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'study duration (sec)' key. Recovered value = " << Key.kvalue << endl); 02862 return 1; 02863 } 02864 ap_IF->frame_group_pause.push_back(pause_duration); 02865 } 02866 02867 02868 // TODO : this key is not used in the current implementation 02869 // --- number of time windows 02870 // "number of time windows " / Intf_fields.nb_time_windows 02871 else if(IntfCheckKeyMatch(Key, "number of time windows ")) 02872 { 02873 if(ConvertFromString(Key.kvalue ,&ap_IF->nb_time_windows) ) 02874 { 02875 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of time windows ' key. Recovered value = " << Key.kvalue << endl); 02876 return 1; 02877 } 02878 Cerr("***** ReadIntfHeader()-> WARNING : 'number of time windows' has no use in the current implementation !"<< endl); 02879 } 02880 02881 02882 // --- Process status : acquired or reconstructed 02883 // "process status" / Intf_fields.process_status 02884 else if(IntfCheckKeyMatch(Key, "process status")) 02885 { 02886 if(ConvertFromString(Key.kvalue ,&ap_IF->process_status) ) 02887 { 02888 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'process status' key. Recovered value = " << Key.kvalue << endl); 02889 return 1; 02890 } 02891 } 02892 02893 02894 // --- Number of energy windows 02895 // "number of energy windows" / Intf_fields.nb_energy_windows 02896 // TODO : not implemented yet. 02897 else if(IntfCheckKeyMatch(Key, "number of energy windows")) 02898 { 02899 if(ConvertFromString(Key.kvalue , &ap_IF->nb_energy_windows) ) 02900 { 02901 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of energy windows' key. Recovered value = " << Key.kvalue << endl); 02902 return 1; 02903 } 02904 //Cerr("***** ReadIntfHeader()-> WARNING : 'number of energy windows' has no use in the current implementation !"<< endl); 02905 } 02906 02907 02908 // --- Number of detector heads (SPECT systems) 02909 // "number of detector heads" / Intf_fields.nb_detector_heads 02910 else if(IntfCheckKeyMatch(Key, "number of detector heads")) 02911 { 02912 if(ConvertFromString(Key.kvalue , &ap_IF->nb_detector_heads) ) 02913 { 02914 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of detector heads' key. Recovered value = " << Key.kvalue << endl); 02915 return 1; 02916 } 02917 } 02918 02919 02920 // --- Number of projections (for one head in SPECT) 02921 // "number of projections" / Intf_fields.nb_energy_windows 02922 else if(IntfCheckKeyMatch(Key, "number of projections")) 02923 { 02924 if(ConvertFromString(Key.kvalue , &ap_IF->nb_projections) ) 02925 { 02926 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of projections' key. Recovered value = " << Key.kvalue << endl); 02927 return 1; 02928 } 02929 } 02930 02931 02932 // --- Angular span of the projections ex: 180, 360 02933 // "extent of rotation " / Intf_fields.extent_rotation 02934 else if(IntfCheckKeyMatch(Key, "extent of rotation")) 02935 { 02936 if(ConvertFromString(Key.kvalue , &ap_IF->extent_rotation) ) 02937 { 02938 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'extent of rotation' key. Recovered value = " << Key.kvalue << endl); 02939 return 1; 02940 } 02941 } 02942 02943 // --- Direction of rotation : CCW (counter-clockwise), CW (clockwise) 02944 // "direction of rotation" / Intf_fields.direction_rotation 02945 else if(IntfCheckKeyMatch(Key, "direction of rotation")) 02946 { 02947 if(ConvertFromString(Key.kvalue , &ap_IF->direction_rotation) ) 02948 { 02949 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'direction_rotation' key. Recovered value = " << Key.kvalue << endl); 02950 return 1; 02951 } 02952 } 02953 02954 // --- Angle of the first view 02955 // "start angle" / Intf_fields.first_angle 02956 // TODO : This may be wrong. First angle may be sum of values of : 02957 // "start angle" + "first projection angle in data set" 02958 else if(IntfCheckKeyMatch(Key, "start angle")) 02959 { 02960 if(ConvertFromString(Key.kvalue , &ap_IF->first_angle) ) 02961 { 02962 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'start angle' key. Recovered value = " << Key.kvalue << endl); 02963 return 1; 02964 } 02965 } 02966 02967 02968 // --- All projection angles as an array key 02969 // "projection_angles" / Intf_fields.projection_angles 02970 else if(IntfCheckKeyMatch(Key, "projection_angles")) 02971 { 02972 if(ConvertFromString(Key.kvalue , &ap_IF->projection_angles) ) 02973 { 02974 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'projection_angles' key. Recovered value = " << Key.kvalue << endl); 02975 return 1; 02976 } 02977 } 02978 02979 02980 // --- All distance between center of rotation and detectors, as unique value similar for each angles 02981 // "Center of rotation to detector distance" / Intf_fields.radius 02982 else if(IntfCheckKeyMatch(Key, "Center of rotation to detector distance")) 02983 { 02984 if(ConvertFromString(Key.kvalue , &ap_IF->radius) ) 02985 { 02986 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'Center of rotation to detector distance' key."); 02987 Cerr(" Recovered value = " << Key.kvalue << endl); 02988 return 1; 02989 } 02990 } 02991 02992 02993 // --- All distance between center of rotation and detectors, as an array key 02994 // "Radius" / Intf_fields.radius 02995 else if(IntfCheckKeyMatch(Key, "Radius")) 02996 { 02997 if(ConvertFromString(Key.kvalue , &ap_IF->radius) ) 02998 { 02999 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'Radius' key. Recovered value = " << Key.kvalue << endl); 03000 return 1; 03001 } 03002 } 03003 03004 // Return error if data compression 03005 else if(IntfCheckKeyMatch(Key, "data compression")) 03006 { 03007 if(! (Key.kvalue == "none" || Key.kvalue.empty() ) ) 03008 { 03009 Cerr("***** ReadIntfHeader()-> Error : compressed interfile images not handled by the current implementation !" << endl); 03010 return 1; 03011 } 03012 } 03013 03014 // Return error if data encoded 03015 else if(IntfCheckKeyMatch(Key, "data encode")) 03016 { 03017 if(! (Key.kvalue == "none" || Key.kvalue.empty() ) ) 03018 { 03019 Cerr("***** ReadIntfHeader()-> Error : encoded interfile images not handled by the current implementation !" << endl); 03020 return 1; 03021 } 03022 } 03023 03024 03025 if(IntfCheckKeyMatch(Key, "number of slices") ) 03026 { 03027 if(!(ap_IF->mtx_size[2] >0)) // Check if not been initialized elsewhere using other field (matrix size [3]) 03028 { 03029 if(ConvertFromString(Key.kvalue , &ap_IF->mtx_size[2]) ) 03030 { 03031 Cerr("***** ReadIntfHeader()-> Error when trying to read data from 'number of slices' key. Recovered value = " << Key.kvalue << endl); 03032 return 1; 03033 } 03034 ap_IF->nb_dims++; 03035 } 03036 continue; 03037 } 03038 03039 } 03040 } 03041 else 03042 { 03043 Cerr("***** ReadIntfHeader()-> Error : couldn't find or read header interfile '"<< a_pathToHeaderFile << "' !" << endl); 03044 return 1; 03045 } 03046 03047 // Recover slice number in ap_IF->mtx_size[2] if provided by "total number of images" or "number of images" tags, and not "matrix size [3]" tag. 03048 if(ap_IF->mtx_size[2]<0 && ap_IF->nb_total_imgs>0) 03049 ap_IF->mtx_size[2] = ap_IF->nb_total_imgs; 03050 03051 03052 // Compute slice thickness if not initialized 03053 ap_IF->vox_size[2] = (ap_IF->vox_size[2] < 0) ? 03054 ((ap_IF->vox_size[0]+ap_IF->vox_size[1])/2.) * ap_IF->slice_thickness_mm: 03055 ap_IF->vox_size[2]; 03056 03057 // Compute nb dimensions for the image 03058 ap_IF->nb_dims = 3; 03059 03060 if(ap_IF->mtx_size[3] >1 || 03061 ap_IF->nb_time_frames>1) 03062 ap_IF->nb_dims++; 03063 if( ap_IF->nb_resp_gates >1) 03064 ap_IF->nb_dims++; 03065 if( ap_IF->nb_card_gates >1) 03066 ap_IF->nb_dims++; 03067 03068 // If there are any frames, check nb frame durations matches nb (frames) group pauses 03069 if(ap_IF->image_duration.size() > 0 && 03070 ap_IF->frame_group_pause.size() > 0 && 03071 ap_IF->image_duration.size() != ap_IF->frame_group_pause.size()) 03072 { 03073 Cerr("***** ReadIntfHeader()-> Error : nb of recovered frame duration ('"<< ap_IF->image_duration.size() 03074 << ") does not match the nb of recovered pauses between frames ('"<< ap_IF->frame_group_pause.size() << ") !" << endl); 03075 return 1; 03076 } 03077 03078 return 0; 03079 } 03080 03081 03082 03083 03084 // ===================================================================== 03085 // --------------------------------------------------------------------- 03086 // --------------------------------------------------------------------- 03087 // ===================================================================== 03088 /* 03089 \fn IntfWriteHeaderMainData 03090 \param a_path : path to the interfile header to write 03091 \param ap_IntfF : Structure containing interfile keys of the image to write 03092 \param vb : verbosity 03093 \brief Write main infos of an Interfile header at the path provided in parameter, 03094 using the interfile keys provided by the field structure parameter. 03095 \todo Recover in 'patient name' the name of the datafile used to reconstruct the images [ bed ]? 03096 \return 0 if success, positive value otherwise. 03097 */ 03098 int IntfWriteHeaderMainData(const string& a_path, const Intf_fields& ap_IntfF, int vb) 03099 { 03100 if (vb >= 4) 03101 { 03102 Cout("--------------------------------------------------------------- " << endl); 03103 Cout("IntfWriteHeaderMainData()-> Start writing header interfile " << a_path << endl); 03104 Cout("--------------------------------------------------------------- " << endl); 03105 } 03106 03107 string path_to_header, path_to_image; 03108 path_to_header = a_path; 03109 03110 path_to_header.append(".hdr"); 03111 path_to_image = GetFileFromPath(a_path).append(".img"); 03112 03113 // Get file and check existence 03114 ofstream ofile(path_to_header.c_str(), ios::out); 03115 ofile << setprecision(std::numeric_limits<FLTNB>::digits10 +1); 03116 03117 if(ofile) 03118 { 03119 sScannerManager* p_scanMgr = sScannerManager::GetInstance(); 03120 03121 ofile << "!INTERFILE := " << endl; 03122 03123 // PET/SPECT/CT according to the scanner type defined in scanner manager 03124 string imaging_modality; 03125 if (p_scanMgr->GetScannerType() == SCANNER_PET) 03126 imaging_modality = "PET"; 03127 else if (p_scanMgr->GetScannerType() == SCANNER_SPECT_CONVERGENT || p_scanMgr->GetScannerType() == SCANNER_SPECT_PINHOLE) 03128 imaging_modality = "SPECT"; 03129 else if(p_scanMgr->GetScannerType() == SCANNER_CT) 03130 imaging_modality = "CT"; 03131 else 03132 imaging_modality = "UNKOWN"; 03133 03134 ofile << "!imaging modality := " << imaging_modality << endl; 03135 03136 // 3D SPECT use 3.3 03137 if(ap_IntfF.data_type == INTF_IMG_SPECT) 03138 ofile << "!version of keys := " << "3.3" << endl; 03139 else 03140 ofile << "!version of keys := " << "CASToRv1.0" << endl; 03141 03142 ofile << endl << "!GENERAL DATA := " << endl; 03143 ofile << "!originating system := " << p_scanMgr->GetScannerName() << endl; 03144 ofile << "!data offset in bytes := " << 0 << endl; 03145 ofile << "!name of data file := " << path_to_image << endl; 03146 //ofile << "patient name:= " << "unknown" << endl; // TODO Recover here the name of the datafile used to reconstruct the images [ bed ]? 03147 03148 ofile << endl << "!GENERAL IMAGE DATA " << endl; 03149 03150 string data_type_str = (ap_IntfF.data_type == INTF_IMG_SPECT) ? 03151 "Tomographic" : 03152 "Static"; 03153 ofile << "!type of data := " << data_type_str << endl; // such as Static|Dynamic|Gated|Tomographic|Curve|ROI|GSPECT|Other 03154 03155 uint32_t nb_imgs = (ap_IntfF.process_status == "acquired") ? 03156 ap_IntfF.nb_projections : // projeted data 03157 ap_IntfF.nb_total_imgs ; // reconstructed data 03158 ofile << "!total number of images := " << nb_imgs << endl; 03159 03160 ofile << "imagedata byte order := " << IntfKeyGetEndianStr(ap_IntfF.endianness) << endl; 03161 03162 // Depending on the mode of output writing, we should write a header for each dynamic image, or append everything in one single header 03163 // output dynamic files as one single file instead of 3D images with a meta header 03164 if( ap_IntfF.nb_time_frames>1 ) 03165 ofile << "number of time frames := " << ap_IntfF.nb_time_frames << endl; 03166 if( ap_IntfF.nb_energy_windows>1 ) 03167 ofile << "number of time windows := " << ap_IntfF.nb_resp_gates * 03168 ap_IntfF.nb_card_gates << endl; 03169 if(ap_IntfF.nb_resp_gates > 1) 03170 ofile << "number of respiratory gates := " << ap_IntfF.nb_resp_gates << endl; 03171 if(ap_IntfF.nb_card_gates > 1) 03172 ofile << "number of cardiac gates := " << ap_IntfF.nb_card_gates << endl; 03173 03174 if(ap_IntfF.process_status == "acquired") 03175 { 03176 ofile << "!number of projections := " << ap_IntfF.nb_projections << endl; 03177 ofile << "!extent of rotation := " << ap_IntfF.extent_rotation << endl; 03178 } 03179 ofile << "process status := " << ap_IntfF.process_status << endl; 03180 03181 03182 if (ap_IntfF.data_type == INTF_IMG_DYNAMIC) // Dynamic study 03183 ofile << endl << "!DYNAMIC STUDY (General) :=" << endl; 03184 else if(p_scanMgr->GetScannerType() == 2) // SPECT study 03185 { 03186 if( ap_IntfF.data_type == INTF_IMG_GATED) 03187 ofile << endl << "!GSPECT STUDY (General) :=" << endl; 03188 else 03189 { 03190 ofile << endl << "!SPECT STUDY (General) :=" << endl; 03191 ofile << endl << "!SPECT STUDY ("<< ap_IntfF.process_status <<" data) :="<< endl; 03192 } 03193 } 03194 else if (ap_IntfF.data_type == INTF_IMG_GATED ) 03195 ofile << endl << "!GATED STUDY (General) :=" << endl; 03196 else // Standard study 03197 ofile << endl << "!STATIC STUDY (General) :=" << endl; 03198 03199 03200 // Patient position 03201 // (not implemented as these informations are most of the time missing, and could be misleading) 03202 // output_header << "patient orientation := " << ap_IntfF.pat_orientation << endl; 03203 // output_header << "patient rotation := " << ap_IntfF.pat_rotation << endl; 03204 // output_header << "slice orientation := " << ap_IntfF.slice_orientation << endl; 03205 03206 03207 // Loop on dynamic images (if any) and write data 03208 03209 if (ap_IntfF.nb_resp_gates>1) 03210 ofile << "!number of frame groups :=" << ap_IntfF.nb_time_frames << endl; 03211 03212 // loop on frames 03213 for(int fr=0 ; fr<ap_IntfF.nb_time_frames ; fr++) 03214 { 03215 if( ap_IntfF.nb_time_frames>1 ) 03216 { 03217 ofile << "!Dynamic Study (each frame group) :=" << endl; 03218 ofile << "!frame group number := " << fr+1 << endl; 03219 } 03220 03221 // loop on respiratory gates 03222 for(int rg=0 ; rg<ap_IntfF.nb_resp_gates ; rg++) 03223 { 03224 if( ap_IntfF.nb_resp_gates>1 ) 03225 { 03226 ofile << "!Respiratory Gated Study (each time window) :=" << endl; 03227 ofile << "!time window number := " << rg+1 << endl; 03228 } 03229 03230 // loop on cardiac gates 03231 for(int cg=0 ; cg<ap_IntfF.nb_card_gates ; cg++) 03232 { 03233 if( ap_IntfF.nb_card_gates>1 ) 03234 { 03235 ofile << "!Cardiac Gated Study (each time window) :=" << endl; 03236 ofile << "!time window number := " << cg+1 << endl; 03237 } 03238 03239 // Write image specific data 03240 if(IntfWriteHeaderImgData(ofile, ap_IntfF, vb) ) 03241 { 03242 Cerr("***** IntfWriteHeaderMainData()-> Error : while trying to write the interfile header '"<< path_to_header << "' !" << endl); 03243 return 1; 03244 } 03245 03246 if( ap_IntfF.nb_card_gates>1 ) 03247 ofile << "!number of images in this time window :=" 03248 << ap_IntfF.mtx_size[2] << endl; 03249 03250 } 03251 03252 if( ap_IntfF.nb_resp_gates>1 ) 03253 ofile << "!number of images in this time window :=" 03254 << ap_IntfF.nb_card_gates * 03255 ap_IntfF.mtx_size[2] << endl; 03256 03257 } 03258 03259 if (ap_IntfF.nb_time_frames>1) 03260 { 03261 ofile << "!number of images this frame group :=" 03262 << ap_IntfF.nb_resp_gates * 03263 ap_IntfF.nb_card_gates * 03264 ap_IntfF.mtx_size[2] << endl; 03265 03266 ofile << "!image duration (sec) := " << ap_IntfF.image_duration[fr] << endl; 03267 if(fr+1 == ap_IntfF.nb_time_frames) // Check if we reached the last frame 03268 ofile << "pause between frame groups (sec) " << 0.0 << endl; 03269 else 03270 ofile << "pause between frame groups (sec) " << ap_IntfF.frame_group_pause[fr] << endl; 03271 } 03272 } 03273 03274 ofile << "!END OF INTERFILE := " << endl; 03275 } 03276 else 03277 { 03278 Cerr("***** IntfWriteHeaderMainData()-> Error : couldn't find output header interfile '"<< a_path << "' !" << endl); 03279 return 1; 03280 } 03281 03282 return 0; 03283 } 03284 03285 03286 03287 03288 // ===================================================================== 03289 // --------------------------------------------------------------------- 03290 // --------------------------------------------------------------------- 03291 // ===================================================================== 03292 /* 03293 \fn IntfWriteHeaderImgData 03294 \param ap_ofile : pointer to the ofstream linked to the image header file 03295 \param ap_IntfF : reference to a structure containing interfile keys of the image to write 03296 \param verbosity : verbosity 03297 \brief Write basic image data info in the file provided in parameter 03298 \todo Data about positionning/offsets ? 03299 \todo : include dectector heads (a set of fields for each head) 03300 \return 0 if success, positive value otherwise. 03301 */ 03302 int IntfWriteHeaderImgData(ofstream &ap_ofile, const Intf_fields& ap_IntfF, int vb) 03303 { 03304 if(ap_ofile) 03305 { 03306 if(ap_IntfF.process_status == "acquired") // projection data (SPECT) 03307 { 03308 // Todo : include dectector heads (a set of fields for each head) 03309 // That means a lot of modifications as angles, radius, ... should be provided for each head 03310 ap_ofile << "!direction of rotation := " << ap_IntfF.direction_rotation << endl; 03311 ap_ofile << "start angle := " << ap_IntfF.first_angle << endl; 03312 ap_ofile << "projection angles := " << ap_IntfF.projection_angles << endl; 03313 03314 // If one common radius for each projection (no "{}"), write it in the 'Radius' key 03315 if( ap_IntfF.radius.find("{}") != string::npos ) 03316 ap_ofile << "Center of rotation to detector distance := " << ap_IntfF.radius << endl; 03317 else 03318 ap_ofile << "Radius := " << ap_IntfF.radius << endl; 03319 03320 ap_ofile << "!matrix size [1] := " << ap_IntfF.mtx_size[0] << endl; 03321 ap_ofile << "!matrix size [2] := " << ap_IntfF.mtx_size[1] << endl; 03322 ap_ofile << "!number format := " << ap_IntfF.nb_format << endl; 03323 ap_ofile << "!number of bytes per pixel := " << sizeof(FLTNB) << endl; 03324 ap_ofile << "scaling factor (mm/pixel) [1] := " << ap_IntfF.vox_size[0] << endl; 03325 ap_ofile << "scaling factor (mm/pixel) [2] := " << ap_IntfF.vox_size[1] << endl; 03326 ap_ofile << "!data offset in bytes := " << 0 << endl; 03327 } 03328 else 03329 { 03330 ap_ofile << "number of dimensions := " << 3 << endl; 03331 ap_ofile << "!matrix size [1] := " << ap_IntfF.mtx_size[0] << endl; 03332 ap_ofile << "!matrix size [2] := " << ap_IntfF.mtx_size[1] << endl; 03333 ap_ofile << "!matrix size [3] := " << ap_IntfF.mtx_size[2] << endl; 03334 ap_ofile << "!number format := " << ap_IntfF.nb_format << endl; 03335 ap_ofile << "!number of bytes per pixel := " << sizeof(FLTNB) << endl; 03336 ap_ofile << "scaling factor (mm/pixel) [1] := " << ap_IntfF.vox_size[0] << endl; 03337 ap_ofile << "scaling factor (mm/pixel) [2] := " << ap_IntfF.vox_size[1] << endl; 03338 ap_ofile << "scaling factor (mm/pixel) [3] := " << ap_IntfF.vox_size[2] << endl; 03339 03340 if(ap_IntfF.rescale_intercept != 1 || ap_IntfF.rescale_slope != 0) 03341 { 03342 ap_ofile << "data rescale offset := " << ap_IntfF.rescale_intercept << endl; 03343 ap_ofile << "data rescale slope := " << ap_IntfF.rescale_slope << endl; 03344 } 03345 // Todo, we may need keys for image positionning in space 03346 //ap_ofile << "origin (mm) [1] := " << ap_ID->Get??? << endl; 03347 //ap_ofile << "origin (mm) [2] := " << ap_ID->Get??? << endl; 03348 //ap_ofile << "origin (mm) [3] := " << ap_ID->Get??? << endl; 03349 03350 03351 // The following keys were used in interfiles format 3.3 03352 // Seem to have been replaced by matrix size [3]/scaling factor (mm/pixel) [3] 03353 // but we may need them for some compatibilities issues 03354 // ap_ofile << "!number of slices := " << ap_IntfF.mtx_size[2] << endl; 03355 // ap_ofile << "slice thickness (pixels) := " << 03356 // ap_IntfF.vox_size[2] /( (ap_IntfF.vox_size[0]+ap_IntfF.vox_size[1])/2.) << endl; 03357 // ap_ofile << "centre-centre slice separation (pixels) := " << endl; 03358 } 03359 } 03360 else 03361 { 03362 Cerr("***** IntfWriteHeaderImgData()-> Error : couldn't open provided interfile header file !" << endl); 03363 return 1; 03364 } 03365 03366 return 0; 03367 } 03368 03369 03370 03371 03372 // ===================================================================== 03373 // --------------------------------------------------------------------- 03374 // --------------------------------------------------------------------- 03375 // ===================================================================== 03376 /* 03377 \fn IntfWriteImage() 03378 \param a_pathToImg : th to the directory where the image will be written 03379 \param ap_outImgMtx : Matrix containing the image to write 03380 \param a_dim : Number of voxels in the 3D image 03381 \param vb : Verbosity level 03382 \brief Write Interfile raw data whose path is provided in parameter, using image matrix provided in parameter. 03383 \brief For 1 dimensional image matrices 03384 \return 0 if success, positive value otherwise. 03385 */ 03386 int IntfWriteImage(const string& a_pathToImg, FLTNB* ap_outImgMtx, uint32_t a_dim, int vb) 03387 { 03388 if(vb >= 5) Cout("IntfWriteImage()*" << endl); 03389 03390 ofstream img_file(a_pathToImg.c_str(), ios::binary | ios::out); 03391 03392 if(img_file) 03393 { 03394 // Call to writing function. 03395 if(IntfWriteData(&img_file, ap_outImgMtx, a_dim, vb) ) 03396 { 03397 Cerr("***** IntfWriteImage()-> Error occurred when writing the image file '"<< a_pathToImg << "' !" << endl); 03398 return 1; 03399 } 03400 } 03401 else 03402 { 03403 Cerr("***** IntfWriteImage()-> Error occurred while trying to write the image file '"<< a_pathToImg << "' !" << endl); 03404 return 1; 03405 } 03406 03407 return 0; 03408 } 03409 03410 03411 03412 03413 // ===================================================================== 03414 // --------------------------------------------------------------------- 03415 // --------------------------------------------------------------------- 03416 // ===================================================================== 03417 /* 03418 \fn IntfWriteImage 03419 \param ap_pathToImgs: List of string containing the paths to the image to write 03420 \param a2p_outImgMtx : 2 dimensional image matrix (1D temporal + 3D) 03421 \param ap_imgDim[2] : Dimensions of ap_outImgMtx 03422 \param vb : Verbosity level 03423 \brief Write Interfile raw data whose path is provided in parameter, using the image matrix provided in parameter. 03424 \brief For 2 dimensional image matrices 03425 \return 0 if success, positive value otherwise. 03426 */ 03427 int IntfWriteImage(vector<string> ap_pathToImgs, FLTNB** a2p_outImgMtx, uint32_t ap_dim[2], int vb) 03428 { 03429 if(vb >= 5) Cout("IntfWriteImage()**" << endl); 03430 03431 if(ap_pathToImgs.size() == 1) // We have one single data file containing all the data 03432 { 03433 // As the content will be appended, make sure there is no existing file with such name 03434 // remove it otherwise 03435 ifstream fcheck(ap_pathToImgs.at(0).append(".img").c_str()); 03436 if(fcheck.good()) 03437 { 03438 fcheck.close(); 03439 #ifdef _WIN32 03440 string dos_instruction = "del " + ap_pathToImgs.at(0); 03441 system(dos_instruction.c_str()); 03442 #else 03443 remove(ap_pathToImgs.at(0).c_str()); 03444 #endif 03445 } 03446 03447 // Add app to the options as we will write the image in several step. 03448 ofstream img_file(ap_pathToImgs.at(0).c_str(), ios::binary | ios::out | ios::app); 03449 03450 if(img_file) 03451 { 03452 // Call to writing function for each image matrix 03453 for(uint32_t d1=0 ; d1<ap_dim[0] ; d1++) 03454 if( IntfWriteData(&img_file, a2p_outImgMtx[d1], ap_dim[1], vb) ) 03455 { 03456 Cerr("***** IntfWriteImage()-> Error occurred when writing the image file '"<< ap_pathToImgs.at(d1) << "' !" << endl); 03457 return 1; 03458 } 03459 } 03460 else 03461 { 03462 Cerr("***** IntfWriteImage()-> Error occurred while trying to write the image file at the path: '"<< ap_pathToImgs.at(0) << "' !" << endl); 03463 return 1; 03464 } 03465 } 03466 else // We have a number of data file equal to the number of image to load 03467 { 03468 if(! ap_pathToImgs.size()!=ap_dim[0]) // Check we have a number of file corresponding to the number of images to load 03469 { 03470 Cerr("***** IntfWriteImage()-> Error : number of interfile images ("<< ap_pathToImgs.size() << ") not consistent with the number of images to load (" << ap_dim[0] << ") !" << endl); 03471 return 1; 03472 } 03473 03474 for(uint32_t d1=0 ; d1<ap_dim[0] ; d1++) 03475 { 03476 ofstream img_file(ap_pathToImgs.at(d1).c_str(), ios::binary | ios::out); // Should be one different path by file 03477 03478 if(img_file) 03479 { 03480 // Call to writing function. 03481 if(IntfWriteData(&img_file, a2p_outImgMtx[d1], ap_dim[1], vb) ) 03482 { 03483 Cerr("***** IntfWriteImage()-> Error occurred when writing the image file '"<< ap_pathToImgs.at(d1) << "' !" << endl); 03484 return 1; 03485 } 03486 } 03487 else 03488 { 03489 Cerr("***** IntfWriteImage()-> Error occurred while trying to write the image file at the path: '"<< ap_pathToImgs.at(d1) << "' !" << endl); 03490 return 1; 03491 } 03492 } 03493 } 03494 03495 return 0; 03496 } 03497 03498 03499 03500 03501 // ===================================================================== 03502 // --------------------------------------------------------------------- 03503 // --------------------------------------------------------------------- 03504 // ===================================================================== 03505 /* 03506 \fn IntfWriteImage() 03507 \param vector<string> ap_pathToImgs: List of string containing the paths to the image to write 03508 \param FLTNB**** a4p_outImgMtx : 4 dimensional image matrix (3D temporal + 3D) 03509 \param int ap_imgDim[4] : Dimensions of ap_outImgMtx 03510 \param int vb : Verbosity level 03511 \brief Write Interfile raw data whose path is provided in parameter, using the image matrix provided in parameter. 03512 \brief For 4 dimensional image matrices 03513 \return 0 if success, positive value otherwise. 03514 */ 03515 int IntfWriteImage(vector<string> ap_pathToImgs, FLTNB**** a4p_outImgMtx, uint32_t ap_dim[4], int vb) 03516 { 03517 if(vb >= 5) Cout("IntfWriteImage()" << endl); 03518 03519 if(ap_pathToImgs.size() == 1) // We have one single data file containing all the data 03520 { 03521 // As the content will be appended, make sure there is no existing file with such name 03522 // remove it otherwise 03523 ifstream fcheck(ap_pathToImgs.at(0).append(".img").c_str()); 03524 if(fcheck.good()) 03525 { 03526 fcheck.close(); 03527 #ifdef _WIN32 03528 string dos_instruction = "del " + ap_pathToImgs.at(0); 03529 system(dos_instruction.c_str()); 03530 #else 03531 remove(ap_pathToImgs.at(0).c_str()); 03532 #endif 03533 } 03534 03535 // Add app to the options as we will write the image in several step. 03536 ofstream img_file(ap_pathToImgs.at(0).c_str(), ios::binary | ios::out | ios::app); 03537 03538 if(img_file) 03539 { 03540 // Call to writing function for each image matrix 03541 for(uint32_t d1=0 ; d1<ap_dim[0] ; d1++) 03542 for(uint32_t d2=0 ; d2<ap_dim[1] ; d2++) 03543 for(uint32_t d3=0 ; d3<ap_dim[2] ; d3++) 03544 if(IntfWriteData(&img_file, a4p_outImgMtx[d1][d2][d3], ap_dim[3], vb) ) 03545 { 03546 int idx_img = d1*ap_dim[1]*ap_dim[2] + d2*ap_dim[2] + d3; 03547 Cerr("***** IntfWriteImage()-> Error occurred when writing the image file '"<< ap_pathToImgs.at(idx_img) << "' !" << endl); 03548 return 1; 03549 } 03550 } 03551 else 03552 { 03553 Cerr("***** IntfWriteImage()-> Error occurred while trying to write the image file '"<< ap_pathToImgs.at(0) << "' !" << endl); 03554 return 1; 03555 } 03556 03557 } 03558 else // We have a number of data file equal to the number of image to load 03559 { 03560 // Check we have a number of file corresponding to the number of images to load 03561 if(ap_pathToImgs.size() != ap_dim[0]*ap_dim[1]*ap_dim[2]) 03562 { 03563 Cerr("***** IntfWriteImage()-> Error : number of interfile images (="<< ap_pathToImgs.size() 03564 << ") not consistent with the number of images to load (=" 03565 << ap_dim[0]*ap_dim[1]*ap_dim[2] << ") !" << endl); 03566 return 1; 03567 } 03568 03569 for(uint32_t d1=0 ; d1<ap_dim[0] ; d1++) 03570 for(uint32_t d2=0 ; d2<ap_dim[1] ; d2++) 03571 for(uint32_t d3=0 ; d3<ap_dim[2] ; d3++) 03572 { 03573 int idx_img = d1*ap_dim[1]*ap_dim[2] + d2*ap_dim[2] + d3; 03574 ofstream img_file(ap_pathToImgs.at(idx_img).append(".img").c_str(), ios::binary | ios::out); // Should be one different path by file 03575 03576 if(img_file) 03577 { 03578 // Call to writing function. 03579 if(IntfWriteData(&img_file, a4p_outImgMtx[d1][d2][d3], ap_dim[3], vb) ) 03580 { 03581 Cerr("***** IntfWriteImage()-> Error occurred when writing the image file '"<< ap_pathToImgs.at(idx_img) << "' !" << endl); 03582 return 1; 03583 } 03584 } 03585 else 03586 { 03587 Cerr("***** IntfWriteImage()-> Error occurred while trying to write the image file at the path: '"<< ap_pathToImgs.at(idx_img) << "' !" << endl); 03588 return 1; 03589 } 03590 } 03591 } 03592 03593 return 0; 03594 } 03595 03596 03597 03598 03599 // ===================================================================== 03600 // --------------------------------------------------------------------- 03601 // --------------------------------------------------------------------- 03602 // ===================================================================== 03603 /* 03604 \fn IntfWriteData 03605 \param ap_oFile : Ofstream pointing to an image file 03606 \param ap_outImgMtx : 3D image matrix with reconstruction dimensions/voxel size containing the image data 03607 \param a_nbVox : A number of voxels in the 3D image matrix with reconstruction dimensions/voxel size 03608 \param vb : Verbosity level 03609 \brief Write the content of the image matrix in the file pointed by ofstream 03610 \todo keep original orientations ? Would require a loop on voxels and a reindexing of voxel indexes 03611 \todo check writing error 03612 \return 0 if success, positive value otherwise. 03613 */ 03614 int IntfWriteData(ofstream* ap_oFile, FLTNB* ap_outImgMatrix, int a_nbVox, int vb) 03615 { 03616 if(vb >= 5) Cout("IntfWriteData() " << endl); 03617 03618 // TODO : Keep original orientations ? Would require a loop on voxels and a reindexing of voxel indexes 03619 // TODO : Perhaps check here if nothing wrong happened during the writing (deletion of objects or writing to full disk). Not sure about how this is done for ofstream 03620 if (ap_oFile->write(reinterpret_cast<char*>(ap_outImgMatrix), a_nbVox*sizeof(FLTNB)) ) // implement error check here 03621 { 03622 //Cerr("***** IntfWriteData()-> Error occurred when trying to write the image file !" << endl); 03623 //return 1; 03624 } 03625 03626 return 0; 03627 } 03628 03629 03630 03631 03632 // ===================================================================== 03633 // --------------------------------------------------------------------- 03634 // --------------------------------------------------------------------- 03635 // ===================================================================== 03636 /* 03637 \fn IntfKeyInitFields 03638 \param ap_IF : Structure containing Interfile keys 03639 \brief Init the keys of an the Interfile structure passed in parameter to default values 03640 */ 03641 void IntfKeyInitFields(Intf_fields* ap_IF) 03642 { 03643 ap_IF->path_to_image = ""; 03644 ap_IF->endianness = User_Endianness; 03645 ap_IF->data_offset = 0; 03646 ap_IF->nb_format = IntfKeyGetPixTypeStr(); 03647 ap_IF->nb_dims = 0; 03648 for(int d=0 ; d<7 ; d++) 03649 ap_IF->mtx_size[d] = 0; 03650 for(int d=0 ; d<3 ; d++) 03651 ap_IF->vox_size[d] = -1.; 03652 03653 ap_IF->slice_thickness_mm = -1.; 03654 ap_IF->ctr_to_ctr_separation = 0; 03655 ap_IF->nb_time_frames = 1; 03656 ap_IF->nb_resp_gates = 1; 03657 ap_IF->nb_card_gates = 1; 03658 ap_IF->nb_total_imgs = 1; 03659 ap_IF->nb_bytes_pixel = 0; 03660 ap_IF->slice_orientation = 0; 03661 ap_IF->pat_rotation = 0; 03662 ap_IF->pat_orientation = 0; 03663 ap_IF->rescale_slope = 1.; 03664 ap_IF->rescale_intercept = 0.; 03665 03666 for(int d=0 ; d<3 ; d++) 03667 { 03668 ap_IF->cmtx_size[d] = 0; 03669 ap_IF->cvox_size[d] = -1; 03670 } 03671 ap_IF->is_mtx_size_different = false; 03672 ap_IF->data_type = -1; 03673 ap_IF->study_duration = -1.; 03674 // ap_IF->image_duration = vector<float>; 03675 // ap_IF->frame_group_pause = vector<float>; 03676 ap_IF->nb_time_windows = 0; 03677 ap_IF->process_status = "reconstructed"; 03678 03679 ap_IF->nb_detector_heads = 0; 03680 ap_IF->nb_energy_windows = 0; 03681 ap_IF->nb_projections = 1; 03682 ap_IF->extent_rotation = -1.; 03683 ap_IF->direction_rotation = "CW"; 03684 ap_IF->first_angle = -1.; 03685 ap_IF->projection_angles = ""; 03686 ap_IF->radius = ""; 03687 }; 03688 03689 03690 03691 03692 // ===================================================================== 03693 // --------------------------------------------------------------------- 03694 // --------------------------------------------------------------------- 03695 // ===================================================================== 03696 /* 03697 \fn IntfKeySetFieldsOutput 03698 \param ap_IF : Structure containing Interfile keys 03699 \param ap_ID : oImageDimensionsAndQuantification object containing additional infos about reconstruction 03700 \brief Init the keys of the Interfile header of an image to be written on disk 03701 \details Init the keys of the Interfile structure passed in parameter for output writing \n 03702 using the ImageDimensions object containing information about the reconstruction 03703 */ 03704 void IntfKeySetFieldsOutput(Intf_fields* ap_IF, oImageDimensionsAndQuantification* ap_ID) 03705 { 03706 // Set header metadata using Image Dimensions object 03707 ap_IF->endianness = User_Endianness; 03708 ap_IF->nb_format = IntfKeyGetPixTypeStr(); 03709 ap_IF->nb_dims = 3; 03710 if(ap_ID->GetNbTimeFrames()>1) ap_IF->nb_dims++; 03711 if(ap_ID->GetNbRespGates()>1) ap_IF->nb_dims++; 03712 if(ap_ID->GetNbCardGates()>1) ap_IF->nb_dims++; 03713 ap_IF->mtx_size[0] = ap_ID->GetNbVoxX(); 03714 ap_IF->mtx_size[1] = ap_ID->GetNbVoxY(); 03715 ap_IF->mtx_size[2] = ap_ID->GetNbVoxZ(); 03716 ap_IF->vox_size[0] = ap_ID->GetVoxSizeX(); 03717 ap_IF->vox_size[1] = ap_ID->GetVoxSizeY(); 03718 ap_IF->vox_size[2] = ap_ID->GetVoxSizeZ(); 03719 03720 ap_IF->slice_thickness_mm = ap_ID->GetVoxSizeZ(); 03721 ap_IF->ctr_to_ctr_separation = 0; 03722 ap_IF->nb_time_frames = ap_ID->GetNbTimeFrames(); 03723 ap_IF->nb_resp_gates = ap_ID->GetNbRespGates(); 03724 ap_IF->nb_card_gates = ap_ID->GetNbCardGates(); 03725 ap_IF->nb_total_imgs = ap_IF->nb_time_frames * 03726 ap_IF->nb_resp_gates * 03727 ap_IF->nb_card_gates* 03728 ap_IF->mtx_size[2]; 03729 03730 ap_IF->nb_bytes_pixel = sizeof(FLTNB); 03731 //ap_IF->slice_orientation = 0; // slice orientation : (=0, default) 03732 //ap_IF->pat_rotation = 0; // patient rotation : supine (=0, default) 03733 //ap_IF->pat_orientation = 0; // slice orientation : head-in (=0, default) 03734 ap_IF->rescale_slope=1.; // multiplicative calibration values. 03735 ap_IF->rescale_intercept=0.; // additive calibration values. 03736 03737 // Just required for reading, not for writing 03738 //ap_IF->cmtx_size[0] = ap_ID->GetNbVoxX(); 03739 //ap_IF->cmtx_size[1] = ap_ID->GetNbVoxY(); 03740 //ap_IF->cmtx_size[2] = ap_ID->GetNbVoxZ(); 03741 //ap_IF->cvox_size[0] = ap_ID->GetVoxSizeX(); 03742 //ap_IF->cvox_size[1] = ap_ID->GetVoxSizeY(); 03743 //ap_IF->cvox_size[2] = ap_ID->GetVoxSizeZ(); 03744 //ap_IF->is_mtx_size_different = false; 03745 ap_IF->data_type = IntfKeyGetOutputImgDataType(ap_ID); 03746 ap_IF->study_duration = ap_ID->GetFinalTimeStopInSec(0) - 03747 ap_ID->GetFrameTimeStartInSec(0,0); 03748 for(int fr=0 ; fr<ap_ID->GetNbTimeFrames() ; fr++) 03749 { 03750 ap_IF->image_duration.push_back(ap_ID->GetFrameDurationInSec(0, fr)); 03751 ap_IF->frame_group_pause.push_back((fr == 0) ? 0 03752 : ap_ID->GetFrameTimeStartInSec(0,fr) - ap_ID->GetFrameTimeStopInSec(0,fr-1)); 03753 } 03754 ap_IF->nb_time_windows = ap_ID->GetNbRespGates() * 03755 ap_ID->GetNbCardGates(); 03756 //ap_IF->process_status; 03757 03758 //ap_IF->detector_heads = 0; 03759 ap_IF->nb_energy_windows = ap_IF->nb_resp_gates * 03760 ap_IF->nb_card_gates; 03761 //ap_IF->nb_projections 03762 //ap_IF->extent_rotation; 03763 //ap_IF->extent_rotation; 03764 //ap_IF->first_angle; 03765 //ap_IF->projection_angles; 03766 //ap_IF->radius; 03767 } 03768 03769 03770 03771 03772 // ===================================================================== 03773 // --------------------------------------------------------------------- 03774 // --------------------------------------------------------------------- 03775 // ===================================================================== 03776 /* 03777 \fn IntfKeyPrintFields 03778 \param ap_IF 03779 \brief Print all the keys of the Intf_fields structure passed in parameter, as well as their values for debugging purposes 03780 */ 03781 void IntfKeyPrintFields(Intf_fields a_IF) 03782 { 03783 Cout("// ------ IntfKeyPrintFields ------ // " << endl << endl); 03784 Cout("path_to_image : " << a_IF.path_to_image << endl); 03785 Cout("endianness : " << IntfKeyGetEndianStr(a_IF.endianness) << endl); 03786 Cout("data_offset : " << a_IF.data_offset << endl); 03787 Cout("nb_format : " << a_IF.nb_format << endl); 03788 Cout("nb_dims : " << unsigned(a_IF.nb_dims) << endl); // uint_8t doesn't output with cout, hence the (us) cast 03789 for(int i=0 ; i<7 ; i++) 03790 Cout("mtx_size["<<i<<"] : " << a_IF.mtx_size[i] << endl); 03791 for(int i=0 ; i<3 ; i++) 03792 Cout("vox_size["<<i<<"] : " << a_IF.vox_size[i] << endl); 03793 03794 Cout("slice_thickness_mm : " << a_IF.slice_thickness_mm << endl); 03795 Cout("ctr_to_ctr_separation : " << a_IF.ctr_to_ctr_separation << endl); 03796 Cout("nb_time_frames : " << a_IF.nb_time_frames << endl); 03797 Cout("nb_resp_gates : " << a_IF.nb_resp_gates << endl); 03798 Cout("nb_card_gates : " << a_IF.nb_card_gates << endl); 03799 Cout("nb_total_imgs : " << a_IF.nb_total_imgs << endl); 03800 Cout("nb_bytes_pixel : " << unsigned(a_IF.nb_bytes_pixel) << endl); // uint_8t doesn't output with cout, hence the (us) cast 03801 Cout("slice_orientation : " << a_IF.slice_orientation << endl); 03802 Cout("pat_rotation : " << a_IF.pat_rotation << endl); 03803 Cout("pat_orientation : " << a_IF.pat_orientation << endl); 03804 Cout("rescale_slope : " << a_IF.rescale_slope << endl); 03805 Cout("rescale_intercept : " << a_IF.rescale_intercept << endl); 03806 for(int i=0 ; i<3 ; i++) 03807 Cout("cmtx_size["<<i<<"] : " << a_IF.cmtx_size[i] << endl); 03808 for(int i=0 ; i<3 ; i++) 03809 Cout("cvox_size["<<i<<"] : " << a_IF.cvox_size[i] << endl); 03810 Cout("is_mtx_size_different : " << a_IF.is_mtx_size_different << endl); 03811 Cout("data_type : " << a_IF.data_type << endl); 03812 Cout("study_duration : " << a_IF.study_duration << endl); 03813 Cout("image_duration(fr) : " << endl); 03814 for(uint32_t fr=0 ; fr<a_IF.image_duration.size() ; fr++) 03815 Cout("image_duration("<<fr+1<<") : " << a_IF.image_duration.at(fr) << endl); 03816 Cout("pause_duration(fr) : " << endl); 03817 for(uint32_t fr=0 ; fr<a_IF.frame_group_pause.size() ; fr++) 03818 Cout("pause_duration("<<fr+1<<") : " << a_IF.frame_group_pause.at(fr) << endl); 03819 03820 //ap_IF->image_duration = vector<float>; 03821 Cout("nb_time_windows : " << a_IF.nb_time_windows << endl); 03822 Cout("process_status : " << a_IF.process_status << endl); 03823 03824 // SPECT and projection related data 03825 Cout("nb_detector_heads : " << a_IF.nb_detector_heads << endl); 03826 Cout("nb_energy_windows : " << a_IF.nb_energy_windows << endl); 03827 Cout("nb_projections : " << a_IF.nb_projections << endl); 03828 Cout("extent_rotation : " << a_IF.extent_rotation << endl); 03829 Cout("direction_rotation : " << a_IF.direction_rotation << endl); 03830 Cout("first_angle : " << a_IF.first_angle << endl); 03831 Cout("projection_angles : " << a_IF.projection_angles << endl); 03832 Cout("radius : " << a_IF.radius << endl); 03833 Cout("// ------ ------------------ ------ // " << endl << endl); 03834 } 03835 03836 03837 03838 03839 // ===================================================================== 03840 // --------------------------------------------------------------------- 03841 // --------------------------------------------------------------------- 03842 // ===================================================================== 03843 /* 03844 \fn IntfRecoverKey 03845 \param ap_Key : Structure to recover the parsed key components (key, value,..) 03846 \param a_line : String to process 03847 \brief Process the line passed in parameter and write the key information in the ap_Key Intf_key member structure 03848 \details .korig : Get original line without comments 03849 .kcase : Get key without spaces and without comments 03850 .klcase: Same as kcase, in lower case 03851 .kvalue: Value of the key, without spaces 03852 \todo Check that IntfToLowerCase() doesn't cause issue with some characters or some ASCII file format (unicode, etc..) 03853 \return 0 if success, positive value otherwise. 03854 */ 03855 int IntfRecoverKey(Intf_key* ap_Key, const string& a_line) 03856 { 03857 string intf_sep = ":="; 03858 03859 // Remove any comment from the original key line 03860 int pos = a_line.find_first_of(';',0); 03861 ap_Key->korig = a_line.substr(0, pos); 03862 03863 // check for interfile separator. 03864 pos = ap_Key->korig.find_first_of(intf_sep); 03865 03866 // If interfile key is not found (not an interfile key or wrong), just add it at the end of the key and proceed 03867 if (ap_Key->korig.find(intf_sep) == string::npos) 03868 ap_Key->korig.append(":="); 03869 03870 // Get key title 03871 ap_Key->kcase = ap_Key->korig.substr(0,pos); 03872 // Get key value 03873 ap_Key->kvalue = ap_Key->korig.substr(pos+2); 03874 03875 // Get case insensitive key, and without space 03876 ap_Key->klcase = ap_Key->kcase; 03877 IntfToLowerCase(&ap_Key->klcase); // TODO Safe ? 03878 03879 // Get key value 03880 ap_Key->kvalue = ap_Key->korig.substr(pos+2); 03881 03882 // Clear the space before the first and after the last element in the key and value; 03883 IntfEraseSpaces(&ap_Key->klcase); 03884 IntfEraseSpaces(&ap_Key->kvalue); 03885 03886 return 0; 03887 } 03888 03889 03890 03891 03892 // ===================================================================== 03893 // --------------------------------------------------------------------- 03894 // --------------------------------------------------------------------- 03895 // ===================================================================== 03896 /* 03897 \fn IntfCheckKeyMatch(Intf_key ap_Key, const string& a_field) 03898 \param ap_Key : Structure containing the parsed key components (key, value,..) 03899 \param a_line : String containing an interfile key 03900 \brief Check if the key matches the string passed in parameter 03901 \todo Be sure it is appropriate to use Intf_key.klcase for each key comparison 03902 \return 1 if success, 0 otherwise (not found). 03903 */ 03904 int IntfCheckKeyMatch(Intf_key ap_Key, const string& a_field) 03905 { 03906 // TODO Check if appropriate to use .klcase in any situation for key comparison 03907 if(ap_Key.klcase == a_field) 03908 return 1; 03909 else 03910 return 0; 03911 } 03912 03913 03914 03915 03916 // ===================================================================== 03917 // --------------------------------------------------------------------- 03918 // --------------------------------------------------------------------- 03919 // ===================================================================== 03920 /* 03921 \fn IntfKeyIsArray() 03922 \param ap_Key 03923 \brief Check if the key passed in parameter is an array (contains brackets '{' and '}' ) 03924 \return 1 if success, 0 otherwise (not array). 03925 */ 03926 int IntfKeyIsArray(Intf_key ap_Key) 03927 { 03928 if(ap_Key.kvalue.find("{") != string::npos && 03929 ap_Key.kvalue.find("}") != string::npos) 03930 return 1; 03931 03932 return 0; 03933 } 03934 03935 03936 03937 03938 // ===================================================================== 03939 // --------------------------------------------------------------------- 03940 // --------------------------------------------------------------------- 03941 // ===================================================================== 03942 /* 03943 \fn IntfKeyGetArrayNbElts 03944 \param a_Key 03945 \brief Return the number of elts in an Interfile array Key 03946 \return the number of elements in the array key, or negative value if error 03947 */ 03948 int IntfKeyGetArrayNbElts(Intf_key a_Key) 03949 { 03950 int nb_elts = 0; 03951 string val_str = a_Key.kvalue; 03952 03953 size_t pos = val_str.find_first_of('{')+1; 03954 03955 while (pos < val_str.length()) 03956 { 03957 size_t pos_c = val_str.find_first_of(",", pos); 03958 03959 // no comma found, looking for closing bracket 03960 if(pos_c == string::npos) 03961 { 03962 pos_c = val_str.find_first_of("}", pos); 03963 if(! (IntfKeyIsArray(a_Key)) ) 03964 { 03965 Cerr("***** IntfKeyGetArrayNbElts-> Error : closing bracket not found in interfile array key : "<< a_Key.korig << " !" << endl); 03966 return -1; 03967 } 03968 else 03969 return nb_elts+1; // add last elt before the end bracket 03970 } 03971 pos = pos_c+1; 03972 nb_elts++; 03973 } 03974 03975 // return error if we end of key if reached without closing bracket 03976 Cerr("***** IntfKeyGetArrayNbElts-> Error : closing bracket not found in interfile array key : "<< a_Key.korig << " !" << endl); 03977 return -1; 03978 } 03979 03980 03981 03982 03983 // ===================================================================== 03984 // --------------------------------------------------------------------- 03985 // --------------------------------------------------------------------- 03986 // ===================================================================== 03987 /* 03988 \fn IntfKeyGetMaxArrayKey 03989 \param ap_Key 03990 \brief Return the maximum value in an array (key value contains brackets '{,,}' ) 03991 \return the max value in the array key. 03992 */ 03993 int IntfKeyGetMaxArrayKey(Intf_key ap_Key) 03994 { 03995 int value, max=0; 03996 string val_str = ap_Key.kvalue; 03997 if (val_str == "") return(max); 03998 03999 size_t pos = val_str.find_first_of('{')+1; 04000 04001 while (pos < val_str.length()) 04002 { 04003 size_t pos_c = val_str.find_first_of(",", pos); 04004 04005 // no comma found, looking for closing bracket 04006 if(pos_c == string::npos) pos_c = val_str.find_first_of("}", pos); 04007 04008 if(ConvertFromString(val_str.substr(pos,pos_c-pos), &value) ) 04009 { 04010 Cerr("***** IntfKeyGetMaxArrayKey()-> An error occured when trying to recover the following value from the array key : "<< val_str.substr(pos,pos_c-pos) << " !" << endl); 04011 return 1; 04012 } 04013 04014 if (value > max) max = value; 04015 04016 pos = pos_c+1; 04017 } 04018 04019 return max; 04020 } 04021 04022 04023 04024 04025 // ===================================================================== 04026 // --------------------------------------------------------------------- 04027 // --------------------------------------------------------------------- 04028 // ===================================================================== 04029 /* 04030 \fn int IntfKeyGetArrayElts() 04031 \param ap_Key 04032 \param T* ap_return : Templated parameter in which the elts will be returned 04033 \brief Get all the elements in an array key in a templated array passed in parameter. 04034 It assumes the return variable has been instanciated with a correct number of elements. 04035 \return 0 if success, positive value otherwise 04036 */ 04037 template<class T> 04038 int IntfKeyGetArrayElts(Intf_key a_Key, T* ap_return) 04039 { 04040 string val_str = a_Key.kvalue; 04041 04042 // Check if interfile key 04043 if(! (IntfKeyIsArray(a_Key)) ) 04044 { 04045 Cerr("***** IntfKeyGetArrayElts-> Error : Problem reading the following interfile array key : "<< a_Key.korig << " !" << endl); 04046 return 1; 04047 } 04048 04049 if (val_str == "") 04050 { 04051 Cerr("***** IntfKeyGetArrayElts-> Error : no elements in the array key : "<< a_Key.korig << " !" << endl); 04052 return 1; 04053 } 04054 04055 size_t pos = val_str.find_first_of('{')+1; 04056 04057 int elt = 0; 04058 04059 while (pos < val_str.length()) 04060 { 04061 size_t pos_c = val_str.find_first_of(",", pos); 04062 04063 // no comma found, looking for closing bracket 04064 if(pos_c == string::npos) pos_c = val_str.find_first_of("}", pos); 04065 04066 if(ConvertFromString(val_str.substr(pos,pos_c-pos), &ap_return[elt]) ) 04067 { 04068 Cerr("***** IntfKeyGetMaxArrayKey()-> An error occured when trying to recover the following value from the array key : "<< val_str.substr(pos,pos_c-pos) << " !" << endl); 04069 return 1; 04070 } 04071 04072 pos = pos_c+1; 04073 elt++; 04074 } 04075 04076 return 0; 04077 } 04078 04079 04080 04081 04082 // ===================================================================== 04083 // --------------------------------------------------------------------- 04084 // --------------------------------------------------------------------- 04085 // ===================================================================== 04086 /* 04087 \fn IntfGetVoxIdxSHTOrientation(Intf_fields a_IF, int a_voxId) 04088 \param ap_IF 04089 \param a_voxId : index of the voxel in a 1-D image vector 04090 \brief Compute a voxel index corresponding to the default orientation (Sup/Hin/Trans) \n 04091 from the orientation informations contained in the Intf_fields passed in parameter 04092 \todo NOT CURRENTLY IMPLEMENTED 04093 \return new voxel index 04094 */ 04095 int IntfGetVoxIdxSHTOrientation(Intf_fields a_IF, int a_voxId) 04096 { 04097 int dimX=a_IF.vox_size[0], dimY=a_IF.vox_size[1], dimZ=a_IF.vox_size[2]; 04098 int dimXY=dimX*dimY; 04099 04100 // Get x,y,z original image indexes 04101 int z = a_voxId/dimXY; 04102 int y = (a_voxId - z*dimXY) / dimX; 04103 int x = a_voxId - z*dimXY - y*dimX; 04104 04105 // new image indexes 04106 int X = x; 04107 int Y = y; 04108 int Z = z; 04109 04110 // Convert to default orientations (TRANSVERSE / SUPINE / HEADIN) 04111 if(a_IF.slice_orientation == INTF_TRANSVERSE) // X=x, Y=y, Z=z 04112 { 04113 if(a_IF.pat_orientation == INTF_HEADIN) 04114 { 04115 if(a_IF.pat_rotation == INTF_SUPINE) 04116 { 04117 // nothing to change 04118 return a_voxId; 04119 } 04120 else // a_IF.pat_rotation == INTF_PRONE 04121 { 04122 X = dimX-x; 04123 Y = dimY-y; 04124 } 04125 } 04126 else // a_IF.pat_orientation == INTF_FEETIN 04127 { 04128 if(a_IF.pat_rotation == INTF_SUPINE) 04129 { 04130 X = dimX-x; 04131 Z = dimZ-z; 04132 } 04133 else // a_IF.pat_rotation == INTF_PRONE 04134 { 04135 Y = dimY-y; 04136 Z = dimZ-z; 04137 } 04138 } 04139 } 04140 else if(a_IF.slice_orientation == INTF_CORONAL) // X=x, Y=z, Z=y 04141 { 04142 if(a_IF.pat_orientation == INTF_HEADIN) 04143 { 04144 if(a_IF.pat_rotation == INTF_SUPINE) 04145 { 04146 Y = z; 04147 Z = y; 04148 } 04149 else // a_IF.pat_rotation == INTF_PRONE 04150 { 04151 X = dimX-x; 04152 Y = dimZ-z; 04153 Z = y; 04154 } 04155 } 04156 else // a_IF.pat_orientation == INTF_FEETIN 04157 { 04158 if(a_IF.pat_rotation == INTF_SUPINE) 04159 { 04160 X = dimX-x; 04161 Y = z; 04162 Z = dimY-y; 04163 } 04164 else // a_IF.pat_rotation == INTF_PRONE 04165 { 04166 X = dimX-x; 04167 Y = dimZ-z; 04168 Z = dimY-y; 04169 } 04170 } 04171 } 04172 else // a_IF.slice_orientation == INTF_SAGITTAL // X=z, Y=x, Z=y 04173 { 04174 if(a_IF.pat_orientation == INTF_HEADIN) 04175 { 04176 if(a_IF.pat_rotation == INTF_SUPINE) 04177 { 04178 X = dimZ-z; 04179 Y = dimX-x; 04180 Z = y; 04181 } 04182 else // a_IF.pat_rotation == INTF_PRONE 04183 { 04184 X = z; 04185 Y = x; 04186 Z = y; 04187 } 04188 } 04189 else // a_IF.pat_orientation == INTF_FEETIN 04190 { 04191 if(a_IF.pat_rotation == INTF_SUPINE) 04192 { 04193 X = z; 04194 Y = dimX-x; 04195 Z = dimY-y; 04196 } 04197 else // a_IF.pat_rotation == INTF_PRONE 04198 { 04199 X = dimZ-z; 04200 Y = x; 04201 Z = dimY-y; 04202 } 04203 } 04204 } 04205 04206 a_voxId = X + Y*dimX + Z*dimX*dimY; 04207 return a_voxId; 04208 } 04209 04210 04211 04212 04213 // ===================================================================== 04214 // --------------------------------------------------------------------- 04215 // --------------------------------------------------------------------- 04216 // ===================================================================== 04217 /* 04218 \fn IntfKeyGetEndianStr() 04219 \param a_val : 04220 \brief return the endian string corresponding to the value passed in parameter (see module INTF_ENDIANNESS). 04221 \return "BIG_ENDIAN" if 0, \n 04222 "LITTLE_ENDIAN" if 1, \n 04223 "UNKNOWN" otherwise 04224 */ 04225 string IntfKeyGetEndianStr(int a_val) 04226 { 04227 if(a_val == 0) return "BIGENDIAN"; 04228 if(a_val == 1) return "LITTLEENDIAN"; 04229 return "UNKNOWN"; 04230 } 04231 04232 04233 /* 04234 \fn IntfKeyGetModalityStr 04235 \param a_modalityIdx 04236 \brief Convert the integer provided in parameter to the string related 04237 to the corresponding modality as defined by the scanner objects 04238 \todo Add other modalities as we implement them 04239 \return string corresponding to the modality 04240 */ 04241 string IntfKeyGetModalityStr(int a_modalityIdx) 04242 { 04243 if(a_modalityIdx == 0) 04244 return "PET"; 04245 else if(a_modalityIdx == 1) 04246 return "SPECT"; 04247 /* 04248 else if(a_modalityIdx == 2) 04249 return "CT"; 04250 else if(a_modalityIdx == 3) // Pinhole 04251 return "SPECT"; 04252 */ 04253 else 04254 return "UNKNOWN"; 04255 } 04256 04257 04258 04259 04260 // ===================================================================== 04261 // --------------------------------------------------------------------- 04262 // --------------------------------------------------------------------- 04263 // ===================================================================== 04264 /* 04265 \fn int IntfKeyGetInputImgDataType() 04266 \param a_str : original string 04267 \brief Get the image data type corresponding to the the input string 04268 \return int value corresponding to the image data type (see module INTF_IMG_TYPE). 04269 */ 04270 int IntfKeyGetInputImgDataType(const string& a_str) 04271 { 04272 if (a_str == "static") 04273 return INTF_IMG_STATIC; 04274 if (a_str == "dynamic") 04275 return INTF_IMG_DYNAMIC; 04276 if (a_str == "gated") 04277 return INTF_IMG_GATED; 04278 if (a_str == "tomographic") 04279 return INTF_IMG_SPECT; 04280 if (a_str == "gspect") 04281 return INTF_IMG_GSPECT; 04282 if (a_str == "pet") 04283 return INTF_IMG_PET; 04284 04285 return INTF_IMG_UNKNOWN; 04286 } 04287 04288 04289 04290 04291 // ===================================================================== 04292 // --------------------------------------------------------------------- 04293 // --------------------------------------------------------------------- 04294 // ===================================================================== 04295 /* 04296 \fn IntfKeyGetOutImgDataType 04297 \param ap_ID 04298 \brief Get the image data type corresponding to the image metadata passed in parameter 04299 \return int value corresponding to the image data type (see module INTF_IMG_TYPE). 04300 */ 04301 int IntfKeyGetOutputImgDataType(oImageDimensionsAndQuantification* ap_ID) 04302 { 04303 sScannerManager* p_scanMgr = sScannerManager::GetInstance(); 04304 04305 int data_type = 0; 04306 04307 // Update data type key 04308 if(ap_ID->GetNbTimeFrames() > 1) 04309 data_type = INTF_IMG_DYNAMIC; 04310 else if(p_scanMgr->GetScannerType() == 2 && 04311 (ap_ID->GetNbRespGates() > 1 || 04312 ap_ID->GetNbCardGates() > 1 ) ) 04313 data_type = INTF_IMG_GSPECT; 04314 else if(ap_ID->GetNbRespGates() > 1 || 04315 ap_ID->GetNbCardGates() > 1 ) 04316 data_type = INTF_IMG_GATED; 04317 else if(p_scanMgr->GetScannerType() == 2) 04318 data_type = INTF_IMG_SPECT; 04319 else 04320 data_type = INTF_IMG_STATIC; 04321 04322 return data_type; 04323 } 04324 04325 04326 04327 04328 // ===================================================================== 04329 // --------------------------------------------------------------------- 04330 // --------------------------------------------------------------------- 04331 // ===================================================================== 04332 /* 04333 \fn IntfKeyGetPixTypeStr() 04334 \brief Return the string corresponding to the nb of bytes in the type FLTNB 04335 \return string corresponding to the pixel data type 04336 */ 04337 string IntfKeyGetPixTypeStr() 04338 { 04339 // Return either "long long float" (long double type), "long float" (doubletype ) or "short float" (float type) 04340 if (sizeof(FLTNB) == 4) return "short float"; 04341 else if (sizeof(FLTNB) == 8) return "long float"; 04342 else if (sizeof(FLTNB) == 16) return "long long float"; 04343 else 04344 { 04345 Cerr("***** oInterfileIO::IntfKeyGetPixTypeStr() -> Size of current float type (" << sizeof(FLTNB) << ") does not correspond to a known type !" << endl); 04346 Exit(EXIT_FAILURE); 04347 } 04348 } 04349 04350 04351 04352 04353 // ===================================================================== 04354 // --------------------------------------------------------------------- 04355 // --------------------------------------------------------------------- 04356 // ===================================================================== 04357 /* 04358 \fn IntfKeyGetPatOrientation() 04359 \param ap_IF 04360 \brief Get the complete patient orientation from an Intf_fields structure according to 04361 the values of keys 'slice orientation', 'patient rotation', and 'patient orientation' 04362 \return int value corresponding to the patient orientation (see module PATIENT_ORIENTATION). 04363 */ 04364 int IntfKeyGetPatOrientation(Intf_fields ap_IF) 04365 { 04366 switch (ap_IF.pat_rotation) 04367 { 04368 case INTF_SUPINE: 04369 switch (ap_IF.pat_orientation) 04370 { 04371 case INTF_HEADIN: 04372 switch (ap_IF.slice_orientation) 04373 { 04374 case INTF_TRANSVERSE: 04375 return(INTF_SUPINE_HEADIN_TRANSAXIAL); break; 04376 case INTF_SAGITTAL : 04377 return(INTF_SUPINE_HEADIN_SAGITTAL); break; 04378 case INTF_CORONAL : 04379 return(INTF_SUPINE_HEADIN_CORONAL); break; 04380 } 04381 break; 04382 04383 case INTF_FEETIN: 04384 switch(ap_IF.slice_orientation) 04385 { 04386 case INTF_TRANSVERSE: 04387 return(INTF_SUPINE_FEETIN_TRANSAXIAL); break; 04388 case INTF_SAGITTAL : 04389 return(INTF_SUPINE_FEETIN_SAGITTAL); break; 04390 case INTF_CORONAL : 04391 return(INTF_SUPINE_FEETIN_CORONAL); break; 04392 } 04393 break; 04394 } 04395 break; 04396 04397 case INTF_PRONE : 04398 switch (ap_IF.pat_orientation) 04399 { 04400 case INTF_HEADIN: 04401 switch (ap_IF.slice_orientation) 04402 { 04403 case INTF_TRANSVERSE: 04404 return(INTF_PRONE_HEADIN_TRANSAXIAL); break; 04405 case INTF_SAGITTAL : 04406 return(INTF_PRONE_HEADIN_SAGITTAL); break; 04407 case INTF_CORONAL : 04408 return(INTF_PRONE_HEADIN_CORONAL); break; 04409 } 04410 break; 04411 case INTF_FEETIN: 04412 switch (ap_IF.slice_orientation) 04413 { 04414 case INTF_TRANSVERSE : 04415 return(INTF_PRONE_FEETIN_TRANSAXIAL); break; 04416 case INTF_SAGITTAL : 04417 return(INTF_PRONE_FEETIN_SAGITTAL); break; 04418 case INTF_CORONAL : 04419 return(INTF_PRONE_FEETIN_CORONAL); break; 04420 } 04421 break; 04422 } 04423 break; 04424 } 04425 04426 return(INTF_SUPINE_HEADIN_TRANSAXIAL); // default 04427 } 04428 04429 04430 04431 04432 // ===================================================================== 04433 // --------------------------------------------------------------------- 04434 // --------------------------------------------------------------------- 04435 // ===================================================================== 04436 /* 04437 \fn GetUserEndianness() 04438 \brief Check user/host computer endianness and write it to the global variable User_Endianness 04439 \details This function should be called once during the initialization step of the algorithm (currently, the singleton initialization) 04440 \todo Maybe better compute it in a preprocessor macro 04441 */ 04442 void GetUserEndianness() 04443 { 04444 int n = 1; 04445 User_Endianness = (*(char *)&n) == 1 ? 1 : 0 ; 04446 } 04447 04448 04449 04450 04451 // ===================================================================== 04452 // --------------------------------------------------------------------- 04453 // --------------------------------------------------------------------- 04454 // ===================================================================== 04455 /* 04456 \fn IntfEraseSpaces() 04457 \param string* input_str 04458 \brief Erase space, blank characters (\t\r\n), and '!' before and after the characters in the string passed in parameter 04459 */ 04460 void IntfEraseSpaces(string* input_str) 04461 { 04462 input_str->erase(0, input_str->find_first_not_of(" !\t\r\n")); // Erase all blank stuff before first character 04463 input_str->erase(input_str->find_last_not_of(" \t\r\n")+1 , input_str->length()); // Erase all blank stuff after last character 04464 } 04465 04466 04467 04468 04469 // ===================================================================== 04470 // --------------------------------------------------------------------- 04471 // --------------------------------------------------------------------- 04472 // ===================================================================== 04473 /* 04474 \fn IntfToLowerCase() 04475 \param string* ap_str : original string 04476 \brief Set all characters of the string passed in parameter to lower case 04477 \todo May have issue with non ASCII characters, file decoding, etc.. 04478 */ 04479 void IntfToLowerCase(string* ap_str) 04480 { 04481 std::transform(ap_str->begin(), ap_str->end(), ap_str->begin(), ::tolower); 04482 } 04483 04484 04485 04486 04487 // ===================================================================== 04488 // --------------------------------------------------------------------- 04489 // --------------------------------------------------------------------- 04490 // ===================================================================== 04491 /* 04492 \fn SwapBytes() 04493 \param T *ap_type : Variable of type T 04494 \brief Use std::reverse to swap the bits of a variable of any type 04495 \details Used for Little/Big Endian conversion 04496 */ 04497 template <class T> 04498 void SwapBytes(T *ap_type) 04499 { 04500 char *buffer = reinterpret_cast<char*>(ap_type); 04501 std::reverse(buffer, buffer + sizeof(T)); 04502 04503 //int i, j; 04504 //for (i=0,j=bytes-1;i < (bytes/2); i++, j--) 04505 //{ 04506 // ptr[i]^=ptr[j]; ptr[j]^=ptr[i]; ptr[i]^=ptr[j]; 04507 //} 04508 04509 //for (int i = 0; i < size; ++i) 04510 // buffer[i] = ((buffer[i]>> 8) | (buffer[i] << 8)); 04511 }