Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2006/07/24 16:36:45 $
   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 #include <fstream>
  20 #include <iostream>
  21 //#include <values.h>
  22
  23 #include "gdcmFile.h"
  24 #include "gdcmFileHelper.h"
  25 #include "gdcmDebug.h"
  26 #include "gdcmDirList.h"
  27 #include "gdcmUtil.h"
  28
  29 #include "gdcmArgMgr.h"
  30
  31 /**
  32   * \brief
  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
  37   */
  38
  39
  40 void Load(std::ifstream &from, std::string imageName, std::string patName,
  41           std::string strStudyUID, int serieNumber);
  42
  43 //std::ifstream& eatwhite(std::ifstream& is);
  44
  45 /*
  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):
  48 >
  49 >"These Attributes are applied according to the following pseudo-code,
  50 >where :
  51 x is the input value,
  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):
  55 >
  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
  59
  60 */
  61 /*
  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
  66
  67 The value of x is the output of the preceding step, the so-called
  68 "modality LUT", which may either be:
  69
  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)
  73
  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
  77
  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
  80   value
  81
  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.
  88
  89 The terms brightness and contrast are not used in radiology imaging
  90 -the window center and width are used instead.
  91 */
  92
  93 int main(int argc, char *argv[])
  94 {
  95    START_USAGE(usage)
  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                         ",
 101    " usage:                                                                   ",
 102    " DenseMultiframeToDicom dirin=rootDirectoryName                           ",
 103    "              [studyUID = ] [patName = ]                                  ",
 104    "              [listonly] [verbose] [debug]                                ",
 105    "                                                                          ",
 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'         ",
 110    FINISH_USAGE
 111
 112    // ----- Initialize Arguments Manager ------
 113
 114    gdcm::ArgMgr *am = new gdcm::ArgMgr(argc, argv);
 115
 116    if (argc == 1 || am->ArgMgrDefined("usage"))
 117    {
 118       am->ArgMgrUsage(usage); // Display 'usage'
 119       delete am;
 120       return 0;
 121    }
 122
 123    const char *dirNamein;
 124    dirNamein  = am->ArgMgrGetString("dirin",".");
 125
 126    if (am->ArgMgrDefined("debug"))
 127       gdcm::Debug::DebugOn();
 128
 129    int verbose  = am->ArgMgrDefined("verbose");
 130    int listonly = am->ArgMgrDefined("listonly");
 131    std::string patName = am->ArgMgrGetString("patname", "g^PatientName");
 132
 133    bool userDefinedStudy = am->ArgMgrDefined("studyUID");
 134    const char *studyUID  = am->ArgMgrGetString("studyUID");
 135
 136 // not described *on purpose* in the Usage !
 137    bool userDefinedSerie = am->ArgMgrDefined("serieUID");
 138    const char *serieUID  = am->ArgMgrGetString("serieUID");
 139
 140
 141    // if unused Param we give up
 142    if ( am->ArgMgrPrintUnusedLabels() )
 143    {
 144       am->ArgMgrUsage(usage);
 145       delete am;
 146       return 0;
 147    }
 148    delete am;  // we don't need Argument Manager any longer
 149
 150    // ----- Begin Processing -----
 151
 152    if ( ! gdcm::DirList::IsDirectory(dirNamein) )
 153    {
 154       std::cout << "KO : [" << dirNamein << "] is not a Directory."
 155                 << std::endl;
 156       return 0;
 157    }
 158    else
 159    {
 160       std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
 161    }
 162
 163    std::string strDirNamein(dirNamein);
 164    gdcm::DirList dirList(strDirNamein, true); // (recursively) the list of files
 165
 166    if (listonly)
 167    {
 168       std::cout << "------------List of found files ------------" << std::endl;
 169       dirList.Print();
 170       std::cout << std::endl;
 171       return 0;
 172     }
 173
 174    std::string filenameout;
 175
 176    std::string strStudyUID;
 177    strStudyUID =  gdcm::Util::CreateUniqueUID();
 178    int serieNumber =0;
 179    gdcm::DirListType fileNames;
 180    fileNames = dirList.GetFilenames();
 181    for (gdcm::DirListType::iterator it = fileNames.begin();
 182                                     it != fileNames.end();
 183                                   ++it)
 184    {
 185       std::ifstream from( (*it).c_str() );
 186       if ( !from )
 187       {
 188          std::cout << "Can't open file" << *it << std::endl;
 189          //return 0;
 190       }
 191       else
 192       {
 193          std::cout << "Success in open file" << *it << std::endl;
 194          Load(from, *it, patName, strStudyUID, serieNumber);
 195          serieNumber+=2;
 196          //return 0;
 197       }
 198    }
 199 }
 200
 201
 202 void Load(std::ifstream &from, std::string imageName, std::string patName,
 203           std::string strStudyUID, int serieNumber)
 204 {
 205    if (!from)
 206       return;
 207
 208 /* was OK for single frames
 209 eg :
 210 ---------------------------
 211  Array dimensions = 58 x 56
 212 The following is the 2D array of peak circumferential strain map, all zero value
 213 pixels are outside the mask
 214      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 215 -----------------------------
 216    std::string str1;
 217    int nx, ny;
 218    // Get nx, ny
 219    from >> str1;
 220    from >> str1;
 221    from >> str1;
 222    from >> nx;
 223    from >> str1;
 224    from >> ny;
 225
 226    std::cout << "nx " << nx << " ny " << ny << std::endl;
 227
 228    // Skip 2 lines.
 229    std::getline(from, str1);
 230    std::cout << "[" << str1 << "]" << 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  */
 236
 237  /* now, we deal with multiframes
 238  eg :
 239  ------------------------------------------
 240 X dim, Y dim, N of frames = 52x59x14
 241 Temporal resolution = 31.9200 ms
 242 First frame starts at 47.9600 ms
 243 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
 244 All pixels with zero strain values are outside the masks.
 245      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 246 --------------------------------------------
 247 */
 248    std::string str1;
 249    int nx, ny, nf;
 250     from >> str1; // X dim,
 251     from >> str1;
 252     from >> str1; // Y dim,
 253     from >> str1;
 254     from >> str1; // N of frames =
 255     from >> str1;
 256     from >> str1;
 257     from >> str1;
 258     from >> str1; // 52x59x14
 259
 260     std::cout << "[" << str1 << "]" << std::endl;
 261
 262     sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
 263     std::cout << nx << " " << ny << " " << nf << std::endl;
 264
 265        // Skip 4 lines.
 266     for (int k=0; k<5; k++) // 5
 267     {  std::cout << "Comment line number : " << k << std::endl;
 268        std::getline(from, str1);
 269        std::cout << str1 << std::endl;
 270     }
 271
 272    //float *f = new float(nx*ny);
 273    float *f = (float *) malloc(nx*ny*nf*sizeof(float));
 274   // float mini = FLT_MAX, maxi = FLT_MIN;
 275    float val;
 276
 277
 278    std::string strSerieUID;
 279    strSerieUID =  gdcm::Util::CreateUniqueUID();
 280
 281   for (int l=0; l<nf; l++)  // Loop on the frames
 282   {
 283
 284       std::cout << "Frame nb " << l << std::endl;
 285      for( int j=0; j<ny; j++)
 286      {
 287      std::cout << " "<< std::endl;
 288       //int l = 0;
 289       for (int i=0; i<nx; i++)
 290       {
 291      // std::cout << "--------------------------" << std::endl;
 292          str1="";
 293          //eatwhite(from);
 294          char c;
 295          for (;;) //eatwhite(from);
 296          {
 297             if (!from.get(c))
 298             {
 299                std::cout << " !from.get(c) ";
 300                break;
 301             }
 302             if (!isspace(c) )
 303             {
 304                //std::cout << " !isspace(c) ";
 305                from.putback(c);
 306                break;
 307             }
 308          } //end eatwhite(from);
 309     //  std::cout << "------end eatwhite-----" << std::endl;
 310
 311         // from >> str1;
 312         // trouble : whe space is missing "-0.0990263-8.8778"
 313         // is not interpreted as TWO values  :-(
 314
 315          int first = 1;
 316          char previous = 'z';
 317          for(;;)
 318          {
 319             from.get(c);
 320             //std::cout << "[" << c << "]" << std::endl;
 321             if ( c == ' ')
 322                break;
 323             if ( first != 1 && c == '-' && previous != 'e')
 324             {
 325                from.putback(c);
 326                std::cout << " One more gauffre in frame:" << std::dec << l
 327                          << ", line : " << j << " element " << i << std::endl;
 328                break;
 329              }
 330
 331              first=0;
 332              previous = c;
 333              str1=str1+c;
 334          }
 335
 336          val = (float)atof(str1.c_str());
 337          std::cout << "  " << val;
 338          *(f+ /*l*nx*ny + */j*nx+i) = val;
 339
 340         if(from.eof())
 341         {
 342             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 343                       << std::endl;
 344            break;
 345          }
 346          //l++;
 347       }
 348
 349       //std::cout << std::endl;
 350       //std::cout << std::endl << " line nb : " << j
 351       //          << " line length : " << l << std::endl;
 352
 353     }
 354
 355    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 356 /*
 357 // values are expressed as % as a fraction, actually!)
 358 // It's useless to rescale them as uint16_t : just *100
 359     uint16_t *img = new uint16_t[ny*nx];
 360     uint16_t *ptr = img;
 361     float *tmp = f;
 362     float div = maxi-mini;
 363     std::cout << "div : " << div << " mini : " << mini << std::endl;
 364     for( int k=0; k<ny*nx; k++)
 365     {
 366        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 367        tmp ++;
 368        ptr++;
 369     }
 370 */
 371     int16_t *img = new int16_t[ny*nx];
 372     int16_t *ptr = img;
 373     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 374     for( int k=0; k<ny*nx; k++)
 375     {
 376        if(*tmp > 1.0) // artefacted pixel
 377           *ptr = 0;
 378        else        /// \todo FIXME : what about max threshold ?
 379         *ptr = (int16_t)(*tmp *100);
 380
 381       // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 382        tmp ++;
 383        ptr++;
 384     }
 385
 386  // gdcm::Debug::DebugOn();
 387
 388         std::ostringstream str;
 389         gdcm::File *file;
 390         file = gdcm::File::New();
 391
 392   // Set the image size
 393         str.str("");
 394         str << nx;
 395         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 396         str.str("");
 397         str << ny;
 398         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 399
 400   // Set the pixel type
 401   //      16; //8, 16, 32
 402         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 403         str.str("");
 404         str << 16; // may be 12 or 16 if componentSize =16
 405         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 406
 407         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 408
 409   // Set the pixel representation // 0/1 1 : signed
 410         file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 411
 412   // Set the samples per pixel // 1:Grey level, 3:RGB
 413         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 414
 415 /*
 416   // Set Rescale Intercept
 417         str.str("");
 418         str << div;
 419         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 420
 421   // Set Rescale Slope
 422         str.str("");
 423         str << mini;
 424         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 425 */
 426
 427 // 0020 0037 DS 6 Image Orientation (Patient)
 428          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 429
 430 // 0020 0032 DS 3 Image Position (Patient)
 431         char charImagePosition[256];
 432         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 433
 434 // 0020 0x1041 DS 1 Slice Location
 435         sprintf(charImagePosition,"%f",float(l%nf));
 436         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 437
 438 //0008 103e LO 1 Series Description
 439         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 440
 441         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 442         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 443         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 444
 445 //0020 0011 "IS" Series Number
 446          sprintf(charImagePosition,"%04d",serieNumber);
 447          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 448
 449    // file->Print();
 450
 451     gdcm::FileHelper *fh;
 452     fh = gdcm::FileHelper::New(file);
 453     // cast is just to avoid warnings (*no* conversion)
 454     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 455     fh->SetWriteModeToRaw();
 456     fh->SetWriteTypeToDcmExplVR();
 457
 458     fh->SetWriteTypeToDcmExplVR();
 459
 460     char numero[10];
 461     sprintf(numero, "%02d", l);
 462     std::string fileName = imageName + "." + numero + ".dcm";
 463     std::cout << "fileName " << fileName << std::endl;
 464
 465     if( !fh->Write(fileName))
 466        std::cout << "Failed for [" << fileName << "]\n"
 467                  << "           File is unwrittable" << std::endl;
 468
 469     delete img;
 470   } // end loop on frames
 471
 472     int pos = (long)(from.tellg());
 473  std::cout << std::hex << "==============================================="
 474            << pos<< std::endl;
 475
 476
 477    // Anatomical Images.
 478 std::cout << "Anatomical images" << std::endl;
 479
 480    strSerieUID =  gdcm::Util::CreateUniqueUID();
 481
 482   for (int l=0; l<nf; l++)  // Loop on the frames
 483   {
 484 std::cout << "Frame nb " << l << std::endl;
 485      for( int j=0; j<ny; j++)
 486      {
 487       int l =0;
 488       for (int i=0; i<nx; i++)
 489       {
 490          //eatwhite(from);
 491          char c;
 492          for (;;)
 493          {
 494             if (!from.get(c))
 495                break;
 496             if (!isspace(c))
 497             {
 498                from.putback(c);
 499                break;
 500             }
 501          }
 502          from >> str1;
 503          val = (float)atof(str1.c_str());
 504         // std::cout << "  " << val;
 505          *(f+ /*l*nx*ny + */j*nx+i) = val;
 506
 507         if(from.eof())
 508         {
 509             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 510                       << std::endl;
 511            break;
 512          }
 513          l++;
 514       }
 515      }
 516
 517     uint16_t *img = new uint16_t[ny*nx];
 518     uint16_t *ptr = img;
 519     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 520     for( int k=0; k<ny*nx; k++)
 521     {
 522        *ptr = (int16_t)(*tmp *100);
 523        tmp ++;
 524        ptr++;
 525     }
 526
 527         std::ostringstream str;
 528         gdcm::File *file;
 529         file = gdcm::File::New();
 530
 531   // Set the image size
 532         str.str("");
 533         str << nx;
 534         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 535         str.str("");
 536         str << ny;
 537         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 538
 539   // Set the pixel type
 540   //      16; //8, 16, 32
 541         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 542         str.str("");
 543         str << 16; // may be 12 or 16 if componentSize =16
 544         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 545
 546         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 547
 548   // Set the pixel representation // 0/1 1 : signed
 549         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 550
 551   // Set the samples per pixel // 1:Grey level, 3:RGB
 552         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 553
 554 /*
 555   // Set Rescale Intercept
 556         str.str("");
 557         str << div;
 558         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 559
 560   // Set Rescale Slope
 561         str.str("");
 562         str << mini;
 563         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 564 */
 565
 566 // 0020 0037 DS 6 Image Orientation (Patient)
 567          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 568
 569 // 0020 0032 DS 3 Image Position (Patient)
 570         char charImagePosition[256];
 571         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 572
 573 // 0020 0x1041 DS 1 Slice Location
 574         sprintf(charImagePosition,"%f",float(l%nf));
 575         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 576
 577 //0008 103e LO 1 Series Description
 578         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 579
 580         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 581         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 582         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 583
 584 //0020 0011 "IS" Series Number
 585          sprintf(charImagePosition,"%04d",serieNumber+1);
 586          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 587
 588    // file->Print();
 589
 590     gdcm::FileHelper *fh;
 591     fh = gdcm::FileHelper::New(file);
 592     // cast is just to avoid warnings (*no* conversion)
 593     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 594     fh->SetWriteModeToRaw();
 595     fh->SetWriteTypeToDcmExplVR();
 596
 597     fh->SetWriteTypeToDcmExplVR();
 598
 599     char numero[10];
 600     sprintf(numero, "%02d", l);
 601     std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
 602     std::cout << "fileName " << fileName << std::endl;
 603
 604     if( !fh->Write(fileName))
 605        std::cout << "Failed for [" << fileName << "]\n"
 606                  << "           File is unwrittable" << std::endl;
 607
 608     delete img;
 609
 610
 611
 612   } // end loop on frames
 613
 614    from.close();
 615 }