Hi Dan,<br>Thanks for your offer to help.<br><br>I got past my previous point where nothing was happening.<br>Now I realized it is probably a combination between my<br>speed function and the FastMarchingImageFilter implementation.<br>
<br>Background: I had this application in matlab working very well.<br>The FastMarching was actually a toolbox contributed by Gabriel Peyre'<br>whose engine is in C (perform_front_propagation_3d)<br>That works with a speed image made from the original subimage<br>
rescaled to [0,1] with 1 being fastest.<br>(no gradient pre-computed).<br><br>If I use the same with itk FM, it won't work. The arrival values at the<br>start point+1 would be infinity.<br><br>If I take my subimage (which is a crop out of the original)<br>
and use the magnitudegradient filter, itk FM works paritially;i.e.,<br>some cases result in a reasonable arrival and path, others do not.<br><br>In all cases so far, the matlab implementation seems to result in better<br>
paths. <br><br>If you can take a look, please visit <br><a href="http://www.rpi.edu/~sharah/testimages/">http://www.rpi.edu/~sharah/testimages/</a><br>You will find 1 working example (*1_1_*)<br>and another failing (*6_1_*)<br>
<br>each has a subregion image, <br>
with an rgb max projection with the start point in red <br>and some endpoints in green (currently I use<br>only one endpoint for now).<br><br>The speed function is the one generated by<br>GradientMagnitudeRecursiveGaussianImageFilter<br>
<div id=":15u" class="ii gt">
with sigma=2. (also the results are sensitive<br>to this sigma because the flow region is very thin.)<br><br>The arrival image and paths are extracted using your <br>filter.<br><br>I would like to know your perspective on this and whether you <br>
can suggest a better speed function generator given the tiny size<br>
of this application.<br><br>I will be testing other gradient options and perhaps will port<br>Peyre's implementation over to see how that arrival function works<br>with your path extraction pipeline.<br><br>Thanks.<br>
Hussein</div><br><br><div class="gmail_quote"><br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div><div class="h5"><br><br><br><br><div class="gmail_quote">
On Mon, Jul 20, 2009 at 8:53 AM, Dan Mueller <span dir="ltr"><<a href="mailto:dan.muel@gmail.com" target="_blank">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>
<br>
I really need the image data as well. Please upload these somewhere<br>
and send the link via email.<br>
<br>
Regards, Dan<br>
<br>
2009/7/18 asamwm <<a href="mailto:asamwm@gmail.com" target="_blank">asamwm@gmail.com</a>>:<br>
<div><div></div><div>> 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>
> )<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<br>
> 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,<br>
> PointType start, PointSetContainerType::Pointer endptscont)<br>
> {<br>
> //typedef float PixelType;//speed<br>
> typedef unsigned char OutputPixelType;<br>
> //typedef itk::Image< PixelType, spr_SPIMGDIMENSION ><br>
> ImageType;//speed<br>
> typedef itk::Image< OutputPixelType, spr_SPIMGDIMENSION ><br>
> OutputImageType;<br>
> //typedef itk::ImageFileReader< ImageType > ReaderType;<br>
> typedef itk::ImageFileWriter< OutputImageType > WriterType;<br>
> //typedef itk::SpeedFunctionToPathFilter< SpineImageType,<br>
> PathType > PathFilterType;<br>
> typedef itk::SpeedFunctionToPathFilter< SpeedImageType,<br>
> PathType > PathFilterType;<br>
> typedef PathFilterType::CostFunctionType::CoordRepType<br>
> CoordRepType;<br>
> typedef itk::PathIterator< OutputImageType, PathType ><br>
> PathIteratorType;<br>
><br>
> char* OutputFilename = "TestPath.tif";<br>
><br>
> speed->DisconnectPipeline();<br>
><br>
> // Create interpolator<br>
> typedef itk::LinearInterpolateImageFunction<SpeedImageType,<br>
> 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 <<<br>
> " 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) <<<br>
> std::endl;<br>
> //itk::ImageRegionIterator<SpeedImageType> arrit(Arrival,<br>
> 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<br>
> 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<br>
> 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>
> On Fri, Jul 17, 2009 at 3:57 AM, Dan Mueller <<a href="mailto:dan.muel@gmail.com" target="_blank">dan.muel@gmail.com</a>> wrote:<br>
>><br>
>> Hi Hussein,<br>
>><br>
>> > In the ArrivalFunctionToPathFilter superclass, you had AddEndPoint,<br>
>> > but in SpeedFunctionToPathFilter, it is invalidated with a warning!<br>
>> > Why is that?<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<br>
>> > their<br>
>> > geodesics, pick the least, and use it for the path.<br>
>><br>
>> 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" target="_blank">asamwm@gmail.com</a>>:<br>
>> > Hi Dan,<br>
>> > I am towards the end of connecting things together however, it looks<br>
>> > 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<br>
>> > 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<br>
>> > 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" target="_blank">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" target="_blank">dan.muel@gmail.com</a>><br>
>> >> 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" target="_blank">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>
><br>
><br>
</div></div></blockquote></div><br>
</div></div></blockquote></div><br>