1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2007/09/18 11:01:55 $
7 Version: $Revision: 1.8 $
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;
184 std::string strStudyUID;
185 std::string strSerieUID;
187 if (userDefinedStudy)
188 strSerieUID = studyUID;
190 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
192 if (userDefinedStudy)
193 strSerieUID = serieUID;
195 strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
199 GDCM_NAME_SPACE::DirListType fileNames;
200 fileNames = dirList.GetFilenames();
201 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
202 it != fileNames.end();
205 if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
210 int sz = (*it).size();
211 if ( (*it).c_str()[sz-1] != 't')
213 // skip non .txt files
216 std::ifstream from( (*it).c_str() );
219 std::cout << "Can't open file" << *it << std::endl;
225 std::cout << "Success in open file" << *it << std::endl;
226 Load(from, *it, patName, strStudyUID, strSerieUID, serieNumber, verbose);
234 void Load(std::ifstream &from, std::string imageName, std::string patName,
235 std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose)
238 std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
242 std::cout << " ========= Create Parametric images" << std::endl;
243 /* was OK for single frames
245 ---------------------------
246 Array dimensions = 58 x 56
247 The following is the 2D array of peak circumferential strain map, all zero value
248 pixels are outside the mask
249 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
250 -----------------------------
261 std::cout << "nx " << nx << " ny " << ny << std::endl;
264 std::getline(from, str1);
265 std::cout << "[" << str1 << "]" << std::endl;
266 std::getline(from, str1);
267 std::cout << "[" << str1 << "]" << std::endl;
268 std::getline(from, str1);
269 std::cout << "[" << str1 << "]" << std::endl;
272 /* now, we deal with multiframes
274 ------------------------------------------
275 X dim, Y dim, N of frames = 52x59x14
276 Temporal resolution = 31.9200 ms
277 First frame starts at 47.9600 ms
278 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
279 All pixels with zero strain values are outside the masks.
280 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
281 --------------------------------------------
285 from >> str1; // X dim,
287 from >> str1; // Y dim,
289 from >> str1; // N of frames =
293 from >> str1; // 52x59x14
296 std::cout << "[" << str1 << "]" << std::endl;
298 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
299 std::cout << nx << " " << ny << " " << nf << std::endl;
301 std::getline(from, str1);
303 from >> str1; // Temporal
304 from >> str1; // Resolution
309 float temporalResolution;
310 sscanf( str1.c_str(),"%f",&temporalResolution);
312 std::cout << "temporal Resolution = " << temporalResolution << std::endl;
313 std::getline(from, str1);
315 from >> str1; // First
316 from >> str1; // frame
317 from >> str1; // starts
322 sscanf( str1.c_str(),"%f",&timeStart);
323 std::cout << "time Start = " << timeStart << std::endl;
324 std::getline(from, str1);
327 for (int k=0; k<2; k++) //
329 std::getline(from, str1);
330 std::cout << str1 << std::endl;
333 //float *f = new float(nx*ny);
334 // -->float *f = new float[nx*ny]; // Would be better !
335 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
336 // float mini = FLT_MAX, maxi = FLT_MIN;
341 currentTime = timeStart;
343 for (l1=0; l1<nf; l1++) // Loop on the frames
345 //std::cout << "Frame nb " << l1 << std::endl;
346 for( int j=0; j<ny; j++)
348 for (int i=0; i<nx; i++)
353 for (;;) //eatwhite(from);
357 std::cout << " !from.get(c) ";
362 //std::cout << " !isspace(c) ";
366 } //end eatwhite(from);
368 // trouble : when space is missing "-0.0990263-8.8778"
369 // is not interpreted as TWO values :-(
379 if ( first != 1 && c == '-' && previous != 'e')
382 //std::cout << " One more gauffre in frame:" << std::dec << l
383 // << ", line : " << j << " element " << i << std::endl;
392 val = (float)atof(str1.c_str());
393 //std::cout << " " << val;
394 *(f+ /*l*nx*ny + */j*nx+i) = val;
398 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
404 //std::cout << std::endl;
405 //std::cout << std::endl << " line nb : " << j
406 // << " line length : " << l << std::endl;
409 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
411 // values are expressed as % as a fraction, actually!)
412 // It's useless to rescale them as uint16_t : just *100
413 uint16_t *img = new uint16_t[ny*nx];
416 float div = maxi-mini;
417 std::cout << "div : " << div << " mini : " << mini << std::endl;
418 for( int k=0; k<ny*nx; k++)
420 *ptr = (uint16_t)((*tmp * 65535.0) / div);
425 int16_t *img = new int16_t[ny*nx];
427 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
428 for( int k=0; k<ny*nx; k++)
430 if(*tmp > 1.0) // artefacted pixel
432 else /// \todo FIXME : what about max threshold ?
433 *ptr = (int16_t)(*tmp *100);
435 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
440 // GDCM_NAME_SPACE::Debug::DebugOn();
442 std::ostringstream str;
443 GDCM_NAME_SPACE::File *file;
444 file = GDCM_NAME_SPACE::File::New();
446 // Set the image size
449 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
452 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
454 // Set the pixel type
456 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
458 str << 16; // may be 12 or 16 if componentSize =16
459 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
461 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
463 // Set the pixel representation // 0/1 1 : signed
464 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
466 // Set the samples per pixel // 1:Grey level, 3:RGB
467 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
470 // Set Rescale Intercept
473 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
478 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
481 // 0020 0037 DS 6 Image Orientation (Patient)
482 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
484 // 0020 0032 DS 3 Image Position (Patient)
485 char charImagePosition[256];
486 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
488 // 0020 0x1041 DS 1 Slice Location
489 sprintf(charImagePosition,"%f",float(l1%nf));
490 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
492 //0008 103e LO 1 Series Description
493 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
495 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
496 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
497 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
499 //0020 0011 "IS" Series Number
500 sprintf(charImagePosition,"%04d",serieNumber);
501 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
503 //0020 0011 "IS" Instance Number
504 sprintf(charImagePosition,"%04d",imageNumber);
505 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
507 //0018 1060 "DS" Time Trigger
508 sprintf(charImagePosition,"%f",currentTime);
509 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
513 GDCM_NAME_SPACE::FileHelper *fh;
514 fh = GDCM_NAME_SPACE::FileHelper::New(file);
515 // cast is just to avoid warnings (*no* conversion)
516 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
517 fh->SetWriteModeToRaw();
518 fh->SetWriteTypeToDcmExplVR();
520 fh->SetWriteTypeToDcmExplVR();
523 sprintf(numero, "%02d", l1);
524 std::string fileName = imageName + "." + numero + ".dcm";
526 std::cout << "fileName " << fileName << std::endl;
528 if( !fh->Write(fileName))
529 std::cout << "Failed for [" << fileName << "]\n"
530 << " File is unwrittable" << std::endl;
533 currentTime += temporalResolution;
536 } // end loop on frames
539 // Anatomical Images.
540 std::cout << " ========= Create Anatomical images" << std::endl;
542 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
544 currentTime = timeStart;
546 for (int fr=0; fr<nf; fr++) // Loop on the frames
548 //std::cout << "Frame nb " << l << std::endl;
549 for( int j=0; j<ny; j++)
552 for (int i=0; i<nx; i++)
567 val = (float)atof(str1.c_str());
568 // std::cout << " " << val;
569 *(f+ /*l*nx*ny + */j*nx+i) = val;
573 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
581 uint16_t *img = new uint16_t[ny*nx];
583 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
584 for( int k=0; k<ny*nx; k++)
586 *ptr = (int16_t)(*tmp);
590 std::ostringstream str;
591 GDCM_NAME_SPACE::File *file;
592 file = GDCM_NAME_SPACE::File::New();
594 // Set the image size
597 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
600 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
602 // Set the pixel type
604 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
606 str << 16; // may be 12 or 16 if componentSize =16
607 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
609 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
611 // Set the pixel representation // 0/1 1 : signed
612 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
614 // Set the samples per pixel // 1:Grey level, 3:RGB
615 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
618 // Set Rescale Intercept
621 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
626 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
629 // 0020 0037 DS 6 Image Orientation (Patient)
630 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
632 // 0020 0032 DS 3 Image Position (Patient)
633 char charImagePosition[256];
634 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
636 // 0020 0x1041 DS 1 Slice Location
637 sprintf(charImagePosition,"%f",float(l1%nf));
638 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
640 //0008 103e LO 1 Series Description
641 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
643 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
644 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
645 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
647 //0020 0011 "IS" Series Number
648 sprintf(charImagePosition,"%04d",serieNumber+1);
649 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
651 //0020 0011 "IS" Instance Number
652 sprintf(charImagePosition,"%04d",imageNumber);
653 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
655 //0018 1060 "DS" Time Trigger
656 sprintf(charImagePosition,"%f",currentTime);
657 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
660 GDCM_NAME_SPACE::FileHelper *fh;
661 fh = GDCM_NAME_SPACE::FileHelper::New(file);
662 // cast is just to avoid warnings (*no* conversion)
663 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
664 fh->SetWriteModeToRaw();
665 fh->SetWriteTypeToDcmExplVR();
667 fh->SetWriteTypeToDcmExplVR();
670 sprintf(numero, "%02d", l1);
671 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
672 std::cout << "fileName " << fileName << std::endl;
674 if( !fh->Write(fileName))
675 std::cout << "Failed for [" << fileName << "]\n"
676 << " File is unwrittable" << std::endl;
679 currentTime += temporalResolution;
682 } // end loop on frames