[Insight-users] itkBSplineScatteredDataImageFilter exception thrown

[Ian Malone]

Hi,

Thanks for the quick response (and the filter...). Just to be clear, 
when you say 'the main routine is attempting to perform the calculation 
on the boundary', you do mean at the filter -> Update (); statement, 
which is where the exception is thrown?

imalone

Nicholas Tustison wrote:
> HI Ian,
>
> The problem is that the main routine is attempting to perform the calculation on the boundary of the open-ended domain of the B-spline object.  I've revised the code quite a bit since the IJ article and recently Matt McCormick has offered a further fix in ITKv4.  You're right in that changing the spacing such that the evaluation is not performed at 1.0 but rather 0.9999 or something like that should make it work.
>
> Nick
>
>
>
>
> On Dec 14, 2010, at 11:12 AM, Ian Malone wrote:
>
>   
>> Hi,
>>
>> I'm trying to use itkBSplineScatteredDataImageFilter as a curve interpolator and I'm following
>> N-D C^k B-Spline Scattered Data Approximation, Tustison and Gee http://www.insight-journal.org/browse/publication/57
>>
>> Just copying the 3D curve example there, wrapping it in an int main(void) and adding the right #includes (as attached), I made one change which was to output any itk::ExceptionObjects caught. It crashes in ITK 3.16 and 3.20 as follows:
>>
>> 3.16:
>> Test 2: itkBSplineScatteredDataImageFilter exception thrown
>>
>> itk::ExceptionObject (0x807d448)
>> Location: "void itk::BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>::GenerateControlLattice() [with TInputPointSet = itk::PointSet<itk::Vector<double, 3u>, 1u, itk::DefaultStaticMeshTraits<itk::Vector<double, 3u>, 1u, 1u, float, float, itk::Vector<double, 3u> > >, TOutputImage = itk::Image<itk::Vector<double, 3u>, 1u>]"
>> File: /var/drc/software/32bit/itk/include/InsightToolkit/Review/itkBSplineScatteredDataPointSetToImageFilter.txx
>> Line: 633
>> Description: itk::ERROR: PointSetToImageFilter(0x807a8a0): The reparameterized point component 1.002 is outside the corresponding parametric domain of [0, 1].
>>
>>
>> 3.20:
>> Test 2: itkBSplineScatteredDataImageFilter exception thrown
>>
>> itk::ExceptionObject (0x8153f38)
>> Location: "void itk::MultiThreader::SingleMethodExecute()"
>> File: /home/samba/user/malone/testing/itk/InsightToolkit-3.20.0/Code/Common/itkMultiThreader.cxx
>> Line: 471
>> Description: itk::ERROR: MultiThreader(0x8149ac0): Exception occurred during SingleMethodExecute
>>
>> By comparison with Testing/Code/Review/itkBSplineScatteredDataPointSetToImageFilterTest2.cxx I guess that this is because itkBSplineScatteredDataImageFilter is trying to do calculations on the range [0:1] with a spacing 0.001, which requires size to be 1001 (when in fact the size Fill ( 1.0/ spacing [0]) results in 999). However this isn't really mentioned in the DOxygen for it. Also at the input end of things if pointSet -> SetPoint ( i , point ); sets points beyond the range [0,1] then I get a similar error. So, is anyone who knows about itkBSplineScatteredDataPointSetToImageFilter able to tell me whether it's supposed to only operate on [0:1]?
>>
>> Thanks for your time.
>>
>> -- 
>> imalone
>> #include <itkBSplineScatteredDataPointSetToImageFilter.h>
>> #include <itkPointSet.h>
>> #include <itkVectorImage.h>
>> #include <itkImage.h>
>>
>> int main (void) {
>> const unsigned int ParametricDimension = 1;
>> const unsigned int DataDimension = 3;
>> typedef double RealType ;
>> typedef itk :: Vector < RealType , DataDimension > VectorType ;
>> typedef itk :: Image < VectorType , ParametricDimension > ImageType ;
>> typedef itk :: PointSet < VectorType , ParametricDimension > PointSetType ;
>> PointSetType :: Pointer pointSet = PointSetType :: New ();
>> // Sample the helix.
>> for ( RealType t = 0.0; t <= 1.0+1e-10; t += 0.05 )
>>  {
>>  unsigned long i = pointSet -> GetNumberOfPoints ();
>>  PointSetType :: PointType point ;
>>  point [0] = t;
>>  pointSet -> SetPoint ( i , point );
>>  VectorType V;
>>  V [0] = 0.25* cos (t *6.0*3.141);
>>  V [1] = 0.25* sin (t *6.0*3.141);
>>  V [2] = 4.0* t;
>>  pointSet -> SetPointData ( i , V );
>>  }
>>
>> // Instantiate the filter and set the parameters
>> typedef itk :: BSplineScatteredDataPointSetToImageFilter
>>    < PointSetType , ImageType > FilterType ;
>> FilterType :: Pointer filter = FilterType :: New ();
>> // Define the parametric domain
>> ImageType :: SpacingType spacing ; spacing . Fill ( 0.001 );
>> ImageType :: SizeType size ; size . Fill ( 1.0/ spacing [0] );
>> ImageType :: PointType origin ; origin . Fill ( 0.0 );
>> filter -> SetSize ( size );
>> filter -> SetOrigin ( origin );
>> filter -> SetSpacing ( spacing );
>> filter -> SetInput ( pointSet );
>> filter -> SetSplineOrder ( 3 );
>> FilterType :: ArrayType ncps ;
>> ncps . Fill ( 4 );
>> filter -> SetNumberOfControlPoints ( ncps );
>> filter -> SetNumberOfLevels ( 5 );
>> filter -> SetGenerateOutputImage ( false );
>> try
>>  {
>>  filter -> Update ();
>>  }
>> catch (itk::ExceptionObject &err)
>>  {
>>  std :: cerr << " Test 2: itkBSplineScatteredDataImageFilter exception thrown " << std :: endl ;
>>  std::cerr << err << std::endl;
>>  return EXIT_FAILURE ;
>>  }
>> for ( RealType t = 0.0; t <= 1.0+1e-10; t += 0.01 )
>>  {
>>  PointSetType :: PointType point ;
>>  point [0] = t;
>>  VectorType V;
>>  filter -> Evaluate ( point , V );
>> }
>>
>> return 0 ;
>> }
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.html
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>     
>
>