]> Creatis software - gdcm.git/blobdiff - vtk/vtkGdcmReader.cxx
* Modification to compile with the distutils.
[gdcm.git] / vtk / vtkGdcmReader.cxx
index 4a18c5ac1cbd56a76ffe16987b3231e950434952..1e0852bdbab52c64923cf2d9c62321b20455b80e 100644 (file)
@@ -1,21 +1,76 @@
-// $Header: /cvs/public/gdcm/vtk/vtkGdcmReader.cxx,v 1.11 2003/06/12 16:58:31 malaterre Exp $
-//CLEANME#include <vtkByteSwap.h>
-#include <stdio.h>
+// $Header: /cvs/public/gdcm/vtk/vtkGdcmReader.cxx,v 1.29 2003/12/22 12:46:19 regrain Exp $
+// //////////////////////////////////////////////////////////////
+// WARNING TODO CLENAME 
+// Actual limitations of this code:
+//
+// /////// Redundant and unnecessary header parsing
+// In it's current state this code actually parses three times the Dicom
+// header of a file before the corresponding image gets loaded in the
+// ad-hoc vtkData !
+// Here is the process:
+//  1/ First loading happens in ExecuteInformation which in order to
+//     positionate the vtk extents calls CheckFileCoherence. The purpose
+//     of CheckFileCoherence is to make sure all the images in the future
+//     stack are "homogenous" (same size, same representation...). This
+//     can only be achieved by parsing all the Dicom headers...
+//  2/ ExecuteData is then responsible for the next two loadings:
+//  2a/ ExecuteData calls AllocateOutputData that in turn seems to 
+//      (indirectely call) ExecuteInformation which ends up in a second
+//      header parsing
+//  2b/ the core of ExecuteData then needs gdcmFile (which in turns
+//      initialises gdcmHeader in the constructor) in order to access
+//      the data-image.
+//
+// Possible solution:
+// maintain a list of gdcmFiles (created by say ExecuteInformation) created
+// once and for all accross the life of vtkGdcmHeader (it would only load
+// new gdcmFile if the user changes the list). ExecuteData would then use 
+// those gdcmFile and hence avoid calling the construtor:
+//  - advantage: the header of the files would only be parser once.
+//  - drawback: once execute information is called (i.e. on creation of
+//              a vtkGdcmHeader) the gdcmFile structure is loaded in memory.
+//              The average size of a gdcmHeader being of 100Ko, is one
+//              loads 10 stacks of images with say 200 images each, you
+//              end-up with a loss of 200Mo...
+//
+// /////// Never unallocated memory:
+// ExecuteData allocates space for the pixel data [which will get pointed
+// by the vtkPointData() through the call
+// data->GetPointData()->GetScalars()->SetVoidArray(mem, StackNumPixels, 0);]
+// This data is never "freed" neither in the destructor nor when the
+// filename list is extended, ExecuteData is called a second (or third)
+// time...
+// //////////////////////////////////////////////////////////////
+
+#include "gdcmFile.h"
+#include "gdcmHeaderHelper.h"
+#include "vtkGdcmReader.h"
+
+//#include <stdio.h>
 #include <vtkObjectFactory.h>
 #include <vtkImageData.h>
 #include <vtkPointData.h>
-#include "vtkGdcmReader.h"
-#include "gdcm.h"
+#include <vtkLookupTable.h>
+
 
 vtkGdcmReader::vtkGdcmReader()
 {
   // Constructor
+  this->LookupTable = NULL;
 }
 
 //----------------------------------------------------------------------------
 vtkGdcmReader::~vtkGdcmReader()
 { 
-  // FIXME free memory
+  this->RemoveAllFileName();
+  this->InternalFileNameList.clear();
+  if(this->LookupTable) this->LookupTable->Delete();
+}
+
+//----------------------------------------------------------------------------
+// Remove all files from the list of images to read.
+void vtkGdcmReader::RemoveAllFileName(void)
+{
   this->FileNameList.clear();
 }
 
@@ -38,11 +93,28 @@ void vtkGdcmReader::SetFileName(const char *name) {
   vtkImageReader2::SetFileName(name);
   // Since we maintain a list of filenames, when building a volume,
   // (see vtkGdcmReader::AddFileName), we additionaly need to purge
-  // this list when we manually positionate the filename:
+  // this list when we manually positionate the filename.
   this->FileNameList.clear();
   this->Modified();
 }
 
