[Vv] Fwd: Bug in clitkDicom2Image

Benoît Presles benoit.presles at u-bourgogne.fr
Thu Jul 20 20:21:28 CEST 2017


Dear vv users,

Please find attached the fix I did.

Best regards,
Ben


-------- Message transféré --------
Sujet : 	Bug in clitkDicom2Image
Date : 	Wed, 12 Jul 2017 09:39:24 +0200
De : 	Benoît Presles <benoit.presles at u-bourgogne.fr>
Pour : 	vv at creatis.insa-lyon.fr



Dear Thomas,

I think there is a bug in the clitkDicom2Image tool. clitkDicom2Image does not take into account the transformation in the dicom header (Image Orientation (Patient) tag in the dicom header).
Indeed, the same two images, the first one in "dicom format" (vv-->Open Dicom) and the other in "image format" generated thanks to clitkDicom2Image, do not overlap.

You can download my test data here: https://cloud.u-bourgogne.fr/index.php/s/dNkQJuvXepdmx7e. You will find a dicom image with an orientation and the corresponding image generated thanks to clitkDicom2Image.


Thanks for your help,
Best regards,
Ben

-------------- next part --------------
/*=========================================================================
  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv

  Authors belong to:
  - University of LYON              http://www.universite-lyon.fr/
  - Léon Bérard cancer center       http://www.centreleonberard.fr
  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the copyright notices for more information.

  It is distributed under dual licence

  - BSD        See included LICENSE.txt file
  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/

// clitk includes
#include "clitkDicom2Image_ggo.h"
#include "clitkCommon.h"
#include "clitkImageCommon.h"
#include "vvImageReader.h"
#include "vvImageWriter.h"
#include <gdcmFile.h>
#include <vtkVersion.h>
#include <vtkImageChangeInformation.h>
#if GDCM_MAJOR_VERSION == 2
  #include <gdcmImageHelper.h>
  #include <gdcmAttribute.h>
  #include <gdcmReader.h>
#endif

#include <set>

//====================================================================
int main(int argc, char * argv[])
{
  // init command line
  GGO(clitkDicom2Image, args_info);
  std::vector<std::string> input_files;
  ///if std_input is given, read the input files from stdin
  if (args_info.std_input_given) {
    while (true) {
      std::string tmp;
      std::cin >> tmp;
      if(std::cin.good())
        input_files.push_back(tmp);
      else break;
    }
    args_info.inputs_num=input_files.size();
  } else for (unsigned int i=0; i<args_info.inputs_num; i++)
      input_files.push_back(args_info.inputs[i]);


  //===========================================
  /// Get slices locations ...
  int series_number = -1;
  std::set<int> series_numbers;
  std::map< int, std::vector<double> > theorigin;
  std::map< int, std::vector<double> > theorientation;
  std::map< int, std::vector<double> > sliceLocations;
  std::map< int, std::vector<std::string> > seriesFiles;
#if GDCM_MAJOR_VERSION == 2
  if (args_info.verbose_flag)
    std::cout << "Using GDCM-2.x" << std::endl;
#else
  if (args_info.verbose_flag) {
    std::cout << "Not using GDCM-2.x" << std::endl;
    std::cout<< "The image orientation is not supported with this version of GDCM" <<std::endl;
  }
#endif
  for(unsigned int i=0; i<args_info.inputs_num; i++) {
    if (args_info.verbose_flag)
        std::cout << "Reading <" << input_files[i] << std::endl;
#if GDCM_MAJOR_VERSION == 2
    gdcm::Reader hreader;
    hreader.SetFileName(input_files[i].c_str());
    hreader.Read();
    gdcm::DataSet& ds = hreader.GetFile().GetDataSet();

    if (args_info.extract_series_flag) {
      gdcm::Attribute<0x20,0x11> series_number_att;
      series_number_att.SetFromDataSet(ds);
      series_number = series_number_att.GetValue();
    }
    
    series_numbers.insert(series_number);
    theorigin[series_number] = gdcm::ImageHelper::GetOriginValue(hreader.GetFile());
    theorientation[series_number] = gdcm::ImageHelper::GetDirectionCosinesValue(hreader.GetFile());
    double n1 = theorientation[series_number][1]*theorientation[series_number][5]-
                theorientation[series_number][2]*theorientation[series_number][4];
    double n2 = theorientation[series_number][3]*theorientation[series_number][2]-
                theorientation[series_number][5]*theorientation[series_number][0];
    double n3 = theorientation[series_number][0]*theorientation[series_number][4]-
                theorientation[series_number][1]*theorientation[series_number][3];
    double sloc = theorigin[series_number][0]*n1+
                  theorigin[series_number][1]*n2+
                  theorigin[series_number][2]*n3;
    sliceLocations[series_number].push_back(sloc);
    seriesFiles[series_number].push_back(input_files[i]);
    
    gdcm::Attribute<0x28, 0x100> pixel_size;
    pixel_size.SetFromDataSet(ds);
    /* if (pixel_size.GetValue() != 16)
       {
       std::cerr << "Pixel type not 2 bytes ! " << std::endl;
       std::cerr << "In file " << input_files[i] << std::endl;
       exit(0);
       }
    */
#else
  gdcm::File *header = new gdcm::File();
  header->SetFileName(input_files[i]);
  header->SetMaxSizeLoadEntry(16384); // required ?
  header->Load();

  if (args_info.extract_series_flag) {
    series_number = atoi(header->GetEntryValue(0x20,0x11).c_str());
  }
  
  series_numbers.insert(series_number);
  theorigin[series_number].resize(3);
  theorigin[series_number][0] = header->GetXOrigin();
  theorigin[series_number][1] = header->GetYOrigin();
  theorigin[series_number][2] = header->GetZOrigin();
  sliceLocations[series_number].push_back(theorigin[series_number][2]);
  seriesFiles[series_number].push_back(input_files[i]);
  /*if (header->GetPixelSize() != 2) {
    std::cerr << "Pixel type 2 bytes ! " << std::endl;
    std::cerr << "In file " << input_files[i] << std::endl;
    exit(0);
  }
  */
