[Insight-users] Question regarding BSplineScatteredDataPointSetToImageFilter
Kerstin Müller
kerstin.mueller612 at googlemail.com
Tue Jan 10 07:49:10 EST 2012
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/4623c51f/attachment-0001.htm>
More information about the Insight-users
mailing list