Line data Source code
1 : //# MultiTermMatrixCleaner.h: Minor Cycle for MSMFS deconvolution
2 : //# Copyright (C) 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: 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: MultiTermMatrixCleaner Urvashi R.V. 2010-12-04 <rurvashi@aoc.nrao.edu$
27 :
28 : #ifndef SYNTHESIS_MULTITERMLATTICECLEANER_H
29 : #define SYNTHESIS_MULTITERMLATTICECLEANER_H
30 :
31 : #include <casacore/scimath/Mathematics/FFTServer.h>
32 : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
33 :
34 : namespace casa { //# NAMESPACE CASA - BEGIN
35 :
36 : class MultiTermMatrixCleaner : public MatrixCleaner
37 : {
38 : public:
39 : // Create a cleaner
40 : MultiTermMatrixCleaner();
41 :
42 : // The copy constructor uses reference semantics
43 : // MultiTermMatrixCleaner(const MultiTermMatrixCleaner & other);
44 :
45 : // The assignment operator also uses reference semantics
46 : // MultiTermMatrixCleaner & operator=(const MultiTermMatrixCleaner & other);
47 :
48 : // The destructor resizes arrays to empty before destruction.
49 : ~MultiTermMatrixCleaner();
50 :
51 : // casacore::Input : number of Taylor terms
52 : // Reshapes PtrBlocks to hold the correct number of PSFs and Residual images
53 : casacore::Bool setntaylorterms(const int & nterms);
54 :
55 : // casacore::Input : scales
56 : casacore::Bool setscales(const casacore::Vector<casacore::Float> & scales);
57 :
58 : // Initialize all the memory being used.
59 : casacore::Bool initialise(casacore::Int nx,casacore::Int ny);
60 :
61 : // Calculate Hessian elements and check for invertibility
62 : // Does not have to be called externally, but can be. Either way, it executes only once.
63 : casacore::Int computeHessianPeak();
64 :
65 : // casacore::Input : psfs and dirty images
66 : casacore::Bool setpsf(int order, casacore::Matrix<casacore::Float> & psf);
67 :
68 : // casacore::Input : psfs and dirty images
69 : casacore::Bool setresidual(int order, casacore::Matrix<casacore::Float> & dirty);
70 :
71 : // casacore::Input : model images
72 : casacore::Bool setmodel(int order, casacore::Matrix<casacore::Float> & model);
73 :
74 : // casacore::Input : mask
75 : casacore::Bool setmask(casacore::Matrix<casacore::Float> & mask);
76 :
77 : // Run the minor cycle
78 : casacore::Int mtclean(casacore::Int maxniter, casacore::Float stopfraction, casacore::Float inputgain, casacore::Float userthreshold);
79 :
80 : // Output : Model images
81 : casacore::Bool getmodel(int order, casacore::Matrix<casacore::Float> & model);
82 :
83 : // Output : psfs and dirty images
84 : casacore::Bool getresidual(int order, casacore::Matrix<casacore::Float> & residual);
85 :
86 : // Compute principal solution - in-place on the residual images in vecDirty.
87 : casacore::Bool computeprincipalsolution();
88 :
89 : // Output : Hessian matrix
90 : casacore::Bool getinvhessian(casacore::Matrix<casacore::Double> & invhessian);
91 :
92 : // Output : Peak residual computed from matR_p (residual convolved with PSF).
93 137 : casacore::Float getpeakresidual(){return rmaxval_p;}
94 :
95 :
96 : private:
97 : casacore::LogIO os;
98 :
99 : using MatrixCleaner::itsCleanType;
100 : using MatrixCleaner::itsMaxNiter;
101 : using MatrixCleaner::itsGain;
102 : using MatrixCleaner::itsThreshold;
103 : using MatrixCleaner::itsMask;
104 : using MatrixCleaner::itsPositionPeakPsf;
105 : // using MatrixCleaner::itsScaleMasks;
106 : // using MatrixCleaner::itsScaleXfrs;
107 : // using MatrixCleaner::itsNscales;
108 : // using MatrixCleaner::itsScalesValid;
109 : using MatrixCleaner::itsMaskThreshold;
110 :
111 : using MatrixCleaner::findMaxAbs;
112 : using MatrixCleaner::findMaxAbsMask;
113 : using MatrixCleaner::makeScale;
114 : // using MatrixCleaner::addTo;
115 : using MatrixCleaner::makeBoxesSameSize;
116 : using MatrixCleaner::validatePsf;
117 : //using MatrixCleaner::makeScaleMasks;
118 :
119 : casacore::Int ntaylor_p; // Number of terms in the Taylor expansion to use.
120 : casacore::Int psfntaylor_p; // Number of terms in the Taylor expansion for PSF.
121 : casacore::Int nscales_p; // Number of scales to use for the multiscale part.
122 : casacore::Int nx_p;
123 : casacore::Int ny_p;
124 : casacore::Int totalIters_p; // Total number of minor-cycle iterations
125 : casacore::Float globalmaxval_p;
126 : casacore::Int maxscaleindex_p;
127 : casacore::IPosition globalmaxpos_p;
128 : casacore::Int itercount_p; // Number of minor cycle iterations
129 : casacore::Int maxniter_p;
130 : casacore::Float stopfraction_p;
131 : casacore::Float inputgain_p;
132 : casacore::Float userthreshold_p;
133 : casacore::Float prev_max_p;
134 : casacore::Float min_max_p;
135 : casacore::Float rmaxval_p;
136 :
137 : casacore::IPosition psfsupport_p;
138 : casacore::IPosition psfpeak_p;
139 : casacore::IPosition blc_p, trc_p, blcPsf_p, trcPsf_p;
140 :
141 : casacore::Vector<casacore::Float> scaleSizes_p; // casacore::Vector of scale sizes in pixels.
142 : casacore::Vector<casacore::Float> scaleBias_p; // casacore::Vector of scale biases !!
143 : casacore::Vector<casacore::Float> totalScaleFlux_p; // casacore::Vector of total scale fluxes.
144 : casacore::Vector<casacore::Float> totalTaylorFlux_p; // casacore::Vector of total flux in each taylor term.
145 : casacore::Vector<casacore::Float> maxScaleVal_p; // casacore::Vector for peaks at each scale size
146 : casacore::Vector<casacore::IPosition> maxScalePos_p; // casacore::Vector of peak positions at each scale size.
147 :
148 : casacore::IPosition gip;
149 : //casacore::Int nx,ny;
150 : casacore::Bool /*donePSF_p,*/ donePSP_p, doneCONV_p;
151 :
152 : casacore::Matrix<casacore::Complex> dirtyFT_p;
153 : casacore::Block<casacore::Matrix<casacore::Float> > vecScaleMasks_p;
154 :
155 : casacore::Matrix<casacore::Complex> cWork_p;
156 : casacore::Block<casacore::Matrix<casacore::Float> > vecWork_p;
157 :
158 : // h(s) [nx,ny,nscales]
159 : casacore::Block<casacore::Matrix<casacore::Float> > vecScales_p;
160 : casacore::Block<casacore::Matrix<casacore::Complex> > vecScalesFT_p;
161 :
162 : // B_k [nx,ny,ntaylor]
163 : // casacore::Block<casacore::Matrix<casacore::Float> > vecPsf_p;
164 : casacore::Block<casacore::Matrix<casacore::Complex> > vecPsfFT_p;
165 :
166 : // I_D : Residual/Dirty Images [nx,ny,ntaylor]
167 : casacore::Block<casacore::Matrix<casacore::Float> > vecDirty_p;
168 :
169 : // I_M : Model Images [nx,ny,ntaylor]
170 : casacore::Block<casacore::Matrix<casacore::Float> > vecModel_p;
171 :
172 : // Previous I_M : Model images at the start of mtclean(). Used for residual calculations.
173 : casacore::Block<casacore::Matrix<casacore::Float> > vecInitialModel_p;
174 :
175 : // casacore::Block <casacore::Matrix<casacore::Float> > vecScaleModel_p;
176 :
177 : // A_{smn} = B_{sm} * B{sn} [nx,ny,ntaylor,ntaylor,nscales,nscales]
178 : // A_{s1s2mn} = B_{s1m} * B{s2n} [nx,ny,ntaylor,ntaylor,nscales,nscales]
179 : casacore::Block<casacore::Matrix<casacore::Float> > cubeA_p;
180 :
181 : // R_{sk} = I_D * B_{sk} [nx,ny,ntaylor,nscales]
182 : casacore::Block<casacore::Matrix<casacore::Float> > matR_p;
183 :
184 : // a_{sk} = Solution vectors. [nx,ny,ntaylor,nscales]
185 : casacore::Block<casacore::Matrix<casacore::Float> > matCoeffs_p;
186 :
187 : // casacore::Memory to be allocated per Matrix
188 : casacore::Double memoryMB_p;
189 :
190 : // Solve [A][Coeffs] = [I_D * B]
191 : // Shape of A : [ntaylor,ntaylor]
192 : casacore::Block<casacore::Matrix<casacore::Double> > matA_p; // 2D matrix to be inverted.
193 : casacore::Block<casacore::Matrix<casacore::Double> > invMatA_p; // Inverse of matA_p;
194 :
195 : // FFTserver
196 : casacore::FFTServer<casacore::Float,casacore::Complex> fftcomplex;
197 :
198 : // Initial setup functions
199 : casacore::Int verifyScaleSizes();
200 : casacore::Int allocateMemory();
201 : casacore::Int setupScaleFunctions();
202 :
203 : // Setup per major cycle
204 : casacore::Int setupUserMask();
205 : casacore::Int computeFluxLimit(casacore::Float &fluxlimit, casacore::Float threshold);
206 : casacore::Int computeRHS();
207 :
208 : // Solver functions : minor-cycle iterations. Need to be efficient.
209 : casacore::Int solveMatrixEqn(casacore::Int ntaylor,casacore::Int scale, casacore::IPosition blc, casacore::IPosition trc);
210 : casacore::Int chooseComponent(casacore::Int ntaylor,casacore::Int scale, casacore::Int criterion, casacore::IPosition blc, casacore::IPosition trc);
211 : casacore::Int updateModelAndRHS(casacore::Float loopgain);
212 : casacore::Int updateRHS(casacore::Int ntaylor, casacore::Int scale, casacore::Float loopgain,casacore::Vector<casacore::Float> coeffs, casacore::IPosition blc, casacore::IPosition trc, casacore::IPosition blcPsf, casacore::IPosition trcPsf);
213 : casacore::Int checkConvergence(casacore::Int updatetype, casacore::Float &fluxlimit, casacore::Float &loopgain);
214 : casacore::Bool buildImagePatches();
215 :
216 : // Helper functions
217 : casacore::Int writeMatrixToDisk(casacore::String imagename, casacore::Matrix<casacore::Float> &themat);
218 : casacore::Int IND2(casacore::Int taylor,casacore::Int scale);
219 : casacore::Int IND4(casacore::Int taylor1, casacore::Int taylor2, casacore::Int scale1, casacore::Int scale2);
220 :
221 : casacore::Bool adbg;
222 : };
223 :
224 : } //# NAMESPACE CASA - END
225 :
226 : #endif
227 :
|
