[Insight-users] Fwd: Question regarding BSplineScatteredDataPointSetToImageFilter
Nicholas Tustison
ntustison at gmail.com
Tue Jan 10 10:24:29 EST 2012
Presumably your points over all time exist within a bounding
box. Do you know what that bounding box is?
If not, you might want to look at the itk::PointSet class and
it has a bounding box
On Jan 10, 2012, at 5:49 AM, Kerstin Müller wrote:
> Hi,
>
>
> than may be I set my origin wrong. Here the sample code with comments.
>
> typedef float RealType;
>
> // It's a volume
> typedef itk::Vector<RealType, 3> VectorType;
> // The deformation is over time
> typedef itk ::Image <VectorType , 4 > TimeVaryingDeformationFieldType;
> typedef itk ::PointSet <VectorType , 4 > TimeVaryingDeformationFieldPointSetType;
>
> TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints = TimeVaryingDeformationFieldPointSetType:: New ();
> fieldPoints->Initialize();
>
> TimeVaryingDeformationFieldType::PointType origin;
> TimeVaryingDeformationFieldType::SpacingType spacing;
> TimeVaryingDeformationFieldType::SizeType size;
> // I want to have a sapcing of 10 mm in each dimension except in the temporal there it is 2
> // Define the parametric domain of BSpline Field
> spacing[0] = 10.0f;
> spacing[1] = 10.0f;
> spacing[2] = 10.0f;
> spacing[3] = 2.0f;
> // Here the size of the BSpline field, I have 133 numberOfTimePoints
> size[0] = 128;
> size[1] = 128;
> size[2] = 128;
> size[3] = 67;
>
>
> // The origin ( where I'm not sure, I also tried 0,0,0)
> float O[3] ={0.0f,0.0f,0.0f};
> for (int i=0; i<3; i++)
> O[i] = m_configuration.origin[i] - spacing[i] * (size[i] - 1) *0.5f;
>
> origin[0] = O[0];
> origin[1] = O[1];
> origin[2] = O[2];
> origin[3] = 0.0f;
>
> // Set the sample points
> for ( int t = 0; t < numberOfTimePoints ; t++ )
> {
> int num = 0;
> for( int idx = 0; idx < endo_points.at(t).size()-3; idx+=3 )
> {
>
> VectorType v;
> // The sampling points range from -70 to +70 for example
> v[0] = endo_points.at(t)[idx] - endo_points.at(0)[idx] ;
> v[1] = endo_points.at(t)[idx+1] - endo_points.at(0)[idx+1];
> v[2] = endo_points.at(t)[idx+2] - endo_points.at(0)[idx+2];
>
>
> TimeVaryingDeformationFieldPointSetType::PointType point;
> point[0] = endo_points.at(0)[idx];
> point[1] = endo_points.at(0)[idx+1];
> point[2] = endo_points.at(0)[idx+2];
> point[3] = t;
>
> fieldPoints->SetPoint( num, point );
> fieldPoints->SetPointData(num, v);
> num++;
>
> }
> }
>
>
>
>
> // Define the filter and set the parameters
> typedef itk::myBSplineScatteredDataPointSetToImageFilter<TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;
> BSplineFilterType::Pointer filter = BSplineFilterType::New();
>
> filter ->SetSize( size );
> filter ->SetOrigin( origin );
> filter ->SetSpacing( spacing );
> filter ->SetGenerateOutputImage ( false );
> filter ->SetNumberOfLevels ( 2 );
> filter ->SetSplineOrder ( 3 );
>
> BSplineFilterType :: ArrayType ncps;
> ncps.SetElement(0,8);
> ncps.SetElement(1,8);
> ncps.SetElement(2,8);
> ncps.SetElement(3,5);
> filter -> SetNumberOfControlPoints ( ncps );
>
> filter ->SetInput( fieldPoints );
>
> filter->SetDebug(true);
>
> try
> {
> filter ->Update ();
> }
> catch ( itk::ExceptionObject & excp )
> {
> std :: cerr << "Test 2: itkBSplineScatteredDataImageFilter exception thrown" << std:: endl;
> std::cerr<< excp <<std::endl;
> }
>
> That's the Spline definition and afterwards I wanna resample the BSpline at arbitrary positions:
>
> float fPoint[3] = {0.0f,0.0f,0.0f};
>
> float sample_O[3] ={0.0f,0.0f,0.0f};
> for (int i=0; i<3; i++)
> sample_O[i] = m_configuration.origin[i] - m_configuration.vs * (m_configuration.size[i] - 1) *0.5f;
>
> int sample_size[3];
> sample_size[0] = m_configuration.size[0];
> sample_size[1] = m_configuration.size[1];
> sample_size[2] = m_configuration.size[2];
>
> for (int iz=0, idx=0; iz<sample_size[2]; iz++)
> {
> fPoint[2] = sample_O[2] + iz*m_configuration.vs;
> for (int iy=0; iy<sample_size[1]; iy++)
> {
> fPoint[1] = sample_O[1] + iy*m_configuration.vs;
> for (int ix=0; ix<sample_size[0]; ix++)
> {
> fPoint[0] = sample_O[0] + ix*m_configuration.vs;
>
> TimeVaryingDeformationFieldPointSetType :: PointType point;
> point [0] = fPoint[0];
> point [1] = fPoint[1];
> point [2] = fPoint[2];
>
> VectorType V;
> filter->EvaluateAtPoint(point, V);
> }
> //cout << iy << endl;
> }
> cout << iz << endl;
> }
>
> Thans a lot for your help!
>
> Best,
>
> Kerstin
> 2012/1/10 Nicholas Tustison <ntustison at gmail.com>
>> That's fine. I'm assuming, though, that they reside
>> in a bounding box of some sort, correct? Just
>> assign those points to their continuous coordinates
>> inside that bounding box.
>>
>>
>>
>>
>> On Jan 10, 2012, at 5:31 AM, Kerstin Müller wrote:
>>
>>
>> But my sample points are on continous world coordinates between for example -70.0 and 70.0.
>>
>> Best,
>>
>> Kerstin
>>
>> 2012/1/10 Nicholas Tustison <ntustison at gmail.com>
>>> Since your points reside in an image, an easy
>>> parameterization is to assign their location to
>>> their x,y,z location in the image. Does that make
>>> sense?
>>>
>>>
>>> On Jan 10, 2012, at 5:19 AM, Kerstin Müller wrote:
>>>
>>>
>>>
>>> Hi,
>>>
>>> yes I also used your sample code, but still I got a heap corruption error. That's why I just wanted to try if my sample points are the problem.
>>> Do I need to scale them to a range between [0 1]?
>>>
>>> Best Regards,
>>>
>>> Kerstin
>>>
>>> 2012/1/10 Nicholas Tustison <ntustison at gmail.com>
>>>> Hi Kerstin,
>>>>
>>>> Please keep the correspondence on the user's list.
>>>>
>>>> Why are you setting your vector to all zeros?
>>>>
>>>>
>>>> VectorType v;
>>>> v[0] = 0.0;
>>>> v[1] = 0.0;
>>>> v[2] = 0.0;
>>>> and why are you setting the location of your points to
>>>> (idx,idx,idx,t)?
>>>>
>>>>
>>>> point[0] = idx;
>>>> point[1] = idx;
>>>> point[2] = idx;
>>>> point[3] = t;
>>>> That can't possibly be right. Did you look at the sample
>>>> code I gave you?
>>>>
>>>> Nick
>>>>
>>>>
>>>>
>>>> On Jan 9, 2012, at 2:04 AM, Kerstin Müller wrote:
>>>>
>>>>
>>>> Hi,
>>>>
>>>> thank you very much!
>>>> However I'm still struggeling with a heap corruption problem. I'm not sure if I set my parameters correct.
>>>>
>>>> typedef float RealType;
>>>>
>>>> typedef itk::Vector<RealType, 3> VectorType;
>>>> typedef itk ::Image <VectorType , 4 > TimeVaryingDeformationFieldType;
>>>> typedef itk :: PointSet <VectorType , 4 > TimeVaryingDeformationFieldPointSetType;
>>>>
>>>> TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints = TimeVaryingDeformationFieldPointSetType:: New ();
>>>> fieldPoints->Initialize();
>>>>
>>>> // Set the sample points just for testing 10 linear points no movement
>>>> for ( int t = 0; t < m_configuration.num_projections ; t++ )
>>>> {
>>>> int num = 0;
>>>> for( int idx = 0; idx <= 10; idx++ )
>>>> {
>>>>
>>>> VectorType v;
>>>> v[0] = 0.0;
>>>> v[1] = 0.0;
>>>> v[2] = 0.0;
>>>>
>>>>
>>>> TimeVaryingDeformationFieldPointSetType::PointType point;
>>>> point[0] = idx;
>>>> point[1] = idx;
>>>> point[2] = idx;
>>>> point[3] = t;
>>>>
>>>> fieldPoints->SetPoint( num, point );
>>>> fieldPoints->SetPointData(num, v);
>>>> num++;
>>>>
>>>> }
>>>> }
>>>>
>>>>
>>>> TimeVaryingDeformationFieldType::PointType origin;
>>>> TimeVaryingDeformationFieldType::SpacingType spacing;
>>>> TimeVaryingDeformationFieldType::SizeType size;
>>>>
>>>> // Define the parametric domain
>>>> spacing[0] = 1.0f;
>>>> spacing[1] = 1.0f;
>>>> spacing[2] = 1.0f;
>>>> spacing[3] = 2.0f;
>>>>
>>>> size[0] = 128;
>>>> size[1] = 128;
>>>> size[2] = 128;
>>>> size[3] = 66;
>>>>
>>>> origin[0] = 0.0f;
>>>> origin[1] = 0.0f;
>>>> origin[2] = 0.0f;
>>>> origin[3] = 0.0f;
>>>>
>>>> // Define the filter and set the parameters
>>>> typedef itk::myBSplineScatteredDataPointSetToImageFilter<TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;
>>>> BSplineFilterType::Pointer filter = BSplineFilterType::New();
>>>>
>>>> filter ->SetSize( size );
>>>> filter ->SetOrigin( origin );
>>>> filter ->SetSpacing( spacing );
>>>> filter -> SetGenerateOutputImage ( false );
>>>> filter -> SetNumberOfLevels ( 2 );
>>>> filter -> SetSplineOrder ( 3 );
>>>>
>>>>
>>>> BSplineFilterType :: ArrayType ncps;
>>>> ncps.SetElement(0,8);
>>>> ncps.SetElement(1,8);
>>>> ncps.SetElement(2,8);
>>>> ncps.SetElement(3,34);
>>>> filter -> SetNumberOfControlPoints ( ncps );
>>>>
>>>> filter ->SetInput( fieldPoints );
>>>>
>>>> filter->SetDebug(true);
>>>>
>>>> try
>>>> {
>>>> filter ->Update ();
>>>> }
>>>> catch ( itk::ExceptionObject & excp )
>>>> {
>>>> std :: cerr << "Test 2: itkBSplineScatteredDataImageFilter exception thrown" << std:: endl;
>>>> std::cerr<< excp <<std::endl;
>>>> }
>>>>
>>>> Thank you very much!
>>>>
>>>> Best,
>>>>
>>>> Kerstin
>>>> 2012/1/4 Nicholas Tustison <ntustison at gmail.com>
>>>>> Okay, here's some pseudocode which hasn't been tested
>>>>> since you have to provide it the interface to your values
>>>>> but you should get the general idea.
>>>>>
>>>>> typedef float RealType;
>>>>> typedef itk::Image<RealType, ImageDimension> RealImageType;
>>>>>
>>>>> // Read in your reference image which will define the domain of
>>>>> // your B-spline displacement field. Assume, for simplicity, that it
>>>>> // has identity direction
>>>>>
>>>>> typedef itk::ImageFileReader<RealImageType> ImageReaderType;
>>>>> typename ImageReaderType::Pointer reader = ImageReaderType::New();
>>>>> reader->SetFileName( XXXX );
>>>>> reader->Update();
>>>>>
>>>>> typedef itk::Vector<RealType, ImageDimension> VectorType;
>>>>> typedef itk::Image<VectorType, ImageDimension+1> TimeVaryingDeformationFieldType;
>>>>>
>>>>> typedef itk::PointSet<VectorType, ImageDimension+1> TimeVaryingDeformationFieldPointSetType;
>>>>> typename TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints =
>>>>> TimeVaryingDeformationFieldPointSetType::New();
>>>>> fieldPoints->Initialize();
>>>>> unsigned long count = 0;
>>>>>
>>>>> // Assume points are stored in a 2-D matrix called 'mySamplePoints' where the column
>>>>> // is the time point and the row is the point
>>>>>
>>>>> for( unsigned int t = 0; t < numberOfTimePoints; t++ )
>>>>> {
>>>>> for( unsigned int n = 0; n < numberOfPointsPerTimePoint; n++ )
>>>>> {
>>>>> VectorType displacement = mySamplePoints[n][t] - mySamplePoints[n][0];
>>>>>
>>>>> TimeVaryingDeformationFieldPointSetType::PointType parametric point;
>>>>> for( unsigned int d = 0; d < ImageDimension; d++ )
>>>>> {
>>>>> point[d] = ( mySamplePoints[n][0] )[d];
>>>>> }
>>>>> point[ImageDimension] = t;
>>>>>
>>>>> fieldPoints->SetPoint( count, point );
>>>>> fieldPoints->SetPointData( count, displacement );
>>>>> count++;
>>>>> }
>>>>> }
>>>>>
>>>>> TimeVaryingDeformationFieldType::PointType origin;
>>>>> TimeVaryingDeformationFieldType::SpacingType spacing;
>>>>> TimeVaryingDeformationFieldType::SizeType size;
>>>>>
>>>>> for( unsigned int d = 0; d < ImageDimension; d++ )
>>>>> {
>>>>> origin[d] = reader->GetOutput()->GetOrigin()[d];
>>>>> size[d] = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[d];
>>>>> spacing[d] = reader->GetOutput()->GetSpacing()[d];
>>>>> }
>>>>> // Now include the temporal information. You can change this to whatever
>>>>> // resolution you like. You just need to make sure that
>>>>> // ( size[ImageDimension] - 1 ) * spacing[ImageDimension] = ( numberOfTimePoints - 1)
>>>>> //
>>>>> origin[ImageDimension] = 0;
>>>>> size[ImageDimension] = numberOfTimePoints - 1;
>>>>> spacing[ImageDimension] = 1;
>>>>>
>>>>> typedef itk::BSplineScatteredDataPointSetToImageFilter
>>>>> <TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;
>>>>> typename BSplineFilterType::Pointer bspliner = BSplineFilterType::New();
>>>>> bspliner->SetOrigin( origin );
>>>>> bspliner->SetSpacing( spacing );
>>>>> bspliner->SetSize( size );
>>>>> bspliner->SetGenerateOutputImage( XXXX );
>>>>> bspliner->SetNumberOfLevels( XXXX );
>>>>> bspliner->SetSplineOrder( XXXX );
>>>>> bspliner->SetNumberOfControlPoints( XXXX );
>>>>> bspliner->SetInput( fieldPoints );
>>>>> bspliner->Update();
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Jan 4, 2012, at 2:52 AM, Kerstin Müller wrote:
>>>>>
>>>>>
>>>>> Hi,
>>>>>
>>>>> yes correctly.
>>>>>
>>>>> 2012/1/3 Nicholas Tustison <ntustison at gmail.com>
>>>>>> Okay, I'm assuming that these points correspond in
>>>>>> time, correct? For example, point 937 in time point 0
>>>>>> corresponds to point 937 in time point 1, 2, 3, ...133,
>>>>>> right?
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Jan 3, 2012, at 11:55 AM, Kerstin Müller wrote:
>>>>>>
>>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> the basic problem I wanna solve:
>>>>>>
>>>>>> I'll have scattered points in 3D which vary over time, they describe a motion of a surface over time. They are not ordered in x-y-and z-dimension.
>>>>>> Now I want to represent that motion by BSplines in order to generate a dense motion vector field defined on specific 3D voxel positions.
>>>>>> I'll have 133 timesteps and ~960 sample points in each time step. Now I want to fit the BSpline to it and resample it.
>>>>>>
>>>>>> All the best,
>>>>>>
>>>>>> Kerstin
>>>>>>
>>>>>> 2012/1/3 Nicholas Tustison <ntustison at gmail.com>
>>>>>>> In that case, let's start with the basic problem set-up.
>>>>>>> Before, you described your problem as follows:
>>>>>>>
>>>>>>>
>>>>>>>> I'm using the BSplineScatteredDataPointSetToImageFilter from the Insight Journal. However,
>>>>>>>> I cannot find the correct parameter in order to make the filter work and does not crash during the evaluation.
>>>>>>>> I'll have scattered 3-D points over time and I want to fit a 4-D Bspline field to that points. Afterwards I want to evaluate the Bspline on a dense volume grid to one time step.
>>>>>>>
>>>>>>> What do these scattered 3D+t points represent?
>>>>>>> Is it a curve, a time-varying scalar field, or something
>>>>>>> else?
>>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120110/2ff06a98/attachment.htm>
More information about the Insight-users
mailing list