SourceXtractorPlusPlus  0.11
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
OpenCvMatImageTraits.h
Go to the documentation of this file.
1 
23 #ifndef MODELFITTING_OPENCVMATIMAGETRAITS_H
24 #define MODELFITTING_OPENCVMATIMAGETRAITS_H
25 
26 #define INTERP_MAXKERNELWIDTH 8 // Max. range of kernel (pixels)
27 
28 #include <boost/math/constants/constants.hpp>
29 #include <cmath>
30 #include <algorithm>
31 #include <chrono>
32 #include <utility>
33 #include <opencv2/opencv.hpp>
35 
36 // Interpolation types
39 
40 namespace ModelFitting {
41 
42 template <>
43 struct ImageTraits<cv::Mat> {
44 
45  using iterator = decltype(std::declval<cv::Mat>().begin<double>());
46 
48  return cv::Mat::zeros((int)height, (int)width, CV_64F);
49  }
50 
51  static std::size_t width(const cv::Mat& image) {
52  return image.cols;
53  }
54 
55  static std::size_t height(const cv::Mat& image) {
56  return image.rows;
57  }
58 
59  static double& at(cv::Mat& image, std::size_t x, std::size_t y) {
60  return image.at<double>((int)y, (int)x);
61  }
62 
63  static double at(const cv::Mat& image, std::size_t x, std::size_t y) {
64  return image.at<double>((int)y, (int)x);
65  }
66 
67  static iterator begin(cv::Mat& image) {
68  return image.begin<double>();
69  }
70 
71  static iterator end(cv::Mat& image) {
72  return image.end<double>();
73  }
74 
75  static void addImageToImage(cv::Mat& image1, const cv::Mat& image2,
76  double scale_factor, double x, double y) {
77  // Calculate the size in pixels of the image2 after in the scale of image1
78  double scaled_width = width(image2) * scale_factor;
79  double scaled_height = height(image2) * scale_factor;
80  // Calculate the window of the image1 which is affected
81  int x_min = std::floor(x - scaled_width / 2.);
82  int x_max = std::ceil(x + scaled_width / 2.);
83  int window_width = x_max - x_min;
84  int y_min = std::floor(y - scaled_height / 2.);
85  int y_max= std::ceil(y + scaled_height / 2.);
86  int window_height = y_max - y_min;
87  // Calculate the shift of the image2 inside the window
88  double x_shift = x - scaled_width / 2. - x_min;
89  double y_shift = y - scaled_height / 2. - y_min;
90  // Create the scaled and shifted window
91  auto window = factory(window_width, window_height);
92 
93 // cv::Mat affine_trans = (cv::Mat_<double>(2, 3) <<
94 // scale_factor, 0., x_shift,
95 // 0., scale_factor, y_shift
96 // );
97  //cv::warpAffine(image2, window, affine_trans, {window_width, window_height});
98  //cv::warpAffine(image2, window, affine_trans, {window_width, window_height}, cv::INTER_LANCZOS4);
99 
100  shiftResize(image2, window, scale_factor, x_shift, y_shift);
101 
102  // We need to correct the window for the scaling, so it has the same integral
103  // with the image2
104  double corr_factor = 1. / (scale_factor * scale_factor);
105  window = corr_factor * window;
106  // Add the window to the image1
107  for(int x_im=std::max(x_min,0); x_im<std::min(x_max,image1.cols); ++x_im) {
108  for (int y_im=std::max(y_min,0); y_im<std::min(y_max,image1.rows); ++y_im) {
109  int x_win = x_im - x_min;
110  int y_win = y_im - y_min;
111  at(image1, x_im, y_im) += at(window, x_win, y_win);
112  }
113  }
114  }
115 
116  static double getClamped(const cv::Mat& image, int x, int y) {
117  return at(image, std::max(0, std::min(x, (int) width(image))), std::max(0, std::min(y, (int) height(image))));
118  }
119 
120  static void shiftResize(const cv::Mat& source, cv::Mat& window, double scale_factor, double x_shift, double y_shift) {
121  int window_width = width(window);
122  int window_height = height(window);
123  for(int x_win=0; x_win < window_width; x_win++) {
124  for(int y_win=0; y_win < window_height; y_win++) {
125  double x = (x_win - 0.5 - x_shift) / scale_factor + 0.5;
126  double y = (y_win - 0.5 - y_shift) / scale_factor + 0.5;
127 
128  int xi = std::floor(x);
129  int yi = std::floor(y);
130 
131  double x_delta = x - xi;
132  double y_delta = y - yi;
133 
134  double v00 = getClamped(source, xi, yi);
135  double v01 = getClamped(source, xi, yi+1);
136  double v10 = getClamped(source, xi+1, yi);
137  double v11 = getClamped(source, xi+1, yi+1);
138 
139  at(window, x_win, y_win) = (1.0 - y_delta) * ((1.0 - x_delta) * v00 + x_delta * v10) +
140  y_delta * ((1.0 - x_delta) * v01 + x_delta * v11);
141  }
142  }
143 
144  }
145 
146  /*
147  static void shiftResizeLancszos(const cv::Mat& source, cv::Mat& window, double scale_factor, double x_shift, double y_shift) {
148  int window_width = width(window);
149  int window_height = height(window);
150  for(int x_win=0; x_win < window_width; x_win++) {
151  for(int y_win=0; y_win < window_height; y_win++) {
152  float x = (x_win - x_shift) / scale_factor;
153  float y = (y_win - y_shift) / scale_factor;
154 
155  //at(window, x_win, y_win) = interpolate_pix(&source.getData()[0], x, y, sourcexsize, sourceysize);
156  }
157  }
158  }
159 */
160 
161  static float interpolate_pix(float *pix, float x, float y,
162  int xsize, int ysize, interpenum interptype) {
163 
164  static const int interp_kernwidth[5]={1,2,4,6,8};
165 
166  float buffer[INTERP_MAXKERNELWIDTH],
167  kernel[INTERP_MAXKERNELWIDTH],
168  *kvector, *pixin, *pixout,
169  dx, dy, val;
170  int i, j, ix, iy, kwidth, step;
171 
172  kwidth = interp_kernwidth[interptype];
173 
174 //-- Get the integer part of the current coordinate or nearest neighbour
175  if (interptype == INTERP_NEARESTNEIGHBOUR) {
176  ix = (int)(x-0.50001);
177  iy = (int)(y-0.50001);
178  } else {
179  ix = (int)x;
180  iy = (int)y;
181  }
182 
183 //-- Store the fractional part of the current coordinate
184  dx = x - ix;
185  dy = y - iy;
186 //-- Check if interpolation start/end exceed image boundary
187  ix -= kwidth / 2;
188  iy -= kwidth / 2;
189  if (ix < 0 || ix + kwidth <= 0 || ix + kwidth > xsize ||
190  iy < 0 || iy + kwidth <= 0 || iy + kwidth > ysize)
191  return 0.0;
192 
193 //-- First step: interpolate along NAXIS1 from the data themselves
194  make_kernel(dx, kernel, interptype);
195  step = xsize - kwidth;
196  pixin = pix + iy * xsize + ix ; // Set starting pointer
197  pixout = buffer;
198  for (j = kwidth; j--;) {
199  val = 0.0;
200  kvector = kernel;
201  for (i = kwidth; i--;)
202  val += *(kvector++) * *(pixin++);
203  *(pixout++) = val;
204  pixin += step;
205  }
206 
207 //-- Second step: interpolate along NAXIS2 from the interpolation buffer
208  make_kernel(dy, kernel, interptype);
209  pixin = buffer;
210  val = 0.0;
211  kvector = kernel;
212  for (i = kwidth; i--;)
213  val += *(kvector++) * *(pixin++);
214 
215  return val;
216  }
217 
218 
219  static void make_kernel(float pos, float *kernel, interpenum interptype) {
220  const float pi = boost::math::constants::pi<float>();
221  const float threshold = 1e-6;
222  float x, val, sinx1,sinx2,sinx3,cosx1;
223 
224  if (interptype == INTERP_NEARESTNEIGHBOUR)
225  *kernel = 1.0;
226  else if (interptype == INTERP_BILINEAR) {
227  *(kernel++) = 1.0 - pos;
228  *kernel = pos;
229  } else if (interptype == INTERP_LANCZOS2) {
230  if (pos < threshold && pos > - threshold) {
231  *(kernel++) = 0.0;
232  *(kernel++) = 1.0;
233  *(kernel++) = 0.0;
234  *kernel = 0.0;
235  } else {
236  x = - pi / 2.0 * (pos + 1.0);
237  sincosf(x, &sinx1, &cosx1);
238  val = (*(kernel++) = sinx1 / (x * x));
239  x += pi / 2.0;
240  val += (*(kernel++) = -cosx1 / (x * x));
241  x += pi / 2.0;
242  val += (*(kernel++) = -sinx1 / (x * x));
243  x += pi / 2.0;
244  val += (*kernel = cosx1 / (x * x));
245  val = 1.0/val;
246  *(kernel--) *= val;
247  *(kernel--) *= val;
248  *(kernel--) *= val;
249  *kernel *= val;
250  }
251  } else if (interptype == INTERP_LANCZOS3) {
252  if (pos < threshold && pos > - threshold) {
253  *(kernel++) = 0.0;
254  *(kernel++) = 0.0;
255  *(kernel++) = 1.0;
256  *(kernel++) = 0.0;
257  *(kernel++) = 0.0;
258  *kernel = 0.0;
259  } else {
260  x = - pi / 3.0 * (pos + 2.0);
261  sincosf(x, &sinx1, &cosx1);
262  val = (*(kernel++) = sinx1 / (x * x));
263  x += pi / 3.0;
264  val += (*(kernel++) = (sinx2 = -0.5 * sinx1 - 0.866025403785 * cosx1)
265  / (x*x));
266  x += pi / 3.0;
267  val += (*(kernel++) = (sinx3 = - 0.5 * sinx1 + 0.866025403785 * cosx1)
268  / (x * x));
269  x += pi / 3.0;
270  val += (*(kernel++) = sinx1 / (x * x));
271  x += pi / 3.0;
272  val += (*(kernel++) = sinx2 / (x * x));
273  x += pi / 3.0;
274  val += (*kernel = sinx3 / (x * x));
275  val = 1.0 / val;
276  *(kernel--) *= val;
277  *(kernel--) *= val;
278  *(kernel--) *= val;
279  *(kernel--) *= val;
280  *(kernel--) *= val;
281  *kernel *= val;
282  }
283  } else if (interptype == INTERP_LANCZOS4) {
284  if (pos < threshold && pos > - threshold) {
285  *(kernel++) = 0.0;
286  *(kernel++) = 0.0;
287  *(kernel++) = 0.0;
288  *(kernel++) = 1.0;
289  *(kernel++) = 0.0;
290  *(kernel++) = 0.0;
291  *(kernel++) = 0.0;
292  *kernel = 0.0;
293  } else {
294  x = - pi / 4.0 * (pos + 3.0);
295  sincosf(x, &sinx1, &cosx1);
296  val = (*(kernel++) = sinx1 / (x * x));
297  x += pi / 4.0;
298  val += (*(kernel++) = - (sinx2 = 0.707106781186 * (sinx1 + cosx1))
299  / (x * x));
300  x += pi / 4.0;
301  val += (*(kernel++) = cosx1 / (x * x));
302  x += pi / 4.0;
303  val += (*(kernel++) = - (sinx3 = 0.707106781186 * (cosx1 - sinx1))
304  /(x * x));
305  x += pi / 4.0;
306  val += (*(kernel++) = -sinx1 / (x * x));
307  x += pi / 4.0;
308  val += (*(kernel++) = sinx2 / (x * x));
309  x += pi / 4.0;
310  val += (*(kernel++) = -cosx1 / (x * x));
311  x += pi / 4.0;
312  val += (*kernel = sinx3 / (x * x));
313  val = 1.0 / val;
314  *(kernel--) *= val;
315  *(kernel--) *= val;
316  *(kernel--) *= val;
317  *(kernel--) *= val;
318  *(kernel--) *= val;
319  *(kernel--) *= val;
320  *(kernel--) *= val;
321  *kernel *= val;
322  }
323  }
324 
325  }
326 
327 }; // end of class ImageTraits<cv::Mat>
328 
329 } // end of namespace ModelFitting
330 
331 #endif /* MODELFITTING_OPENCVMATIMAGETRAITS_H */
332 
constexpr double pi
double getClamped(const ImageInterfaceTypePtr &image, int x, int y)
#define INTERP_MAXKERNELWIDTH
static ImageType factory(std::size_t width, std::size_t height)
static void addImageToImage(cv::Mat &image1, const cv::Mat &image2, double scale_factor, double x, double y)
void shiftResize(const ImageInterfaceTypePtr &source, ImageInterfaceTypePtr &window, double scale_factor, double x_shift, double y_shift)
T ceil(T...args)
static std::size_t height(ImageType &image)
static void make_kernel(float pos, float *kernel, interpenum interptype)
static double & at(ImageType &image, std::size_t x, std::size_t y)
constexpr double e
decltype(std::declval< cv::Mat >().begin< double >()) iterator
T floor(T...args)
static double & at(cv::Mat &image, std::size_t x, std::size_t y)
static std::size_t height(const cv::Mat &image)
T min(T...args)
void make_kernel(float pos, float *kernel, interpenum interptype)
static float interpolate_pix(float *pix, float x, float y, int xsize, int ysize, interpenum interptype)
static iterator end(cv::Mat &image)
T max(T...args)
static iterator begin(cv::Mat &image)
static cv::Mat factory(std::size_t width, std::size_t height)
static double at(const cv::Mat &image, std::size_t x, std::size_t y)
static void shiftResize(const cv::Mat &source, cv::Mat &window, double scale_factor, double x_shift, double y_shift)
std::shared_ptr< EngineParameter > dx
static std::size_t width(ImageType &image)
static std::size_t width(const cv::Mat &image)
std::shared_ptr< EngineParameter > dy
static double getClamped(const cv::Mat &image, int x, int y)