Example/DenseMultiFramesToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseMultiFramesToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2007/06/26 15:42:14 $
   7   Version:   $Revision: 1.6 $
   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 DenseMultiFramesToDicom :                                           \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_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::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_NAME_SPACE::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_NAME_SPACE::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_NAME_SPACE::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_NAME_SPACE::Util::CreateUniqueUID();
 178    int serieNumber =0;
 179    GDCM_NAME_SPACE::DirListType fileNames;
 180    fileNames = dirList.GetFilenames();
 181    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 182                                     it != fileNames.end();
 183                                   ++it)
 184    {
 185       if ( GDCM_NAME_SPACE::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 = new float[nx*ny]; // Would be better !
 310    float *f = (float *) malloc(nx*ny*nf*sizeof(float));
 311   // float mini = FLT_MAX, maxi = FLT_MIN;
 312    float val;
 313
 314    std::string strSerieUID;
 315    strSerieUID =  GDCM_NAME_SPACE::Util::CreateUniqueUID();
 316    int imageNumber = 0;
 317    float currentTime;
 318    currentTime = timeStart;
 319   for (int l=0; l<nf; l++)  // Loop on the frames
 320   {
 321      //std::cout << "Frame nb " << l << std::endl;
 322      for( int j=0; j<ny; j++)
 323      {
 324       for (int i=0; i<nx; i++)
 325       {
 326          str1="";
 327          //eatwhite(from);
 328          char c;
 329          for (;;) //eatwhite(from);
 330          {
 331             if (!from.get(c))
 332             {
 333                std::cout << " !from.get(c) ";
 334                break;
 335             }
 336             if (!isspace(c) )
 337             {
 338                //std::cout << " !isspace(c) ";
 339                from.putback(c);
 340                break;
 341             }
 342          } //end eatwhite(from);
 343
 344         // trouble : when space is missing "-0.0990263-8.8778"
 345         // is not interpreted as TWO values  :-(
 346         // from >> str1;
 347
 348          int first = 1;
 349          char previous = 'z';
 350          for(;;)
 351          {
 352             from.get(c);
 353             if ( c == ' ')
 354                break;
 355             if ( first != 1 && c == '-' && previous != 'e')
 356             {
 357                from.putback(c);
 358                //std::cout << " One more gauffre in frame:" << std::dec << l
 359                //         << ", line : " << j << " element " << i << std::endl;
 360                break;
 361              }
 362
 363              first=0;
 364              previous = c;
 365              str1=str1+c;
 366          }
 367
 368          val = (float)atof(str1.c_str());
 369          //std::cout << "  " << val;
 370          *(f+ /*l*nx*ny + */j*nx+i) = val;
 371
 372         if(from.eof())
 373         {
 374             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 375                       << std::endl;
 376            break;
 377          }
 378       }
 379
 380       //std::cout << std::endl;
 381       //std::cout << std::endl << " line nb : " << j
 382       //          << " line length : " << l << std::endl;
 383     }
 384
 385    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 386 /*
 387 // values are expressed as % as a fraction, actually!)
 388 // It's useless to rescale them as uint16_t : just *100
 389     uint16_t *img = new uint16_t[ny*nx];
 390     uint16_t *ptr = img;
 391     float *tmp = f;
 392     float div = maxi-mini;
 393     std::cout << "div : " << div << " mini : " << mini << std::endl;
 394     for( int k=0; k<ny*nx; k++)
 395     {
 396        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 397        tmp ++;
 398        ptr++;
 399     }
 400 */
 401     int16_t *img = new int16_t[ny*nx];
 402     int16_t *ptr = img;
 403     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 404     for( int k=0; k<ny*nx; k++)
 405     {
 406        if(*tmp > 1.0) // artefacted pixel
 407           *ptr = 0;
 408        else        /// \todo FIXME : what about max threshold ?
 409         *ptr = (int16_t)(*tmp *100);
 410
 411       // std::cout << std::dec << "[" << *tmp <<" : " << *ptr << "] ";
 412        tmp ++;
 413        ptr++;
 414     }
 415
 416  // GDCM_NAME_SPACE::Debug::DebugOn();
 417
 418         std::ostringstream str;
 419         GDCM_NAME_SPACE::File *file;
 420         file = GDCM_NAME_SPACE::File::New();
 421
 422   // Set the image size
 423         str.str("");
 424         str << nx;
 425         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 426         str.str("");
 427         str << ny;
 428         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 429
 430   // Set the pixel type
 431   //      16; //8, 16, 32
 432         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 433         str.str("");
 434         str << 16; // may be 12 or 16 if componentSize =16
 435         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 436
 437         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 438
 439   // Set the pixel representation // 0/1 1 : signed
 440         file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 441
 442   // Set the samples per pixel // 1:Grey level, 3:RGB
 443         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 444
 445 /*
 446   // Set Rescale Intercept
 447         str.str("");
 448         str << div;
 449         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 450
 451   // Set Rescale Slope
 452         str.str("");
 453         str << mini;
 454         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 455 */
 456
 457 // 0020 0037 DS 6 Image Orientation (Patient)
 458          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 459
 460 // 0020 0032 DS 3 Image Position (Patient)
 461         char charImagePosition[256];
 462         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 463
 464 // 0020 0x1041 DS 1 Slice Location
 465         sprintf(charImagePosition,"%f",float(l%nf));
 466         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 467
 468 //0008 103e LO 1 Series Description
 469         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 470
 471         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 472         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 473         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 474
 475 //0020 0011 "IS" Series Number
 476          sprintf(charImagePosition,"%04d",serieNumber);
 477          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 478
 479 //0020 0011 "IS" Instance Number
 480          sprintf(charImagePosition,"%04d",imageNumber);
 481          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 482
 483 //0018 1060 "DS" Time Trigger
 484          sprintf(charImagePosition,"%f",currentTime);
 485          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 486
 487    // file->Print();
 488
 489     GDCM_NAME_SPACE::FileHelper *fh;
 490     fh = GDCM_NAME_SPACE::FileHelper::New(file);
 491     // cast is just to avoid warnings (*no* conversion)
 492     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 493     fh->SetWriteModeToRaw();
 494     fh->SetWriteTypeToDcmExplVR();
 495
 496     fh->SetWriteTypeToDcmExplVR();
 497
 498     char numero[10];
 499     sprintf(numero, "%02d", l);
 500     std::string fileName = imageName + "." + numero + ".dcm";
 501     std::cout << "fileName " << fileName << std::endl;
 502
 503     if( !fh->Write(fileName))
 504        std::cout << "Failed for [" << fileName << "]\n"
 505                  << "           File is unwrittable" << std::endl;
 506
 507     delete img;
 508     currentTime += temporalResolution;
 509     imageNumber ++;
 510
 511   } // end loop on frames
 512
 513
 514    // Anatomical Images.
 515   std::cout << " ========= Create Anatomical images" << std::endl;
 516
 517   strSerieUID =  GDCM_NAME_SPACE::Util::CreateUniqueUID();
 518   imageNumber = 0;
 519   currentTime = timeStart;
 520
 521   for (int l=0; l<nf; l++)  // Loop on the frames
 522   {
 523    //std::cout << "Frame nb " << l << std::endl;
 524      for( int j=0; j<ny; j++)
 525      {
 526       int l =0;
 527       for (int i=0; i<nx; i++)
 528       {
 529          //eatwhite(from);
 530          char c;
 531          for (;;)
 532          {
 533             if (!from.get(c))
 534                break;
 535             if (!isspace(c))
 536             {
 537                from.putback(c);
 538                break;
 539             }
 540          }
 541          from >> str1;
 542          val = (float)atof(str1.c_str());
 543         // std::cout << "  " << val;
 544          *(f+ /*l*nx*ny + */j*nx+i) = val;
 545
 546         if(from.eof())
 547         {
 548             std::cout << "Missing values at [" << std::dec << j <<"," << i << "]"
 549                       << std::endl;
 550            break;
 551          }
 552          l++;
 553       }
 554      }
 555
 556     uint16_t *img = new uint16_t[ny*nx];
 557     uint16_t *ptr = img;
 558     float *tmp = f /* + l*ny*nx */ ; // start on image nb l.
 559     for( int k=0; k<ny*nx; k++)
 560     {
 561        *ptr = (int16_t)(*tmp);
 562        tmp ++;
 563        ptr++;
 564     }
 565         std::ostringstream str;
 566         GDCM_NAME_SPACE::File *file;
 567         file = GDCM_NAME_SPACE::File::New();
 568
 569   // Set the image size
 570         str.str("");
 571         str << nx;
 572         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 573         str.str("");
 574         str << ny;
 575         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 576
 577   // Set the pixel type
 578   //      16; //8, 16, 32
 579         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 580         str.str("");
 581         str << 16; // may be 12 or 16 if componentSize =16
 582         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 583
 584         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 585
 586   // Set the pixel representation // 0/1 1 : signed
 587         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 588
 589   // Set the samples per pixel // 1:Grey level, 3:RGB
 590         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 591
 592 /*
 593   // Set Rescale Intercept
 594         str.str("");
 595         str << div;
 596         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 597
 598   // Set Rescale Slope
 599         str.str("");
 600         str << mini;
 601         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 602 */
 603
 604 // 0020 0037 DS 6 Image Orientation (Patient)
 605          file->InsertEntryString("1.0\\0.0\\0.0\\0.0\\1.0\\0.0",0x0020,0x0037, "DS"); //[1\0\0\0\1\0] : Axial
 606
 607 // 0020 0032 DS 3 Image Position (Patient)
 608         char charImagePosition[256];
 609         sprintf(charImagePosition,"0.0\\0.0\\%f",float(l%nf));
 610
 611 // 0020 0x1041 DS 1 Slice Location
 612         sprintf(charImagePosition,"%f",float(l%nf));
 613         file->InsertEntryString(charImagePosition,0x0020,0x1041, "DS");
 614
 615 //0008 103e LO 1 Series Description
 616         file->InsertEntryString(imageName,0x0008,0x103e, "LO");
 617
 618         file->InsertEntryString(strStudyUID,0x0020,0x000d,"UI");
 619         file->InsertEntryString(strSerieUID,0x0020,0x000e,"UI");
 620         file->InsertEntryString(patName,0x0010,0x0010, "PN");   // Patient's Name
 621
 622 //0020 0011 "IS" Series Number
 623          sprintf(charImagePosition,"%04d",serieNumber+1);
 624          file->InsertEntryString(charImagePosition,0x0020,0x0011, "IS");
 625
 626 //0020 0011 "IS" Instance Number
 627          sprintf(charImagePosition,"%04d",imageNumber);
 628          file->InsertEntryString(charImagePosition,0x0020,0x0013, "IS");
 629
 630 //0018 1060 "DS" Time Trigger
 631          sprintf(charImagePosition,"%f",currentTime);
 632          file->InsertEntryString(charImagePosition,0x0018,0x1060, "DS");
 633    // file->Print();
 634
 635     GDCM_NAME_SPACE::FileHelper *fh;
 636     fh = GDCM_NAME_SPACE::FileHelper::New(file);
 637     // cast is just to avoid warnings (*no* conversion)
 638     fh->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 639     fh->SetWriteModeToRaw();
 640     fh->SetWriteTypeToDcmExplVR();
 641
 642     fh->SetWriteTypeToDcmExplVR();
 643
 644     char numero[10];
 645     sprintf(numero, "%02d", l);
 646     std::string fileName = imageName + ".Anatomical." + numero + ".dcm";
 647     std::cout << "fileName " << fileName << std::endl;
 648
 649     if( !fh->Write(fileName))
 650        std::cout << "Failed for [" << fileName << "]\n"
 651                  << "           File is unwrittable" << std::endl;
 652
 653     delete img;
 654     currentTime += temporalResolution;
 655     imageNumber ++;
 656
 657   } // end loop on frames
 658
 659    from.close();
 660 }