Re: [ITK-users] DVFs interpolation by BSpline
Dženan Zukić <[email protected]> Wed, 3 Oct 2018 14:10:44 -0400
| Newsgroups | gmane.comp.lib.itk.user |
|---|---|
| Message-ID | <CAPf2UMS9gEwOP1BWnu5ctuFG_MdJYFA6JMWLH1zhmEYYPGwc8w@mail.gmail.com> |
--===============0995632974== Content-Type: multipart/alternative; boundary="000000000000b86178057756f5e1" --000000000000b86178057756f5e1 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable Hi JongWon, welcome to ITK! To reproduce your problem we will need the input file too. If it crashes with different input files, please provide the smallest one. Also, we have migrated to the forum <https://discourse.itk.org/>, so please post the updated question there. Regards, D=C5=BEenan On Wed, Oct 3, 2018 at 11:30 AM wonnybro <[email protected]> wrote: > Hi, ITK users > > I really really need your help, please help me... > I am trying to do BSpline interpolation of DVFs obtained by image > registrations of 10 images. > I'm trying to modify the BSplineScatteredDataPointSetToImageFilter from t= he > Insight Journal. However, > I am having trouble to find the correct parameter in order to make the > filter work. > I'll have scattered 3-D points of several DVFs over time and I want to fi= t > these DVFs by a 4-D Bspline. Afterwards I want to evaluate output DVF to > one > time step. > > Here you can find the code, unfortunately it crashes. > > Please help me...please. > > Thanks in advance. > > #include "itkBSplineScatteredDataPointSetToImageFilter.h" > #include "itkPointSet.h" > #include "itkImage.h" > #include "itkVectorImage.h" > #include "itkImageFileReader.h" > #include "itkImageFileWriter.h" > > int main (int argc, char* argv[]) > { > const unsigned int ImageDimension =3D 3; > > using VectorType =3D itk::Vector<float, ImageDimension>; > using InputImageType =3D itk::Image<VectorType, ImageDimension>; > using InputPointType =3D InputImageType::PointType; > using OutputImageType =3D itk::Image<VectorType, ImageDimension + 1>; > using PointSetType =3D itk::PointSet<VectorType, ImageDimension + 1>; > using OutputPointType =3D PointSetType::PointType; > using ReaderType =3D itk::ImageFileReader<InputImageType>; > using WriterType =3D itk::ImageFileWriter<OutputImageType>; > > PointSetType::Pointer pointSet =3D PointSetType::New(); > unsigned long pointId =3D 0; > InputImageType::PointType origin; > InputImageType::SpacingType spacing; > InputImageType::SizeType size; > for (int i =3D 1; i < argc; ++i) > { > ReaderType::Pointer reader =3D ReaderType::New(); > std::cout << "Reading file " << argv[i] << std::endl; > reader->SetFileName(argv[i]); > try > { > reader->Update(); > } > catch (itk::ExceptionObject & err) > { > std::cout << "ExceptionObject caught !" << std::endl; > std::cout << err << std::endl; > return EXIT_FAILURE; > } > using IteratorType =3D itk::ImageRegionConstIterator<InputImageType>; > const InputImageType * image =3D reader->GetOutput(); > IteratorType it(image, image->GetBufferedRegion()); > it.GoToBegin(); > InputPointType inPoint; > OutputPointType outPoint; > int iCount =3D 0; > while (!it.IsAtEnd()) > { > inPoint =3D image->GetPixel(it.GetIndex()); > for (int j =3D 0; j < ImageDimension; ++j) > outPoint[j] =3D inPoint[i]; > outPoint[ImageDimension] =3D i - 1; > pointSet->SetPoint(pointId, outPoint); > // Transfer the pixel data to the value associated with the point. > pointSet->SetPointData(pointId, it.Get()); > ++it; > ++pointId; > ++iCount; > } > std::cout << "Number of points in " << argv[i] << " =3D " << iCount <= < > std::endl; > std::cout << "Total number of points =3D " << > pointSet->GetNumberOfPoints() << std::endl; > origin =3D image->GetOrigin(); > spacing =3D image->GetSpacing(); > size =3D image->GetLargestPossibleRegion().GetSize(); > } > > typedef itk::BSplineScatteredDataPointSetToImageFilter < PointSetType, > OutputImageType > SplineFilterType; > SplineFilterType::Pointer splineFilter =3D SplineFilterType::New(); > > int splineorder=3D3; // complexity of the spline > > SplineFilterType::ArrayType ncontrol; > ncontrol[0]=3Dsplineorder + 1; > SplineFilterType::ArrayType closedim; > closedim[0]=3D 0; > > OutputImageType::PointType parametricDomainOrigin; > OutputImageType::SpacingType parametricDomainSpacing; > OutputImageType::SizeType parametricDomainSize; > for (int i =3D 0; i < ImageDimension; ++i) > { > parametricDomainOrigin[i] =3D origin[i]; > parametricDomainSpacing[i] =3D spacing[i]; > parametricDomainSize[i] =3D size[i]; > } > parametricDomainOrigin[ImageDimension] =3D 0; > parametricDomainSize[ImageDimension] =3D argc - 2; > parametricDomainSpacing[ImageDimension] =3D 10.0; > > splineFilter->SetGenerateOutputImage( true ); // the only reason to > turn > this off is if one only wants to use the control point lattice for furthe= r > processing > splineFilter->SetInput ( pointSet ); > splineFilter->SetSplineOrder ( splineorder ); > splineFilter->SetNumberOfControlPoints ( ncontrol ); > splineFilter->SetNumberOfLevels( 3 ); > splineFilter->SetCloseDimension ( closedim ); > splineFilter->SetSize( parametricDomainSize ); > splineFilter->SetSpacing( parametricDomainSpacing ); > splineFilter->SetOrigin( parametricDomainOrigin ); > std::cout << "Before update spline filter" << std::endl; > splineFilter->Update(); > std::cout << "After update spline filter" << std::endl; > > WriterType::Pointer writer =3D WriterType::New(); > writer->SetInput(splineFilter->GetOutput()); > writer->SetFileName("Output.mhd"); > writer->Update(); > std::cout << "After write image filter" << std::endl; > > return EXIT_SUCCESS; > }; > > > > -- > Sent from: http://itk-users.7.n7.nabble.com/ > The ITK community is transitioning from this mailing list to > discourse.itk.org. Please join us there! > ________________________________ > 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.php > > Please keep messages on-topic and check the ITK FAQ at: > http://www.itk.org/Wiki/ITK_FAQ > > Follow this link to subscribe/unsubscribe: > https://itk.org/mailman/listinfo/insight-users > --000000000000b86178057756f5e1 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div class=3D"gmail_default" style=3D"font-family:verdana,= sans-serif;font-size:small">Hi JongWon,</div><div class=3D"gmail_default" s= tyle=3D"font-family:verdana,sans-serif;font-size:small"><br></div><div clas= s=3D"gmail_default" style=3D"font-family:verdana,sans-serif;font-size:small= ">welcome to ITK! To reproduce your problem we will need the input file too= . If it crashes with different input files, please provide the smallest one= . Also, we have migrated to the <a href=3D"https://discourse.itk.org/">foru= m</a>, so please post the updated question there.</div><div class=3D"gmail_= default" style=3D"font-family:verdana,sans-serif;font-size:small"><br></div= ><div class=3D"gmail_default" style=3D"font-family:verdana,sans-serif;font-= size:small">Regards,</div><div class=3D"gmail_default" style=3D"font-family= :verdana,sans-serif;font-size:small">D=C5=BEenan</div></div><br><div class= =3D"gmail_quote"><div dir=3D"ltr">On Wed, Oct 3, 2018 at 11:30 AM wonnybro = <<a href=3D"mailto:[email protected]">[email protected]</a>> wrote:= <br></div><blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;bord= er-left:1px #ccc solid;padding-left:1ex">Hi, ITK users<br> <br> I really really need your help, please help me...<br> I am trying to do BSpline interpolation of DVFs obtained by image<br> registrations of 10 images. <br> I'm trying to modify the BSplineScatteredDataPointSetToImageFilter from= the<br> Insight Journal. However,<br> I am having trouble to find the correct parameter in order to make the<br> filter work. <br> I'll have scattered 3-D points of several DVFs over time and I want to = fit<br> these DVFs by a 4-D Bspline. Afterwards I want to evaluate output DVF to on= e<br> time step.<br> <br> Here you can find the code, unfortunately it crashes.<br> <br> Please help me...please.<br> <br> Thanks in advance.<br> <br> #include "itkBSplineScatteredDataPointSetToImageFilter.h"<br> #include "itkPointSet.h"<br> #include "itkImage.h"<br> #include "itkVectorImage.h"<br> #include "itkImageFileReader.h"<br> #include "itkImageFileWriter.h"<br> <br> int main (int argc, char* argv[])<br> {<br> =C2=A0 const unsigned int ImageDimension =3D 3;<br> <br> =C2=A0 using VectorType =3D itk::Vector<float, ImageDimension>;<br> =C2=A0 using InputImageType =3D itk::Image<VectorType, ImageDimension>= ;;<br> =C2=A0 using InputPointType =3D InputImageType::PointType;<br> =C2=A0 using OutputImageType =3D itk::Image<VectorType, ImageDimension += 1>;<br> =C2=A0 using PointSetType =3D itk::PointSet<VectorType, ImageDimension += 1>;<br> =C2=A0 using OutputPointType =3D PointSetType::PointType;<br> =C2=A0 using ReaderType =3D itk::ImageFileReader<InputImageType>;<br> =C2=A0 using WriterType =3D itk::ImageFileWriter<OutputImageType>;<br= > <br> =C2=A0 PointSetType::Pointer pointSet =3D PointSetType::New();<br> =C2=A0 unsigned long pointId =3D 0;<br> =C2=A0 InputImageType::PointType origin;<br> =C2=A0 InputImageType::SpacingType spacing;<br> =C2=A0 InputImageType::SizeType size;<br> =C2=A0 for (int i =3D 1; i < argc; ++i)<br> =C2=A0 {<br> =C2=A0 =C2=A0 ReaderType::Pointer reader =3D ReaderType::New();<br> =C2=A0 =C2=A0 std::cout << "Reading file " << argv[i]= << std::endl;<br> =C2=A0 =C2=A0 reader->SetFileName(argv[i]);<br> =C2=A0 =C2=A0 try<br> =C2=A0 =C2=A0 {<br> =C2=A0 =C2=A0 =C2=A0 reader->Update();<br> =C2=A0 =C2=A0 }<br> =C2=A0 =C2=A0 catch (itk::ExceptionObject & err)<br> =C2=A0 =C2=A0 {<br> =C2=A0 =C2=A0 =C2=A0 std::cout << "ExceptionObject caught !"= ; << std::endl;<br> =C2=A0 =C2=A0 =C2=A0 std::cout << err << std::endl;<br> =C2=A0 =C2=A0 =C2=A0 return EXIT_FAILURE;<br> =C2=A0 =C2=A0 }<br> =C2=A0 =C2=A0 using IteratorType =3D itk::ImageRegionConstIterator<Input= ImageType>;<br> =C2=A0 =C2=A0 const InputImageType * image =3D reader->GetOutput();<br> =C2=A0 =C2=A0 IteratorType it(image, image->GetBufferedRegion());<br> =C2=A0 =C2=A0 it.GoToBegin();<br> =C2=A0 =C2=A0 InputPointType inPoint;<br> =C2=A0 =C2=A0 OutputPointType outPoint;<br> =C2=A0 =C2=A0 int iCount =3D 0;<br> =C2=A0 =C2=A0 while (!it.IsAtEnd())<br> =C2=A0 =C2=A0 {<br> =C2=A0 =C2=A0 =C2=A0 inPoint =3D image->GetPixel(it.GetIndex());<br> =C2=A0 =C2=A0 =C2=A0 for (int j =3D 0; j < ImageDimension; ++j)<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 outPoint[j] =3D inPoint[i];<br> =C2=A0 =C2=A0 =C2=A0 outPoint[ImageDimension] =3D i - 1;<br> =C2=A0 =C2=A0 =C2=A0 pointSet->SetPoint(pointId, outPoint);<br> =C2=A0 =C2=A0 =C2=A0 // Transfer the pixel data to the value associated wit= h the point.<br> =C2=A0 =C2=A0 =C2=A0 pointSet->SetPointData(pointId, it.Get());<br> =C2=A0 =C2=A0 =C2=A0 ++it;<br> =C2=A0 =C2=A0 =C2=A0 ++pointId;<br> =C2=A0 =C2=A0 =C2=A0 ++iCount;<br> =C2=A0 =C2=A0 }<br> =C2=A0 =C2=A0 std::cout << "Number of points in " << = argv[i] << " =3D " << iCount <<<br> std::endl;<br> =C2=A0 =C2=A0 std::cout << "Total number of points =3D " &l= t;<<br> pointSet->GetNumberOfPoints() << std::endl;<br> =C2=A0 =C2=A0 origin =3D image->GetOrigin();<br> =C2=A0 =C2=A0 spacing =3D image->GetSpacing();<br> =C2=A0 =C2=A0 size =3D image->GetLargestPossibleRegion().GetSize();<br> =C2=A0 }<br> <br> =C2=A0 typedef itk::BSplineScatteredDataPointSetToImageFilter < PointSet= Type,<br> OutputImageType > SplineFilterType;<br> =C2=A0 SplineFilterType::Pointer splineFilter =3D SplineFilterType::New();<= br> <br> =C2=A0 int splineorder=3D3; // complexity of the spline<br> <br> =C2=A0 SplineFilterType::ArrayType ncontrol;<br> =C2=A0 ncontrol[0]=3Dsplineorder + 1;<br> =C2=A0 SplineFilterType::ArrayType closedim;<br> =C2=A0 closedim[0]=3D 0;<br> <br> =C2=A0 OutputImageType::PointType parametricDomainOrigin;<br> =C2=A0 OutputImageType::SpacingType parametricDomainSpacing;<br> =C2=A0 OutputImageType::SizeType parametricDomainSize;<br> =C2=A0 for (int i =3D 0; i < ImageDimension; ++i)<br> =C2=A0 {<br> =C2=A0 =C2=A0 parametricDomainOrigin[i] =3D origin[i];<br> =C2=A0 =C2=A0 parametricDomainSpacing[i] =3D spacing[i];<br> =C2=A0 =C2=A0 parametricDomainSize[i] =3D size[i];<br> =C2=A0 }<br> =C2=A0 parametricDomainOrigin[ImageDimension] =3D 0;<br> =C2=A0 parametricDomainSize[ImageDimension] =3D argc - 2;<br> =C2=A0 parametricDomainSpacing[ImageDimension] =3D 10.0;<br> <br> =C2=A0 splineFilter->SetGenerateOutputImage( true );=C2=A0 =C2=A0// the = only reason to turn<br> this off is if one only wants to use the control point lattice for further<= br> processing<br> =C2=A0 splineFilter->SetInput ( pointSet );<br> =C2=A0 splineFilter->SetSplineOrder ( splineorder );<br> =C2=A0 splineFilter->SetNumberOfControlPoints ( ncontrol );<br> =C2=A0 splineFilter->SetNumberOfLevels( 3 );<br> =C2=A0 splineFilter->SetCloseDimension ( closedim );<br> =C2=A0 splineFilter->SetSize( parametricDomainSize );<br> =C2=A0 splineFilter->SetSpacing( parametricDomainSpacing );<br> =C2=A0 splineFilter->SetOrigin( parametricDomainOrigin );<br> =C2=A0 std::cout << "Before update spline filter" << = std::endl;<br> =C2=A0 splineFilter->Update();<br> =C2=A0 std::cout << "After update spline filter" << s= td::endl;<br> <br> =C2=A0 WriterType::Pointer writer =3D WriterType::New();<br> =C2=A0 writer->SetInput(splineFilter->GetOutput());<br> =C2=A0 writer->SetFileName("Output.mhd");<br> =C2=A0 writer->Update();<br> =C2=A0 std::cout << "After write image filter" << std= ::endl;<br> <br> =C2=A0 return EXIT_SUCCESS;<br> };<br> <br> <br> <br> --<br> Sent from: <a href=3D"http://itk-users.7.n7.nabble.com/" rel=3D"noreferrer"= target=3D"_blank">http://itk-users.7.n7.nabble.com/</a><br> The ITK community is transitioning from this mailing list to <a href=3D"htt= p://discourse.itk.org" rel=3D"noreferrer" target=3D"_blank">discourse.itk.o= rg</a>. Please join us there!<br> ________________________________<br> Powered by <a href=3D"http://www.kitware.com" rel=3D"noreferrer" target=3D"= _blank">www.kitware.com</a><br> <br> Visit other Kitware open-source projects at<br> <a href=3D"http://www.kitware.com/opensource/opensource.html" rel=3D"norefe= rrer" target=3D"_blank">http://www.kitware.com/opensource/opensource.html</= a><br> <br> Kitware offers ITK Training Courses, for more information visit:<br> <a href=3D"http://www.kitware.com/products/protraining.php" rel=3D"noreferr= er" target=3D"_blank">http://www.kitware.com/products/protraining.php</a><b= r> <br> Please keep messages on-topic and check the ITK FAQ at:<br> <a href=3D"http://www.itk.org/Wiki/ITK_FAQ" rel=3D"noreferrer" target=3D"_b= lank">http://www.itk.org/Wiki/ITK_FAQ</a><br> <br> Follow this link to subscribe/unsubscribe:<br> <a href=3D"https://itk.org/mailman/listinfo/insight-users" rel=3D"noreferre= r" target=3D"_blank">https://itk.org/mailman/listinfo/insight-users</a><br> </blockquote></div> --000000000000b86178057756f5e1-- --===============0995632974== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline The ITK community is transitioning from this mailing list to discourse.itk.org. Please join us there! ________________________________ 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.php Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ Follow this link to subscribe/unsubscribe: https://itk.org/mailman/listinfo/insight-users --===============0995632974==--