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
|
#include <boost/test/unit_test.hpp>
#include "../../math/renderer.h"
#include "../../model/model.h"
#include "../../model/powerlawsed.h"
#include <aocommon/image.h>
#include <aocommon/threadpool.h>
#include <schaapcommon/fitters/gaussianfitter.h>
#include <schaapcommon/math/restoreimage.h>
namespace wsclean {
namespace {
constexpr size_t kThreadCount = 1;
constexpr size_t kWidth = 64;
constexpr size_t kHeight = 64;
constexpr long double kPixelSize = 1 /*amin*/ * (M_PI / 180.0 / 60.0);
struct RendererFixture {
RendererFixture() : restored(kWidth, kHeight, 0.0) {
imageSettings.ra = 0.0;
imageSettings.dec = 0.0;
imageSettings.pixel_scale_l = kPixelSize;
imageSettings.pixel_scale_m = kPixelSize;
imageSettings.l_shift = 0.0;
imageSettings.m_shift = 0.0;
}
aocommon::Image restored;
renderer::ImageCoordinateSettings imageSettings;
};
} // namespace
BOOST_AUTO_TEST_SUITE(render_and_fit)
BOOST_FIXTURE_TEST_CASE(fit_with_bad_initial_value, RendererFixture) {
PowerLawSED sed(150.0e6, 1.0);
ModelComponent component;
component.SetPosDec(0.0);
component.SetPosRA(0.0);
component.SetSED(sed);
ModelSource source;
source.AddComponent(component);
Model model;
model.AddSource(source);
const long double beamMaj = 4.0L * kPixelSize;
const long double beamMin = 4.0L * kPixelSize;
const long double beamPA = 0.0;
const long double estimatedBeamPx = 1.0; // this is on purpose way off
aocommon::ThreadPool::GetInstance().SetNThreads(kThreadCount);
renderer::RestoreWithEllipticalBeam(restored, imageSettings, model, beamMaj,
beamMin, beamPA, 100e6, 200e6,
aocommon::Polarization::StokesI);
const schaapcommon::math::Ellipse ellipse =
schaapcommon::fitters::Fit2DGaussianCentred(
restored.Data(), restored.Width(), restored.Height(), estimatedBeamPx,
10.0, false);
BOOST_CHECK_CLOSE_FRACTION(ellipse.major, 4.0, 1e-4);
BOOST_CHECK_CLOSE_FRACTION(ellipse.minor, 4.0, 1e-4);
}
BOOST_FIXTURE_TEST_CASE(fit_circular, RendererFixture) {
PowerLawSED sed(150.0e6, 1.0);
ModelComponent component;
component.SetPosDec(0.0);
component.SetPosRA(0.0);
component.SetSED(sed);
ModelSource source;
source.AddComponent(component);
Model model;
model.AddSource(source);
const long double beamMaj = 4.0L * kPixelSize;
const long double beamMin = 4.0L * kPixelSize;
const long double beamPA = 0.0;
const long double estimatedBeamPx = 1.0; // this is on purpose way off
aocommon::ThreadPool::GetInstance().SetNThreads(kThreadCount);
renderer::RestoreWithEllipticalBeam(restored, imageSettings, model, beamMaj,
beamMin, beamPA, 100e6, 200e6,
aocommon::Polarization::StokesI);
double fitMajor = estimatedBeamPx;
schaapcommon::fitters::Fit2DCircularGaussianCentred(
restored.Data(), restored.Width(), restored.Height(), fitMajor);
BOOST_CHECK_CLOSE_FRACTION(fitMajor, 4.0, 1e-4);
}
BOOST_FIXTURE_TEST_CASE(fit_small_beam, RendererFixture) {
PowerLawSED sed(150.0e6, 1.0);
ModelComponent component;
component.SetPosDec(0.0);
component.SetPosRA(0.0);
component.SetSED(sed);
ModelSource source;
source.AddComponent(component);
Model model;
model.AddSource(source);
const long double beamMaj = 4.0L * kPixelSize;
const long double beamMin = 0.5L * kPixelSize;
const long double beamPA = 0.0;
const long double estimatedBeamPx = 1.0; // this is on purpose way off
aocommon::ThreadPool::GetInstance().SetNThreads(kThreadCount);
renderer::RestoreWithEllipticalBeam(restored, imageSettings, model, beamMaj,
beamMin, beamPA, 100e6, 200e6,
aocommon::Polarization::StokesI);
const schaapcommon::math::Ellipse ellipse =
schaapcommon::fitters::Fit2DGaussianCentred(
restored.Data(), restored.Width(), restored.Height(), estimatedBeamPx,
10.0, false);
BOOST_CHECK_CLOSE_FRACTION(ellipse.minor, 0.5, 1e-4);
}
BOOST_AUTO_TEST_SUITE_END()
} // namespace wsclean
|
