Line data Source code
1 : //# Copyright (C) 1993,1994,1995,1996,1999,2001
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: aips2-request@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 <casacore/coordinates/Coordinates/DirectionCoordinate.h>
27 : #include <components/ComponentModels/SkyComponentFactory.h>
28 : #include <casacore/casa/Quanta/MVAngle.h>
29 : #include <components/ComponentModels/GaussianDeconvolver.h>
30 : #include <components/ComponentModels/GaussianShape.h>
31 : #include <components/ComponentModels/ComponentType.h>
32 : #include <casacore/images/Images/ImageUtilities.h>
33 :
34 : using namespace casacore;
35 : namespace casa {
36 :
37 393 : SkyComponent SkyComponentFactory::encodeSkyComponent(
38 : LogIO& logIO, Double& facToJy,
39 : const CoordinateSystem& cSys, const Unit& brightnessUnit,
40 : ComponentType::Shape type, const Vector<Double>& parameters,
41 : Stokes::StokesTypes stokes, Bool xIsLong, const GaussianBeam& beam
42 : ) {
43 : // Input:
44 : // pars(0) = FLux image units (e.g. peak flux in Jy/beam)
45 : // pars(1) = x cen abs pix
46 : // pars(2) = y cen abs pix
47 : // pars(3) = major pix
48 : // pars(4) = minor pix
49 : // pars(5) = pa radians (pos +x -> +y)
50 :
51 393 : SkyComponent sky;
52 :
53 : // Account for the fact that 'x' could be longitude or latitude. Urk.
54 :
55 786 : Vector<Double> pars = parameters.copy();
56 393 : if (!xIsLong) {
57 0 : Double tmp = pars(0);
58 :
59 0 : pars(0) = pars(1);
60 0 : pars(1) = tmp;
61 :
62 0 : Double pa0 = pars(5);
63 0 : MVAngle pa(pa0 + C::pi_2);
64 0 : pa(); // +/- pi
65 0 : pars(5) = pa.radian();
66 : }
67 786 : GaussianBeam cbeam = beam;
68 393 : if (brightnessUnit.getName().contains("beam") && beam.isNull()) {
69 0 : cbeam = ImageUtilities::makeFakeBeam(logIO, cSys);
70 : }
71 :
72 393 : sky.fromPixel(facToJy, pars, brightnessUnit, cbeam, cSys, type, stokes);
73 786 : return sky;
74 : }
75 :
76 : /*
77 : // moved from ImageAnalysis. See comments in ImageUtilities.h
78 : // TODO the only thing that uses this is ImageFitter. So move it there
79 : SkyComponent SkyComponentFactory::encodeSkyComponent(
80 : LogIO& os, Double& facToJy,
81 : const ImageInterface<Float>& subIm, ComponentType::Shape model,
82 : const Vector<Double>& parameters, Stokes::StokesTypes stokes,
83 : Bool xIsLong, Bool deconvolveIt, const GaussianBeam& beam
84 : ) {
85 : //
86 : // This function takes a vector of doubles and converts them to
87 : // a SkyComponent. These doubles are in the 'x' and 'y' frames
88 : // (e.g. result from Fit2D). It is possible that the
89 : // x and y axes of the pixel array are lat/long rather than
90 : // long/lat if the CoordinateSystem has been reordered. So we have
91 : // to take this into account before making the SkyComponent as it
92 : // needs to know long/lat values. The subImage holds only the sky
93 :
94 : // Input
95 : // pars(0) = Flux image units
96 : // pars(1) = x cen abs pix
97 : // pars(2) = y cen abs pix
98 : // pars(3) = major pix
99 : // pars(4) = minor pix
100 : // pars(5) = pa radians (pos +x -> +y)
101 : // Output
102 : // facToJy = converts brightness units to Jy
103 : //
104 :
105 : const CoordinateSystem& cSys = subIm.coordinates();
106 : const Unit& bU = subIm.units();
107 : SkyComponent sky = SkyComponentFactory::encodeSkyComponent(
108 : os, facToJy, cSys, bU, model,
109 : parameters, stokes, xIsLong, beam
110 : );
111 : if (!deconvolveIt) {
112 : return sky;
113 : }
114 :
115 : if (beam.isNull()) {
116 : os << LogIO::WARN
117 : << "This image does not have a restoring beam so no deconvolution possible"
118 : << LogIO::POST;
119 : return sky;
120 : }
121 : Int dirCoordinate = cSys.findCoordinate(Coordinate::DIRECTION);
122 : if (dirCoordinate == -1) {
123 : os << LogIO::WARN
124 : << "This image does not have a DirectionCoordinate so no deconvolution possible"
125 : << LogIO::POST;
126 : return sky;
127 : }
128 : return SkyComponentFactory::deconvolveSkyComponent(os, sky, beam);
129 : }
130 : */
131 :
132 : // moved from ImageAnalysis. See comments in ImageUtilities.h
133 0 : SkyComponent SkyComponentFactory::deconvolveSkyComponent(
134 : LogIO& os,
135 : const SkyComponent& skyIn, const GaussianBeam& beam
136 : ) {
137 0 : const ComponentShape& shapeIn = skyIn.shape();
138 0 : ComponentType::Shape type = shapeIn.type();
139 0 : if (type == ComponentType::POINT) {
140 0 : return skyIn;
141 : }
142 0 : SkyComponent skyOut = skyIn.copy();
143 0 : if (type == ComponentType::GAUSSIAN) {
144 : // Recover shape
145 0 : const TwoSidedShape& ts = dynamic_cast<const TwoSidedShape&> (shapeIn);
146 0 : Quantity major = ts.majorAxis();
147 0 : Quantity minor = ts.minorAxis();
148 0 : Quantity pa = ts.positionAngle();
149 0 : Angular2DGaussian source(major, minor, pa);
150 0 : Angular2DGaussian deconvolvedSize;
151 0 : GaussianDeconvolver::deconvolve(deconvolvedSize, source, beam);
152 0 : const MDirection dirRefIn = shapeIn.refDirection();
153 : GaussianShape shapeOut(
154 0 : dirRefIn, deconvolvedSize.getMajor(),
155 0 : deconvolvedSize.getMinor(),
156 0 : deconvolvedSize.getPA(true)
157 0 : );
158 0 : skyOut.setShape(shapeOut);
159 : }
160 : else {
161 0 : os << "Cannot deconvolve components of type " << shapeIn.ident()
162 0 : << LogIO::EXCEPTION;
163 : }
164 0 : return skyOut;
165 : }
166 :
167 58 : Vector<Double> SkyComponentFactory::decodeSkyComponent (
168 : const SkyComponent& sky,
169 : const ImageInfo& ii,
170 : const CoordinateSystem& cSys,
171 : const Unit& brightnessUnit,
172 : Stokes::StokesTypes stokes,
173 : Bool xIsLong
174 : ) {
175 : //
176 : // The decomposition of the SkyComponent gives things as longitide
177 : // and latitude. But it is possible that the x and y axes of the
178 : // pixel array are lat/long rather than long/lat if the CoordinateSystem
179 : // has been reordered. So we have to take this into account.
180 : //
181 : // Output:
182 : // pars(0) = FLux image units (e.g. peak flux in Jy/beam)
183 : // pars(1) = x cen abs pix
184 : // pars(2) = y cen abs pix
185 : // pars(3) = major pix
186 : // pars(4) = minor pix
187 : // pars(5) = pa radians (pos +x -> +y)
188 : //
189 116 : GaussianBeam beam = ii.restoringBeam();
190 :
191 : // pars(1,2) = longitude, latitude centre
192 :
193 116 : Vector<Double> pars = sky.toPixel (brightnessUnit, beam, cSys, stokes).copy();
194 :
195 : // Now account for the fact that 'x' (horizontally displayed axis) could be
196 : // longitude or latitude. Urk.
197 :
198 58 : Double pa0 = pars(5);
199 58 : if (!xIsLong) {
200 0 : Double tmp = pars(0);
201 0 : pars(0) = pars(1);
202 0 : pars(1) = tmp;
203 : //
204 0 : MVAngle pa(pa0 - C::pi_2);
205 0 : pa(); // +/- pi
206 0 : pa0 = pa.radian();
207 : }
208 58 : pars(5) = pa0;
209 : //
210 116 : return pars;
211 : }
212 :
213 681 : void SkyComponentFactory::worldWidthsToPixel(
214 : Vector<Double>& dParameters,
215 : const Vector<Quantum<Double> >& wParameters,
216 : const CoordinateSystem& cSys,
217 : const IPosition& pixelAxes,
218 : Bool doRef
219 : )
220 : //
221 : // world parameters: x, y, major, minor, pa
222 : // pixel parameters: major, minor, pa (rad)
223 : //
224 : {
225 681 : ThrowIf(
226 : pixelAxes.nelements()!=2,
227 : "You must give two pixel axes"
228 : );
229 681 : ThrowIf(
230 : wParameters.nelements() != 5,
231 : "The world parameters vector must be of length 5."
232 : );
233 :
234 681 : dParameters.resize(3);
235 : Int c0, c1, axisInCoordinate0, axisInCoordinate1;
236 681 : cSys.findPixelAxis(c0, axisInCoordinate0, pixelAxes(0));
237 681 : cSys.findPixelAxis(c1, axisInCoordinate1, pixelAxes(1));
238 :
239 : // Find units
240 :
241 682 : String majorUnit = wParameters(2).getFullUnit().getName();
242 682 : String minorUnit = wParameters(3).getFullUnit().getName();
243 :
244 : // This saves me trying to handle mixed pixel/world units which is a pain for coupled coordinates
245 :
246 681 : ThrowIf(
247 : (majorUnit==String("pix") && minorUnit!=String("pix"))
248 : || (majorUnit!=String("pix") && minorUnit==String("pix")),
249 : "If pixel units are used, both major and minor axes must have pixel units"
250 : );
251 :
252 : // Some checks
253 :
254 681 : Coordinate::Type type0 = cSys.type(c0);
255 681 : Coordinate::Type type1 = cSys.type(c1);
256 681 : ThrowIf(
257 : type0 != type1
258 : && (majorUnit!=String("pix") || minorUnit!=String("pix")),
259 : "The coordinate types for the convolution axes are different. "
260 : "Therefore the units of the major and minor axes of "
261 : "the convolution kernel widths must both be pixels."
262 : );
263 681 : ThrowIf(
264 : type0 == Coordinate::DIRECTION && type1 == Coordinate::DIRECTION && c0 != c1,
265 : "The given axes do not come from the same Direction coordinate. "
266 : "This situation requires further code development."
267 : );
268 681 : ThrowIf(
269 : type0 == Coordinate::STOKES || type1 == Coordinate::STOKES,
270 : "Cannot convolve Stokes axes."
271 : );
272 :
273 : // Deal with pixel units separately. Both are in pixels if either is in pixels.
274 :
275 681 : if (majorUnit==String("pix")) {
276 1 : dParameters(0) = max(wParameters(2).getValue(), wParameters(3).getValue());
277 1 : dParameters(1) = min(wParameters(2).getValue(), wParameters(3).getValue());
278 :
279 1 : if (type0==Coordinate::DIRECTION && type1==Coordinate::DIRECTION) {
280 1 : const DirectionCoordinate& dCoord = cSys.directionCoordinate (c0);
281 :
282 : // Use GaussianShape to get the position angle right. Use the specified
283 : // direction or the reference direction
284 :
285 2 : MDirection world;
286 1 : if (doRef) {
287 0 : dCoord.toWorld(world, dCoord.referencePixel());
288 : }
289 : else {
290 1 : world = MDirection(wParameters(0), wParameters(1), dCoord.directionType());
291 : }
292 :
293 2 : Quantity tmpMaj(1.0, Unit("arcsec"));
294 1 : GaussianShape gaussShape(world, tmpMaj, dParameters(1)/dParameters(0),
295 3 : wParameters(4)); // pa is N->E
296 1 : Vector<Double> pars = gaussShape.toPixel (dCoord);
297 2 : dParameters(2) = pars(4); // pa: +x -> +y
298 : }
299 : else {
300 :
301 : // Some 'mixed' plane; the pa is already +x -> +y
302 :
303 0 : dParameters(2) = wParameters(4).getValue(Unit("rad")); // pa
304 : }
305 1 : return;
306 : }
307 :
308 : // Continue on if non-pixel units
309 :
310 680 : if (type0==Coordinate::DIRECTION && type1==Coordinate::DIRECTION) {
311 :
312 : // Check units are angular
313 :
314 1360 : Unit rad(String("rad"));
315 680 : ThrowIf(
316 : ! wParameters(2).check(rad.getValue()),
317 : "The units of the major axis must be angular"
318 : );
319 680 : ThrowIf(
320 : ! wParameters(3).check(rad.getValue()),
321 : "The units of the minor axis must be angular"
322 : );
323 :
324 : // Make a Gaussian shape to convert to pixels at specified location
325 :
326 680 : const DirectionCoordinate& dCoord = cSys.directionCoordinate (c0);
327 :
328 1360 : MDirection world;
329 680 : if (doRef) {
330 0 : dCoord.toWorld(world, dCoord.referencePixel());
331 : }
332 : else {
333 680 : world = MDirection(wParameters(0), wParameters(1), dCoord.directionType());
334 : }
335 1359 : GaussianShape gaussShape(world, wParameters(2), wParameters(3), wParameters(4));
336 679 : Vector<Double> pars = gaussShape.toPixel (dCoord);
337 679 : dParameters(0) = pars(2);
338 679 : dParameters(1) = pars(3);
339 1358 : dParameters(2) = pars(4); // radians; +x -> +y
340 : }
341 : else {
342 :
343 : // The only other coordinates currently available are non-coupled
344 : // ones and linear except for Tabular, which can be non-regular.
345 : // Urk.
346 :
347 : // Find major and minor axes in pixels
348 :
349 0 : dParameters(0) = _worldWidthToPixel (dParameters(2), wParameters(2),
350 : cSys, pixelAxes);
351 0 : dParameters(1) = _worldWidthToPixel (dParameters(2), wParameters(3),
352 : cSys, pixelAxes);
353 0 : dParameters(2) = wParameters(4).getValue(Unit("rad")); // radians; +x -> +y
354 : }
355 :
356 : // Make sure major > minor
357 :
358 679 : Double tmp = dParameters(0);
359 679 : dParameters(0) = max(tmp, dParameters(1));
360 679 : dParameters(1) = min(tmp, dParameters(1));
361 : }
362 :
363 133 : Bool SkyComponentFactory::pixelWidthsToWorld(
364 : GaussianBeam& wParameters,
365 : const Vector<Double>& pParameters, const CoordinateSystem& cSys,
366 : const IPosition& pixelAxes, Bool doRef
367 : ) {
368 : // pixel parameters: x, y, major, minor, pa (rad)
369 : // world parameters: major, minor, pa
370 133 : ThrowIf(
371 : pixelAxes.nelements() != 2,
372 : "You must give two pixel axes"
373 : );
374 133 : ThrowIf(
375 : pParameters.nelements() != 5,
376 : "The parameters vector must be of length 5"
377 : );
378 : Int c0, axis0, c1, axis1;
379 133 : cSys.findPixelAxis(c0, axis0, pixelAxes(0));
380 133 : cSys.findPixelAxis(c1, axis1, pixelAxes(1));
381 133 : Bool flipped = false;
382 133 : if (
383 133 : cSys.type(c1) == Coordinate::DIRECTION
384 133 : && cSys.type(c0) == Coordinate::DIRECTION
385 : ) {
386 133 : ThrowIf(
387 : c0 != c1,
388 : "Cannot handle axes from different DirectionCoordinates"
389 : );
390 133 : flipped = _skyPixelWidthsToWorld(wParameters, cSys, pParameters, pixelAxes, doRef);
391 : }
392 : else {
393 0 : wParameters = GaussianBeam();
394 : // Major/minor
395 : Quantity q0 = _pixelWidthToWorld(
396 0 : pParameters(4), pParameters(2),
397 : cSys, pixelAxes
398 0 : );
399 : Quantity q1 = _pixelWidthToWorld(
400 0 : pParameters(4), pParameters(3),
401 : cSys, pixelAxes
402 0 : );
403 : // Position angle; radians; +x -> +y
404 0 : if (q0.getValue() < q1.getValue(q0.getFullUnit())) {
405 0 : flipped = true;
406 0 : wParameters = GaussianBeam(q1, q0, Quantity(pParameters(4), "rad"));
407 :
408 : }
409 : else {
410 0 : wParameters = GaussianBeam(q0, q1, Quantity(pParameters(4), "rad"));
411 : }
412 : }
413 133 : return flipped;
414 : }
415 :
416 :
417 133 : Bool SkyComponentFactory::_skyPixelWidthsToWorld (
418 : Angular2DGaussian& gauss2d,
419 : const CoordinateSystem& cSys,
420 : const Vector<Double>& pParameters,
421 : const IPosition& pixelAxes, Bool doRef
422 : )
423 : //
424 : // pixel parameters: x, y, major, minor, pa (rad)
425 : // world parameters: major, minor, pa
426 : //
427 : {
428 : // What coordinates are these axes ?
429 :
430 : Int c0, c1, axisInCoordinate0, axisInCoordinate1;
431 133 : cSys.findPixelAxis(c0, axisInCoordinate0, pixelAxes(0));
432 133 : cSys.findPixelAxis(c1, axisInCoordinate1, pixelAxes(1));
433 : // See what sort of coordinates we have. Make sure it is called
434 : // only for the Sky. More development needed otherwise.
435 :
436 133 : Coordinate::Type type0 = cSys.type(c0);
437 133 : Coordinate::Type type1 = cSys.type(c1);
438 133 : ThrowIf(
439 : type0!=Coordinate::DIRECTION || type1!=Coordinate::DIRECTION,
440 : "Can only be called for axes holding the sky"
441 : );
442 133 : ThrowIf(
443 : c0 != c1,
444 : "The given axes do not come from the same Direction coordinate. "
445 : "This situation requires further code development"
446 : );
447 : // Is the 'x' (first axis) the Longitude or Latitude ?
448 :
449 266 : Vector<Int> dirPixelAxes = cSys.pixelAxes(c0);
450 133 : Bool xIsLong = dirPixelAxes(0)==pixelAxes(0) && dirPixelAxes(1)==pixelAxes(1);
451 133 : uInt whereIsX = 0;
452 133 : uInt whereIsY = 1;
453 133 : if (!xIsLong) {
454 0 : whereIsX = 1;
455 0 : whereIsY = 0;
456 : }
457 : // Encode a pretend GaussianShape from these values as a means
458 : // of converting to world.
459 :
460 133 : const DirectionCoordinate& dCoord = cSys.directionCoordinate(c0);
461 266 : GaussianShape gaussShape;
462 133 : Vector<Double> cParameters(pParameters.copy());
463 133 : if (doRef) {
464 0 : cParameters(0) = dCoord.referencePixel()(whereIsX); // x centre
465 0 : cParameters(1) = dCoord.referencePixel()(whereIsY); // y centre
466 : }
467 : else {
468 133 : if (xIsLong) {
469 133 : cParameters(0) = pParameters(0);
470 133 : cParameters(1) = pParameters(1);
471 : } else {
472 0 : cParameters(0) = pParameters(1);
473 0 : cParameters(1) = pParameters(0);
474 : }
475 : }
476 133 : Bool flipped = gaussShape.fromPixel (cParameters, dCoord);
477 266 : gauss2d = Angular2DGaussian(
478 266 : gaussShape.majorAxis(), gaussShape.minorAxis(),
479 266 : gaussShape.positionAngle()
480 133 : );
481 266 : return flipped;
482 : }
483 :
484 0 : Double SkyComponentFactory::_worldWidthToPixel (
485 : Double positionAngle,
486 : const Quantum<Double>& length,
487 : const CoordinateSystem& cSys,
488 : const IPosition& pixelAxes
489 : ) {
490 0 : Int worldAxis0 = cSys.pixelAxisToWorldAxis(pixelAxes(0));
491 0 : Int worldAxis1 = cSys.pixelAxisToWorldAxis(pixelAxes(1));
492 :
493 : // Units of the axes must be consistent for now.
494 : // I will be able to relax this criterion when I get the time
495 :
496 0 : Vector<String> units = cSys.worldAxisUnits();
497 0 : Unit unit0(units(worldAxis0));
498 0 : Unit unit1(units(worldAxis1));
499 0 : ThrowIf(
500 : unit0 != unit1,
501 : "Units of the two axes must be conformant"
502 : );
503 0 : Unit unit(unit0);
504 :
505 : // Check units are ok
506 :
507 0 : if (!length.check(unit.getValue())) {
508 0 : ostringstream oss;
509 0 : oss << "The units of the world length (" << length.getFullUnit().getName()
510 0 : << ") are not consistent with those of Coordinate System ("
511 0 : << unit.getName() << ")";
512 0 : ThrowCc(oss.str());
513 : }
514 :
515 0 : Double w0 = cos(positionAngle) * length.getValue(unit);
516 0 : Double w1 = sin(positionAngle) * length.getValue(unit);
517 :
518 : // Find pixel coordinate of tip of axis relative to reference pixel
519 :
520 0 : Vector<Double> world = cSys.referenceValue().copy();
521 0 : world(worldAxis0) += w0;
522 0 : world(worldAxis1) += w1;
523 :
524 0 : Vector<Double> pixel;
525 0 : ThrowIf(
526 : ! cSys.toPixel (pixel, world),
527 : cSys.errorMessage()
528 : );
529 :
530 0 : return hypot(pixel(pixelAxes(0)), pixel(pixelAxes(1)));
531 : }
532 :
533 0 : Quantum<Double> SkyComponentFactory::_pixelWidthToWorld (
534 : Double positionAngle, Double length,
535 : const CoordinateSystem& cSys2,
536 : const IPosition& pixelAxes
537 : ) {
538 0 : CoordinateSystem cSys(cSys2);
539 0 : Int worldAxis0 = cSys.pixelAxisToWorldAxis(pixelAxes(0));
540 0 : Int worldAxis1 = cSys.pixelAxisToWorldAxis(pixelAxes(1));
541 :
542 : // Units of the axes must be consistent for now.
543 : // I will be able to relax this criterion when I get the time
544 :
545 0 : Vector<String> units = cSys.worldAxisUnits().copy();
546 0 : Unit unit0(units(worldAxis0));
547 0 : Unit unit1(units(worldAxis1));
548 0 : ThrowIf(
549 : unit0 != unit1,
550 : "Units of the axes must be conformant"
551 : );
552 :
553 : // Set units to be the same for both axes
554 :
555 0 : units(worldAxis1) = units(worldAxis0);
556 0 : ThrowIf(
557 : !cSys.setWorldAxisUnits(units),
558 : cSys.errorMessage()
559 : );
560 :
561 0 : Double p0 = cos(positionAngle) * length;
562 0 : Double p1 = sin(positionAngle) * length;
563 :
564 : // Find world coordinate of tip of length relative to reference pixel
565 :
566 0 : Vector<Double> pixel= cSys.referencePixel().copy();
567 0 : pixel(pixelAxes(0)) += p0;
568 0 : pixel(pixelAxes(1)) += p1;
569 :
570 0 : Vector<Double> world;
571 0 : ThrowIf(
572 : ! cSys.toWorld(world, pixel),
573 : cSys.errorMessage()
574 : );
575 0 : Double lengthInWorld = hypot(world(worldAxis0), world(worldAxis1));
576 0 : return Quantum<Double>(lengthInWorld, Unit(units(worldAxis0)));
577 : }
578 :
579 :
580 :
581 : using namespace casacore;
582 : } // end namespace casa
|
