Line data Source code
1 : //# ConvolutionEquation.cc: this defines ConvolutionEquation
2 : //# Copyright (C) 1996,1997,1998,1999,2000,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$
27 :
28 : #include <synthesis/MeasurementEquations/ConvolutionEquation.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/BasicMath/Math.h>
31 : #include <casacore/casa/Arrays/Vector.h>
32 :
33 : using namespace casacore;
34 : namespace casa { //# NAMESPACE CASA - BEGIN
35 :
36 0 : ConvolutionEquation::ConvolutionEquation(){
37 0 : };
38 :
39 0 : ConvolutionEquation::ConvolutionEquation(const Array<Float> & psf,
40 0 : const Array<Float> & dirtyImage){
41 0 : theConv.setPsf(psf, dirtyImage.shape());
42 0 : theConv.setFastConvolve();
43 0 : theRealPsfSize = psf.shape();
44 0 : theMeas.resize(dirtyImage.shape());
45 0 : theMeas = dirtyImage;
46 0 : };
47 :
48 0 : ConvolutionEquation::ConvolutionEquation(const Array<Float>& psf,
49 0 : const MaskedArray<Float>& dirtyImage){
50 0 : theConv.setPsf(psf, dirtyImage.shape());
51 0 : theConv.setFastConvolve();
52 0 : theRealPsfSize = psf.shape();
53 0 : theMeas.resize(dirtyImage.shape());
54 0 : theMeas = dirtyImage.getArray();
55 0 : };
56 :
57 0 : ConvolutionEquation::~ConvolutionEquation(){
58 0 : };
59 :
60 0 : Bool ConvolutionEquation::evaluate(Array<Float> & result,
61 : const LinearModel<Array<Float> > & model){
62 0 : const Array<Float> modelArray = model.getModel();
63 0 : theConv.linearConv(result, modelArray, false);
64 0 : return true;
65 0 : };
66 :
67 0 : Bool ConvolutionEquation::evaluate(Array<Float> & result,
68 : const LinearModel< MaskedArray<Float> > & model) {
69 0 : const MaskedArray<Float> modelArray = model.getModel();
70 0 : theConv.linearConv(result, modelArray.getArray(), false);
71 0 : return true;
72 0 : }
73 :
74 0 : Bool ConvolutionEquation::evaluate(Array<Float> & result,
75 : const IPosition & position,
76 : const Float amplitude,
77 : const IPosition & modelSize){
78 0 : if (thePsf.nelements() == 0){
79 0 : thePsf = theConv.getPsf(false);
80 0 : thePsfOrigin = thePsf.shape()/2;
81 : }
82 0 : IPosition psfSize = thePsf.shape();
83 0 : IPosition blc = thePsfOrigin - position;
84 0 : IPosition trc = blc + modelSize - 1;
85 : // Check if the required bounds are outside the psf. If they are then
86 : // resize the psf (with zero padding) to encompass the requested
87 : // convolution. Another way to do this is to just return the required
88 : // portion of the psf (of size modelSize) suitably padded. This will be
89 : // quicker and more memory efficient for just one evaluation, but if the
90 : // user is cleaning near the edges of their image, then it will have to be
91 : // done for each iteration that is near the edge. By resizing the whole
92 : // psf when necessary it will after a few resizes be at the required size
93 : // for all the iterations.
94 0 : if ((min(blc.asVector()) < 0) ||
95 0 : (max((trc-psfSize).asVector()) >= 0))
96 : {
97 0 : IPosition newSize(thePsf.ndim(), 0);
98 0 : IPosition newtrc(thePsf.ndim(), 0);
99 0 : IPosition newblc(thePsf.ndim(), 0);
100 0 : for(uInt i = 0; i < thePsf.ndim(); i++){
101 0 : newblc(i) = -std::min(blc(i), (ssize_t)0);
102 0 : newSize(i) = std::max(trc(i)+1, psfSize(i)) + newblc(i);
103 0 : newtrc(i) = newblc(i) + psfSize(i) - 1;
104 : }
105 : {
106 0 : Array<Float> newPsf(newSize);
107 0 : newPsf = 0;
108 0 : newPsf(newblc,newtrc) = thePsf;
109 0 : thePsf.reference(newPsf);
110 0 : }
111 0 : thePsfOrigin = thePsfOrigin + newblc;
112 0 : result = thePsf(blc+newblc, trc+newblc);
113 0 : if (!nearAbs(Double(amplitude),Double(1.0)))
114 0 : result = thePsf(blc+newblc, trc+newblc)*amplitude;
115 0 : return true;
116 0 : }
117 : else {
118 0 : result = thePsf(blc, trc);
119 0 : if (!nearAbs(Double(amplitude),1.0))
120 0 : result = result * amplitude;
121 0 : return true;
122 : }
123 0 : };
124 :
125 0 : Bool ConvolutionEquation::
126 : residual(Array<Float> & result,
127 : Float & chisq,
128 : const LinearModel< Array<Float> > & model) {
129 0 : if (residual(result, model)) {
130 0 : chisq = sum(result*result);
131 0 : return true;
132 : }
133 : else
134 0 : return false;
135 : }
136 :
137 :
138 0 : Bool ConvolutionEquation::
139 : residual(Array<Float> & result,
140 : Float & chisq,
141 : Array<Float> & mask,
142 : const LinearModel< Array<Float> > & model) {
143 0 : if (residual(result, model)) {
144 0 : result = result * mask;
145 0 : chisq = sum(result*result);
146 0 : return true;
147 : }
148 : else
149 0 : return false;
150 : }
151 :
152 :
153 0 : Bool ConvolutionEquation::
154 : residual(Array<Float> & result,
155 : const LinearModel< Array<Float> > & model) {
156 0 : if (evaluate(result, model)) {
157 0 : result = theMeas - result;
158 0 : return true;
159 : }
160 : else
161 0 : return false;
162 : }
163 :
164 :
165 :
166 :
167 0 : Bool ConvolutionEquation::residual(Array<Float> & result,
168 : const LinearModel< MaskedArray<Float> > & model) {
169 0 : if (evaluate(result, model)) {
170 0 : result = theMeas - result;
171 0 : return true;
172 : }
173 : else
174 0 : return false;
175 : }
176 :
177 0 : Bool ConvolutionEquation::residual(MaskedArray<Float> & result,
178 : const LinearModel< MaskedArray<Float> > & model) {
179 0 : Array<Float> farray;
180 0 : if (residual(farray, model)) {
181 0 : result.setData(farray, model.getModel().getMask());
182 0 : return true;
183 : }
184 : else
185 0 : return false;
186 0 : }
187 :
188 0 : IPosition ConvolutionEquation::psfSize(){
189 0 : return theRealPsfSize;
190 : }
191 :
192 0 : void ConvolutionEquation::flushPsf(){
193 0 : thePsf.resize(IPosition(0));
194 0 : }
195 :
196 : } //# NAMESPACE CASA - END
197 :
|