+//----------------------------------------------------------------------------
+// Adds a file name to the internal list of images to read.
+void vtkGdcmReader::RemoveAllInternalFileName(void)
+{
+  this->InternalFileNameList.clear();
+}
+
+//----------------------------------------------------------------------------
+// Adds a file name to the internal list of images to read.
+void vtkGdcmReader::AddInternalFileName(const char* name)
+{
+  char * LocalName = new char[strlen(name) + 1];
+  strcpy(LocalName, name);
+  this->InternalFileNameList.push_back(LocalName);
+  delete[] LocalName;
+}
+
 //----------------------------------------------------------------------------
 // vtkGdcmReader can have the file names specified through two ways:
 // (1) by calling the vtkImageReader2::SetFileName(), SetFilePrefix() and
@@ -62,6 +134,7 @@ void vtkGdcmReader::BuildFileListFromPattern()
    if (! this->FileNameList.empty()  )
      {
      vtkDebugMacro("Using the AddFileName specified files");
+         this->InternalFileNameList=this->FileNameList;
      return;
      }
 
@@ -72,12 +145,24 @@ void vtkGdcmReader::BuildFileListFromPattern()
      return;
      }
 
-   for (int idx = this->DataExtent[4]; idx <= this->DataExtent[5]; ++idx)
+   this->RemoveAllInternalFileName();
+   if( this->FileNameList.empty() )
      {
-     this->ComputeInternalFileName(idx);
+     //Multiframe case:
+     this->ComputeInternalFileName(this->DataExtent[4]);
      vtkDebugMacro("Adding file " << this->InternalFileName);
-     this->AddFileName(this->InternalFileName);
+     this->AddInternalFileName(this->InternalFileName);
      }
+   else
+     {
+     //stack of 2D dicom case:
+     for (int idx = this->DataExtent[4]; idx <= this->DataExtent[5]; ++idx)
+       {
+       this->ComputeInternalFileName(idx);
+       vtkDebugMacro("Adding file " << this->InternalFileName);
+       this->AddInternalFileName(this->InternalFileName);
+       }
+    }
 }
 
 //----------------------------------------------------------------------------
@@ -100,7 +185,7 @@ int vtkGdcmReader::CheckFileCoherence()
    int ReturnedTotalNumberOfPlanes = 0;   // The returned value.
 
    this->BuildFileListFromPattern();
