Jerome Piovano

Research in Medical Imaging and Computer Vision

• Home
• Research
• Publications
• Teaching
• Biography
• Shared code and data

2D_FastMarching.C

2D Fast Marching

#include <Image.H>
#include <FastMarching.h>
#include <Globals.h>
#include <Schemes.h>

#include <cmath>
#include <sstream>

int main(int argc, char* argv[])
{       
        using namespace Images;
        using namespace levelset;

#if 1 // random seeds/forbidde circles

        const Dimension nx = 500;
        const Dimension ny = 500;
        const Dimension dims[] = {nx, ny};
        Shape<2> shape(dims);

        BaseImage<2,float> image(shape);    // image of size nx*ny
        Eikonal<2,float>   march(&image);   // fast marching eikonal 

        // initialization 
        image = 0.0;

        srand(10000);

        // 20 random seeds
        for (int k = 0 ; k < 20 ; k++){

                int x = static_cast<int>(nx * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
                int y = static_cast<int>(ny * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
                BaseImage<2,float>::Index ind(x, y);

                march.add_trial_point(ind);
        }

        // 20 circles of forbidden points of radius <= 30
        for (int k = 0 ; k < 20 ; k++){

                int x = static_cast<int>(30 + (nx-60) * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
                int y = static_cast<int>(30 + (ny-60) * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
                int r = static_cast<int>(          30 * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
        
                for (int i = x-r ; i < x+r ; i++)
                        for (int j = y-r ; j< y+r ; j++)
                                if(sqrt( (double) Maths::sqr(i - x) +  Maths::sqr( j - y)) < r){
                                        BaseImage<2,float>::Index ind(i, j);
                                        march.add_forbidden_point(ind);
                                }
        }

        timeval _time;
        Globals::initChrono(_time); 

        // fast marching
        march.run();
  
        Globals::printChrono("\n        => FastMarching : ", _time);
        Globals::initChrono(_time); 

        std::ofstream f("2D_FastMarching.inr");
        f << format("inrimage") << image;
        f.close();

        return 0;
#endif

#if 0 // maze + backtrack

        BaseImage<2, float> maze;
        std::ifstream f("image/maze.inr");
        f >> format("inrimage") >> maze;

        BaseImage<2,float> image(maze.shape());    // image
        Eikonal<2,float>   march(&image);          // fast marching eikonal 

        image = 0.0;

        // forbidden points
        for (BaseImage<2, float>::iterator<domain> ind=maze.begin() ; ind!=maze.end() ; ++ind)
                if (maze(ind) > 0) march.add_forbidden_point(ind.position());

        // fast_marching
        BaseImage<2,float>::Index ind(25, 5);
        march.add_trial_point(ind);
        march.run();

        std::ofstream f_fm("2D_FastMarching.inr");
        f_fm<< format("inrimage") << image;
        f_fm.close();

        #if 1
                // back tracking
                BaseImage<2,float> backtrack = maze.copy();
                Globals::normalize<2,float>(backtrack, 0, 200);
        
                BaseImage<2,float>::Index ind2(240, 495);
                BaseImage<2,float>::Index min;
                int inf_coords[2];
                int sup_coords[2];
                while (ind2 != ind){
                        backtrack(ind2) = 255;
                        Schemes<2,float>::find_neighbours(image, ind2, inf_coords, sup_coords);
                        min = ind2;
                        for (int i=0 ; i<2 ; i++){
                                BaseImage<2,float>::Index inf = ind2;
                                BaseImage<2,float>::Index sup = ind2;
                                inf(i+1) = inf_coords[i];
                                sup(i+1) = sup_coords[i];
        
                                if (march(inf) != Eikonal<2,float>::eForbidden && image(inf)<image(min)) min = inf;
                                if (march(sup) != Eikonal<2,float>::eForbidden && image(sup)<image(min)) min = sup;
                        }
                        ind2 = min;
                }
        
                std::ofstream f_backtrack("backtrack.inr"); 
                f_backtrack << format("inrimage") << backtrack;
                f_backtrack.close();
        #endif

        return 0;
#endif
}

For further information, please contact Jerome Piovano - Last update 2008-02-08