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