[vtkusers] Patch to sample voxels corresponding to input points from gaussian_splatter

[lynx.abraxas at freenet.de]

Hello!


I tried to add a function to vtkGaussianSplatter that probes the output at the
given point. A patch is attached but so far the returned  result  in  my  test
code is not correct. What am I doing wrong?

I also noticed:
-In  the  function  RequestData  i  is  of vtkIdType but j, k are int. Is that
correct?
-ComputeModelBounds does not work with just one point (SetModelBounds is
necessary then)

Thanks for any help or hints
Lynx


On 03/02/10 22:44:22, lynx.abraxas at freenet.de wrote:
> Hello!
> 
> 
> Having  used  GaussianSplatter  to splat some points, I'd now need to read the
> intensity values from the voxels  that  are  closest  to  the  original  splat
> centres  (basicly  the greyvalue at the position of the input point). Is there
> an elegant way to do that with vtk?
> Is there also a way in vtk to get the overall voxel sum of the splats image or
> would I have to use itkImageStatistics for that?

-------------- next part --------------
//doing a weighted spalt
//calculating image/spalt stats


#include <iostream>
//#include <fstream>
//#include <vector>
#include <string>
#include <sstream> //std::stringstream
#include <cmath>
 
#include <vtkGaussianSplatter.h>
#include <vtkFloatArray.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkImageCast.h>
//#include <vtkTIFFWriter.h>
#include <vtkMetaImageWriter.h> 

#include <itkVTKImageToImageFilter.h>
#include <itkStatisticsImageFilter.h>
 
int main(int argc, char **argv) {
 

    vtkPoints* Points = vtkPoints::New();
    vtkFloatArray* weights = vtkFloatArray::New();
    weights->SetNumberOfComponents(1);
    weights->SetName ("weight");
    Points->InsertNextPoint(0, 0, 0);
    weights->InsertNextValue(atof(argv[1]));

    vtkPolyData* polydata = vtkPolyData::New();
    polydata->SetPoints(Points);
    polydata->GetPointData()->SetScalars(weights);

    vtkGaussianSplatter *Splatter = vtkGaussianSplatter::New();
    Splatter->SetInput(polydata);
    Splatter->SetSampleDimensions(atoi(argv[5]),atoi(argv[6]),atoi(argv[7])); //set the resolution of the final! volume
    double r= atof(argv[2]);
    //doesn't work with just one point //let the bonds be computed since r is relative not absolut!
    Splatter->SetModelBounds(-r,r, -r,r, -r,r); //above this no change of the sum!
    //Splatter->SetModelBounds(-2*r,2*r, -2*r,2*r, -2*r,2*r);
    Splatter->SetRadius(r); //GaussianSplat truncated outside Radius
    Splatter->SetExponentFactor(atof(argv[3])); //GaussianSplat decay value
    Splatter->SetAccumulationModeToSum();
    Splatter->ScalarWarpingOn();
    Splatter->CappingOff(); //don't pad with 0
    //Splatter->SetScaleFactor(atof(argv[4]));//scale comes from each point
 
    vtkImageCast* cast = vtkImageCast::New();
    cast->SetInputConnection(Splatter->GetOutputPort());
    //cast->SetOutputScalarTypeToUnsignedChar();
    cast->SetOutputScalarTypeToFloat();


    //vtkTIFFWriter *writer = vtkTIFFWriter::New();
    vtkMetaImageWriter *writer = vtkMetaImageWriter::New();
    std::cout << "Writinig image... "; 
    writer->SetFileName(argv[4]);
    writer->SetFileDimensionality(3);
    //writer->SetInput(cast->GetOutputPort());
    writer->SetInputConnection(cast->GetOutputPort());
    //writer->SetCompressionToNoCompression();
    writer->SetCompression(0);
    writer->Write();
    std::cout << "done." << std::endl; 


    const int dim = 3;
    typedef float PixelType;
    typedef itk::Image< PixelType, dim >    ImageType;

    typedef itk::VTKImageToImageFilter<ImageType> ConnectorType;

    ConnectorType::Pointer vtkitkf = ConnectorType::New();
    vtkitkf->SetInput(cast->GetOutput()); //NOT GetOutputPort()!!!

    typedef itk::StatisticsImageFilter<ImageType> SFType;
    SFType::Pointer stats= SFType::New();

    stats->SetInput(vtkitkf->GetOutput());
    stats->Update();

    cout << "mean: " << stats->GetMean()
         << " std: " << stats->GetSigma()
         << " var: " << stats->GetVariance()
         << " sum: " << stats->GetSum()
         << " max: " << stats->GetMaximum()
         << " min: " << stats->GetMinimum()
         <<  endl;
                                
    double pp[3], v;

    pp[0]= atof(argv[8]);
    pp[1]= atof(argv[9]);
    pp[2]= atof(argv[10]);

    int res= Splatter->ProbePoint(Splatter->GetOutput(), pp, &v);
    
    if (res)
        cout << "Prope Point: " << pp[0] << ";" << pp[1] << ";" << pp[2] << " not insied sample data" <<  endl;
    else
        cout << "Prope Point: " << pp[0] << ";" << pp[1] << ";" << pp[2] << "; value: " << v << endl;

    return 0;
 
    }
-------------- next part --------------
diff -aur orig//vtkGaussianSplatter.cxx probe_point//vtkGaussianSplatter.cxx
--- orig//vtkGaussianSplatter.cxx	2010-02-08 13:16:35.938535000 +0100
+++ probe_point//vtkGaussianSplatter.cxx	2010-02-08 15:38:53.164765000 +0100
@@ -596,3 +596,31 @@
   info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
   return 1;
 }
+
+
+
+int vtkGaussianSplatter::ProbePoint(vtkImageData *output, const double *pp, double *value){
+    int i;
+    double x[3];
+
+
+    // Determine the voxel that the point is in
+    for (i=0; i<3; i++)  
+        {
+        x[i] = round((pp[i] - this->Origin[i]) / this->Spacing[i]); //round(), int(), ceil() or floor()???
+        //check if point is inside output data
+        cout << x[i] << endl;
+
+        if ( x[i] < 0 || x[i] >= this->SampleDimensions[i] )
+            {
+            return(-1);
+            }
+        }
+
+    value= output->GetPointData()->GetScalars()->GetTuple(x[0] + x[1]*this->SampleDimensions[0] + x[2]*this->SampleDimensions[0]*this->SampleDimensions[1]);
+
+    cout << x[0] + x[1]*this->SampleDimensions[0] + x[2]*this->SampleDimensions[0]*this->SampleDimensions[1]<< endl;
+
+    return(0);
+    
+    }
diff -aur orig//vtkGaussianSplatter.h probe_point//vtkGaussianSplatter.h
--- orig//vtkGaussianSplatter.h	2010-02-08 13:16:24.998764000 +0100
+++ probe_point//vtkGaussianSplatter.h	2010-02-08 14:46:25.224003000 +0100
@@ -180,6 +180,10 @@
   vtkGetMacro(NullValue,double);
 
   // Description:
+  // Return the splat value closest to the probe point
+  int ProbePoint(vtkImageData *output, const double *pp, double *value);
+
+  // Description:
   // Compute the size of the sample bounding box automatically from the
   // input data. This is an internal helper function.
   void ComputeModelBounds(vtkDataSet *input, vtkImageData *output,