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