Line data    Source code

       1             : //# SepImageConvolver.cc:  separable convolution of an image
       2             : //# Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003
       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: aips2-request@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: SepImageConvolver.tcc 20615 2009-06-09 02:16:01Z Malte.Marquarding $
      27             : //   
      28             : 
      29             : #include <imageanalysis/ImageAnalysis/SepImageConvolver.h>
      30             : 
      31             : #include <casacore/casa/aips.h>
      32             : 
      33             : #include <casacore/casa/Arrays/Vector.h>
      34             : #include <casacore/casa/Arrays/ArrayMath.h>
      35             : #include <casacore/casa/Containers/Block.h>
      36             : #include <casacore/casa/Exceptions/Error.h>
      37             : #include <casacore/casa/Logging/LogIO.h>
      38             : #include <casacore/scimath/Mathematics/Convolver.h>
      39             : #include <casacore/casa/Quanta/UnitMap.h>
      40             : #include <casacore/casa/Quanta/Quantum.h>
      41             : #include <casacore/casa/BasicSL/String.h>
      42             : 
      43             : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
      44             : #include <casacore/images/Images/PagedImage.h>
      45             : #include <casacore/images/Regions/ImageRegion.h>
      46             : #include <casacore/images/Images/SubImage.h>
      47             : #include <casacore/images/Images/TempImage.h>
      48             : #include <casacore/lattices/Lattices/LatticeIterator.h>
      49             : #include <casacore/lattices/Lattices/LatticeUtilities.h>
      50             : #include <casacore/lattices/Lattices/TiledLineStepper.h>
      51             : #include <casacore/scimath/Mathematics/VectorKernel.h>
      52             : #include <casacore/casa/System/ProgressMeter.h>
      53             : 
      54             : #include <sstream>
      55             : #include <iomanip>
      56             : 
      57             : 
      58             : namespace casa { //# NAMESPACE CASA - BEGIN
      59             : 
      60             : template <class T> 
      61           0 : SepImageConvolver<T>::SepImageConvolver(
      62             :         const casacore::ImageInterface<T>& image, casacore::LogIO &os, casacore::Bool showProgress
      63             : )
      64             : : itsImagePtr(image.cloneII()),
      65             :   itsOs(os),
      66             :   itsAxes(0),
      67             :   itsVectorKernels(0),
      68           0 :   itsShowProgress(showProgress) {}
      69             : 
      70             : template <class T>
      71             : SepImageConvolver<T>::SepImageConvolver(const SepImageConvolver<T> &other)
      72             : : itsImagePtr(0) {
      73             :    operator=(other);
      74             : }
      75             : 
      76             : 
      77             : template <class T> 
      78           0 : SepImageConvolver<T>::~SepImageConvolver () {
      79           0 :    delete itsImagePtr;
      80           0 :    itsImagePtr = 0;
      81           0 :    const casacore::uInt n = itsVectorKernels.nelements();
      82           0 :    for (casacore::uInt i=0; i<n; i++) {
      83           0 :       delete itsVectorKernels[i];
      84           0 :       itsVectorKernels[i] = 0;
      85             :    }
      86           0 : }
      87             : 
      88             : template <class T>
      89             : SepImageConvolver<T> &SepImageConvolver<T>::operator=(
      90             :         const SepImageConvolver<T> &other
      91             : ) {
      92             :    if (this != &other) {
      93             :       if (itsImagePtr!=0) delete itsImagePtr;
      94             :       itsImagePtr = other.itsImagePtr->cloneII();
      95             : 
      96             :       itsOs = other.itsOs;
      97             :       itsAxes = other.itsAxes.copy();
      98             :       itsShowProgress = other.itsShowProgress;
      99             : 
     100             :       const casacore::uInt n = itsVectorKernels.nelements();
     101             :       for (casacore::uInt i=0; i<n; i++) {
     102             :          delete itsVectorKernels[i];
     103             :          itsVectorKernels[i] = 0;
     104             :       }
     105             : 
     106             :       const casacore::uInt n2 = other.itsVectorKernels.nelements();
     107             :       itsVectorKernels.resize(n2);
     108             :       for (casacore::uInt i=0; i<n2; i++) {      
     109             :          itsVectorKernels[i] = new casacore::Vector<T>((other.itsVectorKernels[i])->copy());
     110             :       }
     111             :    }
     112             :    return *this;
     113             : }
     114             : 
     115             : 
     116             : 
     117             : template <class T>
     118             : void SepImageConvolver<T>::setKernel(
     119             :         casacore::uInt axis, const casacore::Vector<T>& kernel
     120             : ) {
     121             :    _checkAxis(axis);
     122             : 
     123             :    casacore::uInt n = itsVectorKernels.nelements() + 1;
     124             :    itsVectorKernels.resize(n, true);
     125             :    itsVectorKernels[n-1] = new casacore::Vector<T>(kernel.copy());
     126             :    itsAxes.resize(n, true);
     127             :    itsAxes(n-1) = axis;
     128             : }
     129             : 
     130             : template <class T>
     131           0 : void SepImageConvolver<T>::setKernel(
     132             :         casacore::uInt axis, casacore::VectorKernel::KernelTypes kernelType,
     133             :     const casacore::Quantum<casacore::Double>& width, casacore::Bool autoScale,
     134             :     casacore::Bool useImageShapeExactly, casacore::Double scale
     135             : ) {
     136             :         // Catch pixel units
     137             : 
     138           0 :    casacore::UnitMap::putUser("pix",casacore::UnitVal(1.0), "pixel units");
     139           0 :    casacore::String sunit = width.getFullUnit().getName();
     140           0 :    if (sunit==casacore::String("pix")) {
     141           0 :       setKernel (axis, kernelType, width.getValue(), autoScale, useImageShapeExactly, scale);
     142           0 :       itsOs.output() << "Axis " << axis+1<< " : setting width "  << width << endl;
     143           0 :       itsOs << casacore::LogIO::NORMAL;
     144           0 :       return;
     145             :    }
     146           0 :    _checkAxis(axis);
     147             : 
     148             :    // Convert width to pixels
     149             : 
     150           0 :    casacore::CoordinateSystem cSys = itsImagePtr->coordinates();
     151           0 :    casacore::Int worldAxis = cSys.pixelAxisToWorldAxis(axis);
     152           0 :    casacore::Double inc = cSys.increment()(worldAxis);
     153             : 
     154           0 :    casacore::Unit unit = casacore::Unit(cSys.worldAxisUnits()(worldAxis));
     155           0 :    if (width.getFullUnit()!=unit) {
     156           0 :       itsOs << "Specified width units (" << width.getUnit() 
     157             :             << ") are inconsistent with image axis unit (" 
     158           0 :             << unit.getName() << casacore::LogIO::EXCEPTION;
     159             :    }
     160           0 :    casacore::Double width2 = abs(width.getValue(unit)/inc);
     161             : 
     162           0 :    itsOs.output() << "Axis " << axis+1<< " : setting width "
     163           0 :          << width << " = " << width2 << " pixels" << endl;
     164           0 :    itsOs <<  casacore::LogIO::NORMAL;
     165           0 :    setKernel(axis, kernelType, width2, autoScale, useImageShapeExactly, scale);
     166             : }
     167             : 
     168             : template <class T>
     169           0 : void SepImageConvolver<T>::setKernel(
     170             :         casacore::uInt axis, casacore::VectorKernel::KernelTypes kernelType,
     171             :     casacore::Double width, casacore::Bool autoScale,
     172             :     casacore::Bool useImageShapeExactly, casacore::Double scale
     173             : ) {
     174           0 :    _checkAxis(axis);
     175             : 
     176             :    // T can only be casacore::Float or Double
     177             : 
     178           0 :    casacore::Bool peakIsUnity = !autoScale;
     179           0 :    casacore::uInt shape = itsImagePtr->shape()(axis);
     180           0 :    casacore::Vector<T> x = casacore::VectorKernel::make(kernelType, T(width), shape, 
     181             :                                     useImageShapeExactly, peakIsUnity);
     182           0 :    if (!autoScale && !casacore::near(scale,1.0)) x *= casacore::Float(scale);
     183           0 :    casacore::uInt n = itsVectorKernels.nelements() + 1;
     184           0 :    itsVectorKernels.resize(n, true);
     185           0 :    itsVectorKernels[n-1] = new casacore::Vector<T>(x.copy());
     186             : 
     187           0 :    itsAxes.resize(n, true);
     188           0 :    itsAxes(n-1) = axis;
     189           0 : }
     190             : 
     191             : 
     192             : template <class T>
     193             : Vector<T> SepImageConvolver<T>::getKernel(casacore::uInt axis) {
     194             :    for (casacore::uInt i=0; i<itsAxes.nelements(); i++) {
     195             :       if (axis==itsAxes(i)) {
     196             :          return *(itsVectorKernels[i]);
     197             :       }
     198             :    }
     199             :    itsOs << "There is no kernel for the specified axis" << casacore::LogIO::EXCEPTION;
     200             :    return casacore::Vector<T>(0);
     201             : }
     202             : 
     203             : template <class T>
     204             : uInt SepImageConvolver<T>::getKernelShape(casacore::uInt axis) {
     205             :    for (casacore::uInt i=0; i<itsAxes.nelements(); i++) {
     206             :       if (axis==itsAxes(i)) {
     207             :          return itsVectorKernels[i]->nelements();
     208             :       }
     209             :    }
     210             :    itsOs << "There is no kernel for the specified axis" << casacore::LogIO::EXCEPTION;
     211             :    return 0;
     212             : }
     213             : 
     214             : template <class T>
     215           0 : void SepImageConvolver<T>::convolve(casacore::ImageInterface<T>& imageOut) {
     216           0 :    const casacore::uInt nAxes = itsAxes.nelements();
     217           0 :    if (nAxes==0) {
     218           0 :       itsOs << "You haven't specified any axes to convolve" << casacore::LogIO::EXCEPTION;
     219             :    }
     220             : 
     221             :    // Some checks
     222             : 
     223           0 :    casacore::IPosition shape = itsImagePtr->shape();
     224           0 :    if (!shape.isEqual(imageOut.shape())) {
     225           0 :       itsOs << "Image shapes are different" << casacore::LogIO::EXCEPTION;
     226             :    }
     227           0 :    casacore::CoordinateSystem cSys = itsImagePtr->coordinates();
     228           0 :    if (!cSys.near(imageOut.coordinates())) {
     229           0 :       itsOs << casacore::LogIO::WARN << "Image CoordinateSystems differ - this may be unwise"
     230           0 :             << casacore::LogIO::POST;
     231             :    }
     232             : 
     233             :    // Give the output image a mask if needed and make it the default
     234             : 
     235           0 :    if (itsImagePtr->isMasked() && !imageOut.isMasked()) {
     236           0 :       if (imageOut.canDefineRegion()) {
     237           0 :          casacore::String maskName = imageOut.makeUniqueRegionName (casacore::String("mask"), 0);
     238           0 :          imageOut.makeMask(maskName, true, true);
     239           0 :          itsOs << casacore::LogIO::NORMAL << "Created mask " << maskName 
     240           0 :                << " and make it the default" << casacore::LogIO::POST;
     241             :       } else {
     242           0 :          itsOs << casacore::LogIO::NORMAL << "Cannot create a mask for this output image" << casacore::LogIO::POST;
     243             :       }
     244             :    }
     245             : 
     246             :    // First copy input to output. We must replace masked pixels by zeros.  These reflect
     247             :    // both the pixel mask and the region mask.  We also set the output mask to the input mask
     248             :  
     249           0 :    casacore::LatticeUtilities::copyDataAndMask(itsOs, imageOut, *itsImagePtr, true);
     250             : 
     251             :    // casacore::Smooth in situ.
     252             :       
     253           0 :    casacore::IPosition niceShape = imageOut.niceCursorShape();
     254           0 :    casacore::uInt axis = 0;
     255           0 :    for (casacore::uInt i=0; i<nAxes; i++) {
     256           0 :       axis = itsAxes(i);
     257           0 :       itsOs << casacore::LogIO::NORMAL << "Convolving axis " << axis+1 << casacore::LogIO::POST;
     258           0 :       const casacore::Int n = shape(axis)/niceShape(axis);
     259           0 :       if (n*niceShape(axis)!=shape(axis)) {
     260           0 :          itsOs << casacore::LogIO::WARN 
     261             :                << "The tile shape is not integral along this axis, performance may degrade" 
     262           0 :                << casacore::LogIO::POST;
     263             :       }
     264           0 :       _smoothProfiles (imageOut, axis, *(itsVectorKernels[i]));
     265             :    }  
     266           0 : }
     267             :  
     268             : template <class T>
     269             : void SepImageConvolver<T>::_zero() {
     270             :    if (itsImagePtr->isMasked()) {
     271             :       itsOs << casacore::LogIO::NORMAL << "Zero masked pixels" << casacore::LogIO::POST;
     272             : 
     273             :       casacore::LatticeIterator<T> iter(*itsImagePtr);
     274             :       casacore::Bool deleteData, deleteMask;
     275             :       casacore::IPosition shape = iter.rwCursor().shape();
     276             :       casacore::Array<T> data(shape);
     277             :       casacore::Array<casacore::Bool> mask(shape);
     278             : 
     279             :       for (iter.reset(); !iter.atEnd(); iter++) {
     280             :          shape = iter.rwCursor().shape();
     281             :          if (!data.shape().isEqual(shape)) data.resize(shape);
     282             :          if (!mask.shape().isEqual(shape)) mask.resize(shape);
     283             : 
     284             :          itsImagePtr->getSlice(data, iter.position(), shape);
     285             :          itsImagePtr->getMaskSlice(mask, iter.position(), shape);
     286             : 
     287             :          T* pData = data.getStorage(deleteData);
     288             :          const casacore::Bool* pMask = mask.getStorage(deleteMask);
     289             : 
     290             :          for (casacore::Int i=0; i<shape.product(); i++) {
     291             :             if (!pMask[i]) pData[i] = 0.0;
     292             :          }
     293             : 
     294             :          data.putStorage(pData, deleteData);
     295             :          mask.freeStorage(pMask, deleteMask);
     296             :       }
     297             :    }     
     298             : }        
     299             :  
     300             : template <class T>
     301           0 : void SepImageConvolver<T>::_smoothProfiles (
     302             :         casacore::ImageInterface<T>& in, const casacore::Int& axis,
     303             :         const casacore::Vector<T>& psf
     304             : ) {
     305           0 :   casacore::ProgressMeter* pProgressMeter = 0;
     306           0 :   if (itsShowProgress) {
     307           0 :      casacore::Double nMin = 0.0;
     308           0 :      casacore::Double nMax = 1.0;
     309           0 :      for (casacore::Int i=0; i<casacore::Int(in.shape().nelements()); i++) {
     310           0 :         if (i!=axis) {
     311           0 :            nMax *= in.shape()(i);
     312             :         }
     313             :      }
     314           0 :      ostringstream oss;
     315           0 :      oss << "Convolve Image Axis " << axis+1;
     316           0 :      pProgressMeter = new casacore::ProgressMeter(nMin, nMax, casacore::String(oss),
     317             :                                         casacore::String("Spectrum Convolutions"),
     318             :                                         casacore::String(""), casacore::String(""),
     319           0 :                                         true, max(1,casacore::Int(nMax/20)));   
     320             :   }
     321             : 
     322           0 :   casacore::TiledLineStepper navIn(in.shape(),
     323             :                          in.niceCursorShape(),
     324             :                          axis);
     325           0 :   casacore::LatticeIterator<T> inIt(in, navIn);
     326           0 :   casacore::Vector<T> result(in.shape()(axis));
     327           0 :   casacore::IPosition shape(1, in.shape()(axis));
     328           0 :   casacore::Convolver<T> conv(psf, shape);
     329             : 
     330           0 :   casacore::uInt i = 0;
     331           0 :   while (!inIt.atEnd()) {
     332           0 :     conv.linearConv(result, inIt.vectorCursor());
     333           0 :     inIt.woVectorCursor() = result;
     334             : 
     335           0 :     if (itsShowProgress) pProgressMeter->update(casacore::Double(i));
     336           0 :     inIt++;
     337           0 :     i++;
     338             :   }
     339           0 :   if (itsShowProgress) delete pProgressMeter;
     340           0 : }
     341             : 
     342             : template <class T>
     343           0 : void SepImageConvolver<T>::_checkAxis(casacore::uInt axis) {
     344           0 :    if (axis>itsImagePtr->ndim()-1) {
     345           0 :       itsOs << "Given pixel axis " << axis 
     346           0 :             << " is greater than the number of axes in the image" << casacore::LogIO::EXCEPTION;
     347             :    }
     348           0 :    const casacore::uInt n = itsAxes.nelements();
     349           0 :    for (casacore::uInt i=0; i<n; i++) {
     350           0 :       if (axis==itsAxes(i)) {
     351           0 :          itsOs << "You have already given this axis to be convolved" << casacore::LogIO::EXCEPTION;
     352             :       }
     353             :    }
     354           0 : }
     355             : 
     356             : } //# NAMESPACE CASA - END
     357             :