Example/RawToInTagDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: RawToInTagDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2009/01/19 17:03:38 $
   7   Version:   $Revision: 1.1 $
   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 an old style Bruker Dicom file, InTag compliant, 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 "gdcmDirList.h"
  29
  30 #include "gdcmArgMgr.h"
  31
  32 #include <iostream>
  33 #include <sstream>
  34
  35
  36 /**
  37   * \brief
  38   *          - explores  the given directory
  39   *          - orders the Raw Files
  40   *            according to their names
  41   *          - fills a single level Directory with *all* the files,
  42   *            converted into a Brucker-like Dicom, Intags compliant
  43   */
  44
  45 typedef char               * PS8;
  46 typedef unsigned char      * PU8;
  47 typedef short int          * PS16;
  48 typedef unsigned short int * PU16;
  49 typedef int                * PS32;
  50 typedef unsigned int       * PU32;
  51 typedef float              * PF32;
  52 typedef double             * PD64;
  53
  54 #define CRR(t1,t2)   { for(int l=0;l<nX*nY*/* nZ* */samplesPerPixel;l++)           \
  55                                   {*((t2)planePixelsOut + l) = *(((t1)pixels)+ l);}\
  56                      }
  57
  58 #define CFR(PPt)  switch ( pixelTypeOutCode ) {     \
  59           case -8   : CRR(PPt,PS8);  break;         \
  60           case 8    : CRR(PPt,PU8);  break;         \
  61           case -16  : CRR(PPt,PS16); break;         \
  62           case 16   : CRR(PPt,PU16); break;         \
  63           case -32  : CRR(PPt,PS32); break;         \
  64           case 32   : CRR(PPt,PU32); break;         \
  65           case 33   : CRR(PPt,PF32); break;         \
  66           case 64   : CRR(PPt,PD64); break;         \
  67           }
  68
  69
  70 void ConvertSwapZone(int pixelSize, void *Raw, size_t RawSize);
  71
  72 void ConvertSwapZone(int pixelSize, void *Raw, size_t RawSize)
  73 {
  74    unsigned int i;
  75    if ( pixelSize == 2 )
  76    {
  77       uint16_t *im16 = (uint16_t*)Raw;
  78       for( i = 0; i < RawSize / 2; i++ )
  79       {
  80          im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
  81       }
  82    }
  83    else if ( pixelSize == 4 )
  84    {
  85       uint32_t s32;
  86       uint16_t high;
  87       uint16_t low;
  88       uint32_t *im32 = (uint32_t*)Raw;
  89
  90       for( i = 0; i < RawSize / 4; i++ )
  91       {
  92          low     = im32[i] & 0x0000ffff; // 3412
  93          high    = im32[i] >> 16;
  94          s32     = low;
  95          im32[i] = ( s32 << 16 ) | high;
  96       }
  97
  98    }
  99 }
 100
 101 int main(int argc, char *argv[])
 102 {
 103    START_USAGE(usage)
 104    " \n RawToDicom : \n                                                       ",
 105    " - explores recursively the given directory,                              ",
 106    " - orders the Raw Files according to their names                          ",
 107    " - fills a single level (*) Directory with *all* the files,               ",
 108    "           converted into a Brucker-like Dicom, InTags compliant          ",
 109    " usage: RawToDicom dirin=rootDirectoryName                                ",
 110    "                   dirout=outputDirectoryName                             ",
 111    "                   rows=nb of Rows                                        ",
 112    "                   lines=nb of Lines,                                     ",
 113    "                   [frames = nb of Frames] //defaulted to 1               ",
 114    "                   pixeltype={8U|8S|16U|16S|32U|32S|32F|64D}              ",
 115    "                   pixeltypeout={8U|8S|16U|16S|32U|32S}                   ",
 116    "                   [{b|l}] b:BigEndian,l:LittleEndian default : l         ",
 117    "                   [samples = {1|3}}       //(1:Gray,3:RGB) defaulted to 1",
 118    "                   [monochrome1]                                          ",
 119    "                   [studyid = ] [patientname = Patient's name]            ",
 120    "                   [debug] [verbose] [listonly]                           ",
 121    "                                                                          ",
 122    "  monochrome1 = user wants MONOCHROME1 photom. interp. (0=white)          ",
 123    "  studyUID   : *aware* user wants to add the serie                        ",
 124    "                                             to an already existing study ",
 125    "  debug      : developper wants to run the program in 'debug mode'        ",
 126    FINISH_USAGE
 127
 128
 129    // ------------ Initialize Arguments Manager ----------------
 130    GDCM_NAME_SPACE::ArgMgr *am= new GDCM_NAME_SPACE::ArgMgr(argc, argv);
 131
 132    if (argc == 1 || am->ArgMgrDefined("usage") )
 133    {
 134       am->ArgMgrUsage(usage); // Display 'usage'
 135       delete am;
 136       return 1;
 137    }
 138
 139    const char *dirNamein;
 140    dirNamein  = am->ArgMgrGetString("dirin",".");
 141
 142    const char *dirNameout;
 143    dirNameout  = am->ArgMgrGetString("dirout",".");
 144
 145    const char *patientName = am->ArgMgrGetString("patientname", "g^Fantomas");
 146
 147    int nX = am->ArgMgrWantInt("rows", usage);
 148    int nY = am->ArgMgrWantInt("lines", usage);
 149    int nZ = am->ArgMgrGetInt("frames", 1);
 150    int samplesPerPixel = am->ArgMgrGetInt("samples", 1);
 151
 152    int b = am->ArgMgrDefined("b");
 153    int l = am->ArgMgrDefined("l");
 154
 155    char *pixelType = am->ArgMgrWantString("pixeltype", usage);
 156    const char *pixelTypeOut = am->ArgMgrGetString("pixeltypeout", pixelType);
 157
 158    bool monochrome1 = ( 0 != am->ArgMgrDefined("monochrome1") );
 159
 160    if (am->ArgMgrDefined("debug"))
 161       GDCM_NAME_SPACE::Debug::DebugOn();
 162
 163    int verbose  = am->ArgMgrDefined("verbose");
 164    int listonly = am->ArgMgrDefined("listonly");
 165
 166    bool userDefinedStudy = ( 0 != am->ArgMgrDefined("studyUID") );
 167    const char *studyUID;
 168    if (userDefinedStudy)
 169       studyUID  = am->ArgMgrGetString("studyUID");
 170
 171    // not described *on purpose* in the Usage !
 172    bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );
 173    const char *serieUID;
 174    if(userDefinedSerie)
 175       serieUID = am->ArgMgrGetString("serieUID");
 176
 177    /* if unused Param we give up */
 178    if ( am->ArgMgrPrintUnusedLabels() )
 179    {
 180       am->ArgMgrUsage(usage);
 181       delete am;
 182       return 1;
 183    }
 184
 185    delete am;  // we don't need Argument Manager any longer
 186
 187    // ----------- End Arguments Manager ---------
 188
 189
 190    // ----- Begin Processing -----
 191
 192
 193    bool bigEndian = GDCM_NAME_SPACE::Util::IsCurrentProcessorBigEndian();
 194
 195    std::string strPixelType(pixelType);
 196    int pixelSign;
 197    int pixelSize;
 198    int pixelTypeCode; // for the switch case
 199
 200    if (strPixelType == "8S")
 201    {
 202       pixelSize = 1;
 203       pixelSign = 1;
 204       pixelTypeCode = -8;
 205    }
 206    else if (strPixelType == "8U")
 207    {
 208       pixelSize = 1;
 209       pixelSign = 0;
 210       pixelTypeCode = 8;
 211    }
 212    else if (strPixelType == "16S")
 213    {
 214       pixelSize = 2;
 215       pixelSign = 1;
 216       pixelTypeCode = -16;
 217    }
 218    else if (strPixelType == "16U")
 219    {
 220       pixelSize = 2;
 221       pixelSign = 0;
 222       pixelTypeCode = 16;
 223    }
 224    else if (strPixelType == "32S")
 225    {
 226       pixelSize = 4;
 227       pixelSign = 1;
 228       pixelTypeCode = -32;
 229    }
 230    else if (strPixelType == "32U")
 231    {
 232       pixelSize = 4;
 233       pixelSign = 0;
 234       pixelTypeCode = 32;
 235    }
 236    else if (strPixelType == "32F")
 237    {
 238       pixelSize = 4;
 239       pixelSign = 0;
 240       pixelTypeCode = 33;
 241    }
 242
 243    else if (strPixelType == "64D")
 244    {
 245       pixelSize = 8;
 246       pixelSign = 0;
 247       pixelTypeCode = 64;
 248    }
 249    else
 250    {
 251       std::cout << "Wrong 'pixeltype' (" << strPixelType << ")" << std::endl;
 252       return 1;
 253    }
 254
 255    std::string strPixelTypeOut(pixelTypeOut);
 256    int pixelSignOut;
 257    int pixelSizeOut;
 258    int pixelTypeOutCode; // for the switch case
 259
 260    if (strPixelTypeOut == "8S")
 261    {
 262       pixelSizeOut = 1;
 263       pixelSignOut = 1;
 264       pixelTypeOutCode = -8;
 265    }
 266    else if (strPixelTypeOut == "8U")
 267    {
 268       pixelSizeOut = 1;
 269       pixelSignOut = 0;
 270       pixelTypeOutCode = 8;
 271    }
 272    else if (strPixelTypeOut == "16S")
 273    {
 274       pixelSizeOut = 2;
 275       pixelSignOut = 1;
 276       pixelTypeOutCode = -16;
 277    }
 278    else if (strPixelTypeOut == "16U")
 279    {
 280       pixelSizeOut = 2;
 281       pixelSignOut = 0;
 282       pixelTypeOutCode = 16;
 283    }
 284    else if (strPixelTypeOut == "32S")
 285    {
 286       pixelSizeOut = 4;
 287       pixelSignOut = 1;
 288       pixelTypeOutCode = -32;
 289    }
 290    else if (strPixelTypeOut == "32U")
 291    {
 292       pixelSizeOut = 4;
 293       pixelSignOut = 0;
 294       pixelTypeOutCode = 32;
 295    }
 296    else
 297    {
 298       std::cout << "Wrong 'pixeltypeout' (" << strPixelTypeOut << ")" << std::endl;
 299       return 1;
 300    }
 301
 302
 303
 304    if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
 305    {
 306       std::cout << "KO : [" << dirNamein << "] is not a Directory." << std::endl;
 307       return 0;
 308
 309    }
 310    else
 311    {
 312       std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
 313    }
 314
 315    std::string systemCommand;
 316
 317    std::cout << "Check for output directory :[" << dirNameout << "]."
 318              <<std::endl;
 319    if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNameout) )    // dirout not found
 320    {
 321       std::string strDirNameout(dirNameout);          // to please gcc 4
 322       systemCommand = "mkdir " +strDirNameout;        // create it!
 323       if (verbose)
 324          std::cout << systemCommand << std::endl;
 325       system (systemCommand.c_str());
 326       if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNameout) ) // be sure it worked
 327       {
 328           std::cout << "KO : not a dir : [" << dirNameout << "] (creation failure ?)" << std::endl;
 329       return 0;
 330
 331       }
 332       else
 333       {
 334         std::cout << "Directory [" << dirNameout << "] created." << std::endl;
 335       }
 336    }
 337    else
 338    {
 339        std::cout << "Output Directory [" << dirNameout << "] already exists; Used as is." << std::endl;
 340    }
 341
 342    std::string strDirNamein(dirNamein);
 343    GDCM_NAME_SPACE::DirList dirList(strDirNamein, false); // DON'T get recursively the list of files
 344    std::string strDirNameout(dirNameout);
 345
 346    if (listonly)
 347    {
 348       std::cout << "------------List of found files ------------" << std::endl;
 349       dirList.Print();
 350       std::cout << std::endl;
 351       return 1;
 352    }
 353
 354    std::string strStudyUID;
 355    std::string strSerieUID;
 356
 357    if (userDefinedStudy)
 358       strSerieUID = studyUID;
 359    else
 360       strStudyUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
 361
 362    if (userDefinedSerie)
 363      strSerieUID = serieUID;
 364    else
 365       strSerieUID = GDCM_NAME_SPACE::Util::CreateUniqueUID();
 366
 367    GDCM_NAME_SPACE::DirListType fileNames;
 368    fileNames = dirList.GetFilenames();
 369
 370    GDCM_NAME_SPACE::File *f;
 371    GDCM_NAME_SPACE::FileHelper *fh;
 372
 373    std::string outputFileName;
 374
 375    int sliceIndex  = 1;  // IMPORTANT : must start with non zero value
 376    int frameIndex  = 1;  // IMPORTANT : must start with non zero value
 377    std::string chSessionIndex("1");
 378
 379    // -----------------------------------------------------
 380    // Iterate fo ALL the files found in the input directory
 381    // -----------------------------------------------------
 382
 383    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 384                                     it != fileNames.end();
 385                                   ++it)
 386    {
 387
 388  // Open the Raw file
 389    std::ifstream *Fp = new std::ifstream((*it).c_str(), std::ios::in | std::ios::binary);
 390    if ( ! *Fp )
 391    {
 392       std::cout << "Cannot open file: [" << *it << "]" << std::endl;
 393       delete Fp;
 394       Fp = 0;
 395       return 0;
 396    }
 397    if (verbose)
 398       std::cout << "Success to open file: [" << *it << "]" << std::endl;
 399
 400    // Read the pixels
 401
 402    int singlePlaneDataSize =  nX*nY*samplesPerPixel*pixelSizeOut;
 403    int dataSizeIn          =  nX*nY*samplesPerPixel*pixelSize*nZ;
 404
 405    uint8_t *pixels         = new uint8_t[dataSizeIn];
 406    uint8_t *planePixelsOut = new uint8_t[singlePlaneDataSize];
 407
 408    Fp->read((char*)pixels, (size_t)dataSizeIn);
 409
 410    if ( pixelSize !=1 && ( (l && bigEndian) || (b && ! bigEndian) ) )
 411    {
 412       ConvertSwapZone(pixelSize, pixels, dataSizeIn);
 413    }
 414
 415 if (verbose)
 416    std::cout << "After ConvertSwapZone" << std::endl;
 417
 418 // Copy (and convert) pixels of a single plane
 419
 420      switch ( pixelTypeCode )
 421      {
 422        case 8    : CFR(PU8);  break;
 423        case -8   : CFR(PS8);  break;
 424        case -16  : CFR(PU16); break;
 425        case 16   : CFR(PS16); break;
 426        case -32  : CFR(PS32); break;
 427        case 32   : CFR(PU32); break;
 428        case 33   : CFR(PF32); break;
 429        case 64   : CFR(PD64); break;
 430      }
 431
 432 if (verbose)
 433    std::cout << "After CFR" << std::endl;
 434
 435 // Create an empty FileHelper
 436
 437    GDCM_NAME_SPACE::FileHelper *fileH = GDCM_NAME_SPACE::FileHelper::New();
 438
 439  // Get the (empty) image header.
 440    GDCM_NAME_SPACE::File *fileToBuild = fileH->GetFile();
 441
 442    // 'Study Instance UID'
 443    // The user is allowed to create his own Study,
 444    //          keeping the same 'Study Instance UID' for various images
 445    // The user may add images to a 'Manufacturer Study',
 446    //          adding new Series to an already existing Study
 447
 448    fileToBuild->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");  //  Study UID
 449
 450    // 'Serie Instance UID'
 451    // The user is allowed to create his own Series,
 452    // keeping the same 'Serie Instance UID' for various images
 453    // The user shouldn't add any image to a 'Manufacturer Serie'
 454    // but there is no way no to prevent him for doing that
 455
 456    fileToBuild->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");  //  Serie UID
 457
 458    std::ostringstream str;
 459
 460    // Set the image size
 461    str.str("");
 462    str << nX;
 463    fileToBuild->InsertEntryString(str.str(),0x0028,0x0011, "US"); // Columns
 464    str.str("");
 465    str << nY;
 466    fileToBuild->InsertEntryString(str.str(),0x0028,0x0010, "US"); // Rows
 467
 468        if (verbose)
 469          std::cout << "before  debut des choses serieuses2"<< std::endl;
 470
 471 //   str.str("");
 472 //   str << nZ;
 473 //   fileToBuild->InsertEntryString(str.str(),0x0028,0x0008, "IS"); // Number of Frames
 474
 475    // Set the pixel type
 476
 477    str.str("");
 478    str << pixelSizeOut*8;
 479    fileToBuild->InsertEntryString(str.str(),0x0028,0x0100, "US"); // Bits Allocated
 480
 481    str.str("");
 482    str << pixelSizeOut*8;
 483    fileToBuild->InsertEntryString(str.str(),0x0028,0x0101, "US"); // Bits Stored
 484
 485    str.str("");
 486    str << ( pixelSizeOut*8 - 1 );
 487    fileToBuild->InsertEntryString(str.str(),0x0028,0x0102, "US"); // High Bit
 488
 489    str.str("");
 490    str << pixelSign;
 491    fileToBuild->InsertEntryString(str.str(),0x0028,0x0103, "US"); // Pixel Representation
 492
 493 // If you deal with a Serie of images, as slices of a volume,
 494 // it's up to you to tell gdcm, for each image, what are the values of :
 495 //
 496 // 0020 0032 DS 3 Image Position (Patient)
 497 // 0020 0037 DS 6 Image Orientation (Patient)
 498
 499    str.str("");
 500    str << "0.0\\0.0\\0.0";
 501
 502 // take Frame Index as position (we don't know anything else) NO!
 503
 504    fileToBuild->InsertEntryString(str.str(),0x0020,0x0032, "DS");
 505
 506    fileToBuild->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 507
 508
 509    str.str("");
 510    str << samplesPerPixel;
 511    fileToBuild->InsertEntryString(str.str(),0x0028,0x0002, "US"); // Samples per Pixel
 512
 513    if (strlen(patientName) != 0)
 514       fileToBuild->InsertEntryString(patientName,0x0010,0x0010, "PN"); // Patient's Name
 515
 516    //  0=white
 517    if(monochrome1)
 518       fileH->SetPhotometricInterpretationToMonochrome1();
 519
 520 // Special InTag, now.
 521
 522       fileH->InsertEntryString(chSessionIndex, 0x0020, 0x0012, "IS");
 523
 524       // Deal with  0x0021, 0x1020 : 'SLICE INDEX'
 525       // will stay to 0, since the stuff deals with single slice directories
 526       char chSliceIndex[5];
 527       sprintf(chSliceIndex, "%04d", sliceIndex);
 528       std::string strChSliceIndex(chSliceIndex);
 529
 530       // Deal with  0x0021, 0x1040 : 'FRAME INDEX'
 531
 532       str.str("");
 533       str << frameIndex;
 534       frameIndex++;  // be ready for next one
 535
 536       fileH->InsertEntryString(strChSliceIndex, 0x0021, 0x1020, "IS");
 537       fileH->InsertEntryString(str.str(),    0x0021, 0x1040, "IS");
 538
 539       // Pixel Size
 540       /// \TODO Ask user to supply 'Pixel Size' value
 541       float pxSzX = 1.0;
 542       float pxSzY = 1.0;
 543
 544       // Field of view
 545       char fov[64];
 546       sprintf(fov, "%f\\%f",nX*pxSzX, nY*pxSzY);
 547       fileH->InsertEntryString(fov, 0x0019, 0x1000, "DS");
 548
 549
 550 // Set the image Pixel Data
 551    fileH->SetImageData(planePixelsOut,singlePlaneDataSize);
 552
 553 // Set the writting mode and write the image
 554    fileH->SetWriteModeToRaw();
 555
 556
 557
 558  // Write a DICOM Explicit VR file
 559    fileH->SetWriteTypeToDcmExplVR();
 560
 561
 562    outputFileName = strDirNameout +  GDCM_NAME_SPACE::GDCM_FILESEPARATOR + *it + "_ForInTag.dcm";
 563    if( !fileH->Write(outputFileName ) )
 564    {
 565       std::cout << "Failed for [" << outputFileName << "]\n"
 566                 << "           File is unwrittable\n";
 567    }
 568
 569 // End of :    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin() ...
 570
 571
 572
 573    fileH->Delete();
 574
 575    delete[] pixels;
 576    delete[] planePixelsOut;
 577
 578   }
 579    return 1;
 580  } // end of : for (GDCM_NAME_SPACE::DirListType::iterator it