[Insight-users] itkBSplineScatteredDataImageFilter exception thrown

[Ian Malone]

Thanks for your help, think I understand what's happening now.

imalone

Nicholas Tustison wrote:
> Sorry, let me clarify.  So when Update() is called the filter does the fitting and the result is a set of B-spline control points which define a B-spline object over the continuous domain which the user defines via the SetSpacing, SetOrigin, and SetSize functions.  The exception is thrown typically because one of your points maps outside that domain which is what I think you were getting at with the first comment.  So what you want to do is make your continuous domain a bit larger by increasing the spacing or size or both.  
>
> I apologize for the code as it is probably a bit outdated.  I'm sure the set-up is still the same, it's just that many of these issues have been fixed or addressed in the code that is actually in ITK without propagating those changes back to the IJ source code.
>
> Nick
>
>
>
> On Dec 14, 2010, at 12:01 PM, Ian Malone wrote:
>
>   
>> 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
>>>>    
>>>>         
>>>  
>>>       
>
>   

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20101215/721106a5/attachment.htm>