-   if (this->FileNameList.empty())
+   if (this->InternalFileNameList.empty())
      {
      vtkErrorMacro("FileNames are not set.");
      return 0;
@@ -110,12 +195,23 @@ int vtkGdcmReader::CheckFileCoherence()
    int  ReferenceNZ = 0;
 
    // Loop on the filenames:
-   // - check for their existence and gdcm "parasability"
+   // - check for their existence and gdcm "parsability"
    // - get the coherence check done:
-   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
-                                        FileName != FileNameList.end();
+   for (std::list<std::string>::iterator FileName  = InternalFileNameList.begin();
+                                        FileName != InternalFileNameList.end();
                                       ++FileName)
      {
+     // The file is always added in the number of planes
+     //  - If file doesn't exist, it will be replaced by a black plane in the 
+     //    ExecuteData method
+     //  - If file has more than 1 plane, other planes will be added later to
+     //    to the ReturnedTotalNumberOfPlanes variable counter
+     ReturnedTotalNumberOfPlanes += 1;
+
+     /////// Stage 0: check for file name:
+         if(*FileName==std::string("GDCM_UNREADABLE"))
+                 continue;
+
      /////// Stage 1: check for file readability:
      // Stage 1.1: check for file existence.
      FILE *fp;
@@ -131,7 +227,7 @@ int vtkGdcmReader::CheckFileCoherence()
      fclose(fp);
    
      // Stage 1.2: check for Gdcm parsability
-     gdcmHeader GdcmHeader(FileName->c_str());
+     gdcmHeaderHelper GdcmHeader(FileName->c_str());
      if (!GdcmHeader.IsReadable())
        {
        vtkErrorMacro("Gdcm cannot parse file " << FileName->c_str());
@@ -148,6 +244,7 @@ int vtkGdcmReader::CheckFileCoherence()
          && (type != "32U") && (type != "32S") )
        {
        vtkErrorMacro("Bad File Type for file" << FileName->c_str());
+       vtkErrorMacro("                      " << type.c_str());
        vtkErrorMacro("Removing this file from readed files "
                      << FileName->c_str());
        *FileName = "GDCM_UNREADABLE";
@@ -192,14 +289,12 @@ int vtkGdcmReader::CheckFileCoherence()
          {
          vtkErrorMacro("This file contains multiple planes (images)"
                        << FileName->c_str());
-         vtkErrorMacro("Removing this file from readed files " 
-                       << FileName->c_str());
          }
 
        // Eventually, this file can be added on the stack. Update the
        // full size of the stack
        vtkDebugMacro("Number of planes added to the stack: " << NZ);
-       ReturnedTotalNumberOfPlanes += NZ;
+       ReturnedTotalNumberOfPlanes += NZ - 1; // First plane already added
        continue;
 
        } else {
@@ -215,20 +310,40 @@ int vtkGdcmReader::CheckFileCoherence()
        this->NumColumns = NX;
        this->NumLines   = NY;
        ReferenceNZ      = NZ;
-       ReturnedTotalNumberOfPlanes += NZ;
+       ReturnedTotalNumberOfPlanes += NZ - 1; // First plane already added
        this->ImageType = type;
        this->PixelSize = GdcmHeader.GetPixelSize();
+       
+       if( GdcmHeader.HasLUT() )
+         {
+         this->NumComponents = GdcmHeader.GetNumberOfScalarComponentsRaw();
+         }
+       else
+         {
+         this->NumComponents = GdcmHeader.GetNumberOfScalarComponents(); //rgb or mono
+         }
+       
+       //Set image spacing
+       this->DataSpacing[0] = GdcmHeader.GetXSpacing();
+       this->DataSpacing[1] = GdcmHeader.GetYSpacing();
+       this->DataSpacing[2] = GdcmHeader.GetZSpacing();
+
+       //Set image origin
+       this->DataOrigin[0] = GdcmHeader.GetXOrigin();
+       this->DataOrigin[1] = GdcmHeader.GetYOrigin();
+       this->DataOrigin[2] = GdcmHeader.GetZOrigin();
+       
        }
      } // End of loop on FileName
 
    ///////// The files we CANNOT load are flaged. On debugging purposes
-   // count the loadable number of files and display thir number:
+   // count the loadable number of files and display their number:
    int NumberCoherentFiles = 0;
-   for (std::list<std::string>::iterator FileName  = FileNameList.begin();
-                                        FileName != FileNameList.end();
-                                      ++FileName)
+   for (std::list<std::string>::iterator Filename  = InternalFileNameList.begin();
+                                        Filename != InternalFileNameList.end();
+                                      ++Filename)
      {
-     if (*FileName != "GDCM_UNREADABLE")
+     if (*Filename != "GDCM_UNREADABLE")
         NumberCoherentFiles++;    
      }
    vtkDebugMacro("Number of coherent files: " << NumberCoherentFiles);
@@ -240,14 +355,14 @@ int vtkGdcmReader::CheckFileCoherence()
 
    vtkDebugMacro("Total number of planes on the stack: "
                  << ReturnedTotalNumberOfPlanes);
-   return ReturnedTotalNumberOfPlanes;
+   
+       return ReturnedTotalNumberOfPlanes;
 }
 
 //----------------------------------------------------------------------------
 // Configure the output e.g. WholeExtent, spacing, origin, scalar type...
 void vtkGdcmReader::ExecuteInformation()
 {
-  //FIXME free any old memory
   this->TotalNumberOfPlanes = this->CheckFileCoherence();
   if ( this->TotalNumberOfPlanes == 0)
     {
@@ -256,7 +371,7 @@ void vtkGdcmReader::ExecuteInformation()
     }
       
   // if the user has not set the extent, but has set the VOI
-  // set the zaxis extent to the VOI z axis
+  // set the z axis extent to the VOI z axis
   if (this->DataExtent[4]==0 && this->DataExtent[5] == 0 &&
      (this->DataVOI[4] || this->DataVOI[5]))
     {
@@ -291,11 +406,8 @@ void vtkGdcmReader::ExecuteInformation()
   this->DataExtent[1] = this->NumColumns - 1;
   this->DataExtent[2] = 0;
   this->DataExtent[3] = this->NumLines - 1;
-  if(this->FileNameList.size() > 1)
-    {
-    this->DataExtent[4] = 0;
-    this->DataExtent[5] = this->TotalNumberOfPlanes - 1;
-    }
+  this->DataExtent[4] = 0;
+  this->DataExtent[5] = this->TotalNumberOfPlanes - 1;
   
   // We don't need to positionate the Endian related stuff (by using
   // this->SetDataByteOrderToBigEndian() or SetDataByteOrderToLittleEndian()
@@ -334,6 +446,9 @@ void vtkGdcmReader::ExecuteInformation()
     this->SetDataScalarTypeToInt();
     }
 
+  //Set number of scalar components:
+  this->SetNumberOfScalarComponents(this->NumComponents);
+
   vtkImageReader::ExecuteInformation();
 }
 
