Example/DenseToDicom.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: DenseToDicom.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2007/06/21 15:06:13 $
   7   Version:   $Revision: 1.3 $
   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
  28 #include "gdcmArgMgr.h"
  29
  30 /**
  31   * \brief
  32   *          - explores recursively the given directory
  33   *          - examines the ".txt" files
  34   *          - Converts the files into 16 bits Dicom Files,
  35   */
  36
  37
  38 void Load(std::ifstream &from, std::string imageName);
  39 //std::ifstream& eatwhite(std::ifstream& is);
  40
  41 /*
  42 window center (level) and width, are defined in the DICOM
  43 standard very precisely as follows (see PS 3.3 C.11.2.1.2):
  44 >
  45 >"These Attributes are applied according to the following pseudo-code,
  46 >where :
  47 x is the input value,
  48 y is an output value with a range from ymin to ymax,
  49 c is Window Center (0028,1050)
  50 w is Window Width (0028,1051):
  51 >
  52 >           if      (x  <= c - 0.5 - (w-1)/2), then y = ymin
  53 >           else if (x > c - 0.5 + (w-1)/2), then y = ymax,
  54 >           else    y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin
  55
  56 */
  57 /*
  58 From:   David Clunie - view profile
  59 Date:   Thurs, Jun 1 2006 3:03 pm
  60 Email:  David Clunie <dclu...@dclunie.com>
  61 Groups: comp.protocols.dicom
  62
  63 The value of x is the output of the preceding step, the so-called
  64 "modality LUT", which may either be:
  65
  66 - identity (no rescale values or Modality LUT, or the SOP Class is
  67   PET and rescale values are ignored), in which case x is the stored
  68   pixel value in (7FE0,0010)
  69
  70 - Rescale Slope and Intercept (as typically used in CT), in which
  71   case x is the value obtained from applying the rescale values to
  72   the stored pixel value
  73
  74 - an actual LUT, in which case x is the value stored in the LUT
  75   corresponding to the LUT index value that is the stored pixel
  76   value
  77
  78 The ymin and ymax are intended to represent the output range; for
  79 example, if the hypothetical Presentation LUT step that follows
  80 the VOI LUT (window) stage is an identity operation, then the
  81 ymin and ymax represent P-Values, which might be the range of
  82 digital driving levels for your display (calibrated to the GSDF),
  83 in the 8-bit wide output case ranging from 0 to 255, for example.
  84
  85 The terms brightness and contrast are not used in radiology imaging
  86 - the window center and width are used instead.
  87 */
  88
  89 int main(int argc, char *argv[])
  90 {
  91    START_USAGE(usage)
  92    " \n DenseToDicom :\n                                                      ",
  93    " - explores recursively the given directory,                              ",
  94    "         - examines the '.txt' files                                      ",
  95    "         - Converts the files into 16 bits Dicom Files,                   ",
  96    " usage:                                                                   ",
  97    " DenseToDicom dirin=rootDirectoryName                                     ",
  98    "                [listonly] [verbose] [debug]                              ",
  99    "                                                                          ",
 100    " verbose  : user wants to run the program in 'verbose mode'               ",
 101    " debug    : *developer*  wants to run the program in 'debug mode'         ",
 102    FINISH_USAGE
 103
 104    // ----- Initialize Arguments Manager ------
 105
 106    GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
 107
 108    if (argc == 1 || am->ArgMgrDefined("usage"))
 109    {
 110       am->ArgMgrUsage(usage); // Display 'usage'
 111       delete am;
 112       return 0;
 113    }
 114
 115    const char *dirNamein;
 116    dirNamein  = am->ArgMgrGetString("dirin",".");
 117
 118    if (am->ArgMgrDefined("debug"))
 119       GDCM_NAME_SPACE::Debug::DebugOn();
 120
 121    int verbose  = am->ArgMgrDefined("verbose");
 122    int listonly = am->ArgMgrDefined("listonly");
 123
 124    // if unused Param we give up
 125    if ( am->ArgMgrPrintUnusedLabels() )
 126    {
 127       am->ArgMgrUsage(usage);
 128       delete am;
 129       return 0;
 130    }
 131    delete am;  // we don't need Argument Manager any longer
 132
 133    // ----- Begin Processing -----
 134
 135    if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirNamein) )
 136    {
 137       std::cout << "KO : [" << dirNamein << "] is not a Directory." << std::endl;
 138       return 0;
 139    }
 140    else
 141    {
 142       std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl;
 143    }
 144
 145    std::string strDirNamein(dirNamein);
 146    GDCM_NAME_SPACE::DirList dirList(strDirNamein, true); // get recursively the list of files
 147
 148    if (listonly)
 149    {
 150       std::cout << "------------List of found files ------------" << std::endl;
 151       dirList.Print();
 152       std::cout << std::endl;
 153       return 0;
 154     }
 155
 156    std::string filenameout;
 157
 158    GDCM_NAME_SPACE::DirListType fileNames;
 159    fileNames = dirList.GetFilenames();
 160    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 161                                     it != fileNames.end();
 162                                   ++it)
 163    {
 164       std::ifstream from( (*it).c_str() );
 165       if ( !from )
 166       {
 167          std::cout << "Can't open file" << *it << std::endl;
 168          //return 0;
 169       }
 170       else
 171       {
 172          filenameout = *it + ".dcm";
 173          std::cout << "Success in open file" << *it << std::endl;
 174          Load(from, filenameout);
 175          //return 0;
 176       }
 177    }
 178 }
 179
 180
 181 void Load(std::ifstream &from, std::string imageName)
 182 {
 183    if (!from)
 184       return;
 185  /*
 186    uint16_t group;
 187    uint16_t elem;
 188    VRKey vr;
 189    TagName vm;
 190    TagName name;
 191 */
 192    std::string str1;
 193    int nx, ny;
 194    // Get nx, ny
 195    from >> str1;
 196    from >> str1;
 197    from >> str1;
 198    from >> nx;
 199    from >> str1;
 200    from >> ny;
 201
 202    std::cout << "nx " << nx << " ny " << ny << std::endl;
 203
 204    // Skip 2 lines.
 205    std::getline(from, str1);
 206    std::cout << "[" << str1 << "]" << std::endl;
 207    std::getline(from, str1);
 208    std::cout << "[" << str1 << "]" << std::endl;
 209    std::getline(from, str1);
 210    std::cout << "[" << str1 << "]" << std::endl;
 211
 212    //float *f = new float(nx*ny);
 213    float *f = (float *) malloc(nx*ny*sizeof(float));
 214   // float mini = FLT_MAX, maxi = FLT_MIN;
 215    float val;
 216    for( int j=0; j<ny; j++)
 217    {
 218       int l =0;
 219       for (int i=0; i<nx; i++)
 220       {
 221          //eatwhite(from);
 222
 223     char c;
 224     for (;;) {
 225       if (!from.get(c))
 226         break;
 227       if (!isspace(c)) {
 228         from.putback(c);
 229         break;
 230       }
 231     }
 232          from >> str1;
 233          val = (float)atof(str1.c_str());
 234          *(f+j*nx+i) = val;
 235
 236          // In our concern, just *100 all the values is enough!
 237          /*
 238         if (val < mini)
 239            mini = val;
 240         if (val > maxi)
 241             maxi = val;
 242          //std::cout << val << " ";
 243         */
 244         if(from.eof())
 245         {
 246             std::cout << "Missing values at [" << j <<"," << i << "]"
 247                       << std::endl;
 248            break;
 249          }
 250          l++;
 251       }
 252       //std::cout << std::endl << " line nb : " << j
 253       //          << " line length : " << l << std::endl;
 254
 255     }
 256
 257    // std::cout << "mini : "<< mini  << " maxi : " << maxi << std::endl;
 258 /*
 259 // values are expressed as %.
 260 // It's useless to rescale them as uint16_t : just *100
 261     uint16_t *img = new uint16_t[ny*nx];
 262     uint16_t *ptr = img;
 263     float *tmp = f;
 264     float div = maxi-mini;
 265     std::cout << "div : " << div << " mini : " << mini << std::endl;
 266     for( int k=0; k<ny*nx; k++)
 267     {
 268        *ptr = (uint16_t)((*tmp * 65535.0) / div);
 269        tmp ++;
 270        ptr++;
 271     }
 272 */
 273
 274     uint16_t *img = new uint16_t[ny*nx];
 275     uint16_t *ptr = img;
 276     float *tmp = f;
 277     for( int k=0; k<ny*nx; k++)
 278     {
 279        *ptr = (uint16_t)(*tmp *100);
 280        tmp ++;
 281        ptr++;
 282     }
 283
 284  /*
 285  Vtk/Python Theralys code for Rescale Slope, Rescale Intercept
 286
 287 mini,maxi=data.GetScalarRange()
 288                      minPrec=2
 289            maxPow=16-minPrec
 290            if mini<=-pow(2,maxPow):mini=-pow(2,maxPow)+1
 291            if maxi>pow(2,maxPow):maxi=pow(2,maxPow)
 292                      maxi=math.ceil(maxi)+1
 293            mini=math.ceil(mini-.5)
 294                      slope,intercept=(pow(2,16)-1)/(maxi-mini),-mini
 295            print mini,maxi,"->",slope,intercept
 296                      sl=math.ceil(1e6/slope)
 297            ok=1
 298            while(ok):
 299                val=1/(sl/1e6)
 300                print sl
 301                if val==int(val):ok=0
 302                else:sl+=1
 303                      print sl/1e6, 1/(sl/1e6)
 304            slope=sl/1e6
 305                      invSlope=1/slope
 306            interc=-math.ceil(intercept*1e6)/1e6
 307                      print invSlope,slope,interc
 308                                imageSC=vtkImageShiftScale()
 309            imageSC.SetOutputScalarTypeToUnsignedShort()
 310            imageSC.SetShift(-interc)
 311            imageSC.SetScale(invSlope)
 312            imageSC.SetInput(data)
 313            imageSC.Update()
 314  */
 315  // GDCM_NAME_SPACE::Debug::DebugOn();
 316
 317        std::ostringstream str;
 318        GDCM_NAME_SPACE::File *file;
 319        file = GDCM_NAME_SPACE::File::New();
 320
 321   // Set the image size
 322         str.str("");
 323         str << nx;
 324         file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
 325         str.str("");
 326         str << ny;
 327         file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
 328   // Set the pixel type
 329   //      16; //8, 16, 32
 330         file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
 331         str.str("");
 332         str << 16; // may be 12 or 16 if componentSize =16
 333         file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
 334
 335         file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit
 336   // Set the pixel representation // 0/1
 337
 338         file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation
 339   // Set the samples per pixel // 1:Grey level, 3:RGB
 340         file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel
 341
 342 /*
 343   // Set Rescale Intercept
 344         str.str("");
 345         str << div;
 346         file->InsertEntryString(str.str(),0x0028,0x1052,"DS");
 347
 348   // Set Rescale Slope
 349         str.str("");
 350         str << mini;
 351         file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
 352 */
 353
 354     file->Print();
 355
 356     GDCM_NAME_SPACE::FileHelper *fileH;
 357     fileH = GDCM_NAME_SPACE::FileHelper::New(file);
 358     // cast is just to avoid warnings (*no* conversion)
 359     fileH->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t));
 360     fileH->SetWriteModeToRaw();
 361     fileH->SetWriteTypeToDcmExplVR();
 362
 363     fileH->Write(imageName);
 364
 365       if( !fileH->Write(imageName))
 366              std::cout << "Failed for [" << imageName << "]\n"
 367                 << "           File is unwrittable" << std::endl;
 368
 369    delete img;
 370    from.close();
 371 }