efros_freeman  0.1
patch_search.c
Go to the documentation of this file.
1 /*
2  Copyright (c) 2016 Lara Raad <lara.raad@cmla.ens-cachan.fr>,
3  Bruno Galerne <bruno.galerne@parisdescartes.fr>
4 
5  Quilting is free software: you can redistribute it and/or modify
6  it under the terms of the GNU Affero General Public License as
7  published by the Free Software Foundation, either version 3 of the
8  License, or (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Affero General Public License for more details.
14 
15  You should have received a copy of the GNU Affero General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18 
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include <math.h>
22 #include <fftw3.h>
23 #include <assert.h>
24 #include <float.h>
25 
26 #include "patch_search.h"
27 #include "random_number.h"
28 
29 #define MAX_DIST FLT_MAX
30 #define MIN_DIST 1.0
31 #define MAX_INT 10000000
32 #define SMOOTH 1
33 
34 /**
35  * @file distances.c
36  * This file contains all the functions for the patch search procedure that
37  * computes the distance of a patch to all input patches and then randomly
38  * pick on of the patches that are close to the minimal distance.
39  */
40 
41 /**
42  * Convention for image sizes
43  * w_img = nx = width of image : Indexes x or j
44  * h_img = ny = height of image : Indexes y or i
45  *
46  */
47 
48 /**
49  * @brief This function computes the cross-correlation between img1 and img2
50  * using FFT, given the fft of img2 as input.
51  * The result is a long double image having the same size (already allocated).
52  * @param img1, img2, two pointers to float images
53  * @param cross_corr image of the cross-correlation
54  * @param nx, ny image sizes
55  */
56 static void cross_correlation_fft2_known(float *img1, fftwl_complex *fft2,
57  long double *cross_corr, size_t nx, size_t ny)
58 {
59  long double *ldimg;
60  fftwl_complex *fft1;
61  fftwl_plan plan_forward, plan_backward;
62  int i, j, si, sj;
63  size_t size = nx*ny;
64  int fftsize = (int) (ny*((int) (nx/2)+1));
65  long double re, im;
66  int imsize = (int) (nx*ny);
67  long double invimsize = 1./(long double) imsize;
68 
69  /* Memory allocation */
70  ldimg = (long double *) malloc(size*sizeof(long double));
71  #pragma omp critical (make_plan)
72  fft1 = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex)*(fftsize));
73 
74  /* Copy symmetric version of img1 into a long double array
75  and compute FFT */
76  for(i=0; i < (int) ny; i++)
77  {
78  // Define (si, sj) = (-i,-j) modulo (nx,ny)
79  si = (i==0) ? 0 : (ny-i);
80  for(j=0; j < (int) nx; j++)
81  {
82  sj = (j==0) ? 0 : (nx-j);
83  assert(nx*si+sj<size);
84  ldimg[nx*si+sj] = (long double) img1[nx*i+j];
85  }
86  }
87  #pragma omp critical (make_plan)
88  plan_forward = fftwl_plan_dft_r2c_2d(ny,nx,ldimg,fft1,FFTW_ESTIMATE);
89  fftwl_execute(plan_forward);
90 
91  /* Pointwisze multiplication of Fourier transforms */
92  for(i=0; i<fftsize; i++)
93  {
94  re = fft1[i][0];
95  im = fft1[i][1];
96  fft1[i][0] = re*fft2[i][0] - im*fft2[i][1];
97  fft1[i][1] = re*fft2[i][1] + im*fft2[i][0];
98  }
99  /* Backward FFT */
100  #pragma omp critical (make_plan)
101  plan_backward = fftwl_plan_dft_c2r_2d(ny,nx,fft1,cross_corr,FFTW_ESTIMATE);
102  fftwl_execute(plan_backward);
103 
104  /* normalization of the convolution */
105  for(i=0; i<(int) size; i++)
106  {
107  cross_corr[i] *= invimsize;
108  }
109 
110  /* free memory */
111  free(ldimg);
112  #pragma omp critical (make_plan)
113  fftwl_destroy_plan(plan_forward);
114  #pragma omp critical (make_plan)
115  fftwl_destroy_plan(plan_backward);
116  #pragma omp critical (make_plan)
117  fftwl_free(fft1);
118 }
119 
120 /**
121  * @brief Computes the squared norms of each RGB patches of the input according
122  * to the 3 possible overlap regions.The computation is done in applying the
123  * cross-correlation between the extended mask and the sum over the 3 channels
124  * of the squares of the pixel values.
125  * @param src_im input RGB Image
126  * @param patch_normsH long double image for the horizontal overlap
127  * (must be allocated)
128  * @param patch_normsV long double image for the vertical overlap
129  * (must be allocated)
130  * @param patch_normsL long double image for the L-shape overlap
131  * (must be allocated)
132  * @param patch_sz patch size
133  * @param overlap_sz size of the overlap at the edge of the patch
134  * @author Lara Raad, Bruno Galerne
135  */
137  long double *patch_normsH, long double *patch_normsV,
138  long double *patch_normsL, int patch_sz, int overlap_sz)
139 {
140  float *mask_ext;
141  long double *sq_rgb_src_im;
142  fftwl_complex *fft_sq_rgb_src_im;
143  fftwl_plan plan_forward;
144  int i,j;
145  int size = (int) (src_im.h*src_im.w);
146  int fftsize = (int) src_im.h*((int) (src_im.w/2)+1);
147  // allocate memory
148  mask_ext = (float *) calloc(size, sizeof(float));
149  for(i=0; i< size; i++) mask_ext[i] = 0.;
150  sq_rgb_src_im = (long double *) malloc(size*sizeof(long double));
151  #pragma omp critical (make_plan)
152  fft_sq_rgb_src_im =
153  (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex)*(fftsize));
154 
155  // Compute sq_rgb_src_im : each pixel is the sum of the squares of
156  // the pixels of the 3 RGB channels
157  for(i=0; i<size; i++)
158  {
159  sq_rgb_src_im[i] = (long double) src_im.img[i]
160  * (long double) src_im.img[i]
161  + (long double) src_im.img[i+size]
162  * (long double) src_im.img[i+size]
163  + (long double) src_im.img[i+2*size]
164  * (long double) src_im.img[i+2*size];
165  }
166  // Compute the fft of this image
167  #pragma omp critical (make_plan)
168  plan_forward = fftwl_plan_dft_r2c_2d(src_im.h,src_im.w,sq_rgb_src_im,
169  fft_sq_rgb_src_im,FFTW_ESTIMATE);
170  fftwl_execute(plan_forward);
171 
172  // Horizontal overlap
173  for(i=0;i<patch_sz;i++)
174  {
175  for(j=0;j<patch_sz;j++)
176  {
177  if(i<overlap_sz)
178  mask_ext[src_im.w*i+j] = 1.;
179  else
180  mask_ext[src_im.w*i+j] = 0.;
181  }
182  }
183  // compute the norm of all patches with cross-correlation
184  cross_correlation_fft2_known(mask_ext, fft_sq_rgb_src_im,
185  patch_normsH, src_im.w, src_im.h);
186 
187  // Vertical overlap
188  for(i=0;i<patch_sz;i++)
189  {
190  for(j=0;j<patch_sz;j++)
191  {
192  if(j<overlap_sz)
193  mask_ext[src_im.w*i+j] = 1.;
194  else
195  mask_ext[src_im.w*i+j] = 0.;
196  }
197  }
198  // compute the norm of all patches with cross-correlation
199  cross_correlation_fft2_known(mask_ext, fft_sq_rgb_src_im,
200  patch_normsV, src_im.w, src_im.h);
201 
202  // L-shape overlap
203  for(i=0;i<patch_sz;i++)
204  {
205  for(j=0;j<patch_sz;j++)
206  {
207  if(i<overlap_sz)
208  mask_ext[src_im.w*i+j] = 1.;
209  else if (j<overlap_sz)
210  mask_ext[src_im.w*i+j] = 1.;
211  else
212  mask_ext[src_im.w*i+j] = 0.;
213  }
214  }
215  // compute the norm of all patches with cross-correlation
216  cross_correlation_fft2_known(mask_ext, fft_sq_rgb_src_im,
217  patch_normsL, src_im.w, src_im.h);
218 
219  // free memory
220  free(mask_ext);
221  free(sq_rgb_src_im);
222  #pragma omp critical (make_plan)
223  fftwl_destroy_plan(plan_forward);
224  #pragma omp critical (make_plan)
225  fftwl_free(fft_sq_rgb_src_im);
226 }
227 
228 
229 
230 /** @brief This function computes the squared distance between a given patch
231  * and all patches of another image. The squared distance is (implicitly)
232  * weighted by a binary mask on the patch domain.It does not appear because
233  * pixels that are not yet defined (outside the L-shape corner generally) are
234  * set to zero, so that multiplying by the binary mask does nàt change the
235  * patch.
236  * @param src_im The input texture sample
237  * @param out_im The currently synthesized texture image
238  * (unknown pixels must be set to zero)
239  * @param temp The patch from the output image to compare with those in
240  * the texture sample (patch to synthesize)
241  * @param patch_sz The pathc size
242  * @param src_patch_normsH, src_patch_normsV, src_patch_normsL long double
243  * arrays corresponding to the norm of input patches for
244  * the 3 possible overlap regions
245  * @param fft_srcR, fft_srcG,fft_srcB long double Fourier transforms of
246  * the 3 input channels
247  * @param distances Image of distances, WARNING: suppposed to be set to 0
248  * at initialization (result)
249  * @author Lara Raad, Bruno Galerne
250  * @date 2014/06/10
251  */
252 static void compute_distance_to_input_patches(Image src_im, Image out_im,
253  Corner temp, int patch_sz, long double *src_patch_normsH,
254  long double *src_patch_normsV, long double *src_patch_normsL,
255  fftwl_complex *fft_srcR,fftwl_complex *fft_srcG,
256  fftwl_complex *fft_srcB, Image_LD *distances)
257 {
258  long double *cross_corr;
259  float *img1;
260  int i,j,k;
261  float patchij;
262  long double normpatch;
263  long double *src_patch_norms;
264  fftwl_complex *fft2;
265 
266  // allocate memory
267  img1 = (float *) calloc(src_im.h*src_im.w, sizeof(float));
268  cross_corr = (long double *) malloc(src_im.h*src_im.w*sizeof(long double));
269 
270  for(k=0; k<(int)src_im.c; k++)
271  {
272  normpatch = 0.;
273  // extend patch channel times mask into a large image
274  // and compute squared norm of patch
275  for(i=0;i<patch_sz;i++)
276  {
277  for(j=0;j<patch_sz;j++)
278  {
279  patchij = out_im.img[out_im.w*(temp.x+i)+temp.y+j+
280  k*out_im.w*out_im.h];
281  img1[src_im.w*i+j] = patchij;
282  normpatch += (long double) patchij * (long double) patchij;
283  }
284  }
285 
286  // compute cross-correlation
287  // change input fft into the concatenation of the 3 arrays ?
288  if (k==0) fft2 = fft_srcR;
289  else if (k==1) fft2 = fft_srcG;
290  else fft2 = fft_srcB;
291  cross_correlation_fft2_known(img1, fft2, cross_corr,
292  src_im.w, src_im.h);
293 
294  // add to distances
295  for(i=0; i< (int) (src_im.h*src_im.w); i++)
296  {
297  // WARNING: distances is supposed to be allocated and be set to 0
298  // before calling this function
299  distances->img[i] += normpatch - 2.0*cross_corr[i];
300  }
301  }
302  // Add norm of input RGB patches
303  // choose overlap case
304  // L-shape overlap : default case
305  src_patch_norms = src_patch_normsL;
306  if (temp.y == 0)
307  src_patch_norms = src_patch_normsH;// horizontal overlap (first column)
308  if (temp.x == 0)
309  src_patch_norms = src_patch_normsV;// vertical overlap (first row)
310 
311  for(i=0; i<(int) (src_im.h*src_im.w); i++)
312  {
313  distances->img[i] += src_patch_norms[i];
314  }
315 
316  /* free memory */
317  free(img1);
318  free(cross_corr);
319 
320 }
321 
322 
323 /**
324  * This function picks a random patch within all the patches of the source
325  * image that verify dist(p, current) < dist_min*(1+tolerance) with dist_min
326  * the minimum distance between the current patch and the all the patches
327  * from the input image
328  * @author Lara Raad, Bruno Galerne
329  * @param distances Image containing the distance between the current
330  * patch and the patch of source image (distance has
331  * the same size as the input image, but the bottom-right
332  * border must be discarded)
333  * @param h, w size of the domain of admissible patches
334  * @param patch The position of the picked patch (Corner type)
335  * @param tol The threshold used to select patches from the input
336  * @date 2014/06/10
337  */
338 static void pick_patch(Image_LD *distances, int h, int w,
339  Corner *patch, float tol)
340 {
341  Corner *cands;//patches veryfing (dist <= min_dist*(1+tolerance))
342  int rand_num, pos_dist;
343  int count_cands;//amount of patches veryfing dist <= min_dist*(1+tolerance)
344  long double min_dist; // minimal distance between the current patch and
345  // those from the input image
346 
347  // compute the minimal distance
348  min_dist = MAX_DIST;//distances->img[0];
349  for(int i=0; i<h; i++)
350  {
351  for (int j=0; j<w; j++)
352  {
353  pos_dist = i*distances->w + j;
354  assert(distances->img[pos_dist]==distances->img[pos_dist]);
355  if(distances->img[pos_dist] <
356  min_dist && distances->img[pos_dist] > 1.)
357  {
358  min_dist = distances->img[pos_dist];
359  }
360  }
361  }
362 
363  // create the list of patches verifying dist <= min_dist*(1+tolerance)
364  cands = (Corner *) calloc(h*w, sizeof(Corner));
365  count_cands = 0;
366  for(int i=0; i<h; i++)
367  {
368  for (int j=0; j<w; j++)
369  {
370  pos_dist = i*distances->w + j;
371  assert(pos_dist < (int) (distances->w*distances->h));
372  if(distances->img[pos_dist] <= min_dist*(1+tol))
373  {
374  cands[count_cands].x = i;
375  cands[count_cands].y = j;
376  count_cands++;
377  }
378  }
379  }
380  assert(count_cands>0);
381  assert(count_cands<=h*w);
382  // choose a random patch within the list of candidates
383  rand_num = random_number(NULL) % count_cands;
384  assert(rand_num<count_cands);
385  patch->x = cands[rand_num].x;
386  patch->y = cands[rand_num].y;
387 
388  //free memory
389  free(cands);
390 }
391 
392 /**
393  * @brief This function executes the two steps of the patch search procedure:
394  * 1. Compute the distance of the current patch
395  * to all the admissible input patches
396  * 2. Randomly pick one of the patch that is close to the minimal distance
397  * @param src_im The input texture sample
398  * @param out_im The currently synthesized texture image
399  * (unknown pixels must be set to zero)
400  * @param temp The patch from the output image to compare with those
401  * in the texture sample (patch to synthesize)
402  * @param patch_sz The pathc size
403  * @param src_patch_normsH, src_patch_normsV, src_patch_normsL long double
404  * arrays corresponding to the norm of input patches for
405  * the 3 possible overlap regions
406  * @param fft_srcR, fft_srcG,fft_srcB long double Fourier transforms
407  * of the 3 input channels
408  * @param h, w size of the domain of admissible patches
409  * @param patch The position of the picked patch (Corner type)
410  * @param tol The threshold used to select patches from the input image
411  * @date 2014/06/10
412  */
413 void patch_search(Image src_im, Image out_im, Corner temp, int patch_sz,
414  long double *src_patch_normsH, long double *src_patch_normsV,
415  long double *src_patch_normsL, fftwl_complex *fft_srcR,
416  fftwl_complex *fft_srcG, fftwl_complex *fft_srcB, int h, int w,
417  Corner *patch, float tol)
418 {
419  Image_LD distances;
420 
421  /* memory allocation*/
422  initialize_image_LD(&distances, src_im.h, src_im.w, 1);
423 
424  /* compute distance to all input patches */
425  compute_distance_to_input_patches(src_im, out_im, temp, patch_sz,
426  src_patch_normsH, src_patch_normsV, src_patch_normsL,
427  fft_srcR, fft_srcG, fft_srcB, &distances);
428 
429  /* pick a random patch */
430  pick_patch(&distances, h, w, patch, tol);
431 
432  /* free memory */
433  free(distances.img);
434 }