vtk/vtkGdcmReader.cxx

   1 // $Header: /cvs/public/gdcm/vtk/vtkGdcmReader.cxx,v 1.7 2003/05/30 18:48:36 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);
 112        vtkErrorMacro("Removing this file from readed files " << *FileName);
 113        FileNameList.remove(*FileName);
 114        continue;
 115        }
 116      fclose(fp);
 117
 118      // Check for Gdcm parsability
 119      gdcmHeader GdcmHeader(FileName->c_str());
 120      if (!GdcmHeader.IsReadable())
 121        {
 122        vtkErrorMacro("Gdcm cannot parse file " << *FileName);
 123        vtkErrorMacro("Removing this file from readed files " << *FileName);
 124        FileNameList.remove(*FileName);
 125        continue;
 126        }
 127
 128      // Coherence stage:
 129      int NX = GdcmHeader.GetXSize();
 130      int NY = GdcmHeader.GetYSize();
 131      int NZ = GdcmHeader.GetZSize();
 132      std::string type = GdcmHeader.GetPixelType();
 133      if (FoundReferenceFile)
 134        {
 135        if (   ( NX != ReferenceNX )
 136            || ( NY != ReferenceNY )
 137            || ( NZ != ReferenceNZ )
 138            || ( type != ReferenceType ) )
 139          {
 140             vtkErrorMacro("This file is not coherent with previous ones"
 141                           << *FileName);
 142             vtkErrorMacro("Removing this file from readed files " << *FileName);
 143             FileNameList.remove(*FileName);
 144             continue;
 145          } else {
 146             vtkDebugMacro("File is coherent with previous ones" << *FileName);
 147          }
 148        } else {
 149          // This file shall be the reference:
 150          FoundReferenceFile = true;
 151          ReferenceNX = NX;
 152          ReferenceNY = NY;
 153          ReferenceNZ = NZ;
 154          ReferenceType = type;
 155          vtkDebugMacro("This file taken as coherence reference:" << *FileName);
 156        }
 157      } // End of loop on FileName
 158
 159    if (this->FileNameList.empty())
 160      {
 161      vtkDebugMacro("No gdcm parsable file.");
 162      return false;
 163      }
 164    if (this->FileNameList.size() == 1)
 165      {
 166      vtkDebugMacro("Single parsable file left after coherence test.");
 167      return true;
 168      }
 169    return true;
 170 }
 171
 172 //----------------------------------------------------------------------------
 173 // Configure the output e.g. WholeExtent, spacing, origin, scalar type...
 174 void vtkGdcmReader::ExecuteInformation()
 175 {
 176   //FIXME free any old memory
 177
 178   // if the user has not set the extent, but has set the VOI
 179   // set the zaxis extent to the VOI z axis
 180   if (this->DataExtent[4]==0 && this->DataExtent[5] == 0 &&
 181      (this->DataVOI[4] || this->DataVOI[5]))
 182     {
 183     this->DataExtent[4] = this->DataVOI[4];
 184     this->DataExtent[5] = this->DataVOI[5];
 185     }
 186   if ( ! this->CheckFileCoherence() )
 187     {
 188        vtkErrorMacro("File set is not coherent. Exiting...");
 189        return;
 190     }
 191   string ReferenceFile = this->FileNameList.front();
 192   gdcmHeader GdcmHeader(ReferenceFile.c_str());
 193
 194   int NX = GdcmHeader.GetXSize();
 195   int NY = GdcmHeader.GetYSize();
 196   int NZ = GdcmHeader.GetZSize();
 197   vtkDebugMacro("Image dimension as read from Gdcm:" <<
 198                 NX << " " << NY << " " << NZ);
 199
 200   if(NZ>1) this->SetFileDimensionality(3);
 201
 202   // When the user has set the VOI, check it's coherence with the file content.
 203   if (this->DataVOI[0] || this->DataVOI[1] ||
 204       this->DataVOI[2] || this->DataVOI[3] ||
 205       this->DataVOI[4] || this->DataVOI[5])
 206     {
 207     if ((this->DataVOI[0] < 0) ||
 208         (this->DataVOI[1] >= NX) ||
 209         (this->DataVOI[2] < 0) ||
 210         (this->DataVOI[3] >= NY) ||
 211         (this->DataVOI[4] < 0) ||
 212         (this->DataVOI[5] >= NZ))
 213       {
 214       vtkWarningMacro("The requested VOI is larger than the file's ("
 215                       << ReferenceFile << ") extent ");
 216       this->DataVOI[0] = 0;
 217       this->DataVOI[1] = NX - 1;
 218       this->DataVOI[2] = 0;
 219       this->DataVOI[3] = NY - 1;
 220       this->DataVOI[4] = 0;
 221       this->DataVOI[5] = NZ - 1;
 222       }
 223     }
 224
 225   // Positionate the Extent.
 226   this->DataExtent[0] = 0;
 227   this->DataExtent[1] = NX - 1;
 228   this->DataExtent[2] = 0;
 229   this->DataExtent[3] = NY - 1;
 230   if(this->GetFileDimensionality()==3)
 231     {
 232       this->DataExtent[4] = 0;
 233       this->DataExtent[5] = NZ - 1;
 234     }
 235
 236   // We don't need to positionate the Endian related stuff (by using
 237   // this->SetDataByteOrderToBigEndian() or SetDataByteOrderToLittleEndian()
 238   // since the reading of the file is done by gdcm.
 239   // But we do need to set up the data type for downstream filters:
 240   std::string type = GdcmHeader.GetPixelType();
 241   if      ( type == "8U" )
 242     {
 243     vtkDebugMacro("8 bits unsigned image");
 244     this->SetDataScalarTypeToUnsignedChar();
 245     }
 246   else if ( type == "8S" )
 247     {
 248     vtkErrorMacro("Cannot handle 8 bit signed files");
 249     return;
 250     }
 251   else if ( type == "16U" )
 252     {
 253     vtkDebugMacro("16 bits unsigned image");
 254     this->SetDataScalarTypeToUnsignedShort();
 255     }
 256   else if ( type == "16S" )
 257     {
 258     vtkDebugMacro("16 bits signed image");
 259     this->SetDataScalarTypeToShort();
 260     //vtkErrorMacro("Cannot handle 16 bit signed files");
 261     }
 262   else if ( type == "32U" )
 263     {
 264     vtkDebugMacro("32 bits unsigned image");
 265     vtkDebugMacro("WARNING: forced to signed int !");
 266     this->SetDataScalarTypeToInt();
 267     }
 268   else if ( type == "32S" )
 269     {
 270     vtkDebugMacro("32 bits signed image");
 271     this->SetDataScalarTypeToInt();
 272     }
 273   else
 274     {
 275     vtkErrorMacro("Bad File Type " << ReferenceFile
 276                                    << "Type " << type.c_str());
 277     return;
 278     }
 279
 280   vtkImageReader::ExecuteInformation();
 281 }
 282
 283 //----------------------------------------------------------------------------
 284 void vtkGdcmReader::LoadImageInMemory(string FileName,
 285                                       unsigned char * Dest,
 286                                       size_t size)
 287 {
 288   vtkDebugMacro("Copying to memmory image" << FileName.c_str());
 289   gdcmFile GdcmFile(FileName.c_str());
 290
 291   if (GdcmFile.GetZSize() != 1 )
 292     vtkErrorMacro("Cannot handle images with multiple planes");
 293
 294   // First check the expected size of the image is the one found by gdcm.
 295   if ( size != GdcmFile.GetImageDataSize() )
 296     {
 297     vtkErrorMacro("Inconsistency with GetImageDataSize for file"
 298                   << FileName.c_str());
 299     vtkErrorMacro("Number of scalar components"
 300                   << this->NumberOfScalarComponents);
 301     }
 302
 303   // If the data structure of vtk for image/volume representation
 304   // were straigthforwards the following would suffice:
 305   //    GdcmFile.GetImageDataIntoVector((void*)Dest, size);
 306   // But vtk chose to invert the lines of an image, that is the last
 307   // line comes first (for some axis related reasons?). Hence we need
 308   // to load the image line by line, starting from the end:
 309   int NumColumns = GdcmFile.GetXSize();
 310   int NumLines   = GdcmFile.GetYSize();
 311   int LineSize   = NumColumns * GdcmFile.GetPixelSize();
 312   unsigned char * Source      = (unsigned char*)GdcmFile.GetImageData();
 313   unsigned char * Destination = Dest + size - LineSize;
 314   for (int i = 0; i < NumLines; i++)
 315     {
 316     memcpy((void*)Destination, (void*)Source, LineSize);
 317     Source      += LineSize;
 318     Destination -= LineSize;
 319     }
 320 }
 321
 322 //----------------------------------------------------------------------------
 323 // Update -> UpdateData -> Execute -> ExecuteData (see vtkSource.cxx for
 324 // last step.
 325 // This function (redefinition of vtkImageReader::ExecuteData, see
 326 // VTK/IO/vtkImageReader.cxx) reads a data from a file. The datas
 327 // extent/axes are assumed to be the
 328 // same as the file extent/order.
 329 void vtkGdcmReader::ExecuteData(vtkDataObject *output)
 330 {
 331   if (this->FileNameList.empty())
 332     {
 333     vtkErrorMacro("A least a valid FileName must be specified.");
 334     return;
 335     }
 336
 337   vtkImageData *data = this->AllocateOutputData(output);
 338   data->SetExtent(this->DataExtent);
 339   data->GetPointData()->GetScalars()->SetName("DicomImage-Volume");
 340
 341   // First check the coherence between the DataExtent and the
 342   // size of the pixel data as annouced by gdcm (looks a bit paranoid)
 343   // for the reference file (i.e. the first one in the list):
 344   string ReferenceFile = this->FileNameList.front();
 345   gdcmFile GdcmFile(ReferenceFile.c_str());
 346   int NumColumns = this->DataExtent[1] - this->DataExtent[0] + 1;
 347   int NumLines   = this->DataExtent[3] - this->DataExtent[2] + 1;
 348   int NumPlanes  = this->DataExtent[5] - this->DataExtent[4] + 1;
 349   size_t size = NumColumns * NumLines * NumPlanes * GdcmFile.GetPixelSize();
 350   if ( size != GdcmFile.GetImageDataSize() )
 351     {
 352     vtkDebugMacro("Inconsistency with GetImageDataSize");
 353     vtkDebugMacro("Number of scalar components"
 354                   << this->NumberOfScalarComponents);
 355     }
 356
 357   // The memory size for a full stack of images of course depends
 358   // on the number of images:
 359   size_t stack_size = size * this->FileNameList.size();
 360   // Allocate pixel data space itself.
 361   unsigned char *mem = new unsigned char [stack_size];
 362
 363   unsigned char * Dest = mem;
 364   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 365                                         FileName != FileNameList.end();
 366                                       ++FileName)
 367     {
 368        this->LoadImageInMemory(*FileName, Dest, size);
 369        Dest += size;
 370     }
 371
 372
 373   // The "size" of the vtkScalars data is expressed in number of points,
 374   // and is not the memory size representing those points:
 375   stack_size = stack_size / GdcmFile.GetPixelSize();
 376   data->GetPointData()->GetScalars()->SetVoidArray(mem, stack_size, 0);
 377   this->Modified();
 378
 379 }
 380
 381 //----------------------------------------------------------------------------
 382 void vtkGdcmReader::PrintSelf(ostream& os, vtkIndent indent)
 383 {
 384   vtkImageReader::PrintSelf(os,indent);
 385   os << indent << "Filenames  : " << endl;
 386   vtkIndent nextIndent = indent.GetNextIndent();
 387   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
 388                                         FileName != FileNameList.end();
 389                                       ++FileName)
 390     {
 391     os << nextIndent << *FileName << endl ;
 392     }
 393 }