Jerome Piovano

Research in Medical Imaging and Computer Vision

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

3D_FastMarching.C

3D Fast Marching

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

#include <cmath>
#include <sstream>

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

#if 1 // sphere

        const Dimension nx = 200;
        const Dimension ny = 200;
        const Dimension nz = 200;
        const Dimension dims[] = {nx, ny, nz};
        Shape<3> shape(dims);

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

        // initialization  
        image = 0.0;
        BaseImage<3,float>::Index ind(100, 100, 100);
        march.add_trial_point(ind);

        march.run();

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

        return 0;

#endif

#if 0 // misc

        const Dimension nx = 200;
        const Dimension ny = 200;
        const Dimension nz = 200; 
        const Dimension dims[] = {nx, ny, nz}; 
        Shape<3> shape(dims);

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

        // initialization  
        image = 0.0;

        srand(10000);

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

                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)));
                int z = static_cast<int>(nz * (double)rand() / ((double)(RAND_MAX)+(double)(1)));
                BaseImage<3,float>::Index ind(x, y, z);

                march.add_trial_point(ind);
        }

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

                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 z = static_cast<int>(30 + (nz-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++)
                                for (int k = z-r ; k< z+r ; k++)
                                        if(sqrt( (double) Maths::sqr(i - x) +  Maths::sqr( j - y) +  Maths::sqr( k - z)) < r){
                                                BaseImage<3,float>::Index ind(i, j, k);
                                                march.add_forbidden_point(ind);
                                        }
        }

        march.run();

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

        return 0;

#endif

}

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