Hi,<br><br><br>than may be I set my origin wrong. Here the sample code with comments.<br><br> typedef float RealType;<br> <br>// It's a volume<br> typedef itk::Vector<RealType, 3> VectorType;<br>// The deformation is over time<br>
typedef itk ::Image <VectorType , 4 > TimeVaryingDeformationFieldType;<br> typedef itk ::PointSet <VectorType , 4 > TimeVaryingDeformationFieldPointSetType;<br><br> TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints = TimeVaryingDeformationFieldPointSetType:: New ();<br>
fieldPoints->Initialize();<br><br> TimeVaryingDeformationFieldType::PointType origin;<br> TimeVaryingDeformationFieldType::SpacingType spacing;<br> TimeVaryingDeformationFieldType::SizeType size;<br>
// I want to have a sapcing of 10 mm in each dimension except in the temporal there it is 2 <br> // Define the parametric domain of BSpline Field<br> spacing[0] = 10.0f;<br> spacing[1] = 10.0f;<br> spacing[2] = 10.0f;<br>
spacing[3] = 2.0f;<br>// Here the size of the BSpline field, I have 133 numberOfTimePoints<br> size[0] = 128;<br> size[1] = 128;<br> size[2] = 128;<br> size[3] = 67;<br> <br> <br>// The origin ( where I'm not sure, I also tried 0,0,0)<br>
float O[3] ={0.0f,0.0f,0.0f};<br> for (int i=0; i<3; i++)<br> O[i] = m_configuration.origin[i] - spacing[i] * (size[i] - 1) *0.5f;<br><br> origin[0] = O[0];<br> origin[1] = O[1];<br> origin[2] = O[2];<br>
origin[3] = 0.0f;<br><br> // Set the sample points<br> for ( int t = 0; t < numberOfTimePoints ; t++ ) <br> {<br> int num = 0;<br> for( int idx = 0; idx < <a href="http://endo_points.at">endo_points.at</a>(t).size()-3; idx+=3 )<br>
{<br> <br> VectorType v;<br>// The sampling points range from -70 to +70 for example<br> v[0] = <a href="http://endo_points.at">endo_points.at</a>(t)[idx] - <a href="http://endo_points.at">endo_points.at</a>(0)[idx] ;<br>
v[1] = <a href="http://endo_points.at">endo_points.at</a>(t)[idx+1] - <a href="http://endo_points.at">endo_points.at</a>(0)[idx+1];<br> v[2] = <a href="http://endo_points.at">endo_points.at</a>(t)[idx+2] - <a href="http://endo_points.at">endo_points.at</a>(0)[idx+2];<br>
<br> <br> TimeVaryingDeformationFieldPointSetType::PointType point;<br> point[0] = <a href="http://endo_points.at">endo_points.at</a>(0)[idx];<br> point[1] = <a href="http://endo_points.at">endo_points.at</a>(0)[idx+1];<br>
point[2] = <a href="http://endo_points.at">endo_points.at</a>(0)[idx+2];<br> point[3] = t;<br><br> fieldPoints->SetPoint( num, point );<br> fieldPoints->SetPointData(num, v);<br>
num++;<br><br> }<br> }<br><br><br> <br><br> // Define the filter and set the parameters<br> typedef itk::myBSplineScatteredDataPointSetToImageFilter<TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;<br>
BSplineFilterType::Pointer filter = BSplineFilterType::New();<br><br> filter ->SetSize( size );<br> filter ->SetOrigin( origin );<br> filter ->SetSpacing( spacing );<br> filter ->SetGenerateOutputImage ( false );<br>
filter ->SetNumberOfLevels ( 2 );<br> filter ->SetSplineOrder ( 3 );<br><br> BSplineFilterType :: ArrayType ncps;<br> ncps.SetElement(0,8);<br> ncps.SetElement(1,8);<br> ncps.SetElement(2,8);<br> ncps.SetElement(3,5);<br>
filter -> SetNumberOfControlPoints ( ncps );<br> <br> filter ->SetInput( fieldPoints );<br><br> filter->SetDebug(true);<br><br> try<br> {<br> filter ->Update ();<br> }<br> catch ( itk::ExceptionObject & excp )<br>
{<br> std :: cerr << "Test 2: itkBSplineScatteredDataImageFilter exception thrown" << std:: endl;<br> std::cerr<< excp <<std::endl;<br> }<br><br>That's the Spline definition and afterwards I wanna resample the BSpline at arbitrary positions:<br>
<br>float fPoint[3] = {0.0f,0.0f,0.0f};<br><br> float sample_O[3] ={0.0f,0.0f,0.0f};<br> for (int i=0; i<3; i++)<br> sample_O[i] = m_configuration.origin[i] - m_configuration.vs * (m_configuration.size[i] - 1) *0.5f;<br>
<br> int sample_size[3];<br> sample_size[0] = m_configuration.size[0];<br> sample_size[1] = m_configuration.size[1];<br> sample_size[2] = m_configuration.size[2];<br> <br> for (int iz=0, idx=0; iz<sample_size[2]; iz++)<br>
{<br> fPoint[2] = sample_O[2] + iz*m_configuration.vs;<br> for (int iy=0; iy<sample_size[1]; iy++)<br> {<br> fPoint[1] = sample_O[1] + iy*m_configuration.vs;<br> for (int ix=0; ix<sample_size[0]; ix++)<br>
{<br> fPoint[0] = sample_O[0] + ix*m_configuration.vs;<br> <br> TimeVaryingDeformationFieldPointSetType :: PointType point;<br> point [0] = fPoint[0];<br>
point [1] = fPoint[1];<br> point [2] = fPoint[2];<br><br> VectorType V;<br> filter->EvaluateAtPoint(point, V);<br> }<br> //cout << iy << endl;<br>
}<br> cout << iz << endl;<br> }<br><br>Thans a lot for your help!<br><br>Best,<br><br>Kerstin<br><div class="gmail_quote">2012/1/10 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com">ntustison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">That's fine. I'm assuming, though, that they reside<div>in a bounding box of some sort, correct? Just </div>
<div>assign those points to their continuous coordinates</div><div>inside that bounding box.</div><div><div class="h5"><div><br></div><div><br><div><div><br></div><div><br><div><div>On Jan 10, 2012, at 5:31 AM, Kerstin Müller wrote:</div>
<br><blockquote type="cite">But my sample points are on continous world coordinates between for example -70.0 and 70.0.<br><br>Best,<br><br>Kerstin<br><br><div class="gmail_quote">2012/1/10 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">Since your points reside in an image, an easy<div>parameterization is to assign their location to </div>
<div>their x,y,z location in the image. Does that make</div><div>sense?<div><div><br><div><div><br><div><div>On Jan 10, 2012, at 5:19 AM, Kerstin Müller wrote:</div><br><blockquote type="cite"><br>Hi,<br><br>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.<br>
Do I need to scale them to a range between [0 1]?<br><br>
Best Regards,<br><br>Kerstin<br><br><div class="gmail_quote">2012/1/10 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word">Hi Kerstin,<div><br></div><div>Please keep the correspondence on the user's list. </div><div><br></div><div>Why are you setting your vector to all zeros? </div><div><div><br>
</div><div><blockquote type="cite"> VectorType v;<br> v[0] = 0.0;<br> v[1] = 0.0;<br> v[2] = 0.0;<br></blockquote><br></div></div><div>and why are you setting the location of your points to </div>
<div>(idx,idx,idx,t)?</div><div><div><br></div><div><blockquote type="cite"> point[0] = idx;<br> point[1] = idx;<br> point[2] = idx;<br> point[3] = t;<br></blockquote>
<br></div></div><div>That can't possibly be right. Did you look at the sample</div><div>code I gave you?</div><div><br></div><div>Nick</div><div><div><div><br></div><div><br></div><div><br><div><div>On Jan 9, 2012, at 2:04 AM, Kerstin Müller wrote:</div>
<br><blockquote type="cite">Hi,<br><br>thank you very much!<br>However I'm still struggeling with a heap corruption problem. I'm not sure if I set my parameters correct.<br><br> typedef float RealType;<br> <br>
typedef itk::Vector<RealType, 3> VectorType;<br>
typedef itk ::Image <VectorType , 4 > TimeVaryingDeformationFieldType;<br> typedef itk :: PointSet <VectorType , 4 > TimeVaryingDeformationFieldPointSetType;<br><br> TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints = TimeVaryingDeformationFieldPointSetType:: New ();<br>
fieldPoints->Initialize();<br><br> // Set the sample points just for testing 10 linear points no movement<br> for ( int t = 0; t < m_configuration.num_projections ; t++ ) <br> {<br> int num = 0;<br>
for( int idx = 0; idx <= 10; idx++ )<br> {<br> <br> VectorType v;<br> v[0] = 0.0;<br> v[1] = 0.0;<br> v[2] = 0.0;<br><br> <br> TimeVaryingDeformationFieldPointSetType::PointType point;<br>
point[0] = idx;<br> point[1] = idx;<br> point[2] = idx;<br> point[3] = t;<br><br> fieldPoints->SetPoint( num, point );<br> fieldPoints->SetPointData(num, v);<br>
num++;<br><br> }<br> }<br><br><br> TimeVaryingDeformationFieldType::PointType origin;<br> TimeVaryingDeformationFieldType::SpacingType spacing;<br> TimeVaryingDeformationFieldType::SizeType size;<br>
<br> // Define the parametric domain<br> spacing[0] = 1.0f;<br> spacing[1] = 1.0f;<br> spacing[2] = 1.0f;<br> spacing[3] = 2.0f;<br><br> size[0] = 128;<br> size[1] = 128;<br> size[2] = 128;<br> size[3] = 66;<br>
<br> origin[0] = 0.0f;<br> origin[1] = 0.0f;<br> origin[2] = 0.0f;<br> origin[3] = 0.0f;<br><br> // Define the filter and set the parameters<br> typedef itk::myBSplineScatteredDataPointSetToImageFilter<TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;<br>
BSplineFilterType::Pointer filter = BSplineFilterType::New();<br><br> filter ->SetSize( size );<br> filter ->SetOrigin( origin );<br> filter ->SetSpacing( spacing );<br> filter -> SetGenerateOutputImage ( false );<br>
filter -> SetNumberOfLevels ( 2 );<br> filter -> SetSplineOrder ( 3 );<br><br> <br> BSplineFilterType :: ArrayType ncps;<br> ncps.SetElement(0,8);<br> ncps.SetElement(1,8);<br> ncps.SetElement(2,8);<br>
ncps.SetElement(3,34);<br> filter -> SetNumberOfControlPoints ( ncps );<br> <br> filter ->SetInput( fieldPoints );<br><br> filter->SetDebug(true);<br><br> try<br> {<br> filter ->Update ();<br>
}<br> catch ( itk::ExceptionObject & excp )<br> {<br> std :: cerr << "Test 2: itkBSplineScatteredDataImageFilter exception thrown" << std:: endl;<br> std::cerr<< excp <<std::endl;<br>
}<br><br>Thank you very much!<br><br>Best,<br><br>Kerstin<br><div class="gmail_quote">2012/1/4 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div>Okay, here's some pseudocode which hasn't been tested</div><div>since you have to provide it the interface to your values </div>
<div>but you should get the general idea.</div><div><br></div><div> typedef float RealType;</div><div><div> typedef itk::Image<RealType, ImageDimension> RealImageType;</div><div><br></div><div> // Read in your reference image which will define the domain of </div>
<div> // your B-spline displacement field. Assume, for simplicity, that it</div><div> // has identity direction</div><div><br></div><div> typedef itk::ImageFileReader<RealImageType> ImageReaderType;</div><div> typename ImageReaderType::Pointer reader = ImageReaderType::New();</div>
<div> reader->SetFileName( XXXX );</div><div> reader->Update();</div><div><br></div><div> typedef itk::Vector<RealType, ImageDimension> VectorType;</div><div> typedef itk::Image<VectorType, ImageDimension+1> TimeVaryingDeformationFieldType;</div>
<div><br></div><div> typedef itk::PointSet<VectorType, ImageDimension+1> TimeVaryingDeformationFieldPointSetType;</div></div><div> typename TimeVaryingDeformationFieldPointSetType::Pointer fieldPoints =</div><div>
TimeVaryingDeformationFieldPointSetType::New();</div><div> fieldPoints->Initialize();</div><div> unsigned long count = 0;</div><div><br></div><div> // Assume points are stored in a 2-D matrix called 'mySamplePoints' where the column</div>
<div> // is the time point and the row is the point</div><div><br></div><div> for( unsigned int t = 0; t < numberOfTimePoints; t++ )</div><div> {</div><div> for( unsigned int n = 0; n < numberOfPointsPerTimePoint; n++ )</div>
<div> {</div><div> VectorType displacement = mySamplePoints[n][t] - mySamplePoints[n][0];</div><div> </div><div> TimeVaryingDeformationFieldPointSetType::PointType parametric point;</div><div> for( unsigned int d = 0; d < ImageDimension; d++ )</div>
<div> {</div><div> point[d] = ( mySamplePoints[n][0] )[d];</div><div> }</div><div> point[ImageDimension] = t;</div><div><br></div><div> fieldPoints->SetPoint( count, point );</div><div> fieldPoints->SetPointData( count, displacement );</div>
<div> count++;</div><div> }</div><div> } </div><div><br></div><div><div> TimeVaryingDeformationFieldType::PointType origin;</div></div><div><div> TimeVaryingDeformationFieldType::SpacingType spacing;</div>
</div><div><div> TimeVaryingDeformationFieldType::SizeType size;</div></div><div><br></div><div> for( unsigned int d = 0; d < ImageDimension; d++ )</div><div> {</div><div> origin[d] = reader->GetOutput()->GetOrigin()[d];</div>
<div> size[d] = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[d];</div><div> spacing[d] = reader->GetOutput()->GetSpacing()[d];</div><div> }</div><div> // Now include the temporal information. You can change this to whatever</div>
<div> // resolution you like. You just need to make sure that </div><div> // ( size[ImageDimension] - 1 ) * spacing[ImageDimension] = ( numberOfTimePoints - 1)</div><div> // </div><div> origin[ImageDimension] = 0;</div>
<div> size[ImageDimension] = numberOfTimePoints - 1;</div><div> spacing[ImageDimension] = 1;</div><div><br></div><div><div> typedef itk::BSplineScatteredDataPointSetToImageFilter</div><div> <TimeVaryingDeformationFieldPointSetType, TimeVaryingDeformationFieldType> BSplineFilterType;</div>
<div> typename BSplineFilterType::Pointer bspliner = BSplineFilterType::New();</div></div><div><div> bspliner->SetOrigin( origin );</div><div> bspliner->SetSpacing( spacing );</div><div> bspliner->SetSize( size );</div>
<div> bspliner->SetGenerateOutputImage( XXXX );</div><div> bspliner->SetNumberOfLevels( XXXX );</div><div> bspliner->SetSplineOrder( XXXX );</div><div> bspliner->SetNumberOfControlPoints( XXXX );</div><div>
bspliner->SetInput( fieldPoints );</div></div><div> bspliner->Update();</div><div><div><div><br></div><div><br></div><div><br></div><br><div><div>On Jan 4, 2012, at 2:52 AM, Kerstin Müller wrote:</div>
<br><blockquote type="cite">Hi,<br><br>yes correctly.<br><br><div class="gmail_quote">2012/1/3 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word"><div>Okay, I'm assuming that these points correspond in </div><div>time, correct? For example, point 937 in time point 0</div><div>corresponds to point 937 in time point 1, 2, 3, ...133,</div>
<div>right? </div><div><div><div><br></div><div><br></div><br><div><div>On Jan 3, 2012, at 11:55 AM, Kerstin Müller wrote:</div><br><blockquote type="cite">Hi,<br><br>the basic problem I wanna solve:<br><br>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.<br>
Now I want to represent that motion by BSplines in order to generate a dense motion vector field defined on specific 3D voxel positions.<br>
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.<br><br>All the best,<br><br>Kerstin<br><br><div class="gmail_quote">2012/1/3 Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">In that case, let's start with the basic problem set-up.<div>Before, you described your problem as follows:<div>
<br><div><div><br><blockquote type="cite"><span style="border-collapse:separate;font-family:Helvetica;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:normal;text-align:-webkit-auto;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;font-size:medium"><blockquote type="cite">
<div><div>I'm using the<span> </span><b>BSplineScatteredDataPointSetToImageFilter<span> </span></b>from the Insight Journal. However,<br>I cannot find the correct parameter in order to make the filter work and does not crash during the evaluation.<br>
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.</div></div></blockquote></span></blockquote>
</div></div></div></div><div><br></div><div>What do these scattered 3D+t points represent? </div><div>Is it a curve, a time-varying scalar field, or something</div><div>else?</div><div><br></div></div></blockquote></div>
<br>
</blockquote></div><br></div></div></div></blockquote></div><br>
</blockquote></div><br></div></div></div></blockquote></div><br>
</blockquote></div><br></div></div></div></div></blockquote></div><br>
</blockquote></div><br></div></div></div></div></div></div></blockquote></div><br>
</blockquote></div><br></div></div></div></div></div></div></blockquote></div><br>