#endif
  }

  //===========================================
  // Sort slices locations ...
  std::set<int>::iterator sn = series_numbers.begin();
  while ( sn != series_numbers.end() ) {
    std::vector<double> locs = sliceLocations[*sn];
    std::vector<double> origin = theorigin[*sn];
    std::vector<std::string> files = seriesFiles[*sn];
    std::vector<int> sliceIndex;
    clitk::GetSortedIndex(locs, sliceIndex);
    if (args_info.verbose_flag) {
      std::cout << locs[sliceIndex[0]] << " -> "
                << sliceIndex[0] << " / " << 0 << " => "
                << "0 mm "
                << files[sliceIndex[0]]
                << std::endl;
      for(unsigned int i=1; i<sliceIndex.size(); i++) {
        std::cout << locs[sliceIndex[i]] << " -> "
                  << sliceIndex[i] << " / " << i << " => "
                  << locs[sliceIndex[i]] - locs[sliceIndex[i-1]]
                  << "mm "
                  << files[sliceIndex[i]]
                  << std::endl;
      }
    }

    //===========================================
    // Analyze slices locations ...
    double currentDist;
    double dist=0;
    double tolerance = args_info.tolerance_arg;
    double previous = locs[sliceIndex[0]];
    for(unsigned int i=1; i<sliceIndex.size(); i++) {
      currentDist = locs[sliceIndex[i]]-previous;
      if (i!=1) {
        if (fabs(dist-currentDist) > tolerance) {
          std::cout << "ERROR : " << std::endl
                    << "Current slice pos is  = " << locs[sliceIndex[i]] << std::endl
                    << "Previous slice pos is = " << previous << std::endl
                    << "Current file is       = " << files[sliceIndex[i]] << std::endl
                    << "Current index is      = " << i << std::endl
                    << "Current sortindex is  = " << sliceIndex[i] << std::endl
                    << "Current slice diff    = " << dist << std::endl
                    << "Current error         = " << fabs(dist-currentDist) << std::endl;
          exit(1);
        }
      } else dist = currentDist;
      previous = locs[sliceIndex[i]];
    }

    //===========================================
    // Create ordered vector of filenames
    std::vector<std::string> sorted_files;
    sorted_files.resize(sliceIndex.size());
    for(unsigned int i=0; i<sliceIndex.size(); i++)
      sorted_files[i] = files[ sliceIndex[i] ];

    //===========================================
    // Read write serie
    vvImageReader::Pointer reader = vvImageReader::New();
    reader->SetInputFilenames(sorted_files);
    reader->Update(vvImageReader::DICOM);
    if (reader->GetLastError().size() != 0) {
      std::cerr << reader->GetLastError() << std::endl;
      return 1;
    }
    
    vvImage::Pointer image = reader->GetOutput();
    vtkImageData* vtk_image = image->GetFirstVTKImageData();
    vtkImageChangeInformation* modifier = vtkImageChangeInformation::New();
    if  (args_info.focal_origin_given) {
      std::vector<double> spacing = image->GetSpacing();
      std::vector<int> size = image->GetSize();
      origin[0] = -spacing[0]*size[0]/2.0;
      origin[1] = -spacing[1]*size[1]/2.0;
#if VTK_MAJOR_VERSION <= 5
      modifier->SetInput(vtk_image);
#else
      modifier->SetInputData(vtk_image);
#endif
      modifier->SetOutputOrigin(origin[0], origin[1], locs[sliceIndex[0]]);
      modifier->Update();
      vvImage::Pointer focal_image = vvImage::New();
      focal_image->AddVtkImage(modifier->GetOutput());
      image = focal_image;
    }

    std::string outfile;
    if (series_numbers.size() == 1)
      outfile = args_info.output_arg;
    else {
      std::ostringstream name;
      name << *sn << "_" << args_info.output_arg;
      outfile = name.str();
    }
    //Check on transform
    bool bId = true;
    for(unsigned int i=0; i<4; i++) {
      for(unsigned int j=0; j<4; j++) {
        double elt = image->GetTransform()[0]->GetMatrix()->GetElement(i,j);
        if(i==j && elt!=1.) {
          bId = false;
        }
        if(i!=j && elt!=0.) {
          bId = false;
        }
      }
    }
    vvImageWriter::Pointer writer = vvImageWriter::New();
    writer->SetInput(image);
    if(!bId) {
        writer->SetSaveTransform(true);
    }
    writer->SetOutputFileName(outfile);
    writer->Update();

    modifier->Delete();
    
    sn++;
  }
  
  return 0;
}
-------------- next part --------------
# file clitkDicom2Image.ggo
package "clitk"
version "Try to convert a DICOM into an image (.hdr, .vox...)"

option "config"	  	 -	"Config file"		 string  	no
option "verbose"	 v "Verbose"		 	 flag off
option "tolerance"	 t "Tolerance for slice position"	 double default="0" no
option "output"      o "Output image filename"		string 	yes
option "std_input"   - "Take the like of input file from stdin, to support huge lists of filenames" flag off
option "focal_origin" - "Output files with FOCAL-like origin, instead of the origin present in the dicom files" flag off
option "extract_series" s "Identify different series in the file list and create the MHDs accordinly" flag off


More information about the vv mailing list