1 /*=========================================================================
4 Module: $RCSfile: DenseMultiFramesToDicom.cxx,v $
6 Date: $Date: 2006/07/20 17:15:28 $
7 Version: $Revision: 1.1 $
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 DenseMultiframeToDicom :\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 " DenseMultiframeToDicom 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", "g^PatientName");
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();
187 std::ifstream from( (*it).c_str() );
190 std::cout << "Can't open file" << *it << std::endl;
195 std::cout << "Success in open file" << *it << std::endl;
196 Load(from, *it, patName, strStudyUID, serieNumber);
204 void Load(std::ifstream &from, std::string imageName, std::string patName,
205 std::string strStudyUID, int serieNumber)
210 /* was OK for single frames
212 ---------------------------
213 Array dimensions = 58 x 56
214 The following is the 2D array of peak circumferential strain map, all zero value
215 pixels are outside the mask
216 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
217 -----------------------------
228 std::cout << "nx " << nx << " ny " << ny << std::endl;
231 std::getline(from, str1);
232 std::cout << "[" << str1 << "]" << std::endl;
233 std::getline(from, str1);
234 std::cout << "[" << str1 << "]" << std::endl;
235 std::getline(from, str1);
236 std::cout << "[" << str1 << "]" << std::endl;
239 /* now, we deal with multiframes
241 ------------------------------------------
242 X dim, Y dim, N of frames = 52x59x14
243 Temporal resolution = 31.9200 ms
244 First frame starts at 47.9600 ms
245 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
246 All pixels with zero strain values are outside the masks.
247 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
248 --------------------------------------------
252 from >> str1; // X dim,
254 from >> str1; // Y dim,
256 from >> str1; // N of frames =
260 from >> str1; // 52x59x14
262 std::cout << "[" << str1 << "]" << std::endl;
264 sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
265 std::cout << nx << " " << ny << " " << nf << std::endl;
268 for (int k=0; k<4; k++)
270 std::getline(from, str1);
271 std::cout << str1 << std::endl;
274 //float *f = new float(nx*ny);
275 float *f = (float *) malloc(nx*ny*nf*sizeof(float));
276 // float mini = FLT_MAX, maxi = FLT_MIN;
280 std::string strSerieUID;
281 strSerieUID = gdcm::Util::CreateUniqueUID();
283 for (int l=0; l<nf; l++) // Loop on the frames
285 for( int j=0; j<ny; j++)
288 for (int i=0; i<nx; i++)
303 val = (float)atof(str1.c_str());
304 // std::cout << " " << val;
305 *(f+ /*l*nx*ny + */j*nx+i) = val;
309 std::cout << "Missing values at [" << j <<"," << i << "]"
316 std::cout << std::endl;
317 //std::cout << std::endl << " line nb : " << j
318 // << " line length : " << l << std::endl;
323 // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl;
325 // values are expressed as % as a fraction, actually!)
326 // It's useless to rescale them as uint16_t : just *100
327 uint16_t *img = new uint16_t[ny*nx];
330 float div = maxi-mini;
331 std::cout << "div : " << div << " mini : " << mini << std::endl;
332 for( int k=0; k<ny*nx; k++)
334 *ptr = (uint16_t)((*tmp * 65535.0) / div);
339 uint16_t *img = new uint16_t[ny*nx];
341 float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
342 for( int k=0; k<ny*nx; k++)
344 if(*tmp < 0) // artefacted pixel
346 else /// \todo FIXME : what about max threshold ?
347 *ptr = (int16_t)(*tmp *100);
349 //std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
354 // gdcm::Debug::DebugOn();
356 std::ostringstream str;
358 file = gdcm::File::New();
360 // Set the image size
363 file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
366 file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
368 // Set the pixel type
370 file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
372 str << 16; // may be 12 or 16 if componentSize =16
373 file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
375 file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
377 // Set the pixel representation // 0/1 1 : signed
378 file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
380 // Set the samples per pixel // 1:Grey level, 3:RGB
381 file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
384 // Set Rescale Intercept
387 file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
392 file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
395 // 0020 0037 DS 6 Image Orientation (Patient)
396 file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
398 // 0020 0032 DS 3 Image Position (Patient)
399 char charImagePosition[256];
400 sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
402 // 0020 0x1041 DS 1 Slice Location
403 sprintf(charImagePosition,"%f",float(l%nf));
404 file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
406 //0008 103e LO 1 Series Description
407 file->InsertEntryString(imageName,0x0008,0x103e, "LO");
409 file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
410 file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
411 file->InsertEntryString(patName,0x0010,0x0010, "PN"); // Patient's Name
413 //0020 0011 "IS" Series Number
414 sprintf(charImagePosition,"%04d",serieNumber);
415 file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
419 gdcm::FileHelper *fh;
420 fh = gdcm::FileHelper::New(file);
421 // cast is just to avoid warnings (*no* conversion)
422 fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
423 fh->SetWriteModeToRaw();
424 fh->SetWriteTypeToDcmExplVR();
426 fh->SetWriteTypeToDcmExplVR();
430 sprintf(numero, "%02d", l);
431 std::string fileName = imageName + "." + numero + ".dcm";
432 std::cout << "fileName " << fileName << std::endl;
434 if( !fh->Write(fileName))
435 std::cout << "Failed for [" << fileName << "]\n"
436 << " File is unwrittable" << std::endl;
439 } // end loop on frames