[Insight-users] 4D Bspline approximation

[DE CRAENE, MATHIEU]

Hi Nick et al.




I am trying to use the itkBSplineScatteredDataPointSetToImageFilter to regularize the displacements of volumetric mesh nodes using a 4D Bspline field. The meshes represent the right and left ventricles of the heart over one cycle.

The idea is to enter velocity as "frame to frame displacements" for a subset of all mesh nodes. Then to find the best b spline fit using the filter.

To check that the result is a correct approximation, I integrated the velocity field using TimeVaryingBSplineVelocityFieldTransform and apply the resulting transform to the first mesh.

My problem is that I obtain a zero displacement i.e. an identity transform at every frame.
It could come from an incorrect way of using the "projection" filter or the final transform.

My input data is here: https://dl.dropbox.com/u/57469885/meshes-stacom.zip
And the source code is below.

It would be very helpful if you can spot any obvious misuse of the filters!

Some short questions :

- Is it correct to use the same size as the number of control points in itkBSplineScatteredDataPointSetToImageFilter ?
- Should I take care of having enough border control points around the mesh nodes (padding) ? How many nodes ?


Thanks for any hint !!



Math


==========================================================================

// Image types and other standard things
#include <iostream>
#include <cstdlib>
#include "vtkUnstructuredGrid.h"
#include <vtkUnstructuredGridReader.h>
#include <vtkUnstructuredGridWriter.h>
#include <vtkSmartPointer.h>
#include <vtkCellData.h>
#include "itkTimeVaryingBSplineVelocityFieldTransform.h"
#include "itkBSplineScatteredDataPointSetToImageFilter.h"
#include <itkPointSet.h>
#include <itkImageRegionIterator.h>