@@ -346,21 +461,49 @@ size_t vtkGdcmReader::LoadImageInMemory(
              const unsigned long UpdateProgressTarget,
              unsigned long & UpdateProgressCount)
 {
-  vtkDebugMacro("Copying to memmory image" << FileName.c_str());
+  vtkDebugMacro("Copying to memory image" << FileName.c_str());
   gdcmFile GdcmFile(FileName.c_str());
-  size_t size = GdcmFile.GetImageDataSize();
+  size_t size;
 
   // If the data structure of vtk for image/volume representation
   // were straigthforwards the following would suffice:
   //    GdcmFile.GetImageDataIntoVector((void*)Dest, size);
-  // But vtk chose to invert the lines of an image, that is the last
+  // But vtk chooses to invert the lines of an image, that is the last
   // line comes first (for some axis related reasons?). Hence we need
-  // to load the image line by line, starting from the end:
-  int NumColumns = GdcmFile.GetXSize();
-  int NumLines   = GdcmFile.GetYSize();
-  int NumPlanes  = GdcmFile.GetZSize();
-  int LineSize   = NumColumns * GdcmFile.GetPixelSize();
-  unsigned char * Source      = (unsigned char*)GdcmFile.GetImageData();
+  // to load the image line by line, starting from the end.
+  int NumColumns = GdcmFile.GetHeader()->GetXSize();
+  int NumLines   = GdcmFile.GetHeader()->GetYSize();
+  int NumPlanes  = GdcmFile.GetHeader()->GetZSize();
+  int LineSize   = NumComponents * NumColumns * GdcmFile.GetHeader()->GetPixelSize();
+
+  unsigned char * Source;
+  if( GdcmFile.GetHeader()->HasLUT() )
+    {
+    size               = GdcmFile.GetImageDataSizeRaw();
+    Source             = (unsigned char*) GdcmFile.GetImageDataRaw();
+    unsigned char *Lut =                  GdcmFile.GetHeader()->GetLUTRGBA();
+  
+    if(!this->LookupTable) this->LookupTable = vtkLookupTable::New();
+    this->LookupTable->SetNumberOfTableValues(256);
+    for (int tmp=0; tmp<256; tmp++)
+      {
+      this->LookupTable->SetTableValue(tmp,
+        (float)Lut[4*tmp+0]/255.0,
+        (float)Lut[4*tmp+1]/255.0,
+        (float)Lut[4*tmp+2]/255.0,
+        1);
+      }
+    this->LookupTable->SetRange(0,255);
+    vtkDataSetAttributes *a=this->GetOutput()->GetPointData();
+    a->GetScalars()->SetLookupTable(this->LookupTable);
+    free(Lut);
+    }
+  else
+    {
+    size        = GdcmFile.GetImageDataSize();
+    Source      = (unsigned char*)GdcmFile.GetImageData();
+    }
+  unsigned char * pSource     = Source; //pointer for later deletion
   unsigned char * Destination = Dest + size - LineSize;
 
   for (int plane = 0; plane < NumPlanes; plane++)
@@ -379,76 +522,92 @@ size_t vtkGdcmReader::LoadImageInMemory(
       UpdateProgressCount++;
       }
     }
+  //GetImageData allocate a (void*)malloc, remove it:
+  free(pSource);
+
   return size;
 }
 
 //----------------------------------------------------------------------------
-// Update -> UpdateData -> Execute -> ExecuteData (see vtkSource.cxx for
-// last step.
+// Update => ouput->Update => UpdateData => Execute => ExecuteData 
+// (see vtkSource.cxx for last step).
 // This function (redefinition of vtkImageReader::ExecuteData, see 
 // VTK/IO/vtkImageReader.cxx) reads a data from a file. The datas
