mcm/FDS_MCM.c

Go to the documentation of this file.

/*
 * Copyright (c) 2010 Adina Ciomaga <ciomaga@cmla.ens-cachan.fr>
 *                    Marco Mondelli <m.mondelli@sssup.it>
 * All rights reserved.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Pre processing
 */

#include <stdlib.h>
#include <tiffio.h>
#include "io_tiff_routine.h"
#include <time.h>
#include <math.h>

/*
 * Print error messages
 */


#define FATAL(MSG)\
        do {\
          fprintf(stderr, MSG "\n");\
        abort();\
        } while (0);

#define INFO(MSG)\
    do {\
        fprintf(stderr, MSG "\n");\
    } while (0);

#define DEBUG(MSG)\
    do {\
        if (debug_flag)\
            fprintf(stderr, __FILE__ ":%03i " MSG "\n", __LINE__);\
    } while (0);


/*
 * Compute one iteration of Heat Equation
 */

static void heat (float *ptr_out, float *ptr_in, size_t n_col, int i, 
                  int i_prev, int i_next, int j, int j_prev, int j_next, 
                  float dt)
{
    float laplacian;
    laplacian = -4.0 * (*(ptr_in + i*n_col + j)) +
                *(ptr_in + i_prev*n_col + j) + *(ptr_in + i*n_col + j_prev) +
                *(ptr_in + i*n_col + j_next) + *(ptr_in + i_next*n_col + j);

    *(ptr_out+i*n_col + j)= *(ptr_in+i*n_col+j) + 0.5 * dt * laplacian;
}

/*
 * Compute one iteration of the Mean Curvature Motion
 */

static void mcm(float *ptr_out, float *ptr_in, size_t n_row, size_t n_col,
                float dt, float t_g)
{
    unsigned int  i, j, i_next, i_prev, j_next, j_prev;
    float         grad, s_x, s_y, lambda0, lambda1, lambda2, lambda3, lambda4;

    /* iterate on j and i, following the array order */
    for (i=0; i<n_row; i++ )  {

        /* symmetry for borders */
        i_next = (i<n_row-1?i+1:i);
        i_prev = (i>0?i-1:i);

        for (j=0; j<n_col; j++)  {

            /* symmetry for borders */
            j_next = (j<n_col-1?j+1:j);
            j_prev = (j>0?j-1:j);

            /* computation of s_x, s_y and gradient */
            s_x= 2 * ( *(ptr_in+i*n_col+j_next) - *(ptr_in+i*n_col+j_prev) ) +
                *(ptr_in+i_prev*n_col+j_next) - *(ptr_in+i_prev*n_col+j_prev) +
                 *(ptr_in+i_next*n_col+j_next) - *(ptr_in+i_next*n_col+j_prev);
            s_y= 2 * ( *(ptr_in+i_next*n_col+j) - *(ptr_in+i_prev*n_col+j) ) +
                *(ptr_in+i_next*n_col+j_next) - *(ptr_in+i_prev*n_col+j_next) +
                *(ptr_in+i_next*n_col+j_prev) - *(ptr_in+i_prev*n_col+j_prev);
            grad = 0.125 * sqrt( (s_x*s_x) + (s_y*s_y) );

            if (grad>t_g) {
                /* MCM diffusion */
                grad*=8;   /* i.e. grad=\sqrt{s_x^2 + s_y^2} */
                lambda0 = 0.5 - ( (s_x*s_x*s_y*s_y) / ( pow(grad,4) ) );
                lambda1= 2 * lambda0 - ( (s_x*s_x) / (grad*grad) );
                lambda2= 2 * lambda0 - ( (s_y*s_y) / (grad*grad) );
                lambda3= -lambda0 + 0.5 * ( 1 - ( (s_x*s_y)/(grad*grad) ) );
                lambda4= -lambda0 + 0.5 * ( 1 + ( (s_x*s_y)/(grad*grad) ) );

                *(ptr_out+i*n_col + j)= *(ptr_in+i*n_col+j) + dt * (
                                         -4 * lambda0 * (*(ptr_in+i*n_col+j)) +
                                         lambda1 * (*(ptr_in+i*n_col+j_next)  +
                                                    *(ptr_in+i*n_col+j_prev)) +
                                         lambda2 * (*(ptr_in+i_next*n_col+j)  +
                                                    *(ptr_in+i_prev*n_col+j)) +
                                         lambda3 * (*(ptr_in+i_prev*n_col
                                                                     +j_prev) +
                                                    *(ptr_in+i_next*n_col
                                                                     +j_next))+
                                         lambda4 * (*(ptr_in+i_prev*n_col
                                                                     +j_next) +
                                                    *(ptr_in+i_next*n_col
                                                                    +j_prev)));
            }

            else
                /* Heat Equation */
                heat (ptr_out, ptr_in, n_col, i, i_prev, i_next, 
                      j, j_prev, j_next, dt);
        }
    }
}


