1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2007/10/29 17:13:59 $
7 Version: $Revision: 1.10 $
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 =========================================================================*/
23 #if defined(__BORLANDC__)
28 #include "gdcmFileHelper.h"
29 #include "gdcmDebug.h"
30 #include "gdcmDirList.h"
33 #include "gdcmArgMgr.h"
37 * - explores recursively the (single Patient, single Study) directory
38 * - examines the ".txt" Dense multiframe files
39 * - Converts the files into 16 bits Dicom Files
40 * WARNING : directory must contain ONLY .txt files
43 void Load(std::ifstream &from, std::string imageName, std::string patName,
44 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose);
46 //std::ifstream& eatwhite(std::ifstream& is);
49 window center (level) and width, are defined in the DICOM
50 standard very precisely as follows (see PS 3.3 C.11.2.1.2):
52 >"These Attributes are applied according to the following pseudo-code,
55 y is an output value with a range from ymin to ymax,
56 c is Window Center (0028,1050)
57 w is Window Width (0028,1051):
59 > if (x <= c - 0.5 - (w-1)/2), then y = ymin
60 > else if (x > c - 0.5 + (w-1)/2), then y = ymax,
61 > else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin
65 From: David Clunie - view profile
66 Date: Thurs, Jun 1 2006 3:03 pm
67 Email: David Clunie <dclu...@dclunie.com>
68 Groups: comp.protocols.dicom
70 The value of x is the output of the preceding step, the so-called
71 "modality LUT", which may either be:
73 - identity (no rescale values or Modality LUT, or the SOP Class is
74 PET and rescale values are ignored), in which case x is the stored
75 pixel value in (7FE0,0010)
77 - Rescale Slope and Intercept (as typically used in CT), in which
78 case x is the value obtained from applying the rescale values to
79 the stored pixel value
81 - an actual LUT, in which case x is the value stored in the LUT
82 corresponding to the LUT index value that is the stored pixel
85 The ymin and ymax are intended to represent the output range; for
86 example, if the hypothetical Presentation LUT step that follows
87 the VOI LUT (window) stage is an identity operation, then the
88 ymin and ymax represent P-Values, which might be the range of
89 digital driving levels for your display (calibrated to the GSDF),
90 in the 8-bit wide output case ranging from 0 to 255, for example.
92 The terms brightness and contrast are not used in radiology imaging
93 -the window center and width are used instead.
96 int main(int argc, char *argv[])
99 " \n DenseMultiFramesToDicom : \n",
100 " - explores recursively the given (single Patient, single Study) directory",
101 " - examines the '.txt' files ",
102 " - Converts the files into 16 bits Dicom files, ",
103 " WARNING : directory must contain ONLY .txt files ",
105 " DenseMultiFramesToDicom dirin=rootDirectoryName ",
106 " [studyUID = ] [patName = ] ",
107 " [listonly] [verbose] [debug] ",
109 " studyUID : *aware* user wants to add the serie ",
110 " to an already existing study ",
111 " verbose : user wants to run the program in 'verbose mode' ",
112 " debug : *developer* wants to run the program in 'debug mode' ",
115 // ----- Initialize Arguments Manager ------
117 GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
119 if (argc == 1 || am->ArgMgrDefined("usage"))
121 am->ArgMgrUsage(usage); // Display 'usage'
126 const char *dirNamein;
127 dirNamein = am->ArgMgrGetString("dirin",".");
129 if (am->ArgMgrDefined("debug"))
130 GDCM_NAME_SPACE::Debug::DebugOn();
132 int verbose = am->ArgMgrDefined("verbose");
133 int listonly = am->ArgMgrDefined("listonly");
134 std::string patName = am->ArgMgrGetString("patname", dirNamein);
136 bool userDefinedStudy = ( 0 != am->ArgMgrDefined("studyUID") );
138 const char *studyUID;
139 if (userDefinedStudy)
140 studyUID = am->ArgMgrGetString("studyUID");
142 // not described *on purpose* in the Usage !
143 bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );
145 const char *serieUID;
147 serieUID = am->ArgMgrGetString("serieUID");
149 // if unused Param we give up
150 if ( am->ArgMgrPrintUnusedLabels() )
152 am->ArgMgrUsage(usage);
156 delete am; // we don't need Argument Manager any longer
158 // ----- Begin Processing -----
160 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
162 std::cout << "KO : [" << dirNamein << "] is not a Directory."
169 std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
172 std::string strDirNamein(dirNamein);
173 GDCM_NAME_SPACE::DirList dirList(strDirNamein, true); // (recursively) the list of files
177 std::cout << "------------List of found files ------------" << std::endl;
179 std::cout << std::endl;
183 std::string filenameout;
185 std::string strStudyUID;
186 std::string strSerieUID;
188 if (userDefinedStudy)
189 strSerieUID = studyUID;
191 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
193 if (userDefinedStudy)
194 strSerieUID = serieUID;
196 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
200 GDCM_NAME_SPACE::DirListType fileNames;
201 fileNames = dirList.GetFilenames();
202 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
203 it != fileNames.end();
206 if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
211 int sz = (*it).size();
212 if ( (*it).c_str()[sz-1] != 't')
214 // skip non .txt files
217 std::ifstream from( (*it).c_str() );
220 std::cout << "Can't open file" << *it << std::endl;
226 std::cout << "Success in open file" << *it << std::endl;
227 Load(from, *it, patName, strStudyUID, strSerieUID, serieNumber, verbose);
236 void Load(std::ifstream &from, std::string imageName, std::string patName,
237 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose)
240 std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
244 std::cout << " ========= Create Parametric images" << std::endl;
245 /* was OK for single frames
247 ---------------------------
248 Array dimensions = 58 x 56
249 The following is the 2D array of peak circumferential strain map, all zero value
250 pixels are outside the mask
251 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
252 -----------------------------
263 std::cout << "nx " << nx << " ny " << ny << std::endl;
266 std::getline(from, str1);
267 std::cout << "[" << str1 << "]" << std::endl;
268 std::getline(from, str1);
269 std::cout << "[" << str1 << "]" << std::endl;
270 std::getline(from, str1);
271 std::cout << "[" << str1 << "]" << std::endl;
274 /* now, we deal with multiframes
276 ------------------------------------------
277 X dim, Y dim, N of frames = 52x59x14
278 Temporal resolution = 31.9200 ms
279 First frame starts at 47.9600 ms
280 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
281 All pixels with zero strain values are outside the masks.
282 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
283 --------------------------------------------
287 from >> str1; // X dim,
289 from >> str1; // Y dim,
291 from >> str1; // N of frames =
295 from >> str1; // 52x59x14
298 std::cout << "[" << str1 << "]" << std::endl;
300 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
301 std::cout << nx << " " << ny << " " << nf << std::endl;
303 std::getline(from, str1);
305 from >> str1; // Temporal
306 from >> str1; // Resolution
311 float temporalResolution;
312 sscanf( str1.c_str(),"%f",&temporalResolution);
314 std::cout << "temporal Resolution = " << temporalResolution << std::endl;
315 std::getline(from, str1);
317 from >> str1; // First
318 from >> str1; // frame
319 from >> str1; // starts
324 sscanf( str1.c_str(),"%f",&timeStart);
325 std::cout << "time Start = " << timeStart << std::endl;
326 std::getline(from, str1);
329 for (int k=0; k<2; k++) //
331 std::getline(from, str1);
332 std::cout << str1 << std::endl;
335 //float *f = new float(nx*ny);
336 // -->float *f = new float[nx*ny]; // Would be better !
337 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
338 // float mini = FLT_MAX, maxi = FLT_MIN;
343 currentTime = timeStart;
345 for (l1=0; l1<nf; l1++) // Loop on the frames
347 //std::cout << "Frame nb " << l1 << std::endl;
348 for( int j=0; j<ny; j++)
350 for (int i=0; i<nx; i++)
355 for (;;) //eatwhite(from);
359 std::cout << " !from.get(c) ";
364 //std::cout << " !isspace(c) ";
368 } //end eatwhite(from);
370 // trouble : when space is missing "-0.0990263-8.8778"
371 // is not interpreted as TWO values :-(
381 if ( first != 1 && c == '-' && previous != 'e')
384 //std::cout << " One more gauffre in frame:" << std::dec << l
385 // << ", line : " << j << " element " << i << std::endl;
394 val = (float)atof(str1.c_str());
395 //std::cout << " " << val;
396 *(f+ /*l*nx*ny + */j*nx+i) = val;
400 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
406 //std::cout << std::endl;
407 //std::cout << std::endl << " line nb : " << j
408 // << " line length : " << l << std::endl;
411 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
413 // values are expressed as % as a fraction, actually!)
414 // It's useless to rescale them as uint16_t : just *100
415 uint16_t *img = new uint16_t[ny*nx];
418 float div = maxi-mini;
419 std::cout << "div : " << div << " mini : " << mini << std::endl;
420 for( int k=0; k<ny*nx; k++)
422 *ptr = (uint16_t)((*tmp * 65535.0) / div);
427 int16_t *img = new int16_t[ny*nx];
429 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
430 for( int k=0; k<ny*nx; k++)
432 if(*tmp > 1.0) // artefacted pixel
434 else /// \todo FIXME : what about max threshold ?
435 *ptr = (int16_t)(*tmp *100);
437 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
442 // GDCM_NAME_SPACE::Debug::DebugOn();
444 std::ostringstream str;
445 GDCM_NAME_SPACE::File *file;
446 file = GDCM_NAME_SPACE::File::New();
448 // Set the image size
451 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
454 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
456 // Set the pixel type
458 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
460 str << 16; // may be 12 or 16 if componentSize =16
461 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
463 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
465 // Set the pixel representation // 0/1 1 : signed
466 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
468 // Set the samples per pixel // 1:Grey level, 3:RGB
469 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
472 // Set Rescale Intercept
475 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
480 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
483 // 0020 0037 DS 6 Image Orientation (Patient)
484 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
486 // 0020 0032 DS 3 Image Position (Patient)
487 char charImagePosition[256];
488 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
490 // 0020 0x1041 DS 1 Slice Location
491 sprintf(charImagePosition,"%f",float(l1%nf));
492 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
494 //0008 103e LO 1 Series Description
495 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
497 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
498 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
499 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
501 //0020 0011 "IS" Series Number
502 sprintf(charImagePosition,"%04d",serieNumber);
503 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
505 //0020 0011 "IS" Instance Number
506 sprintf(charImagePosition,"%04d",imageNumber);
507 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
509 //0018 1060 "DS" Time Trigger
510 sprintf(charImagePosition,"%f",currentTime);
511 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
515 GDCM_NAME_SPACE::FileHelper *fh;
516 fh = GDCM_NAME_SPACE::FileHelper::New(file);
517 // cast is just to avoid warnings (*no* conversion)
518 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
519 fh->SetWriteModeToRaw();
520 fh->SetWriteTypeToDcmExplVR();
522 fh->SetWriteTypeToDcmExplVR();
525 sprintf(numero, "%02d", l1);
526 std::string fileName = imageName + "." + numero + ".dcm";
528 std::cout << "fileName " << fileName << std::endl;
530 if( !fh->Write(fileName))
531 std::cout << "Failed for [" << fileName << "]\n"
532 << " File is unwrittable" << std::endl;
535 currentTime += temporalResolution;
538 } // end loop on frames
541 // Anatomical Images.
542 std::cout << " ========= Create Anatomical images" << std::endl;
544 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
546 currentTime = timeStart;
548 for (int fr=0; fr<nf; fr++) // Loop on the frames
550 //std::cout << "Frame nb " << l << std::endl;
551 for( int j=0; j<ny; j++)
554 for (int i=0; i<nx; i++)
569 val = (float)atof(str1.c_str());
570 // std::cout << " " << val;
571 *(f+ /*l*nx*ny + */j*nx+i) = val;
575 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
583 uint16_t *img = new uint16_t[ny*nx];
585 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
586 for( int k=0; k<ny*nx; k++)
588 *ptr = (int16_t)(*tmp);
592 std::ostringstream str;
593 GDCM_NAME_SPACE::File *file;
594 file = GDCM_NAME_SPACE::File::New();
596 // Set the image size
599 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
602 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
604 // Set the pixel type
606 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
608 str << 16; // may be 12 or 16 if componentSize =16
609 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
611 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
613 // Set the pixel representation // 0/1 1 : signed
614 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
616 // Set the samples per pixel // 1:Grey level, 3:RGB
617 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
620 // Set Rescale Intercept
623 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
628 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
631 // 0020 0037 DS 6 Image Orientation (Patient)
632 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
634 // 0020 0032 DS 3 Image Position (Patient)
635 char charImagePosition[256];
636 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
638 // 0020 0x1041 DS 1 Slice Location
639 sprintf(charImagePosition,"%f",float(l1%nf));
640 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
642 //0008 103e LO 1 Series Description
643 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
645 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
646 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
647 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
649 //0020 0011 "IS" Series Number
650 sprintf(charImagePosition,"%04d",serieNumber+1);
651 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
653 //0020 0011 "IS" Instance Number
654 sprintf(charImagePosition,"%04d",imageNumber);
655 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
657 //0018 1060 "DS" Time Trigger
658 sprintf(charImagePosition,"%f",currentTime);
659 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
662 GDCM_NAME_SPACE::FileHelper *fh;
663 fh = GDCM_NAME_SPACE::FileHelper::New(file);
664 // cast is just to avoid warnings (*no* conversion)
665 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
666 fh->SetWriteModeToRaw();
667 fh->SetWriteTypeToDcmExplVR();
669 fh->SetWriteTypeToDcmExplVR();
672 sprintf(numero, "%02d", l1);
673 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
674 std::cout << "fileName " << fileName << std::endl;
676 if( !fh->Write(fileName))
677 std::cout << "Failed for [" << fileName << "]\n"
678 << " File is unwrittable" << std::endl;
681 currentTime += temporalResolution;
684 } // end loop on frames