1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2007/06/26 15:42:14 $
7 Version: $Revision: 1.6 $
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 DenseMultiFramesToDicom : \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 " DenseMultiFramesToDicom 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_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::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_NAME_SPACE::Debug::DebugOn();
129 int verbose = am->ArgMgrDefined("verbose");
130 int listonly = am->ArgMgrDefined("listonly");
131 std::string patName = am->ArgMgrGetString("patname", dirNamein);
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_NAME_SPACE::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_NAME_SPACE::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_NAME_SPACE::Util::CreateUniqueUID();
179 GDCM_NAME_SPACE::DirListType fileNames;
180 fileNames = dirList.GetFilenames();
181 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
182 it != fileNames.end();
185 if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
190 int sz = (*it).size();
191 if ( (*it).c_str()[sz-1] != 't')
193 // skip non .txt files
196 std::ifstream from( (*it).c_str() );
199 std::cout << "Can't open file" << *it << std::endl;
204 std::cout << "Success in open file" << *it << std::endl;
205 Load(from, *it, patName, strStudyUID, serieNumber);
213 void Load(std::ifstream &from, std::string imageName, std::string patName,
214 std::string strStudyUID, int serieNumber)
216 std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
219 std::cout << " ========= Create Parametric images" << std::endl;
220 /* was OK for single frames
222 ---------------------------
223 Array dimensions = 58 x 56
224 The following is the 2D array of peak circumferential strain map, all zero value
225 pixels are outside the mask
226 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
227 -----------------------------
238 std::cout << "nx " << nx << " ny " << ny << std::endl;
241 std::getline(from, str1);
242 std::cout << "[" << str1 << "]" << std::endl;
243 std::getline(from, str1);
244 std::cout << "[" << str1 << "]" << std::endl;
245 std::getline(from, str1);
246 std::cout << "[" << str1 << "]" << std::endl;
249 /* now, we deal with multiframes
251 ------------------------------------------
252 X dim, Y dim, N of frames = 52x59x14
253 Temporal resolution = 31.9200 ms
254 First frame starts at 47.9600 ms
255 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
256 All pixels with zero strain values are outside the masks.
257 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
258 --------------------------------------------
262 from >> str1; // X dim,
264 from >> str1; // Y dim,
266 from >> str1; // N of frames =
270 from >> str1; // 52x59x14
272 std::cout << "[" << str1 << "]" << std::endl;
274 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
275 std::cout << nx << " " << ny << " " << nf << std::endl;
277 std::getline(from, str1);
279 from >> str1; // Temporal
280 from >> str1; // Resolution
285 float temporalResolution;
286 sscanf( str1.c_str(),"%f",&temporalResolution);
287 std::cout << "temporal Resolution = " << temporalResolution << std::endl;
288 std::getline(from, str1);
290 from >> str1; // First
291 from >> str1; // frame
292 from >> str1; // starts
297 sscanf( str1.c_str(),"%f",&timeStart);
298 std::cout << "time Start = " << timeStart << std::endl;
299 std::getline(from, str1);
302 for (int k=0; k<2; k++) //
304 std::getline(from, str1);
305 std::cout << str1 << std::endl;
308 //float *f = new float(nx*ny);
309 // -->float *f = new float[nx*ny]; // Would be better !
310 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
311 // float mini = FLT_MAX, maxi = FLT_MIN;
314 std::string strSerieUID;
315 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
318 currentTime = timeStart;
319 for (int l=0; l<nf; l++) // Loop on the frames
321 //std::cout << "Frame nb " << l << std::endl;
322 for( int j=0; j<ny; j++)
324 for (int i=0; i<nx; i++)
329 for (;;) //eatwhite(from);
333 std::cout << " !from.get(c) ";
338 //std::cout << " !isspace(c) ";
342 } //end eatwhite(from);
344 // trouble : when space is missing "-0.0990263-8.8778"
345 // is not interpreted as TWO values :-(
355 if ( first != 1 && c == '-' && previous != 'e')
358 //std::cout << " One more gauffre in frame:" << std::dec << l
359 // << ", line : " << j << " element " << i << std::endl;
368 val = (float)atof(str1.c_str());
369 //std::cout << " " << val;
370 *(f+ /*l*nx*ny + */j*nx+i) = val;
374 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
380 //std::cout << std::endl;
381 //std::cout << std::endl << " line nb : " << j
382 // << " line length : " << l << std::endl;
385 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
387 // values are expressed as % as a fraction, actually!)
388 // It's useless to rescale them as uint16_t : just *100
389 uint16_t *img = new uint16_t[ny*nx];
392 float div = maxi-mini;
393 std::cout << "div : " << div << " mini : " << mini << std::endl;
394 for( int k=0; k<ny*nx; k++)
396 *ptr = (uint16_t)((*tmp * 65535.0) / div);
401 int16_t *img = new int16_t[ny*nx];
403 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
404 for( int k=0; k<ny*nx; k++)
406 if(*tmp > 1.0) // artefacted pixel
408 else /// \todo FIXME : what about max threshold ?
409 *ptr = (int16_t)(*tmp *100);
411 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
416 // GDCM_NAME_SPACE::Debug::DebugOn();
418 std::ostringstream str;
419 GDCM_NAME_SPACE::File *file;
420 file = GDCM_NAME_SPACE::File::New();
422 // Set the image size
425 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
428 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
430 // Set the pixel type
432 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
434 str << 16; // may be 12 or 16 if componentSize =16
435 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
437 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
439 // Set the pixel representation // 0/1 1 : signed
440 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
442 // Set the samples per pixel // 1:Grey level, 3:RGB
443 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
446 // Set Rescale Intercept
449 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
454 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
457 // 0020 0037 DS 6 Image Orientation (Patient)
458 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
460 // 0020 0032 DS 3 Image Position (Patient)
461 char charImagePosition[256];
462 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
464 // 0020 0x1041 DS 1 Slice Location
465 sprintf(charImagePosition,"%f",float(l%nf));
466 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
468 //0008 103e LO 1 Series Description
469 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
471 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
472 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
473 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
475 //0020 0011 "IS" Series Number
476 sprintf(charImagePosition,"%04d",serieNumber);
477 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
479 //0020 0011 "IS" Instance Number
480 sprintf(charImagePosition,"%04d",imageNumber);
481 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
483 //0018 1060 "DS" Time Trigger
484 sprintf(charImagePosition,"%f",currentTime);
485 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
489 GDCM_NAME_SPACE::FileHelper *fh;
490 fh = GDCM_NAME_SPACE::FileHelper::New(file);
491 // cast is just to avoid warnings (*no* conversion)
492 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
493 fh->SetWriteModeToRaw();
494 fh->SetWriteTypeToDcmExplVR();
496 fh->SetWriteTypeToDcmExplVR();
499 sprintf(numero, "%02d", l);
500 std::string fileName = imageName + "." + numero + ".dcm";
501 std::cout << "fileName " << fileName << std::endl;
503 if( !fh->Write(fileName))
504 std::cout << "Failed for [" << fileName << "]\n"
505 << " File is unwrittable" << std::endl;
508 currentTime += temporalResolution;
511 } // end loop on frames
514 // Anatomical Images.
515 std::cout << " ========= Create Anatomical images" << std::endl;
517 strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
519 currentTime = timeStart;
521 for (int l=0; l<nf; l++) // Loop on the frames
523 //std::cout << "Frame nb " << l << std::endl;
524 for( int j=0; j<ny; j++)
527 for (int i=0; i<nx; i++)
542 val = (float)atof(str1.c_str());
543 // std::cout << " " << val;
544 *(f+ /*l*nx*ny + */j*nx+i) = val;
548 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
556 uint16_t *img = new uint16_t[ny*nx];
558 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
559 for( int k=0; k<ny*nx; k++)
561 *ptr = (int16_t)(*tmp);
565 std::ostringstream str;
566 GDCM_NAME_SPACE::File *file;
567 file = GDCM_NAME_SPACE::File::New();
569 // Set the image size
572 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
575 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
577 // Set the pixel type
579 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
581 str << 16; // may be 12 or 16 if componentSize =16
582 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
584 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
586 // Set the pixel representation // 0/1 1 : signed
587 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
589 // Set the samples per pixel // 1:Grey level, 3:RGB
590 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
593 // Set Rescale Intercept
596 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
601 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
604 // 0020 0037 DS 6 Image Orientation (Patient)
605 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
607 // 0020 0032 DS 3 Image Position (Patient)
608 char charImagePosition[256];
609 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
611 // 0020 0x1041 DS 1 Slice Location
612 sprintf(charImagePosition,"%f",float(l%nf));
613 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
615 //0008 103e LO 1 Series Description
616 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
618 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
619 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
620 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
622 //0020 0011 "IS" Series Number
623 sprintf(charImagePosition,"%04d",serieNumber+1);
624 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
626 //0020 0011 "IS" Instance Number
627 sprintf(charImagePosition,"%04d",imageNumber);
628 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
630 //0018 1060 "DS" Time Trigger
631 sprintf(charImagePosition,"%f",currentTime);
632 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
635 GDCM_NAME_SPACE::FileHelper *fh;
636 fh = GDCM_NAME_SPACE::FileHelper::New(file);
637 // cast is just to avoid warnings (*no* conversion)
638 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
639 fh->SetWriteModeToRaw();
640 fh->SetWriteTypeToDcmExplVR();
642 fh->SetWriteTypeToDcmExplVR();
645 sprintf(numero, "%02d", l);
646 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
647 std::cout << "fileName " << fileName << std::endl;
649 if( !fh->Write(fileName))
650 std::cout << "Failed for [" << fileName << "]\n"
651 << " File is unwrittable" << std::endl;
654 currentTime += temporalResolution;
657 } // end loop on frames