tools/clitkDicom2Image.cxx

   1 /*=========================================================================
   2
   3 Program:   clitk
   4 Module:    $RCSfile: clitkDicom2Image.cxx,v $
   5 Language:  C++
   6 Date:      $Date: 2010/11/05 00:05:00 $
   7 Version:   $Revision: 1.1 $
   8
   9 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10 l'Image). All rights reserved. See Doc/License.txt or
  11 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13 This software is distributed WITHOUT ANY WARRANTY; without even
  14 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15 PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 /**
  20    =================================================
  21    * @file   clitkDicom2Image.cxx
  22    * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
  23    * @date   26 Oct 2006
  24    *
  25    * @brief
  26    *
  27    =================================================*/
  28
  29 // clitk includes
  30 #include "clitkDicom2Image_ggo.h"
  31 #include <clitkCommon.h>
  32 #include "clitkImageCommon.h"
  33
  34 // itk include
  35 #include "itkImageRegionIterator.h"
  36 #include "itkImageRegionConstIterator.h"
  37
  38 // Why this is needed ???
  39 //const double itk::NumericTraits<double>::Zero = 0.0;
  40
  41 //====================================================================
  42 int main(int argc, char * argv[]) {
  43
  44   // init command line
  45   GGO(clitkDicom2Image, args_info);
  46   std::vector<std::string> input_files;
  47   ///if std_input is given, read the input files from stdin
  48   if (args_info.std_input_given)
  49   {
  50       while (true)
  51       {
  52           std::string tmp;
  53           std::cin >> tmp;
  54           if(std::cin.good())
  55               input_files.push_back(tmp);
  56           else break;
  57       }
  58       args_info.inputs_num=input_files.size();
  59   }
  60   else for (unsigned int i=0;i<args_info.inputs_num;i++)
  61       input_files.push_back(args_info.inputs[i]);
  62
  63
  64   //===========================================
  65   /// Get slices locations ...
  66     std::vector<double> sliceLocations;
  67     for(unsigned int i=0; i<args_info.inputs_num; i++) {
  68       //std::cout << "Reading <" << input_files[i] << std::endl;
  69       gdcm::File * header = clitk::readDicomHeader(input_files[i]);
  70       sliceLocations.push_back(header->GetZOrigin());
  71       // check
  72       /*if (!header->IsSignedPixelData()) {
  73         std::cerr << "Pixel type not signed ! " << std::endl;
  74         std::cerr << "In file " << input_files[i] << std::endl;
  75         exit(0);
  76         }*/
  77       if (header->GetPixelSize() != 2) {
  78         std::cerr << "Pixel type 2 bytes ! " << std::endl;
  79         std::cerr << "In file " << input_files[i] << std::endl;
  80         exit(0);
  81       }
  82     }
  83
  84     //===========================================
  85     // Sort slices locations ...
  86     std::vector<int> sliceIndex;
  87     clitk::GetSortedIndex(sliceLocations, sliceIndex);
  88     if (args_info.verboseSliceLocation_flag) {
  89       std::cout << sliceLocations[sliceIndex[0]] << " -> "
  90                 << sliceIndex[0] << " / " << 0 << " => "
  91                 << "0 mm "
  92                 << input_files[sliceIndex[0]]
  93                 << std::endl;
  94       for(unsigned int i=1; i<sliceIndex.size(); i++) {
  95         std::cout << sliceLocations[sliceIndex[i]] << " -> "
  96                   << sliceIndex[i] << " / " << i << " => "
  97                   << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]]
  98                   << "mm "
  99                   << input_files[sliceIndex[i]]
 100                   << std::endl;
 101       }
 102     }
 103
 104     //===========================================
 105     // Analyze slices locations ...
 106     double currentDist;
 107     double dist=0;
 108     double tolerance = args_info.tolerance_arg;
 109     double previous = sliceLocations[sliceIndex[0]];
 110     for(unsigned int i=1; i<sliceIndex.size(); i++) {
 111       currentDist = sliceLocations[sliceIndex[i]]-previous;
 112       if (i!=1) {
 113         if (fabs(dist-currentDist) > tolerance) {
 114           std::cout << "ERROR : " << std::endl
 115                     << "Current slice pos is  = " << sliceLocations[sliceIndex[i]] << std::endl
 116                     << "Previous slice pos is = " << previous << std::endl
 117                     << "Current file is       = " << input_files[sliceIndex[i]] << std::endl
 118                     << "Current index is      = " << i << std::endl
 119                     << "Current sortindex is  = " << sliceIndex[i] << std::endl
 120                     << "Current slice diff    = " << dist << std::endl
 121                     << "Current error         = " << fabs(dist-currentDist) << std::endl;
 122           exit(1);
 123         }
 124       }
 125       else dist = currentDist;
 126       previous = sliceLocations[sliceIndex[i]];
 127     }
 128
 129     //===========================================
 130     // Create image
 131     gdcm::File * header = clitk::readDicomHeader(input_files[sliceIndex.front()]);
 132     typedef itk::Image<float, 3> ImageType;
 133     typedef itk::Image<signed short, 3> ImageTypeOutput;
 134     ImageType::SpacingType spacing;
 135     ImageType::SizeType    size;
 136     ImageType::PointType    origin;
 137     spacing[0] = header->GetXSpacing();
 138     spacing[1] = header->GetYSpacing();
 139     spacing[2] = dist;
 140     size[0]    = header->GetXSize();
 141     size[1]    = header->GetYSize();
 142     size[2]    = sliceIndex.size();
 143     ///Optional use of special origin scheme used at CLB
 144     if (args_info.focal_origin_flag)
 145     {
 146         origin[0]  = -spacing[0]*size[0]/2;
 147         origin[1]  = -spacing[1]*size[1]/2;
 148     }
 149     else
 150     {
 151         origin[0]  = header->GetXOrigin();
 152         origin[1]  = header->GetYOrigin();
 153     }
 154     origin[2]  = header->GetZOrigin();
 155     ImageTypeOutput::Pointer output = ImageTypeOutput::New();
 156     itk::ImageRegion<3> region;
 157     region.SetSize(size);
 158     output->SetRegions(region);
 159     output->SetOrigin(origin);
 160     output->SetSpacing(spacing);
 161     output->Allocate();
 162
 163     //===========================================
 164     // Fill image
 165     ImageType::Pointer slice;
 166     typedef itk::ImageRegionIterator<ImageTypeOutput> IteratorType;
 167     typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
 168     IteratorType po(output, output->GetLargestPossibleRegion());
 169     po.Begin();
 170     for(unsigned int i=0; i<sliceIndex.size(); i++) {
 171       slice = clitk::readImage<ImageType>(input_files[sliceIndex[i]]);
 172       ConstIteratorType pi(slice, slice->GetLargestPossibleRegion());
 173       while ( !pi.IsAtEnd() ) {
 174         po.Set((signed short)pi.Get());
 175         ++po;
 176         ++pi;
 177       }
 178     }
 179
 180     //===========================================
 181     // Write image
 182     clitk::writeImage<ImageTypeOutput>(output,
 183                                        args_info.output_arg,
 184                                        args_info.verbose_flag);
 185
 186     // this is the end my friend
 187     return 0;
 188 }