Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2006/07/21 17:18:27 $
   7   Version:   $Revision: 1.2 $
   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<4; k++)
 267     {
 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      for( int j=0; j<ny; j++)
 284      {
 285       int l =0;
 286       for (int i=0; i<nx; i++)
 287       {
 288          //eatwhite(from);
 289          char c;
 290          for (;;)
 291          {
 292             if (!from.get(c))
 293             {
 294                std::cout << " !from.get(c) ";
 295                break;
 296             }
 297             if (!isspace(c))
 298             {
 299                //std::cout << " !isspace(c) ";
 300                from.putback(c);
 301                break;
 302             }
 303          }
 304          from >> str1;
 305          val = (float)atof(str1.c_str());
 306         std::cout << "  " << val;
 307          *(f+ /*l*nx*ny + */j*nx+i) = val;
 308
 309         if(from.eof())
 310         {
 311             std::cout << "Missing values at [" << j <<"," << i << "]"
 312                       << std::endl;
 313            break;
 314          }
 315          l++;
 316       }
 317
 318       std::cout << std::endl;
 319       //std::cout << std::endl << " line nb : " << j
 320       //          << " line length : " << l << std::endl;
 321
 322     }
 323
 324    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 325 /*
 326 // values are expressed as % as a fraction, actually!)
 327 // It's useless to rescale them as uint16_t : just *100
 328     uint16_t *img = new uint16_t[ny*nx];
 329     uint16_t *ptr = img;
 330     float *tmp = f;
 331     float div = maxi-mini;
 332     std::cout << "div : " << div << " mini : " << mini << std::endl;
 333     for( int k=0; k<ny*nx; k++)
 334     {
 335        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 336        tmp ++;
 337        ptr++;
 338     }
 339 */
 340     int16_t *img = new int16_t[ny*nx];
 341     int16_t *ptr = img;
 342     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 343     for( int k=0; k<ny*nx; k++)
 344     {
 345        if(*tmp > 1.0) // artefacted pixel
 346           *ptr = 0;
 347        else        /// \todo FIXME : what about max threshold ?
 348         *ptr = (int16_t)(*tmp *100);
 349
 350       // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 351        tmp ++;
 352        ptr++;
 353     }
 354
 355  // gdcm::Debug::DebugOn();
 356
 357         std::ostringstream str;
 358         gdcm::File *file;
 359         file = gdcm::File::New();
 360
 361   // Set the image size
 362         str.str("");
 363         str << nx;
 364         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 365         str.str("");
 366         str << ny;
 367         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 368
 369   // Set the pixel type
 370   //      16; //8, 16, 32
 371         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 372         str.str("");
 373         str << 16; // may be 12 or 16 if componentSize =16
 374         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 375
 376         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 377
 378   // Set the pixel representation // 0/1 1 : signed
 379         file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 380
 381   // Set the samples per pixel // 1:Grey level, 3:RGB
 382         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 383
 384 /*
 385   // Set Rescale Intercept
 386         str.str("");
 387         str << div;
 388         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 389
 390   // Set Rescale Slope
 391         str.str("");
 392         str << mini;
 393         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 394 */
 395
 396 // 0020 0037 DS 6 Image Orientation (Patient)
 397          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 398
 399 // 0020 0032 DS 3 Image Position (Patient)
 400         char charImagePosition[256];
 401         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 402
 403 // 0020 0x1041 DS 1 Slice Location
 404         sprintf(charImagePosition,"%f",float(l%nf));
 405         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 406
 407 //0008 103e LO 1 Series Description
 408         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 409
 410         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 411         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 412         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 413
 414 //0020 0011 "IS" Series Number
 415          sprintf(charImagePosition,"%04d",serieNumber);
 416          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 417
 418    // file->Print();
 419
 420     gdcm::FileHelper *fh;
 421     fh = gdcm::FileHelper::New(file);
 422     // cast is just to avoid warnings (*no* conversion)
 423     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 424     fh->SetWriteModeToRaw();
 425     fh->SetWriteTypeToDcmExplVR();
 426
 427     fh->SetWriteTypeToDcmExplVR();
 428
 429     char numero[10];
 430     sprintf(numero, "%02d", l);
 431     std::string fileName = imageName + "." + numero + ".dcm";
 432     std::cout << "fileName " << fileName << std::endl;
 433
 434     if( !fh->Write(fileName))
 435        std::cout << "Failed for [" << fileName << "]\n"
 436                  << "           File is unwrittable" << std::endl;
 437
 438     delete img;
 439   } // end loop on frames
 440
 441     int pos = (long)(from.tellg());
 442  std::cout << std::hex << "==============================================="
 443            << pos<< std::endl;
 444
 445
 446    // Anatomical Images.
 447 std::cout << "Anatomical images" << std::endl;
 448
 449    strSerieUID =  gdcm::Util::CreateUniqueUID();
 450
 451   for (int l=0; l<nf; l++)  // Loop on the frames
 452   {
 453 std::cout << "Frame nb " << l << std::endl;
 454      for( int j=0; j<ny; j++)
 455      {
 456       int l =0;
 457       for (int i=0; i<nx; i++)
 458       {
 459          //eatwhite(from);
 460          char c;
 461          for (;;)
 462          {
 463             if (!from.get(c))
 464                break;
 465             if (!isspace(c))
 466             {
 467                from.putback(c);
 468                break;
 469             }
 470          }
 471          from >> str1;
 472          val = (float)atof(str1.c_str());
 473         // std::cout << "  " << val;
 474          *(f+ /*l*nx*ny + */j*nx+i) = val;
 475
 476         if(from.eof())
 477         {
 478             std::cout << "Missing values at [" << j <<"," << i << "]"
 479                       << std::endl;
 480            break;
 481          }
 482          l++;
 483       }
 484      }
 485
 486     uint16_t *img = new uint16_t[ny*nx];
 487     uint16_t *ptr = img;
 488     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 489     for( int k=0; k<ny*nx; k++)
 490     {
 491        *ptr = (int16_t)(*tmp *100);
 492        tmp ++;
 493        ptr++;
 494     }
 495
 496         std::ostringstream str;
 497         gdcm::File *file;
 498         file = gdcm::File::New();
 499
 500   // Set the image size
 501         str.str("");
 502         str << nx;
 503         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 504         str.str("");
 505         str << ny;
 506         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 507
 508   // Set the pixel type
 509   //      16; //8, 16, 32
 510         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 511         str.str("");
 512         str << 16; // may be 12 or 16 if componentSize =16
 513         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 514
 515         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 516
 517   // Set the pixel representation // 0/1 1 : signed
 518         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 519
 520   // Set the samples per pixel // 1:Grey level, 3:RGB
 521         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 522
 523 /*
 524   // Set Rescale Intercept
 525         str.str("");
 526         str << div;
 527         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 528
 529   // Set Rescale Slope
 530         str.str("");
 531         str << mini;
 532         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 533 */
 534
 535 // 0020 0037 DS 6 Image Orientation (Patient)
 536          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 537
 538 // 0020 0032 DS 3 Image Position (Patient)
 539         char charImagePosition[256];
 540         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 541
 542 // 0020 0x1041 DS 1 Slice Location
 543         sprintf(charImagePosition,"%f",float(l%nf));
 544         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 545
 546 //0008 103e LO 1 Series Description
 547         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 548
 549         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 550         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 551         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 552
 553 //0020 0011 "IS" Series Number
 554          sprintf(charImagePosition,"%04d",serieNumber+1);
 555          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 556
 557    // file->Print();
 558
 559     gdcm::FileHelper *fh;
 560     fh = gdcm::FileHelper::New(file);
 561     // cast is just to avoid warnings (*no* conversion)
 562     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 563     fh->SetWriteModeToRaw();
 564     fh->SetWriteTypeToDcmExplVR();
 565
 566     fh->SetWriteTypeToDcmExplVR();
 567
 568     char numero[10];
 569     sprintf(numero, "%02d", l);
 570     std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
 571     std::cout << "fileName " << fileName << std::endl;
 572
 573     if( !fh->Write(fileName))
 574        std::cout << "Failed for [" << fileName << "]\n"
 575                  << "           File is unwrittable" << std::endl;
 576
 577     delete img;
 578
 579
 580
 581   } // end loop on frames
 582
 583    from.close();
 584 }