StitchWatershed.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 #include <vigra/seededregiongrowing.hxx>
00029 #include <vigra/convolution.hxx>
00030 #include "vigra_ext/BlendPoisson.h"
00031 #ifdef HAVE_OPENMP
00032 #include <omp.h>
00033 #endif
00034 #include "openmp_vigra.h"
00035 
00036 namespace vigra_ext
00037 {
00038     namespace detail
00039     {
00040         // some helper functions 
00041         struct BuildSeed
00042         {
00043             template <class PixelType>
00044             PixelType operator()(PixelType const& v1, PixelType const& v2) const
00045             {
00046                 return (v1 & 1) | (v2 & 2);
00047             }
00048         };
00049 
00050         template <typename t>
00051         inline double square(t x)
00052         {
00053             return x * x;
00054         }
00055 
00056         struct BuildDiff
00057         {
00058             template <class PixelType>
00059             double operator()(PixelType const& v1, PixelType const& v2) const
00060             {
00061                 return std::abs(static_cast<double>(v1 - v2));
00062             };
00063             template <class PixelType>
00064             double operator()(vigra::RGBValue<PixelType> const& v1, vigra::RGBValue<PixelType> const& v2) const
00065             {
00066                 return sqrt(square(v1.red() - v2.red()) + square(v1.green() - v2.green()) + square(v1.blue() - v2.blue()));
00067             };
00068         };
00069 
00070         struct CombineMasks
00071         {
00072             template <class PixelType>
00073             PixelType operator()(PixelType const& v1, PixelType const& v2) const
00074             {
00075                 if ((v1 & 2) & v2)
00076                 {
00077                     return vigra::NumericTraits<PixelType>::max();
00078                 }
00079                 else
00080                 {
00081                     return vigra::NumericTraits<PixelType>::zero();
00082                 };
00083             };
00084         };
00085 
00086         struct CombineMasksForPoisson
00087         {
00088             template <class PixelType>
00089             PixelType operator()(PixelType const& v1, PixelType const& v2) const
00090             {
00091                 if ((v2 & 2) & v1)
00092                 {
00093                     return 5;
00094                 }
00095                 else
00096                 {
00097                     if (v1 > 0)
00098                     {
00099                         return v2;
00100                     }
00101                     else
00102                     {
00103                         return vigra::NumericTraits<PixelType>::zero();
00104                     };
00105                 };
00106             };
00107         };
00108 
00109         template <class ImageType>
00110         ImageType ResizeImage(const ImageType& image, const vigra::Size2D& newSize)
00111         {
00112             ImageType newImage(std::max(image.size().width(), newSize.width()), std::max(image.size().height(), newSize.height()));
00113             vigra::omp::copyImage(vigra::srcImageRange(image), vigra::destImage(newImage));
00114             return newImage;
00115         };
00116     }; // namespace detail
00117     
00118     template <class ImageType, class MaskType>
00119     void MergeImages(ImageType& image1, MaskType& mask1, const ImageType& image2, const MaskType& mask2, const vigra::Diff2D offset, const bool wrap, const bool hardSeam)
00120     {
00121         const vigra::Point2D offsetPoint(offset);
00122         const vigra::Rect2D offsetRect(offsetPoint, mask2.size());
00123         //increase image size if necessary
00124         if (image1.width() < offsetRect.lowerRight().x || image1.height() < offsetRect.lowerRight().y)
00125         {
00126             image1 = detail::ResizeImage(image1, vigra::Size2D(offsetRect.lowerRight()));
00127             mask1 = detail::ResizeImage(mask1, image1.size());
00128         }
00129         // generate seed mask
00130         vigra::BImage labels(image2.size());
00131         // create a seed mask
00132         // value 0: pixel is not contained in image 1 or 2
00133         // value 1: pixel contains only information from image 1
00134         // value 2: pixel contains only information from image 2
00135         // value 3: pixel contains information from image 1 and 2
00136         vigra::omp::combineTwoImages(vigra::srcImageRange(mask1, offsetRect), vigra::srcImage(mask2), vigra::destImage(labels), detail::BuildSeed());
00137         // find bounding rectangles for all values
00138         vigra::ArrayOfRegionStatistics<vigra::FindBoundingRectangle> roi(3);
00139         vigra::inspectTwoImages(vigra::srcIterRange<vigra::Diff2D>(vigra::Diff2D(0, 0), labels.size()), vigra::srcImage(labels), roi);
00140         // handle some special cases
00141         if (roi.regions[3].size().area() == 0)
00142         {
00143             // images do not overlap, simply copy image2 into image1
00144             vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(mask2), vigra::destImage(image1, offsetPoint, image1.accessor()));
00145             // now merge masks
00146             vigra::copyImageIf(vigra::srcImageRange(mask2), vigra::srcImage(mask2), vigra::destImage(mask1, offsetPoint, mask1.accessor()));
00147             return;
00148         };
00149         if (roi.regions[2].size().area() == 0)
00150         {
00151             // image 2 is fully overlapped by image 1
00152             // we don't need to do anything
00153             return;
00154         };
00155         if (roi.regions[1].size().area() == 0)
00156         {
00157             // image 1 is fully overlapped by image 2
00158             // copy image 2 into output
00159             vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(mask2), vigra::destImage(image1, offsetPoint, image1.accessor()));
00160             // now merge masks
00161             vigra::copyImageIf(vigra::srcImageRange(mask2), vigra::srcImage(mask2), vigra::destImage(mask1, offsetPoint, mask1.accessor()));
00162             return;
00163         }
00164         const double smoothRadius = std::max(1.0, std::max(roi.regions[3].size().width(), roi.regions[3].size().height()) / 1000.0);
00165         const bool doWrap = wrap && (roi.regions[3].size().width() == image1.width());
00166         // build seed map
00167         vigra::omp::transformImage(vigra::srcImageRange(labels), vigra::destImage(labels), vigra::functor::Arg1() % vigra::functor::Param(3));
00168         // build difference, only consider overlapping area
00169         // increase size by 1 pixel in each direction if possible
00170         vigra::Point2D p1(roi.regions[3].upperLeft);
00171         if (p1.x > 0)
00172         {
00173             --(p1.x);
00174         };
00175         if (p1.y > 0)
00176         {
00177             --(p1.y);
00178         };
00179         vigra::Point2D p2(roi.regions[3].lowerRight);
00180         if (p2.x + 1 < image2.width())
00181         {
00182             ++(p2.x);
00183         };
00184         if (p2.y + 1 < image2.height())
00185         {
00186             ++(p2.y);
00187         };
00188         vigra::DImage diff(p2 - p1);
00189         const vigra::Rect2D rect1(offsetPoint + p1, diff.size());
00190         // build difference map
00191         vigra::omp::combineTwoImages(vigra::srcImageRange(image1, rect1), vigra::srcImage(image2, p1), vigra::destImage(diff), detail::BuildDiff());
00192         // scale to 0..255 to faster watershed
00193         vigra::FindMinMax<double> diffMinMax;
00194         vigra::inspectImage(vigra::srcImageRange(diff), diffMinMax);
00195         diffMinMax.max = std::min<double>(diffMinMax.max, 0.25f * vigra::NumericTraits<typename vigra::NumericTraits<typename ImageType::PixelType>::ValueType>::max());
00196         vigra::BImage diffByte(diff.size());
00197         vigra::omp::transformImage(vigra::srcImageRange(diff), vigra::destImage(diffByte), vigra::functor::Param(255) - vigra::functor::Param(255.0f / diffMinMax.max)*vigra::functor::Arg1());
00198         diff.resize(0, 0);
00199         // run watershed algorithm
00200         vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<vigra::UInt8> > stats(3);
00201         if (doWrap)
00202         {
00203             // handle wrapping
00204             const int oldWidth = labels.width();
00205             const int oldHeight = labels.height();
00206             vigra::BImage labelsWrapped(oldWidth * 2, oldHeight);
00207             vigra::omp::copyImage(vigra::srcImageRange(labels), vigra::destImage(labelsWrapped));
00208             vigra::omp::copyImage(labels.upperLeft(), labels.lowerRight(), labels.accessor(), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), labelsWrapped.accessor());
00209             vigra::BImage diffWrapped(oldWidth * 2, diffByte.height());
00210             vigra::omp::copyImage(vigra::srcImageRange(diffByte), vigra::destImage(diffWrapped));
00211             vigra::omp::copyImage(diffByte.upperLeft(), diffByte.lowerRight(), diffByte.accessor(), diffWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), diffWrapped.accessor());
00212             // apply gaussian smoothing with size depending radius
00213             // we need a minimum size of window to apply gaussianSmoothing
00214             if (diffWrapped.width() > 3 * smoothRadius && diffWrapped.height() > 3 * smoothRadius)
00215             {
00216                 vigra::gaussianSmoothing(vigra::srcImageRange(diffWrapped), vigra::destImage(diffWrapped), smoothRadius);
00217             };
00218             vigra::fastSeededRegionGrowing(vigra::srcImageRange(diffWrapped), vigra::destImage(labelsWrapped, p1), stats, vigra::CompleteGrow, vigra::FourNeighborCode(), 255);
00219             vigra::omp::copyImage(labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth / 2, 0), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, oldHeight), labelsWrapped.accessor(),
00220                 labels.upperLeft() + vigra::Diff2D(oldWidth / 2, 0), labels.accessor());
00221             vigra::omp::copyImage(labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth + oldWidth / 2, oldHeight), labelsWrapped.accessor(),
00222                 labels.upperLeft(), labels.accessor());
00223         }
00224         else
00225         {
00226             // apply gaussian smoothing with size depending radius
00227             // we need a minimum size of window to apply gaussianSmoothing
00228             if (diffByte.width() > 3 * smoothRadius && diffByte.height() > 3 * smoothRadius)
00229             {
00230                 vigra::gaussianSmoothing(vigra::srcImageRange(diffByte), vigra::destImage(diffByte), smoothRadius);
00231             };
00232             vigra::fastSeededRegionGrowing(vigra::srcImageRange(diffByte), vigra::destImage(labels, p1), stats, vigra::CompleteGrow, vigra::FourNeighborCode(), 255);
00233         };
00234         // now we can merge the images
00235         // merging the mask is straightforward
00236         vigra::initImageIf(vigra::destImageRange(mask1, offsetRect), vigra::srcImage(mask2), vigra::NumericTraits<typename MaskType::value_type>::max());
00237         if (hardSeam)
00238         {
00239             // the watershed algorithm could also reached area where no informations are available
00240             vigra::omp::combineTwoImages(vigra::srcImageRange(labels), vigra::srcImage(mask2), vigra::destImage(labels), detail::CombineMasks());
00241             // now we can merge the images
00242             vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(labels), vigra::destImage(image1, offsetPoint));
00243         }
00244         else
00245         {
00246             // find all boundaries in new mask
00247             // first filter out unused pixel the watershed algorithm has also processed
00248             vigra::omp::combineTwoImages(vigra::srcImageRange(mask1, offsetRect), vigra::srcImage(labels), vigra::destImage(labels), detail::CombineMasksForPoisson());
00249             // labels has now the following values:
00250             // 0: no image here
00251             // 1: use information from image 1
00252             // 5: use information from image 2
00253             // mark edges in labels for solving Poisson equation with different boundary conditions
00254             vigra::ImagePyramid<vigra::Int8Image> seams;
00255             const int minLength = 8;
00256             vigra_ext::poisson::BuildSeamPyramid(labels, seams, minLength);
00257             // create gradient map
00258             typedef typename vigra::NumericTraits<typename ImageType::PixelType>::RealPromote ImageRealPixelType;
00259             vigra::BasicImage<ImageRealPixelType> gradient(image2.size());
00260             vigra::BasicImage<ImageRealPixelType> target(image2.size());
00261             // build gradient map with special handling of both boundary conditions
00262             vigra_ext::poisson::BuildGradientMap(image1, image2, mask2, seams[0], gradient, offsetPoint, doWrap);
00263             // we start with the values of the image2 as begin
00264             vigra::omp::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(seams[0], vigra_ext::poisson::MaskGreaterAccessor<vigra::Int8>(2)), vigra::destImage(target));
00265             // solve poisson equation
00266             vigra_ext::poisson::Multigrid(target, gradient, seams, minLength, 0.01f, 500, doWrap);
00267             // copy result back into output
00268             vigra::omp::copyImageIf(vigra::srcImageRange(target), vigra::srcImage(seams[0], vigra_ext::poisson::MaskGreaterAccessor<vigra::Int8>(2)), vigra::destImage(image1, offsetPoint));
00269         };
00270     };
00271 
00272 }

Generated on 22 Jan 2018 for Hugintrunk by  doxygen 1.4.7