vtk/vtkGdcmReader.cxx

   1 // $Header: /cvs/public/gdcm/vtk/vtkGdcmReader.cxx,v 1.13 2003/07/01 10:04:37 frog Exp $
   2 #include <stdio.h>
   3 #include <vtkObjectFactory.h>
   4 #include <vtkImageData.h>
   5 #include <vtkPointData.h>
   6 #include "vtkGdcmReader.h"
   7 #include "gdcm.h"
   8
   9 vtkGdcmReader::vtkGdcmReader()
  10 {
  11   // Constructor
  12 }
  13
  14 //----------------------------------------------------------------------------
  15 vtkGdcmReader::~vtkGdcmReader()
  16 {
  17   // FIXME free memory
  18   this->FileNameList.clear();
  19 }
  20
  21 //----------------------------------------------------------------------------
  22 // Adds a file name to the list of images to read.
  23 void vtkGdcmReader::AddFileName(const char* name)
  24 {
  25   // We need to bypass the const pointer [since list<>.push_bash() only
  26   // takes a char* (but not a const char*)] by making a local copy:
  27   char * LocalName = new char[strlen(name) + 1];
  28   strcpy(LocalName, name);
  29   this->FileNameList.push_back(LocalName);
  30   this->Modified();
  31   delete[] LocalName;
  32 }
  33
  34 //----------------------------------------------------------------------------
  35 // Sets up a filename to be read.
  36 void vtkGdcmReader::SetFileName(const char *name) {
  37   vtkImageReader2::SetFileName(name);
  38   // Since we maintain a list of filenames, when building a volume,
  39   // (see vtkGdcmReader::AddFileName), we additionaly need to purge
  40   // this list when we manually positionate the filename:
  41   this->FileNameList.clear();
  42   this->Modified();
  43 }
  44
  45 //----------------------------------------------------------------------------
  46 // vtkGdcmReader can have the file names specified through two ways:
  47 // (1) by calling the vtkImageReader2::SetFileName(), SetFilePrefix() and
  48 //     SetFilePattern()
  49 // (2) By successive calls to vtkGdcmReader::AddFileName()
  50 // When the first method was used by caller we need to update the local
  51 // filename list
  52 void vtkGdcmReader::BuildFileListFromPattern()
  53 {
  54    if ((! this->FileNameList.empty()) && this->FileName )
  55      {
  56      vtkErrorMacro("Both file patterns and AddFileName schemes were used");
  57      vtkErrorMacro("Only the files specified with AddFileName shall be used");
  58      return;
  59      }
  60
  61    if (! this->FileNameList.empty()  )
  62      {
  63      vtkDebugMacro("Using the AddFileName specified files");
  64      return;
  65      }
  66
  67    if (!this->FileName && !this->FilePattern)
  68      {
  69      vtkErrorMacro("FileNames are not set. Either use AddFileName() or");
  70      vtkErrorMacro("specify a FileName or FilePattern.");
  71      return;
  72      }
  73
  74    for (int idx = this->DataExtent[4]; idx <= this->DataExtent[5]; ++idx)
  75      {
  76      this->ComputeInternalFileName(idx);
  77      vtkDebugMacro("Adding file " << this->InternalFileName);
  78      this->AddFileName(this->InternalFileName);
  79      }
  80 }
  81
  82 //----------------------------------------------------------------------------
  83 // When more than one filename is specified (i.e. we expect loading
  84 // a stack or volume) we need to check that the corresponding images/volumes
  85 // to be loaded are coherent i.e. to make sure:
  86 //     - they all share the same X dimensions
  87 //     - they all share the same Y dimensions
  88 //     - they all share the same ImageType ( 8 bit signed, or unsigned...)
  89 //
  90 // Eventually, we emit a warning when all the files do NOT share the
  91 // Z dimension, since we can still build a stack but the
  92 // files are not coherent in Z, which is probably a source a trouble...
  93 //   When files are not readable (either the file cannot be opened or
  94 // because gdcm cannot parse it), they are flagged as "GDCM_UNREADABLE".
  95 //   This method returns the total number of planar images to be loaded
  96 // (i.e. an image represents one plane, but a volume represents many planes)
  97 int vtkGdcmReader::CheckFileCoherence()
  98 {
  99    int ReturnedTotalNumberOfPlanes = 0;   // The returned value.
 100
 101    this->BuildFileListFromPattern();
 102    if (this->FileNameList.empty())
 103      {
 104      vtkErrorMacro("FileNames are not set.");
 105      return 0;
 106      }
 107
 108    bool FoundReferenceFile = false;
 109    int  ReferenceNZ = 0;
 110
 111    // Loop on the filenames:
 112    // - check for their existence and gdcm "parasability"
 113    // - get the coherence check done:
 114    for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 115                                         FileName != FileNameList.end();
 116                                       ++FileName)
 117      {
 118      /////// Stage 1: check for file readability:
 119      // Stage 1.1: check for file existence.
 120      FILE *fp;
 121      fp = fopen(FileName->c_str(),"rb");
 122      if (!fp)
 123        {
 124        vtkErrorMacro("Unable to open file " << FileName->c_str());
 125        vtkErrorMacro("Removing this file from readed files "
 126                      << FileName->c_str());
 127        *FileName = "GDCM_UNREADABLE";
 128        continue;
 129        }
 130      fclose(fp);
 131
 132      // Stage 1.2: check for Gdcm parsability
 133      gdcmHeader GdcmHeader(FileName->c_str());
 134      if (!GdcmHeader.IsReadable())
 135        {
 136        vtkErrorMacro("Gdcm cannot parse file " << FileName->c_str());
 137        vtkErrorMacro("Removing this file from readed files "
 138                      << FileName->c_str());
 139        *FileName = "GDCM_UNREADABLE";
 140        continue;
 141        }
 142
 143      // Stage 1.3: further gdcm compatibility on PixelType
 144      std::string type = GdcmHeader.GetPixelType();
 145      if (   (type !=  "8U") && (type !=  "8S")
 146          && (type != "16U") && (type != "16S")
 147          && (type != "32U") && (type != "32S") )
 148        {
 149        vtkErrorMacro("Bad File Type for file" << FileName->c_str());
 150        vtkErrorMacro("Removing this file from readed files "
 151                      << FileName->c_str());
 152        *FileName = "GDCM_UNREADABLE";
 153        continue;
 154        }
 155
 156      /////// Stage 2: check coherence of the set of files
 157      int NX = GdcmHeader.GetXSize();
 158      int NY = GdcmHeader.GetYSize();
 159      int NZ = GdcmHeader.GetZSize();
 160      if (FoundReferenceFile)
 161        {
 162
 163        // Stage 2.1: mandatory coherence stage:
 164        if (   ( NX   != this->NumColumns )
 165            || ( NY   != this->NumLines )
 166            || ( type != this->ImageType ) )
 167          {
 168          vtkErrorMacro("This file is not coherent with previous ones"
 169                        << FileName->c_str());
 170          vtkErrorMacro("Removing this file from readed files "
 171                        << FileName->c_str());
 172          *FileName = "GDCM_UNREADABLE";
 173          continue;
 174          }
 175
 176        // Stage 2.2: optional coherence stage
 177        if ( NZ != ReferenceNZ )
 178          {
 179          vtkErrorMacro("File is not coherent in Z with previous ones"
 180                        << FileName->c_str());
 181          }
 182        else
 183          {
 184          vtkDebugMacro("File is coherent with previous ones"
 185                        << FileName->c_str());
 186          }
 187
 188        // Stage 2.3: when the file contains a volume (as opposed to an image),
 189        // notify the caller.
 190        if (NZ > 1)
 191          {
 192          vtkErrorMacro("This file contains multiple planes (images)"
 193                        << FileName->c_str());
 194          vtkErrorMacro("Removing this file from readed files "
 195                        << FileName->c_str());
 196          }
 197
 198        // Eventually, this file can be added on the stack. Update the
 199        // full size of the stack
 200        vtkDebugMacro("Number of planes added to the stack: " << NZ);
 201        ReturnedTotalNumberOfPlanes += NZ;
 202        continue;
 203
 204        } else {
 205        // We didn't have a workable reference file yet. Set this one
 206        // as the reference.
 207        FoundReferenceFile = true;
 208        vtkDebugMacro("This file taken as coherence reference:"
 209                      << FileName->c_str());
 210        vtkDebugMacro("Image dimension of reference file as read from Gdcm:" <<
 211                      NX << " " << NY << " " << NZ);
 212        vtkDebugMacro("Number of planes added to the stack: " << NZ);
 213        // Set aside the size of the image
 214        this->NumColumns = NX;
 215        this->NumLines   = NY;
 216        ReferenceNZ      = NZ;
 217        ReturnedTotalNumberOfPlanes += NZ;
 218        this->ImageType = type;
 219        this->PixelSize = GdcmHeader.GetPixelSize();
 220        }
 221      } // End of loop on FileName
 222
 223    ///////// The files we CANNOT load are flaged. On debugging purposes
 224    // count the loadable number of files and display thir number:
 225    int NumberCoherentFiles = 0;
 226    for (std::list<std::string>::iterator Filename  = FileNameList.begin();
 227                                         Filename != FileNameList.end();
 228                                       ++Filename)
 229      {
 230      if (*Filename != "GDCM_UNREADABLE")
 231         NumberCoherentFiles++;
 232      }
 233    vtkDebugMacro("Number of coherent files: " << NumberCoherentFiles);
 234
 235    if (ReturnedTotalNumberOfPlanes == 0)
 236      {
 237      vtkErrorMacro("No loadable file.");
 238      }
 239
 240    vtkDebugMacro("Total number of planes on the stack: "
 241                  << ReturnedTotalNumberOfPlanes);
 242    return ReturnedTotalNumberOfPlanes;
 243 }
 244
 245 //----------------------------------------------------------------------------
 246 // Configure the output e.g. WholeExtent, spacing, origin, scalar type...
 247 void vtkGdcmReader::ExecuteInformation()
 248 {
 249   //FIXME free any old memory
 250   this->TotalNumberOfPlanes = this->CheckFileCoherence();
 251   if ( this->TotalNumberOfPlanes == 0)
 252     {
 253        vtkErrorMacro("File set is not coherent. Exiting...");
 254        return;
 255     }
 256
 257   // if the user has not set the extent, but has set the VOI
 258   // set the zaxis extent to the VOI z axis
 259   if (this->DataExtent[4]==0 && this->DataExtent[5] == 0 &&
 260      (this->DataVOI[4] || this->DataVOI[5]))
 261     {
 262     this->DataExtent[4] = this->DataVOI[4];
 263     this->DataExtent[5] = this->DataVOI[5];
 264     }
 265
 266   // When the user has set the VOI, check it's coherence with the file content.
 267   if (this->DataVOI[0] || this->DataVOI[1] ||
 268       this->DataVOI[2] || this->DataVOI[3] ||
 269       this->DataVOI[4] || this->DataVOI[5])
 270     {
 271     if ((this->DataVOI[0] < 0) ||
 272         (this->DataVOI[1] >= this->NumColumns) ||
 273         (this->DataVOI[2] < 0) ||
 274         (this->DataVOI[3] >= this->NumLines) ||
 275         (this->DataVOI[4] < 0) ||
 276         (this->DataVOI[5] >= this->TotalNumberOfPlanes ))
 277       {
 278       vtkWarningMacro("The requested VOI is larger than expected extent.");
 279       this->DataVOI[0] = 0;
 280       this->DataVOI[1] = this->NumColumns - 1;
 281       this->DataVOI[2] = 0;
 282       this->DataVOI[3] = this->NumLines - 1;
 283       this->DataVOI[4] = 0;
 284       this->DataVOI[5] = this->TotalNumberOfPlanes - 1;
 285       }
 286     }
 287
 288   // Positionate the Extent.
 289   this->DataExtent[0] = 0;
 290   this->DataExtent[1] = this->NumColumns - 1;
 291   this->DataExtent[2] = 0;
 292   this->DataExtent[3] = this->NumLines - 1;
 293   if(this->FileNameList.size() > 1)
 294     {
 295     this->DataExtent[4] = 0;
 296     this->DataExtent[5] = this->TotalNumberOfPlanes - 1;
 297     }
 298
 299   // We don't need to positionate the Endian related stuff (by using
 300   // this->SetDataByteOrderToBigEndian() or SetDataByteOrderToLittleEndian()
 301   // since the reading of the file is done by gdcm.
 302   // But we do need to set up the data type for downstream filters:
 303   if      ( ImageType == "8U" )
 304     {
 305     vtkDebugMacro("8 bits unsigned image");
 306     this->SetDataScalarTypeToUnsignedChar();
 307     }
 308   else if ( ImageType == "8S" )
 309     {
 310     vtkErrorMacro("Cannot handle 8 bit signed files");
 311     return;
 312     }
 313   else if ( ImageType == "16U" )
 314     {
 315     vtkDebugMacro("16 bits unsigned image");
 316     this->SetDataScalarTypeToUnsignedShort();
 317     }
 318   else if ( ImageType == "16S" )
 319     {
 320     vtkDebugMacro("16 bits signed image");
 321     this->SetDataScalarTypeToShort();
 322     //vtkErrorMacro("Cannot handle 16 bit signed files");
 323     }
 324   else if ( ImageType == "32U" )
 325     {
 326     vtkDebugMacro("32 bits unsigned image");
 327     vtkDebugMacro("WARNING: forced to signed int !");
 328     this->SetDataScalarTypeToInt();
 329     }
 330   else if ( ImageType == "32S" )
 331     {
 332     vtkDebugMacro("32 bits signed image");
 333     this->SetDataScalarTypeToInt();
 334     }
 335
 336   vtkImageReader::ExecuteInformation();
 337 }
 338
 339 //----------------------------------------------------------------------------
 340 // Loads the contents of the image/volume contained by Filename at
 341 // the Dest memory address. Returns the size of the data loaded.
 342 size_t vtkGdcmReader::LoadImageInMemory(
 343              std::string FileName,
 344              unsigned char * Dest,
 345              const unsigned long UpdateProgressTarget,
 346              unsigned long & UpdateProgressCount)
 347 {
 348   vtkDebugMacro("Copying to memmory image" << FileName.c_str());
 349   gdcmFile GdcmFile(FileName.c_str());
 350   size_t size = GdcmFile.GetImageDataSize();
 351
 352   // If the data structure of vtk for image/volume representation
 353   // were straigthforwards the following would suffice:
 354   //    GdcmFile.GetImageDataIntoVector((void*)Dest, size);
 355   // But vtk chose to invert the lines of an image, that is the last
 356   // line comes first (for some axis related reasons?). Hence we need
 357   // to load the image line by line, starting from the end:
 358   int NumColumns = GdcmFile.GetXSize();
 359   int NumLines   = GdcmFile.GetYSize();
 360   int NumPlanes  = GdcmFile.GetZSize();
 361   int LineSize   = NumColumns * GdcmFile.GetPixelSize();
 362   unsigned char * Source      = (unsigned char*)GdcmFile.GetImageData();
 363   unsigned char * Destination = Dest + size - LineSize;
 364
 365   for (int plane = 0; plane < NumPlanes; plane++)
 366     {
 367     for (int line = 0; line < NumLines; line++)
 368       {
 369       // Copy one line at proper destination:
 370       memcpy((void*)Destination, (void*)Source, LineSize);
 371       Source      += LineSize;
 372       Destination -= LineSize;
 373       // Update progress related:
 374       if (!(UpdateProgressCount%UpdateProgressTarget))
 375         {
 376         this->UpdateProgress(UpdateProgressCount/(50.0*UpdateProgressTarget));
 377         }
 378       UpdateProgressCount++;
 379       }
 380     }
 381   return size;
 382 }
 383
 384 //----------------------------------------------------------------------------
 385 // Update -> UpdateData -> Execute -> ExecuteData (see vtkSource.cxx for
 386 // last step.
 387 // This function (redefinition of vtkImageReader::ExecuteData, see
 388 // VTK/IO/vtkImageReader.cxx) reads a data from a file. The datas
 389 // extent/axes are assumed to be the
 390 // same as the file extent/order.
 391 void vtkGdcmReader::ExecuteData(vtkDataObject *output)
 392 {
 393   if (this->FileNameList.empty())
 394     {
 395     vtkErrorMacro("A least a valid FileName must be specified.");
 396     return;
 397     }
 398
 399   vtkImageData *data = this->AllocateOutputData(output);
 400   data->SetExtent(this->DataExtent);
 401   data->GetPointData()->GetScalars()->SetName("DicomImage-Volume");
 402
 403   // The memory size for a full stack of images of course depends
 404   // on the number of planes and the size of each image:
 405   size_t StackNumPixels = this->NumColumns * this->NumLines
 406                         * this->TotalNumberOfPlanes;
 407   size_t stack_size = StackNumPixels * this->PixelSize;
 408   // Allocate pixel data space itself.
 409   unsigned char *mem = new unsigned char [stack_size];
 410
 411   // Variables for the UpdateProgress. We shall use 50 steps to signify
 412   // the advance of the process:
 413   unsigned long UpdateProgressTarget = (unsigned long) ceil (this->NumLines
 414                                      * this->TotalNumberOfPlanes
 415                                      / 50.0);
 416   // The actual advance measure:
 417   unsigned long UpdateProgressCount = 0;
 418
 419   // Feeling the allocated memory space with each image/volume:
 420   unsigned char * Dest = mem;
 421   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 422                                         FileName != FileNameList.end();
 423                                       ++FileName)
 424     {
 425     // Images that were tagged as unreadable in CheckFileCoherence()
 426     // are substituted with a black image to let the caller visually
 427     // notice something wrong is going on:
 428     if (*FileName != "GDCM_UNREADABLE")
 429       {
 430       Dest += this->LoadImageInMemory(*FileName, Dest,
 431                                       UpdateProgressTarget,
 432                                       UpdateProgressCount);
 433       } else {
 434       // We insert a black image in the stack for the user to be aware that
 435       // this image/volume couldn't be loaded. We simply skip one image
 436       // size:
 437       Dest += this->NumColumns * this->NumLines * this->PixelSize;
 438       // Update progress related:
 439       UpdateProgressCount += this->NumLines;
 440       if (!(UpdateProgressCount%UpdateProgressTarget))
 441         {
 442         this->UpdateProgress(UpdateProgressCount/(50.0*UpdateProgressTarget));
 443         }
 444       } // Else, file not loadable
 445     } // Loop on files
 446
 447   // The "size" of the vtkScalars data is expressed in number of points,
 448   // and is not the memory size representing those points:
 449   data->GetPointData()->GetScalars()->SetVoidArray(mem, StackNumPixels, 0);
 450   this->Modified();
 451 }
 452
 453 //----------------------------------------------------------------------------
 454 void vtkGdcmReader::PrintSelf(ostream& os, vtkIndent indent)
 455 {
 456   vtkImageReader::PrintSelf(os,indent);
 457   os << indent << "Filenames  : " << endl;
 458   vtkIndent nextIndent = indent.GetNextIndent();
 459   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 460                                         FileName != FileNameList.end();
 461                                       ++FileName)
 462     {
 463     os << nextIndent << FileName->c_str() << endl ;
 464     }
 465 }