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