1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2007/09/20 11:15:54 $
7 Version: $Revision: 1.9 $
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 = ( 0 != am->ArgMgrDefined("studyUID") );
134 const char *studyUID;
135 if (userDefinedStudy)
136 studyUID = am->ArgMgrGetString("studyUID");
138 // not described *on purpose* in the Usage !
139 bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );
141 const char *serieUID;
143 serieUID = am->ArgMgrGetString("serieUID");
145 // if unused Param we give up
146 if ( am->ArgMgrPrintUnusedLabels() )
148 am->ArgMgrUsage(usage);
152 delete am; // we don't need Argument Manager any longer
154 // ----- Begin Processing -----
156 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
158 std::cout << "KO : [" << dirNamein << "] is not a Directory."
165 std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
168 std::string strDirNamein(dirNamein);
169 GDCM_NAME_SPACE::DirList dirList(strDirNamein, true); // (recursively) the list of files
173 std::cout << "------------List of found files ------------" << std::endl;
175 std::cout << std::endl;
179 std::string filenameout;
181 std::string strStudyUID;
182 std::string strSerieUID;
184 if (userDefinedStudy)
185 strSerieUID = studyUID;
187 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
189 if (userDefinedStudy)
190 strSerieUID = serieUID;
192 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
196 GDCM_NAME_SPACE::DirListType fileNames;
197 fileNames = dirList.GetFilenames();
198 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
199 it != fileNames.end();
202 if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
207 int sz = (*it).size();
208 if ( (*it).c_str()[sz-1] != 't')
210 // skip non .txt files
213 std::ifstream from( (*it).c_str() );
216 std::cout << "Can't open file" << *it << std::endl;
222 std::cout << "Success in open file" << *it << std::endl;
223 Load(from, *it, patName, strStudyUID, strSerieUID, serieNumber, verbose);
232 void Load(std::ifstream &from, std::string imageName, std::string patName,
233 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose)
236 std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
240 std::cout << " ========= Create Parametric images" << std::endl;
241 /* was OK for single frames
243 ---------------------------
244 Array dimensions = 58 x 56
245 The following is the 2D array of peak circumferential strain map, all zero value
246 pixels are outside the mask
247 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
248 -----------------------------
259 std::cout << "nx " << nx << " ny " << ny << std::endl;
262 std::getline(from, str1);
263 std::cout << "[" << str1 << "]" << std::endl;
264 std::getline(from, str1);
265 std::cout << "[" << str1 << "]" << std::endl;
266 std::getline(from, str1);
267 std::cout << "[" << str1 << "]" << std::endl;
270 /* now, we deal with multiframes
272 ------------------------------------------
273 X dim, Y dim, N of frames = 52x59x14
274 Temporal resolution = 31.9200 ms
275 First frame starts at 47.9600 ms
276 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
277 All pixels with zero strain values are outside the masks.
278 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
279 --------------------------------------------
283 from >> str1; // X dim,
285 from >> str1; // Y dim,
287 from >> str1; // N of frames =
291 from >> str1; // 52x59x14
294 std::cout << "[" << str1 << "]" << std::endl;
296 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
297 std::cout << nx << " " << ny << " " << nf << std::endl;
299 std::getline(from, str1);
301 from >> str1; // Temporal
302 from >> str1; // Resolution
307 float temporalResolution;
308 sscanf( str1.c_str(),"%f",&temporalResolution);
310 std::cout << "temporal Resolution = " << temporalResolution << std::endl;
311 std::getline(from, str1);
313 from >> str1; // First
314 from >> str1; // frame
315 from >> str1; // starts
320 sscanf( str1.c_str(),"%f",&timeStart);
321 std::cout << "time Start = " << timeStart << std::endl;
322 std::getline(from, str1);
325 for (int k=0; k<2; k++) //
327 std::getline(from, str1);
328 std::cout << str1 << std::endl;
331 //float *f = new float(nx*ny);
332 // -->float *f = new float[nx*ny]; // Would be better !
333 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
334 // float mini = FLT_MAX, maxi = FLT_MIN;
339 currentTime = timeStart;
341 for (l1=0; l1<nf; l1++) // Loop on the frames
343 //std::cout << "Frame nb " << l1 << std::endl;
344 for( int j=0; j<ny; j++)
346 for (int i=0; i<nx; i++)
351 for (;;) //eatwhite(from);
355 std::cout << " !from.get(c) ";
360 //std::cout << " !isspace(c) ";
364 } //end eatwhite(from);
366 // trouble : when space is missing "-0.0990263-8.8778"
367 // is not interpreted as TWO values :-(
377 if ( first != 1 && c == '-' && previous != 'e')
380 //std::cout << " One more gauffre in frame:" << std::dec << l
381 // << ", line : " << j << " element " << i << std::endl;
390 val = (float)atof(str1.c_str());
391 //std::cout << " " << val;
392 *(f+ /*l*nx*ny + */j*nx+i) = val;
396 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
402 //std::cout << std::endl;
403 //std::cout << std::endl << " line nb : " << j
404 // << " line length : " << l << std::endl;
407 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
409 // values are expressed as % as a fraction, actually!)
410 // It's useless to rescale them as uint16_t : just *100
411 uint16_t *img = new uint16_t[ny*nx];
414 float div = maxi-mini;
415 std::cout << "div : " << div << " mini : " << mini << std::endl;
416 for( int k=0; k<ny*nx; k++)
418 *ptr = (uint16_t)((*tmp * 65535.0) / div);
423 int16_t *img = new int16_t[ny*nx];
425 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
426 for( int k=0; k<ny*nx; k++)
428 if(*tmp > 1.0) // artefacted pixel
430 else /// \todo FIXME : what about max threshold ?
431 *ptr = (int16_t)(*tmp *100);
433 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
438 // GDCM_NAME_SPACE::Debug::DebugOn();
440 std::ostringstream str;
441 GDCM_NAME_SPACE::File *file;
442 file = GDCM_NAME_SPACE::File::New();
444 // Set the image size
447 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
450 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
452 // Set the pixel type
454 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
456 str << 16; // may be 12 or 16 if componentSize =16
457 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
459 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
461 // Set the pixel representation // 0/1 1 : signed
462 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
464 // Set the samples per pixel // 1:Grey level, 3:RGB
465 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
468 // Set Rescale Intercept
471 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
476 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
479 // 0020 0037 DS 6 Image Orientation (Patient)
480 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
482 // 0020 0032 DS 3 Image Position (Patient)
483 char charImagePosition[256];
484 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
486 // 0020 0x1041 DS 1 Slice Location
487 sprintf(charImagePosition,"%f",float(l1%nf));
488 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
490 //0008 103e LO 1 Series Description
491 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
493 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
494 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
495 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
497 //0020 0011 "IS" Series Number
498 sprintf(charImagePosition,"%04d",serieNumber);
499 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
501 //0020 0011 "IS" Instance Number
502 sprintf(charImagePosition,"%04d",imageNumber);
503 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
505 //0018 1060 "DS" Time Trigger
506 sprintf(charImagePosition,"%f",currentTime);
507 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
511 GDCM_NAME_SPACE::FileHelper *fh;
512 fh = GDCM_NAME_SPACE::FileHelper::New(file);
513 // cast is just to avoid warnings (*no* conversion)
514 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
515 fh->SetWriteModeToRaw();
516 fh->SetWriteTypeToDcmExplVR();
518 fh->SetWriteTypeToDcmExplVR();
521 sprintf(numero, "%02d", l1);
522 std::string fileName = imageName + "." + numero + ".dcm";
524 std::cout << "fileName " << fileName << std::endl;
526 if( !fh->Write(fileName))
527 std::cout << "Failed for [" << fileName << "]\n"
528 << " File is unwrittable" << std::endl;
531 currentTime += temporalResolution;
534 } // end loop on frames
537 // Anatomical Images.
538 std::cout << " ========= Create Anatomical images" << std::endl;
540 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
542 currentTime = timeStart;
544 for (int fr=0; fr<nf; fr++) // Loop on the frames
546 //std::cout << "Frame nb " << l << std::endl;
547 for( int j=0; j<ny; j++)
550 for (int i=0; i<nx; i++)
565 val = (float)atof(str1.c_str());
566 // std::cout << " " << val;
567 *(f+ /*l*nx*ny + */j*nx+i) = val;
571 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
579 uint16_t *img = new uint16_t[ny*nx];
581 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
582 for( int k=0; k<ny*nx; k++)
584 *ptr = (int16_t)(*tmp);
588 std::ostringstream str;
589 GDCM_NAME_SPACE::File *file;
590 file = GDCM_NAME_SPACE::File::New();
592 // Set the image size
595 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
598 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
600 // Set the pixel type
602 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
604 str << 16; // may be 12 or 16 if componentSize =16
605 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
607 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
609 // Set the pixel representation // 0/1 1 : signed
610 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
612 // Set the samples per pixel // 1:Grey level, 3:RGB
613 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
616 // Set Rescale Intercept
619 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
624 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
627 // 0020 0037 DS 6 Image Orientation (Patient)
628 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
630 // 0020 0032 DS 3 Image Position (Patient)
631 char charImagePosition[256];
632 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
634 // 0020 0x1041 DS 1 Slice Location
635 sprintf(charImagePosition,"%f",float(l1%nf));
636 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
638 //0008 103e LO 1 Series Description
639 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
641 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
642 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
643 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
645 //0020 0011 "IS" Series Number
646 sprintf(charImagePosition,"%04d",serieNumber+1);
647 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
649 //0020 0011 "IS" Instance Number
650 sprintf(charImagePosition,"%04d",imageNumber);
651 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
653 //0018 1060 "DS" Time Trigger
654 sprintf(charImagePosition,"%f",currentTime);
655 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
658 GDCM_NAME_SPACE::FileHelper *fh;
659 fh = GDCM_NAME_SPACE::FileHelper::New(file);
660 // cast is just to avoid warnings (*no* conversion)
661 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
662 fh->SetWriteModeToRaw();
663 fh->SetWriteTypeToDcmExplVR();
665 fh->SetWriteTypeToDcmExplVR();
668 sprintf(numero, "%02d", l1);
669 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
670 std::cout << "fileName " << fileName << std::endl;
672 if( !fh->Write(fileName))
673 std::cout << "Failed for [" << fileName << "]\n"
674 << " File is unwrittable" << std::endl;
677 currentTime += temporalResolution;
680 } // end loop on frames