-// extent/axes are assumed to be the
-// same as the file extent/order.
+// extent/axes are assumed to be the same as the file extent/order.
 void vtkGdcmReader::ExecuteData(vtkDataObject *output)
 {
-  if (this->FileNameList.empty())
+  if (this->InternalFileNameList.empty())
     {
     vtkErrorMacro("A least a valid FileName must be specified.");
     return;
     }
 
+  // FIXME : extraneous parsing of header is made when allocating OuputData
   vtkImageData *data = this->AllocateOutputData(output);
   data->SetExtent(this->DataExtent);
   data->GetPointData()->GetScalars()->SetName("DicomImage-Volume");
 
-  // The memory size for a full stack of images of course depends
-  // on the number of planes and the size of each image:
-  size_t StackNumPixels = this->NumColumns * this->NumLines
-                        * this->TotalNumberOfPlanes;
-  size_t stack_size = StackNumPixels * this->PixelSize;
-  // Allocate pixel data space itself.
-  unsigned char *mem = new unsigned char [stack_size];
-
-  // Variables for the UpdateProgress. We shall use 50 steps to signify
-  // the advance of the process:
-  unsigned long UpdateProgressTarget = (unsigned long) this->NumLines
-                                     * this->TotalNumberOfPlanes
-                                     / 50.0;
-  // The actual advance measure:
-  unsigned long UpdateProgressCount = 0;
-
-  // Feeling the allocated memory space with each image/volume:
-  unsigned char * Dest = mem;
-  for (std::list<std::string>::iterator FileName  = FileNameList.begin();
-                                        FileName != FileNameList.end();
-                                      ++FileName)
-    { 
-    // Images that were tagged as unreadable in CheckFileCoherence()
-    // are substituted with a black image to let the caller visually
-    // notice something wrong is going on:
-    if (*FileName != "GDCM_UNREADABLE")
-      {
-      Dest += this->LoadImageInMemory(*FileName, Dest,
-                                      UpdateProgressTarget,
-                                      UpdateProgressCount);
-      } else {
-      // We insert a black image in the stack for the user to be aware that
-      // this image/volume couldn't be loaded. We simply skip one image
-      // size:
-      Dest += this->NumColumns * this->NumLines * this->PixelSize;
-      // Update progress related:
-      UpdateProgressCount += this->NumLines;
-      if (!(UpdateProgressCount%UpdateProgressTarget))
+  // Test if output has valid extent
+  // Prevent memory errors
+  if((this->DataExtent[1]-this->DataExtent[0]>=0) &&
+     (this->DataExtent[3]-this->DataExtent[2]>=0) &&
+     (this->DataExtent[5]-this->DataExtent[4]>=0))
+    {
+    // The memory size for a full stack of images of course depends
+    // on the number of planes and the size of each image:
+    size_t StackNumPixels = this->NumColumns * this->NumLines
+                          * this->TotalNumberOfPlanes * this->NumComponents;
+    size_t stack_size = StackNumPixels * this->PixelSize;
+    // Allocate pixel data space itself.
+    unsigned char *mem = new unsigned char [stack_size];
+
+    // Variables for the UpdateProgress. We shall use 50 steps to signify
+    // the advance of the process:
+    unsigned long UpdateProgressTarget = (unsigned long) ceil (this->NumLines
+                                       * this->TotalNumberOfPlanes
+                                       / 50.0);
+    // The actual advance measure:
+    unsigned long UpdateProgressCount = 0;
+
+    // Feeling the allocated memory space with each image/volume:
+    unsigned char * Dest = mem;
+    for (std::list<std::string>::iterator FileName  = InternalFileNameList.begin();
+                                          FileName != InternalFileNameList.end();
+                                        ++FileName)
+      { 
+      // Images that were tagged as unreadable in CheckFileCoherence()
+      // are substituted with a black image to let the caller visually
+      // notice something wrong is going on:
+      if (*FileName != "GDCM_UNREADABLE")
         {
-        this->UpdateProgress(UpdateProgressCount/(50.0*UpdateProgressTarget));
-        }
-      } // Else, file not loadable
-    } // Loop on files
-
-  // The "size" of the vtkScalars data is expressed in number of points,
-  // and is not the memory size representing those points:
-  data->GetPointData()->GetScalars()->SetVoidArray(mem, StackNumPixels, 0);
-  this->Modified();
+        // Update progress related for good files is made in LoadImageInMemory
+        Dest += this->LoadImageInMemory(*FileName, Dest,
+                                        UpdateProgressTarget,
+                                        UpdateProgressCount);
+        } else {
+        // We insert a black image in the stack for the user to be aware that
+        // this image/volume couldn't be loaded. We simply skip one image
+        // size:
+        Dest += this->NumColumns * this->NumLines * this->PixelSize;
+
+        // Update progress related for bad files:
+        UpdateProgressCount += this->NumLines;
+        if (UpdateProgressTarget > 0)
+                     {
+          if (!(UpdateProgressCount%UpdateProgressTarget))
+            {
+                       this->UpdateProgress(UpdateProgressCount/(50.0*UpdateProgressTarget));
+            }
+                     }
+        } // Else, file not loadable
+      } // Loop on files
+
+    // The "size" of the vtkScalars data is expressed in number of points,
+    // and is not the memory size representing those points:
+    data->GetPointData()->GetScalars()->SetVoidArray(mem, StackNumPixels, 0);
+    //don't know why it's here, it's calling one more time ExecuteInformation:
+    //this->Modified();
+    }
 }
 
 //----------------------------------------------------------------------------