Add 2 options to clitkImage2Dicom

[clitk.git] / tools / clitkGammaIndex.cxx

diff --git a/tools/clitkGammaIndex.cxx b/tools/clitkGammaIndex.cxx
index 65031a6c22141cc25086596b9ee1ea966ff5ec76..1f5fd03f3a78c905ecd666dc2376569b674e9fa3 100644 (file)
--- a/tools/clitkGammaIndex.cxx
+++ b/tools/clitkGammaIndex.cxx
@@ -15,6 +15,15 @@
   - BSD        See included LICENSE.txt file
   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
 ===========================================================================*/
+
+#include "clitkGammaIndex_ggo.h"
+#include "clitkIO.h"
+#include "clitkDD.h"
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
 #include <vtkPoints.h>
 #include <vtkCellArray.h>
 #include <vtkPointData.h>
@@ -24,17 +33,18 @@
 #include <vtkPNGReader.h>
 #include <vtkPolyData.h>
 #include <vtkCellLocator.h>
-#include <iostream>
-#include <cmath>
-#include <cassert>
-using std::endl;
-using std::cout;
-
-#include "clitkGammaIndex_ggo.h"
 
-#include <clitkIO.h>
 #include <vvImage.h>
 #include <vvImageReader.h>
+#include <vvImageWriter.h>
+
+#include <itkImage.h>
+#include <itkImageRegionIterator.h>
+typedef itk::Image<double> OutputImageType;
+typedef itk::ImageRegionIterator<OutputImageType> OutputImageIterator;
+
+using std::endl;
+using std::cout;
 
 vtkImageData* loadImage(const std::string& filename)
 {
@@ -49,13 +59,14 @@ vtkImageData* loadImage(const std::string& filename)
     return image;
 }
 
-void saveImage(vtkImageData *image,const std::string &filename) {
-    cout << "saving " << filename << endl;
-    vtkImageWriter *writer = vtkMetaImageWriter::New();
-    writer->SetFileName(filename.c_str());
-    writer->SetInput(image);
-    writer->Write();
-    writer->Delete();
+void saveImage(OutputImageType* image,const std::string &filename) {
+    vvImage::Pointer vvimage = vvImage::New();
+    vvimage->AddItkImage(image);
+
+    vvImageWriter::Pointer writer = vvImageWriter::New();
+    writer->SetOutputFileName(filename.c_str());
+    writer->SetInput(vvimage);
+    writer->Update();
 }
 
 
@@ -218,6 +229,17 @@ int main(int argc,char * argv[])
     vtkImageData* reference = loadImage(reference_filename);
     assert(reference);
 
+    // translate target with arguments values
+    // reference is translated instead of target so that the output space stay the same as target
+    {
+       double reference_origin[3];
+        reference->GetOrigin(reference_origin);
+       reference_origin[0] -= args_info.translation_x_arg;
+       reference_origin[1] -= args_info.translation_y_arg;
+       reference_origin[2] -= args_info.translation_z_arg;
+       reference->SetOrigin(reference_origin);
+    }
+
     // intensity normalisation
     if (!use_dose_margin) {
        dose_margin = getMaximum(reference)*dose_rel_margin;
@@ -230,27 +252,40 @@ int main(int argc,char * argv[])
 
     vtkAbstractCellLocator *locator = vtkCellLocator::New();
     locator->SetDataSet(data);
-    data->Delete();
+    DD("here");
+    //    data->Delete();
     locator->CacheCellBoundsOn();
     locator->AutomaticOn();
+    DD("BuildLocator");
     locator->BuildLocator();
+    DD("end BuildLocator");
 
     // load target
-    vtkImageData *target = loadImage(target_filename);
+    vtkImageData* target = loadImage(target_filename);
     assert(target);
 
+
     // allocate output
-    vtkImageData *output = vtkImageData::New();
-    output->SetExtent(target->GetWholeExtent());
-    output->SetOrigin(target->GetOrigin());
-    output->SetSpacing(target->GetSpacing());
-    output->SetScalarTypeToFloat();
-    output->AllocateScalars();
+    OutputImageType::Pointer output = OutputImageType::New();
+    {
+       OutputImageType::SizeType::SizeValueType output_array_size[2];
+       output_array_size[0] = target->GetDimensions()[0];
+       output_array_size[1] = target->GetDimensions()[1];
+       OutputImageType::SizeType output_size;
+       output_size.SetSize(output_array_size);
+       output->SetRegions(OutputImageType::RegionType(output_size));
+       output->SetOrigin(target->GetOrigin());
+       output->SetSpacing(target->GetSpacing());
+       output->Allocate();
+    }
 
     // fill output
-    unsigned long total = 0;
+    unsigned long  kk = 0;
     unsigned long over_one = 0;
-    for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
+    OutputImageIterator iter(output,output->GetLargestPossibleRegion());
+    iter.GoToBegin();
+    DD("loop");
+    while (!iter.IsAtEnd()) {
        double *point = target->GetPoint(kk);
        double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
        assert(value>=0);
@@ -267,27 +302,24 @@ int main(int argc,char * argv[])
        double squared_distance = 0;
 
        locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
-
        double distance = sqrt(squared_distance);
-       output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
-
-       if (value>1) over_one++;
-       total++;
+       iter.Set(distance);
 
+       if (distance>1) over_one++;
+       kk++;
+       ++iter;
     }
 
     if (verbose) {
-       cout << "total=" << total << endl;
+       cout << "total=" << kk << endl;
        cout << "over_one=" << over_one << endl;
-       cout << "ratio=" << static_cast<float>(over_one)/total << endl;
+       cout << "ratio=" << static_cast<float>(over_one)/kk << endl;
     }
 
     locator->Delete();
     target->Delete();
 
     saveImage(output,gamma_filename);
-    output->Delete();
 
     return 0;
 }
-