https://doi.org/10.5201/ipol.2017.198
gti.c
/**
* @file gti.c
*
* @brief Gaussian Texture Inpainting
*
* This code was written for IPOL and is based on the paper
* "Texture Inpainting Using Efficient Gaussian Conditional Simulation"
* (Bruno Galerne, Arthur Leclaire), SIAM Journal on Imaging Sciences, 2017.
* The algorithm is thoroughly discussed in the IPOL companion paper
* "An Algorithm for Gaussian Texture Inpainting"
* (Bruno Galerne, Arthur Leclaire)
*
* The program performs conditional simulation of a Gaussian texture model on the mask, knowing the values of the texture on a border of the mask.
* The Gaussian texture model is estimated on the mask complement.
* The linear system involved in the Gaussian conditional simulation is solved with a conjugate gradient descent.
*
* For more information consult the corresponding paper.
*
* NB: For the sake of clarity, we added in this code (with square brackets [.])
* references to equations or pseudo-codes given in the IPOL paper.
*
* @author Bruno Galerne, Arthur Leclaire
*
* @version 1.3
*
* @section LICENSE
*
* Copyright (c) 2017, Bruno Galerne, Arthur Leclaire
*
* This program is free software: you can use, modify and/or
* redistribute it under the terms of the simplified BSD
* License. You should have received a copy of this license along
* this program. If not, see
* <http://www.opensource.org/licenses/bsd-license.html>.
*
*/
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <png.h>
#include <fftw3.h>
#include <complex.h>
#include <math.h>
#ifdef _OPENMP
#include <omp.h>
#else
#define omp_get_max_threads() 1
#endif
#include "io_png.h"
/**
*
* @brief Implementation of a malloc that quits when the return of malloc is null.
*
*/
void* xmalloc(size_t size){
void* mem = malloc(size);
if( mem == NULL){
perror("Memory allocation error");
exit(EXIT_FAILURE);
}
return mem;
}
#define M_PI 3.141592653589793
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
/**
* @brief Height of the image.
*/
int height;
/**
* @brief Width of the image.
*/
int width;
/**
* @brief Double Height of the image.
*/
int height2;
/**
* @brief Double Width of the image.
*/
int width2;
/**
* @brief Size of Fourier transforms
*/
int ftsize;
/**
* @brief fftw global variables
*
* fftwf_complex is a type of float-precision complex number arrays
* (output of the Fourier transform)
* fftwf_plan is a type of object that contains the FFT plan, that is,
* the stategy that FFTW will adopt to compute the Fourier
* transform for the current image size.
*/
// output of Fourier transform
static fftwf_complex *fftw_out;
// input of Fourier transform
static float *fftw_in;
// plans for direct and inverse FFT
static fftwf_plan fftw_plan_direct,fftw_plan_inverse;
// float arrays to store separately real and imaginary parts of FFT
float *re,*im,*ret,*imt,*retmp,*imtmp;
int numth; // number of threads
static char *wisdomfilename = "/var/tmp/.wisdom";
static char wisdomfilenamestatic[1000];
/**
* @brief Load wisdom file for FFT computation.
*
* This function loads a wisdom file, which contains previously
* computed plans that serve to compute the FFT in an optimal way
* for a given size of image.
* If the wisdom file does not contain an FFT plan that matches
* the current image size, then the first FFT routine will take time
* to compute an optimal plan for this size, and will save it
* in the wisdom file.
* This function creates a wisdom file if it does not exist yet.
*/
void evoke_wisdom(void)
{
if (!wisdomfilename) {
char *homedir = getenv("HOME");
if (!homedir)
return; // no wisdom for you
wisdomfilename = wisdomfilenamestatic;
size_t n = strlen(homedir);
// assert(n + 10 < 1000);
memcpy(wisdomfilename, homedir, n);
memcpy(wisdomfilename + n, "/.wisdomf", 9);
}
FILE *f = fopen(wisdomfilename, "r");
if (f) {
fftwf_import_wisdom_from_file(f);
fclose(f);
}
else {
f = fopen(wisdomfilename, "w");
if (!f) {
perror("could not create wisdom file");
exit(EXIT_FAILURE);
}
fclose(f);
fprintf(stderr, "created wisdom file \"%s\"\n", wisdomfilename);
}
}
/**
* @brief Export FFT wisdom into wisdom file
*
* This function saves the optimal FFT plan used for the
* current image size into the wisdom file.
* It may then be used for future FFT computations with
* the same image size.
*/
void bequeath_wisdom(void)
{
// assert(wisdomfilename);
FILE *f = fopen(wisdomfilename, "w");
if (f)
{
fftwf_export_wisdom_to_file(f);
fclose(f);
}
else {
perror("could not bequeath wisdom");
exit(EXIT_FAILURE);
}
}
/**
* @brief Initialize FFT variables and compute the FFT plan.
*
* NB: This function acts on global variables.
*/
static void init_fftw()
{
#ifdef _OPENMP
fftwf_init_threads();
fftwf_plan_with_nthreads(numth);
printf("numth=%d\n",numth);
#endif
fftw_in = (float*) xmalloc(sizeof(float) * width2*height2);
fftw_out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * ftsize);
evoke_wisdom();
// The two following lines compute an optimal FFT plan for the current image size.
fftw_plan_direct = fftwf_plan_dft_r2c_2d(height2,width2, fftw_in, fftw_out, FFTW_MEASURE);
fftw_plan_inverse = fftwf_plan_dft_c2r_2d(height2,width2, fftw_out, fftw_in, FFTW_MEASURE);
// The input image is real so we use r2c, c2r transforms which are more efficient because
// they do not compute redundant Fourier coefficients.
// r2c -> real to complex
// c2r -> complex to real
// The flag FFTW_MEASURE allows to compute an optimal plan for FFT computation with the current image size.
// Warning: When using FFTW_MEASURE in the following, fftw_in must be assigned after plan allocation.
}
/**
* @brief Compute direct FFT
*
* NB: fftw_in, fftw_out, and fftw_plan_direct are global variables
*/
static void fftw_fft2d(in_re, out_re, out_im)
float *in_re, *out_re, *out_im;
{
int i;
//#pragma omp parallel for shared(fftw_in,in_re) private(i)
for (i=0;i<width2*height2;i++)
fftw_in[i] = in_re[i];
fftwf_execute(fftw_plan_direct);
if (out_re)
//#pragma omp parallel for shared(out_re,fftw_out) private(i)
for (i=0;i<ftsize;i++)
out_re[i] = fftw_out[i][0];
if (out_im)
//#pragma omp parallel for shared(out_im,fftw_out) private(i)
for (i=0;i<ftsize;i++)
out_im[i] = fftw_out[i][1];
}
/**
* @brief Compute inverse FFT
*
* NB: fftw_in, fftw_out, and fftw_plan_inverse are global variables
*/
static void fftw_fft2d_inv(in_re, in_im, out_re)
float *in_re, *in_im, *out_re;
{
int i;
float norm;
//#pragma omp parallel for shared(fftw_out,in_re) private(i)
for (i=0;i<ftsize;i++)
fftw_out[i][0] = in_re[i];
if (in_im)
//#pragma omp parallel for shared(fftw_out,in_im) private(i)
for (i=0;i<ftsize;i++)
fftw_out[i][1] = in_im[i];
fftwf_execute(fftw_plan_inverse);
norm = 1./(float)(width2*height2);
#pragma omp parallel for shared(fftw_in) private(i)
for (i=0;i<width2*height2;i++)
fftw_in[i] *= norm;
if (out_re)
//#pragma omp parallel for shared(out_re,fftw_in) private(i)
for (i=0;i<width2*height2;i++)
out_re[i] = fftw_in[i];
}
/**
* @brief Compute inverse FFT
*/
static void term_fftw()
{
bequeath_wisdom();
fftwf_free(fftw_in); fftwf_free(fftw_out);
fftwf_destroy_plan(fftw_plan_direct); fftwf_destroy_plan(fftw_plan_inverse);
}
/**
* @brief Calculates the position on a torus.
*
* @param i The y coordinate
* @param j The x coordinate
* @param c The channel number
* @return The position in an array of the coordinate (i,j) on channel c
*/
inline int Ic(int i, int j, int c){
if(i < 0)
i += height2;
else if( i >= height2)
i = i - height2;
if(j < 0)
j += width2;
else if( j >= width2)
j = j - width2;
return c*height2*width2 + i*width2 + j;
}
/**
* @brief Sets the mask from an image.
*
* The mask will be obtained by looking at the area that has RGB values (254,254,254) or higher.
*
* @param mask The mask that will be extracted.
* @param im The image from which the mask will be extracted
*
* NB: mask must be allocated with a twice larger size than im
* */
void setMask(float *mask, const float *im){
int i, j, c;
int val;
for(i = 0; i < height2; i++){
for(j = 0; j < width2; j++){
if (i<height && j<width){
if(im[i*width + j] > 253 && im[height*width + i*width + j] > 253 && im[2*height*width + i*width + j] > 253)
val = 0; // 0 value for a masked pixel
else
val = 1;
for(c = 0; c < 3; c++)
mask[Ic(i,j,c)] = val;
}
else {
for(c = 0; c < 3; c++)
mask[Ic(i,j,c)] = 0;
}
}
}
}
/**
* @brief Sets the value of an array to a given value.
*
* @param A The array to be written on
* @param value The value to write to the array
* @param size The size of the array
*/
void setValue(float *A, float value, int size){
int i;
//#pragma omp parallel for shared(A) private(i)
for(i = 0; i < size; i++)
A[i] = value;
}
/**
* @brief Draw an ADSN realization
*
* @param meanv 3-channel mean (0 if NULL)
* @param v 3-channel output
*
* See [Algorithm 2], [Equation (2)] of the IPOL paper.
*
* NB: re, im, ret, imt are global variables.
* The texton FFT is precomputed in the global variables ret, imt
*/
void adsn(float *meanv,float *v){
int adr,ch;
float a,b,c,d;
float *w, *rew, *imw;
w = (float*) xmalloc(sizeof(float)*width2*height2);
rew = (float*) xmalloc(sizeof(float)*ftsize);
imw = (float*) xmalloc(sizeof(float)*ftsize);
// Draw a white noise
for (adr=0;adr<width2*height2;adr++) {
a = ((float)rand())/RAND_MAX; b = ((float)rand())/RAND_MAX;
w[adr] = sqrt(-2.0*log(a))*cos(2.0*M_PI*b);
}
// Convolve with texton (whose FFT has been precomputed)
fftw_fft2d(w,rew,imw);
for (ch=0;ch<3;ch++) {
#pragma omp parallel for shared(re,im,rew,imw,ret,imt) private(adr,a,b,c,d)
for (adr=0;adr<ftsize;adr++) {
a = ret[adr+ch*ftsize]; b = imt[adr+ch*ftsize];
c = rew[adr]; d = imw[adr];
re[adr] = a*c-b*d;
im[adr] = b*c+a*d;
}
fftw_fft2d_inv(re,im,v+ch*height2*width2);
}
// Add the mean component
if (meanv!=NULL)
for (ch=0;ch<3;ch++)
//#pragma omp parallel for shared(v) private(adr)
for (adr=0;adr<width2*height2;adr++)
v[adr+ch*width2*height2] += meanv[ch];
// free(re); free(im);
free(w); free(rew); free(imw);
}
/**
* @brief Convolve single-channel u with monochannel w
*
* @param u 1-channel input
* @param w 1-channel input
* @param v 1-channel output
* NB: re, im are global variables.
*/
void convol_single_channel(float *u,float *w,float *v){
int adr;
float a,b,c,d;
float *rew, *imw;
rew = (float*) xmalloc(sizeof(float)*ftsize);
imw = (float*) xmalloc(sizeof(float)*ftsize);
fftw_fft2d(w,rew,imw);
fftw_fft2d(u,re,im);
#pragma omp parallel for shared(re,im,rew,imw) private(adr,a,b,c,d)
for (adr=0;adr<ftsize;adr++) {
a = re[adr]; b = im[adr];
c = rew[adr]; d = imw[adr];
re[adr] = a*c-b*d;
im[adr] = b*c+a*d;
}
fftw_fft2d_inv(re,im,v);
// printf("a=%f,b=%f,c=%f,d=%f\n",a,b,c,d);
free(rew); free(imw);
}
/**
* @brief Convolve with covariance (with possible restriction)
*
* @param u 3-channel input image
* @param v 3-channel output
* @param dom 1-channel indicator function for restriction
*
* This function allows to convolve the input image u with
* the matrix texton associated to the Gaussian model
* [see equations (4) and (6)]
* and then restrict the result on the domain {dom=1}
* (i.e. set the values to 0 outside this domain)
*
* In the inpainting context, {dom=1} is the set of conditioning points.
*
* The texton FFT is precomputed in the global variables ret, imt
* and thus is not given as an argument of this function.
*
* See [Algorithm 4] of the IPOL paper.
*
* NB: re, im, ret, imt, retmp, imtmp are global variables.
*/
void convcov(float *u,float *v,int *dom){
// NB: can be used with u=v
int adr,ch;
float a,b,c,d;
setValue(retmp,0.,ftsize);
setValue(imtmp,0.,ftsize);
/* The covariance operator is given by the double convolution of [Equation (6)].
See also [Algorithm 3] */
// First convolution with transposed texton [ \tilde{t}_v^T ]
for (ch=0;ch<3;ch++) {
fftw_fft2d(u+ch*height2*width2,re,im);
#pragma omp parallel for shared(re,im,ret,imt,retmp,imtmp) private(adr,a,b,c,d)
for (adr=0;adr<ftsize;adr++) {
a = re[adr]; b = im[adr];
c = ret[adr+ch*ftsize]; d = imt[adr+ch*ftsize];
retmp[adr] += a*c+b*d;
imtmp[adr] += b*c-a*d;
}
}
// Second convolution with texton [ t_v ]
for (ch=0;ch<3;ch++) {
#pragma omp parallel for shared(re,im,retmp,imtmp,ret,imt) private(adr,a,b,c,d)
for (adr=0;adr<ftsize;adr++) {
a = retmp[adr]; b = imtmp[adr];
c = ret[adr+ch*ftsize]; d = imt[adr+ch*ftsize];
re[adr] = a*c-b*d;
im[adr] = b*c+a*d;
}
fftw_fft2d_inv(re,im,v+ch*height2*width2);
}
// Extract the values on the domain dom.
if (dom != NULL)
for (ch=0;ch<3;ch++)
for (adr=0;adr<width2*height2;adr++)
if (dom[adr]==0)
v[adr+ch*height2*width2] = 0.;
// printf("a=%f,b=%f,c=%f,d=%f\n",a,b,c,d);
}
/**
* @brief Get conditioning points
*
* @param mask Input Mask (0 for masked pixels)
* @param cond Output Conditioning points (1 for conditioning points)
* @param w Thickness of conditioning border
*/
void get_conditioning_points(const float *mask,int *cond,int condw){
int i,j;
float *maskc,*tmp;
tmp = (float *) xmalloc(sizeof(float)*width2*height2);
maskc = (float *) xmalloc(sizeof(float)*width2*height2);
for (i=0;i<height2;i++){
for (j=0;j<width2;j++){
tmp[Ic(i,j,0)] = 0.;
if (i<height && j<width)
maskc[Ic(i,j,0)] = 1-mask[Ic(i,j,0)];
else
maskc[Ic(i,j,0)] = 0;
}
}
for (i=-condw;i<=condw;i++)
for (j=-condw;j<=condw;j++)
tmp[Ic(i,j,0)] = 1.;
convol_single_channel(maskc,tmp,tmp);
for (i=0;i<height;i++)
for (j=0;j<width;j++)
if (i<height && j<width)
cond[Ic(i,j,0)] = (int)((tmp[Ic(i,j,0)]>0.5)&&(mask[Ic(i,j,0)]>0.5));
else
cond[Ic(i,j,0)] = 0;
}
/**
* @brief Estimate ADSN model
*
* See [Algorithm 1] and [Equation (1)] in the IPOL paper.
*
* @param u 3-channel input
* @param mask Mask (0 for masked pixels)
* @param meanu Image mean color (3 values)
* @param t Output texton
*/
void estimate_adsn_model(float *u, const float *mask,float *meanu, float *t){
int i,j,c;
float cardm,cardcm;
// Number of pixels in the mask
cardm = 0.;
for (i=0;i<height;i++)
for (j=0;j<width;j++)
if (mask[Ic(i,j,0)]<0.5)
cardm++;
cardcm = width*height-cardm;
printf("nb masked pixels=%f\n",cardm);
// Mean value
for (c=0;c<3;c++) {
meanu[c] = 0.;
for (i=0;i<height;i++)
for (j=0;j<width;j++)
if (mask[Ic(i,j,c)]>0.5)
meanu[c] += u[Ic(i,j,c)];
meanu[c] /= cardcm;
}
printf("meanu = (%f,%f,%f)\n",meanu[0],meanu[1],meanu[2]);
// Texton
for (c=0;c<3;c++)
for (i=0;i<height2;i++)
for (j=0;j<width2;j++)
if (mask[Ic(i,j,c)]>0.5)
t[Ic(i,j,c)] = (u[Ic(i,j,c)]-meanu[c])/sqrt(cardcm);
else
t[Ic(i,j,c)] = 0.;
}
/**
* @brief Conjugate gradient descent on normal equations for kriging system
*
* @param cond Conditioning points (1 for conditioning points)
* @param rhs Right-hand side of the system
* @param v Output
* @param ep Stopping criterion on L^2 norm of residual
* @param imax Maximum number of iterations
*
* See [Algorithm 5] of the IPOL paper.
*
*/
void cgd(int *cond,float *rhs,float *x,float ep,int imax){
int iter = 0, adr;
float *p,*q,*r;
float rn,rn2,rn2old,alpha,beta;
p = (float*) xmalloc(sizeof(float)*width2*height2*3);
q = (float*) xmalloc(sizeof(float)*width2*height2*3);
r = (float*) xmalloc(sizeof(float)*width2*height2*3);
// initialization
setValue(x,0.,width2*height2*3);
convcov(rhs,r,cond);
rn2 = 0.;
for (adr=0;adr<width2*height2*3;adr++) {
rn2 += r[adr]*r[adr];
p[adr] = r[adr];
}
rn = sqrt(rn2);
while (iter++<imax && rn>ep) {
printf("Iteration %d, res norm = %f\n",iter,rn);
// [compute A^T A d_k ]
convcov(p,q,cond);
convcov(q,q,cond);
// [compute \alpha_k ]
alpha = 0;
#pragma omp parallel for shared(p,q) private(adr) reduction(+:alpha)
for (adr=0;adr<width2*height2*3;adr++)
alpha += p[adr]*q[adr];
alpha = rn2/alpha;
/* Update variables. Correspondence with the notation of IPOL paper:
[ \psi_k ] <-> x
[ r_k ] <-> r
[ \|r_k\|^2 ] <-> rn2
*/
rn2old = rn2; rn2 = 0.;
#pragma omp parallel for shared(x,r) private(adr) reduction(+:rn2)
for (adr=0;adr<width2*height2*3;adr++) {
x[adr] += alpha*p[adr];
r[adr] -= alpha*q[adr];
rn2 += r[adr]*r[adr];
}
rn = sqrt(rn2);
// Update [d_{k+1}] <-> p
beta = rn2/rn2old;
#pragma omp parallel for shared(p,r) private(adr)
for (adr=0;adr<width2*height2*3;adr++)
p[adr] = r[adr] + beta*p[adr];
}
free(p);free(q);free(r);
}
/**
* @brief Gaussian texture inpainting
*
* @param image The image to inpaint
* @param mask The mask which indicates the inpainting domain (0 for masked values)
* @return The inpainted image for which a memory has been dynamically allocated.
* @sa THRESHOLD MAXITER
*
* See [Algorithm 6] of the IPOL paper.
*
* NB: re,im,ret,imt,retmp,imtmp are global variables that contain
* the output of Fourier transforms (with separate real and imaginary parts)
*/
float* gausstexinpaint(const float *image, const float *mask,float ep, int imax, int condw){
float *u, *v, *z, *meanu, *t, *rhs;
int *cond;
int i,j,c;
u = (float*) xmalloc(sizeof(float)*width2*height2*3);
v = (float*) xmalloc(sizeof(float)*width2*height2*3);
z = (float*) xmalloc(sizeof(float)*width2*height2*3);
meanu = (float*) xmalloc(sizeof(float)*3);
t = (float*) xmalloc(sizeof(float)*width2*height2*3);
ret = (float*) xmalloc(sizeof(float)*ftsize*3);
imt = (float*) xmalloc(sizeof(float)*ftsize*3);
re = (float*) xmalloc(sizeof(float)*ftsize);
im = (float*) xmalloc(sizeof(float)*ftsize);
retmp = (float*) xmalloc(sizeof(float)*ftsize);
imtmp = (float*) xmalloc(sizeof(float)*ftsize);
rhs = (float*) xmalloc(sizeof(float)*width2*height2*3);
cond = (int*) xmalloc(sizeof(int)*width2*height2);
// initialize FFT structures
init_fftw();
// copy initial image
for (i=0;i<height;i++)
for (j=0;j<width;j++)
for (c=0;c<3;c++)
u[Ic(i,j,c)] = image[c*height*width + i*width + j];
/* Estimate Gaussian model [Algorithm 1] */
estimate_adsn_model(u,mask,meanu,t);
/* Compute Fourier transform of the texton (once and for all)
It is stored in the global variables ret, imt */
for (c=0;c<3;c++)
fftw_fft2d(t+c*height2*width2,ret+c*ftsize,imt+c*ftsize);
/* Get Conditioning Points */
get_conditioning_points(mask,cond,condw);
/* Compute ADSN realization (the texton FFT must be precomputed!)
[Algorithm 2], [Equation (2)] */
adsn(NULL,z);
/* Right-hand side of kriging system */
setValue(rhs,0.,width2*height2*3);
for (c=0;c<3;c++)
for (i=0;i<height2;i++)
for (j=0;j<width2;j++)
if (cond[Ic(i,j,0)]==1)
rhs[Ic(i,j,c)] = u[Ic(i,j,c)] - meanu[c] - z[Ic(i,j,c)];
/* Conjugate gradient descent [Algorithm 5] */
cgd(cond,rhs,v,ep,imax);
/* Final steps:
apply covariance operator, add mean value,
and reimpose initial known values outside the mask. */
convcov(v,v,NULL);
for (c=0;c<3;c++)
for (i=0;i<height2;i++)
for (j=0;j<width2;j++)
if (mask[Ic(i,j,0)]<0.5)
v[Ic(i,j,c)] += meanu[c] + z[Ic(i,j,c)];
else
v[Ic(i,j,c)] = u[Ic(i,j,c)];
/* Cleanup */
free(u); free(z); free(meanu); free(t);
free(ret); free(imt); free(re); free(im); free(retmp); free(imtmp);
free(rhs); free(cond);
// Terminate FFT structures
term_fftw();
return v; // WARNING! OUTPUT must be of size width2*height2*3
}
int main(int argc, char **argv){
if(argc > 1 && (!strcmp(argv[1],"--help") || !strcmp(argv[1],"-h"))){
//printHelp();
return EXIT_FAILURE;
}
if(argc < 4){
printf("Too few arguments.\n\n");
printf("Syntax: gti [image] [mask] [output] [ep] [niter] [w]\n");
// printf(" try --help for more\n");
return EXIT_FAILURE;
}
numth = MIN(8,omp_get_max_threads());
//atoi(getenv("OMP_NUM_THREADS"));
float ep = 1e-3;
int imax = 100;
int condw = 3;
if (argc>4)
ep = (float) atof(argv[4]);
if (argc>5)
imax = (int) atoi(argv[5]);
if (argc>6)
condw = (int) atoi(argv[6]);
float *im = NULL;
size_t swidth, sheight;
im = io_png_read_f32_rgb(argv[1], &swidth, &sheight);
if(im == NULL){
printf("Are you sure %s is a png image?\n", argv[1]);
return 0;
}
/* The following global variables are set once and for all
and contain the dimensions of the image domain,
the dimensions of the duplicated image domain (that serves
for convolution), and the size of Fourier transforms. */
width = (int) swidth;
height = (int) sheight;
width2 = width*2;
height2 = height*2;
ftsize = (width+1)*height2;
float *m,*mask;
m = io_png_read_f32_rgb(argv[2], &swidth, &sheight);
if(m == NULL){
printf("Are you sure %s is a png image?\n", argv[1]);
free(im);
return EXIT_FAILURE;
}
if(width != (int) swidth || height != (int) sheight){
printf("Dimensions of %s and %s don't agree!\n",argv[1],argv[2]);
free(im);
free(m);
return EXIT_FAILURE;
}
mask = (float*) xmalloc(sizeof(float)*width2*height2*3);
setMask(mask,m);
int i,j,c;
for(i = 0; i < width*height*3; i++)
im[i] /= 255;
float *bigres = gausstexinpaint(im,mask,ep,imax,condw);
float *res;
res = (float*) xmalloc(sizeof(float)*width*height*3);
// crop the result on the initial image domain
for (i=0;i<height;i++)
for (j=0;j<width;j++)
for (c=0;c<3;c++)
res[c*height*width + i*width + j] = 255*bigres[Ic(i,j,c)];
io_png_write_f32(argv[3], res, swidth, sheight, 3);
free(im);
free(m); free(mask);
free(res);
return EXIT_SUCCESS;
}
/**
* @mainpage Gaussian Texture Inpainting
* @section intro Introduction
*
* This code was written for IPOL and is based on the paper
* "Texture Inpainting Using Efficient Gaussian Conditional Simulation"
* (Bruno Galerne, Arthur Leclaire), SIAM Journal on Imaging Sciences, 2017.
* The algorithm is thoroughly discussed in the IPOL companion paper
* "An Algorithm for Gaussian Texture Inpainting"
* (Bruno Galerne, Arthur Leclaire)
*
* The program performs conditional simulation of a Gaussian texture model on the mask, knowing the values of the texture on a border of the mask.
* The Gaussian texture model is estimated on the mask complement.
* The linear system involved in the Gaussian conditional simulation is solved with a conjugate gradient descent.
*
* For more information see @ref usage "usage".
*
* @section install Installation
* Run
* @code make clean; make @endcode
* to compile the code. This requires the following libraries:
* @li fftw3
* @li libpng
* @li openmp
* If OpenMP library is not available, then you should edit the makefile
* and remove libraries containing "omp" or "openmp" before compilation.
*
* the latter can be commented out in the Makefile.
*
* @section usage Usage
* To print out all the options and explanations use
* @code gti --help @endcode
* The program takes inputs in the following way:
* @code gti [image] [mask] [output] [options] @endcode
*
* @li @code [image] @endcode PNG image to be inpainted
* @li @code [mask] @endcode PNG image with white region specifying the inpainting domain
* @li @code [output] @endcode filename of the inpainted PNG
* @li @code [options] @endcode @code ep [default 0.001] niter [default 1000] w [default 3] @endcode
*
* The [image] and [mask] images must have the same dimensions.
*
* The software comes with two test images and the inpainting can be done as follows:
* @code gti discharge_print_128x128.png mask.png out.png @endcode
* or if you want to specify parameters ep=0.01, niter=10000, and w=5 as
* @code gti discharge_print_128x128.png mask.png out.png 0.01 10000 5 @endcode
*
* @section licence Licence
* @subsection gti gti.c
*
* Copyright (c) 2017, Arthur Leclaire <Arthur.Leclaire@cmla.ens-cachan.fr>
*
* This program is free software: you can use, modify and/or
* redistribute it under the terms of the simplified BSD License. You
* should have received a copy of this license along this program. If
* not, see <http://opensource.org/licenses/BSD-3-Clause>.
*
* @subsection io_png io_png.c
*
* Copyright (c) 2010-2011, Nicolas Limare <nicolas.limare@cmla.ens-cachan.fr>
* All rights reserved.
*
* This program is free software: you can redistribute it and/or modify
* it under, at your option, 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, or
* the terms of the simplified BSD license.
*
* You should have received a copy of these licenses along this
* program. If not, see <http://www.gnu.org/licenses/> and
* <http://www.opensource.org/licenses/bsd-license.html>.
*
*/
