Example/RawToDicomStack.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: RawToDicomStack.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2009/01/14 14:07:40 $
   7   Version:   $Revision: 1.3 $
   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
 142    const char *pixelTypeOut = am->ArgMgrGetString("pixeltypeout", pixelType);
 143
 144    bool monochrome1 = ( 0 != am->ArgMgrDefined("monochrome1") );
 145
 146    if (am->ArgMgrDefined("debug"))
 147       GDCM_NAME_SPACE::Debug::DebugOn();
 148
 149    bool userDefinedStudy = ( 0 != am->ArgMgrDefined("studyUID") );
 150    const char *studyUID;
 151    if (userDefinedStudy)
 152       studyUID  = am->ArgMgrGetString("studyUID");
 153
 154    // not described *on purpose* in the Usage !
 155    bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );
 156    const char *serieUID;
 157    if(userDefinedSerie)
 158       serieUID = am->ArgMgrGetString("serieUID");
 159
 160    /* if unused Param we give up */
 161    if ( am->ArgMgrPrintUnusedLabels() )
 162    {
 163       am->ArgMgrUsage(usage);
 164       delete am;
 165       return 1;
 166    }
 167
 168    delete am;  // we don't need Argument Manager any longer
 169
 170    // ----------- End Arguments Manager ---------
 171
 172  /// \TODO Deal with all the images of a directory
 173
 174  // Read the Raw file
 175    std::ifstream *Fp = new std::ifstream(inputFileName, std::ios::in | std::ios::binary);
 176    if ( ! *Fp )
 177    {
 178       std::cout << "Cannot open file: " << inputFileName;
 179       delete Fp;
 180       Fp = 0;
 181       return 0;
 182    }
 183
 184    bool bigEndian = GDCM_NAME_SPACE::Util::IsCurrentProcessorBigEndian();
 185
 186    std::string strPixelType(pixelType);
 187    int pixelSign;
 188    int pixelSize;
 189    int pixelTypeCode; // for the switch case
 190
 191    if (strPixelType == "8S")
 192    {
 193       pixelSize = 1;
 194       pixelSign = 1;
 195       pixelTypeCode = -8;
 196    }
 197    else if (strPixelType == "8U")
 198    {
 199       pixelSize = 1;
 200       pixelSign = 0;
 201       pixelTypeCode = 8;
 202    }
 203    else if (strPixelType == "16S")
 204    {
 205       pixelSize = 2;
 206       pixelSign = 1;
 207       pixelTypeCode = -16;
 208    }
 209    else if (strPixelType == "16U")
 210    {
 211       pixelSize = 2;
 212       pixelSign = 0;
 213       pixelTypeCode = 16;
 214    }
 215    else if (strPixelType == "32S")
 216    {
 217       pixelSize = 4;
 218       pixelSign = 1;
 219       pixelTypeCode = -32;
 220    }
 221    else if (strPixelType == "32U")
 222    {
 223       pixelSize = 4;
 224       pixelSign = 0;
 225       pixelTypeCode = 32;
 226    }
 227    else if (strPixelType == "32F")
 228    {
 229       pixelSize = 4;
 230       pixelSign = 0;
 231       pixelTypeCode = 33;
 232    }
 233
 234    else if (strPixelType == "64D")
 235    {
 236       pixelSize = 8;
 237       pixelSign = 0;
 238       pixelTypeCode = 64;
 239    }
 240    else
 241    {
 242       std::cout << "Wrong 'pixeltype' (" << strPixelType << ")" << std::endl;
 243       return 1;
 244    }
 245
 246    std::string strPixelTypeOut(pixelTypeOut);
 247    int pixelSignOut;
 248    int pixelSizeOut;
 249    int pixelTypeOutCode; // for the switch case
 250
 251    if (strPixelTypeOut == "8S")
 252    {
 253       pixelSizeOut = 1;
 254       pixelSignOut = 1;
 255       pixelTypeOutCode = -8;
 256    }
 257    else if (strPixelTypeOut == "8U")
 258    {
 259       pixelSizeOut = 1;
 260       pixelSignOut = 0;
 261       pixelTypeOutCode = 8;
 262    }
 263    else if (strPixelTypeOut == "16S")
 264    {
 265       pixelSizeOut = 2;
 266       pixelSignOut = 1;
 267       pixelTypeOutCode = -16;
 268    }
 269    else if (strPixelTypeOut == "16U")
 270    {
 271       pixelSizeOut = 2;
 272       pixelSignOut = 0;
 273       pixelTypeOutCode = 16;
 274    }
 275    else if (strPixelTypeOut == "32S")
 276    {
 277       pixelSizeOut = 4;
 278       pixelSignOut = 1;
 279       pixelTypeOutCode = -32;
 280    }
 281    else if (strPixelTypeOut == "32U")
 282    {
 283       pixelSizeOut = 4;
 284       pixelSignOut = 0;
 285       pixelTypeOutCode = 32;
 286    }
 287    else
 288    {
 289       std::cout << "Wrong 'pixeltypeout' (" << strPixelTypeOut << ")" << std::endl;
 290       return 1;
 291    }
 292
 293    std::string strStudyUID;
 294    std::string strSerieUID;
 295
 296    if (userDefinedStudy)
 297       strSerieUID = studyUID;
 298    else
 299       strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
 300
 301    if (userDefinedSerie)
 302      strSerieUID = serieUID;
 303    else
 304       strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
 305
 306    // Read the pixels
 307
 308    int singlePlaneDataSize =  nX*nY*samplesPerPixel*pixelSizeOut;
 309    int dataSizeIn          =  nX*nY*samplesPerPixel*pixelSize*nZ;
 310
 311    uint8_t *pixels         = new uint8_t[dataSizeIn];
 312    uint8_t *planePixelsOut = new uint8_t[singlePlaneDataSize];
 313
 314    Fp->read((char*)pixels, (size_t)dataSizeIn);
 315
 316    if ( pixelSize !=1 && ( (l && bigEndian) || (b && ! bigEndian) ) )
 317    {
 318       ConvertSwapZone(pixelSize, pixels, dataSizeIn);
 319    }
 320
 321 // iterate on the planes.
 322
 323    char outputFileName[200];
 324
 325    for(int nbPlanes=0; nbPlanes<nZ; nbPlanes++)
 326    {
 327
 328       sprintf(outputFileName, "%s.%04d.dcm", outputSkeletonFileName,nbPlanes);
 329
 330 // Copy (and convert) pixels of a single plane
 331
 332      switch ( pixelTypeCode )
 333      {
 334        case 8    : CFR(PU8);  break;
 335        case -8   : CFR(PS8);  break;
 336        case -16  : CFR(PU16); break;
 337        case 16   : CFR(PS16); break;
 338        case -32  : CFR(PS32); break;
 339        case 32   : CFR(PU32); break;
 340        case 33   : CFR(PF32); break;
 341        case 64   : CFR(PD64); break;
 342      }
 343
 344 // Create an empty FileHelper
 345
 346    GDCM_NAME_SPACE::FileHelper *fileH = GDCM_NAME_SPACE::FileHelper::New();
 347
 348  // Get the (empty) image header.
 349    GDCM_NAME_SPACE::File *fileToBuild = fileH->GetFile();
 350
 351    // 'Study Instance UID'
 352    // The user is allowed to create his own Study,
 353    //          keeping the same 'Study Instance UID' for various images
 354    // The user may add images to a 'Manufacturer Study',
 355    //          adding new Series to an already existing Study
 356
 357    fileToBuild->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");  //  Study UID
 358
 359    // 'Serie Instance UID'
 360    // The user is allowed to create his own Series,
 361    // keeping the same 'Serie Instance UID' for various images
 362    // The user shouldn't add any image to a 'Manufacturer Serie'
 363    // but there is no way no to prevent him for doing that
 364
 365    fileToBuild->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");  //  Serie UID
 366
 367    std::ostringstream str;
 368
 369    // Set the image size
 370    str.str("");
 371    str << nX;
 372    fileToBuild->InsertEntryString(str.str(),0x0028,0x0011, "US"); // Columns
 373    str.str("");
 374    str << nY;
 375    fileToBuild->InsertEntryString(str.str(),0x0028,0x0010, "US"); // Rows
 376
 377    // Set the pixel type
 378
 379    str.str("");
 380    str << pixelSizeOut*8;
 381    fileToBuild->InsertEntryString(str.str(),0x0028,0x0100, "US"); // Bits Allocated
 382
 383    str.str("");
 384    str << pixelSizeOut*8;
 385    fileToBuild->InsertEntryString(str.str(),0x0028,0x0101, "US"); // Bits Stored
 386
 387    str.str("");
 388    str << ( pixelSizeOut*8 - 1 );
 389    fileToBuild->InsertEntryString(str.str(),0x0028,0x0102, "US"); // High Bit
 390
 391    str.str("");
 392    str << pixelSign;
 393    fileToBuild->InsertEntryString(str.str(),0x0028,0x0103, "US"); // Pixel Representation
 394
 395 // If you deal with a Serie of images, as slices of a volume,
 396 // it up to you to tell gdcm, for each image, what are the values of :
 397 //
 398 // 0020 0032 DS 3 Image Position (Patient)
 399 // 0020 0037 DS 6 Image Orientation (Patient)
 400
 401    str.str("");
 402    str << "0.0 \\ 0.0 \\" <<  nbPlanes;
 403    fileToBuild->InsertEntryString(str.str(),0x0020,0x0032, "DS");
 404
 405    str.str("");
 406    str << samplesPerPixel;
 407    fileToBuild->InsertEntryString(str.str(),0x0028,0x0002, "US"); // Samples per Pixel
 408
 409    if (strlen(patientName) != 0)
 410       fileToBuild->InsertEntryString(patientName,0x0010,0x0010, "PN"); // Patient's Name
 411
 412    //  0=white
 413    if(monochrome1)
 414       fileH->SetPhotometricInterpretationToMonochrome1();
 415
 416 // Set the image Pixel Data (plane)
 417    fileH->SetImageData(planePixelsOut,singlePlaneDataSize);
 418
 419 // Set the writting mode and write the image
 420    fileH->SetWriteModeToRaw();
 421
 422  // Write a DICOM Explicit VR file
 423    fileH->SetWriteTypeToDcmExplVR();
 424
 425    if( !fileH->Write(outputFileName ) )
 426    {
 427       std::cout << "Failed for [" << outputFileName << "]\n"
 428                 << "           File is unwrittable\n";
 429    }
 430
 431    fileH->Delete();
 432
 433 } // end iterate on the planes
 434
 435    delete[] pixels;
 436    delete[] planePixelsOut;
 437    return 1;
 438 }