/*
 * Main function section
 */

int main(int argc, char *argv[])

{
    float        *data_in;                    /* input data */
    float           *ptr_in_rgb[3];
    float        *data_out;                   /* output data */
    float           *ptr_out_rgb[3];
    size_t       n_col, n_row;                /* number of columns and rows */
    size_t       n_channels;                  /* number of channels */
    float        *ptr_in, *ptr_out, *ptr_end;
    int          n_iter=0;                    /* number of iterations */
    float        R;                           /* normalization scale */
    float        t_g=4;                       /* gradient threshold */
    float        dt=0.1;                      /* time step */
    float        n_iter_f;                    /* final number of iterations */
    unsigned int k;
    int          m;
    clock_t      t0, t1;

    /*
     * argv[1] contains the normalized scale
     * argv[2] contains the input file name
     * argv[3] contains the output file name
     */
    if (4!=argc)
        FATAL("Error in the number of arguments!\n");
    if (0==sscanf(argv[1],"%f", &R))
        FATAL("The third argument must be a float!\n");

    /* R<0 is an error because R must be a positive real number */
    if (0 > R) {
        printf("The normalized scale must be a positive real number!\n");
        return 0;
    }

    /* start clock */
    t0 = clock();

    /* number of iterations needed */
    n_iter_f= (R * R) / (2 * dt);

    /* round n_iter_f (must be integer) */
    n_iter= (int) (n_iter_f+0.5);

    /* read the TIFF image */
    if (NULL == (data_in = read_tiff_f32(argv[2], 
                                         &n_col, &n_row, &n_channels)))
        FATAL("error while reading the TIFF image");

    ptr_in_rgb[0] = data_in;
    ptr_in_rgb[1] = data_in + n_row * n_col;
    ptr_in_rgb[2] = data_in + 2 * n_row * n_col;

    /* look for identical RGB channels (check if the image is BW or true RGB) */
    if (1!=n_channels) {
        n_channels = 1;
        for (k=0; k<(n_row*n_col-1); k++)
            if ((*(ptr_in_rgb[0]+k)!= *(ptr_in_rgb[1]+k))||
                (*(ptr_in_rgb[0]+k)!=*(ptr_in_rgb[2]+k))) {
                n_channels=3;
                break;
            }
    }

    /* error if n_iter<0 */
    if (0>n_iter)
        printf("The number of iteration must be a positive number. \n");
    /* input=output if n_iter=0 */
    else if (0==n_iter) {
        write_tiff_f32(argv[3], data_in, n_col, n_row, n_channels);
    } else {
        /* allocate memory for data_out */
        if (NULL == (data_out = (float *) 
                            malloc(n_channels* n_col * n_row * sizeof(float))))
            FATAL("allocation error");

        /* apply MCM to each RGB channel */
        ptr_out_rgb[0] = data_out;
        ptr_out_rgb[1] = data_out +  n_row * n_col;
        ptr_out_rgb[2] = data_out + 2 * n_row * n_col;
        for (k=0; k<n_channels; k++) {
            /* initialize pointers */
            ptr_in  = ptr_in_rgb[k];
            ptr_out = ptr_out_rgb[k];
            ptr_end = ptr_in + n_row*n_col;
            for (m=0; m<n_iter; m++)   {
                /* apply mean curvature motion */
                mcm (ptr_out, ptr_in, n_row, n_col, dt, t_g);
                /* update arrays */
                while (ptr_in<ptr_end) {
                    *ptr_in = *ptr_out;
                    ptr_in++;
                    ptr_out++;
                }
                ptr_in=ptr_in_rgb[k];
                ptr_out=ptr_out_rgb[k];
            }
            /* saturate values to [0,255] */
            while (ptr_in<ptr_end) {
                if (*ptr_in>255) *ptr_in=255;
                if (*ptr_in<0) *ptr_in=0;
                ptr_in++;
            }
        }
        /* write the output image */
        write_tiff_f32(argv[3], data_in, n_col, n_row, n_channels);
        /* free memory */
        free(data_out);
    }

    /* stop clock */
    t1 = clock();
    /* write in the stdout the time required by the algorithm */
    printf("CPU time consumed = %f seconds\n", 
           (double)(t1-t0)/(double)CLOCKS_PER_SEC);
    /* free memory */
    free(data_in);
    return 0;
}

---