Hi Dan,<br>Your code looks great.<br>I have been trying hard to make it work but it <br>has not got to the finish line yet.<br><br>I am carving out small 3D regions for the speed <br>function (on the range of 30x40x30 pixels).<br>
<br>I map the speed pixels to a float range [0.0 - 1.0];<br>and can verify that the start and end points (1 end point for now)<br>are where they should be in the small image.<br><br>Then I pass this to the filter. I added a function to your code<br>
to read the arrival function at the end :<br>after calling pathFilter->Update()<br>I do (InputImagePointer GetArrivalFunc(){return m_CurrentArrivalFunction;}; )<br>in itkSpeedFunctionToPathFilter.h<br><br>This reveals the following statistics for Arrival:<br>
min=0; max=1.70141e+38, mean=4.26e+36 ...<br>I also read the arrival time at the endpoint: 1431.25<br>(reasonable, right?)<br><br>(Even when I use a fake end point that is very close<br>to the start point, the arrival value is 315.8...)<br>
<br>yet, the path is always empty!<br><br>I tried all 3 suggested optimizers, but the result is the same.<br><br>Can you suggest anything?<br>I am attaching my calling method in case you see anything obvious.<br><br>Thanks for your help.<br>
<br><br>#if defined(_MSC_VER)<br>//Warning about: identifier was truncated to '255' characters in the debug information (MVC6.0 Debug)<br>#pragma warning( disable : 4786 )<br>#endif<br><br>// General includes<br>#include <string><br>
#include <iostream><br><br>// ITK includes<br>#include "itkNumericTraits.h"<br>#include "itkImage.h"<br>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"<br>#include "itkPolyLineParametricPath.h"<br>
#include "itkNearestNeighborInterpolateImageFunction.h"<br>#include "itkLinearInterpolateImageFunction.h"<br>#include "itkArrivalFunctionToPathFilter.h"<br>#include "itkSpeedFunctionToPathFilter.h"<br>
#include "itkPathIterator.h"<br>#include "itkGradientDescentOptimizer.h"<br>#include "itkRegularStepGradientDescentOptimizer.h"<br>#include "itkIterateNeighborhoodOptimizer.h"<br>#include "itkStatisticsImageFilter.h"<br>
<br>#include "SpineRing.h"<br><br>PathType::Pointer SpCandidate::FMPath(SpeedImageType::Pointer speed, PointType start, PointSetContainerType::Pointer endptscont)<br>{<br> //typedef float PixelType;//speed<br>
typedef unsigned char OutputPixelType;<br> //typedef itk::Image< PixelType, spr_SPIMGDIMENSION > ImageType;//speed<br> typedef itk::Image< OutputPixelType, spr_SPIMGDIMENSION > OutputImageType;<br>
//typedef itk::ImageFileReader< ImageType > ReaderType;<br> typedef itk::ImageFileWriter< OutputImageType > WriterType;<br> //typedef itk::SpeedFunctionToPathFilter< SpineImageType, PathType > PathFilterType;<br>
typedef itk::SpeedFunctionToPathFilter< SpeedImageType, PathType > PathFilterType;<br> typedef PathFilterType::CostFunctionType::CoordRepType CoordRepType;<br> typedef itk::PathIterator< OutputImageType, PathType > PathIteratorType;<br>
<br> char* OutputFilename = "TestPath.tif";<br><br> speed->DisconnectPipeline();<br><br> // Create interpolator<br> typedef itk::LinearInterpolateImageFunction<SpeedImageType, CoordRepType><br>
InterpolatorType;<br> InterpolatorType::Pointer interp = InterpolatorType::New();<br><br> // Create cost function<br> PathFilterType::CostFunctionType::Pointer cost =<br> PathFilterType::CostFunctionType::New();<br>
cost->SetInterpolator( interp );<br><br> // Create optimizer<br> ///////////////////////////////////////////////////////////////<br> //typedef itk::GradientDescentOptimizer OptimizerType;<br> //OptimizerType::Pointer optimizer = OptimizerType::New();<br>
//optimizer->SetNumberOfIterations( 1000 );<br> //////////////////////////////////////////////////////////////<br><br> ///////////////////////////////////////////////////////////////////<br> //typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br>
//OptimizerType::Pointer optimizer = OptimizerType::New();<br> //optimizer->SetNumberOfIterations( 1000 );<br> //optimizer->SetMaximumStepLength( 0.5 );<br> //optimizer->SetMinimumStepLength( 0.1 );<br>
//optimizer->SetRelaxationFactor( 0.5 );<br> ///////////////////////////////////////////////////////////////////<br><br> ///////////////////////////////////////////////////////////////////<br> // Create IterateNeighborhoodOptimizer<br>
typedef itk::IterateNeighborhoodOptimizer OptimizerType;<br> OptimizerType::Pointer optimizer = OptimizerType::New();<br> optimizer->MinimizeOn( );<br> optimizer->FullyConnectedOn( );<br> OptimizerType::NeighborhoodSizeType nbrsize( 3 );<br>
//float StepLengthFactor = .2;<br> for (unsigned int i=0; i<3; i++)<br> nbrsize[i] = speed->GetSpacing()[i]; //* StepLengthFactor;<br> optimizer->SetNeighborhoodSize( nbrsize );<br><br> // Create path filter<br>
PathFilterType::Pointer pathFilter = PathFilterType::New();<br> pathFilter->SetInput( speed );<br> pathFilter->SetCostFunction( cost );<br> pathFilter->SetOptimizer( optimizer );<br> pathFilter->SetTerminationValue( 2.0 );<br>
<br> // Setup path points<br> PathFilterType::PointType pstart, pend; <br> pstart = start;<br><br> // Add path information<br> PathFilterType::PathInfo info;<br> info.SetStartPoint( pstart );<br> <br>
pathFilter->AddPathInfo( info );<br> <br><br> PointSetType::PointsContainerIterator pciter = endptscont->Begin();<br> PointType pt;<br> while (pciter != endptscont->End()) <br>
{<br> pt = pciter->Value();<br> //for (int j=0; j<spr_SPIMGDIMENSION; j++) <br> //{<br> // pend = pt;<br> //}<br> // just to test a very close end point!<br> for (int ii=0; ii<3; ii++)<br>
pend[ii]=pstart[ii]-5;<br><br> info.SetEndPoint( pend );<br> //////////////////////////////////<br> // Hussein FIXME<br> // for now<br> break;<br> // need to examine endpt arrivals.. for all endpts. Choose best.<br>
/////////////////////////////////<br> pciter++;<br> }<br><br><br> // Compute the path<br> pathFilter->Update( );<br> //////////////////////////////////////////////////////////<br> /////////////// DEBUG STUFF<br>
SpeedImageType::Pointer Arrival = pathFilter->GetArrivalFunc();<br> typedef itk::ImageFileWriter< SpeedImageType > WriterType2;<br> WriterType2::Pointer writer2 = WriterType2::New();<br> writer2->SetFileName( "Arrival.vtk" );<br>
writer2->SetInput( Arrival );<br> writer2->Update();<br> typedef itk::StatisticsImageFilter<SpeedImageType> StatisticsType;<br> StatisticsType::Pointer statistics = StatisticsType::New();<br> statistics->SetInput(Arrival );<br>
statistics->Update();<br> float imin = statistics->GetMinimum();<br> float imax = statistics->GetMaximum();<br> float imean = statistics->GetMean();<br> float istd = statistics->GetSigma();<br>
<br> std::cout << "Input Image Statistics: min:"<< imin << " max:" << imax << " mean:"<< imean << " std:" << istd << std::endl;<br>
SpeedImageType::IndexType arrIndex;<br> for (int ii=0; ii<3; ii++)<br> arrIndex[ii] = pend[ii];<br> std::cout << "End Arrival Val=" << Arrival->GetPixel(arrIndex) << std::endl;<br>
//itk::ImageRegionIterator<SpeedImageType> arrit(Arrival, Arrival->GetRequestedRegion());<br> //SpeedPixelType spdpxl;<br> //for ( arrit.GoToBegin(); !arrit.IsAtEnd(); ++arrit)<br> //{<br> // std::cout << arrit.Get(); //spregion = (spregion-m)/(M-m); %0 to 1<br>
//}<br> /////////// END DEBUG STUFF<br> //////////////////////////////////////////////////////////////<br> ///////////////<br> // Allocate output image<br> OutputImageType::Pointer output = OutputImageType::New();<br>
output->SetRegions( speed->GetLargestPossibleRegion() );<br> output->SetSpacing( speed->GetSpacing() );<br> output->SetOrigin( speed->GetOrigin() );<br> output->Allocate( );<br> output->FillBuffer( itk::NumericTraits<OutputPixelType>::Zero );<br>
<br> PathType::Pointer path;<br> // Rasterize path<br> for (unsigned int i=0; i<pathFilter->GetNumberOfOutputs(); i++)<br> {<br> // Get the path<br> path = pathFilter->GetOutput( i );<br>
<br> // Check path is valid<br> if ( path->GetVertexList()->Size() == 0 )<br> {<br> std::cout << "WARNING: Path " << (i+1) << " contains no points!" << std::endl;<br>
continue;<br> }<br><br> // Iterate path and convert to image<br> PathIteratorType it( output, path );<br> for (it.GoToBegin(); !it.IsAtEnd(); ++it)<br> {<br> it.Set( itk::NumericTraits<OutputPixelType>::max() );<br>
}<br> }<br><br> // Write output<br> WriterType::Pointer writer = WriterType::New();<br> writer->SetFileName( OutputFilename );<br> writer->SetInput( output );<br> writer->Update();<br>
<br> return path;<br>}<br><br><br><div class="gmail_quote">On Fri, Jul 17, 2009 at 3:57 AM, Dan Mueller <span dir="ltr"><<a href="mailto:dan.muel@gmail.com">dan.muel@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi Hussein,<br>
<div class="im"><br>
> In the ArrivalFunctionToPathFilter superclass, you had AddEndPoint,<br>
> but in SpeedFunctionToPathFilter, it is invalidated with a warning!<br>
> Why is that?<br>
</div><div class="im">> In the current implementation, I would have to rerun the wave as many<br>
> times as there are endpoints and choose the one with least geodesic<br>
> distance.<br>
> Or alternatively, set all those end points as way points and compare their<br>
> geodesics, pick the least, and use it for the path.<br>
<br>
</div>Indeed, SpeedFunctionToPathFilter allows only a single end point. This<br>
is because the filter performs the path extraction in sections (eg.<br>
start point to way point 1, way point 1 to way point 2, way point 2 to<br>
end point). I guess it could be extended in the manner you describe,<br>
but to date has not been. Feel free to make this extension and submit<br>
it to the Insight Journal!<br>
<br>
If your path consists of multiple start points and multiple end points<br>
(ie. no way/intermediate points but multiple candidate start points),<br>
then you should instead use the ArrivalFunctionToPathFilter. This<br>
filter expects as input an arrival function, which you can compute<br>
yourself using FastMarchingImageFilter with all of the start points as<br>
trial points. Then simply add the end points via AddPathEndPoint(..)<br>
and the extracted path will select the closest start point.<br>
<br>
HTH<br>
<br>
Cheers, Dan<br>
<br>
2009/7/16 asamwm <<a href="mailto:asamwm@gmail.com">asamwm@gmail.com</a>>:<br>
<div><div></div><div class="h5">> Hi Dan,<br>
> I am towards the end of connecting things together however, it looks like<br>
> you only allow one end point at a time to be set for the path.<br>
> Please correct me if I am wrong, I am assuming the wave starts at<br>
> startpoint (0 geodesic distance) and ends at endpoint.<br>
><br>
> In the ArrivalFunctionToPathFilter superclass, you had AddEndPoint,<br>
> but in SpeedFunctionToPathFilter, it is invalidated with a warning!<br>
> Why is that?<br>
> The reason I think it is needed is when you have multiple end<br>
> point candidates and you want the wave to stop propagation when it reaches<br>
> the first endpoint (that would be the fastest / closest geodesic pt.).<br>
><br>
> In the current implementation, I would have to rerun the wave as many<br>
> times as there are endpoints and choose the one with least geodesic<br>
> distance.<br>
> Or alternatively, set all those end points as way points and compare their<br>
> geodesics, pick the least, and use it for the path.<br>
><br>
> What do you think?<br>
><br>
> Thanks.<br>
> Hussein<br>
><br>
><br>
> On Fri, Jul 10, 2009 at 12:58 AM, asamwm <<a href="mailto:asamwm@gmail.com">asamwm@gmail.com</a>> wrote:<br>
>><br>
>> Thanks Dan.<br>
>> I'll let you know when I get it in place.<br>
>><br>
>> Regards.<br>
>><br>
>><br>
>> On Fri, Jul 10, 2009 at 12:42 AM, Dan Mueller <<a href="mailto:dan.muel@gmail.com">dan.muel@gmail.com</a>> wrote:<br>
>>><br>
>>> Hi,<br>
>>><br>
>>> The filter is not yet available in the main ITK code archive. You must<br>
>>> download the article and code from here:<br>
>>> <a href="http://www.insight-journal.org/browse/publication/213" target="_blank">http://www.insight-journal.org/browse/publication/213</a><br>
>>><br>
>>> Click on the full download ".zip" link and extract the files somewhere<br>
>>> on your hard drive. As explained in the article (section 2.1, page 3)<br>
>>> the main filter you will want to use is<br>
>>> itk::SpeedFunctionToPathFilter, which is located under the "Source"<br>
>>> folder. To use this filter within your own project, you will need to<br>
>>> copy all of the files under the "Source" folder to your own source<br>
>>> code directory.<br>
>>><br>
>>> Make sure you leave a review at the journal! ;P<br>
>>><br>
>>> Hope this helps.<br>
>>><br>
>>> Cheers, Dan<br>
>>><br>
>>> 2009/7/9 asamwm <<a href="mailto:asamwm@gmail.com">asamwm@gmail.com</a>>:<br>
>>> > Hi all,<br>
>>> > I spent a lot of time searching for this but could not find it.<br>
>>> > There is this insight-journal publication:<br>
>>> > <a href="http://www.insight-journal.org/browse/publication/213" target="_blank">http://www.insight-journal.org/browse/publication/213</a><br>
>>> ><br>
>>> > but where is that filter? I searched under the itk code tree for<br>
>>> > the keywords, speed function, minimal path, .. no success<br>
>>> ><br>
>>> > Any successful examples? any alternatives?<br>
>>> ><br>
>>> > Thanks<br>
>>> ><br>
>>> ><br>
>>> > _____________________________________<br>
>>> > Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
>>> ><br>
>>> > Visit other Kitware open-source projects at<br>
>>> > <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
>>> ><br>
>>> > Please keep messages on-topic and check the ITK FAQ at:<br>
>>> > <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
>>> ><br>
>>> > Follow this link to subscribe/unsubscribe:<br>
>>> > <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
>>> ><br>
>>> ><br>
>><br>
><br>
><br>
</div></div></blockquote></div><br>