Line data    Source code

       1             : //# ImageMetaData.cc: Meta information for Images
       2             : //# Copyright (C) 2009
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: ImageMetaData.cc 20886 2010-04-29 14:06:56Z gervandiepen $
      27             : 
      28             : #include <components/ComponentModels/GaussianDeconvolver.h>
      29             : #include <casacore/casa/BasicMath/Math.h>
      30             : #include <casacore/scimath/Mathematics/GaussianBeam.h>
      31             : 
      32             : using namespace casacore;
      33             : namespace casa {
      34             : 
      35          14 : Bool GaussianDeconvolver::deconvolve(
      36             :         Angular2DGaussian& deconvolvedSize,
      37             :         const Angular2DGaussian& convolvedSize,
      38             :         const GaussianBeam& beam
      39             : ) {
      40          14 :         Unit radians(String("rad"));
      41          14 :         Unit positionAngleModelUnit = deconvolvedSize.getPA(false).getFullUnit();
      42          14 :         Unit majorAxisModelUnit = deconvolvedSize.getMajor().getFullUnit();
      43          14 :         Unit minorAxisModelUnit = deconvolvedSize.getMinor().getFullUnit();
      44             : 
      45             :         // Get values in radians
      46          14 :         Double majorSource = convolvedSize.getMajor().getValue(radians);
      47          14 :         Double minorSource = convolvedSize.getMinor().getValue(radians);
      48          14 :         Double thetaSource = convolvedSize.getPA(true).getValue(radians);
      49          14 :         Double majorBeam = beam.getMajor().getValue(radians);
      50          14 :         Double minorBeam = beam.getMinor().getValue(radians);
      51          14 :         Double thetaBeam = beam.getPA(true).getValue(radians);
      52             :         // Do the sums
      53             : 
      54          14 :         Double alpha  = square(majorSource*cos(thetaSource)) +
      55          14 :                 square(minorSource*sin(thetaSource)) -
      56          14 :                 square(majorBeam*cos(thetaBeam)) -
      57          14 :                 square(minorBeam*sin(thetaBeam));
      58          14 :         Double beta   = square(majorSource*sin(thetaSource)) +
      59          14 :                 square(minorSource*cos(thetaSource)) -
      60          14 :                 square(majorBeam*sin(thetaBeam)) -
      61          14 :                 square(minorBeam*cos(thetaBeam));
      62          14 :         Double gamma  = 2 * ( (square(minorSource)-square(majorSource))*sin(thetaSource)*cos(thetaSource) -
      63          14 :                 (square(minorBeam)-square(majorBeam))*sin(thetaBeam)*cos(thetaBeam) );
      64             :         // Set result in radians
      65             : 
      66          14 :         Double s = alpha + beta;
      67          14 :         Double t = sqrt(square(alpha-beta) + square(gamma));
      68          14 :         Double limit = min(majorSource,minorSource);
      69          14 :         limit = min(limit,majorBeam);
      70          14 :         limit = min(limit,minorBeam);
      71          14 :         limit = 0.1*limit*limit;
      72             : 
      73          14 :         if(alpha<0.0 || beta<0.0 || s<t) {
      74           0 :                 if(0.5*(s-t)<limit && alpha>-limit && beta>-limit) {
      75             :                         // Point source. Fill in values of beam
      76           0 :                         deconvolvedSize = GaussianBeam(
      77           0 :                                 Quantity(beam.getMajor().get(majorAxisModelUnit)),
      78           0 :                                 Quantity(beam.getMinor().get(minorAxisModelUnit)),
      79           0 :                                 Quantity(beam.getPA(true).get(positionAngleModelUnit))
      80           0 :                         );
      81             :                         // unwrap
      82           0 :                         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
      83           0 :                         return true;
      84             :                 }
      85             :                 else {
      86           0 :                         throw AipsError("Source may be only (slightly) resolved in one direction");
      87             :                 }
      88             :         }
      89          14 :         Quantity majax(sqrt(0.5*(s+t)), radians);
      90          14 :         majax.convert(majorAxisModelUnit);
      91          14 :         Quantity minax(sqrt(0.5*(s-t)), radians);
      92          14 :         minax.convert(minorAxisModelUnit);
      93             :         Quantity pa(
      94          14 :                 abs(gamma)+abs(alpha-beta) == 0.0
      95          14 :                         ? 0.0
      96           0 :                         : 0.5*atan2(-gamma,alpha-beta),
      97          14 :                 radians);
      98          14 :         pa.convert(positionAngleModelUnit);
      99          14 :         deconvolvedSize = GaussianBeam(majax, minax, pa);
     100             :         // unwrap
     101          14 :         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
     102          14 :         return false;
     103          14 : }
     104             : 
     105             : 
     106             : 
     107             : } //# NAMESPACE CASA - END
     108             :