1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2006/07/24 16:36:45 $
7 Version: $Revision: 1.3 $
9 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
10 l'Image). All rights reserved. See Doc/License.txt or
11 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
13 This software is distributed WITHOUT ANY WARRANTY; without even
14 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 PURPOSE. See the above copyright notices for more information.
17 =========================================================================*/
24 #include "gdcmFileHelper.h"
25 #include "gdcmDebug.h"
26 #include "gdcmDirList.h"
29 #include "gdcmArgMgr.h"
33 * - explores recursively the (single Patient, single Study) directory
34 * - examines the ".txt" Dense multiframe files
35 * - Converts the files into 16 bits Dicom Files
36 * WARNING : directory must contain ONLY .txt files
40 void Load(std::ifstream &from, std::string imageName, std::string patName,
41 std::string strStudyUID, int serieNumber);
43 //std::ifstream& eatwhite(std::ifstream& is);
46 window center (level) and width, are defined in the DICOM
47 standard very precisely as follows (see PS 3.3 C.11.2.1.2):
49 >"These Attributes are applied according to the following pseudo-code,
52 y is an output value with a range from ymin to ymax,
53 c is Window Center (0028,1050)
54 w is Window Width (0028,1051):
56 > if (x <= c - 0.5 - (w-1)/2), then y = ymin
57 > else if (x > c - 0.5 + (w-1)/2), then y = ymax,
58 > else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin
62 From: David Clunie - view profile
63 Date: Thurs, Jun 1 2006 3:03 pm
64 Email: David Clunie <dclu...@dclunie.com>
65 Groups: comp.protocols.dicom
67 The value of x is the output of the preceding step, the so-called
68 "modality LUT", which may either be:
70 - identity (no rescale values or Modality LUT, or the SOP Class is
71 PET and rescale values are ignored), in which case x is the stored
72 pixel value in (7FE0,0010)
74 - Rescale Slope and Intercept (as typically used in CT), in which
75 case x is the value obtained from applying the rescale values to
76 the stored pixel value
78 - an actual LUT, in which case x is the value stored in the LUT
79 corresponding to the LUT index value that is the stored pixel
82 The ymin and ymax are intended to represent the output range; for
83 example, if the hypothetical Presentation LUT step that follows
84 the VOI LUT (window) stage is an identity operation, then the
85 ymin and ymax represent P-Values, which might be the range of
86 digital driving levels for your display (calibrated to the GSDF),
87 in the 8-bit wide output case ranging from 0 to 255, for example.
89 The terms brightness and contrast are not used in radiology imaging
90 -the window center and width are used instead.
93 int main(int argc, char *argv[])
96 " \n DenseMultiframeToDicom :\n ",
97 " - explores recursively the given (single Patient, single Study) directory",
98 " - examines the '.txt' files ",
99 " - Converts the files into 16 bits Dicom files, ",
100 " WARNING : directory must contain ONLY .txt files ",
102 " DenseMultiframeToDicom dirin=rootDirectoryName ",
103 " [studyUID = ] [patName = ] ",
104 " [listonly] [verbose] [debug] ",
106 "studyUID : *aware* user wants to add the serie ",
107 " to an already existing study ",
108 " verbose : user wants to run the program in 'verbose mode' ",
109 " debug : *developer* wants to run the program in 'debug mode' ",
112 // ----- Initialize Arguments Manager ------
114 gdcm::ArgMgr *am = new gdcm::ArgMgr(argc, argv);
116 if (argc == 1 || am->ArgMgrDefined("usage"))
118 am->ArgMgrUsage(usage); // Display 'usage'
123 const char *dirNamein;
124 dirNamein = am->ArgMgrGetString("dirin",".");
126 if (am->ArgMgrDefined("debug"))
127 gdcm::Debug::DebugOn();
129 int verbose = am->ArgMgrDefined("verbose");
130 int listonly = am->ArgMgrDefined("listonly");
131 std::string patName = am->ArgMgrGetString("patname", "g^PatientName");
133 bool userDefinedStudy = am->ArgMgrDefined("studyUID");
134 const char *studyUID = am->ArgMgrGetString("studyUID");
136 // not described *on purpose* in the Usage !
137 bool userDefinedSerie = am->ArgMgrDefined("serieUID");
138 const char *serieUID = am->ArgMgrGetString("serieUID");
141 // if unused Param we give up
142 if ( am->ArgMgrPrintUnusedLabels() )
144 am->ArgMgrUsage(usage);
148 delete am; // we don't need Argument Manager any longer
150 // ----- Begin Processing -----
152 if ( ! gdcm::DirList::IsDirectory(dirNamein) )
154 std::cout << "KO : [" << dirNamein << "] is not a Directory."
160 std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
163 std::string strDirNamein(dirNamein);
164 gdcm::DirList dirList(strDirNamein, true); // (recursively) the list of files
168 std::cout << "------------List of found files ------------" << std::endl;
170 std::cout << std::endl;
174 std::string filenameout;
176 std::string strStudyUID;
177 strStudyUID = gdcm::Util::CreateUniqueUID();
179 gdcm::DirListType fileNames;
180 fileNames = dirList.GetFilenames();
181 for (gdcm::DirListType::iterator it = fileNames.begin();
182 it != fileNames.end();
185 std::ifstream from( (*it).c_str() );
188 std::cout << "Can't open file" << *it << std::endl;
193 std::cout << "Success in open file" << *it << std::endl;
194 Load(from, *it, patName, strStudyUID, serieNumber);
202 void Load(std::ifstream &from, std::string imageName, std::string patName,
203 std::string strStudyUID, int serieNumber)
208 /* was OK for single frames
210 ---------------------------
211 Array dimensions = 58 x 56
212 The following is the 2D array of peak circumferential strain map, all zero value
213 pixels are outside the mask
214 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
215 -----------------------------
226 std::cout << "nx " << nx << " ny " << ny << std::endl;
229 std::getline(from, str1);
230 std::cout << "[" << str1 << "]" << std::endl;
231 std::getline(from, str1);
232 std::cout << "[" << str1 << "]" << std::endl;
233 std::getline(from, str1);
234 std::cout << "[" << str1 << "]" << std::endl;
237 /* now, we deal with multiframes
239 ------------------------------------------
240 X dim, Y dim, N of frames = 52x59x14
241 Temporal resolution = 31.9200 ms
242 First frame starts at 47.9600 ms
243 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
244 All pixels with zero strain values are outside the masks.
245 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
246 --------------------------------------------
250 from >> str1; // X dim,
252 from >> str1; // Y dim,
254 from >> str1; // N of frames =
258 from >> str1; // 52x59x14
260 std::cout << "[" << str1 << "]" << std::endl;
262 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
263 std::cout << nx << " " << ny << " " << nf << std::endl;
266 for (int k=0; k<5; k++) // 5
267 { std::cout << "Comment line number : " << k << std::endl;
268 std::getline(from, str1);
269 std::cout << str1 << std::endl;
272 //float *f = new float(nx*ny);
273 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
274 // float mini = FLT_MAX, maxi = FLT_MIN;
278 std::string strSerieUID;
279 strSerieUID = gdcm::Util::CreateUniqueUID();
281 for (int l=0; l<nf; l++) // Loop on the frames
284 std::cout << "Frame nb " << l << std::endl;
285 for( int j=0; j<ny; j++)
287 std::cout << " "<< std::endl;
289 for (int i=0; i<nx; i++)
291 // std::cout << "--------------------------" << std::endl;
295 for (;;) //eatwhite(from);
299 std::cout << " !from.get(c) ";
304 //std::cout << " !isspace(c) ";
308 } //end eatwhite(from);
309 // std::cout << "------end eatwhite-----" << std::endl;
312 // trouble : whe space is missing "-0.0990263-8.8778"
313 // is not interpreted as TWO values :-(
320 //std::cout << "[" << c << "]" << std::endl;
323 if ( first != 1 && c == '-' && previous != 'e')
326 std::cout << " One more gauffre in frame:" << std::dec << l
327 << ", line : " << j << " element " << i << std::endl;
336 val = (float)atof(str1.c_str());
337 std::cout << " " << val;
338 *(f+ /*l*nx*ny + */j*nx+i) = val;
342 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
349 //std::cout << std::endl;
350 //std::cout << std::endl << " line nb : " << j
351 // << " line length : " << l << std::endl;
355 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
357 // values are expressed as % as a fraction, actually!)
358 // It's useless to rescale them as uint16_t : just *100
359 uint16_t *img = new uint16_t[ny*nx];
362 float div = maxi-mini;
363 std::cout << "div : " << div << " mini : " << mini << std::endl;
364 for( int k=0; k<ny*nx; k++)
366 *ptr = (uint16_t)((*tmp * 65535.0) / div);
371 int16_t *img = new int16_t[ny*nx];
373 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
374 for( int k=0; k<ny*nx; k++)
376 if(*tmp > 1.0) // artefacted pixel
378 else /// \todo FIXME : what about max threshold ?
379 *ptr = (int16_t)(*tmp *100);
381 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
386 // gdcm::Debug::DebugOn();
388 std::ostringstream str;
390 file = gdcm::File::New();
392 // Set the image size
395 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
398 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
400 // Set the pixel type
402 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
404 str << 16; // may be 12 or 16 if componentSize =16
405 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
407 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
409 // Set the pixel representation // 0/1 1 : signed
410 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
412 // Set the samples per pixel // 1:Grey level, 3:RGB
413 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
416 // Set Rescale Intercept
419 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
424 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
427 // 0020 0037 DS 6 Image Orientation (Patient)
428 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
430 // 0020 0032 DS 3 Image Position (Patient)
431 char charImagePosition[256];
432 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
434 // 0020 0x1041 DS 1 Slice Location
435 sprintf(charImagePosition,"%f",float(l%nf));
436 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
438 //0008 103e LO 1 Series Description
439 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
441 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
442 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
443 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
445 //0020 0011 "IS" Series Number
446 sprintf(charImagePosition,"%04d",serieNumber);
447 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
451 gdcm::FileHelper *fh;
452 fh = gdcm::FileHelper::New(file);
453 // cast is just to avoid warnings (*no* conversion)
454 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
455 fh->SetWriteModeToRaw();
456 fh->SetWriteTypeToDcmExplVR();
458 fh->SetWriteTypeToDcmExplVR();
461 sprintf(numero, "%02d", l);
462 std::string fileName = imageName + "." + numero + ".dcm";
463 std::cout << "fileName " << fileName << std::endl;
465 if( !fh->Write(fileName))
466 std::cout << "Failed for [" << fileName << "]\n"
467 << " File is unwrittable" << std::endl;
470 } // end loop on frames
472 int pos = (long)(from.tellg());
473 std::cout << std::hex << "==============================================="
477 // Anatomical Images.
478 std::cout << "Anatomical images" << std::endl;
480 strSerieUID = gdcm::Util::CreateUniqueUID();
482 for (int l=0; l<nf; l++) // Loop on the frames
484 std::cout << "Frame nb " << l << std::endl;
485 for( int j=0; j<ny; j++)
488 for (int i=0; i<nx; i++)
503 val = (float)atof(str1.c_str());
504 // std::cout << " " << val;
505 *(f+ /*l*nx*ny + */j*nx+i) = val;
509 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
517 uint16_t *img = new uint16_t[ny*nx];
519 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
520 for( int k=0; k<ny*nx; k++)
522 *ptr = (int16_t)(*tmp *100);
527 std::ostringstream str;
529 file = gdcm::File::New();
531 // Set the image size
534 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
537 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
539 // Set the pixel type
541 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
543 str << 16; // may be 12 or 16 if componentSize =16
544 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
546 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
548 // Set the pixel representation // 0/1 1 : signed
549 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
551 // Set the samples per pixel // 1:Grey level, 3:RGB
552 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
555 // Set Rescale Intercept
558 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
563 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
566 // 0020 0037 DS 6 Image Orientation (Patient)
567 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
569 // 0020 0032 DS 3 Image Position (Patient)
570 char charImagePosition[256];
571 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
573 // 0020 0x1041 DS 1 Slice Location
574 sprintf(charImagePosition,"%f",float(l%nf));
575 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
577 //0008 103e LO 1 Series Description
578 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
580 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
581 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
582 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
584 //0020 0011 "IS" Series Number
585 sprintf(charImagePosition,"%04d",serieNumber+1);
586 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
590 gdcm::FileHelper *fh;
591 fh = gdcm::FileHelper::New(file);
592 // cast is just to avoid warnings (*no* conversion)
593 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
594 fh->SetWriteModeToRaw();
595 fh->SetWriteTypeToDcmExplVR();
597 fh->SetWriteTypeToDcmExplVR();
600 sprintf(numero, "%02d", l);
601 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
602 std::cout << "fileName " << fileName << std::endl;
604 if( !fh->Write(fileName))
605 std::cout << "Failed for [" << fileName << "]\n"
606 << " File is unwrittable" << std::endl;
612 } // end loop on frames