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