BlendPoisson.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00002 
00012 /*  This program is free software; you can redistribute it and/or
00013 *  modify it under the terms of the GNU General Public
00014 *  License as published by the Free Software Foundation; either
00015 *  version 2 of the License, or (at your option) any later version.
00016 *
00017 *  This software is distributed in the hope that it will be useful,
00018 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00019 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 *  General Public License for more details.
00021 *
00022 *  You should have received a copy of the GNU General Public
00023 *  License along with this software. If not, see
00024 *  <http://www.gnu.org/licenses/>.
00025 *
00026 */
00027 
00028 #ifndef POISSON_BLEND_H
00029 #define POISSON_BLEND_H
00030 
00031 #include <iostream>
00032 #include <vigra/stdimage.hxx>
00033 #include <vigra/convolution.hxx>
00034 #include <vigra/stdconvolution.hxx>
00035 #include <vigra/basicgeometry.hxx>
00036 #include "openmp_vigra.h"
00037 
00038 #define DEBUG_TIMING
00039 
00040 namespace vigra_ext
00041 {
00042 
00043 namespace poisson
00044 {
00045 
00046 namespace detail
00047 {
00048 // helper functions for build gradient map
00049 template <class Image, class Mask>
00050 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)
00051 {
00052     const typename Mask::PixelType m1 = mask[y + dy][x + dx];
00053     const typename Mask::PixelType m2 = mask[y - dy][x - dx];
00054     if (m1 > 0 && m2 > 0)
00055     {
00056         return image[y + dy][x + dx] + image[y - dy][x - dx];
00057     };
00058     // if one of the 2 pixel has no information (outside of image area)
00059     // we extrapolate the value from neighbor pixel
00060     if (m1 > 0)
00061     {
00062         return 2 * image[y + dy][x + dx];
00063     }
00064     else
00065     {
00066         return 2 * image[y - dy][x - dx];
00067     };
00068 };
00069 
00070 template <class Image, class Mask, class SeamMask>
00071 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)
00072 {
00073     const typename SeamMask::PixelType seam1 = seam[y + dy][x + dx];
00074     const typename SeamMask::PixelType seam2 = seam[y - dy][x - dx];
00075     const typename Mask::PixelType mask1 = mask[y + dy][x + dx];
00076     const typename Mask::PixelType mask2 = mask[y - dy][x - dx];
00077     if (seam1 > 0 && seam2 > 0)
00078     {
00079         if (mask1 > 0 && mask2 > 0)
00080         {
00081             return image[y + dy][x + dx] + image[y - dy][x - dx];
00082         };
00083         if (mask1 > 0)
00084         {
00085             return 2 *image[y + dy][x + dx];
00086         }
00087         else
00088         {
00089             return 2 * image[y - dy][x - dx];
00090         };
00091     };
00092     if (seam1 > 0)
00093     {
00094         if (mask1 > 0)
00095         {
00096             return 2 * image[y + dy][x + dx];
00097         }
00098         else
00099         {
00100             return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
00101         };
00102     };
00103     if (seam2 > 0)
00104     {
00105         if (mask2 > 0)
00106         {
00107             return 2 * image[y - dy][x - dx];
00108         }
00109         else
00110         {
00111             return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
00112         };
00113     };
00114     return vigra::NumericTraits<typename vigra::NumericTraits<typename Image::PixelType>::RealPromote>::zero();
00115 };
00116 
00117 template <class Image, class SeamMask>
00118 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)
00119 {
00120     if (seams[y + dy][x + dx] == 1)
00121     {
00122         return image1[offset.y + y + dy][offset.x + x + dx];
00123     };
00124     return vigra::NumericTraits<typename Image::PixelType>::zero();
00125 };
00126 
00127 template <class Image, class SeamMask>
00128 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)
00129 {
00130     const typename SeamMask::PixelType s1 = seams[y + dy][x + dx];
00131     const typename SeamMask::PixelType s2 = seams[y - dy][x - dx];
00132     if (s1 > 1 && s2 > 1)
00133     {
00134         return image[y + dy][x + dx] + image[y - dy][x - dx];
00135     }
00136     else
00137     {
00138         return (2 - std::min<int>(s2, 2))*image[y + dy][x + dx] + (2 - std::min<int>(s1, 2))*image[y - dy][x - dx];
00139     };
00140 };
00141 
00142 template <class Image>
00143 void RestrictErrorToNextLevel(const Image & in, Image & out)
00144 {
00145     typedef typename Image::PixelType ImagePixelType;
00146     const int smallWidth = out.width() - 1;
00147     const int smallHeight = out.height() - 1;
00148 #pragma omp parallel for schedule(dynamic, 100)
00149     for (int y = 0; y < smallHeight; y++)
00150     {
00151         for (int x = 0; x < smallWidth; x++)
00152         {
00153             out[y][x] = in[2 * y][2 * x] + in[2 * y][2 * x + 1] + in[2 * y + 1][2 * x] + in[2 * y + 1][2 * x + 1];
00154         };
00155         // special case of most right pixel
00156         ImagePixelType p2 = 0;
00157         if (2 * smallWidth + 1 < in.width())
00158         {
00159             p2 = in[2 * y][2 * smallWidth + 1];
00160         }
00161         ImagePixelType p3 = 0;
00162         ImagePixelType p4 = 0;
00163         if (2 * y + 1 < in.height())
00164         {
00165             p3 = in[2 * y + 1][2 * smallWidth];
00166             if (2 * smallWidth + 1 < in.width())
00167             {
00168                 p4 = in[2 * y + 1][2 * smallWidth + 1];
00169             };
00170         };
00171         out[y][smallWidth] = in[2 * y][2 * smallWidth] + p2 + p3 + p4;
00172     };
00173     // special case of last line
00174     for (int x = 0; x < smallWidth; x++)
00175     {
00176         ImagePixelType p34 = 0;
00177         if (2 * smallHeight + 1 < in.height())
00178         {
00179             p34 = in[2 * smallHeight + 1][2 * x] + in[2 * smallHeight + 1][2 * x + 1];
00180         };
00181         out[smallHeight][x] = in[2 * smallHeight][2 * x] + in[2 * smallHeight][2 * x + 1] + p34;
00182     };
00183     ImagePixelType p2 = 0;
00184     if (2 * smallWidth + 1 < in.width())
00185     {
00186         p2 = in[2 * smallHeight][2 * smallWidth + 1];
00187     }
00188     ImagePixelType p3 = 0;
00189     ImagePixelType p4 = 0;
00190     if (2 * smallHeight + 1 < in.height())
00191     {
00192         p3 = in[2 * smallHeight + 1][2 * smallWidth];
00193         if (2 * smallWidth + 1 < in.width())
00194         {
00195             p4 = in[2 * smallHeight + 1][2 * smallWidth + 1];
00196         };
00197     };
00198     out[smallHeight][smallWidth] = in[2 * smallHeight][2 * smallWidth] + p2 + p3 + p4;
00199 };
00200 
00201 // internal use by FindEdgesForPoisson
00202 struct FilterEdges
00203 {
00204     vigra::Int8 operator()(float const& vf) const
00205     {
00206         const int v = static_cast<int>(vf);
00207         if (v > 0 && v < 17)
00208         {
00209             return 1;
00210         }
00211         else
00212         {
00213             if (v >= 85)
00214             {
00215                 if (v == 85 || v == 89 || v == 93 || v == 97)
00216                 {
00217                     return 3;
00218                 }
00219                 else
00220                 {
00221                     return 2;
00222                 }
00223             }
00224             else
00225             {
00226                 return 0;
00227             };
00228         };
00229     };
00230 };
00231 
00232 // multi-threaded variant for convolution of images, only 4-neighborhood is considerd
00233 template <class Image1, class Image2>
00234 void SimpleConvolveImage4(const Image1& image1, Image2& image2, const double factor1, const double factor2)
00235 {
00236     vigra_precondition(image1.size() == image2.size(), "ConvolveImage: Image size does not match");
00237     vigra_precondition(image1.width() >= 2 && image1.height() >= 2, "ConvolveImage: Image too small");
00238     const int width = image1.width();
00239     const int height = image1.height();
00240     //special treatment of first line
00241     image2[0][0] = factor1*image1[0][0] + factor2*image1[0][1] + factor2*image1[1][0];
00242     for (int x = 1; x < width - 1; ++x)
00243     {
00244         image2[0][x] = factor1*image1[0][x] + factor2*image1[0][x - 1] + factor2*image1[0][x + 1] + factor2*image1[1][x];
00245     };
00246     image2[0][width - 1] = factor1*image1[0][width - 1] + factor2*image1[0][width - 2] + factor2*image1[1][width - 1];
00247 #pragma omp parallel for
00248     for (int y = 1; y < height - 1; ++y)
00249     {
00250         image2[y][0] = factor1*image1[y][0] + factor2*image1[y - 1][0]
00251             + factor2*image1[y][1] + factor2*image1[y + 1][0];
00252         for (size_t x = 1; x < width - 1; ++x)
00253         {
00254             image2[y][x] = factor1*image1[y][x] + factor2*image1[y - 1][x]
00255                 + factor2*image1[y][x - 1] + factor2*image1[y + 1][x] + factor2*image1[y][x + 1];
00256         }
00257         image2[y][width - 1] = factor1*image1[y][width - 1] + factor2*image1[y - 1][width - 1]
00258             + factor2*image1[y][width - 2] + factor2*image1[y + 1][width - 1];
00259     };
00260     //special treatment of last line
00261     image2[height - 1][0] = factor1*image1[height - 1][0] + factor2*image1[height - 1][1] + factor2*image1[height - 2][0];
00262     for (size_t x = 1; x < width - 1; ++x)
00263     {
00264         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];
00265     };
00266     image2[height - 1][width - 1] = factor1*image1[height - 1][width - 1] + factor2*image1[height - 1][width - 2] + factor2*image1[height - 2][width - 1];
00267 };
00268 
00281 template <class Image>
00282 vigra::Int8Image FindEdgesForPoisson(const Image& input)
00283 {
00284     vigra::Int8Image output(input.size());
00285     SimpleConvolveImage4(input, output, 21, -1);
00286     vigra::omp::transformImage(srcImageRange(output), destImage(output), detail::FilterEdges());
00287     return output;
00288 };
00289 
00290 template <class ComponentType>
00291 double GetRealValue(const ComponentType& val) { return val; }
00292 
00293 template <class ComponentType>
00294 double GetRealValue(const vigra::RGBValue<ComponentType>& val) { return val.magnitude(); }
00295 
00296 template <class Image, class SeamMask>
00297 void SOR(Image& target, const Image& gradient, const SeamMask& seams, const float omega, const float errorThreshold, const int maxIter, const bool doWrap)
00298 {
00299     typedef typename Image::PixelType TargetPixelType;
00300     const int width = target.width();
00301     const int height = target.height();
00302 
00303     // changes in last iteration
00304     double oldError = 0;
00305     for (int j = 0; j < maxIter; j++)
00306     {
00307         // changes in current iteration
00308         double error = 0;
00309         if (seams[0][0]>1)
00310         {
00311             if (doWrap)
00312             {
00313                 const TargetPixelType delta = omega*((gradient[0][0] + target[0][1] + 2 * target[1][0] + target[0][width - 1]) / 4.0f - target[0][0]);
00314                 error += detail::GetRealValue(delta*delta);
00315                 target[0][0] += delta;
00316             }
00317             else
00318             {
00319                 const TargetPixelType delta = omega*((gradient[0][0] + 2 * target[0][1] + 2 * target[1][0]) / 4.0f - target[0][0]);
00320                 error += detail::GetRealValue(delta*delta);
00321                 target[0][0] += delta;
00322             };
00323         };
00324         for (int x = 1; x < width - 1; ++x)
00325         {
00326             if (seams[0][x]>1)
00327             {
00328                 const TargetPixelType delta = omega*((gradient[0][x] + detail::GetBorderValues(x, 0, 1, 0, target, seams) + 2 * target[1][x]) / 4.0f - target[0][x]);
00329                 error += detail::GetRealValue(delta*delta);
00330                 target[0][x] += delta;
00331             };
00332         };
00333         if (seams[0][width - 1] > 1)
00334         {
00335             if (doWrap)
00336             {
00337                 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]);
00338                 error += detail::GetRealValue(delta*delta);
00339                 target[0][width - 1] += delta;
00340             }
00341             else
00342             {
00343                 const TargetPixelType delta = omega*((gradient[0][width - 1] + 2 * target[0][width - 2] + 2 * target[1][width - 1]) / 4.0f - target[0][width - 1]);
00344                 error += detail::GetRealValue(delta*delta);
00345                 target[0][width - 1] += delta;
00346             };
00347         };
00348 #pragma omp parallel for reduction(+: error) schedule(dynamic, 100)
00349         for (int y = 1; y < height - 1; ++y)
00350         {
00351             if (seams[y][0] > 1)
00352             {
00353                 if (doWrap)
00354                 {
00355                     const TargetPixelType delta = omega*((gradient[y][0] + detail::GetBorderValues(0, y, 0, 1, target, seams)
00356                         + target[y][1] + target[y][width - 1]) / 4.0f - target[y][0]);
00357                     error += detail::GetRealValue(delta*delta);
00358                     target[y][0] += delta;
00359                 }
00360                 else
00361                 {
00362                     const TargetPixelType delta = omega*((gradient[y][0] + detail::GetBorderValues(0, y, 0, 1, target, seams) + 2 * target[y][1]) / 4.0f - target[y][0]);
00363                     error += detail::GetRealValue(delta*delta);
00364                     target[y][0] += delta;
00365                 };
00366             };
00367             for (int x = 1; x < width - 1; ++x)
00368             {
00369                 const typename SeamMask::value_type maskValue = seams[y][x];
00370                 if (maskValue > 1)
00371                 {
00372                     if (maskValue == 2)
00373                     {
00374                         // border pixel
00375                         const TargetPixelType sum = detail::GetBorderValues(x, y, 1, 0, target, seams) + detail::GetBorderValues(x, y, 0, 1, target, seams);
00376                         const TargetPixelType delta = omega*((gradient[y][x] + sum) / 4.0f - target[y][x]);
00377                         error += detail::GetRealValue(delta*delta);
00378                         target[y][x] += delta;
00379                     }
00380                     else
00381                     {
00382                         const TargetPixelType sum = target[y + 1][x] + target[y][x + 1] + target[y - 1][x] + target[y][x - 1];
00383                         const TargetPixelType delta = omega*((gradient[y][x] + sum) / 4.0f - target[y][x]);
00384                         error += detail::GetRealValue(delta*delta);
00385                         target[y][x] += delta;
00386                     };
00387                 };
00388             };
00389             if (seams[y][width - 1] > 1)
00390             {
00391                 if (doWrap)
00392                 {
00393                     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]);
00394                     error += detail::GetRealValue(delta*delta);
00395                     target[y][width - 1] += delta;
00396                 }
00397                 else
00398                 {
00399                     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]);
00400                     error += detail::GetRealValue(delta*delta);
00401                     target[y][width - 1] += delta;
00402                 };
00403             };
00404         };
00405         // last line
00406         if (seams[height - 1][0] > 1)
00407         {
00408             if (doWrap)
00409             {
00410                 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]);
00411                 error += detail::GetRealValue(delta*delta);
00412                 target[height - 1][0] += delta;
00413             }
00414             else
00415             {
00416                 const TargetPixelType delta = omega*((gradient[height - 1][0] + 2 * target[height - 2][0] + 2 * target[height - 1][1]) / 4.0f - target[height - 1][0]);
00417                 error += detail::GetRealValue(delta*delta);
00418                 target[height - 1][0] += delta;
00419             };
00420         };
00421         for (int x = 1; x < width - 1; ++x)
00422         {
00423             if (seams[height - 1][x] > 1)
00424             {
00425                 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]);
00426                 error += detail::GetRealValue(delta*delta);
00427                 target[height - 1][x] += delta;
00428             };
00429         };
00430         if (seams[height - 1][width - 1] > 1)
00431         {
00432             if (doWrap)
00433             {
00434                 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]);
00435                 error += detail::GetRealValue(delta*delta);
00436                 target[height - 1][width - 1] += delta;
00437             }
00438             else
00439             {
00440                 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]);
00441                 error += detail::GetRealValue(delta*delta);
00442                 target[height - 1][width - 1] += delta;
00443             };
00444         };
00445 
00446         if (oldError > 0 && log(oldError / error) / log(10.0) < errorThreshold)
00447         {
00448             break;
00449         }
00450         oldError = error;
00451     }
00452 }
00453 
00454 template <class Image, class SeamMask>
00455 void CalcResidualError(Image& error, const Image& target, const Image& gradient, const SeamMask& seam, const bool doWrap)
00456 {
00457     typedef typename Image::PixelType ImagePixelType;
00458     const int width = target.width();
00459     const int height = target.height();
00460 
00461     if (seam[0][0] > 1)
00462     {
00463         if (doWrap)
00464         {
00465             const ImagePixelType sum = 2 * target[1][0] + target[0][1] + target[0][width - 1];
00466             error[0][0] = (4 * target[0][0] - sum - gradient[0][0]);
00467         }
00468         else
00469         {
00470             const ImagePixelType sum = 2 * target[1][0] + 2 * target[0][1];
00471             error[0][0] = (4 * target[0][0] - sum - gradient[0][0]);
00472         };
00473     };
00474     for (int x = 1; x < width - 1; ++x)
00475     {
00476         if (seam[0][x]>1)
00477         {
00478             const ImagePixelType sum = 2 * target[1][x] + detail::GetBorderValues(x, 0, 1, 0, target, seam);
00479             error[0][x] = (4 * target[0][x] - sum - gradient[0][x]);
00480         };
00481     };
00482     if (seam[0][width - 1] > 1)
00483     {
00484         if (doWrap)
00485         {
00486             const ImagePixelType sum = 2 * target[1][width - 1] + target[0][width - 2] + target[0][0];
00487             error[0][width - 1] = (4 * target[0][width - 1] - sum - gradient[0][width - 1]);
00488         }
00489         else
00490         {
00491             const ImagePixelType sum = 2 * target[1][width - 1] + 2 * target[0][width - 2];
00492             error[0][width - 1] = (4 * target[0][width - 1] - sum - gradient[0][width - 1]);
00493         }
00494     };
00495 #pragma omp parallel for schedule(dynamic, 100)
00496     for (int y = 1; y < height - 1; ++y)
00497     {
00498         if (seam[y][0] > 1)
00499         {
00500             if (doWrap)
00501             {
00502                 const ImagePixelType sum = detail::GetBorderValues(0, y, 0, 1, target, seam) + target[y][1] + target[y][width - 1];
00503                 error[y][0] = (4 * target[y][0] - sum - gradient[y][0]);
00504             }
00505             else
00506             {
00507                 const ImagePixelType sum = detail::GetBorderValues(0, y, 0, 1, target, seam) + 2 * target[y][1];
00508                 error[y][0] = (4 * target[y][0] - sum - gradient[y][0]);
00509             };
00510         }
00511         for (int x = 1; x < width - 1; ++x)
00512         {
00513             const typename SeamMask::value_type maskValue = seam[y][x];
00514             if (maskValue > 1)
00515             {
00516                 if (maskValue == 2)
00517                 {
00518                     // border pixel
00519                     const ImagePixelType sum = detail::GetBorderValues(x, y, 1, 0, target, seam) + detail::GetBorderValues(x, y, 0, 1, target, seam);
00520                     error[y][x] = (4 * target[y][x] - sum - gradient[y][x]);
00521                 }
00522                 else
00523                 {
00524                     const ImagePixelType sum = target[y + 1][x] + target[y][x + 1] + target[y - 1][x] + target[y][x - 1];
00525                     error[y][x] = (4 * target[y][x] - sum - gradient[y][x]);
00526                 };
00527             };
00528         };
00529         if (seam[y][width - 1] > 1)
00530         {
00531             if (doWrap)
00532             {
00533                 const ImagePixelType sum = detail::GetBorderValues(width - 1, y, 0, 1, target, seam) + target[y][width - 2] + target[y][0];
00534                 error[y][width - 1] = (4 * target[y][width - 1] - sum - gradient[y][width - 1]);
00535             }
00536             else
00537             {
00538                 const ImagePixelType sum = detail::GetBorderValues(width - 1, y, 0, 1, target, seam) + 2 * target[y][width - 2];
00539                 error[y][width - 1] = (4 * target[y][width - 1] - sum - gradient[y][width - 1]);
00540             };
00541         };
00542     };
00543     // last line
00544     if (seam[height - 1][0] > 1)
00545     {
00546         if (doWrap)
00547         {
00548             const ImagePixelType sum = 2 * target[height - 2][0] + target[height - 1][width - 1] + target[height - 1][1];
00549             error[height - 1][0] = (4 * target[height - 1][0] - sum - gradient[height - 1][0]);
00550         }
00551         else
00552         {
00553             const ImagePixelType sum = 2 * target[height - 2][0] + 2 * target[height - 1][1];
00554             error[height - 1][0] = (4 * target[height - 1][0] - sum - gradient[height - 1][0]);
00555         };
00556     };
00557     for (int x = 1; x < width - 1; ++x)
00558     {
00559         if (seam[height - 1][x]>1)
00560         {
00561             const ImagePixelType sum = 2 * target[height-2][x] + detail::GetBorderValues(x, height - 1, 1, 0, target, seam);
00562             error[height - 1][x] = (4 * target[height - 1][x] - sum - gradient[height - 1][x]);
00563         };
00564     };
00565     if (seam[height - 1][width - 1] > 1)
00566     {
00567         if (doWrap)
00568         {
00569             const ImagePixelType sum = 2 * target[height - 2][width - 1] + target[height - 1][width - 2] + target[height - 1][0];
00570             error[height - 1][width - 1] = (4 * target[height - 1][width - 1] - sum - gradient[height - 1][width - 1]);
00571         }
00572         else
00573         {
00574             const ImagePixelType sum = 2 * target[height - 2][width - 1] + 2 * target[height - 1][width - 2];
00575             error[height - 1][width - 1] = (4 * target[height - 1][width - 1] - sum - gradient[height - 1][width - 1]);
00576         };
00577     };
00578 
00579 }
00580 
00581 } // namespace detail
00582 template <class PixelType>
00583 class MaskGreaterAccessor
00584 {
00585 public:
00586     explicit MaskGreaterAccessor(PixelType val) : m_value(val) {};
00587 
00588     template <class ITERATOR>
00589     PixelType operator()(ITERATOR const & i) const {
00590         if (PixelType(*i) >= m_value)
00591         {
00592             return vigra::NumericTraits<PixelType>::max();
00593         }
00594         else
00595         {
00596             return vigra::NumericTraits<PixelType>::zero();
00597         };
00598     }
00599 
00600     template <class ITERATOR, class DIFFERENCE>
00601     PixelType operator()(ITERATOR const & i, DIFFERENCE d) const
00602     {
00603         if (PixelType(i[d]) >= m_value)
00604         {
00605             return vigra::NumericTraits<PixelType>::max();
00606         }
00607         else
00608         {
00609             return vigra::NumericTraits<PixelType>::zero();
00610         };
00611     }
00612     PixelType m_value;
00613 };
00614 
00615 template <class PixelType>
00616 class MaskSmallerAccessor
00617 {
00618 public:
00619     explicit MaskSmallerAccessor(PixelType val) : m_value(val) {};
00620 
00621     template <class ITERATOR>
00622     PixelType operator()(ITERATOR const & i) const {
00623         if (PixelType(*i) < m_value)
00624         {
00625             return vigra::NumericTraits<PixelType>::max();
00626         }
00627         else
00628         {
00629             return vigra::NumericTraits<PixelType>::zero();
00630         };
00631     }
00632 
00633     template <class ITERATOR, class DIFFERENCE>
00634     PixelType operator()(ITERATOR const & i, DIFFERENCE d) const
00635     {
00636         if (PixelType(i[d]) < m_value)
00637         {
00638             return vigra::NumericTraits<PixelType>::max();
00639         }
00640         else
00641         {
00642             return vigra::NumericTraits<PixelType>::zero();
00643         };
00644     }
00645     PixelType m_value;
00646 };
00647 
00648 template <class Image, class PyramidImage>
00649 void BuildSeamPyramid(const Image& input, vigra::ImagePyramid<PyramidImage>& seams, const int minLength)
00650 {
00651     const int nlevels = (int)ceil(log(std::min(input.width(), input.height()) / (double)minLength) / log(2.0));
00652     seams.resize(0, nlevels, input.size());
00653     Image scaledImage = input;
00654     seams[0] = detail::FindEdgesForPoisson(input);
00655     for (size_t i = 1; i <= seams.highestLevel(); ++i)
00656     {
00657         Image smaller((scaledImage.width() + 1) / 2, (scaledImage.height() + 1) / 2);
00658         vigra::resizeImageNoInterpolation(vigra::srcImageRange(scaledImage), vigra::destImageRange(smaller));
00659         seams[i] = detail::FindEdgesForPoisson(smaller);
00660         scaledImage = smaller;
00661     };
00662 }
00663 
00664 template<class Image, class Mask, class SeamMask, class GradientType>
00665 void BuildGradientMap(const Image& image1, const Image& image2, const Mask& mask2, const SeamMask& seam, GradientType& gradient, const vigra::Point2D& offset, const bool doWrap)
00666 {
00667     typedef typename GradientType::PixelType GradientPixelType;
00668     const int width = image2.width();
00669     const int height = image2.height();
00670     //special treatment of first line
00671     if (seam[0][0] == 2)
00672     {
00673         if (doWrap)
00674         {
00675             GradientPixelType value = 4 * image2[0][0] - image2[0][1] - 2 * image2[1][0] - image2[0][width - 1];
00676             // copy values for Dirichlet boundary condition 
00677             value += detail::GetBorderGradient(0, 0, 1, 0, seam, image1, offset);
00678             value += detail::GetBorderGradient(0, 0, 0, 1, seam, image1, offset);
00679             value += detail::GetBorderGradient(width - 2, 0, 1, 0, seam, image1, offset);
00680             gradient[0][0] = value;
00681         }
00682         else
00683         {
00684             GradientPixelType value = 4 * image2[0][0] - 2 * image2[0][1] - 2 * image2[1][0];
00685             // copy values for Dirichlet boundary condition 
00686             value += detail::GetBorderGradient(0, 0, 1, 0, seam, image1, offset);
00687             value += detail::GetBorderGradient(0, 0, 0, 1, seam, image1, offset);
00688             gradient[0][0] = value;
00689         };
00690     };
00691     for (int x = 1; x < width - 1; ++x)
00692     {
00693         if (seam[0][x] == 2)
00694         {
00695             GradientPixelType value = 4 * image2[0][x] - detail::ProcessBorderPixel(x, 0, 1, 0, image2, mask2, seam) - 2 * image2[1][x];
00696             value += detail::GetBorderGradient(x, 0, -1, 0, seam, image1, offset);
00697             value += detail::GetBorderGradient(x, 0, 1, 0, seam, image1, offset);
00698             value += detail::GetBorderGradient(x, 0, 0, 1, seam, image1, offset);
00699             gradient[0][x] = value;
00700         };
00701     };
00702     if (seam[0][width - 1] == 2)
00703     {
00704         if (doWrap)
00705         {
00706             GradientPixelType value = 4 * image2[0][width - 1] - image2[0][width - 2] - 2 * image2[1][width - 1] - image2[0][0];
00707             value += detail::GetBorderGradient(width - 1, 0, -1, 0, seam, image1, offset);
00708             value += detail::GetBorderGradient(width - 1, 0, 0, 1, seam, image1, offset);
00709             value += detail::GetBorderGradient(1, 0, -1, 0, seam, image1, offset);
00710             gradient[0][width - 1] = value;
00711         }
00712         else
00713         {
00714             GradientPixelType value = 4 * image2[0][width - 1] - 2 * image2[0][width - 2] - 2 * image2[1][width - 1];
00715             value += detail::GetBorderGradient(width - 1, 0, -1, 0, seam, image1, offset);
00716             value += detail::GetBorderGradient(width - 1, 0, 0, 1, seam, image1, offset);
00717             gradient[0][width - 1] = value;
00718         };
00719     };
00720 #pragma omp parallel for
00721     for (int y = 1; y < height - 1; ++y)
00722     {
00723         if (seam[y][0] == 2)
00724         {
00725             if (doWrap)
00726             {
00727                 GradientPixelType value = 4 * image2[y][0] - detail::ProcessBorderPixel(0, y, 0, 1, image2, mask2, seam) - image2[y][1] - image2[y][width - 1];
00728                 value += detail::GetBorderGradient(0, y, 1, 0, seam, image1, offset);
00729                 value += detail::GetBorderGradient(0, y, 0, 1, seam, image1, offset);
00730                 value += detail::GetBorderGradient(0, y, 0, -1, seam, image1, offset);
00731                 value += detail::GetBorderGradient(width - 2, y, 1, 0, seam, image1, offset);
00732                 gradient[y][0] = value;
00733             }
00734             else
00735             {
00736                 GradientPixelType value = 4 * image2[y][0] - detail::ProcessBorderPixel(0, y, 0, 1, image2, mask2, seam) - 2 * image2[y][1];
00737                 value += detail::GetBorderGradient(0, y, 1, 0, seam, image1, offset);
00738                 value += detail::GetBorderGradient(0, y, 0, 1, seam, image1, offset);
00739                 value += detail::GetBorderGradient(0, y, 0, -1, seam, image1, offset);
00740                 gradient[y][0] = value;
00741             };
00742         };
00743         for (int x = 1; x < width - 1; ++x)
00744         {
00745             const typename SeamMask::PixelType seamVal = seam[y][x];
00746             if (seamVal>1)
00747             {
00748                 GradientPixelType value = 4.0 * image2[y][x];
00749                 if (seamVal == 3)
00750                 {
00751                     value -= detail::ProcessNeighborPixels(x, y, 1, 0, image2, mask2);
00752                     value -= detail::ProcessNeighborPixels(x, y, 0, 1, image2, mask2);
00753                 }
00754                 else
00755                 {
00756                     // seamVal==2, border pixel
00757                     value -= detail::ProcessBorderPixel(x, y, 1, 0, image2, mask2, seam);
00758                     value -= detail::ProcessBorderPixel(x, y, 0, 1, image2, mask2, seam);
00759                 };
00760                 // copy values for Dirichlet boundary condition 
00761                 value += detail::GetBorderGradient(x, y, 1, 0, seam, image1, offset);
00762                 value += detail::GetBorderGradient(x, y, 0, 1, seam, image1, offset);
00763                 value += detail::GetBorderGradient(x, y, -1, 0, seam, image1, offset);
00764                 value += detail::GetBorderGradient(x, y, 0, -1, seam, image1, offset);
00765                 gradient[y][x] = value;
00766             }
00767         };
00768         if (seam[y][width - 1] == 2)
00769         {
00770             if (doWrap)
00771             {
00772                 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];
00773                 value += detail::GetBorderGradient(width - 1, y, -1, 0, seam, image1, offset);
00774                 value += detail::GetBorderGradient(width - 1, y, 0, 1, seam, image1, offset);
00775                 value += detail::GetBorderGradient(width - 1, y, 0, -1, seam, image1, offset);
00776                 value += detail::GetBorderGradient(0, y, -1, 0, seam, image1, offset);
00777                 gradient[y][width - 1] = value;
00778             }
00779             else
00780             {
00781                 GradientPixelType value = 4.0 * image2[y][width - 1] - detail::ProcessBorderPixel(width - 1, y, 0, 1, image2, mask2, seam) - 2 * image2[y][width - 2];
00782                 value += detail::GetBorderGradient(width - 1, y, -1, 0, seam, image1, offset);
00783                 value += detail::GetBorderGradient(width - 1, y, 0, 1, seam, image1, offset);
00784                 value += detail::GetBorderGradient(width - 1, y, 0, -1, seam, image1, offset);
00785                 gradient[y][width - 1] = value;
00786             };
00787         };
00788     };
00789     //special treatment of last line
00790     if (seam[height - 1][0] == 2)
00791     {
00792         if (doWrap)
00793         {
00794             GradientPixelType value = 4 * image2[height - 1][0] - image2[height - 1][1] - 2 * image2[height - 2][0] - image2[height - 1][width - 1];
00795             value += detail::GetBorderGradient(0, height - 1, 1, 0, seam, image1, offset);
00796             value += detail::GetBorderGradient(0, height - 1, 0, -1, seam, image1, offset);
00797             value += detail::GetBorderGradient(width - 2, height - 1, 1, 0, seam, image1, offset);
00798             gradient[height - 1][0] = value;
00799         }
00800         else
00801         {
00802             GradientPixelType value = 4 * image2[height - 1][0] - 2 * image2[height - 1][1] - 2 * image2[height - 2][0];
00803             value += detail::GetBorderGradient(0, height - 1, 1, 0, seam, image1, offset);
00804             value += detail::GetBorderGradient(0, height - 1, 0, -1, seam, image1, offset);
00805             gradient[height - 1][0] = value;
00806         };
00807     };
00808     for (size_t x = 1; x < width - 1; ++x)
00809     {
00810         if (seam[height - 1][x] == 2)
00811         {
00812             GradientPixelType value = 4 * image2[height - 1][x] - detail::ProcessBorderPixel(x, height - 1, 1, 0, image2, mask2, seam) - 2 * image2[height - 2][x];
00813             value += detail::GetBorderGradient(x, height - 1, -1, 0, seam, image1, offset);
00814             value += detail::GetBorderGradient(x, height - 1, 1, 0, seam, image1, offset);
00815             value += detail::GetBorderGradient(x, height - 1, 0, -1, seam, image1, offset);
00816             gradient[height - 1][x] = value;
00817         };
00818     };
00819     if (seam[height - 1][width - 1] == 2)
00820     {
00821         if (doWrap)
00822         {
00823             GradientPixelType value = 4 * image2[height - 1][width - 1] - image2[height - 1][width - 2] - 2 * image2[height - 2][width - 1] - image2[height - 1][0];
00824             value += detail::GetBorderGradient(width - 1, height - 1, -1, 0, seam, image1, offset);
00825             value += detail::GetBorderGradient(width - 1, height - 1, 0, -1, seam, image1, offset);
00826             value += detail::GetBorderGradient(1, height - 1, -1, 0, seam, image1, offset);
00827             gradient[height - 1][width - 1] = value;
00828         }
00829         else
00830         {
00831             GradientPixelType value = 4 * image2[height - 1][width - 1] - 2 * image2[height - 1][width - 2] - 2 * image2[height - 2][width - 1];
00832             value += detail::GetBorderGradient(width - 1, height - 1, -1, 0, seam, image1, offset);
00833             value += detail::GetBorderGradient(width - 1, height - 1, 0, -1, seam, image1, offset);
00834             gradient[height - 1][width - 1] = value;
00835         };
00836     };
00837 };
00838 
00839 template <class Image, class SeamMask>
00840 void Multigrid(Image& out, const Image& gradient, const vigra::ImagePyramid<SeamMask>& seamMaskPyramid, int minLen, const float errorThreshold, const int maxIter, const bool doWrap)
00841 {
00842     const int width = out.width();
00843     const int height = out.height();
00844 
00845     if (width < minLen || height < minLen)
00846     {
00847         return;
00848     }
00849     Image err(width, height);
00850     Image err2((width + 1) / 2, (height + 1) / 2);
00851     Image out2(err2.size());
00852     int maskIndex = -1;
00853     for (int i = 0; i <= seamMaskPyramid.highestLevel(); ++i)
00854     {
00855         if (out.size() == seamMaskPyramid[i].size())
00856         {
00857             maskIndex = i;
00858             break;
00859         };
00860     }
00861     if (maskIndex == -1)
00862     {
00863         std::cout << "ERROR: No suitable mask, this should not happen." << std::endl
00864             << "searching " << out.size() << ", finest " << seamMaskPyramid[seamMaskPyramid.highestLevel()].size() << std::endl;
00865         return;
00866     };
00867     // pre-smoothing
00868     detail::SOR(out, gradient, seamMaskPyramid[maskIndex], 1.95f, errorThreshold, maxIter, doWrap);
00869     detail::CalcResidualError(err, out, gradient, seamMaskPyramid[maskIndex], doWrap);    // Fehler berechnen
00870     detail::RestrictErrorToNextLevel(err, err2);
00871     Multigrid(out2, err2, seamMaskPyramid, minLen, errorThreshold, maxIter, doWrap);
00872     // again W cycle
00873     Multigrid(out2, err2, seamMaskPyramid, minLen, errorThreshold, maxIter, doWrap); 
00874     vigra::resizeImageLinearInterpolation(srcImageRange(out2), destImageRange(err)); 
00875     vigra::omp::combineTwoImagesIf(vigra::srcImageRange(out), vigra::srcImage(err),
00876         vigra::srcImage(seamMaskPyramid[maskIndex], MaskGreaterAccessor<typename SeamMask::PixelType>(2)),
00877         vigra::destImage(out),
00878         vigra::functor::Arg1() - vigra::functor::Arg2());
00879     // post smoothing
00880     detail::SOR(out, gradient, seamMaskPyramid[maskIndex], 1.95f, 2.0f * errorThreshold, maxIter, doWrap);
00881     return;
00882 }
00883 
00884 } // namespace poisson
00885 } // namespace vigra_ext
00886 
00887 #endif //POISSON_BLEND_H

Generated on 20 Apr 2018 for Hugintrunk by  doxygen 1.4.7