Hugintrunk  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BlendPoisson.h
Go to the documentation of this file.
1 // -*- c-basic-offset: 4 -*-
2 
12 /* This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public
14 * License as published by the Free Software Foundation; either
15 * version 2 of the License, or (at your option) any later version.
16 *
17 * This software is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public
23 * License along with this software. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 */
27 
28 #ifndef POISSON_BLEND_H
29 #define POISSON_BLEND_H
30 
31 #include <iostream>
32 #include <vigra/stdimage.hxx>
33 #include <vigra/convolution.hxx>
34 #include <vigra/stdconvolution.hxx>
35 #include <vigra/basicgeometry.hxx>
36 #include "openmp_vigra.h"
37 
38 #define DEBUG_TIMING
39 
40 namespace vigra_ext
41 {
42 
43 namespace poisson
44 {
45 
46 namespace detail
47 {
48 // helper functions for build gradient map
49 template <class Image, class Mask>
50 inline typename vigra::NumericTraits<typename Image::PixelType>::RealPromote ProcessNeighborPixels(const int x, const int y, const int dx, const int dy, const Image& image, const Mask& mask)
51 {
52  const typename Mask::PixelType m1 = mask[y + dy][x + dx];
53  const typename Mask::PixelType m2 = mask[y - dy][x - dx];
54  if (m1 > 0 && m2 > 0)
55  {
56  return image[y + dy][x + dx] + image[y - dy][x - dx];
57  };
58  // if one of the 2 pixel has no information (outside of image area)
59  // we extrapolate the value from neighbor pixel
60  if (m1 > 0)
61  {
62  return 2 * image[y + dy][x + dx];
63  }
64  else
65  {
66  return 2 * image[y - dy][x - dx];
67  };
68 };
69 
70 template <class Image, class Mask, class SeamMask>
71 inline typename vigra::NumericTraits<typename Image::PixelType>::RealPromote ProcessBorderPixel(const int x, const int y, const int dx, const int dy, const Image& image, const Mask& mask, const SeamMask& seam)
72 {
73  const typename SeamMask::PixelType seam1 = seam[y + dy][x + dx];
74  const typename SeamMask::PixelType seam2 = seam[y - dy][x - dx];
75  const typename Mask::PixelType mask1 = mask[y + dy][x + dx];
76  const typename Mask::PixelType mask2 = mask[y - dy][x - dx];
77  if (seam1 > 0 && seam2 > 0)
78  {
79  if (mask1 > 0 && mask2 > 0)
80  {
81  return image[y + dy][x + dx] + image[y - dy][x - dx];
82  };
83  if (mask1 > 0)
84  {
85  return 2 *image[y + dy][x + dx];
86  }
87  else
88  {
89  return 2 * image[y - dy][x - dx];
90  };
91  };
92  if (seam1 > 0)
93  {
94  if (mask1 > 0)
95  {
96  return 2 * image[y + dy][x + dx];
97  }
98  else
99  {
100  return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
101  };
102  };
103  if (seam2 > 0)
104  {
105  if (mask2 > 0)
106  {
107  return 2 * image[y - dy][x - dx];
108  }
109  else
110  {
111  return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
112  };
113  };
114  return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
115 };
116 
117 template <class Image, class SeamMask>
118 inline typename Image::PixelType GetBorderGradient(const int x, const int y, const int dx, const int dy, const SeamMask& seams, const Image& image1, const vigra::Point2D& offset)
119 {
120  if (seams[y + dy][x + dx] == 1)
121  {
122  return image1[offset.y + y + dy][offset.x + x + dx];
123  };
124  return vigra::NumericTraits<typename Image::PixelType>::zero();
125 };
126 
127 template <class Image, class SeamMask>
128 inline typename vigra::NumericTraits<typename Image::PixelType>::RealPromote GetBorderValues(const int x, const int y, int dx, int dy, const Image& image, const SeamMask& seams)
129 {
130  const typename SeamMask::PixelType s1 = seams[y + dy][x + dx];
131  const typename SeamMask::PixelType s2 = seams[y - dy][x - dx];
132  if (s1 > 1 && s2 > 1)
133  {
134  return image[y + dy][x + dx] + image[y - dy][x - dx];
135  }
136  else
137  {
138  return (2 - std::min<int>(s2, 2))*image[y + dy][x + dx] + (2 - std::min<int>(s1, 2))*image[y - dy][x - dx];
139  };
140 };
141 
142 template <class Image>
143 void RestrictErrorToNextLevel(const Image & in, Image & out)
144 {
145  // smooth input image
146  vigra::Kernel2D<double> filter2D;
147  filter2D.initExplicitly(vigra::Diff2D(-1, -1), vigra::Diff2D(1, 1)) = 0.25, 0.5, 0.25, 0.5, 1, 0.5, 0.25, 0.5, 0.25;
148  Image smoothImage(in.size());
149  vigra::convolveImage(vigra::srcImageRange(in), vigra::destImage(smoothImage), vigra::kernel2d(filter2D));
150  // downsample smoothed image
151  vigra::resizeImageLinearInterpolation(vigra::srcImageRange(smoothImage), vigra::destImageRange(out));
152 }
153 
154 // internal use by FindEdgesForPoisson
156 {
157  vigra::Int8 operator()(float const& vf) const
158  {
159  const int v = static_cast<int>(vf);
160  if (v > 0 && v < 17)
161  {
162  return 1;
163  }
164  else
165  {
166  if (v >= 85)
167  {
168  if (v == 85 || v == 89 || v == 93 || v == 97)
169  {
170  return 3;
171  }
172  else
173  {
174  return 2;
175  }
176  }
177  else
178  {
179  return 0;
180  };
181  };
182  };
183 };
184 
185 // multi-threaded variant for convolution of images, only 4-neighborhood is considerd
186 template <class Image1, class Image2>
187 void SimpleConvolveImage4(const Image1& image1, Image2& image2, const double factor1, const double factor2)
188 {
189  vigra_precondition(image1.size() == image2.size(), "ConvolveImage: Image size does not match");
190  vigra_precondition(image1.width() >= 2 && image1.height() >= 2, "ConvolveImage: Image too small");
191  const int width = image1.width();
192  const int height = image1.height();
193  //special treatment of first line
194  image2[0][0] = factor1*image1[0][0] + factor2*image1[0][1] + factor2*image1[1][0];
195  for (int x = 1; x < width - 1; ++x)
196  {
197  image2[0][x] = factor1*image1[0][x] + factor2*image1[0][x - 1] + factor2*image1[0][x + 1] + factor2*image1[1][x];
198  };
199  image2[0][width - 1] = factor1*image1[0][width - 1] + factor2*image1[0][width - 2] + factor2*image1[1][width - 1];
200 #pragma omp parallel for
201  for (int y = 1; y < height - 1; ++y)
202  {
203  image2[y][0] = factor1*image1[y][0] + factor2*image1[y - 1][0]
204  + factor2*image1[y][1] + factor2*image1[y + 1][0];
205  for (size_t x = 1; x < width - 1; ++x)
206  {
207  image2[y][x] = factor1*image1[y][x] + factor2*image1[y - 1][x]
208  + factor2*image1[y][x - 1] + factor2*image1[y + 1][x] + factor2*image1[y][x + 1];
209  }
210  image2[y][width - 1] = factor1*image1[y][width - 1] + factor2*image1[y - 1][width - 1]
211  + factor2*image1[y][width - 2] + factor2*image1[y + 1][width - 1];
212  };
213  //special treatment of last line
214  image2[height - 1][0] = factor1*image1[height - 1][0] + factor2*image1[height - 1][1] + factor2*image1[height - 2][0];
215  for (size_t x = 1; x < width - 1; ++x)
216  {
217  image2[height - 1][x] = factor1*image1[height - 1][x] + factor2*image1[height - 1][x - 1] + factor2*image1[height - 1][x + 1] + factor2*image1[height - 2][x];
218  };
219  image2[height - 1][width - 1] = factor1*image1[height - 1][width - 1] + factor2*image1[height - 1][width - 2] + factor2*image1[height - 2][width - 1];
220 };
221 
234 template <class Image>
235 vigra::Int8Image FindEdgesForPoisson(const Image& input)
236 {
237  vigra::Int8Image output(input.size());
238  SimpleConvolveImage4(input, output, 21, -1);
240  return output;
241 };
242 
243 template <class ComponentType>
244 double GetRealValue(const ComponentType& val) { return val; }
245 
246 template <class ComponentType>
247 double GetRealValue(const vigra::RGBValue<ComponentType>& val) { return val.magnitude(); }
248 
249 template <class Image, class SeamMask>
250 void SOR(Image& target, const Image& gradient, const SeamMask& seams, const float omega, const float errorThreshold, const int maxIter, const bool doWrap)
251 {
252  typedef typename Image::PixelType TargetPixelType;
253  const int width = target.width();
254  const int height = target.height();
255 
256  // changes in last iteration
257  double oldError = 0;
258  for (int j = 0; j < maxIter; j++)
259  {
260  // changes in current iteration
261  double error = 0;
262  // special treatment of first line
263  if (seams[0][0] > 1)
264  {
265  if (doWrap)
266  {
267  const TargetPixelType delta = omega * ((gradient[0][0] + target[0][1] + 2 * target[1][0] + target[0][width - 1]) / 4.0f - target[0][0]);
268  error += detail::GetRealValue(delta*delta);
269  target[0][0] += delta;
270  }
271  else
272  {
273  const TargetPixelType delta = omega * ((gradient[0][0] + 2 * target[0][1] + 2 * target[1][0]) / 4.0f - target[0][0]);
274  error += detail::GetRealValue(delta*delta);
275  target[0][0] += delta;
276  };
277  };
278  for (int x = 1; x < width - 1; ++x)
279  {
280  if (seams[0][x] > 1)
281  {
282  const TargetPixelType delta = omega * ((gradient[0][x] + detail::GetBorderValues(x, 0, 1, 0, target, seams) + 2 * target[1][x]) / 4.0f - target[0][x]);
283  error += detail::GetRealValue(delta*delta);
284  target[0][x] += delta;
285  };
286  };
287  if (seams[0][width - 1] > 1)
288  {
289  if (doWrap)
290  {
291  const TargetPixelType delta = omega * ((gradient[0][width - 1] + target[0][width - 2] + 2 * target[1][width - 1] + target[0][0]) / 4.0f - target[0][width - 1]);
292  error += detail::GetRealValue(delta*delta);
293  target[0][width - 1] += delta;
294  }
295  else
296  {
297  const TargetPixelType delta = omega * ((gradient[0][width - 1] + 2 * target[0][width - 2] + 2 * target[1][width - 1]) / 4.0f - target[0][width - 1]);
298  error += detail::GetRealValue(delta*delta);
299  target[0][width - 1] += delta;
300  };
301  };
302 #pragma omp parallel for reduction(+: error) schedule(dynamic, 100)
303  for (int y = 1; y < height - 1; ++y)
304  {
305  if (seams[y][0] > 1)
306  {
307  if (doWrap)
308  {
309  const TargetPixelType delta = omega * ((gradient[y][0] + detail::GetBorderValues(0, y, 0, 1, target, seams)
310  + target[y][1] + target[y][width - 1]) / 4.0f - target[y][0]);
311  error += detail::GetRealValue(delta*delta);
312  target[y][0] += delta;
313  }
314  else
315  {
316  const TargetPixelType delta = omega * ((gradient[y][0] + detail::GetBorderValues(0, y, 0, 1, target, seams) + 2 * target[y][1]) / 4.0f - target[y][0]);
317  error += detail::GetRealValue(delta*delta);
318  target[y][0] += delta;
319  };
320  };
321  for (int x = 1; x < width - 1; ++x)
322  {
323  const typename SeamMask::value_type maskValue = seams[y][x];
324  if (maskValue > 1)
325  {
326  if (maskValue == 2)
327  {
328  // border pixel
329  const TargetPixelType sum = detail::GetBorderValues(x, y, 1, 0, target, seams) + detail::GetBorderValues(x, y, 0, 1, target, seams);
330  const TargetPixelType delta = omega * ((gradient[y][x] + sum) / 4.0f - target[y][x]);
331  error += detail::GetRealValue(delta*delta);
332  target[y][x] += delta;
333  }
334  else
335  {
336  const TargetPixelType sum = target[y + 1][x] + target[y][x + 1] + target[y - 1][x] + target[y][x - 1];
337  const TargetPixelType delta = omega * ((gradient[y][x] + sum) / 4.0f - target[y][x]);
338  error += detail::GetRealValue(delta*delta);
339  target[y][x] += delta;
340  };
341  };
342  };
343  if (seams[y][width - 1] > 1)
344  {
345  if (doWrap)
346  {
347  const TargetPixelType delta = omega * ((gradient[y][width - 1] + detail::GetBorderValues(width - 1, y, 0, 1, target, seams) + target[y][width - 2] + target[y][0]) / 4.0f - target[y][width - 1]);
348  error += detail::GetRealValue(delta*delta);
349  target[y][width - 1] += delta;
350  }
351  else
352  {
353  const TargetPixelType delta = omega * ((gradient[y][width - 1] + detail::GetBorderValues(width - 1, y, 0, 1, target, seams) + 2 * target[y][width - 2]) / 4.0f - target[y][width - 1]);
354  error += detail::GetRealValue(delta*delta);
355  target[y][width - 1] += delta;
356  };
357  };
358  };
359  // last line
360  if (seams[height - 1][0] > 1)
361  {
362  if (doWrap)
363  {
364  const TargetPixelType delta = omega * ((gradient[height - 1][0] + 2 * target[height - 2][0] + target[height - 1][1] + target[height - 1][width - 1]) / 4.0f - target[height - 1][0]);
365  error += detail::GetRealValue(delta*delta);
366  target[height - 1][0] += delta;
367  }
368  else
369  {
370  const TargetPixelType delta = omega * ((gradient[height - 1][0] + 2 * target[height - 2][0] + 2 * target[height - 1][1]) / 4.0f - target[height - 1][0]);
371  error += detail::GetRealValue(delta*delta);
372  target[height - 1][0] += delta;
373  };
374  };
375  for (int x = 1; x < width - 1; ++x)
376  {
377  if (seams[height - 1][x] > 1)
378  {
379  const TargetPixelType delta = omega * ((gradient[height - 1][x] + detail::GetBorderValues(x, height - 1, 1, 0, target, seams) + 2 * target[height - 2][x]) / 4.0f - target[height - 1][x]);
380  error += detail::GetRealValue(delta*delta);
381  target[height - 1][x] += delta;
382  };
383  };
384  if (seams[height - 1][width - 1] > 1)
385  {
386  if (doWrap)
387  {
388  const TargetPixelType delta = omega * ((gradient[height - 1][width - 1] + 2 * target[height - 2][width - 1] + target[height - 1][width - 2] + target[height - 1][0]) / 4.0f - target[height - 1][width - 1]);
389  error += detail::GetRealValue(delta*delta);
390  target[height - 1][width - 1] += delta;
391  }
392  else
393  {
394  const TargetPixelType delta = omega * ((gradient[height - 1][width - 1] + 2 * target[height - 2][width - 1] + 2 * target[height - 1][width - 2]) / 4.0f - target[height - 1][width - 1]);
395  error += detail::GetRealValue(delta*delta);
396  target[height - 1][width - 1] += delta;
397  };
398  };
399 
400  if (oldError > 0 && log(oldError / error) / log(10.0) < errorThreshold)
401  {
402  break;
403  }
404  oldError = error;
405  }
406 }
407 
408 template <class Image, class SeamMask>
409 void CalcResidualError(Image& error, const Image& target, const Image& gradient, const SeamMask& seam, const bool doWrap)
410 {
411  typedef typename Image::PixelType ImagePixelType;
412  const int width = target.width();
413  const int height = target.height();
414 
415  if (seam[0][0] > 1)
416  {
417  if (doWrap)
418  {
419  const ImagePixelType sum = 2 * target[1][0] + target[0][1] + target[0][width - 1];
420  error[0][0] = (4 * target[0][0] - sum - gradient[0][0]);
421  }
422  else
423  {
424  const ImagePixelType sum = 2 * target[1][0] + 2 * target[0][1];
425  error[0][0] = (4 * target[0][0] - sum - gradient[0][0]);
426  };
427  };
428  for (int x = 1; x < width - 1; ++x)
429  {
430  if (seam[0][x]>1)
431  {
432  const ImagePixelType sum = 2 * target[1][x] + detail::GetBorderValues(x, 0, 1, 0, target, seam);
433  error[0][x] = (4 * target[0][x] - sum - gradient[0][x]);
434  };
435  };
436  if (seam[0][width - 1] > 1)
437  {
438  if (doWrap)
439  {
440  const ImagePixelType sum = 2 * target[1][width - 1] + target[0][width - 2] + target[0][0];
441  error[0][width - 1] = (4 * target[0][width - 1] - sum - gradient[0][width - 1]);
442  }
443  else
444  {
445  const ImagePixelType sum = 2 * target[1][width - 1] + 2 * target[0][width - 2];
446  error[0][width - 1] = (4 * target[0][width - 1] - sum - gradient[0][width - 1]);
447  }
448  };
449 #pragma omp parallel for schedule(dynamic, 100)
450  for (int y = 1; y < height - 1; ++y)
451  {
452  if (seam[y][0] > 1)
453  {
454  if (doWrap)
455  {
456  const ImagePixelType sum = detail::GetBorderValues(0, y, 0, 1, target, seam) + target[y][1] + target[y][width - 1];
457  error[y][0] = (4 * target[y][0] - sum - gradient[y][0]);
458  }
459  else
460  {
461  const ImagePixelType sum = detail::GetBorderValues(0, y, 0, 1, target, seam) + 2 * target[y][1];
462  error[y][0] = (4 * target[y][0] - sum - gradient[y][0]);
463  };
464  }
465  for (int x = 1; x < width - 1; ++x)
466  {
467  const typename SeamMask::value_type maskValue = seam[y][x];
468  if (maskValue > 1)
469  {
470  if (maskValue == 2)
471  {
472  // border pixel
473  const ImagePixelType sum = detail::GetBorderValues(x, y, 1, 0, target, seam) + detail::GetBorderValues(x, y, 0, 1, target, seam);
474  error[y][x] = (4 * target[y][x] - sum - gradient[y][x]);
475  }
476  else
477  {
478  const ImagePixelType sum = target[y + 1][x] + target[y][x + 1] + target[y - 1][x] + target[y][x - 1];
479  error[y][x] = (4 * target[y][x] - sum - gradient[y][x]);
480  };
481  };
482  };
483  if (seam[y][width - 1] > 1)
484  {
485  if (doWrap)
486  {
487  const ImagePixelType sum = detail::GetBorderValues(width - 1, y, 0, 1, target, seam) + target[y][width - 2] + target[y][0];
488  error[y][width - 1] = (4 * target[y][width - 1] - sum - gradient[y][width - 1]);
489  }
490  else
491  {
492  const ImagePixelType sum = detail::GetBorderValues(width - 1, y, 0, 1, target, seam) + 2 * target[y][width - 2];
493  error[y][width - 1] = (4 * target[y][width - 1] - sum - gradient[y][width - 1]);
494  };
495  };
496  };
497  // last line
498  if (seam[height - 1][0] > 1)
499  {
500  if (doWrap)
501  {
502  const ImagePixelType sum = 2 * target[height - 2][0] + target[height - 1][width - 1] + target[height - 1][1];
503  error[height - 1][0] = (4 * target[height - 1][0] - sum - gradient[height - 1][0]);
504  }
505  else
506  {
507  const ImagePixelType sum = 2 * target[height - 2][0] + 2 * target[height - 1][1];
508  error[height - 1][0] = (4 * target[height - 1][0] - sum - gradient[height - 1][0]);
509  };
510  };
511  for (int x = 1; x < width - 1; ++x)
512  {
513  if (seam[height - 1][x]>1)
514  {
515  const ImagePixelType sum = 2 * target[height-2][x] + detail::GetBorderValues(x, height - 1, 1, 0, target, seam);
516  error[height - 1][x] = (4 * target[height - 1][x] - sum - gradient[height - 1][x]);
517  };
518  };
519  if (seam[height - 1][width - 1] > 1)
520  {
521  if (doWrap)
522  {
523  const ImagePixelType sum = 2 * target[height - 2][width - 1] + target[height - 1][width - 2] + target[height - 1][0];
524  error[height - 1][width - 1] = (4 * target[height - 1][width - 1] - sum - gradient[height - 1][width - 1]);
525  }
526  else
527  {
528  const ImagePixelType sum = 2 * target[height - 2][width - 1] + 2 * target[height - 1][width - 2];
529  error[height - 1][width - 1] = (4 * target[height - 1][width - 1] - sum - gradient[height - 1][width - 1]);
530  };
531  };
532 
533 }
534 
535 } // namespace detail
536 template <class PixelType>
538 {
539 public:
540  explicit MaskGreaterAccessor(PixelType val) : m_value(val) {};
541 
542  template <class ITERATOR>
543  PixelType operator()(ITERATOR const & i) const {
544  if (PixelType(*i) >= m_value)
545  {
547  }
548  else
549  {
550  return vigra::NumericTraits<PixelType>::zero();
551  };
552  }
553 
554  template <class ITERATOR, class DIFFERENCE>
555  PixelType operator()(ITERATOR const & i, DIFFERENCE d) const
556  {
557  if (PixelType(i[d]) >= m_value)
558  {
560  }
561  else
562  {
563  return vigra::NumericTraits<PixelType>::zero();
564  };
565  }
566  PixelType m_value;
567 };
568 
569 template <class PixelType>
571 {
572 public:
573  explicit MaskSmallerAccessor(PixelType val) : m_value(val) {};
574 
575  template <class ITERATOR>
576  PixelType operator()(ITERATOR const & i) const {
577  if (PixelType(*i) < m_value)
578  {
580  }
581  else
582  {
583  return vigra::NumericTraits<PixelType>::zero();
584  };
585  }
586 
587  template <class ITERATOR, class DIFFERENCE>
588  PixelType operator()(ITERATOR const & i, DIFFERENCE d) const
589  {
590  if (PixelType(i[d]) < m_value)
591  {
593  }
594  else
595  {
596  return vigra::NumericTraits<PixelType>::zero();
597  };
598  }
599  PixelType m_value;
600 };
601 
602 template <class Image, class PyramidImage>
603 void BuildSeamPyramid(const Image& input, vigra::ImagePyramid<PyramidImage>& seams, const int minLength)
604 {
605  const int nlevels = (int)ceil(log(std::min(input.width(), input.height()) / (double)minLength) / log(2.0));
606  seams.resize(0, nlevels, input.size());
607  Image scaledImage = input;
608  seams[0] = detail::FindEdgesForPoisson(input);
609  for (size_t i = 1; i <= seams.highestLevel(); ++i)
610  {
611  Image smaller((scaledImage.width() + 1) / 2, (scaledImage.height() + 1) / 2);
612  vigra::resizeImageNoInterpolation(vigra::srcImageRange(scaledImage), vigra::destImageRange(smaller));
613  seams[i] = detail::FindEdgesForPoisson(smaller);
614  scaledImage = smaller;
615  };
616 }
617 
618 template<class Image, class Mask, class SeamMask, class GradientType>
619 void BuildGradientMap(const Image& image1, const Image& image2, const Mask& mask2, const SeamMask& seam, GradientType& gradient, const vigra::Point2D& offset, const bool doWrap)
620 {
621  typedef typename GradientType::PixelType GradientPixelType;
622  const int width = image2.width();
623  const int height = image2.height();
624  //special treatment of first line
625  if (seam[0][0] == 2)
626  {
627  if (doWrap)
628  {
629  GradientPixelType value = 4 * image2[0][0] - image2[0][1] - 2 * image2[1][0] - image2[0][width - 1];
630  // copy values for Dirichlet boundary condition
631  value += detail::GetBorderGradient(0, 0, 1, 0, seam, image1, offset);
632  value += detail::GetBorderGradient(0, 0, 0, 1, seam, image1, offset);
633  value += detail::GetBorderGradient(width - 2, 0, 1, 0, seam, image1, offset);
634  gradient[0][0] = value;
635  }
636  else
637  {
638  GradientPixelType value = 4 * image2[0][0] - 2 * image2[0][1] - 2 * image2[1][0];
639  // copy values for Dirichlet boundary condition
640  value += detail::GetBorderGradient(0, 0, 1, 0, seam, image1, offset);
641  value += detail::GetBorderGradient(0, 0, 0, 1, seam, image1, offset);
642  gradient[0][0] = value;
643  };
644  };
645  for (int x = 1; x < width - 1; ++x)
646  {
647  if (seam[0][x] == 2)
648  {
649  GradientPixelType value = 4 * image2[0][x] - detail::ProcessBorderPixel(x, 0, 1, 0, image2, mask2, seam) - 2 * image2[1][x];
650  value += detail::GetBorderGradient(x, 0, -1, 0, seam, image1, offset);
651  value += detail::GetBorderGradient(x, 0, 1, 0, seam, image1, offset);
652  value += detail::GetBorderGradient(x, 0, 0, 1, seam, image1, offset);
653  gradient[0][x] = value;
654  };
655  };
656  if (seam[0][width - 1] == 2)
657  {
658  if (doWrap)
659  {
660  GradientPixelType value = 4 * image2[0][width - 1] - image2[0][width - 2] - 2 * image2[1][width - 1] - image2[0][0];
661  value += detail::GetBorderGradient(width - 1, 0, -1, 0, seam, image1, offset);
662  value += detail::GetBorderGradient(width - 1, 0, 0, 1, seam, image1, offset);
663  value += detail::GetBorderGradient(1, 0, -1, 0, seam, image1, offset);
664  gradient[0][width - 1] = value;
665  }
666  else
667  {
668  GradientPixelType value = 4 * image2[0][width - 1] - 2 * image2[0][width - 2] - 2 * image2[1][width - 1];
669  value += detail::GetBorderGradient(width - 1, 0, -1, 0, seam, image1, offset);
670  value += detail::GetBorderGradient(width - 1, 0, 0, 1, seam, image1, offset);
671  gradient[0][width - 1] = value;
672  };
673  };
674 #pragma omp parallel for
675  for (int y = 1; y < height - 1; ++y)
676  {
677  if (seam[y][0] == 2)
678  {
679  if (doWrap)
680  {
681  GradientPixelType value = 4 * image2[y][0] - detail::ProcessBorderPixel(0, y, 0, 1, image2, mask2, seam) - image2[y][1] - image2[y][width - 1];
682  value += detail::GetBorderGradient(0, y, 1, 0, seam, image1, offset);
683  value += detail::GetBorderGradient(0, y, 0, 1, seam, image1, offset);
684  value += detail::GetBorderGradient(0, y, 0, -1, seam, image1, offset);
685  value += detail::GetBorderGradient(width - 2, y, 1, 0, seam, image1, offset);
686  gradient[y][0] = value;
687  }
688  else
689  {
690  GradientPixelType value = 4 * image2[y][0] - detail::ProcessBorderPixel(0, y, 0, 1, image2, mask2, seam) - 2 * image2[y][1];
691  value += detail::GetBorderGradient(0, y, 1, 0, seam, image1, offset);
692  value += detail::GetBorderGradient(0, y, 0, 1, seam, image1, offset);
693  value += detail::GetBorderGradient(0, y, 0, -1, seam, image1, offset);
694  gradient[y][0] = value;
695  };
696  };
697  for (int x = 1; x < width - 1; ++x)
698  {
699  const typename SeamMask::PixelType seamVal = seam[y][x];
700  if (seamVal>1)
701  {
702  GradientPixelType value = 4.0 * image2[y][x];
703  if (seamVal == 3)
704  {
705  value -= detail::ProcessNeighborPixels(x, y, 1, 0, image2, mask2);
706  value -= detail::ProcessNeighborPixels(x, y, 0, 1, image2, mask2);
707  }
708  else
709  {
710  // seamVal==2, border pixel
711  value -= detail::ProcessBorderPixel(x, y, 1, 0, image2, mask2, seam);
712  value -= detail::ProcessBorderPixel(x, y, 0, 1, image2, mask2, seam);
713  };
714  // copy values for Dirichlet boundary condition
715  value += detail::GetBorderGradient(x, y, 1, 0, seam, image1, offset);
716  value += detail::GetBorderGradient(x, y, 0, 1, seam, image1, offset);
717  value += detail::GetBorderGradient(x, y, -1, 0, seam, image1, offset);
718  value += detail::GetBorderGradient(x, y, 0, -1, seam, image1, offset);
719  gradient[y][x] = value;
720  }
721  };
722  if (seam[y][width - 1] == 2)
723  {
724  if (doWrap)
725  {
726  GradientPixelType value = 4.0 * image2[y][width - 1] - detail::ProcessBorderPixel(width - 1, y, 0, 1, image2, mask2, seam) - image2[y][width - 2] - image2[y][0];
727  value += detail::GetBorderGradient(width - 1, y, -1, 0, seam, image1, offset);
728  value += detail::GetBorderGradient(width - 1, y, 0, 1, seam, image1, offset);
729  value += detail::GetBorderGradient(width - 1, y, 0, -1, seam, image1, offset);
730  value += detail::GetBorderGradient(0, y, -1, 0, seam, image1, offset);
731  gradient[y][width - 1] = value;
732  }
733  else
734  {
735  GradientPixelType value = 4.0 * image2[y][width - 1] - detail::ProcessBorderPixel(width - 1, y, 0, 1, image2, mask2, seam) - 2 * image2[y][width - 2];
736  value += detail::GetBorderGradient(width - 1, y, -1, 0, seam, image1, offset);
737  value += detail::GetBorderGradient(width - 1, y, 0, 1, seam, image1, offset);
738  value += detail::GetBorderGradient(width - 1, y, 0, -1, seam, image1, offset);
739  gradient[y][width - 1] = value;
740  };
741  };
742  };
743  //special treatment of last line
744  if (seam[height - 1][0] == 2)
745  {
746  if (doWrap)
747  {
748  GradientPixelType value = 4 * image2[height - 1][0] - image2[height - 1][1] - 2 * image2[height - 2][0] - image2[height - 1][width - 1];
749  value += detail::GetBorderGradient(0, height - 1, 1, 0, seam, image1, offset);
750  value += detail::GetBorderGradient(0, height - 1, 0, -1, seam, image1, offset);
751  value += detail::GetBorderGradient(width - 2, height - 1, 1, 0, seam, image1, offset);
752  gradient[height - 1][0] = value;
753  }
754  else
755  {
756  GradientPixelType value = 4 * image2[height - 1][0] - 2 * image2[height - 1][1] - 2 * image2[height - 2][0];
757  value += detail::GetBorderGradient(0, height - 1, 1, 0, seam, image1, offset);
758  value += detail::GetBorderGradient(0, height - 1, 0, -1, seam, image1, offset);
759  gradient[height - 1][0] = value;
760  };
761  };
762  for (size_t x = 1; x < width - 1; ++x)
763  {
764  if (seam[height - 1][x] == 2)
765  {
766  GradientPixelType value = 4 * image2[height - 1][x] - detail::ProcessBorderPixel(x, height - 1, 1, 0, image2, mask2, seam) - 2 * image2[height - 2][x];
767  value += detail::GetBorderGradient(x, height - 1, -1, 0, seam, image1, offset);
768  value += detail::GetBorderGradient(x, height - 1, 1, 0, seam, image1, offset);
769  value += detail::GetBorderGradient(x, height - 1, 0, -1, seam, image1, offset);
770  gradient[height - 1][x] = value;
771  };
772  };
773  if (seam[height - 1][width - 1] == 2)
774  {
775  if (doWrap)
776  {
777  GradientPixelType value = 4 * image2[height - 1][width - 1] - image2[height - 1][width - 2] - 2 * image2[height - 2][width - 1] - image2[height - 1][0];
778  value += detail::GetBorderGradient(width - 1, height - 1, -1, 0, seam, image1, offset);
779  value += detail::GetBorderGradient(width - 1, height - 1, 0, -1, seam, image1, offset);
780  value += detail::GetBorderGradient(1, height - 1, -1, 0, seam, image1, offset);
781  gradient[height - 1][width - 1] = value;
782  }
783  else
784  {
785  GradientPixelType value = 4 * image2[height - 1][width - 1] - 2 * image2[height - 1][width - 2] - 2 * image2[height - 2][width - 1];
786  value += detail::GetBorderGradient(width - 1, height - 1, -1, 0, seam, image1, offset);
787  value += detail::GetBorderGradient(width - 1, height - 1, 0, -1, seam, image1, offset);
788  gradient[height - 1][width - 1] = value;
789  };
790  };
791 };
792 
793 template <class Image, class SeamMask>
794 void Multigrid(Image& out, const Image& gradient, const vigra::ImagePyramid<SeamMask>& seamMaskPyramid, int minLen, const float errorThreshold, const int maxIter, const bool doWrap)
795 {
796  const int width = out.width();
797  const int height = out.height();
798 
799  if (width < minLen || height < minLen)
800  {
801  return;
802  }
803  Image err(width, height);
804  Image err2((width + 1) / 2, (height + 1) / 2);
805  Image out2(err2.size());
806  int maskIndex = -1;
807  for (int i = 0; i <= seamMaskPyramid.highestLevel(); ++i)
808  {
809  if (out.size() == seamMaskPyramid[i].size())
810  {
811  maskIndex = i;
812  break;
813  };
814  }
815  if (maskIndex == -1)
816  {
817  std::cout << "ERROR: No suitable mask, this should not happen." << std::endl
818  << "searching " << out.size() << ", finest " << seamMaskPyramid[seamMaskPyramid.highestLevel()].size() << std::endl;
819  return;
820  };
821  // pre-smoothing
822  const float omega = 1.6f; // relaxation parameter: 0 < omega < 2
823  detail::SOR(out, gradient, seamMaskPyramid[maskIndex], omega, errorThreshold, maxIter, doWrap);
824  detail::CalcResidualError(err, out, gradient, seamMaskPyramid[maskIndex], doWrap); // Fehler berechnen
826  Multigrid(out2, err2, seamMaskPyramid, minLen, errorThreshold, maxIter, doWrap);
827  // now do a W cycle calculation, another run of multigrid follows
828  Multigrid(out2, err2, seamMaskPyramid, minLen, errorThreshold, maxIter, doWrap);
829  vigra::resizeImageLinearInterpolation(srcImageRange(out2), destImageRange(err));
832  vigra::destImage(out),
833  vigra::functor::Arg1() - vigra::functor::Arg2());
834  // post smoothing
835  detail::SOR(out, gradient, seamMaskPyramid[maskIndex], omega, errorThreshold, maxIter, doWrap);
836  return;
837 }
838 
839 } // namespace poisson
840 } // namespace vigra_ext
841 
842 #endif //POISSON_BLEND_H
void RestrictErrorToNextLevel(const Image &in, Image &out)
Definition: BlendPoisson.h:143
PixelType operator()(ITERATOR const &i) const
Definition: BlendPoisson.h:576
Image::PixelType GetBorderGradient(const int x, const int y, const int dx, const int dy, const SeamMask &seams, const Image &image1, const vigra::Point2D &offset)
Definition: BlendPoisson.h:118
double GetRealValue(const ComponentType &val)
Definition: BlendPoisson.h:244
PixelType operator()(ITERATOR const &i, DIFFERENCE d) const
Definition: BlendPoisson.h:555
vigra::NumericTraits< typename Image::PixelType >::RealPromote GetBorderValues(const int x, const int y, int dx, int dy, const Image &image, const SeamMask &seams)
Definition: BlendPoisson.h:128
vigra::pair< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImage(const ROIImage< Image, Mask > &img)
Definition: ROIImage.h:300
void CalcResidualError(Image &error, const Image &target, const Image &gradient, const SeamMask &seam, const bool doWrap)
Definition: BlendPoisson.h:409
void BuildGradientMap(const Image &image1, const Image &image2, const Mask &mask2, const SeamMask &seam, GradientType &gradient, const vigra::Point2D &offset, const bool doWrap)
Definition: BlendPoisson.h:619
PixelType operator()(ITERATOR const &i) const
Definition: BlendPoisson.h:543
void SOR(Image &target, const Image &gradient, const SeamMask &seams, const float omega, const float errorThreshold, const int maxIter, const bool doWrap)
Definition: BlendPoisson.h:250
vigra::Int8Image FindEdgesForPoisson(const Image &input)
mark edges in input image for poisson blending * input: expected an image with following meanings lab...
Definition: BlendPoisson.h:235
void combineTwoImagesIf(SrcImageIterator1 src1_upperleft, SrcImageIterator1 src1_lowerright, SrcAccessor1 src1_acc, SrcImageIterator2 src2_upperleft, SrcAccessor2 src2_acc, MaskImageIterator mask_upperleft, MaskAccessor mask_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc, const Functor &func)
Definition: openmp_vigra.h:268
vigra::NumericTraits< typename Image::PixelType >::RealPromote ProcessNeighborPixels(const int x, const int y, const int dx, const int dy, const Image &image, const Mask &mask)
Definition: BlendPoisson.h:50
vigra::pair< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImage(ROIImage< Image, Alpha > &img)
Definition: ROIImage.h:324
vigra::triple< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImageRange(const ROIImage< Image, Mask > &img)
helper function for ROIImages
Definition: ROIImage.h:287
vigra::RGBValue< T, RIDX, GIDX, BIDX > log(vigra::RGBValue< T, RIDX, GIDX, BIDX > const &v)
component-wise logarithm
Definition: utils.h:209
PixelType operator()(ITERATOR const &i, DIFFERENCE d) const
Definition: BlendPoisson.h:588
static T max(T x, T y)
Definition: svm.cpp:65
vigra::Int8 operator()(float const &vf) const
Definition: BlendPoisson.h:157
void SimpleConvolveImage4(const Image1 &image1, Image2 &image2, const double factor1, const double factor2)
Definition: BlendPoisson.h:187
vigra::triple< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImageRange(ROIImage< Image, Alpha > &img)
Definition: ROIImage.h:312
void transformImage(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc, const Functor &func)
Definition: openmp_vigra.h:330
void BuildSeamPyramid(const Image &input, vigra::ImagePyramid< PyramidImage > &seams, const int minLength)
Definition: BlendPoisson.h:603
vigra::NumericTraits< typename Image::PixelType >::RealPromote ProcessBorderPixel(const int x, const int y, const int dx, const int dy, const Image &image, const Mask &mask, const SeamMask &seam)
Definition: BlendPoisson.h:71
static T min(T x, T y)
Definition: svm.cpp:62
void Multigrid(Image &out, const Image &gradient, const vigra::ImagePyramid< SeamMask > &seamMaskPyramid, int minLen, const float errorThreshold, const int maxIter, const bool doWrap)
Definition: BlendPoisson.h:794