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 }