[Insight-users] Where is Minimal Path Extraction From FastMarching Speed Function
asamwm
asamwm at gmail.com
Fri Jul 17 21:21:39 EDT 2009
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/20090717/720b2eca/attachment-0001.htm>
More information about the Insight-users
mailing list