vtk/vtkGdcmReader.cxx

   1 // $Header: /cvs/public/gdcm/vtk/vtkGdcmReader.cxx,v 1.9 2003/06/03 10:26:07 frog Exp $
   2 //CLEANME#include <vtkByteSwap.h>
   3 #include <stdio.h>
   4 #include <vtkObjectFactory.h>
   5 #include <vtkImageData.h>
   6 #include <vtkPointData.h>
   7 #include "vtkGdcmReader.h"
   8 #include "gdcm.h"
   9
  10 vtkGdcmReader::vtkGdcmReader()
  11 {
  12   // Constructor
  13 }
  14
  15 //----------------------------------------------------------------------------
  16 vtkGdcmReader::~vtkGdcmReader()
  17 {
  18   // FIXME free memory
  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    // Starting from two files we have a stack of images:
  31    if(this->FileNameList.size() >= 2)
  32       this->SetFileDimensionality(3);
  33   this->Modified();
  34 }
  35
  36 //----------------------------------------------------------------------------
  37 // Sets up a filename to be read.
  38 void vtkGdcmReader::SetFileName(const char *name) {
  39    vtkImageReader2::SetFileName(name);
  40    // Since we maintain a list of filenames (when building a volume)
  41    // we additionaly need to maintain this list. First we clean-up the
  42    // list and then positionate the incoming filename:
  43    this->FileNameList.empty();
  44    this->AddFileName(name);
  45 }
  46
  47 //----------------------------------------------------------------------------
  48 // vtkGdcmReader can have the file names specified through two ways:
  49 // (1) by calling the vtkImageReader2::SetFileName(), SetFilePrefix() and
  50 //     SetFilePattern()
  51 // (2) By successive calls to vtkGdcmReader::SetFileName()
  52 // When the first method was used by caller we need to update the local
  53 // filename list
  54 void vtkGdcmReader::BuilFileListFromPattern()
  55 {
  56    if (! this->FileNameList.empty())
  57       return;
  58    if (!this->FileName && !this->FilePattern)
  59      {
  60      vtkErrorMacro("FileNames are not set. Either use AddFileName() or");
  61      vtkErrorMacro("specify a FileName or FilePattern.");
  62      return;
  63      }
  64    for (int idx = this->DataExtent[4]; idx <= this->DataExtent[5]; ++idx)
  65      {
  66      this->ComputeInternalFileName(idx);
  67      vtkDebugMacro("Adding file " << this->InternalFileName);
  68      this->AddFileName(this->InternalFileName);
  69      }
  70 }
  71
  72 //----------------------------------------------------------------------------
  73 // When more than one filename is specified (i.e. we expect loading
  74 // a stack or volume) we need to check the corresponding images are
  75 // coherent:
  76 //  - they all share the same X dimensions
  77 //  - they all share the same Y dimensions
  78 //  - each file a Z dimension of 1
  79 //  - they all share the same type ( 8 bit signed, or unsigned...)
  80 bool vtkGdcmReader::CheckFileCoherence()
  81 {
  82    this->BuilFileListFromPattern();
  83    if (this->FileNameList.empty())
  84      {
  85      vtkErrorMacro("FileNames are not set.");
  86      return false;
  87      }
  88    if (this->FileNameList.size() == 1)
  89      {
  90      vtkDebugMacro("Single file specified.");
  91      return true;
  92      }
  93
  94    // Loop on the filenames:
  95    // - check for their existence and gdcm "parasability"
  96    // - get the coherence check done:
  97    bool FoundReferenceFile = false;
  98    int ReferenceNX;
  99    int ReferenceNY;
 100    int ReferenceNZ;
 101    std::string ReferenceType;
 102    for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 103                                         FileName != FileNameList.end();
 104                                       ++FileName)
 105      {
 106      // Check for file existence.
 107      FILE *fp;
 108      fp = fopen(FileName->c_str(),"rb");
 109      if (!fp)
 110        {
 111        vtkErrorMacro("Unable to open file " << *FileName->c_str());
 112        vtkErrorMacro("Removing this file from readed files "
 113                      << *FileName->c_str());
 114        FileNameList.remove(*FileName);
 115        continue;
 116        }
 117      fclose(fp);
 118
 119      // Check for Gdcm parsability
 120      gdcmHeader GdcmHeader(FileName->c_str());
 121      if (!GdcmHeader.IsReadable())
 122        {
 123        vtkErrorMacro("Gdcm cannot parse file " << *FileName->c_str());
 124        vtkErrorMacro("Removing this file from readed files "
 125                      << *FileName->c_str());
 126        FileNameList.remove(*FileName);
 127        continue;
 128        }
 129
 130      // We don't know how to handle multiple images in one file yet:
 131      int NZ = GdcmHeader.GetZSize();
 132      if (NZ > 1)
 133        {
 134        vtkErrorMacro("This file contains multiple planes (images)"
 135                      << *FileName->c_str());
 136        vtkErrorMacro("Removing this file from readed files "
 137                      << *FileName->c_str());
 138        FileNameList.remove(*FileName);
 139        continue;
 140        }
 141
 142      // Coherence stage:
 143      int NX = GdcmHeader.GetXSize();
 144      int NY = GdcmHeader.GetYSize();
 145      std::string type = GdcmHeader.GetPixelType();
 146      if (FoundReferenceFile)
 147        {
 148        if (   ( NX != ReferenceNX )
 149            || ( NY != ReferenceNY )
 150            || ( NZ != ReferenceNZ )
 151            || ( type != ReferenceType ) )
 152          {
 153             vtkErrorMacro("This file is not coherent with previous ones"
 154                           << *FileName->c_str());
 155             vtkErrorMacro("Removing this file from readed files "
 156                           << *FileName->c_str());
 157             FileNameList.remove(*FileName);
 158             continue;
 159          } else {
 160             vtkDebugMacro("File is coherent with previous ones"
 161                           << *FileName->c_str());
 162          }
 163        } else {
 164          // This file shall be the reference:
 165          FoundReferenceFile = true;
 166          ReferenceNX = NX;
 167          ReferenceNY = NY;
 168          ReferenceNZ = NZ;
 169          ReferenceType = type;
 170          vtkDebugMacro("This file taken as coherence reference:"
 171                        << *FileName->c_str());
 172        }
 173      } // End of loop on FileName
 174
 175    if (this->FileNameList.empty())
 176      {
 177      vtkDebugMacro("No gdcm parsable file.");
 178      return false;
 179      }
 180    if (this->FileNameList.size() == 1)
 181      {
 182      vtkDebugMacro("Single parsable file left after coherence test.");
 183      return true;
 184      }
 185    return true;
 186 }
 187
 188 //----------------------------------------------------------------------------
 189 // Configure the output e.g. WholeExtent, spacing, origin, scalar type...
 190 void vtkGdcmReader::ExecuteInformation()
 191 {
 192   //FIXME free any old memory
 193
 194   // if the user has not set the extent, but has set the VOI
 195   // set the zaxis extent to the VOI z axis
 196   if (this->DataExtent[4]==0 && this->DataExtent[5] == 0 &&
 197      (this->DataVOI[4] || this->DataVOI[5]))
 198     {
 199     this->DataExtent[4] = this->DataVOI[4];
 200     this->DataExtent[5] = this->DataVOI[5];
 201     }
 202   if ( ! this->CheckFileCoherence() )
 203     {
 204        vtkErrorMacro("File set is not coherent. Exiting...");
 205        return;
 206     }
 207   std::string ReferenceFile = this->FileNameList.front();
 208   gdcmHeader GdcmHeader(ReferenceFile.c_str());
 209
 210   int NX = GdcmHeader.GetXSize();
 211   int NY = GdcmHeader.GetYSize();
 212   int NZ = GdcmHeader.GetZSize();
 213   vtkDebugMacro("Image dimension as read from Gdcm:" <<
 214                 NX << " " << NY << " " << NZ);
 215
 216   // When the user has set the VOI, check it's coherence with the file content.
 217   if (this->DataVOI[0] || this->DataVOI[1] ||
 218       this->DataVOI[2] || this->DataVOI[3] ||
 219       this->DataVOI[4] || this->DataVOI[5])
 220     {
 221     if ((this->DataVOI[0] < 0) ||
 222         (this->DataVOI[1] >= NX) ||
 223         (this->DataVOI[2] < 0) ||
 224         (this->DataVOI[3] >= NY) ||
 225         (this->DataVOI[4] < 0) ||
 226         (this->DataVOI[5] >= this->FileNameList.size()))
 227       {
 228       vtkWarningMacro("The requested VOI is larger than the file's ("
 229                       << ReferenceFile.c_str() << ") extent ");
 230       this->DataVOI[0] = 0;
 231       this->DataVOI[1] = NX - 1;
 232       this->DataVOI[2] = 0;
 233       this->DataVOI[3] = NY - 1;
 234       this->DataVOI[4] = 0;
 235       this->DataVOI[5] = this->FileNameList.size() - 1;
 236       }
 237     }
 238
 239   // Positionate the Extent.
 240   this->DataExtent[0] = 0;
 241   this->DataExtent[1] = NX - 1;
 242   this->DataExtent[2] = 0;
 243   this->DataExtent[3] = NY - 1;
 244   if(this->FileNameList.size() > 1)
 245     {
 246     this->DataExtent[4] = 0;
 247     this->DataExtent[5] = this->FileNameList.size() - 1;
 248     }
 249
 250
 251   // We don't need to positionate the Endian related stuff (by using
 252   // this->SetDataByteOrderToBigEndian() or SetDataByteOrderToLittleEndian()
 253   // since the reading of the file is done by gdcm.
 254   // But we do need to set up the data type for downstream filters:
 255   std::string type = GdcmHeader.GetPixelType();
 256   if      ( type == "8U" )
 257     {
 258     vtkDebugMacro("8 bits unsigned image");
 259     this->SetDataScalarTypeToUnsignedChar();
 260     }
 261   else if ( type == "8S" )
 262     {
 263     vtkErrorMacro("Cannot handle 8 bit signed files");
 264     return;
 265     }
 266   else if ( type == "16U" )
 267     {
 268     vtkDebugMacro("16 bits unsigned image");
 269     this->SetDataScalarTypeToUnsignedShort();
 270     }
 271   else if ( type == "16S" )
 272     {
 273     vtkDebugMacro("16 bits signed image");
 274     this->SetDataScalarTypeToShort();
 275     //vtkErrorMacro("Cannot handle 16 bit signed files");
 276     }
 277   else if ( type == "32U" )
 278     {
 279     vtkDebugMacro("32 bits unsigned image");
 280     vtkDebugMacro("WARNING: forced to signed int !");
 281     this->SetDataScalarTypeToInt();
 282     }
 283   else if ( type == "32S" )
 284     {
 285     vtkDebugMacro("32 bits signed image");
 286     this->SetDataScalarTypeToInt();
 287     }
 288   else
 289     {
 290     vtkErrorMacro("Bad File Type " << ReferenceFile.c_str()
 291                                    << "Type " << type.c_str());
 292     return;
 293     }
 294
 295   vtkImageReader::ExecuteInformation();
 296 }
 297
 298 //----------------------------------------------------------------------------
 299 void vtkGdcmReader::LoadImageInMemory(std::string FileName,
 300                                       unsigned char * Dest,
 301                                       size_t size)
 302 {
 303   vtkDebugMacro("Copying to memmory image" << FileName.c_str());
 304   gdcmFile GdcmFile(FileName.c_str());
 305
 306   if (GdcmFile.GetZSize() != 1 )
 307     vtkErrorMacro("Cannot handle images with multiple planes");
 308
 309   // First check the expected size of the image is the one found by gdcm.
 310   if ( size != GdcmFile.GetImageDataSize() )
 311     {
 312     vtkErrorMacro("Inconsistency with GetImageDataSize for file"
 313                   << FileName.c_str());
 314     vtkErrorMacro("Number of scalar components"
 315                   << this->NumberOfScalarComponents);
 316     }
 317
 318   // If the data structure of vtk for image/volume representation
 319   // were straigthforwards the following would suffice:
 320   //    GdcmFile.GetImageDataIntoVector((void*)Dest, size);
 321   // But vtk chose to invert the lines of an image, that is the last
 322   // line comes first (for some axis related reasons?). Hence we need
 323   // to load the image line by line, starting from the end:
 324   int NumColumns = GdcmFile.GetXSize();
 325   int NumLines   = GdcmFile.GetYSize();
 326   int LineSize   = NumColumns * GdcmFile.GetPixelSize();
 327   unsigned char * Source      = (unsigned char*)GdcmFile.GetImageData();
 328   unsigned char * Destination = Dest + size - LineSize;
 329   for (int i = 0; i < NumLines; i++)
 330     {
 331     memcpy((void*)Destination, (void*)Source, LineSize);
 332     Source      += LineSize;
 333     Destination -= LineSize;
 334     }
 335 }
 336
 337 //----------------------------------------------------------------------------
 338 // Update -> UpdateData -> Execute -> ExecuteData (see vtkSource.cxx for
 339 // last step.
 340 // This function (redefinition of vtkImageReader::ExecuteData, see
 341 // VTK/IO/vtkImageReader.cxx) reads a data from a file. The datas
 342 // extent/axes are assumed to be the
 343 // same as the file extent/order.
 344 void vtkGdcmReader::ExecuteData(vtkDataObject *output)
 345 {
 346   if (this->FileNameList.empty())
 347     {
 348     vtkErrorMacro("A least a valid FileName must be specified.");
 349     return;
 350     }
 351
 352   vtkImageData *data = this->AllocateOutputData(output);
 353   data->SetExtent(this->DataExtent);
 354   data->GetPointData()->GetScalars()->SetName("DicomImage-Volume");
 355
 356   // First check the coherence between the DataExtent and the
 357   // size of the pixel data as annouced by gdcm (looks a bit paranoid)
 358   // for the reference file (i.e. the first one in the list):
 359   std::string ReferenceFile = this->FileNameList.front();
 360   gdcmFile GdcmFile(ReferenceFile.c_str());
 361   int NumColumns = this->DataExtent[1] - this->DataExtent[0] + 1;
 362   int NumLines   = this->DataExtent[3] - this->DataExtent[2] + 1;
 363   int NumPlanes  = 1;   // This has been checked in CheckFileCoherence
 364   size_t size = NumColumns * NumLines * NumPlanes * GdcmFile.GetPixelSize();
 365   if ( size != GdcmFile.GetImageDataSize() )
 366     {
 367     vtkDebugMacro("Inconsistency with GetImageDataSize");
 368     vtkDebugMacro("Number of scalar components"
 369                   << this->NumberOfScalarComponents);
 370     }
 371
 372   // The memory size for a full stack of images of course depends
 373   // on the number of images:
 374   size_t stack_size = size * this->FileNameList.size();
 375   // Allocate pixel data space itself.
 376   unsigned char *mem = new unsigned char [stack_size];
 377
 378   unsigned char * Dest = mem;
 379   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 380                                         FileName != FileNameList.end();
 381                                       ++FileName)
 382     {
 383        this->LoadImageInMemory(*FileName, Dest, size);
 384        Dest += size;
 385     }
 386
 387   // The "size" of the vtkScalars data is expressed in number of points,
 388   // and is not the memory size representing those points:
 389   stack_size = stack_size / GdcmFile.GetPixelSize();
 390   data->GetPointData()->GetScalars()->SetVoidArray(mem, stack_size, 0);
 391   this->Modified();
 392 }
 393
 394 //----------------------------------------------------------------------------
 395 void vtkGdcmReader::PrintSelf(ostream& os, vtkIndent indent)
 396 {
 397   vtkImageReader::PrintSelf(os,indent);
 398   os << indent << "Filenames  : " << endl;
 399   vtkIndent nextIndent = indent.GetNextIndent();
 400   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 401                                         FileName != FileNameList.end();
 402                                       ++FileName)
 403     {
 404     os << nextIndent << *FileName->c_str() << endl ;
 405     }
 406 }