Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2006/07/20 17:15:28 $
   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 #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
 186
 187       std::ifstream from( (*it).c_str() );
 188       if ( !from )
 189       {
 190          std::cout << "Can't open file" << *it << std::endl;
 191          //return 0;
 192       }
 193       else
 194       {
 195          std::cout << "Success in open file" << *it << std::endl;
 196          Load(from, *it, patName, strStudyUID, serieNumber);
 197          serieNumber++;
 198          //return 0;
 199       }
 200    }
 201 }
 202
 203
 204 void Load(std::ifstream &from, std::string imageName, std::string patName,
 205           std::string strStudyUID, int serieNumber)
 206 {
 207    if (!from)
 208       return;
 209
 210 /* was OK for single frames
 211 eg :
 212 ---------------------------
 213  Array dimensions = 58 x 56
 214 The following is the 2D array of peak circumferential strain map, all zero value
 215 pixels are outside the mask
 216      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 217 -----------------------------
 218    std::string str1;
 219    int nx, ny;
 220    // Get nx, ny
 221    from >> str1;
 222    from >> str1;
 223    from >> str1;
 224    from >> nx;
 225    from >> str1;
 226    from >> ny;
 227
 228    std::cout << "nx " << nx << " ny " << ny << std::endl;
 229
 230    // Skip 2 lines.
 231    std::getline(from, str1);
 232    std::cout << "[" << str1 << "]" << std::endl;
 233    std::getline(from, str1);
 234    std::cout << "[" << str1 << "]" << std::endl;
 235    std::getline(from, str1);
 236    std::cout << "[" << str1 << "]" << std::endl;
 237  */
 238
 239  /* now, we deal with multiframes
 240  eg :
 241  ------------------------------------------
 242 X dim, Y dim, N of frames = 52x59x14
 243 Temporal resolution = 31.9200 ms
 244 First frame starts at 47.9600 ms
 245 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
 246 All pixels with zero strain values are outside the masks.
 247      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 248 --------------------------------------------
 249 */
 250    std::string str1;
 251    int nx, ny, nf;
 252     from >> str1; // X dim,
 253     from >> str1;
 254     from >> str1; // Y dim,
 255     from >> str1;
 256     from >> str1; // N of frames =
 257     from >> str1;
 258     from >> str1;
 259     from >> str1;
 260     from >> str1; // 52x59x14
 261
 262     std::cout << "[" << str1 << "]" << std::endl;
 263
 264     sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
 265     std::cout << nx << " " << ny << " " << nf << std::endl;
 266
 267        // Skip 4 lines.
 268     for (int k=0; k<4; k++)
 269     {
 270        std::getline(from, str1);
 271        std::cout << str1 << std::endl;
 272     }
 273
 274    //float *f = new float(nx*ny);
 275    float *f = (float *) malloc(nx*ny*nf*sizeof(float));
 276   // float mini = FLT_MAX, maxi = FLT_MIN;
 277    float val;
 278
 279
 280    std::string strSerieUID;
 281    strSerieUID =  gdcm::Util::CreateUniqueUID();
 282
 283   for (int l=0; l<nf; l++)  // Loop on the frames
 284   {
 285      for( int j=0; j<ny; j++)
 286      {
 287       int l =0;
 288       for (int i=0; i<nx; i++)
 289       {
 290          //eatwhite(from);
 291          char c;
 292          for (;;)
 293          {
 294             if (!from.get(c))
 295                break;
 296             if (!isspace(c))
 297             {
 298                from.putback(c);
 299                break;
 300             }
 301          }
 302          from >> str1;
 303          val = (float)atof(str1.c_str());
 304         // std::cout << "  " << val;
 305          *(f+ /*l*nx*ny + */j*nx+i) = val;
 306
 307         if(from.eof())
 308         {
 309             std::cout << "Missing values at [" << j <<"," << i << "]"
 310                       << std::endl;
 311            break;
 312          }
 313          l++;
 314       }
 315
 316        std::cout << std::endl;
 317       //std::cout << std::endl << " line nb : " << j
 318       //          << " line length : " << l << std::endl;
 319
 320     }
 321
 322
 323    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 324 /*
 325 // values are expressed as % as a fraction, actually!)
 326 // It's useless to rescale them as uint16_t : just *100
 327     uint16_t *img = new uint16_t[ny*nx];
 328     uint16_t *ptr = img;
 329     float *tmp = f;
 330     float div = maxi-mini;
 331     std::cout << "div : " << div << " mini : " << mini << std::endl;
 332     for( int k=0; k<ny*nx; k++)
 333     {
 334        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 335        tmp ++;
 336        ptr++;
 337     }
 338 */
 339     uint16_t *img = new uint16_t[ny*nx];
 340     uint16_t *ptr = img;
 341     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 342     for( int k=0; k<ny*nx; k++)
 343     {
 344        if(*tmp < 0) // artefacted pixel
 345           *ptr = 0;
 346        else        /// \todo FIXME : what about max threshold ?
 347         *ptr = (int16_t)(*tmp *100);
 348
 349        //std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 350        tmp ++;
 351        ptr++;
 352     }
 353
 354  // gdcm::Debug::DebugOn();
 355
 356         std::ostringstream str;
 357         gdcm::File *file;
 358         file = gdcm::File::New();
 359
 360   // Set the image size
 361         str.str("");
 362         str << nx;
 363         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 364         str.str("");
 365         str << ny;
 366         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 367
 368   // Set the pixel type
 369   //      16; //8, 16, 32
 370         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 371         str.str("");
 372         str << 16; // may be 12 or 16 if componentSize =16
 373         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 374
 375         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 376
 377   // Set the pixel representation // 0/1 1 : signed
 378         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 379
 380   // Set the samples per pixel // 1:Grey level, 3:RGB
 381         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 382
 383 /*
 384   // Set Rescale Intercept
 385         str.str("");
 386         str << div;
 387         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 388
 389   // Set Rescale Slope
 390         str.str("");
 391         str << mini;
 392         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 393 */
 394
 395 // 0020 0037 DS 6 Image Orientation (Patient)
 396          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 397
 398 // 0020 0032 DS 3 Image Position (Patient)
 399         char charImagePosition[256];
 400         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 401
 402 // 0020 0x1041 DS 1 Slice Location
 403         sprintf(charImagePosition,"%f",float(l%nf));
 404         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 405
 406 //0008 103e LO 1 Series Description
 407         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 408
 409         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 410         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 411         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 412
 413 //0020 0011 "IS" Series Number
 414          sprintf(charImagePosition,"%04d",serieNumber);
 415          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 416
 417    // file->Print();
 418
 419     gdcm::FileHelper *fh;
 420     fh = gdcm::FileHelper::New(file);
 421     // cast is just to avoid warnings (*no* conversion)
 422     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 423     fh->SetWriteModeToRaw();
 424     fh->SetWriteTypeToDcmExplVR();
 425
 426     fh->SetWriteTypeToDcmExplVR();
 427
 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    from.close();
 442 }