Example/DenseMultiFramesToDicom.cxx

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