[Insight-users] Where is Minimal Path Extraction From FastMarching Speed Function
Dan Mueller
dan.muel at gmail.com
Fri Jul 24 14:48:44 EDT 2009
Hi Hussein,
Please indicate which points (in pixel space) you are trying to
extract paths between.
I tried your datasets with a single start and end point (ie.
SpeedToPathImageFilter) and didn't seem have any issues.
To reproduce my images (on Windows platform):
0. Edit your mhd files with the correct raw files
1. Download SharpImage tool:
http://www.insight-journal.org/browse/publication/161
2. Open C6_1_SpeedFunc_Path.mhd
3. Select "Script" > "Console"
4. Type "SpeedToPath", press enter
5. Select "plus" button to add path points
6. Start point: [5, 1, 6]
7. End point: [10, 24, 16]
8. Select "OK"
9. Type "MaximumProjection", press enter
10. Type "AddPaths", press enter
Cheers, Dan
2009/7/24 asamwm <asamwm at gmail.com>:
> 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 --------------
A non-text attachment was scrubbed...
Name: path_1_1.png
Type: image/png
Size: 6770 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090724/286559bd/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: path_6_1.png
Type: image/png
Size: 5597 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090724/286559bd/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: path_6_1_screen.png
Type: image/png
Size: 23522 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090724/286559bd/attachment-0005.png>
More information about the Insight-users
mailing list