Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2007/09/20 11:15:54 $
   7   Version:   $Revision: 1.9 $
   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 void Load(std::ifstream &from, std::string imageName, std::string patName,
  40           std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose);
  41
  42 //std::ifstream& eatwhite(std::ifstream& is);
  43
  44 /*
  45 window center (level) and width, are defined in the DICOM
  46 standard very precisely as follows (see PS 3.3 C.11.2.1.2):
  47 >
  48 >"These Attributes are applied according to the following pseudo-code,
  49 >where :
  50 x is the input value,
  51 y is an output value with a range from ymin to ymax,
  52 c is Window Center (0028,1050)
  53 w is Window Width (0028,1051):
  54 >
  55 >           if      (x  <= c - 0.5 - (w-1)/2), then y = ymin
  56 >           else if (x > c - 0.5 + (w-1)/2), then y = ymax,
  57 >           else    y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin
  58
  59 */
  60 /*
  61 From:   David Clunie - view profile
  62 Date:   Thurs, Jun 1 2006 3:03 pm
  63 Email:  David Clunie <dclu...@dclunie.com>
  64 Groups: comp.protocols.dicom
  65
  66 The value of x is the output of the preceding step, the so-called
  67 "modality LUT", which may either be:
  68
  69 - identity (no rescale values or Modality LUT, or the SOP Class is
  70   PET and rescale values are ignored), in which case x is the stored
  71   pixel value in (7FE0,0010)
  72
  73 - Rescale Slope and Intercept (as typically used in CT), in which
  74   case x is the value obtained from applying the rescale values to
  75   the stored pixel value
  76
  77 - an actual LUT, in which case x is the value stored in the LUT
  78   corresponding to the LUT index value that is the stored pixel
  79   value
  80
  81 The ymin and ymax are intended to represent the output range; for
  82 example, if the hypothetical Presentation LUT step that follows
  83 the VOI LUT (window) stage is an identity operation, then the
  84 ymin and ymax represent P-Values, which might be the range of
  85 digital driving levels for your display (calibrated to the GSDF),
  86 in the 8-bit wide output case ranging from 0 to 255, for example.
  87
  88 The terms brightness and contrast are not used in radiology imaging
  89 -the window center and width are used instead.
  90 */
  91
  92 int main(int argc, char *argv[])
  93 {
  94    START_USAGE(usage)
  95    " \n DenseMultiFramesToDicom :                                           \n",
  96    " - explores recursively the given (single Patient, single Study) directory",
  97    "         - examines the '.txt' files                                      ",
  98    "         - Converts the files into 16 bits Dicom files,                   ",
  99    " WARNING : directory must contain ONLY .txt files                         ",
 100    " usage:                                                                   ",
 101    " DenseMultiFramesToDicom dirin=rootDirectoryName                          ",
 102    "              [studyUID = ] [patName = ]                                  ",
 103    "              [listonly] [verbose] [debug]                                ",
 104    "                                                                          ",
 105    " studyUID   : *aware* user wants to add the serie                         ",
 106    "                                             to an already existing study ",
 107    " verbose  : user wants to run the program in 'verbose mode'               ",
 108    " debug    : *developer*  wants to run the program in 'debug mode'         ",
 109    FINISH_USAGE
 110
 111    // ----- Initialize Arguments Manager ------
 112
 113    GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
 114
 115    if (argc == 1 || am->ArgMgrDefined("usage"))
 116    {
 117       am->ArgMgrUsage(usage); // Display 'usage'
 118       delete am;
 119       return 0;
 120    }
 121
 122    const char *dirNamein;
 123    dirNamein  = am->ArgMgrGetString("dirin",".");
 124
 125    if (am->ArgMgrDefined("debug"))
 126       GDCM_NAME_SPACE::Debug::DebugOn();
 127
 128    int verbose  = am->ArgMgrDefined("verbose");
 129    int listonly = am->ArgMgrDefined("listonly");
 130    std::string patName = am->ArgMgrGetString("patname", dirNamein);
 131
 132    bool userDefinedStudy = ( 0 != am->ArgMgrDefined("studyUID") );
 133
 134    const char *studyUID;
 135    if (userDefinedStudy)
 136       studyUID  = am->ArgMgrGetString("studyUID");
 137
 138    // not described *on purpose* in the Usage !
 139    bool userDefinedSerie = ( 0 != am->ArgMgrDefined("serieUID") );
 140
 141    const char *serieUID;
 142    if(userDefinedSerie)
 143       serieUID = am->ArgMgrGetString("serieUID");
 144
 145    // if unused Param we give up
 146    if ( am->ArgMgrPrintUnusedLabels() )
 147    {
 148       am->ArgMgrUsage(usage);
 149       delete am;
 150       return 0;
 151    }
 152    delete am;  // we don't need Argument Manager any longer
 153
 154    // ----- Begin Processing -----
 155
 156    if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
 157    {
 158       std::cout << "KO : [" << dirNamein << "] is not a Directory."
 159                 << std::endl;
 160       return 0;
 161    }
 162    else
 163    {
 164       if (verbose)
 165          std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
 166    }
 167
 168    std::string strDirNamein(dirNamein);
 169    GDCM_NAME_SPACE::DirList dirList(strDirNamein, true); // (recursively) the list of files
 170
 171    if (listonly)
 172    {
 173       std::cout << "------------List of found files ------------" << std::endl;
 174       dirList.Print();
 175       std::cout << std::endl;
 176       return 0;
 177     }
 178
 179    std::string filenameout;
 180
 181    std::string strStudyUID;
 182    std::string strSerieUID;
 183
 184    if (userDefinedStudy)
 185       strSerieUID =  studyUID;
 186    else
 187       strStudyUID =  GDCM_NAME_SPACE::Util::CreateUniqueUID();
 188
 189    if (userDefinedStudy)
 190      strSerieUID =  serieUID;
 191    else
 192       strStudyUID =  GDCM_NAME_SPACE::Util::CreateUniqueUID();
 193
 194
 195    int serieNumber =0;
 196    GDCM_NAME_SPACE::DirListType fileNames;
 197    fileNames = dirList.GetFilenames();
 198    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 199                                     it != fileNames.end();
 200                                   ++it)
 201    {
 202       if ( GDCM_NAME_SPACE::Util::GetName((*it)).c_str()[0] == '.' )
 203       {
 204       // skip hidden files
 205          continue;
 206       }
 207       int sz = (*it).size();
 208       if ( (*it).c_str()[sz-1] != 't')
 209       {
 210          // skip non .txt files
 211          continue;
 212       }
 213       std::ifstream from( (*it).c_str() );
 214       if ( !from )
 215       {
 216          std::cout << "Can't open file" << *it << std::endl;
 217          //return 0;
 218       }
 219       else
 220       {
 221          if (verbose)
 222            std::cout << "Success in open file" << *it << std::endl;
 223          Load(from, *it, patName, strStudyUID, strSerieUID, serieNumber, verbose);
 224          serieNumber+=2;
 225          //return 0;
 226       }
 227    }
 228    return 1;
 229 }
 230
 231
 232 void Load(std::ifstream &from, std::string imageName, std::string patName,
 233           std::string strStudyUID, std::string strSerieUID, int serieNumber, int verbose)
 234 {
 235    if (verbose)
 236       std::cout << " ========= Deal with file [" << imageName << "]" << std::endl;
 237    if (!from)
 238       return;
 239    if (verbose)
 240      std::cout << " ========= Create Parametric images" << std::endl;
 241 /* was OK for single frames
 242 eg :
 243 ---------------------------
 244  Array dimensions = 58 x 56
 245 The following is the 2D array of peak circumferential strain map, all zero value
 246 pixels are outside the mask
 247      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 248 -----------------------------
 249    std::string str1;
 250    int nx, ny;
 251    // Get nx, ny
 252    from >> str1;
 253    from >> str1;
 254    from >> str1;
 255    from >> nx;
 256    from >> str1;
 257    from >> ny;
 258
 259    std::cout << "nx " << nx << " ny " << ny << std::endl;
 260
 261    // Skip 2 lines.
 262    std::getline(from, str1);
 263    std::cout << "[" << str1 << "]" << std::endl;
 264    std::getline(from, str1);
 265    std::cout << "[" << str1 << "]" << std::endl;
 266    std::getline(from, str1);
 267    std::cout << "[" << str1 << "]" << std::endl;
 268  */
 269
 270  /* now, we deal with multiframes
 271  eg :
 272  ------------------------------------------
 273 X dim, Y dim, N of frames = 52x59x14
 274 Temporal resolution = 31.9200 ms
 275 First frame starts at 47.9600 ms
 276 The following are the 3D arrays of peak circumferential strain map and the magnitude images,
 277 All pixels with zero strain values are outside the masks.
 278      0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
 279 --------------------------------------------
 280 */
 281    std::string str1;
 282    int nx, ny, nf;
 283     from >> str1; // X dim,
 284     from >> str1;
 285     from >> str1; // Y dim,
 286     from >> str1;
 287     from >> str1; // N of frames =
 288     from >> str1;
 289     from >> str1;
 290     from >> str1;
 291     from >> str1; // 52x59x14
 292
 293     if(verbose)
 294        std::cout << "[" << str1 << "]" << std::endl;
 295
 296     sscanf( str1.c_str(),"%dx%dx%d", &nx,&ny,&nf);
 297     std::cout << nx << " " << ny << " " << nf << std::endl;
 298
 299     std::getline(from, str1);
 300
 301     from >> str1; // Temporal
 302     from >> str1; // Resolution
 303     from >> str1; // =
 304
 305     from >> str1;
 306
 307     float temporalResolution;
 308     sscanf( str1.c_str(),"%f",&temporalResolution);
 309     if(verbose)
 310       std::cout << "temporal Resolution = " << temporalResolution << std::endl;
 311     std::getline(from, str1);
 312
 313     from >> str1; // First
 314     from >> str1; // frame
 315     from >> str1; // starts
 316     from >> str1; // at
 317
 318     from >> str1;
 319     float timeStart;
 320     sscanf( str1.c_str(),"%f",&timeStart);
 321     std::cout << "time Start = " << timeStart << std::endl;
 322     std::getline(from, str1);
 323
 324        // Skip 3 lines.
 325     for (int k=0; k<2; k++) //
 326     {
 327        std::getline(from, str1);
 328        std::cout << str1 << std::endl;
 329     }
 330
 331   //float *f = new float(nx*ny);
 332   // -->float *f = new float[nx*ny]; // Would be better !
 333    float *f = (float *) malloc(nx*ny*nf*sizeof(float));
 334   // float mini = FLT_MAX, maxi = FLT_MIN;
 335    float val;
 336
 337    int imageNumber = 0;
 338    float currentTime;
 339    currentTime = timeStart;
 340    int l1;
 341    for (l1=0; l1<nf; l1++)  // Loop on the frames
 342    {
 343      //std::cout << "Frame nb " << l1 << std::endl;
 344      for( int j=0; j<ny; j++)
 345      {
 346       for (int i=0; i<nx; i++)
 347       {
 348          str1="";
 349          //eatwhite(from);
 350          char c;
 351          for (;;) //eatwhite(from);
 352          {
 353             if (!from.get(c))
 354             {
 355                std::cout << " !from.get(c) ";
 356                break;
 357             }
 358             if (!isspace(c) )
 359             {
 360                //std::cout << " !isspace(c) ";
 361                from.putback(c);
 362                break;
 363             }
 364          } //end eatwhite(from);
 365
 366         // trouble : when space is missing "-0.0990263-8.8778"
 367         // is not interpreted as TWO values  :-(
 368         // from >> str1;
 369
 370          int first = 1;
 371          char previous = 'z';
 372          for(;;)
 373          {
 374             from.get(c);
 375             if ( c == ' ')
 376                break;
 377             if ( first != 1 && c == '-' && previous != 'e')
 378             {
 379                from.putback(c);
 380                //std::cout << " One more gauffre in frame:" << std::dec << l
 381                //         << ", line : " << j << " element " << i << std::endl;
 382                break;
 383              }
 384
 385              first=0;
 386              previous = c;
 387              str1=str1+c;
 388          }
 389
 390          val = (float)atof(str1.c_str());
 391          //std::cout << "  " << val;
 392          *(f+ /*l*nx*ny + */j*nx+i) = val;
 393
 394         if(from.eof())
 395         {
 396             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 397                       << std::endl;
 398            break;
 399          }
 400       }
 401
 402       //std::cout << std::endl;
 403       //std::cout << std::endl << " line nb : " << j
 404       //          << " line length : " << l << std::endl;
 405     }
 406
 407    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 408 /*
 409 // values are expressed as % as a fraction, actually!)
 410 // It's useless to rescale them as uint16_t : just *100
 411     uint16_t *img = new uint16_t[ny*nx];
 412     uint16_t *ptr = img;
 413     float *tmp = f;
 414     float div = maxi-mini;
 415     std::cout << "div : " << div << " mini : " << mini << std::endl;
 416     for( int k=0; k<ny*nx; k++)
 417     {
 418        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 419        tmp ++;
 420        ptr++;
 421     }
 422 */
 423     int16_t *img = new int16_t[ny*nx];
 424     int16_t *ptr = img;
 425     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 426     for( int k=0; k<ny*nx; k++)
 427     {
 428        if(*tmp > 1.0) // artefacted pixel
 429           *ptr = 0;
 430        else        /// \todo FIXME : what about max threshold ?
 431         *ptr = (int16_t)(*tmp *100);
 432
 433       // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 434        tmp ++;
 435        ptr++;
 436     }
 437
 438  // GDCM_NAME_SPACE::Debug::DebugOn();
 439
 440         std::ostringstream str;
 441         GDCM_NAME_SPACE::File *file;
 442         file = GDCM_NAME_SPACE::File::New();
 443
 444   // Set the image size
 445         str.str("");
 446         str << nx;
 447         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 448         str.str("");
 449         str << ny;
 450         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 451
 452   // Set the pixel type
 453   //      16; //8, 16, 32
 454         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 455         str.str("");
 456         str << 16; // may be 12 or 16 if componentSize =16
 457         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 458
 459         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 460
 461   // Set the pixel representation // 0/1 1 : signed
 462         file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 463
 464   // Set the samples per pixel // 1:Grey level, 3:RGB
 465         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 466
 467 /*
 468   // Set Rescale Intercept
 469         str.str("");
 470         str << div;
 471         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 472
 473   // Set Rescale Slope
 474         str.str("");
 475         str << mini;
 476         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 477 */
 478
 479 // 0020 0037 DS 6 Image Orientation (Patient)
 480          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 481
 482 // 0020 0032 DS 3 Image Position (Patient)
 483         char charImagePosition[256];
 484         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
 485
 486 // 0020 0x1041 DS 1 Slice Location
 487         sprintf(charImagePosition,"%f",float(l1%nf));
 488         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 489
 490 //0008 103e LO 1 Series Description
 491         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 492
 493         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 494         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 495         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 496
 497 //0020 0011 "IS" Series Number
 498          sprintf(charImagePosition,"%04d",serieNumber);
 499          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 500
 501 //0020 0011 "IS" Instance Number
 502          sprintf(charImagePosition,"%04d",imageNumber);
 503          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 504
 505 //0018 1060 "DS" Time Trigger
 506          sprintf(charImagePosition,"%f",currentTime);
 507          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 508
 509    // file->Print();
 510
 511     GDCM_NAME_SPACE::FileHelper *fh;
 512     fh = GDCM_NAME_SPACE::FileHelper::New(file);
 513     // cast is just to avoid warnings (*no* conversion)
 514     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 515     fh->SetWriteModeToRaw();
 516     fh->SetWriteTypeToDcmExplVR();
 517
 518     fh->SetWriteTypeToDcmExplVR();
 519
 520     char numero[10];
 521     sprintf(numero, "%02d", l1);
 522     std::string fileName = imageName + "." + numero + ".dcm";
 523     if(verbose)
 524       std::cout << "fileName " << fileName << std::endl;
 525
 526     if( !fh->Write(fileName))
 527        std::cout << "Failed for [" << fileName << "]\n"
 528                  << "           File is unwrittable" << std::endl;
 529
 530     delete img;
 531     currentTime += temporalResolution;
 532     imageNumber ++;
 533
 534   } // end loop on frames
 535
 536
 537    // Anatomical Images.
 538   std::cout << " ========= Create Anatomical images" << std::endl;
 539
 540   strSerieUID =  GDCM_NAME_SPACE::Util::CreateUniqueUID();
 541   imageNumber = 0;
 542   currentTime = timeStart;
 543
 544   for (int fr=0; fr<nf; fr++)  // Loop on the frames
 545   {
 546    //std::cout << "Frame nb " << l << std::endl;
 547      for( int j=0; j<ny; j++)
 548      {
 549       int l =0;
 550       for (int i=0; i<nx; i++)
 551       {
 552          //eatwhite(from);
 553          char c;
 554          for (;;)
 555          {
 556             if (!from.get(c))
 557                break;
 558             if (!isspace(c))
 559             {
 560                from.putback(c);
 561                break;
 562             }
 563          }
 564          from >> str1;
 565          val = (float)atof(str1.c_str());
 566         // std::cout << "  " << val;
 567          *(f+ /*l*nx*ny + */j*nx+i) = val;
 568
 569         if(from.eof())
 570         {
 571             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 572                       << std::endl;
 573            break;
 574          }
 575          l++;
 576       }
 577      }
 578
 579     uint16_t *img = new uint16_t[ny*nx];
 580     uint16_t *ptr = img;
 581     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 582     for( int k=0; k<ny*nx; k++)
 583     {
 584        *ptr = (int16_t)(*tmp);
 585        tmp ++;
 586        ptr++;
 587     }
 588         std::ostringstream str;
 589         GDCM_NAME_SPACE::File *file;
 590         file = GDCM_NAME_SPACE::File::New();
 591
 592   // Set the image size
 593         str.str("");
 594         str << nx;
 595         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 596         str.str("");
 597         str << ny;
 598         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 599
 600   // Set the pixel type
 601   //      16; //8, 16, 32
 602         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 603         str.str("");
 604         str << 16; // may be 12 or 16 if componentSize =16
 605         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 606
 607         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 608
 609   // Set the pixel representation // 0/1 1 : signed
 610         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 611
 612   // Set the samples per pixel // 1:Grey level, 3:RGB
 613         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 614
 615 /*
 616   // Set Rescale Intercept
 617         str.str("");
 618         str << div;
 619         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 620
 621   // Set Rescale Slope
 622         str.str("");
 623         str << mini;
 624         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 625 */
 626
 627 // 0020 0037 DS 6 Image Orientation (Patient)
 628          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 629
 630 // 0020 0032 DS 3 Image Position (Patient)
 631         char charImagePosition[256];
 632         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l1%nf));
 633
 634 // 0020 0x1041 DS 1 Slice Location
 635         sprintf(charImagePosition,"%f",float(l1%nf));
 636         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 637
 638 //0008 103e LO 1 Series Description
 639         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 640
 641         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 642         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 643         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 644
 645 //0020 0011 "IS" Series Number
 646          sprintf(charImagePosition,"%04d",serieNumber+1);
 647          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 648
 649 //0020 0011 "IS" Instance Number
 650          sprintf(charImagePosition,"%04d",imageNumber);
 651          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 652
 653 //0018 1060 "DS" Time Trigger
 654          sprintf(charImagePosition,"%f",currentTime);
 655          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 656    // file->Print();
 657
 658     GDCM_NAME_SPACE::FileHelper *fh;
 659     fh = GDCM_NAME_SPACE::FileHelper::New(file);
 660     // cast is just to avoid warnings (*no* conversion)
 661     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 662     fh->SetWriteModeToRaw();
 663     fh->SetWriteTypeToDcmExplVR();
 664
 665     fh->SetWriteTypeToDcmExplVR();
 666
 667     char numero[10];
 668     sprintf(numero, "%02d", l1);
 669     std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
 670     std::cout << "fileName " << fileName << std::endl;
 671
 672     if( !fh->Write(fileName))
 673        std::cout << "Failed for [" << fileName << "]\n"
 674                  << "           File is unwrittable" << std::endl;
 675
 676     delete img;
 677     currentTime += temporalResolution;
 678     imageNumber ++;
 679
 680   } // end loop on frames
 681
 682    from.close();
 683 } // end void Load(