laplacian, DFT and Poisson routines More...

#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <fftw3.h>
#include "retinex_pde_lib.h"

Include dependency graph for retinex_pde_lib.c:

Defines
#define	M_PI 3.14159265358979323846
Functions
static float *	discrete_laplacian_threshold (float data_out, const float data_in, size_t nx, size_t ny, float t)
	compute the discrete laplacian of a 2D array with a threshold
static double *	cos_table (size_t size)
	compute a cosinus table
static float *	retinex_poisson_dct (float *data, size_t nx, size_t ny, double m)
	perform a Poisson PDE in the Fourier DCT space
float *	retinex_pde (float *data, size_t nx, size_t ny, float t)
	retinex PDE implementation

Detailed Description

laplacian, DFT and Poisson routines

Author:: Nicolas Limare <nicolas.limare@cmla.ens-cachan.fr>

Definition in file retinex_pde_lib.c.

Define Documentation

#define M_PI 3.14159265358979323846

macro definition for Pi

Definition at line 36 of file retinex_pde_lib.c.

Function Documentation

static double* cos_table ( size_t size ) [static]

compute a cosinus table

Allocate and fill a table of n values cos(i Pi / n) for i in [0..n[.

Parameters:

size

the table size

Returns:: the table, allocated and filled

Definition at line 147 of file retinex_pde_lib.c.

Here is the caller graph for this function:

static float* discrete_laplacian_threshold	(	float *	data_out,
		const float *	data_in,
		size_t	nx,
		size_t	ny,
		float	t
	)			`[static]`

compute the discrete laplacian of a 2D array with a threshold

This function computes the discrete laplacian, ie $(F_{i - 1, j} - F_{i, j}) + (F_{i + 1, j} - F_{i, j}) + (F_{i, j - 1} - F_{i, j}) + (F_{i, j + 1} - F_{i, j})$ . On the border, differences with "outside of the array" are 0. If the absolute value of difference is < t, 0 is used instead.

This step takes a significant part of the computation time, and needs to be fast. In that case, we observed that (with our compiler and architecture):

pointer arithmetic is faster than data[i]
if() is faster than ( ? : )

Parameters:

	data_out	output array
	data_in	input array
	nx,ny	array size
	t	threshold

Returns:: data_out

Todo:: split corner/border/inner

Definition at line 71 of file retinex_pde_lib.c.

Here is the caller graph for this function:

float* retinex_pde	(	float *	data,
		size_t	nx,
		size_t	ny,
		float	t
	)

retinex PDE implementation

This function solves the Retinex PDE equation with forward and backward DCT.

The input array is processed as follow:

a discrete laplacian is computed with a threshold;
this discrete laplacian array is symmetrised in both directions;
this data is transformed by forward DFT (both steps can be handled by a simple DCT);
the DFT data is modified by $\hat{u}(i, j) = \frac{\hat{F}(i, j)} {4 - 2 \cos(\frac{i \pi}{n_x}) - 2 \cos(\frac{j \pi}{n_y})}$ ;
this data is transformed by backward DFT.

Parameters:

	data	input/output array
	nx,ny	dimension
	t	retinex threshold

Returns:: data, or NULL if an error occured

Definition at line 268 of file retinex_pde_lib.c.

Here is the call graph for this function:

Here is the caller graph for this function:

static float* retinex_poisson_dct	(	float *	data,
		size_t	nx,
		size_t	ny,
		double	m
	)			`[static]`

perform a Poisson PDE in the Fourier DCT space

$ u(i, j) = F(i, j) * m / (4 - 2 cos(i PI / nx) - 2 cos(j PI / ny)) $ if $(i, j) \neq (0, 0)$ , $ u(0, 0) = 0 $

When this function is successively used on arrays of identical size, the trigonometric computation is redundant and could be kept in memory for a faster code. However, in our use case, the speedup is marginal and we prefer to recompute this data and keep the code simple.

Parameters:

	data	the dct complex coefficients, of size nx x ny
	nx,ny	data array size
	m	global multiplication parameter (DCT normalization)

Returns:: the data array, updated

Definition at line 190 of file retinex_pde_lib.c.

Here is the call graph for this function:

Here is the caller graph for this function:

retinex_pde_lib.c File Reference

Defines

Functions

Detailed Description

Define Documentation

Function Documentation