tools/clitkDicom2Image.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18
  19 // clitk includes
  20 #include "clitkDicom2Image_ggo.h"
  21 #include "clitkCommon.h"
  22 #include "clitkImageCommon.h"
  23 #include "vvImageReader.h"
  24 #include "vvImageWriter.h"
  25 #include <gdcmFile.h>
  26 #include <vtkImageChangeInformation.h>
  27 #if GDCM_MAJOR_VERSION == 2
  28   #include <gdcmImageHelper.h>
  29   #include <gdcmAttribute.h>
  30   #include <gdcmReader.h>
  31 #endif
  32
  33 //====================================================================
  34 int main(int argc, char * argv[])
  35 {
  36   // init command line
  37   GGO(clitkDicom2Image, args_info);
  38   std::vector<std::string> input_files;
  39   ///if std_input is given, read the input files from stdin
  40   if (args_info.std_input_given) {
  41     while (true) {
  42       std::string tmp;
  43       std::cin >> tmp;
  44       if(std::cin.good())
  45         input_files.push_back(tmp);
  46       else break;
  47     }
  48     args_info.inputs_num=input_files.size();
  49   } else for (unsigned int i=0; i<args_info.inputs_num; i++)
  50       input_files.push_back(args_info.inputs[i]);
  51
  52
  53   //===========================================
  54   /// Get slices locations ...
  55   std::vector<double> theorigin(3);
  56   std::vector<double> sliceLocations;
  57   for(unsigned int i=0; i<args_info.inputs_num; i++) {
  58     //std::cout << "Reading <" << input_files[i] << std::endl;
  59 #if GDCM_MAJOR_VERSION == 2
  60     if (args_info.verbose_flag)
  61       std::cout << "Using GDCM-2.x" << std::endl;
  62     gdcm::Reader hreader;
  63     hreader.SetFileName(input_files[i].c_str());
  64     hreader.Read();
  65     theorigin = gdcm::ImageHelper::GetOriginValue(hreader.GetFile());
  66     sliceLocations.push_back(theorigin[2]);
  67     gdcm::Attribute<0x28, 0x100> pixel_size;
  68     gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
  69     pixel_size.SetFromDataSet(ds);
  70     if (pixel_size.GetValue() != 16)
  71     {
  72       std::cerr << "Pixel type 2 bytes ! " << std::endl;
  73       std::cerr << "In file " << input_files[i] << std::endl;
  74       exit(0);
  75     }
  76 #else
  77     if (args_info.verbose_flag)
  78       std::cout << "Not using GDCM-2.x" << std::endl;
  79   gdcm::File *header = new gdcm::File();
  80   header->SetFileName(input_files[i]);
  81   header->SetMaxSizeLoadEntry(16384); // required ?
  82   header->Load();
  83   theorigin[0] = header->GetXOrigin();
  84   theorigin[1] = header->GetYOrigin();
  85   theorigin[2] = header->GetZOrigin();
  86   sliceLocations.push_back(theorigin[2]);
  87   if (header->GetPixelSize() != 2) {
  88     std::cerr << "Pixel type 2 bytes ! " << std::endl;
  89     std::cerr << "In file " << input_files[i] << std::endl;
  90     exit(0);
  91   }
  92 #endif
  93   }
  94
  95   //===========================================
  96   // Sort slices locations ...
  97   std::vector<int> sliceIndex;
  98   clitk::GetSortedIndex(sliceLocations, sliceIndex);
  99   if (args_info.verboseSliceLocation_flag) {
 100     std::cout << sliceLocations[sliceIndex[0]] << " -> "
 101               << sliceIndex[0] << " / " << 0 << " => "
 102               << "0 mm "
 103               << input_files[sliceIndex[0]]
 104               << std::endl;
 105     for(unsigned int i=1; i<sliceIndex.size(); i++) {
 106       std::cout << sliceLocations[sliceIndex[i]] << " -> "
 107                 << sliceIndex[i] << " / " << i << " => "
 108                 << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]]
 109                 << "mm "
 110                 << input_files[sliceIndex[i]]
 111                 << std::endl;
 112     }
 113   }
 114
 115   //===========================================
 116   // Analyze slices locations ...
 117   double currentDist;
 118   double dist=0;
 119   double tolerance = args_info.tolerance_arg;
 120   double previous = sliceLocations[sliceIndex[0]];
 121   for(unsigned int i=1; i<sliceIndex.size(); i++) {
 122     currentDist = sliceLocations[sliceIndex[i]]-previous;
 123     if (i!=1) {
 124       if (fabs(dist-currentDist) > tolerance) {
 125         std::cout << "ERROR : " << std::endl
 126                   << "Current slice pos is  = " << sliceLocations[sliceIndex[i]] << std::endl
 127                   << "Previous slice pos is = " << previous << std::endl
 128                   << "Current file is       = " << input_files[sliceIndex[i]] << std::endl
 129                   << "Current index is      = " << i << std::endl
 130                   << "Current sortindex is  = " << sliceIndex[i] << std::endl
 131                   << "Current slice diff    = " << dist << std::endl
 132                   << "Current error         = " << fabs(dist-currentDist) << std::endl;
 133         exit(1);
 134       }
 135     } else dist = currentDist;
 136     previous = sliceLocations[sliceIndex[i]];
 137   }
 138
 139   //===========================================
 140   // Create ordered vector of filenames
 141   std::vector<std::string> sorted_files;
 142   sorted_files.resize(sliceIndex.size());
 143   for(unsigned int i=0; i<sliceIndex.size(); i++)
 144     sorted_files[i] = input_files[ sliceIndex[i] ];
 145
 146   //===========================================
 147   // Read write serie
 148   vvImageReader::Pointer reader = vvImageReader::New();
 149   reader->SetInputFilenames(sorted_files);
 150   reader->Update(vvImageReader::DICOM);
 151   if (reader->GetLastError().size() != 0) {
 152     std::cerr << reader->GetLastError() << std::endl;
 153     return 1;
 154   }
 155
 156   vvImage::Pointer image = reader->GetOutput();
 157   vtkImageData* vtk_image = image->GetFirstVTKImageData();
 158   vtkImageChangeInformation* modifier = vtkImageChangeInformation::New();
 159   if  (args_info.focal_origin_given) {
 160     std::vector<double> spacing = image->GetSpacing();
 161     std::vector<int> size = image->GetSize();
 162     theorigin[0] = -spacing[0]*size[0]/2.0;
 163     theorigin[1] = -spacing[1]*size[1]/2.0;
 164     modifier->SetInput(vtk_image);
 165     modifier->SetOutputOrigin(theorigin[0], theorigin[1], sliceLocations[sliceIndex[0]]);
 166     modifier->Update();
 167     vvImage::Pointer focal_image = vvImage::New();
 168     focal_image->AddVtkImage(modifier->GetOutput());
 169     image = focal_image;
 170   }
 171
 172   vvImageWriter::Pointer writer = vvImageWriter::New();
 173   writer->SetInput(image);
 174   writer->SetOutputFileName(args_info.output_arg);
 175   writer->Update();
 176
 177   modifier->Delete();
 178
 179   return 0;
 180 }