[Insight-users] 4D Bspline approximation
Nicholas Tustison
ntustison at gmail.com
Thu Sep 20 07:37:45 EDT 2012
Hi Mathieu,
Perhaps answering your short questions will help us figure out
what's going on.
> - Is it correct to use the same size as the number of control points in itkBSplineScatteredDataPointSetToImageFilter ?
No, they're distinct. The spacing and size combination define the
sampled parametric domain (I did this primarily since the filter
derives from the PointSetToImageFilter base class). So if I have
a 20x20 control point grid defining a b-spline object over a domain of
100mm x 100mm, there are many ways I can define the domain
which won't matter for the internal calculations of the final "PhiLattice".
For example, a size and spacing combination of
size = 101, spacing = 1
will define the same parametric domain as
size = 1001, spacing = 0.1
The only difference will be the size of the "image", i.e. sampled B-spline
object when you call GetOutput() but the GetPhiLattice() will be the same.
> - Should I take care of having enough border control points around the mesh nodes (padding) ? How many nodes ?
In the B-spline classes that I've written (including a rewrite of the
B-spline transform class), I don't make a distinction between the
internal control points and the border control points. When you
specify the parametric domain for the BSplineScattered...Filter
(by specifying the origin, size, and spacing), the control points get
laid out automatically based on that domain and the spline order(s)
including the border ones.
I hope that helps. Let me know if it doesn't and I can try running things on my
end.
Nick
On Sep 19, 2012, at 1:39 PM, "DE CRAENE, MATHIEU" <mathieu.de_craene at philips.com> wrote:
> 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.
>
>
More information about the Insight-users
mailing list