1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2007/06/27 08:38:44 $
7 Version: $Revision: 1.7 $
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
39 void Load(std::ifstream &from, std::string imageName, std::string patName,
40 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose);
42 //std::ifstream& eatwhite(std::ifstream& is);
45 window center (level) and width, are defined in the DICOM
46 standard very precisely as follows (see PS 3.3 C.11.2.1.2):
48 >"These Attributes are applied according to the following pseudo-code,
51 y is an output value with a range from ymin to ymax,
52 c is Window Center (0028,1050)
53 w is Window Width (0028,1051):
55 > if (x <= c - 0.5 - (w-1)/2), then y = ymin
56 > else if (x > c - 0.5 + (w-1)/2), then y = ymax,
57 > else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin
61 From: David Clunie - view profile
62 Date: Thurs, Jun 1 2006 3:03 pm
63 Email: David Clunie <dclu...@dclunie.com>
64 Groups: comp.protocols.dicom
66 The value of x is the output of the preceding step, the so-called
67 "modality LUT", which may either be:
69 - identity (no rescale values or Modality LUT, or the SOP Class is
70 PET and rescale values are ignored), in which case x is the stored
71 pixel value in (7FE0,0010)
73 - Rescale Slope and Intercept (as typically used in CT), in which
74 case x is the value obtained from applying the rescale values to
75 the stored pixel value
77 - an actual LUT, in which case x is the value stored in the LUT
78 corresponding to the LUT index value that is the stored pixel
81 The ymin and ymax are intended to represent the output range; for
82 example, if the hypothetical Presentation LUT step that follows
83 the VOI LUT (window) stage is an identity operation, then the
84 ymin and ymax represent P-Values, which might be the range of
85 digital driving levels for your display (calibrated to the GSDF),
86 in the 8-bit wide output case ranging from 0 to 255, for example.
88 The terms brightness and contrast are not used in radiology imaging
89 -the window center and width are used instead.
92 int main(int argc, char *argv[])
95 " \n DenseMultiFramesToDicom : \n",
96 " - explores recursively the given (single Patient, single Study) directory",
97 " - examines the '.txt' files ",
98 " - Converts the files into 16 bits Dicom files, ",
99 " WARNING : directory must contain ONLY .txt files ",
101 " DenseMultiFramesToDicom dirin=rootDirectoryName ",
102 " [studyUID = ] [patName = ] ",
103 " [listonly] [verbose] [debug] ",
105 " studyUID : *aware* user wants to add the serie ",
106 " to an already existing study ",
107 " verbose : user wants to run the program in 'verbose mode' ",
108 " debug : *developer* wants to run the program in 'debug mode' ",
111 // ----- Initialize Arguments Manager ------
113 GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
115 if (argc == 1 || am->ArgMgrDefined("usage"))
117 am->ArgMgrUsage(usage); // Display 'usage'
122 const char *dirNamein;
123 dirNamein = am->ArgMgrGetString("dirin",".");
125 if (am->ArgMgrDefined("debug"))
126 GDCM_NAME_SPACE::Debug::DebugOn();
128 int verbose = am->ArgMgrDefined("verbose");
129 int listonly = am->ArgMgrDefined("listonly");
130 std::string patName = am->ArgMgrGetString("patname", dirNamein);
132 bool userDefinedStudy = am->ArgMgrDefined("studyUID");
133 const char *studyUID;
134 if (userDefinedStudy)
135 studyUID = am->ArgMgrGetString("studyUID");
137 // not described *on purpose* in the Usage !
139 bool userDefinedSerie = am->ArgMgrDefined("serieUID");
140 const char *serieUID;
142 serieUID = am->ArgMgrGetString("serieUID");
144 // if unused Param we give up
145 if ( am->ArgMgrPrintUnusedLabels() )
147 am->ArgMgrUsage(usage);
151 delete am; // we don't need Argument Manager any longer
153 // ----- Begin Processing -----
155 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
157 std::cout << "KO : [" << dirNamein << "] is not a Directory."
164 std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
167 std::string strDirNamein(dirNamein);
168 GDCM_NAME_SPACE::DirList dirList(strDirNamein, true); // (recursively) the list of files
172 std::cout << "------------List of found files ------------" << std::endl;
174 std::cout << std::endl;
178 std::string filenameout;
183 std::string strStudyUID;
184 std::string strSerieUID;
186 if (userDefinedStudy)
187 strSerieUID = studyUID;
189 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
191 if (userDefinedStudy)
192 strSerieUID = serieUID;
194 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
198 GDCM_NAME_SPACE::DirListType fileNames;
199 fileNames = dirList.GetFilenames();
200 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
201 it != fileNames.end();
204 if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
209 int sz = (*it).size();
210 if ( (*it).c_str()[sz-1] != 't')
212 // skip non .txt files
215 std::ifstream from( (*it).c_str() );
218 std::cout << "Can't open file" << *it << std::endl;
224 std::cout << "Success in open file" << *it << std::endl;
225 Load(from, *it, patName, strStudyUID, strSerieUID, serieNumber, verbose);
233 void Load(std::ifstream &from, std::string imageName, std::string patName,
234 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose)
237 std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
241 std::cout << " ========= Create Parametric images" << std::endl;
242 /* was OK for single frames
244 ---------------------------
245 Array dimensions = 58 x 56
246 The following is the 2D array of peak circumferential strain map, all zero value
247 pixels are outside the mask
248 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
249 -----------------------------
260 std::cout << "nx " << nx << " ny " << ny << std::endl;
263 std::getline(from, str1);
264 std::cout << "[" << str1 << "]" << std::endl;
265 std::getline(from, str1);
266 std::cout << "[" << str1 << "]" << std::endl;
267 std::getline(from, str1);
268 std::cout << "[" << str1 << "]" << std::endl;
271 /* now, we deal with multiframes
273 ------------------------------------------
274 X dim, Y dim, N of frames = 52x59x14
275 Temporal resolution = 31.9200 ms
276 First frame starts at 47.9600 ms
277 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
278 All pixels with zero strain values are outside the masks.
279 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
280 --------------------------------------------
284 from >> str1; // X dim,
286 from >> str1; // Y dim,
288 from >> str1; // N of frames =
292 from >> str1; // 52x59x14
295 std::cout << "[" << str1 << "]" << std::endl;
297 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
298 std::cout << nx << " " << ny << " " << nf << std::endl;
300 std::getline(from, str1);
302 from >> str1; // Temporal
303 from >> str1; // Resolution
308 float temporalResolution;
309 sscanf( str1.c_str(),"%f",&temporalResolution);
311 std::cout << "temporal Resolution = " << temporalResolution << std::endl;
312 std::getline(from, str1);
314 from >> str1; // First
315 from >> str1; // frame
316 from >> str1; // starts
321 sscanf( str1.c_str(),"%f",&timeStart);
322 std::cout << "time Start = " << timeStart << std::endl;
323 std::getline(from, str1);
326 for (int k=0; k<2; k++) //
328 std::getline(from, str1);
329 std::cout << str1 << std::endl;
332 //float *f = new float(nx*ny);
333 // -->float *f = new float[nx*ny]; // Would be better !
334 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
335 // float mini = FLT_MAX, maxi = FLT_MIN;
340 currentTime = timeStart;
342 for (l1=0; l1<nf; l1++) // Loop on the frames
344 //std::cout << "Frame nb " << l1 << std::endl;
345 for( int j=0; j<ny; j++)
347 for (int i=0; i<nx; i++)
352 for (;;) //eatwhite(from);
356 std::cout << " !from.get(c) ";
361 //std::cout << " !isspace(c) ";
365 } //end eatwhite(from);
367 // trouble : when space is missing "-0.0990263-8.8778"
368 // is not interpreted as TWO values :-(
378 if ( first != 1 && c == '-' && previous != 'e')
381 //std::cout << " One more gauffre in frame:" << std::dec << l
382 // << ", line : " << j << " element " << i << std::endl;
391 val = (float)atof(str1.c_str());
392 //std::cout << " " << val;
393 *(f+ /*l*nx*ny + */j*nx+i) = val;
397 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
403 //std::cout << std::endl;
404 //std::cout << std::endl << " line nb : " << j
405 // << " line length : " << l << std::endl;
408 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
410 // values are expressed as % as a fraction, actually!)
411 // It's useless to rescale them as uint16_t : just *100
412 uint16_t *img = new uint16_t[ny*nx];
415 float div = maxi-mini;
416 std::cout << "div : " << div << " mini : " << mini << std::endl;
417 for( int k=0; k<ny*nx; k++)
419 *ptr = (uint16_t)((*tmp * 65535.0) / div);
424 int16_t *img = new int16_t[ny*nx];
426 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
427 for( int k=0; k<ny*nx; k++)
429 if(*tmp > 1.0) // artefacted pixel
431 else /// \todo FIXME : what about max threshold ?
432 *ptr = (int16_t)(*tmp *100);
434 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
439 // GDCM_NAME_SPACE::Debug::DebugOn();
441 std::ostringstream str;
442 GDCM_NAME_SPACE::File *file;
443 file = GDCM_NAME_SPACE::File::New();
445 // Set the image size
448 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
451 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
453 // Set the pixel type
455 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
457 str << 16; // may be 12 or 16 if componentSize =16
458 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
460 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
462 // Set the pixel representation // 0/1 1 : signed
463 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
465 // Set the samples per pixel // 1:Grey level, 3:RGB
466 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
469 // Set Rescale Intercept
472 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
477 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
480 // 0020 0037 DS 6 Image Orientation (Patient)
481 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
483 // 0020 0032 DS 3 Image Position (Patient)
484 char charImagePosition[256];
485 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
487 // 0020 0x1041 DS 1 Slice Location
488 sprintf(charImagePosition,"%f",float(l1%nf));
489 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
491 //0008 103e LO 1 Series Description
492 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
494 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
495 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
496 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
498 //0020 0011 "IS" Series Number
499 sprintf(charImagePosition,"%04d",serieNumber);
500 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
502 //0020 0011 "IS" Instance Number
503 sprintf(charImagePosition,"%04d",imageNumber);
504 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
506 //0018 1060 "DS" Time Trigger
507 sprintf(charImagePosition,"%f",currentTime);
508 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
512 GDCM_NAME_SPACE::FileHelper *fh;
513 fh = GDCM_NAME_SPACE::FileHelper::New(file);
514 // cast is just to avoid warnings (*no* conversion)
515 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
516 fh->SetWriteModeToRaw();
517 fh->SetWriteTypeToDcmExplVR();
519 fh->SetWriteTypeToDcmExplVR();
522 sprintf(numero, "%02d", l1);
523 std::string fileName = imageName + "." + numero + ".dcm";
525 std::cout << "fileName " << fileName << std::endl;
527 if( !fh->Write(fileName))
528 std::cout << "Failed for [" << fileName << "]\n"
529 << " File is unwrittable" << std::endl;
532 currentTime += temporalResolution;
535 } // end loop on frames
538 // Anatomical Images.
539 std::cout << " ========= Create Anatomical images" << std::endl;
541 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
543 currentTime = timeStart;
545 for (int fr=0; fr<nf; fr++) // Loop on the frames
547 //std::cout << "Frame nb " << l << std::endl;
548 for( int j=0; j<ny; j++)
551 for (int i=0; i<nx; i++)
566 val = (float)atof(str1.c_str());
567 // std::cout << " " << val;
568 *(f+ /*l*nx*ny + */j*nx+i) = val;
572 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
580 uint16_t *img = new uint16_t[ny*nx];
582 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
583 for( int k=0; k<ny*nx; k++)
585 *ptr = (int16_t)(*tmp);
589 std::ostringstream str;
590 GDCM_NAME_SPACE::File *file;
591 file = GDCM_NAME_SPACE::File::New();
593 // Set the image size
596 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
599 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
601 // Set the pixel type
603 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
605 str << 16; // may be 12 or 16 if componentSize =16
606 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
608 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
610 // Set the pixel representation // 0/1 1 : signed
611 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
613 // Set the samples per pixel // 1:Grey level, 3:RGB
614 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
617 // Set Rescale Intercept
620 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
625 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
628 // 0020 0037 DS 6 Image Orientation (Patient)
629 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
631 // 0020 0032 DS 3 Image Position (Patient)
632 char charImagePosition[256];
633 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
635 // 0020 0x1041 DS 1 Slice Location
636 sprintf(charImagePosition,"%f",float(l1%nf));
637 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
639 //0008 103e LO 1 Series Description
640 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
642 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
643 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
644 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
646 //0020 0011 "IS" Series Number
647 sprintf(charImagePosition,"%04d",serieNumber+1);
648 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
650 //0020 0011 "IS" Instance Number
651 sprintf(charImagePosition,"%04d",imageNumber);
652 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
654 //0018 1060 "DS" Time Trigger
655 sprintf(charImagePosition,"%f",currentTime);
656 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
659 GDCM_NAME_SPACE::FileHelper *fh;
660 fh = GDCM_NAME_SPACE::FileHelper::New(file);
661 // cast is just to avoid warnings (*no* conversion)
662 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
663 fh->SetWriteModeToRaw();
664 fh->SetWriteTypeToDcmExplVR();
666 fh->SetWriteTypeToDcmExplVR();
669 sprintf(numero, "%02d", l1);
670 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
671 std::cout << "fileName " << fileName << std::endl;
673 if( !fh->Write(fileName))
674 std::cout << "Failed for [" << fileName << "]\n"
675 << " File is unwrittable" << std::endl;
678 currentTime += temporalResolution;
681 } // end loop on frames