Raw File
young.c
/*
 * Copyright (c) 2016, Anmol Popli <anmol.ap020@gmail.com>
 * All rights reserved.
 *
 * This program is free software: you can use, modify and/or
 * redistribute 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. You should have received a copy of this license along
 * this program. If not, see <http://www.gnu.org/licenses/>.
 */


/**
 * @file young.c
 * @brief Fast IIR approximation of gaussian filter using Young's approach
 *
 * @author ANMOL POPLI <anmol.ap020@gmail.com>
 *         PRAVIN NAIR  <sreehari1390@gmail.com>
 **/

#include "headersreq.h"
void convolve_young2D(int rows, int columns, int sigma, double complex** ip_padded);
void symmetric_padding(int rows,int columns,double complex **in,int w);
double bf[3], bb[3] , B;
static int w;
/**
 * \brief Convolve input array with 1D Causal filter
 *        (Young and van Vliet's algorithm) 
 * \param in        Pointer to input array
 * \param out       Pointer to output array
 * \param datasize  Input array size
 *
 * This routine performs constant time convolution of the
 * 1D input array of complex doubles with 1D Causal filter
 * of Young and van Vliet's algorithm. The 1D filter is an
 * IIR filter. 
 */

void convolve_youngCausal(double complex* in, double complex* out, int datasize) {
    
    int i, j;
    
    /* Compute first 3 output elements */  
    out[0] = B*in[0];
    out[1] = B*in[1] + bf[2]*out[0];
    out[2] = B*in[2] + (bf[1]*out[0]+bf[2]*out[1]);
    
    /* Recursive computation of output in forward direction using filter parameters bf and B */
    for (i=3; i<datasize; i++) {
        out[i] = B*in[i];
        for (j=0; j<3; j++) {
            out[i] += bf[j]*out[i-(3-j)];
        }
    }
         
}

/**
 * \brief Convolve input array with 1D AntiCausal filter
 *        (Young and van Vliet's algorithm) 
 * \param in        Pointer to input array
 * \param out       Pointer to output array
 * \param datasize  Input array size
 *
 * This routine performs constant time convolution of the
 * 1D input array of complex doubles with 1D AntiCausal filter
 * of Young and van Vliet's algorithm. The 1D filter is an
 * IIR filter. 
 */
void convolve_youngAnticausal(double complex* in, double complex* out, int datasize) {
    
    int i, j;

    /* Compute last 3 output elements */
    out[datasize-1] = B*in[datasize-1];
    out[datasize-2] = B*in[datasize-2] + bb[0]*out[datasize-1];
    out[datasize-3] = B*in[datasize-3] + (bb[0]*out[datasize-2]+bb[1]*out[datasize-1]);
    
    /* Recursive computation of output in backward direction using filter parameters bb and B */
    for (i=datasize-4; i>=w; i--) {
        out[i] = B*in[i];
        for (j=0; j<3; j++) {
            out[i] += bb[j]*out[i+(j+1)];
        }
    }
        
}

/**
 * \brief Convolve input array with 1D Gaussian filter
 *        (Young and van Vliet's algorithm) 
 * \param in        Pointer to input array
 * \param out       Pointer to output array
 * \param datasize  Input array size
 *
 * This routine performs constant time convolution of the
 * 1D input array of complex doubles with 1D Gaussian filter
 * using Young and van Vliet's algorithm. The input array is
 * first convolved with 1D Causal filter, the result of
 * which is convolved with 1D AntiCausal filter.
 */
void convolve_young1D(double complex* in, double complex* out, int datasize) {
    /** \brief Array to store output of Causal filter convolution */
    convolve_youngCausal(in, out, datasize);
    convolve_youngAnticausal(out, in, datasize);    
}

/**
 * \brief Apply 2D Gaussian filter to input image
 *        (Young and van Vliet's algorithm) 
 * \param rows      Image height
 * \param columns   Image width
 * \param sigma     Gaussian kernel standard deviation
 * \param ip_padded Pointer to input image
 * \param op_padded Pointer to output image
 *
 * This routine applies 2D Gaussian filter of s.d.
 * sigma to input image ip_padded of dimensions
 * rows x columns and computes output image op_padded.
 * 1D filter is first convolved along rows and then
 * along columns. The 1D convolution is performed using
 * Young and van Vliet's fast recursive algorithm.
 */
