Example/Dense2007ToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: Dense2007ToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2008/03/17 13:16:10 $
   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 #if defined(__BORLANDC__)
  24 #include <ctype.h>
  25 #endif
  26
  27 #include "gdcmFile.h"
  28 #include "gdcmFileHelper.h"
  29 #include "gdcmDebug.h"
  30 #include "gdcmDirList.h"
  31
  32 #include "gdcmArgMgr.h"
  33
  34 /**
  35   * \brief
  36   *  Converts the "Dense" ".txt" (2007 version)  files into 16 bits Dicom Files,
  37   * Hope they don't change soon!
  38   */
  39
  40
  41 void LoadPeakStrain(std::ifstream &from, std::string imageName);
  42 void LoadStrain(std::ifstream &from, std::string imageName);
  43 void MakeDicomImage(float *tabVal, float *X, float *Y, float *Z, int NP, std::string dcmImageName);
  44
  45 int verbose;
  46
  47 int main(int argc, char *argv[])
  48 {
  49    START_USAGE(usage)
  50    " \n Dense2007ToDicom :\n                                                  ",
  51    "        Converts the '.txt' files into 16 bits Dicom Files,               ",
  52    " usage:                                                                   ",
  53    " Dense2007ToDicom strain=...strain.txt  peak_strain=...peak_strain.txt    ",
  54    "                 [verbose] [debug]                                        ",
  55    "                                                                          ",
  56    " verbose  : user wants to run the program in 'verbose mode'               ",
  57    " debug    : *developer*  wants to run the program in 'debug mode'         ",
  58    FINISH_USAGE
  59
  60    // ----- Initialize Arguments Manager ------
  61
  62    GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
  63
  64    if (argc == 1 || am->ArgMgrDefined("usage"))
  65    {
  66       am->ArgMgrUsage(usage); // Display 'usage'
  67       delete am;
  68       return 0;
  69    }
  70    // Seems that ArgMgrWantString doesn't work on MacOS
  71    if(!am->ArgMgrDefined("strain"))
  72    {
  73       std::cout << "strain is mandatory" << std::endl;
  74       exit(0);
  75    }
  76    if(!am->ArgMgrDefined("peak_strain"))
  77    {
  78       std::cout << "peak_strain is mandatory" << std::endl;
  79       exit(0);
  80    }
  81
  82    const char *strain      = am->ArgMgrWantString("strain",usage);
  83    const char *peak_strain = am->ArgMgrWantString("peak_strain",usage);
  84
  85    if (am->ArgMgrDefined("debug"))
  86       GDCM_NAME_SPACE::Debug::DebugOn();
  87
  88    verbose  = am->ArgMgrDefined("verbose");
  89
  90    // if unused Param we give up
  91    if ( am->ArgMgrPrintUnusedLabels() )
  92    {
  93       am->ArgMgrUsage(usage);
  94       delete am;
  95       return 0;
  96    }
  97    delete am;  // we don't need Argument Manager any longer
  98
  99    // ----- Begin Processing -----
 100
 101
 102    std::ifstream fromPeakStrain( peak_strain );
 103    if ( !fromPeakStrain )
 104    {
 105       std::cout << "Can't open file" << peak_strain << std::endl;
 106       exit(0);
 107    }
 108
 109    std::ifstream fromStrain( strain );
 110    if ( !fromStrain )
 111    {
 112       std::cout << "Can't open file" << strain << std::endl;
 113       exit(0);
 114    }
 115
 116    std::cout << "Success in open file" << peak_strain << std::endl;
 117    LoadPeakStrain(fromPeakStrain, peak_strain);
 118    fromPeakStrain.close();
 119
 120    std::cout << "Success in open file" << strain << std::endl;
 121    LoadStrain(fromStrain, strain);
 122    fromStrain.close();
 123    return 1;
 124 }
 125
 126 // =====================================================================================================================
 127
 128 void LoadPeakStrain(std::ifstream &from, std::string textFileName)
 129 {
 130 // in sax_base_slice0_peak_strain.txt :
 131
 132 /*
 133 Number of material points (NP) = 181
 134 Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.3243 3.19392
 135 88.2381
 136 Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
 137 Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
 138 Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
 139 The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
 140 followed by their peak Ecc strain, an array of NP elements,
 141 followed by their peak Err strain, an array of NP elements,
 142 followed by their peak E11 strain, an array of NP elements,
 143 followed by their Peak E22 strain, an array of NP elements,
 144       42.0000      10.0000     0.000000
 145       ...
 146     -0.154905   -0.0840482    -0.157350    -0.221403    -0.168118    -0.131331
 147     -0.153781    -0.148481    -0.166602    -0.232858    -0.222650    -0.213712
 148     ...
 149 */
 150
 151    if (!from)
 152       return;
 153
 154    std::string str1;
 155    int NP;
 156
 157    //Number of material points (NP) = 181
 158     from >> str1;
 159     from >> str1;
 160     from >> str1;
 161     from >> str1;
 162     from >> str1;
 163     from >> str1;
 164     from >> NP;
 165
 166     std::cout << "NP : " << NP << std::endl;
 167
 168    //Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.3243 3.19392 88.2381
 169      from >> str1;
 170      from >> str1;
 171      from >> str1;
 172      from >> str1;
 173      from >> str1;
 174      from >> str1;
 175      from >> str1;
 176      from >> str1;
 177      from >> str1;
 178      from >> str1;
 179      from >> str1;
 180
 181      float readout,  phase_enc, slice_sel;
 182      from >> readout;
 183      from >> phase_enc;
 184      from >> slice_sel;
 185      std::cout << " readout " << readout << " phase_enc " << phase_enc << " slice_sel " << slice_sel << std::endl;
 186
 187     // Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
 188
 189     from >> str1;
 190     from >> str1;
 191     from >> str1;
 192     from >> str1;
 193     from >> str1;
 194
 195     float readoutX, readoutY, readoutZ;
 196     from >> readoutX;
 197     from >> readoutY;
 198     from >> readoutZ;
 199     std::cout << " readoutX " << readoutX <<  " readoutY " << readoutY <<  " readoutZ " << readoutZ << std::endl;
 200
 201 // Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
 202
 203      from >> str1;
 204      from >> str1;
 205      from >> str1;
 206      from >> str1;
 207      from >> str1;
 208      from >> str1;
 209
 210     float phase_encX, phase_encY, phase_encZ;
 211     from >> phase_encX;
 212     from >> phase_encY;
 213     from >> phase_encZ;
 214     std::cout << " phase_encX " << phase_encX <<  " phase_encY " << phase_encY <<  " phase_encZ " << phase_encZ << std::endl;
 215
 216 // Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
 217      from >> str1;
 218      from >> str1;
 219      from >> str1;
 220      from >> str1;
 221      from >> str1;
 222      from >> str1;
 223
 224     float slice_selX, slice_selY, slice_selZ;
 225     from >> slice_selX;
 226     from >> slice_selY;
 227     from >> slice_selZ;
 228     std::cout << " slice_selX " << slice_selX <<  " slice_selY " << slice_selY <<  " slice_selZ " << slice_selZ << std::endl;
 229
 230
 231
 232 // Skip 5 lines :
 233 /*
 234 The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
 235 followed by their peak Ecc strain, an array of NP elements,
 236 followed by their peak Err strain, an array of NP elements,
 237 followed by their peak E11 strain, an array of NP elements,
 238 followed by their Peak E22 strain, an array of NP elements,
 239 */
 240
 241 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 242    std::getline(from, str1);
 243    std::cout << "[" << str1 << "]" << std::endl;
 244 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 245    std::getline(from, str1);
 246    std::cout << "[" << str1 << "]" << std::endl;
 247 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 248    std::getline(from, str1);
 249    std::cout << "[" << str1 << "]" << std::endl;
 250 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 251    std::getline(from, str1);
 252    std::cout << "[" << str1 << "]" << std::endl;
 253 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 254    std::getline(from, str1);
 255    std::cout << "[" << str1 << "]" << std::endl;
 256 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 257    std::getline(from, str1);
 258    std::cout << "[" << str1 << "]" << std::endl;
 259 std::cout << "------------stop skipping ---------------- " << std::endl;
 260
 261    float *X = new float[NP];
 262    float *Y = new float[NP];
 263    float *Z = new float[NP];
 264
 265    char c;
 266    int i;
 267    for (i=0; i<NP; i++) {
 268
 269       from >> X[i];
 270       for (;;) {
 271         if (!from.get(c))
 272           break;
 273         if (!isspace(c)) {
 274           from.putback(c);
 275           break;
 276         }
 277      }
 278
 279       from >> Y[i];
 280       for (;;) {
 281         if (!from.get(c))
 282           break;
 283         if (!isspace(c)) {
 284           from.putback(c);
 285           break;
 286         }
 287      }
 288       from >> Z[i];
 289
 290    } // end for i<NP
 291
 292    std::cout << "--------------- Ecc_strain ------------------" << std::endl;
 293    float *ecc_strain = new float[NP];
 294    for (i=0; i<NP; i++) {
 295        from >> ecc_strain[i];
 296        if (verbose)
 297        std::cout <<  ecc_strain[i] <<  std::endl;
 298    }
 299
 300    std::cout << "--------------- Err_strain ------------------" << std::endl;
 301    float *err_strain = new float[NP];
 302    for (i=0; i<NP; i++) {
 303        from >> err_strain[i];
 304        if (verbose)
 305        std::cout <<  err_strain[i] <<  std::endl;
 306    }
 307
 308    std::cout << "--------------- E11_strain ------------------" << std::endl;
 309    float *e11_strain = new float[NP];
 310    for (i=0; i<NP; i++) {
 311        from >> e11_strain[i];
 312        if (verbose)
 313        std::cout <<  e11_strain[i] <<  std::endl;
 314    }
 315
 316    std::cout << "--------------- E22_strain ------------------" << std::endl;
 317    float *e22_strain = new float[NP];
 318    for (i=0; i<NP; i++) {
 319        from >> e22_strain[i];
 320        if (verbose)
 321        std::cout <<  e22_strain[i] <<  std::endl;
 322    }
 323
 324    std::string dcmImageName;
 325
 326 //followed by their peak Ecc strain, an array of NP elements,
 327    dcmImageName = textFileName + "_peak_Ecc_strain.dcm";
 328    MakeDicomImage(ecc_strain, X, Y, Z, NP, dcmImageName);
 329
 330 //followed by their peak Err strain, an array of NP elements,
 331    dcmImageName = textFileName + "_peak_Err_strain.dcm";
 332    MakeDicomImage(err_strain, X, Y, Z, NP, dcmImageName);
 333
 334 //followed by their peak E11 strain, an array of NP elements,
 335    dcmImageName = textFileName + "_peak_E11_strain.dcm";
 336    MakeDicomImage(e11_strain, X, Y, Z, NP, dcmImageName);
 337
 338 //followed by their Peak E22 strain, an array of NP elements,
 339    dcmImageName = textFileName + "_peak_E22_strain.dcm";
 340    MakeDicomImage(e22_strain, X, Y, Z, NP, dcmImageName);
 341 }
 342
 343 // =====================================================================================================================
 344
 345 void LoadStrain(std::ifstream &from, std::string textFileName)
 346 {
 347
 348 // in sax_base_slice0_strain.txt :
 349 /*
 350 Number of cine frames = 18
 351 Temporal resolution = 32.0000 ms
 352 First frame starts at 48.0000 ms
 353 Number of material points (NP) = 181
 354 Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.324341 3.193918 88.238113
 355 Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
 356 Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
 357 Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
 358 The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
 359 followed by their Ecc strain, an array of dimensions(NP, number of cine frames),
 360 followed by their Err strain, an array of dimensions(NP, number of cine frames),
 361 followed by their E11 strain, an array of dimensions(NP, number of cine frames),
 362 followed by their E22 strain, an array of dimensions(NP, number of cine frames),
 363 Note that RV Err, E11 and E22 strains are not calculated due to the small thickness.
 364       42.0000      10.0000     0.000000
 365       44.0000      10.0000     0.000000
 366       ...
 367    -0.0622793   -0.0840482    -0.157350    -0.196722    -0.105844    -0.131331
 368     -0.153781  -0.00940573   -0.0542236    -0.100403   -0.0369671   -0.0696840
 369 */
 370
 371    if (!from)
 372       return;
 373
 374    std::string str1;
 375    int NP;    // Number of Pints
 376    int NCF;   // Number of cine frames
 377    float TR;  // Temporal resolution
 378    float FFS; // First frame starts
 379
 380    // Number of cine frames = 18
 381     from >> str1;
 382     from >> str1;
 383     from >> str1;
 384     from >> str1;
 385     from >> str1;
 386     from >> NCF;
 387
 388    // Temporal resolution = 32.0000 ms
 389     from >> str1;
 390     from >> str1;
 391     from >> str1;
 392     from >> TR;
 393     from >> str1;
 394
 395    // First frame starts at 48.0000 ms
 396     from >> str1;
 397     from >> str1;
 398     from >> str1;
 399     from >> str1;
 400     from >> FFS;
 401     from >> str1;
 402
 403    //Number of material points (NP) = 181
 404     from >> str1;
 405     from >> str1;
 406     from >> str1;
 407     from >> str1;
 408     from >> str1;
 409     from >> str1;
 410     from >> NP;
 411
 412     std::cout << "NP : " << NP << std::endl;
 413
 414    //Origin of (readout, phase enc, slice sel) coordinates in 3D = 87.324341 3.193918 88.238113
 415      from >> str1;
 416      from >> str1;
 417      from >> str1;
 418      from >> str1;
 419      from >> str1;
 420      from >> str1;
 421      from >> str1;
 422      from >> str1;
 423      from >> str1;
 424      from >> str1;
 425      from >> str1;
 426
 427      float readout,  phase_enc, slice_sel;
 428      from >> readout;
 429      from >> phase_enc;
 430      from >> slice_sel;
 431      std::cout << " readout " << readout << " phase_enc " << phase_enc << " slice_sel " << slice_sel << std::endl;
 432
 433     // Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
 434
 435     from >> str1;
 436     from >> str1;
 437     from >> str1;
 438     from >> str1;
 439     from >> str1;
 440
 441     float readoutX, readoutY, readoutZ;
 442     from >> readoutX;
 443     from >> readoutY;
 444     from >> readoutZ;
 445     std::cout << " readoutX " << readoutX <<  " readoutY " << readoutY <<  " readoutZ " << readoutZ << std::endl;
 446
 447 // Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
 448
 449      from >> str1;
 450      from >> str1;
 451      from >> str1;
 452      from >> str1;
 453      from >> str1;
 454      from >> str1;
 455
 456     float phase_encX, phase_encY, phase_encZ;
 457     from >> phase_encX;
 458     from >> phase_encY;
 459     from >> phase_encZ;
 460     std::cout << " phase_encX " << phase_encX <<  " phase_encY " << phase_encY <<  " phase_encZ " << phase_encZ << std::endl;
 461
 462 // Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
 463      from >> str1;
 464      from >> str1;
 465      from >> str1;
 466      from >> str1;
 467      from >> str1;
 468      from >> str1;
 469
 470     float slice_selX, slice_selY, slice_selZ;
 471     from >> slice_selX;
 472     from >> slice_selY;
 473     from >> slice_selZ;
 474     std::cout << " slice_selX " << slice_selX <<  " slice_selY " << slice_selY <<  " slice_selZ " << slice_selZ << std::endl;
 475
 476
 477
 478 // Skip 6 lines :
 479 /*
 480 The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
 481 followed by their Ecc strain, an array of dimensions(NP, number of cine frames),
 482 followed by their Err strain, an array of dimensions(NP, number of cine frames),
 483 followed by their E11 strain, an array of dimensions(NP, number of cine frames),
 484 followed by their E22 strain, an array of dimensions(NP, number of cine frames),
 485 Note that RV Err, E11 and E22 strains are not calculated due to the small thickness.
 486 */
 487 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 488    std::getline(from, str1);
 489    std::cout << "[" << str1 << "]" << std::endl;
 490 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 491    std::getline(from, str1);
 492    std::cout << "[" << str1 << "]" << std::endl;
 493 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 494    std::getline(from, str1);
 495    std::cout << "[" << str1 << "]" << std::endl;
 496 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 497    std::getline(from, str1);
 498    std::cout << "[" << str1 << "]" << std::endl;
 499 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 500    std::getline(from, str1);
 501    std::cout << "[" << str1 << "]" << std::endl;
 502 std::cout << "------------start skipping 1 line---------------- " << std::endl;
 503    std::getline(from, str1);
 504    std::cout << "[" << str1 << "]" << std::endl;
 505 std::cout << "------------stop skipping ---------------- " << std::endl;
 506    std::getline(from, str1);
 507    std::cout << "[" << str1 << "]" << std::endl;
 508 std::cout << "------------stop skipping ---------------- " << std::endl;
 509
 510    float *X = new float[NP];
 511    float *Y = new float[NP];
 512    float *Z = new float[NP];
 513
 514    char c;
 515    int i;
 516    for (i=0; i<NP; i++) {
 517
 518       from >> X[i];
 519       for (;;) {
 520         if (!from.get(c))
 521           break;
 522         if (!isspace(c)) {
 523           from.putback(c);
 524           break;
 525         }
 526      }
 527       from >> Y[i];
 528       for (;;) {
 529         if (!from.get(c))
 530           break;
 531         if (!isspace(c)) {
 532           from.putback(c);
 533           break;
 534         }
 535      }
 536       from >> Z[i];
 537    }
 538
 539    std::cout << "--------------- Ecc_strain ------------------" << std::endl;
 540    float *ecc_strain = new float[NP];
 541    for (i=0; i<NP; i++) {
 542        from >> ecc_strain[i];
 543        if (verbose)
 544        std::cout <<  ecc_strain[i] <<  std::endl;
 545    }
 546
 547    std::cout << "--------------- Err_strain ------------------" << std::endl;
 548    float *err_strain = new float[NP];
 549    for (i=0; i<NP; i++) {
 550        from >> err_strain[i];
 551        if (verbose)
 552        std::cout <<  err_strain[i] <<  std::endl;
 553    }
 554
 555    std::cout << "--------------- E11_strain ------------------" << std::endl;
 556    float *e11_strain = new float[NP];
 557    for (i=0; i<NP; i++) {
 558        from >> e11_strain[i];
 559        if (verbose)
 560        std::cout <<  e11_strain[i] <<  std::endl;
 561    }
 562
 563    std::cout << "--------------- E22_strain ------------------" << std::endl;
 564    float *e22_strain = new float[NP];
 565    for (i=0; i<NP; i++) {
 566        from >> e22_strain[i];
 567        if (verbose)
 568        std::cout <<  e22_strain[i] <<  std::endl;
 569    }
 570
 571    std::string dcmImageName;
 572
 573 //followed by their Ecc strain, an array of NP elements,
 574    dcmImageName = textFileName + "_Ecc_strain.dcm";
 575    MakeDicomImage(ecc_strain, X, Y, Z, NP, dcmImageName);
 576
 577 //followed by their  Err strain, an array of NP elements,
 578    dcmImageName = textFileName + "_Err_strain.dcm";
 579    MakeDicomImage(err_strain, X, Y, Z, NP, dcmImageName);
 580
 581 //followed by their E11 strain, an array of NP elements,
 582    dcmImageName = textFileName + "_E11_strain.dcm";
 583    MakeDicomImage(e11_strain, X, Y, Z, NP, dcmImageName);
 584
 585 //followed by their E22 strain, an array of NP elements,
 586    dcmImageName = textFileName + "_E22_strain.dcm";
 587    MakeDicomImage(e22_strain, X, Y, Z, NP, dcmImageName);
 588
 589 }
 590
 591
 592 // =====================================================================================================================
 593
 594
 595 void MakeDicomImage(float *tabVal, float *X, float *Y, float *Z, int NP, std::string dcmImageName)
 596 {
 597    float minX = 99999., minY = 99999., minZ = 99999.;
 598    float maxX = 0., maxY = 0., maxZ = 0.;
 599    int i;
 600
 601    for (i=0; i<NP; i++) {
 602       // std::cout << X[i] << " " << Y[i] << " " << Z[i] <<  std::endl;
 603       if(maxX < X[i])
 604          maxX = X[i];
 605       if(maxY < Y[i])
 606          maxY = Y[i];
 607       if(maxZ < Z[i])
 608          maxZ = Z[i];
 609
 610       if(minX > X[i])
 611          minX = X[i];
 612       if(minY > Y[i])
 613          minY = Y[i];
 614       if(minZ > Z[i])
 615          minZ = Z[i];
 616    }
 617    std::cout << "Min X,Y,Z " << minX << " " << minY << " " << minZ <<  std::endl;
 618    std::cout << "Max X,Y,Z " << maxX << " " << maxY << " " << maxZ <<  std::endl;
 619    std::cout << "Size X,Y,Z " << maxX-minX << " " << maxY-minY << " " << maxZ-minZ <<  std::endl;
 620
 621 //   uint16_t *img = new uint16_t[int(maxX+0.5)*int(maxY+0.5)];
 622    uint16_t *img = new uint16_t[int(maxX*4.)*int(maxY*4.)];
 623
 624    // Set whole image to 0
 625    for(int i3=0;i3<int(maxX*4.)*int(maxY*4.);i3++)
 626 //   for(int i3=0;i3<int(maxX+0.5)*int(maxY+0.5);i3++)
 627       img[i3] = 0;
 628
 629    for(int i2=0; i2<NP; i2++) {
 630
 631       int ordX = int(X[i2]*4.-30);
 632       int ordY = int(maxY*4.) - int(Y[i2]*4.)+30;
 633 //   img[ /*int(maxX) -*/ int(X[i2]+0.5-1)   +  (int(maxY+0.5) - int(Y[i2]+0.5))   * int(maxX+0.5) ] = int(tabVal[i2]*100);
 634
 635 //   img[ /*int(maxX) -*/ int(X[i2]*4.-30)   +  (int(maxY*4.) - int(Y[i2]*4.)+30)   * int(maxX*4.) ] = int(tabVal[i2]*100);
 636    img[ /*int(maxX) -*/ ordX   +  ordY   * int(maxX*4.) ] = int(tabVal[i2]*100);
 637
 638       // Try to round up, just to see.
 639        for(int iii=ordY-3; iii<ordY+4; iii++) for(int jjj=ordX-3; jjj<ordX+4; jjj++)
 640           img[  jjj  +  iii   * int(maxX*4.) ] = int(tabVal[i2]*100);
 641       std::cout << int(X[i2]*4.) << " " << int(Y[i2]*4.) << " = " << int(tabVal[i2]*100) << std::endl;
 642    }
 643
 644  // GDCM_NAME_SPACE::Debug::DebugOn();
 645
 646    std::ostringstream str;
 647
 648    GDCM_NAME_SPACE::File *file;
 649    file = GDCM_NAME_SPACE::File::New();
 650
 651   // Set the image size
 652    str.str("");
 653    str << (int)(maxX*4.);
 654    file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 655    str.str("");
 656    str << (int)(maxY*4.);
 657    file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 658
 659   // Set the pixel type
 660   //      16; //8, 16, 32
 661    file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 662    str.str("");
 663    str << 16; // may be 12 or 16 if componentSize =16
 664    file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 665    file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 666
 667   // Set the pixel representation // 0/1 , 0=unsigned
 668    file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
 669   // Set the samples per pixel // 1:Grey level, 3:RGB
 670    file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 671
 672 /*
 673   // Set Rescale Intercept
 674         str.str("");
 675         str << div;
 676         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 677
 678   // Set Rescale Slope
 679         str.str("");
 680         str << mini;
 681         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 682 */
 683
 684    GDCM_NAME_SPACE::FileHelper *fileH;
 685    fileH = GDCM_NAME_SPACE::FileHelper::New(file);
 686    // cast is just to avoid warnings (*no* conversion)
 687    //fileH->SetImageData((uint8_t *)img,int(maxX*maxY)*sizeof(uint16_t)); // troubles when maxX, mayY are *actually* float!
 688    fileH->SetImageData((uint8_t *)img,int(maxX*4.)*int(maxY*4.)*sizeof(uint16_t));
 689    fileH->SetWriteModeToRaw();
 690    fileH->SetWriteTypeToDcmExplVR();
 691
 692    if( !fileH->Write(dcmImageName))
 693       std::cout << "Failed for [" << dcmImageName << "]\n"
 694                 << "           File is unwrittable" << std::endl;
 695
 696    //file->Print();
 697
 698    delete img;
 699    file->Delete();
 700    fileH->Delete();
 701 }