1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2006/07/26 17:02:55 $
7 Version: $Revision: 1.4 $
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 DenseMultiFramessToDicom : \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::ArgMgr *am = new gdcm::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::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::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::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::Util::CreateUniqueUID();
179 gdcm::DirListType fileNames;
180 fileNames = dirList.GetFilenames();
181 for (gdcm::DirListType::iterator it = fileNames.begin();
182 it != fileNames.end();
185 if ( gdcm::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 = (float *) malloc(nx*ny*nf*sizeof(float));
310 // float mini = FLT_MAX, maxi = FLT_MIN;
313 std::string strSerieUID;
314 strSerieUID = gdcm::Util::CreateUniqueUID();
317 currentTime = timeStart;
318 for (int l=0; l<nf; l++) // Loop on the frames
320 //std::cout << "Frame nb " << l << std::endl;
321 for( int j=0; j<ny; j++)
323 for (int i=0; i<nx; i++)
328 for (;;) //eatwhite(from);
332 std::cout << " !from.get(c) ";
337 //std::cout << " !isspace(c) ";
341 } //end eatwhite(from);
343 // trouble : when space is missing "-0.0990263-8.8778"
344 // is not interpreted as TWO values :-(
354 if ( first != 1 && c == '-' && previous != 'e')
357 //std::cout << " One more gauffre in frame:" << std::dec << l
358 // << ", line : " << j << " element " << i << std::endl;
367 val = (float)atof(str1.c_str());
368 //std::cout << " " << val;
369 *(f+ /*l*nx*ny + */j*nx+i) = val;
373 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
379 //std::cout << std::endl;
380 //std::cout << std::endl << " line nb : " << j
381 // << " line length : " << l << std::endl;
384 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
386 // values are expressed as % as a fraction, actually!)
387 // It's useless to rescale them as uint16_t : just *100
388 uint16_t *img = new uint16_t[ny*nx];
391 float div = maxi-mini;
392 std::cout << "div : " << div << " mini : " << mini << std::endl;
393 for( int k=0; k<ny*nx; k++)
395 *ptr = (uint16_t)((*tmp * 65535.0) / div);
400 int16_t *img = new int16_t[ny*nx];
402 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
403 for( int k=0; k<ny*nx; k++)
405 if(*tmp > 1.0) // artefacted pixel
407 else /// \todo FIXME : what about max threshold ?
408 *ptr = (int16_t)(*tmp *100);
410 // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
415 // gdcm::Debug::DebugOn();
417 std::ostringstream str;
419 file = gdcm::File::New();
421 // Set the image size
424 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
427 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
429 // Set the pixel type
431 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
433 str << 16; // may be 12 or 16 if componentSize =16
434 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
436 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
438 // Set the pixel representation // 0/1 1 : signed
439 file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
441 // Set the samples per pixel // 1:Grey level, 3:RGB
442 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
445 // Set Rescale Intercept
448 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
453 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
456 // 0020 0037 DS 6 Image Orientation (Patient)
457 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
459 // 0020 0032 DS 3 Image Position (Patient)
460 char charImagePosition[256];
461 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
463 // 0020 0x1041 DS 1 Slice Location
464 sprintf(charImagePosition,"%f",float(l%nf));
465 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
467 //0008 103e LO 1 Series Description
468 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
470 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
471 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
472 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
474 //0020 0011 "IS" Series Number
475 sprintf(charImagePosition,"%04d",serieNumber);
476 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
478 //0020 0011 "IS" Instance Number
479 sprintf(charImagePosition,"%04d",imageNumber);
480 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
482 //0018 1060 "DS" Time Trigger
483 sprintf(charImagePosition,"%f",currentTime);
484 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
488 gdcm::FileHelper *fh;
489 fh = gdcm::FileHelper::New(file);
490 // cast is just to avoid warnings (*no* conversion)
491 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
492 fh->SetWriteModeToRaw();
493 fh->SetWriteTypeToDcmExplVR();
495 fh->SetWriteTypeToDcmExplVR();
498 sprintf(numero, "%02d", l);
499 std::string fileName = imageName + "." + numero + ".dcm";
500 std::cout << "fileName " << fileName << std::endl;
502 if( !fh->Write(fileName))
503 std::cout << "Failed for [" << fileName << "]\n"
504 << " File is unwrittable" << std::endl;
507 currentTime += temporalResolution;
510 } // end loop on frames
513 // Anatomical Images.
514 std::cout << " ========= Create Anatomical images" << std::endl;
516 strSerieUID = gdcm::Util::CreateUniqueUID();
518 currentTime = timeStart;
520 for (int l=0; l<nf; l++) // Loop on the frames
522 //std::cout << "Frame nb " << l << std::endl;
523 for( int j=0; j<ny; j++)
526 for (int i=0; i<nx; i++)
541 val = (float)atof(str1.c_str());
542 // std::cout << " " << val;
543 *(f+ /*l*nx*ny + */j*nx+i) = val;
547 std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
555 uint16_t *img = new uint16_t[ny*nx];
557 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
558 for( int k=0; k<ny*nx; k++)
560 *ptr = (int16_t)(*tmp);
564 std::ostringstream str;
566 file = gdcm::File::New();
568 // Set the image size
571 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
574 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
576 // Set the pixel type
578 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
580 str << 16; // may be 12 or 16 if componentSize =16
581 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
583 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
585 // Set the pixel representation // 0/1 1 : signed
586 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
588 // Set the samples per pixel // 1:Grey level, 3:RGB
589 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
592 // Set Rescale Intercept
595 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
600 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
603 // 0020 0037 DS 6 Image Orientation (Patient)
604 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
606 // 0020 0032 DS 3 Image Position (Patient)
607 char charImagePosition[256];
608 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
610 // 0020 0x1041 DS 1 Slice Location
611 sprintf(charImagePosition,"%f",float(l%nf));
612 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
614 //0008 103e LO 1 Series Description
615 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
617 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
618 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
619 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
621 //0020 0011 "IS" Series Number
622 sprintf(charImagePosition,"%04d",serieNumber+1);
623 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
625 //0020 0011 "IS" Instance Number
626 sprintf(charImagePosition,"%04d",imageNumber);
627 file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
629 //0018 1060 "DS" Time Trigger
630 sprintf(charImagePosition,"%f",currentTime);
631 file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
634 gdcm::FileHelper *fh;
635 fh = gdcm::FileHelper::New(file);
636 // cast is just to avoid warnings (*no* conversion)
637 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
638 fh->SetWriteModeToRaw();
639 fh->SetWriteTypeToDcmExplVR();
641 fh->SetWriteTypeToDcmExplVR();
644 sprintf(numero, "%02d", l);
645 std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
646 std::cout << "fileName " << fileName << std::endl;
648 if( !fh->Write(fileName))
649 std::cout << "Failed for [" << fileName << "]\n"
650 << " File is unwrittable" << std::endl;
653 currentTime += temporalResolution;
656 } // end loop on frames