12 namespace SourceXtractor {
15 constexpr
float PI = boost::math::constants::pi<float>();
18 float x, val, sinx1, cosx1;
20 if (pos < threshold && pos > -threshold) {
27 x = -
PI / 2.0 * (pos + 1.0);
28 sincosf(x, &sinx1, &cosx1);
29 val = (*(kernel++) = sinx1 / (x * x));
31 val += (*(kernel++) = -cosx1 / (x * x));
33 val += (*(kernel++) = -sinx1 / (x * x));
35 val += (*kernel = cosx1 / (x *
x));
45 float x, val, sinx1, sinx2, sinx3, cosx1;
47 if (pos < threshold && pos > -threshold) {
56 x = -
PI / 3.0 * (pos + 2.0);
57 sincosf(x, &sinx1, &cosx1);
58 val = (*(kernel++) = sinx1 / (x * x));
60 val += (*(kernel++) = (sinx2 = -0.5 * sinx1 - 0.866025403785 * cosx1) / (x *
x));
62 val += (*(kernel++) = (sinx3 = -0.5 * sinx1 + 0.866025403785 * cosx1) / (x *
x));
64 val += (*(kernel++) = sinx1 / (x * x));
66 val += (*(kernel++) = sinx2 / (x * x));
68 val += (*kernel = sinx3 / (x *
x));
80 float x, val, sinx1, sinx2, sinx3, cosx1;
82 if (pos < threshold && pos > -threshold) {
93 x = -
PI / 4.0 * (pos + 3.0);
94 sincosf(x, &sinx1, &cosx1);
95 val = (*(kernel++) = sinx1 / (x * x));
97 val += (*(kernel++) = -(sinx2 = 0.707106781186 * (sinx1 + cosx1)) / (x * x));
99 val += (*(kernel++) = cosx1 / (x * x));
101 val += (*(kernel++) = -(sinx3 = 0.707106781186 * (cosx1 - sinx1)) / (x * x));
103 val += (*(kernel++) = -sinx1 / (x * x));
105 val += (*(kernel++) = sinx2 / (x * x));
107 val += (*(kernel++) = -cosx1 / (x * x));
109 val += (*kernel = sinx3 / (x *
x));
123 const float threshold = 1
e-6;
125 switch (interptype) {
130 *(kernel++) = 1.0 - pos;
148 int xsize,
int ysize,
interpenum interptype) {
150 static const int interp_kernwidth[5] = {1, 2, 4, 6, 8};
154 *kvector, *pixin, *pixout,
156 int i, j, ix, iy, kwidth, step;
158 kwidth = interp_kernwidth[interptype];
162 ix = (int) (x - 0.50001);
163 iy = (int) (y - 0.50001);
176 if (ix < 0 || ix + kwidth <= 0 || ix + kwidth > xsize ||
177 iy < 0 || iy + kwidth <= 0 || iy + kwidth > ysize)
182 step = xsize - kwidth;
183 pixin = pix + iy * xsize + ix;
185 for (j = kwidth; j--;) {
188 for (i = kwidth; i--;)
189 val += *(kvector++) * *(pixin++);
199 for (i = kwidth; i--;)
200 val += *(kvector++) * *(pixin++);
212 double scale_factor,
double x_shift,
double y_shift) {
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;
223 double x_delta = x - xi;
224 double y_delta = y - yi;
229 double v11 =
getClamped(source, xi + 1, yi + 1);
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);
238 double scale_factor,
double x_shift,
double y_shift) {
242 auto data = &source->getData()[0];
243 auto source_width = source->getWidth();
244 auto source_height = source->getHeight();
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;
262 double scale_factor,
double x_shift,
double y_shift) {
267 auto source_width = source->getWidth();
268 auto source_height = source->getHeight();
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;
279 if (ipos < 0 || ipos + 7 >= source_height) {
284 for (
unsigned int buff_y = 0; buff_y <
Traits::height(buffer); buff_y++) {
286 for (
int i = 0; i < 8; i++) {
287 val += kernel[i] *
Traits::at(source, buff_y, ipos + i);
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) {
302 for (
int y_win = 0; y_win < window_height; y_win++) {
304 for (
int i = 0; i < 8; i++) {
305 val += kernel[i] *
Traits::at(buffer, y_win, ipos + i);
315 namespace ModelFitting {
320 double scale_factor,
double x,
double y) {
322 double scaled_width =
width(source_image) * scale_factor;
323 double scaled_height =
height(source_image) * scale_factor;
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;
332 double x_shift = x - scaled_width / 2. - x_min;
333 double y_shift = y - scaled_height / 2. - y_min;
335 auto window = factory(window_width, window_height);
343 double corr_factor = 1. / (scale_factor * scale_factor);
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);
#define INTERP_MAXKERNELWIDTH
static std::size_t height(const ImageInterfaceTypePtr &image)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
static void addImageToImage(ImageType &image1, const ImageType &image2, double scale, double x, double y)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
static ImageInterfaceType::PixelType & at(ImageInterfaceTypePtr &image, std::size_t x, std::size_t y)
static ImageInterfaceTypePtr factory(std::size_t width, std::size_t height)
static std::size_t width(const ImageInterfaceTypePtr &image)
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > dy