Line data    Source code

       1             : //# Matrix Non Amnesiac Cleaner Copyright (C) 2015
       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  General Public License as published by
       6             : //# the Free Software Foundation; either version 3 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             : //# $Id:  $
      26             : 
      27             : #include <casacore/casa/Arrays/Matrix.h>
      28             : #include <casacore/casa/Arrays/Cube.h>
      29             : #include <casacore/casa/Arrays/MaskedArray.h>
      30             : #include <casacore/casa/Arrays/ArrayMath.h>
      31             : #include <casacore/casa/Arrays/MaskArrMath.h>
      32             : #include <casacore/casa/Arrays/MatrixMath.h>
      33             : #include <casacore/casa/IO/ArrayIO.h>
      34             : #include <casacore/casa/BasicMath/Math.h>
      35             : #include <casacore/casa/BasicSL/Complex.h>
      36             : #include <casacore/casa/Logging/LogIO.h>
      37             : #include <casacore/casa/OS/File.h>
      38             : #include <casacore/casa/Containers/Record.h>
      39             : 
      40             : #include <casacore/lattices/LRegions/LCBox.h>
      41             : #include <casacore/casa/Arrays/Slicer.h>
      42             : #include <casacore/scimath/Mathematics/FFTServer.h>
      43             : #include <casacore/casa/OS/HostInfo.h>
      44             : #include <casacore/casa/Arrays/ArrayError.h>
      45             : #include <casacore/casa/Arrays/ArrayIter.h>
      46             : #include <casacore/casa/Arrays/VectorIter.h>
      47             : 
      48             : 
      49             : #include <casacore/casa/Utilities/GenSort.h>
      50             : #include <casacore/casa/BasicSL/String.h>
      51             : #include <casacore/casa/Utilities/Assert.h>
      52             : #include <casacore/casa/Utilities/Fallible.h>
      53             : 
      54             : #include <casacore/casa/BasicSL/Constants.h>
      55             : 
      56             : #include <casacore/casa/Logging/LogSink.h>
      57             : #include <casacore/casa/Logging/LogMessage.h>
      58             : 
      59             : #include <synthesis/MeasurementEquations/MatrixNACleaner.h>
      60             : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
      61             : #include <casacore/coordinates/Coordinates/TabularCoordinate.h>
      62             : #ifdef _OPENMP
      63             : #include <omp.h>
      64             : #endif
      65             : using namespace std;
      66             : using namespace casacore;
      67             : namespace casa { //# NAMESPACE CASA - BEGIN
      68             : 
      69             :  
      70             : 
      71             :   
      72             :  
      73           0 : MatrixNACleaner::MatrixNACleaner():
      74           0 :   itsMask( ),
      75           0 :   itsDirty( ),
      76           0 :   itsBitPix( ),
      77           0 :   itsMaximumResidual(1e100),
      78           0 :   itsTotalFlux(0.0),
      79           0 :   itsSupport(3),
      80           0 :   psfShape_p(0),
      81           0 :   itsPositionPeakPsf(0),
      82           0 :   itsRms(0.0),
      83           0 :   typeOfMemory_p(2),
      84           0 :   numSigma_p(5.0)
      85             : {
      86             :   
      87             :  
      88           0 : }
      89             : 
      90             :  
      91             : 
      92           0 : MatrixNACleaner::MatrixNACleaner(const Matrix<Float> & psf,
      93           0 :                                  const Matrix<Float> &dirty, const Int memType, const Float  numSigma):
      94           0 :   itsBitPix( ), itsSupport(3), itsRms(0.0),typeOfMemory_p(memType), numSigma_p(numSigma)
      95             : {
      96           0 :   psfShape_p.resize(0, false);
      97           0 :   psfShape_p=psf.shape();
      98             :   // Check that everything is the same dimension and that none of the
      99             :   // dimensions is zero length.
     100           0 :   AlwaysAssert(psf.shape().nelements() == dirty.shape().nelements(),
     101             :                AipsError);
     102           0 :   AlwaysAssert(dirty.shape().product() != 0, AipsError);
     103             :  
     104           0 :   itsDirty=make_shared<Matrix<Float> >();
     105           0 :   itsDirty->reference(dirty);
     106           0 :   setPsf(psf);
     107           0 :   itsResidual=std::make_shared<Matrix<Float> >();
     108           0 :   itsResidual->assign(dirty);
     109           0 :   itsMask= std::make_shared<Matrix<Float> >(itsDirty->shape(), Float(0.0));
     110             :   //itsRms=rms(*itsResidual);
     111           0 :   itsMaximumResidual=1e100;
     112           0 : }
     113             : 
     114             : 
     115           0 :   Bool MatrixNACleaner::validatePsf(const Matrix<Float> & psf)
     116             : {
     117           0 :   LogIO os(LogOrigin("MatrixCleaner", "validatePsf()", WHERE));
     118             :   
     119             :   // Find the peak of the raw Psf
     120           0 :   AlwaysAssert(psf.shape().product() != 0, AipsError);
     121           0 :   Float maxPsf=0;
     122           0 :   itsPositionPeakPsf=IPosition(psf.shape().nelements(), 0);
     123           0 :   Int psfSupport = max(psf.shape().asVector())/2;
     124           0 :   MatrixCleaner::findPSFMaxAbs(psf, maxPsf, itsPositionPeakPsf, psfSupport);
     125           0 :   os << "Peak of PSF = " << maxPsf << " at " << itsPositionPeakPsf
     126           0 :      << LogIO::POST;
     127           0 :   return true;
     128           0 : }
     129             :   
     130             : 
     131           0 : void MatrixNACleaner::setPsf(const Matrix<Float>& psf){
     132           0 :   AlwaysAssert(validatePsf(psf), AipsError);
     133           0 :   psfShape_p.resize(0, false);
     134           0 :   psfShape_p=psf.shape();
     135           0 :   itsPsf=make_shared<Matrix<Float> >();
     136           0 :   itsPsf->reference(psf);
     137             :   //cout << "shapes " << itsXfr->shape() << " psf " << psf.shape() << endl;
     138           0 : }
     139             : 
     140           0 : MatrixNACleaner::MatrixNACleaner(const MatrixNACleaner & other)
     141             :    
     142             : {
     143           0 :   operator=(other);
     144           0 : }
     145             : 
     146           0 : MatrixNACleaner & MatrixNACleaner::operator=(const MatrixNACleaner & other) {
     147           0 :   if (this != &other) {
     148             :     
     149             :    
     150           0 :     itsMask = other.itsMask;
     151           0 :     itsDirty = other.itsDirty;
     152           0 :     itsResidual=other.itsResidual;
     153           0 :     itsBitPix=other.itsBitPix;
     154           0 :     itsTotalFlux=other.itsTotalFlux;
     155           0 :     psfShape_p.resize(0, false);
     156           0 :     psfShape_p=other.psfShape_p;
     157           0 :     itsSupport=other.itsSupport;
     158           0 :     itsRms=other.itsRms;
     159           0 :     typeOfMemory_p=other.typeOfMemory_p;
     160             :     
     161             :   }
     162           0 :   return *this;
     163             : }
     164             : 
     165           0 : MatrixNACleaner::~MatrixNACleaner()
     166             : {
     167             :   
     168             :   
     169           0 : }
     170             : 
     171             : 
     172             : 
     173             : 
     174             : // add a mask image
     175           0 :   void MatrixNACleaner::setMask(Matrix<Float> & mask) 
     176             : {
     177             :  
     178             : 
     179             :   //cerr << "Mask Shape " << mask.shape() << endl;
     180             :   // This is not needed after the first steps
     181           0 :   itsMask = std::make_shared<Matrix<Float> >(mask.shape());
     182           0 :   itsMask->reference(mask);
     183             :  
     184             : 
     185           0 : }
     186             :  
     187           0 : void MatrixNACleaner::setPixFlag(const Matrix<Bool> & bitpix) 
     188             : {
     189             :  
     190             : 
     191             :   //cerr << "Mask Shape " << mask.shape() << endl;
     192             :   // This is not needed after the first steps
     193           0 :   itsBitPix = std::make_shared<Matrix<Bool> >(bitpix.shape());
     194           0 :   itsBitPix->reference(bitpix);
     195             :  
     196             : 
     197           0 : }
     198             : 
     199             : 
     200             : // Set up the control parameters
     201           0 : void MatrixNACleaner::setcontrol(
     202             :                                    const Int niter,
     203             :                                    const Float gain,
     204             :                                    const Quantity& aThreshold,
     205             :                                    const Int masksupp, const Int memType, const Float numSigma)
     206             : {
     207             :   
     208           0 :   itsMaxNiter=niter;
     209           0 :   itsGain=gain;
     210           0 :   itsThreshold=aThreshold;
     211           0 :   itsSupport=masksupp;
     212           0 :   typeOfMemory_p=memType;
     213           0 :   numSigma_p=numSigma;
     214             :   
     215           0 : }
     216             : 
     217             : 
     218             : 
     219             : 
     220             : // Do the clean as set up
     221           0 : Int MatrixNACleaner::clean(Matrix<Float>& model)
     222             : {
     223           0 :   AlwaysAssert(model.shape()==itsDirty->shape(), AipsError);
     224             : 
     225           0 :   LogIO os(LogOrigin("MatrixCleaner", "clean()", WHERE));
     226             : 
     227             : 
     228           0 :   if(typeOfMemory_p==0)
     229           0 :     f_p=std::bind(&MatrixNACleaner::amnesiac, this, std::placeholders::_1);
     230           0 :   else if(typeOfMemory_p==1)
     231           0 :     f_p=std::bind(&MatrixNACleaner::weak, this, std::placeholders::_1);
     232           0 :   else if(typeOfMemory_p==3)
     233           0 :     f_p=std::bind(&MatrixNACleaner::strong, this, std::placeholders::_1);
     234             :   else
     235           0 :     f_p=std::bind(&MatrixNACleaner::medium, this, std::placeholders::_1);
     236           0 :   Int converged=0;
     237           0 :   itsTotalFlux=0.0;
     238           0 :    if(!itsBitPix){
     239           0 :      itsRms=rms(*itsResidual);
     240             :    }else{
     241           0 :      itsRms=rms(MaskedArray<Float>(*itsResidual, *itsBitPix));
     242             :    }
     243             : 
     244             :   // Define a subregion for the inner quarter
     245           0 :   IPosition blcDirty(model.shape().nelements(), 0);
     246           0 :   IPosition trcDirty(model.shape()-1);
     247             : 
     248           0 :   if(!itsMask){
     249           0 :     itsMask=make_shared<Matrix<Float> >(model.shape(), Float(0.0));
     250             :   } 
     251             : 
     252             :     
     253             :    
     254             :  
     255             :  
     256           0 :   os << "Starting iteration"<< LogIO::POST;
     257           0 :   Float maxima=0.0;
     258           0 :   IPosition posMaximum;
     259           0 :   itsIteration = 0;
     260           0 :   for (Int ii=0; ii < itsMaxNiter; ii++) {
     261           0 :     itsIteration++;
     262             :  
     263             :         
     264             :         
     265             :         
     266           0 :     if(!findMaxAbsMask(*itsResidual, *itsMask,
     267             :                        maxima, posMaximum, itsSupport))
     268           0 :       break;
     269             : 
     270             :         
     271             : 
     272             :    
     273             :    
     274             : 
     275             :     // Now add to the total flux
     276           0 :     itsTotalFlux += (maxima*itsGain);
     277             :     
     278             : 
     279             : 
     280           0 :     if(ii==0 ) {
     281           0 :       itsMaximumResidual=abs(maxima);
     282           0 :       os << "Initial maximum residual is " << itsMaximumResidual;
     283             :      
     284           0 :       os << LogIO::POST;
     285             :     }
     286             : 
     287             :     // Various ways of stopping:
     288             :     //    1. stop if below threshold
     289           0 :     if(abs(maxima)<threshold() ) {
     290           0 :       os << "Reached stopping threshold " << threshold() << " at iteration "<<
     291           0 :             ii << LogIO::POST;
     292           0 :       os << "Optimum flux is " << abs(maxima) << LogIO::POST;
     293           0 :       converged = 1;
     294           0 :       break;
     295             :     }
     296             :     
     297             :     /*
     298             :     if(progress) {
     299             :       progress->info(false, itsIteration, itsMaxNiter, maxima,
     300             :                      posMaximum, itsStrengthOptimum,
     301             :                      optimumScale, positionOptimum,
     302             :                      totalFlux, totalFluxScale,
     303             :                      itsJustStarting );
     304             :       itsJustStarting = false;
     305             :       } else*/ {
     306           0 :       if (itsIteration ==   1) {
     307           0 :           os << "iteration    MaximumResidual   CleanedFlux" << LogIO::POST;
     308             :       }
     309           0 :       if ((itsIteration % (itsMaxNiter/10 > 0 ? itsMaxNiter/10 : 1)) == 0) {
     310             :         //Good time to redo rms if necessary
     311           0 :         if(abs(maxima) >  3*itsRms){
     312           0 :           if(!itsBitPix){
     313           0 :             itsRms=rms(*itsResidual);
     314             :           }else{
     315           0 :             itsRms=rms(MaskedArray<Float>(*itsResidual, *itsBitPix));
     316             :           }
     317             :         }
     318             :         os << itsIteration <<"      "<< maxima<<"      "
     319           0 :            << itsTotalFlux << " rms " << itsRms << LogIO::POST;
     320             :       }
     321             :     }
     322             : 
     323             :    
     324             : 
     325             :     // Continuing: subtract the peak that we found from all dirty images
     326             :     // Define a subregion so that that the peak is centered
     327           0 :     IPosition support(model.shape());
     328             :     
     329             : 
     330           0 :     IPosition inc(model.shape().nelements(), 1);
     331             :     //cout << "support " << support.asVector()  << endl;
     332             :     //support(0)=1024;
     333             :     //support(1)=1024;
     334             :     //support(0)=min(Int(support(0)), Int(trcDirty(0)-blcDirty(0)));
     335             :     //support(1)=min(Int(support(1)), Int(trcDirty(1)-blcDirty(1)));
     336             :     // support(0)=min(Int(support(0)), (trcDirty(0)-blcDirty(0)+
     337             :     //                          Int(2*abs(positionOptimum(0)-blcDirty(0)/2.0-trcDirty(0)/2.0))));
     338             :     //support(1)=min(Int(support(1)), (trcDirty(1)-blcDirty(1)+
     339             :     //                          Int(2*abs(positionOptimum(1)-blcDirty(1)/2.0-trcDirty(1)/2.0))));
     340             : 
     341           0 :     IPosition blc(posMaximum-support/2);
     342           0 :     IPosition trc(posMaximum+support/2-1);
     343           0 :     LCBox::verify(blc, trc, inc, model.shape());
     344             :     
     345             :     //cout << "blc " << blc.asVector() << " trc " << trc.asVector() << endl;
     346             : 
     347           0 :     IPosition blcPsf(blc+itsPositionPeakPsf-posMaximum);
     348           0 :     IPosition trcPsf(trc+itsPositionPeakPsf-posMaximum);
     349           0 :     LCBox::verify(blcPsf, trcPsf, inc, model.shape());
     350           0 :     makeBoxesSameSize(blc,trc,blcPsf,trcPsf);
     351             :     // cout << "blcPsf " << blcPsf.asVector() << " trcPsf " << trcPsf.asVector() << endl;
     352             :     //cout << "blc " << blc.asVector() << " trc " << trc.asVector() << endl;
     353             :     //    LCBox subRegion(blc, trc, model.shape());
     354             :     //  LCBox subRegionPsf(blcPsf, trcPsf, model.shape());
     355             :     
     356             :    
     357           0 :     Matrix<Float> resSub=(*itsResidual)(blc,trc);
     358             :     
     359             :  
     360             :     // Now do the addition to the model image....
     361           0 :     model(posMaximum) += itsGain*maxima;
     362             : 
     363             :                         
     364             :   
     365             :       
     366             :         
     367           0 :     Matrix<Float> psfSub=(*itsPsf)(blcPsf, trcPsf);
     368           0 :     resSub -= itsGain*maxima*psfSub;
     369             :             
     370             :       
     371             :     
     372           0 :   }
     373             :   // End of iteration
     374             : 
     375             :  
     376             :     os << LogIO::NORMAL
     377           0 :        << "  " << "  Total Flux  " << itsTotalFlux
     378           0 :        << LogIO::POST;
     379             :   
     380             :   // Finish off the plot, etc.
     381             :   /*
     382             :   if(progress) {
     383             :     progress->info(true, itsIteration, itsMaxNiter, maxima, posMaximum,
     384             :                    itsStrengthOptimum,
     385             :                    optimumScale, positionOptimum,
     386             :                    totalFlux, totalFluxScale);
     387             :   }
     388             :   */
     389             : 
     390           0 :   if(!converged) {
     391           0 :     os << "Failed to reach stopping threshold" << LogIO::POST;
     392             :   }
     393             : 
     394           0 :   return converged;
     395           0 : }
     396             : 
     397             : 
     398             : 
     399             : 
     400             : 
     401             : 
     402           0 : Bool MatrixNACleaner::findMaxAbsMask(const Matrix<Float>& lattice,
     403             :                                    Matrix<Float>& mask,
     404             :                                               Float& peakval,
     405             :                                    IPosition& posPeak, const Int support)
     406             : {
     407             : 
     408             :   /*
     409             :     Here is the secret of non masking algorithm
     410             :     find peak of mask*resid -- m1*p1   p1 (value of resid at that peak)
     411             :     find peak of resid -- p2
     412             :     if p2 < m1*p1
     413             :        return p1 
     414             :        modify box (using user given support) at pos of p1 in mask to be f(p1) or m@p1 which ever is higher
     415             :    else 
     416             :         return p2
     417             :          modify box at pos of p2 in mask to be f(p2) or m@p2  which ever is higher
     418             :     
     419             :    stop modifying mask when peak reaches 5 (or user settable) sigma
     420             :    
     421             :    f(p) is a memory function...f(p)=1 implies no memory or normal clean
     422             :    default right now we are using f(p)=p
     423             :     a weak memory can be f(p)=1+ k*p  ( k << 1)  or p ^0.1
     424             :     a strong memory can be f(p)=p^2
     425             :    */
     426             : 
     427             :   //cerr << "SHAPES " << lattice.shape() << "  " << mask.shape() << endl;
     428           0 :   Matrix<Float> wgtresid=lattice*mask;
     429             :   
     430           0 :   IPosition posMaxWgt = IPosition(lattice.shape().nelements(), 0);
     431           0 :   Float maxValWgt=0.0;
     432             :   Float minValWgt;
     433           0 :   IPosition posMinWgt(lattice.shape().nelements(), 0);
     434           0 :   if(!itsBitPix)
     435           0 :     minMax(minValWgt, maxValWgt, posMinWgt, posMaxWgt, wgtresid);
     436             :   else
     437           0 :      minMax(minValWgt, maxValWgt, posMinWgt, posMaxWgt, MaskedArray<Float>(wgtresid, *itsBitPix));
     438           0 :   if(abs(minValWgt) > abs(maxValWgt)){
     439           0 :     posMaxWgt=posMinWgt;
     440           0 :     maxValWgt=minValWgt;
     441             :   }
     442           0 :   IPosition posMaxRes = IPosition(lattice.shape().nelements(), 0);
     443           0 :   Float maxValRes=0.0;
     444             :   Float minValRes;
     445           0 :   IPosition posMinRes(lattice.shape().nelements(), 0);
     446           0 :   if(!itsBitPix)
     447           0 :     minMax(minValRes, maxValRes, posMinRes, posMaxRes, lattice);
     448             :   else
     449           0 :     minMax(minValRes, maxValRes, posMinRes, posMaxRes,  MaskedArray<Float>(lattice, *itsBitPix));
     450           0 :   if(abs(minValRes) > abs(maxValRes)){
     451           0 :     posMaxRes=posMinRes;
     452           0 :     maxValRes=minValRes;
     453             :   }
     454           0 :   if(abs(maxValWgt) > abs(maxValRes)){
     455           0 :     posMaxRes=posMaxWgt;
     456           0 :     maxValRes=lattice(posMaxRes);
     457             :   }
     458           0 :   posPeak=posMaxRes;
     459           0 :   peakval=maxValRes;
     460           0 :   Int nx=lattice.shape()(0);
     461           0 :   Int ny=lattice.shape()(1);
     462           0 :   Float apeakval=abs(peakval);
     463           0 :   if(apeakval > numSigma_p*itsRms){
     464           0 :     for (Int y=-support; y < support+1; ++y){
     465           0 :       if(((posPeak[1]+y) >=0) && ((posPeak[1]+y) < ny)){
     466           0 :         for(Int x=-support; x < support+1; ++x){
     467           0 :           if(((posPeak[0]+x) >=0) && ((posPeak[0]+x) < nx))
     468           0 :             if(mask(x+posPeak[0], y+posPeak[1]) < f_p(apeakval))
     469           0 :               mask(x+posPeak[0], y+posPeak[1])=f_p(apeakval);
     470             :         }
     471             :       }
     472             :     }
     473             :   }
     474           0 :   itsMaximumResidual=peakval;
     475             :   //if(peakval > 2*itsRms)
     476             :   //  return true;
     477             : 
     478           0 :   return true;
     479           0 : }
     480           0 :   Float  MatrixNACleaner::amnesiac(const Float& ){
     481           0 :     return 1.0;
     482             :   }
     483             : 
     484           0 :   Float MatrixNACleaner::weak(const Float& v){
     485           0 :     return 1.0+0.1*v;
     486             :   }
     487             : 
     488           0 :   Float MatrixNACleaner::medium(const Float& v){
     489           0 :     return v;
     490             : 
     491             :   }
     492           0 :   Float MatrixNACleaner::strong(const Float& v){
     493           0 :     return v*v;
     494             :   }
     495             : 
     496             : 
     497           0 : void MatrixNACleaner::setDirty(const Matrix<Float>& dirty){
     498             : 
     499           0 :   itsDirty=make_shared<Matrix<Float> >(dirty.shape());
     500           0 :   itsDirty->reference(dirty);
     501           0 :    itsResidual=make_shared<Matrix<Float> >(dirty.shape());
     502           0 :   itsResidual->assign(dirty);
     503           0 :   itsRms=rms(*itsResidual);
     504             :   
     505             :   
     506             :   
     507             : 
     508           0 : } 
     509             : 
     510             : //Define the scales without doing anything else
     511             : // user will call make makePsfScales and makeDirtyScales like an adult in the know
     512             : 
     513             : 
     514             : 
     515             : 
     516           0 : void MatrixNACleaner::unsetMask()
     517             : {
     518             :  
     519           0 :   if(!itsMask)
     520           0 :     itsMask=nullptr;
     521             :   
     522           0 : }
     523             : 
     524             : 
     525             : 
     526             : 
     527             : 
     528             : 
     529             : 
     530             : 
     531           0 : Float MatrixNACleaner::threshold() const
     532             : {
     533           0 :   return Float(itsThreshold.getValue("Jy"));
     534             : }
     535             : 
     536             : 
     537             : 
     538           0 : void MatrixNACleaner::makeBoxesSameSize(IPosition& blc1, IPosition& trc1, 
     539             :                   IPosition &blc2, IPosition& trc2)
     540             : {
     541           0 :   const IPosition shape1 = trc1 - blc1;
     542           0 :   const IPosition shape2 = trc2 - blc2;
     543             : 
     544           0 :   AlwaysAssert(shape1.nelements() == shape2.nelements(), AipsError);
     545             :   
     546           0 :   if (shape1 == shape2) {
     547           0 :       return;
     548             :   }
     549           0 :   for (uInt i=0;i<shape1.nelements();++i) {
     550           0 :        Int minLength = shape1[i];
     551           0 :        if (shape2[i]<minLength) {
     552           0 :            minLength = shape2[i];
     553             :        }
     554           0 :        AlwaysAssert(minLength>=0, AipsError);
     555             :        //if (minLength % 2 != 0) {
     556             :            // if the number of pixels is odd, ensure that the centre stays 
     557             :            // the same by making this number even
     558             :            //--minLength; // this code is a mistake and should be removed
     559             :        //}
     560           0 :        const Int increment1 = shape1[i] - minLength;
     561           0 :        const Int increment2 = shape2[i] - minLength;
     562           0 :        blc1[i] += increment1/2;
     563           0 :        trc1[i] -= increment1/2 + (increment1 % 2 != 0 ? 1 : 0);
     564           0 :        blc2[i] += increment2/2;
     565           0 :        trc2[i] -= increment2/2 + (increment2 % 2 != 0 ? 1 : 0);
     566             :   }
     567           0 : }
     568             : 
     569             : 
     570             : 
     571             : } //# NAMESPACE CASA - END
     572             :