Line data Source code
1 : //# ConvolutionTVI.h: This file contains the implementation of the ConvolutionTVI class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <mstransform/TVI/ConvolutionTVI.h>
24 :
25 : using namespace casacore;
26 : namespace casa { //# NAMESPACE CASA - BEGIN
27 :
28 : namespace vi { //# NAMESPACE VI - BEGIN
29 :
30 : //////////////////////////////////////////////////////////////////////////
31 : // ConvolutionTVI class
32 : //////////////////////////////////////////////////////////////////////////
33 :
34 : // -----------------------------------------------------------------------
35 : //
36 : // -----------------------------------------------------------------------
37 0 : ConvolutionTVI::ConvolutionTVI( ViImplementation2 * inputVii,
38 0 : const Record &configuration):
39 0 : FreqAxisTVI (inputVii)
40 : {
41 : // Parse and check configuration parameters
42 : // Note: if a constructor finishes by throwing an exception, the memory
43 : // associated with the object itself is cleaned up — there is no memory leak.
44 0 : if (not parseConfiguration(configuration))
45 : {
46 0 : throw AipsError("Error parsing ConvolutionTVI configuration");
47 : }
48 :
49 0 : initialize();
50 :
51 0 : return;
52 : }
53 :
54 : // -----------------------------------------------------------------------
55 : //
56 : // -----------------------------------------------------------------------
57 0 : Bool ConvolutionTVI::parseConfiguration(const Record &configuration)
58 : {
59 0 : int exists = -1;
60 0 : Bool ret = true;
61 :
62 : // Parse kernel parameter (optional)
63 0 : exists = -1;
64 0 : exists = configuration.fieldNumber ("kernel");
65 0 : if (exists >= 0)
66 : {
67 0 : if( configuration.type(exists) == casacore::TpArrayFloat )
68 : {
69 0 : convCoeff_p.resize(0,false);
70 0 : convCoeff_p = configuration.asArrayFloat( exists );
71 0 : logger_p << LogIO::NORMAL << LogOrigin("ConvolutionTVI", __FUNCTION__)
72 0 : << "Kernel is " << convCoeff_p << LogIO::POST;
73 : }
74 : else
75 : {
76 0 : ret = false;
77 0 : logger_p << LogIO::SEVERE << LogOrigin("ConvolutionTVI", __FUNCTION__)
78 : << "Wrong format of kernel parameter (only float/double/int arrays are supported) "
79 0 : << LogIO::POST;
80 : }
81 : }
82 :
83 0 : return ret;
84 : }
85 :
86 : // -----------------------------------------------------------------------
87 : //
88 : // -----------------------------------------------------------------------
89 0 : void ConvolutionTVI::initialize()
90 : {
91 : // Populate nchan input-output maps
92 : Int spw;
93 0 : uInt spw_idx = 0;
94 0 : map<Int,vector<Int> >::iterator iter;
95 0 : for(iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
96 : {
97 0 : spw = iter->first;
98 0 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size();
99 :
100 0 : spw_idx++;
101 : }
102 :
103 0 : return;
104 : }
105 :
106 : // -----------------------------------------------------------------------
107 : //
108 : // -----------------------------------------------------------------------
109 0 : void ConvolutionTVI::flag(Cube<Bool>& flagCube) const
110 : {
111 : // Get input VisBuffer and SPW
112 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
113 :
114 : // Configure Transformation Engine
115 0 : ConvolutionLogicalORKernel<Bool> kernel(&convCoeff_p);
116 0 : ConvolutionTransformEngine<Bool> transformer(&kernel,convCoeff_p.size());
117 :
118 : // Transform data
119 0 : transformFreqAxis(vb->flagCube(),flagCube,transformer);
120 :
121 0 : return;
122 : }
123 :
124 : // -----------------------------------------------------------------------
125 : //
126 : // -----------------------------------------------------------------------
127 0 : void ConvolutionTVI::floatData (Cube<Float> & vis) const
128 : {
129 : // Get input VisBuffer and SPW
130 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
131 :
132 : // Configure Transformation Engine
133 0 : ConvolutionDataKernel<Float> kernel(&convCoeff_p);
134 0 : ConvolutionTransformEngine<Float> transformer(&kernel,convCoeff_p.size());
135 :
136 : // Transform data
137 0 : transformFreqAxis(vb->visCubeFloat(),vis,transformer);
138 :
139 0 : return;
140 : }
141 :
142 : // -----------------------------------------------------------------------
143 : //
144 : // -----------------------------------------------------------------------
145 0 : void ConvolutionTVI::visibilityObserved (Cube<Complex> & vis) const
146 : {
147 : // Get input VisBuffer and SPW
148 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
149 :
150 : // Configure Transformation Engine
151 0 : ConvolutionDataKernel<Complex> kernel(&convCoeff_p);
152 0 : ConvolutionTransformEngine<Complex> transformer(&kernel,convCoeff_p.size());
153 :
154 : // Transform data
155 0 : transformFreqAxis(vb->visCube(),vis,transformer);
156 :
157 0 : return;
158 : }
159 :
160 : // -----------------------------------------------------------------------
161 : //
162 : // -----------------------------------------------------------------------
163 0 : void ConvolutionTVI::visibilityCorrected (Cube<Complex> & vis) const
164 : {
165 : // Get input VisBuffer and SPW
166 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
167 :
168 : // Configure Transformation Engine
169 0 : ConvolutionDataKernel<Complex> kernel(&convCoeff_p);
170 0 : ConvolutionTransformEngine<Complex> transformer(&kernel,convCoeff_p.size());
171 :
172 : // Transform data
173 0 : transformFreqAxis(vb->visCubeCorrected(),vis,transformer);
174 :
175 0 : return;
176 : }
177 :
178 : // -----------------------------------------------------------------------
179 : //
180 : // -----------------------------------------------------------------------
181 0 : void ConvolutionTVI::visibilityModel (Cube<Complex> & vis) const
182 : {
183 : // Get input VisBuffer and SPW
184 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
185 :
186 : // Configure Transformation Engine
187 0 : ConvolutionDataKernel<Complex> kernel(&convCoeff_p);
188 0 : ConvolutionTransformEngine<Complex> transformer(&kernel,convCoeff_p.size());
189 :
190 : // Transform data
191 0 : transformFreqAxis(vb->visCubeModel(),vis,transformer);
192 :
193 0 : return;
194 : }
195 :
196 : // -----------------------------------------------------------------------
197 : //
198 : // -----------------------------------------------------------------------
199 0 : void ConvolutionTVI::weightSpectrum(Cube<Float> &weightSp) const
200 : {
201 : // Get input VisBuffer and SPW
202 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
203 :
204 : // Configure Transformation Engine
205 0 : ConvolutionWeightPropagationKernel<Float> kernel(&convCoeff_p);
206 0 : ConvolutionTransformEngine<Float> transformer(&kernel,convCoeff_p.size());
207 :
208 : // Transform data
209 0 : transformFreqAxis(vb->weightSpectrum(),weightSp,transformer);
210 :
211 0 : return;
212 : }
213 :
214 : // -----------------------------------------------------------------------
215 : //
216 : // -----------------------------------------------------------------------
217 0 : void ConvolutionTVI::sigmaSpectrum(Cube<Float> &sigmaSp) const
218 : {
219 : // Get input VisBuffer and SPW
220 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
221 :
222 : // Configure Transformation Engine
223 0 : ConvolutionWeightPropagationKernel<Float> kernel(&convCoeff_p);
224 0 : ConvolutionTransformEngine<Float> transformer(&kernel,convCoeff_p.size());
225 :
226 : // Get weightSpectrum from sigmaSpectrum
227 0 : Cube<Float> weightSpFromSigmaSp;
228 0 : weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),false);
229 0 : weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
230 0 : arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
231 :
232 : // Transform data
233 0 : transformFreqAxis(weightSpFromSigmaSp,sigmaSp,transformer);
234 :
235 : // Transform back from weight format to sigma format
236 0 : arrayTransformInPlace (sigmaSp,weightToSigma);
237 :
238 0 : return;
239 : }
240 :
241 : //////////////////////////////////////////////////////////////////////////
242 : // ConvolutionTVIFactory class
243 : //////////////////////////////////////////////////////////////////////////
244 :
245 : // -----------------------------------------------------------------------
246 : //
247 : // -----------------------------------------------------------------------
248 0 : ConvolutionTVIFactory::ConvolutionTVIFactory (Record &configuration,
249 0 : ViImplementation2 *inputVii)
250 : {
251 0 : inputVii_p = inputVii;
252 0 : configuration_p = configuration;
253 0 : }
254 :
255 : // -----------------------------------------------------------------------
256 : //
257 : // -----------------------------------------------------------------------
258 0 : vi::ViImplementation2 * ConvolutionTVIFactory::createVi(VisibilityIterator2 *) const
259 : {
260 0 : return new ConvolutionTVI(inputVii_p,configuration_p);
261 : }
262 :
263 : // -----------------------------------------------------------------------
264 : //
265 : // -----------------------------------------------------------------------
266 0 : vi::ViImplementation2 * ConvolutionTVIFactory::createVi() const
267 : {
268 0 : return new ConvolutionTVI(inputVii_p,configuration_p);
269 : }
270 :
271 : //////////////////////////////////////////////////////////////////////////
272 : // ConvolutionTransformEngine class
273 : //////////////////////////////////////////////////////////////////////////
274 :
275 : // -----------------------------------------------------------------------
276 : //
277 : // -----------------------------------------------------------------------
278 0 : template<class T> ConvolutionTransformEngine<T>::ConvolutionTransformEngine
279 : (ConvolutionKernel<T> *kernel,
280 0 : uInt width)
281 : {
282 0 : width_p = width;
283 0 : convolutionKernel_p = kernel;
284 0 : }
285 :
286 : // -----------------------------------------------------------------------
287 : //
288 : // -----------------------------------------------------------------------
289 0 : template<class T> void ConvolutionTransformEngine<T>::transform(Vector<T> &inputVector,
290 : Vector<T> &outputVector)
291 : {
292 0 : uInt startChanIndex = 0;
293 0 : uInt outChanStart = width_p / 2;
294 0 : uInt outChanIndex = outChanStart;
295 0 : uInt outChanStop = inputVector.size() - width_p / 2;
296 0 : while (startChanIndex < outChanStop)
297 : {
298 0 : convolutionKernel_p->kernel(inputVector,outputVector,startChanIndex,outChanIndex);
299 0 : startChanIndex += 1;
300 0 : outChanIndex += 1;
301 : }
302 :
303 : // Process low end
304 0 : for (uInt chanIndex = 0; chanIndex<outChanStart; chanIndex++)
305 : {
306 0 : convolutionKernel_p->kernel(inputVector,outputVector,chanIndex,chanIndex);
307 0 : chanIndex += 1;
308 : }
309 :
310 : // Process high end
311 0 : for (uInt chanIndex = outChanStop; chanIndex<inputVector.size(); chanIndex++)
312 : {
313 0 : convolutionKernel_p->kernel(inputVector,outputVector,chanIndex,chanIndex);
314 0 : chanIndex += 1;
315 : }
316 :
317 0 : return;
318 : }
319 :
320 : //////////////////////////////////////////////////////////////////////////
321 : // ConvolutionKernel class
322 : //////////////////////////////////////////////////////////////////////////
323 0 : template<class T> ConvolutionKernel<T>::ConvolutionKernel(Vector<Float> *convCoeff)
324 : {
325 0 : convCoeff_p = convCoeff;
326 0 : width_p = convCoeff->size();
327 0 : }
328 :
329 : //////////////////////////////////////////////////////////////////////////
330 : // ConvolutionDataKernel class
331 : //////////////////////////////////////////////////////////////////////////
332 :
333 : // -----------------------------------------------------------------------
334 : //
335 : // -----------------------------------------------------------------------
336 0 : template<class T> ConvolutionDataKernel<T>::ConvolutionDataKernel(Vector<Float> *convCoeff):
337 0 : ConvolutionKernel<T>(convCoeff)
338 : {
339 :
340 0 : }
341 :
342 : // -----------------------------------------------------------------------
343 : //
344 : // -----------------------------------------------------------------------
345 0 : template<class T> void ConvolutionDataKernel<T>::kernel( Vector<T> &inputVector,
346 : Vector<T> &outputVector,
347 : uInt startInputPos,
348 : uInt outputPos)
349 : {
350 : // Do not process edges
351 0 : if (startInputPos == outputPos)
352 : {
353 0 : outputVector(outputPos) = inputVector(startInputPos);
354 0 : return;
355 : }
356 :
357 : // Initialization
358 0 : outputVector(outputPos) = (*convCoeff_p)(0)*inputVector(startInputPos);
359 :
360 : // Main loop
361 0 : for (uInt chanIndex = 1; chanIndex<width_p; chanIndex++)
362 : {
363 0 : outputVector(outputPos) += (*convCoeff_p)(chanIndex)*inputVector(startInputPos+chanIndex);
364 : }
365 :
366 0 : return;
367 : }
368 :
369 : //////////////////////////////////////////////////////////////////////////
370 : // ConvolutionLogicalORKernel class
371 : //////////////////////////////////////////////////////////////////////////
372 :
373 : // -----------------------------------------------------------------------
374 : //
375 : // -----------------------------------------------------------------------
376 0 : template<class T> ConvolutionLogicalORKernel<T>::ConvolutionLogicalORKernel(Vector<Float> *convCoeff):
377 0 : ConvolutionKernel<T>(convCoeff)
378 : {
379 :
380 :
381 0 : }
382 :
383 : // -----------------------------------------------------------------------
384 : //
385 : // -----------------------------------------------------------------------
386 0 : template<class T> void ConvolutionLogicalORKernel<T>::kernel( Vector<T> &inputVector,
387 : Vector<T> &outputVector,
388 : uInt startInputPos,
389 : uInt outputPos)
390 : {
391 : // Flag edges
392 0 : if (startInputPos == outputPos)
393 : {
394 0 : outputVector(outputPos) = true;
395 0 : return;
396 : }
397 :
398 0 : Bool outputFlag = false;
399 : // Output sample is flagged if any of the contributors are flagged
400 0 : for (uInt chanIndex = 0; chanIndex<width_p; chanIndex++)
401 : {
402 0 : if (inputVector(startInputPos+chanIndex))
403 : {
404 0 : outputFlag = true;
405 0 : break;
406 : }
407 : }
408 :
409 0 : outputVector(outputPos) = outputFlag;
410 :
411 0 : return;
412 : }
413 :
414 : //////////////////////////////////////////////////////////////////////////
415 : // ConvolutionWeightPropagationKernel class
416 : //////////////////////////////////////////////////////////////////////////
417 :
418 : // -----------------------------------------------------------------------
419 : //
420 : // -----------------------------------------------------------------------
421 0 : template<class T> ConvolutionWeightPropagationKernel<T>::ConvolutionWeightPropagationKernel(Vector<Float> *convCoeff):
422 0 : ConvolutionKernel<T>(convCoeff)
423 : {
424 :
425 0 : }
426 :
427 : // -----------------------------------------------------------------------
428 : //
429 : // -----------------------------------------------------------------------
430 0 : template<class T> void ConvolutionWeightPropagationKernel<T>::kernel( Vector<T> &inputVector,
431 : Vector<T> &outputVector,
432 : uInt startInputPos,
433 : uInt outputPos)
434 : {
435 : // Do not process edges
436 0 : if (startInputPos == outputPos)
437 : {
438 0 : outputVector(outputPos) = inputVector(startInputPos);
439 0 : return;
440 : }
441 :
442 : // Initialization (mind for zeros as there is a division operation)
443 0 : outputVector(outputPos) = 0;
444 0 : if (inputVector(startInputPos) > FLT_MIN)
445 : {
446 0 : outputVector(outputPos) = (*convCoeff_p)(0)*(*convCoeff_p)(0)/inputVector(startInputPos);
447 : }
448 :
449 : // Main accumulation loop
450 0 : for (uInt chanIndex = 1; chanIndex<width_p; chanIndex++)
451 : {
452 : // Mind for zeros as there is a division operation
453 0 : if (inputVector(startInputPos+chanIndex) > FLT_MIN)
454 : {
455 0 : outputVector(outputPos) += (*convCoeff_p)(chanIndex)*(*convCoeff_p)(chanIndex)/inputVector(startInputPos+chanIndex);
456 : }
457 : }
458 :
459 : // Final propagated weight is the inverse of the accumulation
460 0 : if (outputVector(outputPos) > FLT_MIN)
461 : {
462 0 : outputVector(outputPos) = 1/outputVector(outputPos);
463 : }
464 :
465 0 : return;
466 : }
467 :
468 : } //# NAMESPACE VI - END
469 :
470 : } //# NAMESPACE CASA - END
471 :
472 :
|
