Line data Source code
1 : //# Copyright (C) 1995,1997,1998,1999,2000,2001,2002,2003
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU Library General Public License as published by
6 : //# the Free Software Foundation; either version 2 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
12 : //# License for more details.
13 : //#
14 : //# You should have received a copy of the GNU Library General Public License
15 : //# along with this library; if not, write to the Free Software Foundation,
16 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: casa-feedback@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 :
26 : #include <components/ComponentModels/GaussianDeconvolver.h>
27 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
28 :
29 : using namespace casacore;
30 : namespace casa {
31 :
32 0 : CasaImageBeamSet::CasaImageBeamSet() : ImageBeamSet() {}
33 :
34 0 : CasaImageBeamSet::CasaImageBeamSet(const Matrix<GaussianBeam>& beams)
35 0 : : ImageBeamSet(beams) {}
36 :
37 0 : CasaImageBeamSet::CasaImageBeamSet(const GaussianBeam& beam)
38 0 : : ImageBeamSet(beam) {}
39 :
40 0 : CasaImageBeamSet::CasaImageBeamSet(
41 : uInt nchan, uInt nstokes, const GaussianBeam& beam
42 0 : ) : ImageBeamSet(nchan, nstokes, beam) {}
43 :
44 0 : CasaImageBeamSet::CasaImageBeamSet(const CasaImageBeamSet& other)
45 0 : : ImageBeamSet(other) {}
46 :
47 0 : CasaImageBeamSet::CasaImageBeamSet(const ImageBeamSet& other)
48 0 : : ImageBeamSet(other) {}
49 :
50 0 : CasaImageBeamSet::~CasaImageBeamSet() {}
51 :
52 0 : CasaImageBeamSet& CasaImageBeamSet::operator=(const CasaImageBeamSet& other) {
53 0 : ImageBeamSet::operator=(other);
54 0 : return *this;
55 : }
56 :
57 :
58 0 : const String& CasaImageBeamSet::className() {
59 0 : static const String c = "CasaImageBeamSet";
60 0 : return c;
61 : }
62 :
63 0 : GaussianBeam CasaImageBeamSet::getCommonBeam() const {
64 0 : if (empty()) {
65 0 : throw AipsError("This beam set is empty.");
66 : }
67 0 : const Matrix<GaussianBeam> beams = getBeams();
68 0 : if (allTrue(beams == GaussianBeam::NULL_BEAM)) {
69 0 : throw AipsError("All beams are null.");
70 : }
71 0 : if (allTrue(beams == beams(IPosition(2, 0)))) {
72 0 : return beams(IPosition(2, 0));
73 : }
74 0 : BeamIter end = beams.end();
75 0 : Bool largestBeamWorks = true;
76 0 : GaussianBeam junk;
77 0 : GaussianBeam problemBeam;
78 0 : GaussianBeam maxBeam = getMaxAreaBeam();
79 : //Double myMajor = maxBeam.getMajor("arcsec");
80 : // Double myMinor = maxBeam.getMinor("arcsec");
81 :
82 0 : for (
83 0 : BeamIter iBeam = beams.begin();
84 0 : iBeam != end; iBeam++
85 : ) {
86 0 : if (*iBeam != maxBeam && !iBeam->isNull()) {
87 : //myMajor = max(myMajor, iBeam->getMajor("arcsec"));
88 : //myMinor = max(myMinor, iBeam->getMinor("arcsec"));
89 : try {
90 0 : if (GaussianDeconvolver::deconvolve(junk, maxBeam, *iBeam)) {
91 0 : largestBeamWorks = false;
92 0 : problemBeam = *iBeam;
93 : }
94 : }
95 0 : catch (const AipsError& x) {
96 0 : largestBeamWorks = false;
97 0 : problemBeam = *iBeam;
98 0 : }
99 : }
100 0 : }
101 0 : if (largestBeamWorks) {
102 0 : return maxBeam;
103 : }
104 :
105 : // transformation 1, rotate so one of the ellipses' major axis lies
106 : // along the x axis. Ellipse A is _maxBeam, ellipse B is problemBeam,
107 : // ellipse C is our wanted ellipse
108 :
109 0 : Double tB1 = problemBeam.getPA("rad", true) - maxBeam.getPA("rad", true);
110 :
111 0 : if (abs(tB1) == C::pi / 2) {
112 0 : Bool maxHasMajor = maxBeam.getMajor("arcsec")
113 0 : >= problemBeam.getMajor("arcsec");
114 : // handle the situation of right angles explicitly because things blow up otherwise
115 : Quantity major =
116 0 : maxHasMajor ? maxBeam.getMajor() : problemBeam.getMajor();
117 : Quantity minor =
118 0 : maxHasMajor ? problemBeam.getMajor() : maxBeam.getMajor();
119 : Quantity pa =
120 0 : maxHasMajor ? maxBeam.getPA(true) : problemBeam.getPA(true);
121 0 : return GaussianBeam(major, minor, pa);
122 0 : }
123 :
124 0 : Double aA1 = maxBeam.getMajor("arcsec");
125 0 : Double bA1 = maxBeam.getMinor("arcsec");
126 0 : Double aB1 = problemBeam.getMajor("arcsec");
127 0 : Double bB1 = problemBeam.getMinor("arcsec");
128 :
129 : // transformation 2: Squeeze along the x axis and stretch along y axis so
130 : // ellipse A becomes a circle, preserving its area
131 0 : Double aA2 = sqrt(aA1 * bA1);
132 0 : Double bA2 = aA2;
133 0 : Double p = aA2 / aA1;
134 0 : Double q = bA2 / bA1;
135 :
136 : // ellipse B's parameters after transformation 2
137 : Double aB2, bB2, tB2;
138 :
139 0 : _transformEllipseByScaling(aB2, bB2, tB2, aB1, bB1, tB1, p, q);
140 :
141 : // Now the enclosing transformed ellipse, C, has semi-major axis equal to aB2,
142 : // minor axis is aA2 == bA2, and the pa is tB2
143 0 : Double aC2 = aB2;
144 0 : Double bC2 = aA2;
145 0 : Double tC2 = tB2;
146 :
147 : // Now reverse the transformations, first transforming ellipse C by shrinking the coordinate
148 : // system of transformation 2 yaxis and expanding its xaxis to return to transformation 1.
149 :
150 : Double aC1, bC1, tC1;
151 0 : _transformEllipseByScaling(aC1, bC1, tC1, aC2, bC2, tC2, 1 / p, 1 / q);
152 :
153 : // now rotate by _maxBeam.getPA() to get the untransformed enclosing ellipse
154 :
155 0 : Double aC = aC1;
156 0 : Double bC = bC1;
157 0 : Double tC = tC1 + maxBeam.getPA("rad", true);
158 :
159 : // confirm that we can indeed convolve both beams with the enclosing ellipse
160 0 : GaussianBeam newMaxBeam = GaussianBeam(Quantity(aC, "arcsec"),
161 0 : Quantity(bC, "arcsec"), Quantity(tC, "rad"));
162 : // Sometimes (due to precision issues I suspect), the found beam has to be increased slightly
163 : // so our deconvolving method doesn't fail
164 0 : Bool ok = false;
165 0 : while (!ok) {
166 : try {
167 0 : if (GaussianDeconvolver::deconvolve(junk, newMaxBeam, maxBeam)) {
168 0 : throw AipsError();
169 : }
170 0 : if (GaussianDeconvolver::deconvolve(junk, newMaxBeam, problemBeam)) {
171 0 : throw AipsError();
172 : }
173 0 : ok = true;
174 : }
175 0 : catch (const AipsError& x) {
176 : // deconvolution issues, increase the enclosing beam size slightly
177 0 : aC *= 1.001;
178 0 : bC *= 1.001;
179 0 : newMaxBeam = GaussianBeam(Quantity(aC, "arcsec"),
180 0 : Quantity(bC, "arcsec"), Quantity(tC, "rad"));
181 0 : }
182 : }
183 : // create a new beam set to run this method on, replacing _maxBeam with ellipse C
184 :
185 0 : CasaImageBeamSet newBeamSet(*this);
186 0 : Array<GaussianBeam> newBeams = beams.copy();
187 0 : newBeams(getMaxAreaBeamPosition()) = newMaxBeam;
188 0 : newBeamSet.setBeams(newBeams);
189 :
190 0 : return newBeamSet.getCommonBeam();
191 0 : }
192 :
193 0 : void CasaImageBeamSet::_transformEllipseByScaling(
194 : Double& transformedMajor,
195 : Double& transformedMinor, Double& transformedPA, Double major,
196 : Double minor, Double pa, Double xScaleFactor, Double yScaleFactor
197 : ) {
198 0 : Double mycos = cos(pa);
199 0 : Double mysin = sin(pa);
200 0 : Double cos2 = mycos * mycos;
201 0 : Double sin2 = mysin * mysin;
202 0 : Double major2 = major * major;
203 0 : Double minor2 = minor * minor;
204 0 : Double a = cos2 / (major2) + sin2 / (minor2);
205 0 : Double b = -2 * mycos * mysin * (1 / (major2) - 1 / (minor2));
206 0 : Double c = sin2 / (major2) + cos2 / (minor2);
207 :
208 0 : Double xs = xScaleFactor * xScaleFactor;
209 0 : Double ys = yScaleFactor * yScaleFactor;
210 :
211 0 : Double r = a / xs;
212 0 : Double s = b * b / (4 * xs * ys);
213 0 : Double t = c / ys;
214 :
215 0 : Double u = r - t;
216 0 : Double u2 = u * u;
217 :
218 0 : Double f1 = u2 + 4 * s;
219 0 : Double f2 = sqrt(f1) * abs(u);
220 :
221 0 : Double j1 = (f2 + f1) / f1 / 2;
222 0 : Double j2 = (-f2 + f1) / f1 / 2;
223 :
224 0 : Double k1 = (j1 * r + j1 * t - t) / (2 * j1 - 1);
225 0 : Double k2 = (j2 * r + j2 * t - t) / (2 * j2 - 1);
226 :
227 0 : Double c1 = sqrt(1 / k1);
228 0 : Double c2 = sqrt(1 / k2);
229 :
230 0 : if (c1 == c2) {
231 : // the transformed ellipse is a circle
232 0 : transformedMajor = sqrt(k1);
233 0 : transformedMinor = transformedMajor;
234 0 : transformedPA = 0;
235 0 : } else if (c1 > c2) {
236 : // c1 is the major axis and so j1 is the solution for the corresponding pa
237 : // of the transformed ellipse
238 0 : transformedMajor = c1;
239 0 : transformedMinor = c2;
240 0 : transformedPA = (pa >= 0 ? 1 : -1) * acos(sqrt(j1));
241 : } else {
242 : // c2 is the major axis and so j2 is the solution for the corresponding pa
243 : // of the transformed ellipse
244 0 : transformedMajor = c2;
245 0 : transformedMinor = c1;
246 0 : transformedPA = (pa >= 0 ? 1 : -1) * acos(sqrt(j2));
247 : }
248 0 : }
249 :
250 : }
|
