]> Creatis software - gdcm.git/blob - Example/RawToDicomStack.cxx
Add some explanations
[gdcm.git] / Example / RawToDicomStack.cxx
1 /*=========================================================================
2                                                                                 
3   Program:   gdcm
4   Module:    $RCSfile: RawToDicomStack.cxx,v $
5   Language:  C++
6   Date:      $Date: 2009/01/19 17:05:13 $
7   Version:   $Revision: 1.4 $
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.
12                                                                                 
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.
16                                       
17 =========================================================================*/
18
19 /**
20  * Writes a Dicom file from a Raw File
21  * User has to supply parameters. 
22  * 
23  */
24 #include "gdcmFile.h"
25 #include "gdcmFileHelper.h"
26 #include "gdcmDebug.h"
27 #include "gdcmUtil.h"
28 #include "gdcmArgMgr.h"
29
30 #include <iostream>
31 #include <sstream>
32
33 typedef char               * PS8;
34 typedef unsigned char      * PU8;
35 typedef short int          * PS16;
36 typedef unsigned short int * PU16;
37 typedef int                * PS32;
38 typedef unsigned int       * PU32;
39 typedef float              * PF32;
40 typedef double             * PD64;
41
42 #define CRR(t1,t2)   { for(int i=0;i<nY;i++)                   \
43                          for(int j=0;j<nX;j++)                 \
44                            for(int k=0;k<samplesPerPixel;k++)  {\
45                                *((t2)planePixelsOut + i*nX +j+k) = *(((t1)pixels)+ nbPlanes*nX*nY + i*nX +j+k);\
46                            }                                                                                  \
47                      }
48     
49 #define CFR(PPt)  switch ( pixelTypeOutCode ) {     \
50           case -8   : CRR(PPt,PS8);  break;         \
51           case 8    : CRR(PPt,PU8);  break;         \
52           case -16  : CRR(PPt,PS16); break;         \
53           case 16   : CRR(PPt,PU16); break;         \
54           case -32  : CRR(PPt,PS32); break;         \
55           case 32   : CRR(PPt,PU32); break;         \
56           case 33   : CRR(PPt,PF32); break;         \
57           case 64   : CRR(PPt,PD64); break;         \
58           }
59
60
61 void ConvertSwapZone(int pixelSize, void *Raw, size_t RawSize);
62
63 void ConvertSwapZone(int pixelSize, void *Raw, size_t RawSize)
64 {
65    unsigned int i;    
66    if ( pixelSize == 2 )
67    {
68       uint16_t *im16 = (uint16_t*)Raw;
69       for( i = 0; i < RawSize / 2; i++ )
70       {
71          im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
72       }
73    }
74    else if ( pixelSize == 4 )
75    {
76       uint32_t s32;
77       uint16_t high;
78       uint16_t low;
79       uint32_t *im32 = (uint32_t*)Raw;
80
81       for( i = 0; i < RawSize / 4; i++ )
82       {
83          low     = im32[i] & 0x0000ffff; // 3412
84          high    = im32[i] >> 16;
85          s32     = low;
86          im32[i] = ( s32 << 16 ) | high;
87       }
88       
89    }
90 }
91
92 int main(int argc, char *argv[])
93 {
94    START_USAGE(usage)
95    " \n RawToDicom : \n                                                       ",
96    " Writes a Dicom file from a Raw File                                      ",
97    " usage: RawToDicom filein=inputFileName                                   ",
98    "                   fileout=outputSkeletonName                             ",
99    "                   rows=nb of Rows                                        ",
100    "                   lines=nb of Lines,                                     ",
101    "                   [frames = nb of Frames] //defaulted to 1               ",
102    "                   pixeltype={8U|8S|16U|16S|32U|32S|32F|64D}              ",
103    "                   pixeltypeout={8U|8S|16U|16S|32U|32S}                   ",
104    "                   [{b|l}] b:BigEndian,l:LittleEndian default : l         ",
105    "                   [samples = {1|3}}       //(1:Gray,3:RGB) defaulted to 1",
106    "                   [monochrome1]                                          ",
107    "                   [studyid = ] [patientname = Patient's name]            ",
108    "                   [debug]                                                ",
109    "                                                                          ",
110    "  monochrome1 = user wants MONOCHROME1 photom. interp. (0=white)          ",
111    "  studyUID   : *aware* user wants to add the serie                        ",
112    "                                             to an already existing study ",
113    "  debug      : developper wants to run the program in 'debug mode'        ",
114    FINISH_USAGE
115    
116
117    // ------------ Initialize Arguments Manager ----------------  
118    GDCM_NAME_SPACE::ArgMgr *am= new GDCM_NAME_SPACE::ArgMgr(argc, argv);
119   
120    if (argc == 1 || am->ArgMgrDefined("usage") )
121    {
122       am->ArgMgrUsage(usage); // Display 'usage'
123       delete am;
124       return 1;
125    }
126
127    const char *inputFileName          = am->ArgMgrGetString("filein");
128    const char *outputSkeletonFileName = am->ArgMgrGetString("fileout");
129    
130    const char *patientName = am->ArgMgrGetString("patientname", "g^Fantomas");
131    
132    int nX = am->ArgMgrWantInt("rows", usage);
133    int nY = am->ArgMgrWantInt("lines", usage);
134    int nZ = am->ArgMgrGetInt("frames", 1);
135    int samplesPerPixel = am->ArgMgrGetInt("samples", 1);
136    
137    int b = am->ArgMgrDefined("b");
138    int l = am->ArgMgrDefined("l");
139       
140    char *pixelType = am->ArgMgrWantString("pixeltype", usage);
141    const char *pixelTypeOut = am->ArgMgrGetString("pixeltypeout", pixelType);   
142
143    bool monochrome1 = ( 0 != am->ArgMgrDefined("monochrome1") );
144    
145    if (am->ArgMgrDefined("debug"))
146       GDCM_NAME_SPACE::Debug::DebugOn();
147
148    bool userDefinedStudy = ( 0 != am->ArgMgrDefined("studyUID") );
149    const char *studyUID;
150    if (userDefinedStudy)
151       studyUID  = am->ArgMgrGetString("studyUID");  
152
153    // not described *on purpose* in the Usage !
154    bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );       
155    const char *serieUID;
156    if(userDefinedSerie)
157       serieUID = am->ArgMgrGetString("serieUID");
158
159    /* if unused Param we give up */
160    if ( am->ArgMgrPrintUnusedLabels() )
161    {
162       am->ArgMgrUsage(usage);
163       delete am;
164       return 1;
165    } 
166
167    delete am;  // we don't need Argument Manager any longer
168
169    // ----------- End Arguments Manager ---------
170    
171  /// \TODO Deal with all the images of a directory
172
173  // Read the Raw file
174    std::ifstream *Fp = new std::ifstream(inputFileName, std::ios::in | std::ios::binary);
175    if ( ! *Fp )
176    {
177       std::cout << "Cannot open file: " << inputFileName;
178       delete Fp;
179       Fp = 0;
180       return 0;
181    }
182
183    bool bigEndian = GDCM_NAME_SPACE::Util::IsCurrentProcessorBigEndian();
184
185    std::string strPixelType(pixelType);
186    int pixelSign;
187    int pixelSize;
188    int pixelTypeCode; // for the switch case
189    
190    if (strPixelType == "8S")
191    {
192       pixelSize = 1;
193       pixelSign = 1;
194       pixelTypeCode = -8;
195    }
196    else if (strPixelType == "8U")
197    {
198       pixelSize = 1;
199       pixelSign = 0;
200       pixelTypeCode = 8;
201    }
202    else if (strPixelType == "16S")
203    {
204       pixelSize = 2;
205       pixelSign = 1;
206       pixelTypeCode = -16;
207    }   
208    else if (strPixelType == "16U")
209    {
210       pixelSize = 2;
211       pixelSign = 0;
212       pixelTypeCode = 16;
213    }      
214    else if (strPixelType == "32S")
215    {
216       pixelSize = 4;
217       pixelSign = 1;
218       pixelTypeCode = -32;
219    }   
220    else if (strPixelType == "32U")
221    {
222       pixelSize = 4;
223       pixelSign = 0;
224       pixelTypeCode = 32;
225    }
226    else if (strPixelType == "32F")
227    {
228       pixelSize = 4;
229       pixelSign = 0;
230       pixelTypeCode = 33;
231    }
232
233    else if (strPixelType == "64D")
234    {
235       pixelSize = 8;
236       pixelSign = 0;
237       pixelTypeCode = 64;
238    }
239    else
240    {
241       std::cout << "Wrong 'pixeltype' (" << strPixelType << ")" << std::endl;
242       return 1;
243    }
244
245    std::string strPixelTypeOut(pixelTypeOut);
246    int pixelSignOut;
247    int pixelSizeOut;
248    int pixelTypeOutCode; // for the switch case
249     
250    if (strPixelTypeOut == "8S")
251    {
252       pixelSizeOut = 1;
253       pixelSignOut = 1;
254       pixelTypeOutCode = -8;
255    }
256    else if (strPixelTypeOut == "8U")
257    {
258       pixelSizeOut = 1;
259       pixelSignOut = 0;
260       pixelTypeOutCode = 8;
261    }
262    else if (strPixelTypeOut == "16S")
263    {
264       pixelSizeOut = 2;
265       pixelSignOut = 1; 
266       pixelTypeOutCode = -16;
267    }   
268    else if (strPixelTypeOut == "16U")
269    {
270       pixelSizeOut = 2;
271       pixelSignOut = 0;
272       pixelTypeOutCode = 16;
273    }      
274    else if (strPixelTypeOut == "32S")
275    {
276       pixelSizeOut = 4;
277       pixelSignOut = 1;
278       pixelTypeOutCode = -32;
279    }   
280    else if (strPixelTypeOut == "32U")
281    {
282       pixelSizeOut = 4;
283       pixelSignOut = 0;
284       pixelTypeOutCode = 32;
285    }
286    else
287    {
288       std::cout << "Wrong 'pixeltypeout' (" << strPixelTypeOut << ")" << std::endl;
289       return 1;
290    }
291
292    std::string strStudyUID;
293    std::string strSerieUID;
294
295    if (userDefinedStudy)
296       strSerieUID = studyUID;
297    else
298       strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
299    
300    if (userDefinedSerie)
301      strSerieUID = serieUID;
302    else
303       strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
304       
305    // Read the pixels
306
307    int singlePlaneDataSize =  nX*nY*samplesPerPixel*pixelSizeOut;
308    int dataSizeIn          =  nX*nY*samplesPerPixel*pixelSize*nZ;
309
310    uint8_t *pixels         = new uint8_t[dataSizeIn];
311    uint8_t *planePixelsOut = new uint8_t[singlePlaneDataSize];
312
313    Fp->read((char*)pixels, (size_t)dataSizeIn);
314
315    if ( pixelSize !=1 && ( (l && bigEndian) || (b && ! bigEndian) ) )
316    {  
317       ConvertSwapZone(pixelSize, pixels, dataSizeIn);
318    }
319
320 // iterate on the planes.
321
322    char outputFileName[200];
323
324    for(int nbPlanes=0; nbPlanes<nZ; nbPlanes++)
325    {
326
327       sprintf(outputFileName, "%s.%04d.dcm", outputSkeletonFileName,nbPlanes);
328
329 // Copy (and convert) pixels of a single plane
330
331      switch ( pixelTypeCode )
332      {
333        case 8    : CFR(PU8);  break;
334        case -8   : CFR(PS8);  break;
335        case -16  : CFR(PU16); break;
336        case 16   : CFR(PS16); break;
337        case -32  : CFR(PS32); break;
338        case 32   : CFR(PU32); break;
339        case 33   : CFR(PF32); break;
340        case 64   : CFR(PD64); break;
341      }
342
343 // Create an empty FileHelper
344
345    GDCM_NAME_SPACE::FileHelper *fileH = GDCM_NAME_SPACE::FileHelper::New();
346  
347  // Get the (empty) image header.
348    GDCM_NAME_SPACE::File *fileToBuild = fileH->GetFile();
349
350    // 'Study Instance UID'
351    // The user is allowed to create his own Study, 
352    //          keeping the same 'Study Instance UID' for various images
353    // The user may add images to a 'Manufacturer Study',
354    //          adding new Series to an already existing Study
355
356    fileToBuild->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");  //  Study UID   
357
358    // 'Serie Instance UID'
359    // The user is allowed to create his own Series, 
360    // keeping the same 'Serie Instance UID' for various images
361    // The user shouldn't add any image to a 'Manufacturer Serie'
362    // but there is no way no to prevent him for doing that
363    
364    fileToBuild->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");  //  Serie UID
365    
366    std::ostringstream str;
367
368    // Set the image size
369    str.str("");
370    str << nX;
371    fileToBuild->InsertEntryString(str.str(),0x0028,0x0011, "US"); // Columns
372    str.str("");
373    str << nY;
374    fileToBuild->InsertEntryString(str.str(),0x0028,0x0010, "US"); // Rows
375    
376    // Set the pixel type
377    
378    str.str("");
379    str << pixelSizeOut*8;
380    fileToBuild->InsertEntryString(str.str(),0x0028,0x0100, "US"); // Bits Allocated
381
382    str.str("");
383    str << pixelSizeOut*8;
384    fileToBuild->InsertEntryString(str.str(),0x0028,0x0101, "US"); // Bits Stored
385
386    str.str("");
387    str << ( pixelSizeOut*8 - 1 );
388    fileToBuild->InsertEntryString(str.str(),0x0028,0x0102, "US"); // High Bit
389
390    str.str("");
391    str << pixelSign;
392    fileToBuild->InsertEntryString(str.str(),0x0028,0x0103, "US"); // Pixel Representation
393
394 // If you deal with a Serie of images, as slices of a volume,
395 // it up to you to tell gdcm, for each image, what are the values of :
396 // 
397 // 0020 0032 DS 3 Image Position (Patient)
398 // 0020 0037 DS 6 Image Orientation (Patient)
399
400    str.str("");
401    str << "0.0 \\ 0.0 \\" <<  nbPlanes;
402    fileToBuild->InsertEntryString(str.str(),0x0020,0x0032, "DS");
403
404    str.str("");
405    str << samplesPerPixel;
406    fileToBuild->InsertEntryString(str.str(),0x0028,0x0002, "US"); // Samples per Pixel
407
408    if (strlen(patientName) != 0)
409       fileToBuild->InsertEntryString(patientName,0x0010,0x0010, "PN"); // Patient's Name
410     
411    //  0=white  
412    if(monochrome1)
413       fileH->SetPhotometricInterpretationToMonochrome1();
414      
415 // Set the image Pixel Data (plane)
416    fileH->SetImageData(planePixelsOut,singlePlaneDataSize);
417
418 // Set the writting mode and write the image
419    fileH->SetWriteModeToRaw();
420
421  // Write a DICOM Explicit VR file
422    fileH->SetWriteTypeToDcmExplVR();
423
424    if( !fileH->Write(outputFileName ) )
425    {
426       std::cout << "Failed for [" << outputFileName << "]\n"
427                 << "           File is unwrittable\n";
428    }
429
430    fileH->Delete();
431    
432 } // end iterate on the planes
433
434    delete[] pixels;
435    delete[] planePixelsOut;
436    return 1;
437 }