amss/FDS_AMSS.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 <time.h>
#include <math.h>
#include "io_tiff_routine.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) + dt * laplacian;

}


/*
 * Compute one iteration of the Affine Morphological Scale Space
 */

static void amss(float *ptr_out, float *ptr_in, size_t n_row, size_t n_col,
                 float dt, float t_g)
{
    float         grad, s_x, s_y, eta0, eta1, eta2, eta3, eta4, l_comb;
    unsigned int  i, j, i_next, i_prev, j_next, j_prev;

    /* 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 of the 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)
                /* Heat Equation */
                heat (ptr_out, ptr_in, n_col, i, i_prev, i_next, 
                      j, j_prev, j_next, dt);
            else {
                /* AMSS diffusion */
                eta0 =0.5 * (s_x*s_x + s_y*s_y)- 
                            (s_x*s_x*s_y*s_y)/(s_x*s_x + s_y*s_y);
                eta1= 2 * eta0 - (s_x*s_x);
                eta2= 2 * eta0 - (s_y*s_y);
                eta3= - eta0 + 0.5 * (s_x*s_x + s_y*s_y - s_x*s_y);
                eta4= - eta0 + 0.5 * (s_x*s_x + s_y*s_y + s_x*s_y);

                l_comb= ( -4 * eta0 * ( *(ptr_in+i*n_col+j) ) +
                          eta1*( *(ptr_in+i*n_col+j_next) + 
                                 *(ptr_in+i*n_col+j_prev) ) +
                          eta2*( *(ptr_in+i_next*n_col+j) + 
                                 *(ptr_in+i_prev*n_col+j) ) +
                          eta3*( *(ptr_in+i_next*n_col+j_next) + 
                                 *(ptr_in+i_prev*n_col+j_prev) ) +
                          eta4*( *(ptr_in+i_prev*n_col+j_next) + 
                                 *(ptr_in+i_next*n_col+j_prev) ) );

                /* pow(x,y) which returns x^y requires that x>=0 */
                if (l_comb>0)
                    *(ptr_out+i*n_col+j)=
                        *(ptr_in+i*n_col+j)+0.25* dt * pow(l_comb,0.33333333);
                else
                    *(ptr_out+i*n_col+j)=
                        *(ptr_in+i*n_col+j)-0.25* dt * pow(-l_comb,0.33333333);
            }
        }
    }
}


/*
 * 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;                         /* gradient threshold */
    float        dt=0.1;                      /* time step */
    float        n_iter_f;
    int          m;
    unsigned int k;
    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 iteration needed */
    n_iter_f= 0.75 * ( pow(R,1.33333333) / dt );

    /* round n_iter_f (must be an 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 {
        /* memory allocation 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++) {
            t_g=4;
            /* 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 affine morphological scale space */
                amss(ptr_out, ptr_in, n_row, n_col, dt, t_g);
                /* decrease threshold after the first iteration */
                if (m==0) t_g=1;
                /* 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(data_in);
    return 0;
}

---