Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: DicomHandler.cc 84839 2014-10-21 13:44:55Z gcosmo $ 26 // 27 // 27 /// \file medical/DICOM/src/DicomHandler.cc 28 /// \file medical/DICOM/src/DicomHandler.cc 28 /// \brief Implementation of the DicomHandler 29 /// \brief Implementation of the DicomHandler class 29 // 30 // 30 // The code was written by : 31 // The code was written by : 31 // *Louis Archambault louis.archambault@p 32 // *Louis Archambault louis.archambault@phy.ulaval.ca, 32 // *Luc Beaulieu beaulieu@phy.ulaval.ca 33 // *Luc Beaulieu beaulieu@phy.ulaval.ca 33 // +Vincent Hubert-Tremblay at tigre.2@sy 34 // +Vincent Hubert-Tremblay at tigre.2@sympatico.ca 34 // 35 // 35 // 36 // 36 // *Centre Hospitalier Universitaire de Quebec 37 // *Centre Hospitalier Universitaire de Quebec (CHUQ), 37 // Hotel-Dieu de Quebec, departement de Radio- 38 // Hotel-Dieu de Quebec, departement de Radio-oncologie 38 // 11 cote du palais. Quebec, QC, Canada, G1R 39 // 11 cote du palais. Quebec, QC, Canada, G1R 2J6 39 // tel (418) 525-4444 #6720 40 // tel (418) 525-4444 #6720 40 // fax (418) 691 5268 41 // fax (418) 691 5268 41 // 42 // 42 // + University Laval, Quebec (QC) Canada 43 // + University Laval, Quebec (QC) Canada 43 //******************************************** 44 //******************************************************* 44 // 45 // 45 //******************************************** 46 //******************************************************* 46 // 47 // 47 /// DicomHandler.cc : 48 /// DicomHandler.cc : 48 /// - Handling of DICM images 49 /// - Handling of DICM images 49 /// - Reading headers and pixels 50 /// - Reading headers and pixels 50 /// - Transforming pixel to density and 51 /// - Transforming pixel to density and creating *.g4dcm 51 /// files 52 /// files 52 //******************************************** 53 //******************************************************* 53 54 54 #include "DicomHandler.hh" 55 #include "DicomHandler.hh" 55 << 56 #include "DicomPhantomZSliceHeader.hh" << 57 #include "DicomPhantomZSliceMerged.hh" << 58 << 59 #include "G4ios.hh" << 60 #include "globals.hh" 56 #include "globals.hh" >> 57 #include "G4ios.hh" >> 58 #include <fstream> 61 59 62 #include <cctype> 60 #include <cctype> 63 #include <cstring> 61 #include <cstring> 64 #include <fstream> << 62 >> 63 #include "DicomPhantomZSliceHeader.hh" >> 64 #include "DicomPhantomZSliceMerged.hh" 65 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 67 67 68 DicomHandler* DicomHandler::fInstance = 0; << 68 DicomHandler* DicomHandler::fgInstance = 0; 69 69 70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 71 71 72 DicomHandler* DicomHandler::Instance() 72 DicomHandler* DicomHandler::Instance() 73 { 73 { 74 if (fInstance == 0) { << 74 return fgInstance; 75 static DicomHandler dicomhandler; << 76 fInstance = &dicomhandler; << 77 } << 78 return fInstance; << 79 } 75 } 80 76 81 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 78 83 G4String DicomHandler::GetDicomDataPath() << 79 DicomHandler::DicomHandler() >> 80 : DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512), >> 81 fCompression(0), fNFiles(0), fRows(0), fColumns(0), >> 82 fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0), >> 83 fPixelSpacingX(0.), fPixelSpacingY(0.), >> 84 fSliceThickness(0.), fSliceLocation(0.), >> 85 fRescaleIntercept(0), fRescaleSlope(0), >> 86 fLittleEndian(true), fImplicitEndian(false), >> 87 fPixelRepresentation(0), nbrequali(0), >> 88 valuedensity(NULL),valueCT(NULL),readCalibration(false), >> 89 mergedSlices(NULL),driverFile("Data.dat"),ct2densityFile("CT2Density.dat") 84 { 90 { 85 // default is current directory << 91 mergedSlices = new DicomPhantomZSliceMerged; 86 G4String driverPath = "."; << 87 // check environment << 88 const char* path = std::getenv("DICOM_PATH") << 89 << 90 if (path) { << 91 // if is set in environment << 92 return G4String(path); << 93 } << 94 else { << 95 // if DICOM_USE_HEAD, look for data instal << 96 #ifdef DICOM_USE_HEAD << 97 G4cerr << "Warning! DICOM was compiled wit << 98 << "the DICOM_PATH was not set!" << << 99 G4String datadir = G4GetEnv<G4String>("G4E << 100 if (datadir.length() > 0) { << 101 auto _last = datadir.rfind("/"); << 102 if (_last != std::string::npos) datadir. << 103 driverPath = datadir + "/DICOM1.1/DICOM_ << 104 G4int rc = setenv("DICOM_PATH", driverPa << 105 G4cerr << "\t --> Using '" << driverPath << 106 G4ConsumeParameters(rc); << 107 } << 108 #endif << 109 } << 110 return driverPath; << 111 } 92 } 112 93 113 //....oooOO0OOooo........oooOO0OOooo........oo << 114 << 115 G4String DicomHandler::GetDicomDataFile() << 116 { << 117 #if defined(DICOM_USE_HEAD) && defined(G4_DCMT << 118 return GetDicomDataPath() + "/Data.dat.new"; << 119 #else << 120 return GetDicomDataPath() + "/Data.dat"; << 121 #endif << 122 } << 123 94 124 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 125 96 126 #ifdef DICOM_USE_HEAD << 97 DicomHandler::~DicomHandler() 127 DicomHandler::DicomHandler() << 128 : DATABUFFSIZE(8192), << 129 LINEBUFFSIZE(5020), << 130 FILENAMESIZE(512), << 131 fCompression(0), << 132 fNFiles(0), << 133 fRows(0), << 134 fColumns(0), << 135 fBitAllocated(0), << 136 fMaxPixelValue(0), << 137 fMinPixelValue(0), << 138 fPixelSpacingX(0.), << 139 fPixelSpacingY(0.), << 140 fSliceThickness(0.), << 141 fSliceLocation(0.), << 142 fRescaleIntercept(0), << 143 fRescaleSlope(0), << 144 fLittleEndian(true), << 145 fImplicitEndian(false), << 146 fPixelRepresentation(0), << 147 fNbrequali(0), << 148 fValueDensity(NULL), << 149 fValueCT(NULL), << 150 fReadCalibration(false), << 151 fMergedSlices(NULL), << 152 fCt2DensityFile("null.dat") << 153 { 98 { 154 fMergedSlices = new DicomPhantomZSliceMerged << 155 fDriverFile = GetDicomDataFile(); << 156 G4cout << "Reading the DICOM_HEAD project " << 157 } << 158 #else << 159 DicomHandler::DicomHandler() << 160 : DATABUFFSIZE(8192), << 161 LINEBUFFSIZE(5020), << 162 FILENAMESIZE(512), << 163 fCompression(0), << 164 fNFiles(0), << 165 fRows(0), << 166 fColumns(0), << 167 fBitAllocated(0), << 168 fMaxPixelValue(0), << 169 fMinPixelValue(0), << 170 fPixelSpacingX(0.), << 171 fPixelSpacingY(0.), << 172 fSliceThickness(0.), << 173 fSliceLocation(0.), << 174 fRescaleIntercept(0), << 175 fRescaleSlope(0), << 176 fLittleEndian(true), << 177 fImplicitEndian(false), << 178 fPixelRepresentation(0), << 179 fNbrequali(0), << 180 fValueDensity(NULL), << 181 fValueCT(NULL), << 182 fReadCalibration(false), << 183 fMergedSlices(NULL), << 184 fDriverFile("Data.dat"), << 185 fCt2DensityFile("CT2Density.dat") << 186 { << 187 fMergedSlices = new DicomPhantomZSliceMerged << 188 } << 189 #endif << 190 //....oooOO0OOooo........oooOO0OOooo........oo << 191 99 192 DicomHandler::~DicomHandler() {} << 100 } 193 101 194 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 195 103 196 G4int DicomHandler::ReadFile(FILE* dicom, char 104 G4int DicomHandler::ReadFile(FILE* dicom, char* filename2) 197 { 105 { 198 G4cout << " ReadFile " << filename2 << G4end << 106 G4cout << " ReadFile " << filename2 << G4endl; 199 << 107 G4int returnvalue = 0; size_t rflag = 0; 200 G4int returnvalue = 0; << 108 char * buffer = new char[LINEBUFFSIZE]; 201 size_t rflag = 0; << 109 202 char* buffer = new char[LINEBUFFSIZE]; << 110 fImplicitEndian = false; 203 << 111 fLittleEndian = true; 204 fImplicitEndian = false; << 112 205 fLittleEndian = true; << 113 rflag = std::fread( buffer, 1, 128, dicom ); // The first 128 bytes 206 << 114 //are not important 207 rflag = std::fread(buffer, 1, 128, dicom); << 115 // Reads the "DICOM" letters 208 << 116 rflag = std::fread( buffer, 1, 4, dicom ); 209 << 117 // if there is no preamble, the FILE pointer is rewinded. 210 rflag = std::fread(buffer, 1, 4, dicom); << 118 if(std::strncmp("DICM", buffer, 4) != 0) { 211 // if there is no preamble, the FILE pointer << 119 std::fseek(dicom, 0, SEEK_SET); 212 if (std::strncmp("DICM", buffer, 4) != 0) { << 120 fImplicitEndian = true; 213 std::fseek(dicom, 0, SEEK_SET); << 121 } 214 fImplicitEndian = true; << 122 215 } << 123 short readGroupId; // identify the kind of input data 216 << 124 short readElementId; // identify a particular type information 217 short readGroupId; // identify the kind of << 125 short elementLength2; // deal with element length in 2 bytes 218 short readElementId; // identify a particul << 126 //unsigned int elementLength4; 219 short elementLength2; // deal with element << 127 // deal with element length in 4 bytes 220 // unsigned int eleme << 128 unsigned long elementLength4; // deal with element length in 4 bytes 221 // deal with element length in 4 bytes << 129 222 unsigned long elementLength4; // deal with << 130 char * data = new char[DATABUFFSIZE]; 223 << 131 224 char* data = new char[DATABUFFSIZE]; << 132 225 << 133 // Read information up to the pixel data 226 // Read information up to the pixel data << 134 while(true) { 227 while (true) { << 135 228 // Reading groups and elements : << 136 //Reading groups and elements : 229 readGroupId = 0; << 137 readGroupId = 0; 230 readElementId = 0; << 138 readElementId = 0; 231 // group ID << 139 // group ID 232 rflag = std::fread(buffer, 2, 1, dicom); << 140 rflag = std::fread(buffer, 2, 1, dicom); 233 GetValue(buffer, readGroupId); << 141 GetValue(buffer, readGroupId); 234 // element ID << 142 // element ID 235 rflag = std::fread(buffer, 2, 1, dicom); << 143 rflag = std::fread(buffer, 2, 1, dicom); 236 GetValue(buffer, readElementId); << 144 GetValue(buffer, readElementId); 237 << 145 238 // Creating a tag to be identified afterwa << 146 // Creating a tag to be identified afterward 239 G4int tagDictionary = readGroupId * 0x1000 << 147 G4int tagDictionary = readGroupId*0x10000 + readElementId; 240 << 148 241 // beginning of the pixels << 149 // beginning of the pixels 242 if (tagDictionary == 0x7FE00010) { << 150 if(tagDictionary == 0x7FE00010) { 243 // Following 2 fread's are modifications << 151 // Folling 2 fread's are modifications to original DICOM example 244 // (Jonathan Madsen) << 152 //(Jonathan Madsen) 245 if (!fImplicitEndian) rflag = std::fread << 153 rflag = std::fread(buffer,2,1,dicom); // Reserved 2 bytes 246 // (not used for pixels) << 154 // (not used for pixels) 247 rflag = std::fread(buffer, 4, 1, dicom); << 155 rflag = std::fread(buffer,4,1,dicom); // Element Length 248 // (not used for pixels) << 156 // (not used for pixels) 249 break; // Exit to ReadImageData() << 157 break; // Exit to ReadImageData() 250 } << 158 } 251 << 159 252 // VR or element length << 160 // VR or element length 253 rflag = std::fread(buffer, 2, 1, dicom); << 161 rflag = std::fread(buffer,2,1,dicom); 254 GetValue(buffer, elementLength2); << 162 GetValue(buffer, elementLength2); 255 << 163 256 // If value representation (VR) is OB, OW, << 164 // If value representation (VR) is OB, OW, SQ, UN, added OF and UT 257 // the next length is 32 bits << 165 //the next length is 32 bits 258 if ((elementLength2 == 0x424f || // "OB" << 166 if((elementLength2 == 0x424f || // "OB" 259 elementLength2 == 0x574f || // "OW" << 167 elementLength2 == 0x574f || // "OW" 260 elementLength2 == 0x464f || // "OF" << 168 elementLength2 == 0x464f || // "OF" 261 elementLength2 == 0x5455 || // "UT" << 169 elementLength2 == 0x5455 || // "UT" 262 elementLength2 == 0x5153 || // "SQ" << 170 elementLength2 == 0x5153 || // "SQ" 263 elementLength2 == 0x4e55) << 171 elementLength2 == 0x4e55) && // "UN" 264 && // "UN" << 172 !fImplicitEndian ) { // explicit VR 265 !fImplicitEndian) << 173 266 { // explicit VR << 174 rflag = std::fread(buffer, 2, 1, dicom); // Skip 2 reserved bytes 267 << 175 268 rflag = std::fread(buffer, 2, 1, dicom); << 269 << 270 // element length 176 // element length 271 rflag = std::fread(buffer, 4, 1, dicom); 177 rflag = std::fread(buffer, 4, 1, dicom); 272 GetValue(buffer, elementLength4); 178 GetValue(buffer, elementLength4); 273 << 179 274 if (elementLength2 == 0x5153) { << 180 if(elementLength2 == 0x5153) 275 if (elementLength4 == 0xFFFFFFFF) { << 181 { 276 read_undefined_nested(dicom); << 182 if(elementLength4 == 0xFFFFFFFF) 277 elementLength4 = 0; << 183 { 278 } << 184 read_undefined_nested( dicom ); 279 else { << 185 elementLength4=0; 280 if (read_defined_nested(dicom, eleme << 186 } else{ 281 G4Exception("DicomHandler::ReadFil << 187 if(read_defined_nested( dicom, elementLength4 )==0){ 282 "Function read_defined << 188 G4Exception("DicomHandler::ReadFile", >> 189 "DICOM001", >> 190 FatalException, >> 191 "Function read_defined_nested() failed!"); >> 192 } 283 } 193 } 284 } << 194 } else { 285 } << 286 else { << 287 // Reading the information with data 195 // Reading the information with data 288 rflag = std::fread(data, elementLength << 196 rflag = std::fread(data, elementLength4,1,dicom); 289 } 197 } 290 } << 198 291 else { << 199 >> 200 } else { >> 201 292 // explicit with VR different than prev 202 // explicit with VR different than previous ones 293 if (!fImplicitEndian || readGroupId == 2 << 203 if(!fImplicitEndian || readGroupId == 2) { 294 // G4cout << "Reading DICOM files wit << 204 295 // element length (2 bytes) << 205 //G4cout << "Reading DICOM files with Explicit VR"<< G4endl; >> 206 // element length (2 bytes) 296 rflag = std::fread(buffer, 2, 1, dicom 207 rflag = std::fread(buffer, 2, 1, dicom); 297 GetValue(buffer, elementLength2); 208 GetValue(buffer, elementLength2); 298 elementLength4 = elementLength2; 209 elementLength4 = elementLength2; 299 << 210 300 rflag = std::fread(data, elementLength 211 rflag = std::fread(data, elementLength4, 1, dicom); 301 } << 212 302 else { // Implicit VR << 213 } else { // Implicit VR 303 << 214 304 // G4cout << "Reading DICOM files wit << 215 //G4cout << "Reading DICOM files with Implicit VR"<< G4endl; 305 << 216 306 // element length (4 bytes) 217 // element length (4 bytes) 307 if (std::fseek(dicom, -2, SEEK_CUR) != << 218 if(std::fseek(dicom, -2, SEEK_CUR) != 0) { 308 G4Exception("DicomHandler::ReadFile" << 219 G4Exception("DicomHandler::ReadFile", >> 220 "DICOM001", >> 221 FatalException, >> 222 "fseek failed"); 309 } 223 } 310 << 224 311 rflag = std::fread(buffer, 4, 1, dicom 225 rflag = std::fread(buffer, 4, 1, dicom); 312 GetValue(buffer, elementLength4); 226 GetValue(buffer, elementLength4); 313 << 227 314 // G4cout << std::hex<< elementLength << 228 //G4cout << std::hex<< elementLength4 << G4endl; 315 << 229 316 if (elementLength4 == 0xFFFFFFFF) { << 230 if(elementLength4 == 0xFFFFFFFF) 317 read_undefined_nested(dicom); << 231 { 318 elementLength4 = 0; << 232 read_undefined_nested(dicom); 319 } << 233 elementLength4=0; 320 else { << 234 } else{ 321 rflag = std::fread(data, elementLeng 235 rflag = std::fread(data, elementLength4, 1, dicom); 322 } 236 } >> 237 >> 238 } 323 } 239 } >> 240 >> 241 // NULL termination >> 242 data[elementLength4] = '\0'; >> 243 >> 244 // analyzing information >> 245 GetInformation(tagDictionary, data); 324 } 246 } >> 247 >> 248 G4String fnameG4DCM = G4String(filename2) + ".g4dcm"; 325 249 326 // NULL termination << 250 // Perform functions originally written straight to file 327 data[elementLength4] = '\0'; << 251 DicomPhantomZSliceHeader* zslice = new DicomPhantomZSliceHeader(fnameG4DCM); 328 << 329 // analyzing information << 330 GetInformation(tagDictionary, data); << 331 } << 332 252 333 G4String fnameG4DCM = G4String(filename2) + << 253 std::map<G4float,G4String>::const_iterator ite; 334 << 254 for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ++ite){ 335 // Perform functions originally written stra << 255 zslice->AddMaterial(ite->second); 336 DicomPhantomZSliceHeader* zslice = new Dicom << 256 } 337 for (auto ite = fMaterialIndices.cbegin(); i << 257 338 zslice->AddMaterial(ite->second); << 258 zslice->SetNoVoxelX(fColumns/fCompression); 339 } << 259 zslice->SetNoVoxelY(fRows/fCompression); 340 << 260 zslice->SetNoVoxelZ(1); 341 zslice->SetNoVoxelsX(fColumns / fCompression << 342 zslice->SetNoVoxelsY(fRows / fCompression); << 343 zslice->SetNoVoxelsZ(1); << 344 261 345 zslice->SetMinX(-fPixelSpacingX * fColumns / << 262 zslice->SetMinX(-fPixelSpacingX*fColumns/2.); 346 zslice->SetMaxX(fPixelSpacingX * fColumns / << 263 zslice->SetMaxX(fPixelSpacingX*fColumns/2.); 347 264 348 zslice->SetMinY(-fPixelSpacingY * fRows / 2. << 265 zslice->SetMinY(-fPixelSpacingY*fRows/2.); 349 zslice->SetMaxY(fPixelSpacingY * fRows / 2.) << 266 zslice->SetMaxY(fPixelSpacingY*fRows/2.); 350 267 351 zslice->SetMinZ(fSliceLocation - fSliceThick << 268 zslice->SetMinZ(fSliceLocation-fSliceThickness/2.); 352 zslice->SetMaxZ(fSliceLocation + fSliceThick << 269 zslice->SetMaxZ(fSliceLocation+fSliceThickness/2.); 353 270 354 //=== << 271 //=== 355 272 356 ReadData(dicom, filename2); << 273 ReadData( dicom, filename2 ); 357 274 358 // DEPRECIATED << 275 // DEPRECIATED 359 // StoreData( foutG4DCM ); << 276 //StoreData( foutG4DCM ); 360 // foutG4DCM.close(); << 277 //foutG4DCM.close(); 361 278 362 StoreData(zslice); << 279 StoreData( zslice ); 363 280 364 // Dumped 2 file after DicomPhantomZSliceMer << 281 // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency 365 // zslice->DumpToFile(); << 282 //zslice->DumpToFile(); 366 283 367 fMergedSlices->AddZSlice(zslice); << 284 mergedSlices->AddZSlice(zslice); 368 285 369 // << 286 // 370 delete[] buffer; << 287 delete [] buffer; 371 delete[] data; << 288 delete [] data; 372 289 373 if (rflag) return returnvalue; << 290 if (rflag) return returnvalue; 374 return returnvalue; << 291 return returnvalue; 375 } 292 } 376 293 377 //....oooOO0OOooo........oooOO0OOooo........oo 294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 378 295 379 void DicomHandler::GetInformation(G4int& tagDi << 296 void DicomHandler::GetInformation(G4int & tagDictionary, char * data) 380 { 297 { 381 if (tagDictionary == 0x00280010) { // Numbe << 298 if(tagDictionary == 0x00280010 ) { // Number of Rows 382 GetValue(data, fRows); << 299 GetValue(data, fRows); 383 std::printf("[0x00280010] Rows -> %i\n", f << 300 std::printf("[0x00280010] Rows -> %i\n",fRows); 384 } << 301 385 else if (tagDictionary == 0x00280011) { // << 302 } else if(tagDictionary == 0x00280011 ) { // Number of fColumns 386 GetValue(data, fColumns); << 303 GetValue(data, fColumns); 387 std::printf("[0x00280011] Columns -> %i\n" << 304 std::printf("[0x00280011] Columns -> %i\n",fColumns); 388 } << 305 389 else if (tagDictionary == 0x00280102) { // << 306 } else if(tagDictionary == 0x00280102 ) { // High bits ( not used ) 390 short highBits; << 307 short highBits; 391 GetValue(data, highBits); << 308 GetValue(data, highBits); 392 std::printf("[0x00280102] High bits -> %i\ << 309 std::printf("[0x00280102] High bits -> %i\n",highBits); 393 } << 310 394 else if (tagDictionary == 0x00280100) { // << 311 } else if(tagDictionary == 0x00280100 ) { // Bits allocated 395 GetValue(data, fBitAllocated); << 312 GetValue(data, fBitAllocated); 396 std::printf("[0x00280100] Bits allocated - << 313 std::printf("[0x00280100] Bits allocated -> %i\n", fBitAllocated); 397 } << 314 398 else if (tagDictionary == 0x00280101) { // << 315 } else if(tagDictionary == 0x00280101 ) { // Bits stored ( not used ) 399 short bitStored; << 316 short bitStored; 400 GetValue(data, bitStored); << 317 GetValue(data, bitStored); 401 std::printf("[0x00280101] Bits stored -> % << 318 std::printf("[0x00280101] Bits stored -> %i\n",bitStored); 402 } << 319 403 else if (tagDictionary == 0x00280106) { // << 320 } else if(tagDictionary == 0x00280106 ) { // Min. pixel value 404 GetValue(data, fMinPixelValue); << 321 GetValue(data, fMinPixelValue); 405 std::printf("[0x00280106] Min. pixel value << 322 std::printf("[0x00280106] Min. pixel value -> %i\n", fMinPixelValue); 406 } << 323 407 else if (tagDictionary == 0x00280107) { // << 324 } else if(tagDictionary == 0x00280107 ) { // Max. pixel value 408 GetValue(data, fMaxPixelValue); << 325 GetValue(data, fMaxPixelValue); 409 std::printf("[0x00280107] Max. pixel value << 326 std::printf("[0x00280107] Max. pixel value -> %i\n", fMaxPixelValue); 410 } << 327 411 else if (tagDictionary == 0x00281053) { // << 328 } else if(tagDictionary == 0x00281053) { // Rescale slope 412 fRescaleSlope = atoi(data); << 329 fRescaleSlope = atoi(data); 413 std::printf("[0x00281053] Rescale Slope -> << 330 std::printf("[0x00281053] Rescale Slope -> %d\n", fRescaleSlope); 414 } << 331 415 else if (tagDictionary == 0x00281052) { // << 332 } else if(tagDictionary == 0x00281052 ) { // Rescalse intercept 416 fRescaleIntercept = atoi(data); << 333 fRescaleIntercept = atoi(data); 417 std::printf("[0x00281052] Rescale Intercep << 334 std::printf("[0x00281052] Rescale Intercept -> %d\n", 418 } << 335 fRescaleIntercept ); 419 else if (tagDictionary == 0x00280103) { << 336 420 // Pixel representation ( functions not d << 337 } else if(tagDictionary == 0x00280103 ) { 421 fPixelRepresentation = atoi(data); // 0: << 338 // Pixel representation ( functions not design to read signed bits ) 422 std::printf("[0x00280103] Pixel Representa << 339 fPixelRepresentation = atoi(data); // 0: unsigned 1: signed 423 if (fPixelRepresentation == 1) { << 340 std::printf("[0x00280103] Pixel Representation -> %i\n", 424 std::printf("### PIXEL REPRESENTATION = << 341 fPixelRepresentation); 425 std::printf("DICOM READING SCAN FOR UNSI << 342 if(fPixelRepresentation == 1 ) { 426 std::printf("ERROR !!!!!! -> \n"); << 343 std::printf("### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, "); >> 344 std::printf("DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE "); >> 345 std::printf("ERROR !!!!!! -> \n"); >> 346 } >> 347 >> 348 } else if(tagDictionary == 0x00080006 ) { // Modality >> 349 std::printf("[0x00080006] Modality -> %s\n", data); >> 350 >> 351 } else if(tagDictionary == 0x00080070 ) { // Manufacturer >> 352 std::printf("[0x00080070] Manufacturer -> %s\n", data); >> 353 >> 354 } else if(tagDictionary == 0x00080080 ) { // Institution Name >> 355 std::printf("[0x00080080] Institution Name -> %s\n", data); >> 356 >> 357 } else if(tagDictionary == 0x00080081 ) { // Institution Address >> 358 std::printf("[0x00080081] Institution Address -> %s\n", data); >> 359 >> 360 } else if(tagDictionary == 0x00081040 ) { // Institution Department Name >> 361 std::printf("[0x00081040] Institution Department Name -> %s\n", data); >> 362 >> 363 } else if(tagDictionary == 0x00081090 ) { // Manufacturer's Model Name >> 364 std::printf("[0x00081090] Manufacturer's Model Name -> %s\n", data); >> 365 >> 366 } else if(tagDictionary == 0x00181000 ) { // Device Serial Number >> 367 std::printf("[0x00181000] Device Serial Number -> %s\n", data); >> 368 >> 369 } else if(tagDictionary == 0x00080008 ) { // Image type ( not used ) >> 370 std::printf("[0x00080008] Image Types -> %s\n", data); >> 371 >> 372 } else if(tagDictionary == 0x00283000 ) { //Modality LUT Sequence(not used) >> 373 std::printf("[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data); >> 374 >> 375 } else if(tagDictionary == 0x00283002 ) { // LUT Descriptor ( not used ) >> 376 std::printf("[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data); >> 377 >> 378 } else if(tagDictionary == 0x00283003 ) { // LUT Explanation ( not used ) >> 379 std::printf("[0x00283003] LUT Explanation LO 1 -> %s\n", data); >> 380 >> 381 } else if(tagDictionary == 0x00283004 ) { // Modality LUT ( not used ) >> 382 std::printf("[0x00283004] Modality LUT Type LO 1 -> %s\n", data); >> 383 >> 384 } else if(tagDictionary == 0x00283006 ) { // LUT Data ( not used ) >> 385 std::printf("[0x00283006] LUT Data US or SS -> %s\n", data); >> 386 >> 387 } else if(tagDictionary == 0x00283010 ) { // VOI LUT ( not used ) >> 388 std::printf("[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data); >> 389 >> 390 } else if(tagDictionary == 0x00280120 ) { // Pixel Padding Value (not used) >> 391 std::printf("[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data); >> 392 >> 393 } else if(tagDictionary == 0x00280030 ) { // Pixel Spacing >> 394 G4String datas(data); >> 395 int iss = datas.find('\\'); >> 396 fPixelSpacingX = atof( datas.substr(0,iss).c_str() ); >> 397 fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() ); >> 398 >> 399 } else if(tagDictionary == 0x00200037 ) { // Image Orientation ( not used ) >> 400 std::printf("[0x00200037] Image Orientation (Phantom) -> %s\n", data); >> 401 >> 402 } else if(tagDictionary == 0x00200032 ) { // Image Position ( not used ) >> 403 std::printf("[0x00200032] Image Position (Phantom,mm) -> %s\n", data); >> 404 >> 405 } else if(tagDictionary == 0x00180050 ) { // Slice Thickness >> 406 fSliceThickness = atof(data); >> 407 std::printf("[0x00180050] Slice Thickness (mm) -> %f\n", fSliceThickness); >> 408 >> 409 } else if(tagDictionary == 0x00201041 ) { // Slice Location >> 410 fSliceLocation = atof(data); >> 411 std::printf("[0x00201041] Slice Location -> %f\n", fSliceLocation); >> 412 >> 413 } else if(tagDictionary == 0x00280004 ) { // Photometric Interpretation >> 414 // ( not used ) >> 415 std::printf("[0x00280004] Photometric Interpretation -> %s\n", data); >> 416 >> 417 } else if(tagDictionary == 0x00020010) { // Endian >> 418 if(strcmp(data, "1.2.840.10008.1.2") == 0) >> 419 fImplicitEndian = true; >> 420 else if(strncmp(data, "1.2.840.10008.1.2.2", 19) == 0) >> 421 fLittleEndian = false; >> 422 //else 1.2.840..10008.1.2.1 (explicit little endian) >> 423 >> 424 std::printf("[0x00020010] Endian -> %s\n", data); 427 } 425 } 428 } << 429 else if (tagDictionary == 0x00080006) { // << 430 std::printf("[0x00080006] Modality -> %s\n << 431 } << 432 else if (tagDictionary == 0x00080070) { // << 433 std::printf("[0x00080070] Manufacturer -> << 434 } << 435 else if (tagDictionary == 0x00080080) { // << 436 std::printf("[0x00080080] Institution Name << 437 } << 438 else if (tagDictionary == 0x00080081) { // << 439 std::printf("[0x00080081] Institution Addr << 440 } << 441 else if (tagDictionary == 0x00081040) { // << 442 std::printf("[0x00081040] Institution Depa << 443 } << 444 else if (tagDictionary == 0x00081090) { // << 445 std::printf("[0x00081090] Manufacturer's M << 446 } << 447 else if (tagDictionary == 0x00181000) { // << 448 std::printf("[0x00181000] Device Serial Nu << 449 } << 450 else if (tagDictionary == 0x00080008) { // << 451 std::printf("[0x00080008] Image Types -> % << 452 } << 453 else if (tagDictionary == 0x00283000) { // << 454 std::printf("[0x00283000] Modality LUT Seq << 455 } << 456 else if (tagDictionary == 0x00283002) { // << 457 std::printf("[0x00283002] LUT Descriptor U << 458 } << 459 else if (tagDictionary == 0x00283003) { // << 460 std::printf("[0x00283003] LUT Explanation << 461 } << 462 else if (tagDictionary == 0x00283004) { // << 463 std::printf("[0x00283004] Modality LUT Typ << 464 } << 465 else if (tagDictionary == 0x00283006) { // << 466 std::printf("[0x00283006] LUT Data US or S << 467 } << 468 else if (tagDictionary == 0x00283010) { // << 469 std::printf("[0x00283010] VOI LUT Sequence << 470 } << 471 else if (tagDictionary == 0x00280120) { // << 472 std::printf("[0x00280120] Pixel Padding Va << 473 } << 474 else if (tagDictionary == 0x00280030) { // << 475 G4String datas(data); << 476 G4int iss = G4int(datas.find('\\')); << 477 fPixelSpacingX = atof(datas.substr(0, iss) << 478 fPixelSpacingY = atof(datas.substr(iss + 1 << 479 } << 480 else if (tagDictionary == 0x00200037) { // << 481 std::printf("[0x00200037] Image Orientatio << 482 } << 483 else if (tagDictionary == 0x00200032) { // << 484 std::printf("[0x00200032] Image Position ( << 485 } << 486 else if (tagDictionary == 0x00180050) { // << 487 fSliceThickness = atof(data); << 488 std::printf("[0x00180050] Slice Thickness << 489 } << 490 else if (tagDictionary == 0x00201041) { // << 491 fSliceLocation = atof(data); << 492 std::printf("[0x00201041] Slice Location - << 493 } << 494 else if (tagDictionary == 0x00280004) { // << 495 // ( not used ) << 496 std::printf("[0x00280004] Photometric Inte << 497 } << 498 else if (tagDictionary == 0x00020010) { // << 499 if (strcmp(data, "1.2.840.10008.1.2") == 0 << 500 fImplicitEndian = true; << 501 else if (strncmp(data, "1.2.840.10008.1.2. << 502 fLittleEndian = false; << 503 // else 1.2.840..10008.1.2.1 (explicit lit << 504 426 505 std::printf("[0x00020010] Endian -> %s\n", << 427 // others 506 } << 428 else { >> 429 //std::printf("[0x%x] -> %s\n", tagDictionary, data); >> 430 ; >> 431 } 507 432 508 // others << 509 else { << 510 // std::printf("[0x%x] -> %s\n", tagDictio << 511 ; << 512 } << 513 } 433 } 514 434 515 //....oooOO0OOooo........oooOO0OOooo........oo 435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 516 436 517 void DicomHandler::StoreData(DicomPhantomZSlic 437 void DicomHandler::StoreData(DicomPhantomZSliceHeader* dcmPZSH) 518 { 438 { 519 G4int mean; << 439 G4int mean; 520 G4double density; << 440 G4double density; 521 G4bool overflow = false; << 441 G4bool overflow = false; 522 << 442 523 if (!dcmPZSH) { << 443 if(!dcmPZSH) { return; } 524 return; << 444 525 } << 445 dcmPZSH->SetSliceLocation(fSliceLocation); 526 << 446 527 dcmPZSH->SetSliceLocation(fSliceLocation); << 447 //----- Print indices of material 528 << 448 if(fCompression == 1) { // no fCompression: each pixel has a density value) 529 //----- Print indices of material << 449 for( G4int ww = 0; ww < fRows; ww++) { 530 if (fCompression == 1) { // no fCompression << 450 dcmPZSH->AddRow(); 531 for (G4int ww = 0; ww < fRows; ++ww) { << 451 for( G4int xx = 0; xx < fColumns; xx++) { 532 dcmPZSH->AddRow(); << 452 mean = fTab[ww][xx]; 533 for (G4int xx = 0; xx < fColumns; ++xx) << 453 density = Pixel2density(mean); 534 mean = fTab[ww][xx]; << 454 dcmPZSH->AddValue(density); 535 density = Pixel2density(mean); << 455 dcmPZSH->AddMateID(GetMaterialIndex(density)); 536 dcmPZSH->AddValue(density); << 456 } 537 dcmPZSH->AddMateID(GetMaterialIndex(G4 << 538 } << 539 } << 540 } << 541 else { << 542 // density value is the average of a squar << 543 // fCompression*fCompression pixels << 544 for (G4int ww = 0; ww < fRows; ww += fComp << 545 dcmPZSH->AddRow(); << 546 for (G4int xx = 0; xx < fColumns; xx += << 547 overflow = false; << 548 mean = 0; << 549 for (G4int sumx = 0; sumx < fCompressi << 550 for (G4int sumy = 0; sumy < fCompres << 551 if (ww + sumy >= fRows || xx + sum << 552 mean += fTab[ww + sumy][xx + sumx] << 553 } << 554 if (overflow) break; << 555 } 457 } 556 mean /= fCompression * fCompression; << 557 458 558 if (!overflow) { << 459 } else { 559 density = Pixel2density(mean); << 460 // density value is the average of a square region of 560 dcmPZSH->AddValue(density); << 461 // fCompression*fCompression pixels 561 dcmPZSH->AddMateID(GetMaterialIndex( << 462 for(G4int ww = 0; ww < fRows ;ww += fCompression ) { >> 463 dcmPZSH->AddRow(); >> 464 for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) { >> 465 overflow = false; >> 466 mean = 0; >> 467 for(int sumx = 0; sumx < fCompression; sumx++) { >> 468 for(int sumy = 0; sumy < fCompression; sumy++) { >> 469 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true; >> 470 mean += fTab[ww+sumy][xx+sumx]; >> 471 } >> 472 if(overflow) break; >> 473 } >> 474 mean /= fCompression*fCompression; >> 475 >> 476 if(!overflow) { >> 477 density = Pixel2density(mean); >> 478 dcmPZSH->AddValue(density); >> 479 dcmPZSH->AddMateID(GetMaterialIndex(density)); >> 480 } 562 } 481 } 563 } 482 } 564 } 483 } 565 } << 484 566 << 485 dcmPZSH->FlipData(); 567 dcmPZSH->FlipData(); << 568 } 486 } 569 487 570 //....oooOO0OOooo........oooOO0OOooo........oo 488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 571 // This function is depreciated as it is handl << 489 // This function is depreciated as it is handled by 572 // DicomPhantomZSliceHeader::DumpToFile 490 // DicomPhantomZSliceHeader::DumpToFile 573 void DicomHandler::StoreData(std::ofstream& fo 491 void DicomHandler::StoreData(std::ofstream& foutG4DCM) 574 { 492 { 575 G4int mean; 493 G4int mean; 576 G4double density; 494 G4double density; 577 G4bool overflow = false; 495 G4bool overflow = false; 578 << 496 579 //----- Print indices of material 497 //----- Print indices of material 580 if (fCompression == 1) { // no fCompression << 498 if(fCompression == 1) { // no fCompression: each pixel has a density value) 581 for (G4int ww = 0; ww < fRows; ++ww) { << 499 for( G4int ww = 0; ww < fRows; ww++) { 582 for (G4int xx = 0; xx < fColumns; ++xx) << 500 for( G4int xx = 0; xx < fColumns; xx++) { 583 mean = fTab[ww][xx]; 501 mean = fTab[ww][xx]; 584 density = Pixel2density(mean); 502 density = Pixel2density(mean); 585 foutG4DCM << GetMaterialIndex(G4float( << 503 foutG4DCM << GetMaterialIndex( density ) << " "; 586 } 504 } 587 foutG4DCM << G4endl; 505 foutG4DCM << G4endl; 588 } 506 } 589 } << 507 590 else { << 508 } else { 591 // density value is the average of a squar 509 // density value is the average of a square region of 592 // fCompression*fCompression pixels 510 // fCompression*fCompression pixels 593 for (G4int ww = 0; ww < fRows; ww += fComp << 511 for(G4int ww = 0; ww < fRows ;ww += fCompression ) { 594 for (G4int xx = 0; xx < fColumns; xx += << 512 for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) { 595 overflow = false; 513 overflow = false; 596 mean = 0; 514 mean = 0; 597 for (G4int sumx = 0; sumx < fCompressi << 515 for(int sumx = 0; sumx < fCompression; sumx++) { 598 for (G4int sumy = 0; sumy < fCompres << 516 for(int sumy = 0; sumy < fCompression; sumy++) { 599 if (ww + sumy >= fRows || xx + sum << 517 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true; 600 mean += fTab[ww + sumy][xx + sumx] << 518 mean += fTab[ww+sumy][xx+sumx]; 601 } 519 } 602 if (overflow) break; << 520 if(overflow) break; 603 } 521 } 604 mean /= fCompression * fCompression; << 522 mean /= fCompression*fCompression; 605 << 523 606 if (!overflow) { << 524 if(!overflow) { 607 density = Pixel2density(mean); 525 density = Pixel2density(mean); 608 foutG4DCM << GetMaterialIndex(G4floa << 526 foutG4DCM << GetMaterialIndex( density ) << " "; 609 } 527 } 610 } 528 } 611 foutG4DCM << G4endl; 529 foutG4DCM << G4endl; 612 } 530 } >> 531 613 } 532 } 614 << 533 615 //----- Print densities 534 //----- Print densities 616 if (fCompression == 1) { // no fCompression << 535 if(fCompression == 1) { // no fCompression: each pixel has a density value) 617 for (G4int ww = 0; ww < fRows; ww++) { << 536 for( G4int ww = 0; ww < fRows; ww++) { 618 for (G4int xx = 0; xx < fColumns; xx++) << 537 for( G4int xx = 0; xx < fColumns; xx++) { 619 mean = fTab[ww][xx]; 538 mean = fTab[ww][xx]; 620 density = Pixel2density(mean); 539 density = Pixel2density(mean); 621 foutG4DCM << density << " "; 540 foutG4DCM << density << " "; 622 if (xx % 8 == 3) foutG4DCM << G4endl; << 541 if( xx%8 == 3 ) foutG4DCM << G4endl; // just for nicer reading 623 } 542 } 624 } 543 } 625 } << 544 626 else { << 545 } else { 627 // density value is the average of a squar 546 // density value is the average of a square region of 628 // fCompression*fCompression pixels 547 // fCompression*fCompression pixels 629 for (G4int ww = 0; ww < fRows; ww += fComp << 548 for(G4int ww = 0; ww < fRows ;ww += fCompression ) { 630 for (G4int xx = 0; xx < fColumns; xx += << 549 for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) { 631 overflow = false; 550 overflow = false; 632 mean = 0; 551 mean = 0; 633 for (G4int sumx = 0; sumx < fCompressi << 552 for(int sumx = 0; sumx < fCompression; sumx++) { 634 for (G4int sumy = 0; sumy < fCompres << 553 for(int sumy = 0; sumy < fCompression; sumy++) { 635 if (ww + sumy >= fRows || xx + sum << 554 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true; 636 mean += fTab[ww + sumy][xx + sumx] << 555 mean += fTab[ww+sumy][xx+sumx]; 637 } 556 } 638 if (overflow) break; << 557 if(overflow) break; 639 } 558 } 640 mean /= fCompression * fCompression; << 559 mean /= fCompression*fCompression; 641 << 560 642 if (!overflow) { << 561 if(!overflow) { 643 density = Pixel2density(mean); 562 density = Pixel2density(mean); 644 foutG4DCM << density << " "; << 563 foutG4DCM << density << " "; 645 if (xx / fCompression % 8 == 3) fout << 564 if( xx/fCompression%8 == 3 ) foutG4DCM << G4endl; // just for nicer 646 // reading 565 // reading 647 } 566 } 648 } 567 } 649 } 568 } >> 569 650 } 570 } >> 571 651 } 572 } 652 573 653 //....oooOO0OOooo........oooOO0OOooo........oo << 574 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 654 void DicomHandler::ReadMaterialIndices(std::if << 575 void DicomHandler::ReadMaterialIndices( std::ifstream& finData) 655 { 576 { 656 unsigned int nMate; 577 unsigned int nMate; 657 G4String mateName; 578 G4String mateName; 658 G4float densityMax; 579 G4float densityMax; 659 finData >> nMate; 580 finData >> nMate; 660 if (finData.eof()) return; << 581 if( finData.eof() ) return; 661 << 582 662 G4cout << " ReadMaterialIndices " << nMate < 583 G4cout << " ReadMaterialIndices " << nMate << G4endl; 663 for (unsigned int ii = 0; ii < nMate; ++ii) << 584 for( unsigned int ii = 0; ii < nMate; ii++ ){ 664 finData >> mateName >> densityMax; 585 finData >> mateName >> densityMax; 665 fMaterialIndices[densityMax] = mateName; 586 fMaterialIndices[densityMax] = mateName; 666 // G4cout << ii << " ReadMaterialIndice << 587 G4cout << ii << " ReadMaterialIndices " << mateName << " " 667 // << densityMax << G4endl; << 588 << densityMax << G4endl; 668 } 589 } >> 590 669 } 591 } 670 592 671 //....oooOO0OOooo........oooOO0OOooo........oo 593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 672 594 673 unsigned int DicomHandler::GetMaterialIndex(G4 << 595 unsigned int DicomHandler::GetMaterialIndex( G4float density ) 674 { 596 { 675 std::map<G4float, G4String>::const_reverse_i << 597 std::map<G4float,G4String>::reverse_iterator ite; 676 std::size_t ii = fMaterialIndices.size(); << 598 G4int ii = fMaterialIndices.size(); 677 << 599 for( ite = fMaterialIndices.rbegin(); ite != fMaterialIndices.rend(); 678 for (ite = fMaterialIndices.crbegin(); ite ! << 600 ite++, ii-- ) { 679 if (density >= (*ite).first) { << 601 if( density >= (*ite).first ) { 680 break; 602 break; 681 } 603 } 682 } 604 } 683 << 605 //- G4cout << " GetMaterialIndex " << density << " = " << ii << G4endl; 684 if (ii == fMaterialIndices.size()) { << 606 685 ii = fMaterialIndices.size() - 1; << 607 if(static_cast<unsigned int>(ii) == fMaterialIndices.size()) 686 } << 608 { ii = fMaterialIndices.size()-1; } 687 << 609 688 return unsigned(ii); << 610 return ii; >> 611 689 } 612 } 690 613 691 //....oooOO0OOooo........oooOO0OOooo........oo 614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 692 615 693 G4int DicomHandler::ReadData(FILE* dicom, char << 616 G4int DicomHandler::ReadData(FILE *dicom,char * filename2) 694 { 617 { 695 G4int returnvalue = 0; << 618 G4int returnvalue = 0; size_t rflag = 0; 696 size_t rflag = 0; << 619 697 << 620 // READING THE PIXELS : 698 // READING THE PIXELS << 621 G4int w = 0; 699 << 622 700 fTab = new G4int*[fRows]; 623 fTab = new G4int*[fRows]; 701 for (G4int i = 0; i < fRows; ++i) { << 624 for ( G4int i = 0; i < fRows; i ++ ) { 702 fTab[i] = new G4int[fColumns]; << 625 fTab[i] = new G4int[fColumns]; 703 } << 626 } 704 << 627 705 if (fBitAllocated == 8) { // Case 8 bits : << 628 if(fBitAllocated == 8) { // Case 8 bits : 706 << 629 707 std::printf("@@@ Error! Picture != 16 bits 630 std::printf("@@@ Error! Picture != 16 bits...\n"); 708 std::printf("@@@ Error! Picture != 16 bits 631 std::printf("@@@ Error! Picture != 16 bits...\n"); 709 std::printf("@@@ Error! Picture != 16 bits 632 std::printf("@@@ Error! Picture != 16 bits...\n"); 710 << 633 711 unsigned char ch = 0; 634 unsigned char ch = 0; 712 << 635 713 for (G4int j = 0; j < fRows; ++j) { << 636 for(G4int j = 0; j < fRows; j++) { 714 for (G4int i = 0; i < fColumns; ++i) { << 637 for(G4int i = 0; i < fColumns; i++) { 715 rflag = std::fread(&ch, 1, 1, dicom); << 638 w++; 716 fTab[j][i] = ch * fRescaleSlope + fRes << 639 rflag = std::fread( &ch, 1, 1, dicom); >> 640 fTab[j][i] = ch*fRescaleSlope + fRescaleIntercept; 717 } 641 } 718 } 642 } 719 returnvalue = 1; 643 returnvalue = 1; 720 } << 644 721 else { // from 12 to 16 bits : << 645 } else { // from 12 to 16 bits : 722 char sbuff[2]; 646 char sbuff[2]; 723 short pixel; 647 short pixel; 724 for (G4int j = 0; j < fRows; ++j) { << 648 for( G4int j = 0; j < fRows; j++) { 725 for (G4int i = 0; i < fColumns; ++i) { << 649 for( G4int i = 0; i < fColumns; i++) { >> 650 w++; 726 rflag = std::fread(sbuff, 2, 1, dicom) 651 rflag = std::fread(sbuff, 2, 1, dicom); 727 GetValue(sbuff, pixel); 652 GetValue(sbuff, pixel); 728 fTab[j][i] = pixel * fRescaleSlope + f << 653 fTab[j][i] = pixel*fRescaleSlope + fRescaleIntercept; 729 } 654 } 730 } 655 } 731 } 656 } 732 << 657 733 // Creation of .g4 files wich contains avera 658 // Creation of .g4 files wich contains averaged density data 734 G4String nameProcessed = filename2 + G4Strin << 659 char * nameProcessed = new char[FILENAMESIZE]; 735 FILE* fileOut = std::fopen(nameProcessed.c_s << 660 FILE* fileOut; 736 << 661 737 G4cout << "### Writing of " << nameProcessed << 662 std::sprintf(nameProcessed,"%s.g4dcmb",filename2); 738 << 663 fileOut = std::fopen(nameProcessed,"w+b"); >> 664 std::printf("### Writing of %s ###\n",nameProcessed); >> 665 739 unsigned int nMate = fMaterialIndices.size() 666 unsigned int nMate = fMaterialIndices.size(); 740 rflag = std::fwrite(&nMate, sizeof(unsigned 667 rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut); 741 //--- Write materials 668 //--- Write materials 742 for (auto ite = fMaterialIndices.cbegin(); i << 669 std::map<G4float,G4String>::const_iterator ite; >> 670 for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++ ){ 743 G4String mateName = (*ite).second; 671 G4String mateName = (*ite).second; 744 for (std::size_t ii = (*ite).second.length << 672 for( G4int ii = (*ite).second.length(); ii < 40; ii++ ) { 745 mateName += " "; 673 mateName += " "; 746 } // mateName = const_cast<char*>(((*ite) << 674 } //mateName = const_cast<char*>(((*ite).second).c_str()); 747 << 675 748 const char* mateNameC = mateName.c_str(); 676 const char* mateNameC = mateName.c_str(); 749 rflag = std::fwrite(mateNameC, sizeof(char << 677 rflag = std::fwrite(mateNameC, sizeof(char),40, fileOut); 750 } 678 } 751 << 679 752 unsigned int fRowsC = fRows / fCompression; << 680 unsigned int fRowsC = fRows/fCompression; 753 unsigned int fColumnsC = fColumns / fCompres << 681 unsigned int fColumnsC = fColumns/fCompression; 754 unsigned int planesC = 1; 682 unsigned int planesC = 1; 755 G4float pixelLocationXM = -G4float(fPixelSpa << 683 G4float pixelLocationXM = -fPixelSpacingX*fColumns/2.; 756 G4float pixelLocationXP = G4float(fPixelSpac << 684 G4float pixelLocationXP = fPixelSpacingX*fColumns/2.; 757 G4float pixelLocationYM = -G4float(fPixelSpa << 685 G4float pixelLocationYM = -fPixelSpacingY*fRows/2.; 758 G4float pixelLocationYP = G4float(fPixelSpac << 686 G4float pixelLocationYP = fPixelSpacingY*fRows/2.; 759 G4float fSliceLocationZM = G4float(fSliceLoc << 687 G4float fSliceLocationZM = fSliceLocation-fSliceThickness/2.; 760 G4float fSliceLocationZP = G4float(fSliceLoc << 688 G4float fSliceLocationZP = fSliceLocation+fSliceThickness/2.; 761 //--- Write number of voxels (assume only on 689 //--- Write number of voxels (assume only one voxel in Z) 762 rflag = std::fwrite(&fRowsC, sizeof(unsigned 690 rflag = std::fwrite(&fRowsC, sizeof(unsigned int), 1, fileOut); 763 rflag = std::fwrite(&fColumnsC, sizeof(unsig 691 rflag = std::fwrite(&fColumnsC, sizeof(unsigned int), 1, fileOut); 764 rflag = std::fwrite(&planesC, sizeof(unsigne 692 rflag = std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut); 765 //--- Write minimum and maximum extensions 693 //--- Write minimum and maximum extensions 766 rflag = std::fwrite(&pixelLocationXM, sizeof 694 rflag = std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut); 767 rflag = std::fwrite(&pixelLocationXP, sizeof 695 rflag = std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut); 768 rflag = std::fwrite(&pixelLocationYM, sizeof 696 rflag = std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut); 769 rflag = std::fwrite(&pixelLocationYP, sizeof 697 rflag = std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut); 770 rflag = std::fwrite(&fSliceLocationZM, sizeo 698 rflag = std::fwrite(&fSliceLocationZM, sizeof(G4float), 1, fileOut); 771 rflag = std::fwrite(&fSliceLocationZP, sizeo 699 rflag = std::fwrite(&fSliceLocationZP, sizeof(G4float), 1, fileOut); 772 // rflag = std::fwrite(&fCompression, sizeof 700 // rflag = std::fwrite(&fCompression, sizeof(unsigned int), 1, fileOut); 773 << 701 774 std::printf("%8i %8i\n", fRows, fColumns); << 702 std::printf("%8i %8i\n",fRows,fColumns); 775 std::printf("%8f %8f\n", fPixelSpacingX, f << 703 std::printf("%8f %8f\n",fPixelSpacingX,fPixelSpacingY); 776 std::printf("%8f\n", fSliceThickness); 704 std::printf("%8f\n", fSliceThickness); 777 std::printf("%8f\n", fSliceLocation); 705 std::printf("%8f\n", fSliceLocation); 778 std::printf("%8i\n", fCompression); 706 std::printf("%8i\n", fCompression); 779 << 707 780 G4int compSize = fCompression; 708 G4int compSize = fCompression; 781 G4int mean; 709 G4int mean; 782 G4float density; 710 G4float density; 783 G4bool overflow = false; 711 G4bool overflow = false; 784 << 712 785 //----- Write index of material for each pix 713 //----- Write index of material for each pixel 786 if (compSize == 1) { // no fCompression: ea << 714 if(compSize == 1) { // no fCompression: each pixel has a density value) 787 for (G4int ww = 0; ww < fRows; ++ww) { << 715 for( G4int ww = 0; ww < fRows; ww++) { 788 for (G4int xx = 0; xx < fColumns; ++xx) << 716 for( G4int xx = 0; xx < fColumns; xx++) { 789 mean = fTab[ww][xx]; 717 mean = fTab[ww][xx]; 790 density = Pixel2density(mean); 718 density = Pixel2density(mean); 791 unsigned int mateID = GetMaterialIndex << 719 unsigned int mateID = GetMaterialIndex( density ); 792 rflag = std::fwrite(&mateID, sizeof(un 720 rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut); 793 } 721 } 794 } 722 } 795 } << 723 796 else { << 724 } else { 797 // density value is the average of a squar 725 // density value is the average of a square region of 798 // fCompression*fCompression pixels 726 // fCompression*fCompression pixels 799 for (G4int ww = 0; ww < fRows; ww += compS << 727 for(G4int ww = 0; ww < fRows ;ww += compSize ) { 800 for (G4int xx = 0; xx < fColumns; xx += << 728 for(G4int xx = 0; xx < fColumns ;xx +=compSize ) { 801 overflow = false; 729 overflow = false; 802 mean = 0; 730 mean = 0; 803 for (G4int sumx = 0; sumx < compSize; << 731 for(int sumx = 0; sumx < compSize; sumx++) { 804 for (G4int sumy = 0; sumy < compSize << 732 for(int sumy = 0; sumy < compSize; sumy++) { 805 if (ww + sumy >= fRows || xx + sum << 733 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true; 806 mean += fTab[ww + sumy][xx + sumx] << 734 mean += fTab[ww+sumy][xx+sumx]; 807 } 735 } 808 if (overflow) break; << 736 if(overflow) break; 809 } 737 } 810 mean /= compSize * compSize; << 738 mean /= compSize*compSize; 811 << 739 812 if (!overflow) { << 740 if(!overflow) { 813 density = Pixel2density(mean); 741 density = Pixel2density(mean); 814 unsigned int mateID = GetMaterialInd << 742 unsigned int mateID = GetMaterialIndex( density ); 815 rflag = std::fwrite(&mateID, sizeof( 743 rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut); 816 } 744 } 817 } 745 } >> 746 818 } 747 } 819 } 748 } 820 << 749 821 //----- Write density for each pixel 750 //----- Write density for each pixel 822 if (compSize == 1) { // no fCompression: ea << 751 if(compSize == 1) { // no fCompression: each pixel has a density value) 823 for (G4int ww = 0; ww < fRows; ++ww) { << 752 for( G4int ww = 0; ww < fRows; ww++) { 824 for (G4int xx = 0; xx < fColumns; ++xx) << 753 for( G4int xx = 0; xx < fColumns; xx++) { 825 mean = fTab[ww][xx]; 754 mean = fTab[ww][xx]; 826 density = Pixel2density(mean); 755 density = Pixel2density(mean); 827 rflag = std::fwrite(&density, sizeof(G 756 rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut); 828 } 757 } 829 } 758 } 830 } << 759 831 else { << 760 } else { 832 // density value is the average of a squar 761 // density value is the average of a square region of 833 // fCompression*fCompression pixels 762 // fCompression*fCompression pixels 834 for (G4int ww = 0; ww < fRows; ww += compS << 763 for(G4int ww = 0; ww < fRows ;ww += compSize ) { 835 for (G4int xx = 0; xx < fColumns; xx += << 764 for(G4int xx = 0; xx < fColumns ;xx +=compSize ) { 836 overflow = false; 765 overflow = false; 837 mean = 0; 766 mean = 0; 838 for (G4int sumx = 0; sumx < compSize; << 767 for(int sumx = 0; sumx < compSize; sumx++) { 839 for (G4int sumy = 0; sumy < compSize << 768 for(int sumy = 0; sumy < compSize; sumy++) { 840 if (ww + sumy >= fRows || xx + sum << 769 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true; 841 mean += fTab[ww + sumy][xx + sumx] << 770 mean += fTab[ww+sumy][xx+sumx]; 842 } 771 } 843 if (overflow) break; << 772 if(overflow) break; 844 } 773 } 845 mean /= compSize * compSize; << 774 mean /= compSize*compSize; 846 << 775 847 if (!overflow) { << 776 if(!overflow) { 848 density = Pixel2density(mean); 777 density = Pixel2density(mean); 849 rflag = std::fwrite(&density, sizeof 778 rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut); 850 } 779 } 851 } 780 } >> 781 852 } 782 } 853 } 783 } 854 << 784 855 rflag = std::fclose(fileOut); 785 rflag = std::fclose(fileOut); 856 << 786 857 delete[] nameProcessed; << 787 delete [] nameProcessed; 858 << 788 859 /* for ( G4int i = 0; i < fRows; i ++ ) { 789 /* for ( G4int i = 0; i < fRows; i ++ ) { 860 delete [] fTab[i]; 790 delete [] fTab[i]; 861 } 791 } 862 delete [] fTab; 792 delete [] fTab; 863 */ 793 */ 864 << 794 865 if (rflag) return returnvalue; 795 if (rflag) return returnvalue; 866 return returnvalue; 796 return returnvalue; 867 } 797 } 868 798 869 //....oooOO0OOooo........oooOO0OOooo........oo << 870 << 871 // DICOM HEAD does not use the calibration cur << 872 799 873 #ifdef DICOM_USE_HEAD << 800 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 874 void DicomHandler::ReadCalibration() << 875 { << 876 fNbrequali = 0; << 877 fReadCalibration = false; << 878 G4cout << "No calibration curve for the DICO << 879 } << 880 #else << 881 // Separated out of Pixel2density 801 // Separated out of Pixel2density 882 // No need to read in same calibration EVERY t 802 // No need to read in same calibration EVERY time 883 // Increases the speed of reading file by seve 803 // Increases the speed of reading file by several orders of magnitude 884 << 885 void DicomHandler::ReadCalibration() 804 void DicomHandler::ReadCalibration() 886 { 805 { 887 fNbrequali = 0; << 806 nbrequali = 0; >> 807 888 // CT2Density.dat contains the calibration c 808 // CT2Density.dat contains the calibration curve to convert CT (Hounsfield) 889 // number to physical density 809 // number to physical density 890 std::ifstream calibration(fCt2DensityFile.c_ << 810 std::ifstream calibration(ct2densityFile.c_str()); 891 calibration >> fNbrequali; << 811 calibration >> nbrequali; 892 fValueDensity = new G4double[fNbrequali]; << 812 893 fValueCT = new G4double[fNbrequali]; << 813 valuedensity = new G4double[nbrequali]; 894 << 814 valueCT = new G4double[nbrequali]; 895 if (!calibration) { << 815 896 G4Exception("DicomHandler::ReadFile", "DIC << 816 if(!calibration) { >> 817 G4Exception("DicomHandler::ReadFile", >> 818 "DICOM001", >> 819 FatalException, 897 "@@@ No value to transform pix 820 "@@@ No value to transform pixels in density!"); 898 } << 821 899 else { // calibration was successfully open << 822 } else { // calibration was successfully opened 900 for (G4int i = 0; i < fNbrequali; ++i) { << 823 for(G4int i = 0; i < nbrequali; i++) { // Loop to store all the 901 calibration >> fValueCT[i] >> fValueDens << 824 //pts in CT2Density.dat >> 825 calibration >> valueCT[i] >> valuedensity[i]; 902 } 826 } 903 } 827 } 904 calibration.close(); 828 calibration.close(); 905 fReadCalibration = true; << 829 906 } << 830 readCalibration = true; 907 #endif << 908 << 909 #ifdef DICOM_USE_HEAD << 910 G4float DicomHandler::Pixel2density(G4int pixe << 911 { << 912 G4float density = -1; << 913 << 914 // Air << 915 if (pixel == -998.) density = 0.001290; << 916 // Soft Tissue << 917 else if (pixel == 24.) << 918 density = 1.00; << 919 // Brain << 920 else if (pixel == 52.) << 921 density = 1.03; << 922 // Spinal disc << 923 else if (pixel == 92.) << 924 density = 1.10; << 925 // Trabecular bone << 926 else if (pixel == 197.) << 927 density = 1.18; << 928 // Cortical Bone << 929 else if (pixel == 923.) << 930 density = 1.85; << 931 // Tooth dentine << 932 else if (pixel == 1280.) << 933 density = 2.14; << 934 // Tooth enamel << 935 else if (pixel == 2310.) << 936 density = 2.89; << 937 << 938 return density; << 939 } 831 } 940 832 941 #else << 833 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 942 //....oooOO0OOooo........oooOO0OOooo........oo << 943 834 944 G4float DicomHandler::Pixel2density(G4int pixe 835 G4float DicomHandler::Pixel2density(G4int pixel) 945 { 836 { 946 if (!fReadCalibration) { << 837 if(!readCalibration) { ReadCalibration(); } 947 ReadCalibration(); << 838 948 } << 949 << 950 G4float density = -1.; 839 G4float density = -1.; 951 G4double deltaCT = 0; 840 G4double deltaCT = 0; 952 G4double deltaDensity = 0; 841 G4double deltaDensity = 0; 953 << 842 954 for (G4int j = 1; j < fNbrequali; ++j) { << 843 955 if (pixel >= fValueCT[j - 1] && pixel < fV << 844 for(G4int j = 1; j < nbrequali; j++) { 956 deltaCT = fValueCT[j] - fValueCT[j - 1]; << 845 if( pixel >= valueCT[j-1] && pixel < valueCT[j]) { 957 deltaDensity = fValueDensity[j] - fValue << 846 958 << 847 deltaCT = valueCT[j] - valueCT[j-1]; >> 848 deltaDensity = valuedensity[j] - valuedensity[j-1]; >> 849 959 // interpolating linearly 850 // interpolating linearly 960 density = G4float(fValueDensity[j] - ((f << 851 density = valuedensity[j] - ((valueCT[j] - pixel)*deltaDensity/deltaCT ); 961 break; 852 break; 962 } 853 } 963 } 854 } 964 << 855 965 if (density < 0.) { << 856 if(density < 0.) { 966 std::printf( << 857 std::printf("@@@ Error density = %f && Pixel = %i \ 967 "@@@ Error density = %f && Pixel = %i \ << 858 (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel, 968 (0x%x) && deltaDensity/deltaCT = %f\n", << 859 deltaDensity/deltaCT); 969 density, pixel, pixel, deltaDensity / de << 970 } 860 } 971 << 861 972 return density; 862 return density; 973 } 863 } 974 #endif << 864 975 //....oooOO0OOooo........oooOO0OOooo........oo 865 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 976 866 977 void DicomHandler::CheckFileFormat() 867 void DicomHandler::CheckFileFormat() 978 { 868 { 979 std::ifstream checkData(fDriverFile.c_str()) << 869 std::ifstream checkData(driverFile.c_str()); 980 char* oneLine = new char[128]; << 870 char * oneLine = new char[128]; 981 << 871 982 if (!(checkData.is_open())) { // Check exis << 872 if(!(checkData.is_open())) { //Check existance of Data.dat 983 << 873 984 G4String message = "\nDicomG4 needs Data.d << 874 G4String message = >> 875 "\nDicomG4 needs Data.dat (or another driver file specified"; 985 message += " in command line):\n"; 876 message += " in command line):\n"; 986 message += "\tFirst line: number of image 877 message += "\tFirst line: number of image pixel for a voxel (G4Box)\n"; 987 message += "\tSecond line: number of image 878 message += "\tSecond line: number of images (CT slices) to read\n"; 988 message += "\tEach following line contains 879 message += "\tEach following line contains the name of a Dicom image"; 989 message += " except for the .dcm extension 880 message += " except for the .dcm extension"; 990 G4Exception("DicomHandler::ReadFile", "DIC << 881 G4Exception("DicomHandler::ReadFile", >> 882 "DICOM001", >> 883 FatalException, >> 884 message.c_str()); 991 } 885 } 992 << 886 993 checkData >> fCompression; 887 checkData >> fCompression; 994 checkData >> fNFiles; 888 checkData >> fNFiles; 995 G4String oneName; 889 G4String oneName; 996 checkData.getline(oneLine, 100); << 890 checkData.getline(oneLine,100); 997 std::ifstream testExistence; 891 std::ifstream testExistence; 998 G4bool existAlready = true; 892 G4bool existAlready = true; 999 for (G4int rep = 0; rep < fNFiles; ++rep) { << 893 for(G4int rep = 0; rep < fNFiles; rep++) { 1000 checkData.getline(oneLine, 100); << 894 checkData.getline(oneLine,100); 1001 oneName = oneLine; 895 oneName = oneLine; 1002 oneName += ".g4dcm"; // create dicomFile << 896 oneName += ".g4dcm"; // create dicomFile.g4dcm 1003 G4cout << fNFiles << " test file " << one 897 G4cout << fNFiles << " test file " << oneName << G4endl; 1004 testExistence.open(oneName.data()); 898 testExistence.open(oneName.data()); 1005 if (!(testExistence.is_open())) { << 899 if(!(testExistence.is_open())) { 1006 existAlready = false; 900 existAlready = false; 1007 testExistence.clear(); 901 testExistence.clear(); 1008 testExistence.close(); 902 testExistence.close(); 1009 } 903 } 1010 testExistence.clear(); 904 testExistence.clear(); 1011 testExistence.close(); 905 testExistence.close(); 1012 } 906 } 1013 << 907 1014 ReadMaterialIndices(checkData); << 908 ReadMaterialIndices( checkData ); 1015 << 909 1016 checkData.close(); 910 checkData.close(); 1017 delete[] oneLine; << 911 delete [] oneLine; 1018 << 912 1019 if (existAlready == false) { // The files << 913 if( existAlready == false ) { // The files *.g4dcm have to be created 1020 << 914 1021 G4cout << "\nAll the necessary images wer 915 G4cout << "\nAll the necessary images were not found in processed form " 1022 << ", starting with .dcm images\n" 916 << ", starting with .dcm images\n"; 1023 << 917 1024 FILE* dicom; << 918 FILE * dicom; 1025 FILE* lecturePref; << 919 FILE * lecturePref; 1026 char* fCompressionc = new char[LINEBUFFSI << 920 char * fCompressionc = new char[LINEBUFFSIZE]; 1027 char* maxc = new char[LINEBUFFSIZE]; << 921 char * maxc = new char[LINEBUFFSIZE]; 1028 // char name[300], inputFile[300]; << 922 //char name[300], inputFile[300]; 1029 char* inputFile = new char[FILENAMESIZE]; << 923 char * name = new char[FILENAMESIZE]; >> 924 char * inputFile = new char[FILENAMESIZE]; 1030 G4int rflag; 925 G4int rflag; 1031 lecturePref = std::fopen(fDriverFile.c_st << 926 1032 << 927 lecturePref = std::fopen(driverFile.c_str(),"r"); 1033 rflag = std::fscanf(lecturePref, "%s", fC << 928 rflag = std::fscanf(lecturePref,"%s",fCompressionc); 1034 fCompression = atoi(fCompressionc); 929 fCompression = atoi(fCompressionc); 1035 rflag = std::fscanf(lecturePref, "%s", ma << 930 rflag = std::fscanf(lecturePref,"%s",maxc); 1036 fNFiles = atoi(maxc); 931 fNFiles = atoi(maxc); 1037 G4cout << " fNFiles " << fNFiles << G4end 932 G4cout << " fNFiles " << fNFiles << G4endl; 1038 << 933 1039 ///////////////////// << 934 for( G4int i = 1; i <= fNFiles; i++ ) { // Begin loop on filenames 1040 << 935 1041 #ifdef DICOM_USE_HEAD << 936 rflag = std::fscanf(lecturePref,"%s",inputFile); 1042 for (G4int i = 1; i <= fNFiles; ++i) { / << 937 std::sprintf(name,"%s.dcm",inputFile); 1043 rflag = std::fscanf(lecturePref, "%s", << 938 G4cout << "check 1: " << name << G4endl; 1044 G4String path = GetDicomDataPath() + "/ << 1045 << 1046 G4String name = inputFile + G4String(". << 1047 // Writes the results to a character st << 1048 << 1049 G4String name2 = path + name; << 1050 // Open input file and give it to gest 939 // Open input file and give it to gestion_dicom : 1051 G4cout << "### Opening " << name2 << " << 940 std::printf("### Opening %s and reading :\n",name); 1052 dicom = std::fopen(name2.c_str(), "rb") << 941 dicom = std::fopen(name,"rb"); 1053 // Reading the .dcm in two steps: 942 // Reading the .dcm in two steps: 1054 // 1. reading the header 943 // 1. reading the header 1055 // 2. reading the pixel data and s 944 // 2. reading the pixel data and store the density in Moyenne.dat 1056 if (dicom != 0) { << 945 if( dicom != 0 ) { 1057 ReadFile(dicom, inputFile); << 946 ReadFile(dicom,inputFile); 1058 } << 947 } else { 1059 else { << 1060 G4cout << "\nError opening file : " < << 1061 } << 1062 rflag = std::fclose(dicom); << 1063 } << 1064 #else << 1065 << 1066 for (G4int i = 1; i <= fNFiles; ++i) { / << 1067 rflag = std::fscanf(lecturePref, "%s", << 1068 << 1069 G4String name = inputFile + G4String(". << 1070 // Writes the results to a character st << 1071 << 1072 // G4cout << "check: " << name << G4end << 1073 // Open input file and give it to gest << 1074 G4cout << "### Opening " << name << " a << 1075 dicom = std::fopen(name.c_str(), "rb"); << 1076 // Reading the .dcm in two steps: << 1077 // 1. reading the header << 1078 // 2. reading the pixel data and s << 1079 if (dicom != 0) { << 1080 ReadFile(dicom, inputFile); << 1081 } << 1082 else { << 1083 G4cout << "\nError opening file : " < 948 G4cout << "\nError opening file : " << name << G4endl; 1084 } 949 } 1085 rflag = std::fclose(dicom); 950 rflag = std::fclose(dicom); 1086 } 951 } 1087 #endif << 1088 << 1089 rflag = std::fclose(lecturePref); 952 rflag = std::fclose(lecturePref); 1090 << 953 1091 // Checks the spacing is correct. Dumps t 954 // Checks the spacing is correct. Dumps to files 1092 fMergedSlices->CheckSlices(); << 955 mergedSlices->CheckSlices(); 1093 << 956 1094 delete[] fCompressionc; << 957 delete [] fCompressionc; 1095 delete[] maxc; << 958 delete [] maxc; 1096 delete[] inputFile; << 959 delete [] name; >> 960 delete [] inputFile; 1097 if (rflag) return; 961 if (rflag) return; >> 962 1098 } 963 } 1099 << 964 1100 if (fValueDensity) { << 965 if(valuedensity) { delete [] valuedensity; } 1101 delete[] fValueDensity; << 966 if(valueCT) { delete [] valueCT; } 1102 } << 967 if(mergedSlices) { delete mergedSlices; } 1103 if (fValueCT) { << 968 1104 delete[] fValueCT; << 1105 } << 1106 if (fMergedSlices) { << 1107 delete fMergedSlices; << 1108 } << 1109 } 969 } 1110 970 1111 //....oooOO0OOooo........oooOO0OOooo........o << 971 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1112 972 1113 G4int DicomHandler::read_defined_nested(FILE* << 973 G4int DicomHandler::read_defined_nested(FILE * nested,G4int SQ_Length) 1114 { 974 { 1115 // VARIABLES 975 // VARIABLES 1116 unsigned short item_GroupNumber; 976 unsigned short item_GroupNumber; 1117 unsigned short item_ElementNumber; 977 unsigned short item_ElementNumber; 1118 G4int item_Length; 978 G4int item_Length; 1119 G4int items_array_length = 0; << 979 G4int items_array_length=0; 1120 char* buffer = new char[LINEBUFFSIZE]; << 980 char * buffer= new char[LINEBUFFSIZE]; 1121 size_t rflag = 0; 981 size_t rflag = 0; 1122 << 982 1123 while (items_array_length < SQ_Length) { << 983 while(items_array_length < SQ_Length) 1124 rflag = std::fread(buffer, 2, 1, nested); << 984 { 1125 GetValue(buffer, item_GroupNumber); << 985 rflag = std::fread(buffer, 2, 1, nested); 1126 << 986 GetValue(buffer, item_GroupNumber); 1127 rflag = std::fread(buffer, 2, 1, nested); << 987 1128 GetValue(buffer, item_ElementNumber); << 988 rflag = std::fread(buffer, 2, 1, nested); 1129 << 989 GetValue(buffer, item_ElementNumber); 1130 rflag = std::fread(buffer, 4, 1, nested); << 990 1131 GetValue(buffer, item_Length); << 991 rflag = std::fread(buffer, 4, 1, nested); 1132 << 992 GetValue(buffer, item_Length); 1133 rflag = std::fread(buffer, item_Length, 1 << 993 1134 << 994 rflag = std::fread(buffer, item_Length, 1, nested); 1135 items_array_length = items_array_length + << 995 1136 } << 996 items_array_length= items_array_length+8+item_Length; 1137 << 997 } 1138 delete[] buffer; << 998 1139 << 999 delete [] buffer; 1140 if (SQ_Length > items_array_length) << 1000 >> 1001 if( SQ_Length>items_array_length ) 1141 return 0; 1002 return 0; 1142 else 1003 else 1143 return 1; 1004 return 1; 1144 if (rflag) return 1; 1005 if (rflag) return 1; 1145 } 1006 } 1146 1007 1147 //....oooOO0OOooo........oooOO0OOooo........o << 1008 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1148 1009 1149 void DicomHandler::read_undefined_nested(FILE << 1010 void DicomHandler::read_undefined_nested(FILE * nested) 1150 { 1011 { 1151 // VARIABLES 1012 // VARIABLES 1152 unsigned short item_GroupNumber; 1013 unsigned short item_GroupNumber; 1153 unsigned short item_ElementNumber; 1014 unsigned short item_ElementNumber; 1154 unsigned int item_Length; 1015 unsigned int item_Length; 1155 char* buffer = new char[LINEBUFFSIZE]; << 1016 char * buffer= new char[LINEBUFFSIZE]; 1156 size_t rflag = 0; 1017 size_t rflag = 0; 1157 << 1018 1158 do { << 1019 do 1159 rflag = std::fread(buffer, 2, 1, nested); << 1020 { 1160 GetValue(buffer, item_GroupNumber); << 1021 rflag = std::fread(buffer, 2, 1, nested); 1161 << 1022 GetValue(buffer, item_GroupNumber); 1162 rflag = std::fread(buffer, 2, 1, nested); << 1023 1163 GetValue(buffer, item_ElementNumber); << 1024 rflag = std::fread(buffer, 2, 1, nested); 1164 << 1025 GetValue(buffer, item_ElementNumber); 1165 rflag = std::fread(buffer, 4, 1, nested); << 1026 1166 GetValue(buffer, item_Length); << 1027 rflag = std::fread(buffer, 4, 1, nested); 1167 << 1028 GetValue(buffer, item_Length); 1168 if (item_Length != 0xffffffff) << 1029 1169 rflag = std::fread(buffer, item_Length, << 1030 if(item_Length!=0xffffffff) 1170 else << 1031 rflag = std::fread(buffer, item_Length, 1, nested); 1171 read_undefined_item(nested); << 1032 else 1172 << 1033 read_undefined_item(nested); 1173 } while (item_GroupNumber != 0xFFFE || item << 1034 1174 << 1035 1175 delete[] buffer; << 1036 } while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD 1176 if (rflag) return; << 1037 || item_Length!=0); >> 1038 >> 1039 delete [] buffer; >> 1040 if (rflag) return; 1177 } 1041 } 1178 1042 1179 //....oooOO0OOooo........oooOO0OOooo........o << 1043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1180 1044 1181 void DicomHandler::read_undefined_item(FILE* << 1045 void DicomHandler::read_undefined_item(FILE * nested) 1182 { 1046 { 1183 // VARIABLES 1047 // VARIABLES 1184 unsigned short item_GroupNumber; 1048 unsigned short item_GroupNumber; 1185 unsigned short item_ElementNumber; 1049 unsigned short item_ElementNumber; 1186 G4int item_Length; << 1050 G4int item_Length; size_t rflag = 0; 1187 size_t rflag = 0; << 1051 char *buffer= new char[LINEBUFFSIZE]; 1188 char* buffer = new char[LINEBUFFSIZE]; << 1052 1189 << 1053 do 1190 do { << 1054 { 1191 rflag = std::fread(buffer, 2, 1, nested); << 1055 rflag = std::fread(buffer, 2, 1, nested); 1192 GetValue(buffer, item_GroupNumber); << 1056 GetValue(buffer, item_GroupNumber); 1193 << 1057 1194 rflag = std::fread(buffer, 2, 1, nested); << 1058 rflag = std::fread(buffer, 2, 1, nested); 1195 GetValue(buffer, item_ElementNumber); << 1059 GetValue(buffer, item_ElementNumber); 1196 << 1060 1197 rflag = std::fread(buffer, 4, 1, nested); << 1061 rflag = std::fread(buffer, 4, 1, nested); 1198 GetValue(buffer, item_Length); << 1062 GetValue(buffer, item_Length); 1199 << 1063 1200 if (item_Length != 0) rflag = std::fread( << 1064 1201 << 1065 if(item_Length!=0) 1202 } while (item_GroupNumber != 0xFFFE || item << 1066 rflag = std::fread(buffer,item_Length,1,nested); 1203 << 1067 1204 delete[] buffer; << 1068 } >> 1069 while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D >> 1070 || item_Length!=0); >> 1071 >> 1072 delete [] buffer; 1205 if (rflag) return; 1073 if (rflag) return; 1206 } 1074 } 1207 1075 1208 //....oooOO0OOooo........oooOO0OOooo........o << 1076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1209 1077 1210 template<class Type> << 1078 template <class Type> 1211 void DicomHandler::GetValue(char* _val, Type& << 1079 void DicomHandler::GetValue(char * _val, Type & _rval) { 1212 { << 1080 1213 #if BYTE_ORDER == BIG_ENDIAN 1081 #if BYTE_ORDER == BIG_ENDIAN 1214 if (fLittleEndian) { // little endian << 1082 if(fLittleEndian) { // little endian 1215 #else // BYTE_ORDER == LITTLE_ENDIAN << 1083 #else // BYTE_ORDER == LITTLE_ENDIAN 1216 if (!fLittleEndian) { // big endian << 1084 if(!fLittleEndian) { // big endian 1217 #endif 1085 #endif 1218 const G4int SIZE = sizeof(_rval); << 1086 const int SIZE = sizeof(_rval); 1219 char ctemp; << 1087 char ctemp; 1220 for (G4int i = 0; i < SIZE / 2; ++i) { << 1088 for(int i = 0; i < SIZE/2; i++) { 1221 ctemp = _val[i]; << 1089 ctemp = _val[i]; 1222 _val[i] = _val[SIZE - 1 - i]; << 1090 _val[i] = _val[SIZE - 1 - i]; 1223 _val[SIZE - 1 - i] = ctemp; << 1091 _val[SIZE - 1 - i] = ctemp; >> 1092 } 1224 } 1093 } >> 1094 _rval = *(Type *)_val; 1225 } 1095 } 1226 _rval = *(Type*)_val; << 1096 1227 } << 1097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1228 << 1098 1229 //....oooOO0OOooo........oooOO0OOooo........o << 1230 1099