# poisson_lca_lib.c

/* * Copyright 2009, 2010 IPOL Image Processing On Line http://www.ipol.im/ * * 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/>. *//** * @file poisson_lca_lib.c * @brief   gradient, divergence  and Poisson routines using DFT.  * * @author Catalina Sbert <catalina.sbert@uib.es/> * @author Ana BelÃ©n Petro <anabelen.petro@uib.es/> */#include <stdlib.h>#include <math.h>#include <fftw3.h>#include "auxiliary_lib.h"/*#include "poisson_lca_lib.h"*//* M_PI is a POSIX definition */#ifndef M_PI/** macro definition for Pi */#define M_PI 3.14159265358979323846#endif                          /* !M_PI *//** * @brief compute the gradient of a 2D array using the discrete cosinus transform  * * We supose that the input image (N X M array) is extended simmetrically across its sides to a *  2N X 2M array. Then we obtain a function wich is even around (-0.5 , j) ,(N-0.5,j ),  * (i,-0.5) and (i, M-0.5). To compute the gradient we use the following definition *  of cosinus transform, (we write it in dimension 1 for dimension two is similar)*    *    @f$C(X)_k= 2\sum_{j=0}^{N-1} X_j cos( \frac{\pi k (j+1/2)}{N}) \f$* Note that with the symmetry C(X)_N=0.* * With this definition the sinus transform of the derivative is defined by *  * @f$S(X')_k=\frac{-\pi(k+1)}{N } C(X)_{k+1} , k=0,..., N-1 \f$* *  This new function is odd around -1 and even aroun N-1* Then we compute X'_k using the invers sinus  transform whose definition is ** @f$X'_k=\frac{2\sum_{j=0}^{N-1} S(X')_j sin (\frac{\pi(j+1)(k+1/2)}{N})}{2N} \f$** Note that the result is odd around -0.5 and around N-0.5.** @param data_x  partial derivative respect to the first variable* @param data_y  partial derivative respect to the second variable* @param data  input array* @param nx x-the dimensions of the array* @param ny y-the dimensions of the array** @return data_x, data_y  or NULL if an error occured*/ static double *gradient(double *data_x, double *data_y, double *data_in, size_t nx, size_t ny){	int x,y,dim,dim4;	double *ftout, *ftux, *ftuy;	fftw_plan p;	double normx, normy;	/*  pointers to the right and bottom neighbour values */	double *ptr_end, *ptr_in_x1,*ptr_in_y1, *ptr_x, *ptr_y;	/* sanity check */    if (NULL == data_in || NULL == data_x || NULL== data_y)    {        fprintf(stderr, "a pointer is NULL and should not be so\n");        abort();    }	dim=(int)nx*(int)ny;	dim4=4*dim;	/* allocate the memory of  fourier transform*/	ftout=(double*)fftw_malloc(sizeof(double)*dim);	ftux=(double*)fftw_malloc(sizeof(double)*dim);	ftuy=(double*)fftw_malloc(sizeof(double)*dim);	/* compute the cosinus transform of the input array*/	p=fftw_plan_r2r_2d((int) ny,(int) nx, data_in, ftout,FFTW_REDFT10,FFTW_REDFT10,FFTW_ESTIMATE);	fftw_execute(p); 	normx=M_PI/nx; normy=M_PI/ny;	/* compute the Fourier transform of partial derivatives respect to x and y*/	ptr_in_x1=ftout+1;	ptr_in_y1=ftout+nx;	ptr_x=ftux;	ptr_y=ftuy; 	for(y=0;y<(int)ny;y++)		for(x=0;x< (int) nx;x++)		{			if(x < (int) nx-1) *ptr_x=-normx*(x+1)*(*ptr_in_x1);			else *ptr_x=0.;			if(y < (int) ny-1) *ptr_y=-normy*(y+1)*(*ptr_in_y1);			else *ptr_y=0.;			ptr_x++;			ptr_y++;			ptr_in_x1++;			ptr_in_y1++;		}	/* compute the inverse Fourier transform to obtain the gradient of the input array*/	p=fftw_plan_r2r_2d((int) ny,(int) nx, ftux, data_x, FFTW_REDFT01,FFTW_RODFT01,FFTW_ESTIMATE);	fftw_execute(p);	p=fftw_plan_r2r_2d((int) ny, (int) nx, ftuy, data_y, FFTW_RODFT01,FFTW_REDFT01,FFTW_ESTIMATE);	fftw_execute(p);	/* normalize the inverse fourier transform*/         	ptr_end=data_x+dim;	ptr_x=data_x;	ptr_y=data_y;	while(ptr_x < ptr_end){		*ptr_x=*ptr_x/dim4;		*ptr_y=*ptr_y/dim4;		ptr_x++; ptr_y++;	}	/*deallocate the memory*/  	fftw_destroy_plan(p);	fftw_free(ftout);	fftw_free(ftux); fftw_free(ftuy);	return data_x; return data_y;}/** * @brief compute the Fourier transform of the divergence of a vector field using the *  @brief       discrete sinus/cosinus transform* *  We have two  N X M arrays V_1 and V_2. V_1 is  odd around (-0.5, j) and (N-0.5, j) and even around * (i,-0.5) and (i,M-0.5). V_2 is even  around (-0.5, j) and (N-0.5, j) and odd around * (i,-0.5) and (i,M-0.5).* We compute the Fourier transform of  divergence of the vector field V= (V_1, V_2), **     @f$div(V)=\frac{\partial V_1}{\partial x}+\frac{\partial V_2}{\partial y} \f$**   Due to the kind of symmetry to compute the Fourier transform of* @f$\frac{\partial V_1}{\partial x}\f$ ( respectively @f$frac{\partial V_2}{\partial y}\f$)* we use the sinus transform, which in dimension 1 is defined by** @f$S(X)_k=2\sum_{j=0}^{N-1} X_j sin(\frac{\pi(j+1/2)(k+1)}{N}) \f$ ** The cosinus transform of the derivative of X is defined by** @f$C(X')_k=\frac{\pi k}{N} S(X)_{k-1} , k=1,... N-1, C(X')_0=0 \f$ * * @param div  the cosinus transform of the divergence of a vector field* @param data_x the first component of the vector field* @param data_y the second component of the vector field* @param nx x-the dimensions of the array* @param ny y-the dimensions of the array ** @return div or NULL */static double *divergence(double *div, double *data_x, double *data_y, size_t nx, size_t ny){	int x,y, dim;	double *ftout,*ptr_div,*ptr_x1;	fftw_plan p;	double normx, normy;    /* sanity check*/	if (NULL == div || NULL == data_x || NULL== data_y)    {        fprintf(stderr, "a pointer is NULL and should not be so\n");        abort();    }	dim=(int)nx*(int)ny;	normx=M_PI/nx; normy=M_PI/ny;	/* allocate the Fourier transform*/	ftout=(double*)fftw_malloc(sizeof(double)*dim);	/* compute the Fourier transform of the first component of the vector field*/	p=fftw_plan_r2r_2d((int) ny,(int) nx, data_x, ftout,FFTW_REDFT10,FFTW_RODFT10,FFTW_ESTIMATE);	fftw_execute(p);         /*compute the fourier transform of the derivative respect to x of the first component*/    	ptr_x1=ftout;	ptr_div=div;	for(y=0;y< (int) ny; y++)		for(x=0; x< (int) nx; x++){			if(x==0) *ptr_div=0.;			else {				*ptr_div= normx*x* *(ptr_x1-1);			}			ptr_div++;			ptr_x1++;		}    /* compute the Fourier transform of the second component of the vector field*/	p=fftw_plan_r2r_2d((int) ny, (int) nx,data_y, ftout,FFTW_RODFT10,FFTW_REDFT10,FFTW_ESTIMATE);	fftw_execute(p); 	/* compute the fourier transform of the derivative respect to y of the second component and the divergence*/   	ptr_x1=ftout;	ptr_div=div+nx;	for(y=1;y< (int) ny; y++)		for(x=0; x< (int) nx; x++){			*ptr_div+=normy*y*(*ptr_x1);			ptr_div++;			ptr_x1++;		}    /* deallocate the memory*/	fftw_destroy_plan(p);	fftw_free(ftout);	return div;}/*** @brief Solve a Poisson equation using FFTW** Given a vector filed V=(ux,uy), this routine solves the Poisson equation* @f$\Delta u= div(V)\f$ with homogeneous Neumann boundary conditions*** @li compute  the  Fourier transform of the divergence of the vector field V, C(di.* @li compute the cosinus transform of the laplacian*  @f$C(\Delta u)_{ij}=-[(\frac{\pi i}{nx})Â²+(\frac{\pi j}{ny})Â²] C(u)_{ij}\f$* @li The cosinus transform of the solution of the Poisson equation is* @f\$ C(u)_{ij}=-\frac{C(div)}{(\frac{\pi i}{nx})Â²+(\frac{\pi j}{ny})Â²}* @li the solution is computed by the inverse cosinus transform ** @param ux first component of the vector* @param uy second component of the vector* @param nx x-dimension of the array* @param ny y-dimension of the array* * @return data_out*/double *solve_poisson_equation(double *data_out, double *ux, double *uy, size_t nx, size_t ny){	double *div;	double normx, normy, A;	int x,y,l, dim, dim4;	fftw_plan p;	dim=(int)nx*(int)ny;	dim4=4*dim;	div=(double*)fftw_malloc(sizeof(double)*dim);	/* compute the fourier transform of the divergence of the vector (ux,uy)*/	(void) divergence(div, ux,uy, nx,ny);	normx=M_PI/nx; normy=M_PI/ny;	normx*=normx; normy*=normy;	/* compute the fourier transform of the solution of the Poisson equation*/     	for(y=0;y< (int) ny;y++){		l=y*(int)nx;		for(x=0;x< (int) nx;x++){			if( x==0 && y==0) div[0]=0.;			else{                A=(normx*x*x+normy*y*y);                div[l+x]=-div[l+x]/A;			}		}	}    /* compute the inverse fourier transform*/	p=fftw_plan_r2r_2d((int) ny, (int) nx, div, data_out,FFTW_REDFT01, FFTW_REDFT01, FFTW_ESTIMATE);	fftw_execute(p); 	for(l=0; l< dim; l++) data_out[l]/=dim4;	return  data_out;}/*** @brief  Linear option for modifying the gradient, solving a Poisson equation using the Fourier Transform**  The input array is processed as follow:** @li compute the gradient of the input array* @li Modify the gradient: amplifying by a factor the gradient in the interior of the dark zones of the image* @li solve the Poisson equation aplying FFT** @param data_in input array* @param T threshold for the dark regions* @param a amplification factor for the gradient* @param nx x-dimension of the array* @param ny y-dimension of the array* * @return data_out*/ double  *poisson_local(double *data_out, double *data_in, double T, double a, size_t nx, size_t ny){	double *ux,*uy;	double *ptr_x,*ptr_y, *ptr_data,  *ptr_xr, *ptr_xl,*ptr_yt,*ptr_yb;	int x,y, dim;	if (NULL == data_out || NULL == data_in )	{        fprintf(stderr, "a pointer is NULL and should not be so\n");        abort();    }	dim=(int)nx*(int)ny;	/* allocate the memory for the Fourier transform*/      	ux=(double*)fftw_malloc(sizeof(double)*dim);	uy=(double*)fftw_malloc(sizeof(double)*dim);	/* compute the gradient of the input array*/	gradient(ux, uy, data_in, nx,ny);          	/* multiply the gradient by a factor a in the interior of the dark zones which are defined by a threshold T*/   	ptr_x=ux; ptr_y=uy;	ptr_data=data_in;	ptr_xr=data_in+1;	ptr_xl=data_in-1;	ptr_yt=data_in-nx;	ptr_yb=data_in+nx;	for(y=0; y < (int) ny; y++)		for(x=0; x < (int) nx; x++){			if( x >0  && x < (int) nx-1  && y>0 && y < (int) ny-1)			if(*ptr_xl <= T && *ptr_data <= T && *ptr_xr <= T && *ptr_yt <= T && *ptr_yb <= T){				*ptr_x=a*(*ptr_x); *ptr_y=a* (*ptr_y);			}			ptr_data++; ptr_x++; ptr_y++; ptr_xr++; ptr_xl++; ptr_yt++; ptr_yb++;		}	  	solve_poisson_equation(data_out,ux,uy, nx,ny);	fftw_free(ux); fftw_free(uy);    return data_out;}/*** @brief  Power option for modifying the gradient, solving a Poisson equation using the Fourier Transform**  The input array is processed as follow:** @li compute the gradient of the input array* @li Modify the gradient: power by alpha the gradient* @li solve the Poisson equation** @param data_in input array* @param alpha Power function for the gradient* @param nx x-dimension of the array* @param ny y-dimension of the array** @return data_out*/double  *poisson_global(double *data_out, double *data_in, double a, size_t nx, size_t ny){	double *ux,*uy;	int l,dim;	double A;	if (NULL == data_out || NULL == data_in )    {        fprintf(stderr, "a pointer is NULL and should not be so\n");        abort();    }	dim=(int)nx*(int)ny;	/* allocate the memory for the Fourier transform*/	ux=(double*)fftw_malloc(sizeof(double)*dim);	uy=(double*)fftw_malloc(sizeof(double)*dim);	/* compute and modify the gradient of the input array*/	gradient(ux, uy, data_in, nx,ny);	for(l=0;l< dim;l++){		A=fabs(ux[l])+fabs(uy[l])+0.1;		A=pow(A, a-1.);		ux[l]*=A; 		uy[l]*=A;	}	solve_poisson_equation(data_out, ux, uy, nx,ny);	fftw_free(ux); fftw_free(uy);	return data_out;}/** * @brief  The algorithm applied to the gray intensity and the color channels are  computed  *  such that the ratios R/G/B are preserved. * * * @param data input array * @param T the threshold to compute the dark regions * @param a the gradient amplification factor  * @param alpha Power function for the gradient * @param nx x-dimension of the array * @param ny y-dimension of the array * @param channel the number of channels of the input image * @param s the percentage of pixels saturated to minimum/maximum value * * @return data_norm  simplest color balance of the input with s% of saturation * @return local_out the result of the local option of the algorithm * @return global_out the result of the global option of the algoritm */double *poisson_gray(double *data_norm, double *local_out, double *global_out, 					 double *data, double T, double a, double alpha, size_t nx, size_t ny, int channel, float s){	int dim;	double *data_gray,*data_gray1;	dim=(int)nx*(int)ny;	data_gray = (double*) malloc(dim* sizeof(double));	data_gray1 = (double*) malloc(dim * sizeof(double));	if(channel == 1){		simplest_color_balance(data_norm, data, dim, s);		poisson_local(local_out, data_norm, T, a, nx,ny);		simplest_color_balance(local_out, local_out, dim,s);		poisson_global(global_out, data_norm, alpha, nx,ny);		simplest_color_balance(global_out, global_out, dim,s);	}	else{		gray_intensity(data_gray, data, dim);		simplest_color_balance(data_gray1, data_gray, dim,s);		color(data_norm, data, data_gray, data_gray1, dim);		poisson_local(data_gray, data_gray1, T, a, nx,ny);		simplest_color_balance(data_gray, data_gray, dim,s);		color(local_out, data_norm, data_gray1, data_gray, dim);		poisson_global(data_gray, data_gray1, alpha, nx,ny);		simplest_color_balance(data_gray, data_gray, dim,s);		color(global_out, data_norm, data_gray1, data_gray, dim);	}	free(data_gray); free(data_gray1);	return data_norm;	return local_out;	return global_out;}/** * @brief  The algorithm applied to  the color channels independently  * * @param data input array * @param T the threshold to compute the dark regions * @param a the gradient amplification factor  * @param alpha Power function for the gradient * @param nx x-dimension of the array * @param ny y-dimension of the array * @param channel the number of channels of the input image * @param s the percentage of pixels saturated to minimum/maximum value * * @return data_norm  simplest color balance of the input with s% of saturation * @return local_out the result of the local option of the algorithm * @return global_out the result of the global option of the algoritm */double *poisson_rgb(double *data_norm, double *local_out, double *global_out, 					double *data, double T, double a, double alpha, size_t nx, size_t ny, int channel, float s){	int dim, n;	dim=(int)nx*(int)ny;	for(n=0;n<channel;n++) simplest_color_balance(data_norm+n*dim, data+n*dim, dim, s);	for(n=0;n<channel;n++){		poisson_local(local_out+n*dim, data_norm+n*dim, T, a, nx,ny);		simplest_color_balance(local_out+n*dim, local_out+n*dim, dim, s);	}	for(n=0;n<channel;n++){		poisson_global(global_out+n*dim, data_norm+n*dim, alpha, nx,ny);		simplest_color_balance(global_out+n*dim, global_out+n*dim, dim, s);	}	return data_norm;	return local_out;	return global_out;	   }