int main( int argc, char * argv[] )
{
        if (argc < 3)
        {
                std::cerr << "Usage : " << argv[0] << "InputMeshList.txt outputMeshPrefix" << std::endl;
                exit(EXIT_FAILURE);
        }
        const unsigned int MAX_PATH_LENGTH = 1024;
        const unsigned int ImageDimension = 3;

        int arg = 1;
        const char* InputMeshListFileName = argv[arg]; arg++;
        const char* outputMeshPrefixFileName = argv[arg]; arg++;

        // reading input list
        std::ifstream f (InputMeshListFileName, std::ifstream::in);
        if (f==0)
        {
                return EXIT_FAILURE;
        }

        std::string buffer;
        std::vector<std::string> inputUGFileNames;
        f>>buffer;

        while (!f.eof())
        {
                inputUGFileNames.push_back(buffer);
                f>>buffer;
        }

        // reading all meshes
        std::vector<vtkSmartPointer<vtkUnstructuredGrid>> ugs;
        for (unsigned int t=0; t<inputUGFileNames.size(); t++)
        {
                vtkSmartPointer<vtkUnstructuredGridReader> reader
                        = vtkSmartPointer<vtkUnstructuredGridReader>::New();
                reader->SetFileName(inputUGFileNames[t].c_str());
                reader->Update();
                ugs.push_back(reader->GetOutput());
        }

        // check that we have at least to meshes
        if (ugs.size() <=1)
        {
                std::cerr << "Insufficient (" << ugs.size() << ") number of meshes" << std::endl;
                exit(EXIT_FAILURE);
        }

        // Instantiate velocity and adjust bspline field to mesh bounding box
        typedef itk::TimeVaryingBSplineVelocityFieldTransform< double, ImageDimension > TransformType;
        typedef TransformType::DisplacementVectorType DisplacementVectorType;
        typedef TransformType::TimeVaryingVelocityFieldControlPointLatticeType VelocityFieldType;
        typedef itk::PointSet<DisplacementVectorType, ImageDimension + 1> PointSetType;
        typedef itk::BSplineScatteredDataPointSetToImageFilter<PointSetType, VelocityFieldType> BSplineFilterType;

        // Filling point set with velocity (incr. displ.) data
        PointSetType::Pointer velocityFieldPoints = PointSetType::New();
        velocityFieldPoints->Initialize();

        typedef BSplineFilterType::WeightsContainerType WeightsContainerType;
        WeightsContainerType::Pointer velocityFieldWeights = WeightsContainerType::New();

        PointSetType::PointType spatioTemporalPoint;
        DisplacementVectorType velocity;
        PointSetType::PointIdentifier id = 0;
        for (unsigned int t=0; t<ugs.size()-1; t++)
        {
                for (int p=0; p<ugs[t]->GetNumberOfPoints(); p+=50,id++)
                {
                        for (int d=0; d<ImageDimension; d++)
                        {
                                spatioTemporalPoint[d] = ugs[t]->GetPoints()->GetPoint(p)[d];
                                velocity[d] =
                                (   ugs[t+1]->GetPoints()->GetPoint(p)[d]
                                -       ugs[t  ]->GetPoints()->GetPoint(p)[d] )
                                * (double) ugs.size();
                        }
                        spatioTemporalPoint[ImageDimension] = (double)t / (double) ugs.size();

                        velocityFieldPoints->SetPoint(id,  spatioTemporalPoint);
                        velocityFieldPoints->SetPointData(id, velocity );
                        velocityFieldWeights->InsertElement( id, 1. );
                }
        }

        // Setting up Bspline regression of the velocity field
        BSplineFilterType::Pointer bspliner = BSplineFilterType::New();
        unsigned int nocp[4] = {20,20,20,10};
        BSplineFilterType::SizeType domainSize;
        for (int d=0; d<4; d++)
                domainSize[d] = nocp[d];
        BSplineFilterType::ArrayType numberOfControlPoints(nocp);
        BSplineFilterType::PointType origin;
        BSplineFilterType::SpacingType spacing;
        double bb[6];
        ugs[0]->GetBounds(bb);

        // safety margins
        double offset[4]={5., 5., 5., 0.5};

        spacing[0] = ( bb[1] - bb[0] + 2.*offset[0] ) / (double)(domainSize[0]-3);
        spacing[1] = ( bb[3] - bb[2] + 2.*offset[1] ) / (double)(domainSize[1]-3);
        spacing[2] = ( bb[5] - bb[4] + 2.*offset[2] ) / (double)(domainSize[2]-3);
        spacing[3] = ( 1.                        + 2.*offset[3] ) / (double)(domainSize[3]-3);

        origin[0] = bb[0] - offset[0] - spacing[0];
        origin[1] = bb[2] - offset[1] - spacing[1];
        origin[2] = bb[4] - offset[2] - spacing[2];
        origin[3] =       - offset[3] - spacing[3];

        bspliner->SetOrigin( origin );
        bspliner->SetSpacing( spacing );
        bspliner->SetNumberOfControlPoints( numberOfControlPoints );
        bspliner->SetNumberOfLevels( 1 );
        bspliner->SetSplineOrder( 3 );
        bspliner->SetSize( domainSize );
        bspliner->SetInput( velocityFieldPoints );
        bspliner->SetPointWeights( velocityFieldWeights );
        bspliner->GenerateOutputImageOff();

        try
        {
                bspliner->Update();
        }
        catch( itk::ExceptionObject & err )
        {
                std::cerr << "ExceptionObject caught !" << std::endl;
                std::cerr << err << std::endl;
                return EXIT_FAILURE;
        }

        TransformType::TimeVaryingVelocityFieldControlPointLatticePointer lattice = bspliner->GetPhiLattice();

        // Transformation that is supposed to approximate our velocity field
        TransformType::Pointer transform = TransformType::New();
        transform->SetTimeVaryingVelocityFieldControlPointLattice(lattice);
        transform->SetSplineOrder(bspliner->GetSplineOrder()[0]);

        // Allocating
        std::vector<vtkSmartPointer<vtkUnstructuredGrid>> out_ugs;
        for (unsigned int t=0; t<ugs.size(); t++)
        {
                vtkSmartPointer<vtkUnstructuredGrid> oug = vtkSmartPointer<vtkUnstructuredGrid>::New();
                oug->DeepCopy(ugs[0]);
                out_ugs.push_back(oug);
        }

        // Recomputing transform for all points
        TransformType::OutputPointType outp;
        for (unsigned int t=0; t<ugs.size(); t++)
        {
                transform->SetUpperTimeBound((double)t / (double)ugs.size());
                transform->SetLowerTimeBound((double)0);
                transform->IntegrateVelocityField();

                for (int p=0; p<ugs[0]->GetNumberOfPoints(); p++)
                {
                        for (int d=0; d<ImageDimension; d++)
                        {
                                spatioTemporalPoint[d] = ugs[0]->GetPoints()->GetPoint(p)[d];
                        }
                        spatioTemporalPoint[ImageDimension] = 0.;

                        outp = transform->TransformPoint(ugs[0]->GetPoints()->GetPoint(p));
                        out_ugs[t]->GetPoints()->SetPoint(p, outp.GetDataPointer());
                }
        }

        // Writing output
        char outputFileName[MAX_PATH_LENGTH];
        for (unsigned int t=0; t<ugs.size(); t++)
        {
                sprintf(outputFileName, "%s%03d.vtk", outputMeshPrefixFileName, t);
                vtkSmartPointer<vtkUnstructuredGridWriter> wr =                                                                 vtkSmartPointer<vtkUnstructuredGridWriter>::New();
                wr->SetFileName(outputFileName);
                wr->SetInput(out_ugs[t]);
                wr->Update();
        }


        exit(EXIT_SUCCESS);
}



________________________________
The information contained in this message may be confidential and legally protected under applicable law. The message is intended solely for the addressee(s). If you are not the intended recipient, you are hereby notified that any use, forwarding, dissemination, or reproduction of this message is strictly prohibited and may be unlawful. If you are not the intended recipient, please contact the sender by return e-mail and destroy all copies of the original message.