Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
#include "facet.h"
#include <aocommon/imagecoordinates.h>
#include <cassert>
#include <cmath>
#include <stdexcept>
namespace schaapcommon {
namespace facets {
BoundingBox::BoundingBox(const std::vector<Pixel>& pixels, size_t align) {
if (pixels.empty()) {
throw std::invalid_argument("Cannot create boundingbox for 0 pixels");
}
min_ = max_ = pixels.front();
for (auto i = pixels.begin() + 1; i != pixels.end(); ++i) {
min_.x = std::min(min_.x, i->x);
max_.x = std::max(max_.x, i->x);
min_.y = std::min(min_.y, i->y);
max_.y = std::max(max_.y, i->y);
}
if (align > 1) {
min_.x -= min_.x % align;
min_.y -= min_.y % align;
max_.x += align - 1 - (max_.x + align - 1) % align;
max_.y += align - 1 - (max_.y + align - 1) % align;
}
}
std::pair<int, int> Facet::HorizontalIntersections(
const int y_intersect) const {
if (y_intersect < min_y_ || y_intersect >= max_y_) {
return std::make_pair<int, int>(0, 0);
}
bool found_x1 = false;
int x1 = 0;
int x2 = 0;
size_t i;
for (i = 0; i != pixels_.size(); ++i) {
const Pixel& p1 = pixels_[i];
const Pixel& p2 = pixels_[(i + 1) % pixels_.size()];
if ((p1.y <= y_intersect && p2.y > y_intersect) ||
(p2.y <= y_intersect && p1.y > y_intersect)) {
int x;
if (p1.y == y_intersect) {
x = p1.x;
} else if (p2.y == y_intersect) {
x = p2.x;
} else {
const double beta = double(p2.x - p1.x) / double(p2.y - p1.y);
const double xfl = p1.x + beta * (y_intersect - p1.y);
x = round(xfl);
}
if (!found_x1) {
x1 = x;
found_x1 = true;
} else {
x2 = x;
break;
}
}
}
// The loop should have found x1 and x2, and then stopped using 'break'.
assert(i != pixels_.size());
return std::minmax(x1, x2);
}
void Facet::CalculatePixels(double phase_centre_ra, double phase_centre_dec,
double px_scale_x, double px_scale_y,
size_t image_width, size_t image_height,
double shift_l, double shift_m,
bool origin_at_centre, double padding, size_t align,
bool make_square) {
if (padding < 1.0) {
throw std::invalid_argument("Padding factor should be >= 1.0");
}
if ((align > 1) && ((image_width % align) || (image_height % align))) {
throw std::invalid_argument("Image is not aligned");
}
if (coords_.size() < 3) {
throw std::runtime_error("Number of coordinates < 3, facet incomplete!");
}
pixels_.clear();
pixels_.reserve(coords_.size());
for (const Coord& coord : coords_) {
double l, m;
int x, y;
aocommon::ImageCoordinates::RaDecToLM(coord.ra, coord.dec, phase_centre_ra,
phase_centre_dec, l, m);
l -= shift_l;
m -= shift_m;
aocommon::ImageCoordinates::LMToXY(l, m, px_scale_x, px_scale_y,
image_width, image_height, x, y);
// Clip x and y to the image size.
x = std::max(x, 0);
y = std::max(y, 0);
x = std::min(x, static_cast<int>(image_width));
y = std::min(y, static_cast<int>(image_height));
if (origin_at_centre) {
x = static_cast<int>(image_width / 2) - x;
y = y - static_cast<int>(image_height / 2);
}
pixels_.emplace_back(x, y);
}
// Calculate bounding box for the pixels only, and set min_y_ and max_y_.
BoundingBox pixel_box(pixels_);
min_y_ = pixel_box.Min().y;
max_y_ = pixel_box.Max().y;
// Calculate the trimmed_box_.
trimmed_box_ = BoundingBox({pixel_box.Min(), pixel_box.Max()}, align);
if (make_square) {
size_t trimmed_size = std::max(trimmed_box_.Width(), trimmed_box_.Height());
const Pixel square_pad((trimmed_size - trimmed_box_.Width()) / 2,
(trimmed_size - trimmed_box_.Height()) / 2);
trimmed_box_ = BoundingBox(
{trimmed_box_.Min() - square_pad, trimmed_box_.Max() + square_pad});
}
// Calculate the untrimmed_box_.
size_t width = static_cast<size_t>(std::ceil(trimmed_box_.Width() * padding));
size_t height =
static_cast<size_t>(std::ceil(trimmed_box_.Height() * padding));
// Calculate padding. Divide by two since the padding occurs at both sides.
const Pixel pad((width - trimmed_box_.Width()) / 2,
(height - trimmed_box_.Height()) / 2);
// Create the padded and aligned bounding box for the facet.
untrimmed_box_ =
BoundingBox({trimmed_box_.Min() - pad, trimmed_box_.Max() + pad}, align);
}
} // namespace facets
} // namespace schaapcommon