void convolve_young2D(int rows, int columns, int sigma, double complex** ip_padded) {
    
    /** \brief Filter radius */
    w = 3*sigma;
    /** \brief Filter parameter q */
    double q;
    if (sigma < 2.5)
        q = 3.97156 - 4.14554*sqrt(1-0.26891*sigma);
    else
        q = 0.98711*sigma - 0.9633;
    
    /** \brief Filter parameters b0, b1, b2, b3 */
    double b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q;
    double b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q;
    double b2 = -(1.4281*q*q + 1.26661*q*q*q);
    double b3 = 0.422205*q*q*q;
    
    /** \brief Filter parameters bf, bb, B */
    bf[0] = b3/b0; bf[1] = b2/b0; bf[2] = b1/b0;
    bb[0] = b1/b0; bb[1] = b2/b0; bb[2] = b3/b0;
    B = 1 - (b1+b2+b3)/b0;
    
    int i,j;

    symmetric_padding(rows,columns,ip_padded,w);
    /* Convolve each row with 1D Gaussian filter */
    double complex *out_t = calloc(columns+(2*w),sizeof(double complex));
    for (i=0; i<rows+2*w; i++) {
        convolve_young1D(ip_padded[i], out_t, columns+2*w);
    }
    free(out_t);
    double complex *intemp=calloc(rows+(2*w),sizeof(double complex)),*outtemp=calloc(rows+(2*w),sizeof(double complex));
    for (j=w; j<columns+w; j++) 
    {
        /* Convolve each column with 1D Gaussian filter */
        for (i=0;i<rows+(2*w);i++)
        {
            intemp[i]=ip_padded[i][j];      
        }
        convolve_young1D(intemp, outtemp, rows+2*w);
        /* Store the convolved column in row of output matrix*/
        for (i=0;i<rows+(2*w);i++)
        {
            ip_padded[i][j]=intemp[i];
        }
    }
    free(intemp);
    free(outtemp);
}
/**
 * \brief Apply symmetric padding to input image
 *        (Young and van Vliet's algorithm) 
 * \param rows      Image height
 * \param columns   Image width
 * \param in Pointer to input image padded with zeros
 *
 * This routine applies mirror boundary conditions
 * to input image which is zero padded i.e size of
 * input image will be [rows+2*w, columns+2*w]
 */
void symmetric_padding(int rows,int columns,double complex **in,int w){
    int i,j;
    double complex pixval;
    /* Rows 0 - w are duplicated to satisfy mirror bondary condition */
    for (i=0; i<w; i++) {
        for (j=0; j<columns; j++) {
            pixval = in[i+w][j+w];
            if (j < w)      {
                in[i+w][w-1-j] = pixval;
                in[w-1-i][w-1-j] = pixval;
            }
            if (j >= columns-w) {
                in[i+w][columns+w+(columns-j)-1] = pixval;
                in[w-1-i][columns+w+(columns-j)-1] = pixval;
            }
            in[w-1-i][j+w] = pixval;
        }
    }
    /* Rows rows-w - rows  are duplicated to satisfy mirror bondary condition */
    for (i=rows-w; i<rows; i++) {
        for (j=0; j<columns; j++) {
            pixval = in[i+w][j+w];
            if (j < w)      {
                in[i+w][w-1-j] = pixval;
                in[rows+w+(rows-i)-1][w-1-j] = pixval;
            }
            if (j >= columns-w) {
                in[i+w][columns+w+(columns-j)-1] = pixval;
                in[rows+w+(rows-i)-1][columns+w+(columns-j)-1] = pixval;
            }
            in[rows+w+(rows-i)-1][j+w] = pixval;
        }
    }
    /* Remaining rows are duplicated to satisfy mirror bondary condition*/ 
    for (i=w; i<rows-w; i++) {
        for (j=0; j<w; j++) {
            in[i+w][w-1-j]=in[i+w][j+w];
        }
    }
    for (i=w; i<rows-w; i++) {
        for (j=columns-w; j<columns; j++) {
            in[i+w][columns+w+(columns-j)-1] = in[i+w][j+w];
        }
    }

}

back to top