[Insight-users] Where is Minimal Path Extraction From FastMarching Speed Function
Dan Mueller
dan.muel at gmail.com
Mon Jul 20 08:53:32 EDT 2009
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
>> >>> >
>> >>> >
>> >>
>> >
>> >
>
>
More information about the Insight-users
mailing list