Line data Source code
1 : //# LatConvEquation.cc: this defines LatConvEquation
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,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$
27 :
28 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
29 : #include <synthesis/MeasurementEquations/LinearModel.h>
30 : #include <casacore/lattices/LEL/LatticeExpr.h>
31 : #include <casacore/lattices/LEL/LatticeExprNode.h>
32 : #include <casacore/lattices/Lattices/SubLattice.h>
33 : #include <casacore/lattices/LRegions/LCBox.h>
34 : #include <casacore/casa/Utilities/Assert.h>
35 : #include <casacore/casa/Arrays/Vector.h>
36 : #include <casacore/casa/Arrays/Array.h>
37 : #include <casacore/casa/Arrays/ArrayMath.h>
38 : #include <casacore/casa/Arrays/IPosition.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 :
41 : using namespace casacore;
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 :
44 0 : LatConvEquation::LatConvEquation(Lattice<Float> & psf,
45 0 : Lattice<Float> & dirtyImage)
46 : :itsMeas(&dirtyImage),
47 : itsPsf(&psf),
48 0 : itsConv(psf, dirtyImage.shape(), true)
49 : {
50 0 : itsVirgin = true;
51 0 : itsRealPsfSize = psf.shape();
52 0 : itsPsfOrigin = itsPsf->shape()/2;
53 0 : }
54 :
55 0 : LatConvEquation::~LatConvEquation() {
56 0 : if (!itsVirgin) {
57 : // cout << "Deleting nonVirgin PSF" << endl;
58 0 : delete itsPsf;
59 : }
60 0 : }
61 :
62 0 : Bool LatConvEquation::evaluate(Lattice<Float> & result,
63 : const LinearModel<Lattice<Float> > & model){
64 0 : const Lattice<Float> & modelLattice = model.getModel();
65 : // DebugAssert(modelLattice.shape().isEqual(result.shape()), AipsError);
66 0 : itsConv.linear(result, modelLattice);
67 0 : return true;
68 : }
69 :
70 :
71 : Lattice<Float> *
72 0 : LatConvEquation::evaluate(const IPosition & position,
73 : const Float amplitude,
74 : const IPosition & modelSize){
75 0 : if (itsPsf->nelements() == 0){
76 0 : itsConv.getPsf(*itsPsf);
77 0 : itsPsfOrigin = itsPsf->shape()/2;
78 0 : itsRealPsfSize = itsPsf->shape();
79 : }
80 0 : IPosition psfSize = itsPsf->shape();
81 0 : IPosition blc = itsPsfOrigin - position;
82 0 : IPosition trc = blc + modelSize - 1;
83 : // Check if the required bounds are outside the psf. If they are then
84 : // resize the psf (with zero padding) to encompass the requested
85 : // convolution. Another way to do this is to just return the required
86 : // portion of the psf (of size modelSize) suitably padded. This will be
87 : // quicker and more memory efficient for just one evaluation, but if the
88 : // user is cleaning near the edges of their image, then it will have to be
89 : // done for each iteration that is near the edge. By resizing the whole
90 : // psf when necessary it will after a few resizes be at the required size
91 : // for all the iterations.
92 : //
93 : // If we resize the psf and itsVirgin == true, we need to create a new psf and
94 : // make sure itsVirgin = false
95 : //
96 0 : if ((min(blc.asVector()) < 0) ||
97 0 : (max((trc-psfSize).asVector()) >= 0))
98 : {
99 :
100 :
101 : // newtrc and newblc are the points in the newPsf that
102 : // refer to the old Psf (itsPsf).
103 0 : IPosition newSize(itsPsf->ndim(), 0);
104 0 : IPosition newtrc(itsPsf->ndim(), 0);
105 0 : IPosition newblc(itsPsf->ndim(), 0);
106 0 : for(uInt i = 0; i < itsPsf->ndim(); i++){
107 0 : newblc(i) = - min(blc(i), 0);
108 0 : newSize(i) = std::max(trc(i)+1, psfSize(i)) + newblc(i);
109 0 : newtrc(i) = newblc(i) + psfSize(i) - 1;
110 : }
111 : {
112 : // expand/pad itsPsf
113 : Lattice<Float> *newPsf;
114 0 : newPsf = new TempLattice<Float>(newSize);
115 0 : newPsf->set(0.0);
116 0 : LCBox box(newblc,newtrc,newSize);
117 0 : SubLattice<Float> newPsfSub (*newPsf, box, true);
118 0 : newPsfSub.copyData(*itsPsf); // fill the old Psf into the larger new one
119 :
120 0 : if (itsVirgin) {
121 : // just drop old itsPsf; someone else is responcible
122 0 : itsPsf = newPsf;
123 0 : itsVirgin = false;
124 : } else {
125 : // we must clean up old itsPsf
126 0 : delete itsPsf;
127 0 : itsPsf = newPsf;
128 : }
129 : }
130 0 : itsPsfOrigin = itsPsfOrigin + newblc;
131 0 : LCBox box(blc+newblc, trc+newtrc, newSize);
132 0 : Lattice<Float> *result = 0;
133 0 : if (!nearAbs(Double(amplitude),1.0)) {
134 0 : SubLattice<Float> newPsfSub( *itsPsf, box, true);
135 0 : result = new TempLattice<Float>(newPsfSub.shape());
136 0 : result->copyData((LatticeExpr<Float>)((amplitude) * newPsfSub ));
137 : } else {
138 0 : result = new SubLattice<Float>( *itsPsf, box, true);
139 : }
140 0 : return result;
141 : } else {
142 0 : LCBox box(blc,trc,itsPsf->shape());
143 0 : Lattice<Float> *result = 0;
144 0 : if (!nearAbs(Double(amplitude),1.0)) {
145 0 : SubLattice<Float> newPsfSub( *itsPsf, box, true);
146 0 : result = new TempLattice<Float>(newPsfSub.shape());
147 0 : result->copyData((LatticeExpr<Float>)((amplitude) * newPsfSub ));
148 : } else {
149 0 : result = new SubLattice<Float>( *itsPsf, box, true);
150 : }
151 0 : return result;
152 : }
153 : };
154 :
155 :
156 0 : Bool LatConvEquation::evaluate(Array<Float> & result,
157 : const IPosition & position,
158 : const Float amplitude,
159 : const IPosition & modelSize)
160 : {
161 0 : Lattice<Float> *resultLattice = 0;
162 0 : resultLattice = evaluate(position, amplitude, modelSize);
163 0 : if (resultLattice != 0) {
164 0 : result = resultLattice->get();
165 0 : delete resultLattice; resultLattice=0;
166 0 : return true;
167 : } else {
168 0 : return false;
169 : }
170 : };
171 :
172 :
173 0 : Bool LatConvEquation::residual(Lattice<Float> & result,
174 : const LinearModel< Lattice<Float> > & model) {
175 :
176 0 : if (evaluate(result, model)) {
177 0 : LatticeExpr<Float> expr = *itsMeas - result;
178 0 : result.copyData(expr);
179 0 : return true;
180 : }
181 0 : return false;
182 : }
183 :
184 :
185 0 : Bool LatConvEquation::residual(Lattice<Float> & result,
186 : Float & chisq,
187 : const LinearModel<Lattice<Float> > & model) {
188 0 : if (residual(result, model)) {
189 0 : LatticeExprNode myChisq = sum(result * result);
190 0 : chisq = myChisq.getFloat();
191 0 : return true;
192 : }
193 0 : return false;
194 : }
195 :
196 :
197 0 : Bool LatConvEquation::residual(Lattice<Float> & result,
198 : Float & chisq,
199 : Lattice<Float> & mask,
200 : const LinearModel<Lattice<Float> > & model) {
201 0 : if (residual(result, model)) {
202 0 : result.copyData( (LatticeExpr<Float>) (result * mask) );
203 0 : LatticeExprNode myChisq = sum(result * result);
204 0 : chisq = myChisq.getFloat();
205 0 : return true;
206 : }
207 0 : return false;
208 : }
209 :
210 0 : IPosition LatConvEquation::psfSize() {
211 0 : return itsConv.psfShape();
212 : }
213 :
214 : } //# NAMESPACE CASA - END
215 :
|
