Hugintrunk  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
StitchWatershed.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 #include <vigra/seededregiongrowing.hxx>
29 #include <vigra/convolution.hxx>
30 #include "vigra_ext/BlendPoisson.h"
31 #ifdef HAVE_OPENMP
32 #include <omp.h>
33 #endif
34 #include "openmp_vigra.h"
35 
36 namespace vigra_ext
37 {
38  namespace detail
39  {
40  // some helper functions
41  struct BuildSeed
42  {
43  template <class PixelType>
44  PixelType operator()(PixelType const& v1, PixelType const& v2) const
45  {
46  return (v1 & 1) | (v2 & 2);
47  }
48  };
49 
50  template <typename t>
51  inline double square(t x)
52  {
53  return x * x;
54  }
55 
56  struct BuildDiff
57  {
58  template <class PixelType>
59  double operator()(PixelType const& v1, PixelType const& v2) const
60  {
61  return std::abs(static_cast<double>(v1 - v2));
62  };
63  template <class PixelType>
64  double operator()(vigra::RGBValue<PixelType> const& v1, vigra::RGBValue<PixelType> const& v2) const
65  {
66  return sqrt(square(v1.red() - v2.red()) + square(v1.green() - v2.green()) + square(v1.blue() - v2.blue()));
67  };
68  };
69 
70  struct CombineMasks
71  {
72  template <class PixelType>
73  PixelType operator()(PixelType const& v1, PixelType const& v2) const
74  {
75  if ((v1 & 2) & v2)
76  {
78  }
79  else
80  {
81  return vigra::NumericTraits<PixelType>::zero();
82  };
83  };
84  };
85 
87  {
88  template <class PixelType>
89  PixelType operator()(PixelType const& v1, PixelType const& v2) const
90  {
91  if ((v2 & 2) & v1)
92  {
93  return 5;
94  }
95  else
96  {
97  if (v1 > 0)
98  {
99  return v2;
100  }
101  else
102  {
103  return vigra::NumericTraits<PixelType>::zero();
104  };
105  };
106  };
107  };
108 
109  template <class ImageType>
110  ImageType ResizeImage(const ImageType& image, const vigra::Size2D& newSize)
111  {
112  ImageType newImage(std::max(image.size().width(), newSize.width()), std::max(image.size().height(), newSize.height()));
114  return newImage;
115  };
116  }; // namespace detail
117 
118  // blend 2 images with Poisson blending
119  // labels contains the map which images should deliver the input, see comment in MergeImages for the values
120  template <class ImageType, class MaskType>
121  void PoissonBlend(ImageType& image1, const ImageType& image2, const MaskType& mask2, const vigra::BImage& labels, const vigra::Point2D& offsetPoint, const bool doWrap)
122  {
123  // mark edges in labels for solving Poisson equation with different boundary conditions
124  vigra::ImagePyramid<vigra::Int8Image> seams;
125  const int minLength = 8;
126  vigra_ext::poisson::BuildSeamPyramid(labels, seams, minLength);
127  // create gradient map
128  typedef typename vigra::NumericTraits<typename ImageType::PixelType>::RealPromote ImageRealPixelType;
129  vigra::BasicImage<ImageRealPixelType> gradient(image2.size());
130  vigra::BasicImage<ImageRealPixelType> target(image2.size());
131  // build gradient map with special handling of both boundary conditions
132  vigra_ext::poisson::BuildGradientMap(image1, image2, mask2, seams[0], gradient, offsetPoint, doWrap);
133  // we start with the values of the image2 as begin
135  // solve poisson equation
136  vigra_ext::poisson::Multigrid(target, gradient, seams, minLength, 0.1f, 500, doWrap);
137  // copy result back into output
139  }
140 
141  template <class ImageType, class MaskType>
142  void MergeImages(ImageType& image1, MaskType& mask1, const ImageType& image2, const MaskType& mask2, const vigra::Diff2D offset, const bool wrap, const bool hardSeam)
143  {
144  const vigra::Point2D offsetPoint(offset);
145  const vigra::Rect2D offsetRect(offsetPoint, mask2.size());
146  //increase image size if necessary
147  if (image1.width() < offsetRect.lowerRight().x || image1.height() < offsetRect.lowerRight().y)
148  {
149  image1 = detail::ResizeImage(image1, vigra::Size2D(offsetRect.lowerRight()));
150  mask1 = detail::ResizeImage(mask1, image1.size());
151  }
152  // generate seed mask
153  vigra::BImage labels(image2.size());
154  // create a seed mask
155  // value 0: pixel is not contained in image 1 or 2
156  // value 1: pixel contains only information from image 1
157  // value 2: pixel contains only information from image 2
158  // value 3: pixel contains information from image 1 and 2
160  // find bounding rectangles for all values
161  vigra::ArrayOfRegionStatistics<vigra::FindBoundingRectangle> roi(3);
162  vigra::inspectTwoImages(vigra::srcIterRange<vigra::Diff2D>(vigra::Diff2D(0, 0), labels.size()), vigra::srcImage(labels), roi);
163  // handle some special cases
164  const bool haveOverlap = roi.regions[3].size().area() > 0;
165  if (!haveOverlap && hardSeam)
166  {
167  // images do not overlap, simply copy image2 into image1
168  vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(mask2), vigra::destImage(image1, offsetPoint, image1.accessor()));
169  // now merge masks
170  vigra::copyImageIf(vigra::srcImageRange(mask2), vigra::srcImage(mask2), vigra::destImage(mask1, offsetPoint, mask1.accessor()));
171  return;
172  };
173  if (roi.regions[2].size().area() == 0)
174  {
175  // image 2 is fully overlapped by image 1
176  // we don't need to do anything
177  return;
178  };
179  if (roi.regions[1].size().area() == 0)
180  {
181  // image 1 is fully overlapped by image 2
182  // copy image 2 into output
183  vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(mask2), vigra::destImage(image1, offsetPoint, image1.accessor()));
184  // now merge masks
185  vigra::copyImageIf(vigra::srcImageRange(mask2), vigra::srcImage(mask2), vigra::destImage(mask1, offsetPoint, mask1.accessor()));
186  return;
187  }
188  const double smoothRadius = std::max(1.0, std::max(roi.regions[3].size().width(), roi.regions[3].size().height()) / 1000.0);
189  // build seed map
190  vigra::omp::transformImage(vigra::srcImageRange(labels), vigra::destImage(labels), vigra::functor::Arg1() % vigra::functor::Param(3));
191  if (haveOverlap)
192  {
193  // build difference, only consider overlapping area
194  // increase size by 1 pixel in each direction if possible
195  vigra::Point2D p1(roi.regions[3].upperLeft);
196  if (p1.x > 0)
197  {
198  --(p1.x);
199  };
200  if (p1.y > 0)
201  {
202  --(p1.y);
203  };
204  vigra::Point2D p2(roi.regions[3].lowerRight);
205  if (p2.x + 1 < image2.width())
206  {
207  ++(p2.x);
208  };
209  if (p2.y + 1 < image2.height())
210  {
211  ++(p2.y);
212  };
213  vigra::DImage diff(p2 - p1);
214  const vigra::Rect2D rect1(offsetPoint + p1, diff.size());
215  // build difference map
217  // scale to 0..255 to faster watershed
218  vigra::FindMinMax<double> diffMinMax;
219  vigra::inspectImage(vigra::srcImageRange(diff), diffMinMax);
220  diffMinMax.max = std::min<double>(diffMinMax.max, 0.25f * vigra::NumericTraits<typename vigra::NumericTraits<typename ImageType::PixelType>::ValueType>::max());
221  vigra::BImage diffByte(diff.size());
222  vigra::omp::transformImage(vigra::srcImageRange(diff), vigra::destImage(diffByte), vigra::functor::Param(255) - vigra::functor::Param(255.0f / diffMinMax.max) * vigra::functor::Arg1());
223  diff.resize(0, 0);
224  // run watershed algorithm
225  vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<vigra::UInt8> > stats(3);
226  if (wrap && (roi.regions[3].size().width() == image1.width()))
227  {
228  // handle wrapping
229  // only if requested and if overlap extends above 360 deg boundary
230  const int oldWidth = labels.width();
231  const int oldHeight = labels.height();
232  vigra::BImage labelsWrapped(oldWidth * 2, oldHeight);
234  vigra::omp::copyImage(labels.upperLeft(), labels.lowerRight(), labels.accessor(), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), labelsWrapped.accessor());
235  vigra::BImage diffWrapped(oldWidth * 2, diffByte.height());
237  vigra::omp::copyImage(diffByte.upperLeft(), diffByte.lowerRight(), diffByte.accessor(), diffWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), diffWrapped.accessor());
238  // apply gaussian smoothing with size depending radius
239  // we need a minimum size of window to apply gaussianSmoothing
240  if (diffWrapped.width() > 3 * smoothRadius && diffWrapped.height() > 3 * smoothRadius)
241  {
242  vigra::gaussianSmoothing(vigra::srcImageRange(diffWrapped), vigra::destImage(diffWrapped), smoothRadius);
243  };
244  vigra::fastSeededRegionGrowing(vigra::srcImageRange(diffWrapped), vigra::destImage(labelsWrapped, p1), stats, vigra::CompleteGrow, vigra::FourNeighborCode(), 255);
245  vigra::omp::copyImage(labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth / 2, 0), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, oldHeight), labelsWrapped.accessor(),
246  labels.upperLeft() + vigra::Diff2D(oldWidth / 2, 0), labels.accessor());
247  vigra::omp::copyImage(labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth, 0), labelsWrapped.upperLeft() + vigra::Diff2D(oldWidth + oldWidth / 2, oldHeight), labelsWrapped.accessor(),
248  labels.upperLeft(), labels.accessor());
249  }
250  else
251  {
252  // apply gaussian smoothing with size depending radius
253  // we need a minimum size of window to apply gaussianSmoothing
254  if (diffByte.width() > 3 * smoothRadius && diffByte.height() > 3 * smoothRadius)
255  {
256  vigra::gaussianSmoothing(vigra::srcImageRange(diffByte), vigra::destImage(diffByte), smoothRadius);
257  };
258  vigra::fastSeededRegionGrowing(vigra::srcImageRange(diffByte), vigra::destImage(labels, p1), stats, vigra::CompleteGrow, vigra::FourNeighborCode(), 255);
259  };
260  };
261  // now we can merge the images
262  // merging the mask is straightforward
264  if (hardSeam)
265  {
266  // the watershed algorithm could also reached area where no informations are available
268  // now we can merge the images
269  vigra::copyImageIf(vigra::srcImageRange(image2), vigra::srcImage(labels), vigra::destImage(image1, offsetPoint));
270  }
271  else
272  {
273  // find all boundaries in new mask
274  // first filter out unused pixel the watershed algorithm has also processed
276  const bool doWrapBlending = wrap && (
277  // overlap regions wrap around 360 deg border
278  (roi.regions[3].size().width() == image1.width()) ||
279  // touching region wrap around 360 deg border
280  (roi.regions[2].size().width() == image1.width())
281  );
282  // labels has now the following values:
283  // 0: no image here
284  // 1: use information from image 1
285  // 5: use information from image 2
286  PoissonBlend(image1, image2, mask2, labels, offsetPoint, doWrapBlending);
287  };
288  };
289 
290 }
PixelType operator()(PixelType const &v1, PixelType const &v2) const
double square(t x)
void copyImageIf(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, MaskImageIterator mask_upperleft, MaskAccessor mask_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc)
Definition: openmp_vigra.h:316
solve poisson equation for blending images
double operator()(PixelType const &v1, PixelType const &v2) const
PixelType operator()(PixelType const &v1, PixelType const &v2) const
void MergeImages(ImageType &image1, MaskType &mask1, const ImageType &image2, const MaskType &mask2, const vigra::Diff2D offset, const bool wrap, const bool hardSeam)
vigra::pair< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImage(const ROIImage< Image, Mask > &img)
Definition: ROIImage.h:300
double operator()(vigra::RGBValue< PixelType > const &v1, vigra::RGBValue< PixelType > const &v2) const
void combineTwoImages(SrcImageIterator1 src1_upperleft, SrcImageIterator1 src1_lowerright, SrcAccessor1 src1_acc, SrcImageIterator2 src2_upperleft, SrcAccessor2 src2_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc, const Functor &func)
Definition: openmp_vigra.h:249
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
vigra::FRGBImage ImageType
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
ImageType ResizeImage(const ImageType &image, const vigra::Size2D &newSize)
static T max(T x, T y)
Definition: svm.cpp:65
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 copyImage(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc)
Definition: openmp_vigra.h:305
PixelType operator()(PixelType const &v1, PixelType const &v2) const
void BuildSeamPyramid(const Image &input, vigra::ImagePyramid< PyramidImage > &seams, const int minLength)
Definition: BlendPoisson.h:603
void PoissonBlend(ImageType &image1, const ImageType &image2, const MaskType &mask2, const vigra::BImage &labels, const vigra::Point2D &offsetPoint, const bool doWrap)
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