Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2006/07/26 17:02:55 $
   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 #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 DenseMultiFramessToDicom :                                          \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    " DenseMultiFramesToDicom 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", dirNamein);
 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       if ( gdcm::Util::GetName((*it)).c_str()[0] == '.' )
 186       {
 187       // skip hidden files
 188          continue;
 189       }
 190       int sz = (*it).size();
 191       if ( (*it).c_str()[sz-1] != 't')
 192       {
 193          // skip non .txt files
 194          continue;
 195       }
 196       std::ifstream from( (*it).c_str() );
 197       if ( !from )
 198       {
 199          std::cout << "Can't open file" << *it << std::endl;
 200          //return 0;
 201       }
 202       else
 203       {
 204          std::cout << "Success in open file" << *it << std::endl;
 205          Load(from, *it, patName, strStudyUID, serieNumber);
 206          serieNumber+=2;
 207          //return 0;
 208       }
 209    }
 210 }
 211
 212
 213 void Load(std::ifstream &from, std::string imageName, std::string patName,
 214           std::string strStudyUID, int serieNumber)
 215 {
 216   std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
 217    if (!from)
 218       return;
 219   std::cout << " ========= Create Parametric images" << std::endl;
 220 /* was OK for single frames
 221 eg :
 222 ---------------------------
 223  Array dimensions = 58 x 56
 224 The following is the 2D array of peak circumferential strain map, all zero value
 225 pixels are outside the mask
 226      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 227 -----------------------------
 228    std::string str1;
 229    int nx, ny;
 230    // Get nx, ny
 231    from >> str1;
 232    from >> str1;
 233    from >> str1;
 234    from >> nx;
 235    from >> str1;
 236    from >> ny;
 237
 238    std::cout << "nx " << nx << " ny " << ny << std::endl;
 239
 240    // Skip 2 lines.
 241    std::getline(from, str1);
 242    std::cout << "[" << str1 << "]" << std::endl;
 243    std::getline(from, str1);
 244    std::cout << "[" << str1 << "]" << std::endl;
 245    std::getline(from, str1);
 246    std::cout << "[" << str1 << "]" << std::endl;
 247  */
 248
 249  /* now, we deal with multiframes
 250  eg :
 251  ------------------------------------------
 252 X dim, Y dim, N of frames = 52x59x14
 253 Temporal resolution = 31.9200 ms
 254 First frame starts at 47.9600 ms
 255 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
 256 All pixels with zero strain values are outside the masks.
 257      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 258 --------------------------------------------
 259 */
 260    std::string str1;
 261    int nx, ny, nf;
 262     from >> str1; // X dim,
 263     from >> str1;
 264     from >> str1; // Y dim,
 265     from >> str1;
 266     from >> str1; // N of frames =
 267     from >> str1;
 268     from >> str1;
 269     from >> str1;
 270     from >> str1; // 52x59x14
 271
 272     std::cout << "[" << str1 << "]" << std::endl;
 273
 274     sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
 275     std::cout << nx << " " << ny << " " << nf << std::endl;
 276
 277     std::getline(from, str1);
 278
 279     from >> str1; // Temporal
 280     from >> str1; // Resolution
 281     from >> str1; // =
 282
 283     from >> str1;
 284
 285     float temporalResolution;
 286     sscanf( str1.c_str(),"%f",&temporalResolution);
 287     std::cout << "temporal Resolution = " << temporalResolution << std::endl;
 288     std::getline(from, str1);
 289
 290     from >> str1; // First
 291     from >> str1; // frame
 292     from >> str1; // starts
 293     from >> str1; // at
 294
 295     from >> str1;
 296     float timeStart;
 297     sscanf( str1.c_str(),"%f",&timeStart);
 298     std::cout << "time Start = " << timeStart << std::endl;
 299     std::getline(from, str1);
 300
 301        // Skip 3 lines.
 302     for (int k=0; k<2; k++) //
 303     {
 304        std::getline(from, str1);
 305        std::cout << str1 << std::endl;
 306     }
 307
 308   //float *f = new float(nx*ny);
 309    float *f = (float *) malloc(nx*ny*nf*sizeof(float));
 310   // float mini = FLT_MAX, maxi = FLT_MIN;
 311    float val;
 312
 313    std::string strSerieUID;
 314    strSerieUID =  gdcm::Util::CreateUniqueUID();
 315    int imageNumber = 0;
 316    float currentTime;
 317    currentTime = timeStart;
 318   for (int l=0; l<nf; l++)  // Loop on the frames
 319   {
 320      //std::cout << "Frame nb " << l << std::endl;
 321      for( int j=0; j<ny; j++)
 322      {
 323       for (int i=0; i<nx; i++)
 324       {
 325          str1="";
 326          //eatwhite(from);
 327          char c;
 328          for (;;) //eatwhite(from);
 329          {
 330             if (!from.get(c))
 331             {
 332                std::cout << " !from.get(c) ";
 333                break;
 334             }
 335             if (!isspace(c) )
 336             {
 337                //std::cout << " !isspace(c) ";
 338                from.putback(c);
 339                break;
 340             }
 341          } //end eatwhite(from);
 342
 343         // trouble : when space is missing "-0.0990263-8.8778"
 344         // is not interpreted as TWO values  :-(
 345         // from >> str1;
 346
 347          int first = 1;
 348          char previous = 'z';
 349          for(;;)
 350          {
 351             from.get(c);
 352             if ( c == ' ')
 353                break;
 354             if ( first != 1 && c == '-' && previous != 'e')
 355             {
 356                from.putback(c);
 357                //std::cout << " One more gauffre in frame:" << std::dec << l
 358                //         << ", line : " << j << " element " << i << std::endl;
 359                break;
 360              }
 361
 362              first=0;
 363              previous = c;
 364              str1=str1+c;
 365          }
 366
 367          val = (float)atof(str1.c_str());
 368          //std::cout << "  " << val;
 369          *(f+ /*l*nx*ny + */j*nx+i) = val;
 370
 371         if(from.eof())
 372         {
 373             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 374                       << std::endl;
 375            break;
 376          }
 377       }
 378
 379       //std::cout << std::endl;
 380       //std::cout << std::endl << " line nb : " << j
 381       //          << " line length : " << l << std::endl;
 382     }
 383
 384    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 385 /*
 386 // values are expressed as % as a fraction, actually!)
 387 // It's useless to rescale them as uint16_t : just *100
 388     uint16_t *img = new uint16_t[ny*nx];
 389     uint16_t *ptr = img;
 390     float *tmp = f;
 391     float div = maxi-mini;
 392     std::cout << "div : " << div << " mini : " << mini << std::endl;
 393     for( int k=0; k<ny*nx; k++)
 394     {
 395        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 396        tmp ++;
 397        ptr++;
 398     }
 399 */
 400     int16_t *img = new int16_t[ny*nx];
 401     int16_t *ptr = img;
 402     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 403     for( int k=0; k<ny*nx; k++)
 404     {
 405        if(*tmp > 1.0) // artefacted pixel
 406           *ptr = 0;
 407        else        /// \todo FIXME : what about max threshold ?
 408         *ptr = (int16_t)(*tmp *100);
 409
 410       // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 411        tmp ++;
 412        ptr++;
 413     }
 414
 415  // gdcm::Debug::DebugOn();
 416
 417         std::ostringstream str;
 418         gdcm::File *file;
 419         file = gdcm::File::New();
 420
 421   // Set the image size
 422         str.str("");
 423         str << nx;
 424         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 425         str.str("");
 426         str << ny;
 427         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 428
 429   // Set the pixel type
 430   //      16; //8, 16, 32
 431         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 432         str.str("");
 433         str << 16; // may be 12 or 16 if componentSize =16
 434         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 435
 436         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 437
 438   // Set the pixel representation // 0/1 1 : signed
 439         file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 440
 441   // Set the samples per pixel // 1:Grey level, 3:RGB
 442         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 443
 444 /*
 445   // Set Rescale Intercept
 446         str.str("");
 447         str << div;
 448         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 449
 450   // Set Rescale Slope
 451         str.str("");
 452         str << mini;
 453         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 454 */
 455
 456 // 0020 0037 DS 6 Image Orientation (Patient)
 457          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 458
 459 // 0020 0032 DS 3 Image Position (Patient)
 460         char charImagePosition[256];
 461         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 462
 463 // 0020 0x1041 DS 1 Slice Location
 464         sprintf(charImagePosition,"%f",float(l%nf));
 465         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 466
 467 //0008 103e LO 1 Series Description
 468         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 469
 470         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 471         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 472         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 473
 474 //0020 0011 "IS" Series Number
 475          sprintf(charImagePosition,"%04d",serieNumber);
 476          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 477
 478 //0020 0011 "IS" Instance Number
 479          sprintf(charImagePosition,"%04d",imageNumber);
 480          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 481
 482 //0018 1060 "DS" Time Trigger
 483          sprintf(charImagePosition,"%f",currentTime);
 484          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 485
 486    // file->Print();
 487
 488     gdcm::FileHelper *fh;
 489     fh = gdcm::FileHelper::New(file);
 490     // cast is just to avoid warnings (*no* conversion)
 491     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 492     fh->SetWriteModeToRaw();
 493     fh->SetWriteTypeToDcmExplVR();
 494
 495     fh->SetWriteTypeToDcmExplVR();
 496
 497     char numero[10];
 498     sprintf(numero, "%02d", l);
 499     std::string fileName = imageName + "." + numero + ".dcm";
 500     std::cout << "fileName " << fileName << std::endl;
 501
 502     if( !fh->Write(fileName))
 503        std::cout << "Failed for [" << fileName << "]\n"
 504                  << "           File is unwrittable" << std::endl;
 505
 506     delete img;
 507     currentTime += temporalResolution;
 508     imageNumber ++;
 509
 510   } // end loop on frames
 511
 512
 513    // Anatomical Images.
 514   std::cout << " ========= Create Anatomical images" << std::endl;
 515
 516   strSerieUID =  gdcm::Util::CreateUniqueUID();
 517   imageNumber = 0;
 518   currentTime = timeStart;
 519
 520   for (int l=0; l<nf; l++)  // Loop on the frames
 521   {
 522    //std::cout << "Frame nb " << l << std::endl;
 523      for( int j=0; j<ny; j++)
 524      {
 525       int l =0;
 526       for (int i=0; i<nx; i++)
 527       {
 528          //eatwhite(from);
 529          char c;
 530          for (;;)
 531          {
 532             if (!from.get(c))
 533                break;
 534             if (!isspace(c))
 535             {
 536                from.putback(c);
 537                break;
 538             }
 539          }
 540          from >> str1;
 541          val = (float)atof(str1.c_str());
 542         // std::cout << "  " << val;
 543          *(f+ /*l*nx*ny + */j*nx+i) = val;
 544
 545         if(from.eof())
 546         {
 547             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 548                       << std::endl;
 549            break;
 550          }
 551          l++;
 552       }
 553      }
 554
 555     uint16_t *img = new uint16_t[ny*nx];
 556     uint16_t *ptr = img;
 557     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 558     for( int k=0; k<ny*nx; k++)
 559     {
 560        *ptr = (int16_t)(*tmp);
 561        tmp ++;
 562        ptr++;
 563     }
 564         std::ostringstream str;
 565         gdcm::File *file;
 566         file = gdcm::File::New();
 567
 568   // Set the image size
 569         str.str("");
 570         str << nx;
 571         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 572         str.str("");
 573         str << ny;
 574         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 575
 576   // Set the pixel type
 577   //      16; //8, 16, 32
 578         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 579         str.str("");
 580         str << 16; // may be 12 or 16 if componentSize =16
 581         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 582
 583         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 584
 585   // Set the pixel representation // 0/1 1 : signed
 586         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 587
 588   // Set the samples per pixel // 1:Grey level, 3:RGB
 589         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 590
 591 /*
 592   // Set Rescale Intercept
 593         str.str("");
 594         str << div;
 595         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 596
 597   // Set Rescale Slope
 598         str.str("");
 599         str << mini;
 600         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 601 */
 602
 603 // 0020 0037 DS 6 Image Orientation (Patient)
 604          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 605
 606 // 0020 0032 DS 3 Image Position (Patient)
 607         char charImagePosition[256];
 608         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 609
 610 // 0020 0x1041 DS 1 Slice Location
 611         sprintf(charImagePosition,"%f",float(l%nf));
 612         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 613
 614 //0008 103e LO 1 Series Description
 615         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 616
 617         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 618         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 619         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 620
 621 //0020 0011 "IS" Series Number
 622          sprintf(charImagePosition,"%04d",serieNumber+1);
 623          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 624
 625 //0020 0011 "IS" Instance Number
 626          sprintf(charImagePosition,"%04d",imageNumber);
 627          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 628
 629 //0018 1060 "DS" Time Trigger
 630          sprintf(charImagePosition,"%f",currentTime);
 631          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 632    // file->Print();
 633
 634     gdcm::FileHelper *fh;
 635     fh = gdcm::FileHelper::New(file);
 636     // cast is just to avoid warnings (*no* conversion)
 637     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 638     fh->SetWriteModeToRaw();
 639     fh->SetWriteTypeToDcmExplVR();
 640
 641     fh->SetWriteTypeToDcmExplVR();
 642
 643     char numero[10];
 644     sprintf(numero, "%02d", l);
 645     std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
 646     std::cout << "fileName " << fileName << std::endl;
 647
 648     if( !fh->Write(fileName))
 649        std::cout << "Failed for [" << fileName << "]\n"
 650                  << "           File is unwrittable" << std::endl;
 651
 652     delete img;
 653     currentTime += temporalResolution;
 654     imageNumber ++;
 655
 656   } // end loop on frames
 657
 658    from.close();
 659 }