[Insight-users] Where is Minimal Path Extraction From FastMarching Speed Function
asamwm
asamwm at gmail.com
Fri Jul 24 12:55:57 EDT 2009
Hi Dan,
Thanks for your offer to help.
I got past my previous point where nothing was happening.
Now I realized it is probably a combination between my
speed function and the FastMarchingImageFilter implementation.
Background: I had this application in matlab working very well.
The FastMarching was actually a toolbox contributed by Gabriel Peyre'
whose engine is in C (perform_front_propagation_3d)
That works with a speed image made from the original subimage
rescaled to [0,1] with 1 being fastest.
(no gradient pre-computed).
If I use the same with itk FM, it won't work. The arrival values at the
start point+1 would be infinity.
If I take my subimage (which is a crop out of the original)
and use the magnitudegradient filter, itk FM works paritially;i.e.,
some cases result in a reasonable arrival and path, others do not.
In all cases so far, the matlab implementation seems to result in better
paths.
If you can take a look, please visit
http://www.rpi.edu/~sharah/testimages/
You will find 1 working example (*1_1_*)
and another failing (*6_1_*)
each has a subregion image,
with an rgb max projection with the start point in red
and some endpoints in green (currently I use
only one endpoint for now).
The speed function is the one generated by
GradientMagnitudeRecursiveGaussianImageFilter
with sigma=2. (also the results are sensitive
to this sigma because the flow region is very thin.)
The arrival image and paths are extracted using your
filter.
I would like to know your perspective on this and whether you
can suggest a better speed function generator given the tiny size
of this application.
I will be testing other gradient options and perhaps will port
Peyre's implementation over to see how that arrival function works
with your path extraction pipeline.
Thanks.
Hussein
>
>
>
> On Mon, Jul 20, 2009 at 8:53 AM, Dan Mueller <dan.muel at gmail.com> wrote:
>
>> Hi Hussein,
>>
>> I really need the image data as well. Please upload these somewhere
>> and send the link via email.
>>
>> Regards, Dan
>>
>> 2009/7/18 asamwm <asamwm at gmail.com>:
>> > Hi Dan,
>> > Your code looks great.
>> > I have been trying hard to make it work but it
>> > has not got to the finish line yet.
>> >
>> > I am carving out small 3D regions for the speed
>> > function (on the range of 30x40x30 pixels).
>> >
>> > I map the speed pixels to a float range [0.0 - 1.0];
>> > and can verify that the start and end points (1 end point for now)
>> > are where they should be in the small image.
>> >
>> > Then I pass this to the filter. I added a function to your code
>> > to read the arrival function at the end :
>> > after calling pathFilter->Update()
>> > I do (InputImagePointer GetArrivalFunc(){return
>> m_CurrentArrivalFunction;};
>> > )
>> > in itkSpeedFunctionToPathFilter.h
>> >
>> > This reveals the following statistics for Arrival:
>> > min=0; max=1.70141e+38, mean=4.26e+36 ...
>> > I also read the arrival time at the endpoint: 1431.25
>> > (reasonable, right?)
>> >
>> > (Even when I use a fake end point that is very close
>> > to the start point, the arrival value is 315.8...)
>> >
>> > yet, the path is always empty!
>> >
>> > I tried all 3 suggested optimizers, but the result is the same.
>> >
>> > Can you suggest anything?
>> > I am attaching my calling method in case you see anything obvious.
>> >
>> > Thanks for your help.
>> >
>> >
>> > #if defined(_MSC_VER)
>> > //Warning about: identifier was truncated to '255' characters in the
>> debug
>> > information (MVC6.0 Debug)
>> > #pragma warning( disable : 4786 )
>> > #endif
>> >
>> > // General includes
>> > #include <string>
>> > #include <iostream>
>> >
>> > // ITK includes
>> > #include "itkNumericTraits.h"
>> > #include "itkImage.h"
>> > #include "itkImageFileReader.h"
>> > #include "itkImageFileWriter.h"
>> > #include "itkPolyLineParametricPath.h"
>> > #include "itkNearestNeighborInterpolateImageFunction.h"
>> > #include "itkLinearInterpolateImageFunction.h"
>> > #include "itkArrivalFunctionToPathFilter.h"
>> > #include "itkSpeedFunctionToPathFilter.h"
>> > #include "itkPathIterator.h"
>> > #include "itkGradientDescentOptimizer.h"
>> > #include "itkRegularStepGradientDescentOptimizer.h"
>> > #include "itkIterateNeighborhoodOptimizer.h"
>> > #include "itkStatisticsImageFilter.h"
>> >
>> > #include "SpineRing.h"
>> >
>> > PathType::Pointer SpCandidate::FMPath(SpeedImageType::Pointer speed,
>> > PointType start, PointSetContainerType::Pointer endptscont)
>> > {
>> > //typedef float PixelType;//speed
>> > typedef unsigned char OutputPixelType;
>> > //typedef itk::Image< PixelType, spr_SPIMGDIMENSION >
>> > ImageType;//speed
>> > typedef itk::Image< OutputPixelType, spr_SPIMGDIMENSION >
>> > OutputImageType;
>> > //typedef itk::ImageFileReader< ImageType > ReaderType;
>> > typedef itk::ImageFileWriter< OutputImageType > WriterType;
>> > //typedef itk::SpeedFunctionToPathFilter< SpineImageType,
>> > PathType > PathFilterType;
>> > typedef itk::SpeedFunctionToPathFilter< SpeedImageType,
>> > PathType > PathFilterType;
>> > typedef PathFilterType::CostFunctionType::CoordRepType
>> > CoordRepType;
>> > typedef itk::PathIterator< OutputImageType, PathType >
>> > PathIteratorType;
>> >
>> > char* OutputFilename = "TestPath.tif";
>> >
>> > speed->DisconnectPipeline();
>> >
>> > // Create interpolator
>> > typedef itk::LinearInterpolateImageFunction<SpeedImageType,
>> > CoordRepType>
>> > InterpolatorType;
>> > InterpolatorType::Pointer interp = InterpolatorType::New();
>> >
>> > // Create cost function
>> > PathFilterType::CostFunctionType::Pointer cost =
>> > PathFilterType::CostFunctionType::New();
>> > cost->SetInterpolator( interp );
>> >
>> > // Create optimizer
>> > ///////////////////////////////////////////////////////////////
>> > //typedef itk::GradientDescentOptimizer OptimizerType;
>> > //OptimizerType::Pointer optimizer = OptimizerType::New();
>> > //optimizer->SetNumberOfIterations( 1000 );
>> > //////////////////////////////////////////////////////////////
>> >
>> > ///////////////////////////////////////////////////////////////////
>> > //typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>> > //OptimizerType::Pointer optimizer = OptimizerType::New();
>> > //optimizer->SetNumberOfIterations( 1000 );
>> > //optimizer->SetMaximumStepLength( 0.5 );
>> > //optimizer->SetMinimumStepLength( 0.1 );
>> > //optimizer->SetRelaxationFactor( 0.5 );
>> > ///////////////////////////////////////////////////////////////////
>> >
>> > ///////////////////////////////////////////////////////////////////
>> > // Create IterateNeighborhoodOptimizer
>> > typedef itk::IterateNeighborhoodOptimizer OptimizerType;
>> > OptimizerType::Pointer optimizer = OptimizerType::New();
>> > optimizer->MinimizeOn( );
>> > optimizer->FullyConnectedOn( );
>> > OptimizerType::NeighborhoodSizeType nbrsize( 3 );
>> > //float StepLengthFactor = .2;
>> > for (unsigned int i=0; i<3; i++)
>> > nbrsize[i] = speed->GetSpacing()[i]; //* StepLengthFactor;
>> > optimizer->SetNeighborhoodSize( nbrsize );
>> >
>> > // Create path filter
>> > PathFilterType::Pointer pathFilter = PathFilterType::New();
>> > pathFilter->SetInput( speed );
>> > pathFilter->SetCostFunction( cost );
>> > pathFilter->SetOptimizer( optimizer );
>> > pathFilter->SetTerminationValue( 2.0 );
>> >
>> > // Setup path points
>> > PathFilterType::PointType pstart, pend;
>> > pstart = start;
>> >
>> > // Add path information
>> > PathFilterType::PathInfo info;
>> > info.SetStartPoint( pstart );
>> >
>> > pathFilter->AddPathInfo( info );
>> >
>> >
>> > PointSetType::PointsContainerIterator pciter =
>> endptscont->Begin();
>> > PointType pt;
>> > while (pciter != endptscont->End())
>> > {
>> > pt = pciter->Value();
>> > //for (int j=0; j<spr_SPIMGDIMENSION; j++)
>> > //{
>> > // pend = pt;
>> > //}
>> > // just to test a very close end point!
>> > for (int ii=0; ii<3; ii++)
>> > pend[ii]=pstart[ii]-5;
>> >
>> > info.SetEndPoint( pend );
>> > //////////////////////////////////
>> > // Hussein FIXME
>> > // for now
>> > break;
>> > // need to examine endpt arrivals.. for all endpts. Choose
>> best.
>> > /////////////////////////////////
>> > pciter++;
>> > }
>> >
>> >
>> > // Compute the path
>> > pathFilter->Update( );
>> > //////////////////////////////////////////////////////////
>> > /////////////// DEBUG STUFF
>> > SpeedImageType::Pointer Arrival = pathFilter->GetArrivalFunc();
>> > typedef itk::ImageFileWriter< SpeedImageType > WriterType2;
>> > WriterType2::Pointer writer2 = WriterType2::New();
>> > writer2->SetFileName( "Arrival.vtk" );
>> > writer2->SetInput( Arrival );
>> > writer2->Update();
>> > typedef itk::StatisticsImageFilter<SpeedImageType> StatisticsType;
>> > StatisticsType::Pointer statistics = StatisticsType::New();
>> > statistics->SetInput(Arrival );
>> > statistics->Update();
>> > float imin = statistics->GetMinimum();
>> > float imax = statistics->GetMaximum();
>> > float imean = statistics->GetMean();
>> > float istd = statistics->GetSigma();
>> >
>> > std::cout << "Input Image Statistics: min:"<< imin << " max:" <<
>> imax <<
>> > " mean:"<< imean << " std:" << istd << std::endl;
>> > SpeedImageType::IndexType arrIndex;
>> > for (int ii=0; ii<3; ii++)
>> > arrIndex[ii] = pend[ii];
>> > std::cout << "End Arrival Val=" << Arrival->GetPixel(arrIndex) <<
>> > std::endl;
>> > //itk::ImageRegionIterator<SpeedImageType> arrit(Arrival,
>> > Arrival->GetRequestedRegion());
>> > //SpeedPixelType spdpxl;
>> > //for ( arrit.GoToBegin(); !arrit.IsAtEnd(); ++arrit)
>> > //{
>> > // std::cout << arrit.Get(); //spregion = (spregion-m)/(M-m); %0
>> to
>> > 1
>> > //}
>> > /////////// END DEBUG STUFF
>> > //////////////////////////////////////////////////////////////
>> > ///////////////
>> > // Allocate output image
>> > OutputImageType::Pointer output = OutputImageType::New();
>> > output->SetRegions( speed->GetLargestPossibleRegion() );
>> > output->SetSpacing( speed->GetSpacing() );
>> > output->SetOrigin( speed->GetOrigin() );
>> > output->Allocate( );
>> > output->FillBuffer( itk::NumericTraits<OutputPixelType>::Zero );
>> >
>> > PathType::Pointer path;
>> > // Rasterize path
>> > for (unsigned int i=0; i<pathFilter->GetNumberOfOutputs(); i++)
>> > {
>> > // Get the path
>> > path = pathFilter->GetOutput( i );
>> >
>> > // Check path is valid
>> > if ( path->GetVertexList()->Size() == 0 )
>> > {
>> > std::cout << "WARNING: Path " << (i+1) << " contains
>> no
>> > points!" << std::endl;
>> > continue;
>> > }
>> >
>> > // Iterate path and convert to image
>> > PathIteratorType it( output, path );
>> > for (it.GoToBegin(); !it.IsAtEnd(); ++it)
>> > {
>> > it.Set( itk::NumericTraits<OutputPixelType>::max() );
>> > }
>> > }
>> >
>> > // Write output
>> > WriterType::Pointer writer = WriterType::New();
>> > writer->SetFileName( OutputFilename );
>> > writer->SetInput( output );
>> > writer->Update();
>> >
>> > return path;
>> > }
>> >
>> >
>> > On Fri, Jul 17, 2009 at 3:57 AM, Dan Mueller <dan.muel at gmail.com>
>> wrote:
>> >>
>> >> Hi Hussein,
>> >>
>> >> > In the ArrivalFunctionToPathFilter superclass, you had AddEndPoint,
>> >> > but in SpeedFunctionToPathFilter, it is invalidated with a warning!
>> >> > Why is that?
>> >> > In the current implementation, I would have to rerun the wave as many
>> >> > times as there are endpoints and choose the one with least geodesic
>> >> > distance.
>> >> > Or alternatively, set all those end points as way points and compare
>> >> > their
>> >> > geodesics, pick the least, and use it for the path.
>> >>
>> >> Indeed, SpeedFunctionToPathFilter allows only a single end point. This
>> >> is because the filter performs the path extraction in sections (eg.
>> >> start point to way point 1, way point 1 to way point 2, way point 2 to
>> >> end point). I guess it could be extended in the manner you describe,
>> >> but to date has not been. Feel free to make this extension and submit
>> >> it to the Insight Journal!
>> >>
>> >> If your path consists of multiple start points and multiple end points
>> >> (ie. no way/intermediate points but multiple candidate start points),
>> >> then you should instead use the ArrivalFunctionToPathFilter. This
>> >> filter expects as input an arrival function, which you can compute
>> >> yourself using FastMarchingImageFilter with all of the start points as
>> >> trial points. Then simply add the end points via AddPathEndPoint(..)
>> >> and the extracted path will select the closest start point.
>> >>
>> >> HTH
>> >>
>> >> Cheers, Dan
>> >>
>> >> 2009/7/16 asamwm <asamwm at gmail.com>:
>> >> > Hi Dan,
>> >> > I am towards the end of connecting things together however, it looks
>> >> > like
>> >> > you only allow one end point at a time to be set for the path.
>> >> > Please correct me if I am wrong, I am assuming the wave starts at
>> >> > startpoint (0 geodesic distance) and ends at endpoint.
>> >> >
>> >> > In the ArrivalFunctionToPathFilter superclass, you had AddEndPoint,
>> >> > but in SpeedFunctionToPathFilter, it is invalidated with a warning!
>> >> > Why is that?
>> >> > The reason I think it is needed is when you have multiple end
>> >> > point candidates and you want the wave to stop propagation when it
>> >> > reaches
>> >> > the first endpoint (that would be the fastest / closest geodesic
>> pt.).
>> >> >
>> >> > In the current implementation, I would have to rerun the wave as many
>> >> > times as there are endpoints and choose the one with least geodesic
>> >> > distance.
>> >> > Or alternatively, set all those end points as way points and compare
>> >> > their
>> >> > geodesics, pick the least, and use it for the path.
>> >> >
>> >> > What do you think?
>> >> >
>> >> > Thanks.
>> >> > Hussein
>> >> >
>> >> >
>> >> > On Fri, Jul 10, 2009 at 12:58 AM, asamwm <asamwm at gmail.com> wrote:
>> >> >>
>> >> >> Thanks Dan.
>> >> >> I'll let you know when I get it in place.
>> >> >>
>> >> >> Regards.
>> >> >>
>> >> >>
>> >> >> On Fri, Jul 10, 2009 at 12:42 AM, Dan Mueller <dan.muel at gmail.com>
>> >> >> wrote:
>> >> >>>
>> >> >>> Hi,
>> >> >>>
>> >> >>> The filter is not yet available in the main ITK code archive. You
>> must
>> >> >>> download the article and code from here:
>> >> >>> http://www.insight-journal.org/browse/publication/213
>> >> >>>
>> >> >>> Click on the full download ".zip" link and extract the files
>> somewhere
>> >> >>> on your hard drive. As explained in the article (section 2.1, page
>> 3)
>> >> >>> the main filter you will want to use is
>> >> >>> itk::SpeedFunctionToPathFilter, which is located under the "Source"
>> >> >>> folder. To use this filter within your own project, you will need
>> to
>> >> >>> copy all of the files under the "Source" folder to your own source
>> >> >>> code directory.
>> >> >>>
>> >> >>> Make sure you leave a review at the journal! ;P
>> >> >>>
>> >> >>> Hope this helps.
>> >> >>>
>> >> >>> Cheers, Dan
>> >> >>>
>> >> >>> 2009/7/9 asamwm <asamwm at gmail.com>:
>> >> >>> > Hi all,
>> >> >>> > I spent a lot of time searching for this but could not find it.
>> >> >>> > There is this insight-journal publication:
>> >> >>> > http://www.insight-journal.org/browse/publication/213
>> >> >>> >
>> >> >>> > but where is that filter? I searched under the itk code tree for
>> >> >>> > the keywords, speed function, minimal path, .. no success
>> >> >>> >
>> >> >>> > Any successful examples? any alternatives?
>> >> >>> >
>> >> >>> > Thanks
>> >> >>> >
>> >> >>> >
>> >> >>> > _____________________________________
>> >> >>> > Powered by www.kitware.com
>> >> >>> >
>> >> >>> > Visit other Kitware open-source projects at
>> >> >>> > http://www.kitware.com/opensource/opensource.html
>> >> >>> >
>> >> >>> > Please keep messages on-topic and check the ITK FAQ at:
>> >> >>> > http://www.itk.org/Wiki/ITK_FAQ
>> >> >>> >
>> >> >>> > Follow this link to subscribe/unsubscribe:
>> >> >>> > http://www.itk.org/mailman/listinfo/insight-users
>> >> >>> >
>> >> >>> >
>> >> >>
>> >> >
>> >> >
>> >
>> >
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090724/4d943f29/attachment-0001.htm>
More information about the Insight-users
mailing list