Line data Source code
1 : /*
2 : * SidebandSeparator.cc
3 : *
4 : * Created on: 2017/07/19
5 : * Author: kana
6 : */
7 :
8 : // STL
9 : #include <ctype.h>
10 :
11 : // cascore
12 : #include <casacore/casa/OS/File.h>
13 : #include <casacore/casa/Logging/LogIO.h>
14 :
15 : #include <imageanalysis/ImageAnalysis/ImageFactory.h>
16 : #include <imageanalysis/ImageAnalysis/ImageMetaData.h>
17 : #include <imageanalysis/ImageAnalysis/ImageMaskAttacher.h>
18 :
19 : // casa
20 : #include <synthesis/MeasurementEquations/SideBandSeparator.h>
21 :
22 : using namespace std ;
23 : using namespace casacore ;
24 :
25 : namespace casa {
26 :
27 : // constructors
28 0 : SideBandSeparatorBase::SideBandSeparatorBase(const vector<string>& inputname)
29 : {
30 0 : init();
31 0 : setInput(inputname);
32 : {// logging
33 0 : LogIO os(LogOrigin("SideBandSeparatorBase","SideBandSeparatorBase()", WHERE));
34 0 : os << "Found " << inputNames_.size() << " images." << LogIO::POST;
35 0 : os << "Input Data to be processed:" << LogIO::POST;
36 0 : for (size_t i = 0; i < inputNames_.size(); ++i) {
37 0 : os << "\t" << inputNames_[i] << LogIO::POST;
38 : }
39 0 : }
40 0 : };
41 :
42 0 : SideBandSeparatorBase::~SideBandSeparatorBase()
43 : {
44 0 : };
45 :
46 0 : void SideBandSeparatorBase::init()
47 : {
48 : // shifts
49 0 : initshift();
50 : // image list
51 0 : inputNames_.resize(0);
52 : // solution parameters
53 0 : otherside_ = false;
54 0 : doboth_ = false;
55 0 : rejlimit_ = 0.2;
56 0 : };
57 :
58 0 : void SideBandSeparatorBase::initshift()
59 : {
60 : // shifts
61 0 : nshift_ = 0;
62 0 : nchan_ = 0;
63 0 : sigShift_.resize(0);
64 0 : imgShift_.resize(0);
65 0 : };
66 :
67 0 : void SideBandSeparatorBase::setInput(const vector<string>& inputname) {
68 0 : inputNames_.resize(0);
69 : // check existence of images
70 0 : for (size_t i = 0 ; i < inputname.size(); ++i) {
71 0 : if (checkFile(inputname[i], "d")) {
72 0 : inputNames_.push_back(inputname[i]);
73 : } else {
74 0 : inputNames_.resize(0);
75 0 : throw( AipsError("Could not find "+inputname[i]) );
76 : }
77 : }
78 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setImage()", WHERE));
79 0 : if (inputNames_.size() < 2)
80 0 : throw( AipsError("At least two valid input data are required for processing") );
81 0 : }
82 :
83 0 : void SideBandSeparatorBase::setShift(const vector<double> &shift, const bool signal)
84 : {
85 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setShift()", WHERE));
86 0 : if (shift.size() != 0 && shift.size() != inputNames_.size())
87 0 : throw( AipsError("The number of shift should match that of images") );
88 0 : vector<double> *target = signal ? &sigShift_ : &imgShift_;
89 0 : target->resize(shift.size());
90 0 : for (unsigned int i = 0; i < shift.size(); i++)
91 0 : target->at(i) = - shift[i]; /// NOTE if spw shifts +3ch, need to shift back -3ch in operation.
92 :
93 0 : if (target->size() == 0) {
94 0 : os << "Channel shifts of " << (signal ? "SIGNAL" : "IMAGE") << " sideband are cleared." << LogIO::POST;
95 : } else {
96 0 : os << "Channel shifts of " << (signal ? "SIGNAL" : "IMAGE") << " sideband are set: ( ";
97 0 : for (unsigned int i = 0; i < target->size(); i++) {
98 0 : os << (signal ? sigShift_[i] : imgShift_[i]);
99 0 : if (i != target->size()-1) os << " , ";
100 : }
101 0 : os << " ) [channels]" << LogIO::POST;
102 : }
103 0 : };
104 :
105 0 : void SideBandSeparatorBase::setThreshold(const double limit)
106 : {
107 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setThreshold()", WHERE));
108 0 : if (limit <= 0 || limit >= 1.0)
109 0 : throw( AipsError("Rejection limit should be > 0.0 and < 1.0") );
110 :
111 0 : rejlimit_ = limit;
112 0 : os << "Rejection limit is set to " << rejlimit_ << LogIO::POST;
113 0 : };
114 :
115 : /////////////// PROTECTED FUNCTIONS //////////////////////
116 :
117 0 : size_t SideBandSeparatorBase::setupShift()
118 : {
119 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setupShift()", WHERE));
120 0 : if (sigShift_.size() > 0 && imgShift_.size() > 0) {
121 0 : return sigShift_.size();
122 0 : } else if (sigShift_.size() > 0 && imgShift_.size() == 0) {
123 0 : os << "Channel shift set only for signal sideband. Assuming the same shift in the opposite direction." << LogIO::POST;
124 0 : imgShift_.resize(sigShift_.size());
125 0 : for (size_t i = 0; i < sigShift_.size(); ++i) {
126 0 : imgShift_[i] = - sigShift_[i];
127 : }
128 0 : } else if (sigShift_.size() == 0 && imgShift_.size() > 0) {
129 0 : os << "Channel shift set only for image sideband. Assuming the same shift in the opposite direction." << LogIO::POST;
130 0 : sigShift_.resize(imgShift_.size());
131 0 : for (size_t i = 0; i < imgShift_.size(); ++i) {
132 0 : sigShift_[i] = - imgShift_[i];
133 : }
134 : } else {
135 0 : throw( AipsError("Channel shift was not been set.") );
136 : }
137 0 : return sigShift_.size();
138 0 : }
139 :
140 0 : Vector<bool> SideBandSeparatorBase::collapseMask(const Matrix<bool> &flagMat,
141 : const vector<uInt> &inIdvec,
142 : const bool signal)
143 : {
144 0 : LogIO os(LogOrigin("SideBandSeparatorBase","collapseFlag()", WHERE));
145 0 : if (inIdvec.size() == 0)
146 0 : throw(AipsError("Internal error. Table index is not defined."));
147 0 : if (flagMat.ncolumn() != inIdvec.size())
148 0 : throw(AipsError("Internal error. The row number of input matrix is not conformant."));
149 0 : if (flagMat.nrow() != nchan_)
150 0 : throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
151 :
152 0 : const size_t nspec = inIdvec.size();
153 : vector<double> *thisShift;
154 0 : if (signal == otherside_) {
155 : // (solve signal && solveother = T) OR (solve image && solveother = F)
156 0 : thisShift = &imgShift_;
157 : } else {
158 : // (solve signal && solveother = F) OR (solve image && solveother = T)
159 0 : thisShift = &sigShift_;
160 : }
161 :
162 0 : Vector<bool> outflag(nchan_, true);
163 : double tempshift;
164 0 : Vector<bool> shiftvec(nchan_, true);
165 0 : Vector<bool> accflag(nchan_, true);
166 : uInt shiftId;
167 0 : for (uInt i = 0 ; i < nspec; ++i) {
168 0 : shiftId = inIdvec[i];
169 0 : tempshift = - thisShift->at(shiftId);
170 0 : shiftFlag(flagMat.column(i), tempshift, shiftvec);
171 : // Now accumulate Mask (true only if all data is valid)
172 0 : for (uInt j = 0 ; j < nchan_ ; ++j)
173 : // accflag[j] |= shiftvec[j];
174 0 : accflag[j] &= shiftvec[j];
175 : }
176 0 : outflag = accflag;
177 : // Shift back Flag
178 : //cout << "Shifting FLAG back to " << thisShift->at(0) << " channels" << endl;
179 : //shiftFlag(accflag, thisShift->at(0), outflag);
180 :
181 0 : return outflag;
182 0 : }
183 :
184 :
185 0 : Vector<float> SideBandSeparatorBase::solve(const Matrix<float> &specmat,
186 : const vector<uInt> &inIdvec,
187 : const bool signal)
188 : {
189 0 : LogIO os(LogOrigin("SideBandSeparatorBase","solve()", WHERE));
190 0 : if (inIdvec.size() == 0)
191 0 : throw(AipsError("Internal error. Table index is not defined."));
192 0 : if (specmat.ncolumn() != inIdvec.size())
193 0 : throw(AipsError("Internal error. The row number of input matrix is not conformant."));
194 0 : if (specmat.nrow() != nchan_)
195 0 : throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
196 :
197 :
198 : os << LogIO::DEBUG1 << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband."
199 0 : << LogIO::POST;
200 :
201 0 : const size_t nspec = inIdvec.size();
202 : vector<double> *thisShift, *otherShift;
203 0 : if (signal == otherside_) {
204 : // (solve signal && solveother = T) OR (solve image && solveother = F)
205 0 : thisShift = &imgShift_;
206 0 : otherShift = &sigShift_;
207 0 : os << LogIO::DEBUG1 << "Image sideband will be deconvolved." << LogIO::POST;
208 : } else {
209 : // (solve signal && solveother = F) OR (solve image && solveother = T)
210 0 : thisShift = &sigShift_;
211 0 : otherShift = &imgShift_;
212 0 : os << LogIO::DEBUG1 << "Signal sideband will be deconvolved." << LogIO::POST;
213 : }
214 :
215 0 : vector<double> spshift(nspec);
216 0 : Matrix<float> shiftSpecmat(nchan_, nspec, 0.);
217 : double tempshift;
218 0 : Vector<float> shiftspvec;
219 : uInt shiftId;
220 0 : for (uInt i = 0 ; i < nspec; i++) {
221 0 : shiftId = inIdvec[i];
222 0 : spshift[i] = otherShift->at(shiftId) - thisShift->at(shiftId);
223 0 : tempshift = - thisShift->at(shiftId);
224 0 : shiftspvec.reference(shiftSpecmat.column(i));
225 0 : shiftSpectrum(specmat.column(i), tempshift, shiftspvec);
226 : }
227 :
228 0 : Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.);
229 0 : vector<float> thisvec(nchan_, 0.);
230 :
231 : float minval, maxval;
232 : {//Debug
233 0 : minMax(minval, maxval, shiftSpecmat);
234 : os << LogIO::DEBUG1 << "Max/Min of input Matrix = (max: " << maxval
235 0 : << ", min: " << minval << ")" << LogIO::POST;
236 : }
237 :
238 0 : os << LogIO::DEBUG1 << "starting deconvolution" << LogIO::POST;
239 0 : deconvolve(shiftSpecmat, spshift, rejlimit_, convmat);
240 :
241 : {//Debug
242 0 : minMax(minval, maxval, convmat);
243 : os << LogIO::DEBUG1 << "Max/Min of output Matrix = (max: " << maxval
244 0 : << ", min: " << minval << ")" << LogIO::POST;
245 : }
246 :
247 0 : aggregateMat(convmat, thisvec);
248 :
249 0 : if (!otherside_) return Vector<float>(thisvec);
250 :
251 : // subtract from the other side band.
252 0 : vector<float> othervec(nchan_, 0.);
253 0 : subtractFromOther(shiftSpecmat, thisvec, spshift, othervec);
254 0 : return Vector<float>(othervec);
255 0 : };
256 :
257 :
258 0 : void SideBandSeparatorBase::shiftSpectrum(const Vector<float> &invec,
259 : double shift,
260 : Vector<float> &outvec)
261 : {
262 0 : LogIO os(LogOrigin("SideBandSeparatorBase","shiftSpectrum()", WHERE));
263 0 : if (invec.size() != nchan_)
264 0 : throw(AipsError("Internal error. The length of input vector differs from nchan_"));
265 0 : if (outvec.size() != nchan_)
266 0 : throw(AipsError("Internal error. The length of output vector differs from nchan_"));
267 :
268 : // cout << "Start shifting spectrum for " << shift << " channels" << endl;
269 :
270 : // tweak shift to be in 0 ~ nchan_-1
271 0 : if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
272 0 : if (shift < 0.) shift += nchan_;
273 0 : double rweight = fmod(shift, 1.);
274 0 : if (rweight < 0.) rweight += 1.;
275 0 : double lweight = 1. - rweight;
276 : uInt lchan, rchan;
277 :
278 0 : outvec = 0;
279 0 : for (uInt i = 0 ; i < nchan_ ; i++) {
280 0 : lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
281 : if (lchan < 0.) lchan += nchan_;
282 0 : rchan = ( (lchan + 1) % nchan_ );
283 0 : outvec(lchan) += invec(i) * lweight;
284 0 : outvec(rchan) += invec(i) * rweight;
285 :
286 : // if (i == 2350 || i== 2930) {
287 : // cout << "Channel= " << i << " of input vector: " << endl;
288 : // cout << "L channel = " << lchan << endl;
289 : // cout << "R channel = " << rchan << endl;
290 : // cout << "L weight = " << lweight << endl;
291 : // cout << "R weight = " << rweight << endl;
292 : // }
293 : }
294 0 : };
295 :
296 :
297 0 : void SideBandSeparatorBase::shiftFlag(const Vector<bool> &invec,
298 : double shift,
299 : Vector<bool> &outvec)
300 : {
301 0 : LogIO os(LogOrigin("SideBandSeparatorBase","shiftFlag()", WHERE));
302 0 : if (invec.size() != nchan_)
303 0 : throw(AipsError("Internal error. The length of input vector differs from nchan_"));
304 0 : if (outvec.size() != nchan_)
305 0 : throw(AipsError("Internal error. The length of output vector differs from nchan_"));
306 :
307 : // cout << "Start shifting flag for " << shift << "channels" << endl;
308 :
309 : // shift is almost integer think it as int.
310 : // tolerance should be in 0 - 1
311 0 : double tolerance = 0.01;
312 : // tweak shift to be in 0 ~ nchan_-1
313 0 : if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
314 0 : if (shift < 0.) shift += nchan_;
315 0 : double rweight = fmod(shift, 1.);
316 0 : bool ruse(true), luse(true);
317 0 : if (rweight < 0.) rweight += 1.;
318 0 : if (rweight < tolerance){
319 : // the shift is almost lchan
320 0 : ruse = false;
321 0 : luse = true;
322 : }
323 0 : if (rweight > 1-tolerance){
324 : // the shift is almost rchan
325 0 : ruse = true;
326 0 : luse = false;
327 : }
328 : uInt lchan, rchan;
329 :
330 0 : outvec = false;
331 0 : for (uInt i = 0 ; i < nchan_ ; i++) {
332 0 : lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
333 : if (lchan < 0.) lchan += nchan_;
334 0 : rchan = ( (lchan + 1) % nchan_ );
335 0 : outvec(lchan) |= (invec(i) && luse);
336 0 : outvec(rchan) |= (invec(i) && ruse);
337 :
338 : // if (i == 2350 || i == 2930) {
339 : // cout << "Channel= " << i << " of input vector: " << endl;
340 : // cout << "L channel = " << lchan << endl;
341 : // cout << "R channel = " << rchan << endl;
342 : // cout << "L channel will be " << (luse ? "used" : "ignored") << endl;
343 : // cout << "R channel will be " << (ruse ? "used" : "ignored") << endl;
344 : // }
345 : }
346 0 : };
347 :
348 :
349 0 : void SideBandSeparatorBase::deconvolve(Matrix<float> &specmat,
350 : const vector<double> shiftvec,
351 : const double threshold,
352 : Matrix<float> &outmat)
353 : {
354 0 : LogIO os(LogOrigin("SideBandSeparatorBase","deconvolve()", WHERE));
355 0 : if (specmat.nrow() != nchan_)
356 0 : throw(AipsError("Internal error. The length of input matrix differs from nchan_"));
357 0 : if (specmat.ncolumn() != shiftvec.size())
358 0 : throw(AipsError("Internal error. The number of input shifts and spectrum differs."));
359 :
360 : float minval, maxval;
361 : {//Debug
362 0 : minMax(minval, maxval, specmat);
363 : os << LogIO::DEBUG1 << "Max/Min of input Matrix = (max: " << maxval
364 0 : << ", min: " << minval << ")" << LogIO::POST;
365 : }
366 :
367 0 : uInt ninsp = shiftvec.size();
368 0 : outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.);
369 0 : Matrix<Complex> fftspmat(nchan_/2+1, ninsp, Complex(0.));
370 0 : Vector<float> rvecref(nchan_, 0.);
371 0 : Vector<Complex> cvecref(nchan_/2+1, Complex(0.));
372 0 : uInt icol = 0;
373 0 : unsigned int nreject = 0;
374 :
375 0 : os << LogIO::DEBUG1 << "Starting initial FFT. The number of input spectra = " << ninsp << LogIO::POST;
376 0 : os << LogIO::DEBUG1 << "out matrix has ncolumn = " << outmat.ncolumn() << LogIO::POST;
377 :
378 0 : for (uInt isp = 0 ; isp < ninsp ; isp++) {
379 0 : rvecref.reference( specmat.column(isp) );
380 0 : cvecref.reference( fftspmat.column(isp) );
381 :
382 : {//Debug
383 0 : minMax(minval, maxval, rvecref);
384 : os << LogIO::DEBUG1 << "Max/Min of inv FFTed Matrix = (max: " << maxval
385 0 : << ", min: " << minval << ")" << LogIO::POST;
386 : }
387 :
388 0 : fftsf.fft0(cvecref, rvecref, true);
389 : }
390 :
391 : //Liberate from reference
392 0 : rvecref.unique();
393 :
394 0 : Vector<Complex> cspec(nchan_/2+1, Complex(0.));
395 0 : const double PI = 6.0 * asin(0.5);
396 0 : const double nchani = 1. / (float) nchan_;
397 0 : const Complex trans(0., 1.);
398 :
399 0 : os << LogIO::DEBUG1 << "starting actual deconvolution" << LogIO::POST;
400 0 : for (uInt j = 0 ; j < ninsp ; j++) {
401 0 : for (uInt k = j+1 ; k < ninsp ; k++) {
402 0 : const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani;
403 :
404 0 : for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){
405 0 : cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5;
406 0 : double phase = dx*ichan;
407 0 : if ( fabs( sin(phase) ) > threshold){
408 0 : cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5
409 0 : * trans * cos(0.5*phase) / sin(0.5*phase);
410 : // * trans * sin(phase) / ( 1. - cos(phase) );
411 : } else {
412 0 : nreject++;
413 : }
414 : } // end of channel loop
415 : // os << LogIO::DEBUG1 << "done calculation of cspec" << LogIO::POST;
416 :
417 0 : Vector<Float> rspec;
418 0 : rspec.reference( outmat.column(icol) );
419 :
420 : // os << LogIO::DEBUG1 << "Starting inverse FFT. icol = " << icol << LogIO::POST;
421 0 : fftsi.fft0(rspec, cspec, false);
422 :
423 0 : icol++;
424 0 : }
425 : }
426 :
427 : {//Debug
428 0 : minMax(minval, maxval, outmat);
429 0 : os << LogIO::DEBUG1 << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << LogIO::POST;
430 : }
431 :
432 : os << LogIO::DEBUG1 << "Threshold = " << threshold
433 0 : << ", Rejected channels = " << nreject << LogIO::POST;
434 0 : };
435 :
436 :
437 0 : void SideBandSeparatorBase::aggregateMat(const Matrix<float> &inmat,
438 : vector<float> &outvec)
439 : {
440 0 : LogIO os(LogOrigin("SideBandSeparatorBase","aggregateMat()", WHERE));
441 0 : if (inmat.nrow() != nchan_)
442 0 : throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
443 : // if (outvec.size() != nchan_)
444 : // throw(AipsError("Internal error. The size of output vector should be equal to nchan_"));
445 :
446 : os << LogIO::DEBUG1 << "Averaging " << inmat.ncolumn() << " spectra in the input matrix."
447 0 : << LogIO::POST;
448 :
449 0 : const uInt nspec = inmat.ncolumn();
450 0 : const double scale = 1./( (double) nspec );
451 : // initialize values with 0
452 0 : outvec.assign(nchan_, 0);
453 0 : for (uInt isp = 0 ; isp < nspec ; isp++) {
454 0 : for (uInt ich = 0 ; ich < nchan_ ; ich++) {
455 0 : outvec[ich] += inmat(ich, isp);
456 : }
457 : }
458 :
459 0 : vector<float>::iterator iter;
460 0 : for (iter = outvec.begin(); iter != outvec.end(); iter++){
461 0 : *iter *= scale;
462 : }
463 0 : };
464 :
465 0 : void SideBandSeparatorBase::subtractFromOther(const Matrix<float> &shiftmat,
466 : const vector<float> &invec,
467 : const vector<double> &shift,
468 : vector<float> &outvec)
469 : {
470 0 : LogIO os(LogOrigin("SideBandSeparatorBase","subtractFromOther()", WHERE));
471 0 : if (shiftmat.nrow() != nchan_)
472 0 : throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
473 0 : if (invec.size() != nchan_)
474 0 : throw(AipsError("Internal error. The length of input vector should be nchan_"));
475 0 : if (shift.size() != shiftmat.ncolumn())
476 0 : throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift"));
477 :
478 0 : const uInt nspec = shiftmat.ncolumn();
479 0 : Vector<float> subsp(nchan_, 0.), shiftsub;
480 0 : Matrix<float> submat(nchan_, nspec, 0.);
481 0 : for (uInt isp = 0 ; isp < nspec ; isp++) {
482 0 : for (uInt ich = 0; ich < nchan_ ; ich++) {
483 0 : subsp(ich) = shiftmat(ich, isp) - invec[ich];
484 : }
485 0 : shiftsub.reference(submat.column(isp));
486 : // shift back spectrum so that other side stays still.
487 0 : shiftSpectrum(subsp, -shift[isp], shiftsub);
488 : }
489 :
490 0 : aggregateMat(submat, outvec);
491 0 : };
492 :
493 0 : bool SideBandSeparatorBase::interpolateMaskedChannels(Array<float> spectrum,
494 : const Array<bool> mask)
495 : {
496 : // if (spectrum.ndim() != 1 || mask.ndim() != 1)
497 : // throw AipsError("Array arguments must be 1-dimension.");
498 0 : if (spectrum.size() != mask.size())
499 0 : throw AipsError("Size mismatch between spectrum and mask.");
500 0 : size_t num_chan = spectrum.size();
501 0 : float* spectrum_p = spectrum.data();
502 0 : const bool* mask_p = mask.data();
503 : // get the first and the last valid channel IDs.
504 0 : size_t ledge=0, redge=num_chan-1;
505 0 : while (!mask_p[ledge] && ledge < num_chan-1) {
506 0 : ++ledge;
507 : }
508 0 : while (!mask_p[redge] && redge > 0) {
509 0 : --redge;
510 : }
511 : // Return if no valid channel
512 0 : if (redge < ledge) return false;
513 : // Nearest-neighbor for edges;
514 0 : for (size_t i = 0; i < ledge; ++i) spectrum_p[i] = spectrum_p[ledge];
515 0 : for (size_t i = num_chan-1; i > redge; --i) spectrum_p[i] = spectrum_p[redge];
516 : // Linear interpolation for intermediate gaps
517 0 : size_t mstart = -1, mend = -1, i0;
518 : float slope;
519 0 : for (size_t i = ledge+1; i < redge; ++i) {
520 0 : if (!mask_p[i]) {
521 : // the first masked channel
522 0 : mstart=i;
523 : // the last valid channel
524 0 : i0 = mstart-1;
525 : // search for the end of masked channel range
526 0 : while(!mask_p[i] && i < redge) {
527 0 : mend = i;
528 0 : ++i;
529 : }
530 : // 'mend' points to the last masked channel
531 : // while 'i' points to the next valid channel
532 0 : slope = (spectrum_p[i] - spectrum_p[i0])/static_cast<float>(i-i0);
533 0 : for (size_t j = mstart; j < mend+1; ++j) {
534 0 : spectrum_p[j] = slope*static_cast<float>(j-i0) + spectrum_p[i0];
535 : }
536 0 : mstart = -1;
537 0 : mend = -1;
538 : }
539 : }
540 0 : return true;
541 : }
542 :
543 0 : Bool SideBandSeparatorBase::checkFile(const string name, string type)
544 : {
545 0 : File file(name);
546 0 : if (!file.exists()){
547 0 : return false;
548 0 : } else if (type.empty()) {
549 0 : return true;
550 : } else {
551 : // Check for file type
552 0 : switch (std::tolower(type[0])) {
553 0 : case 'f':
554 0 : return file.isRegular(True);
555 0 : case 'd':
556 0 : return file.isDirectory(True);
557 0 : case 's':
558 0 : return file.isSymLink();
559 0 : default:
560 0 : throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
561 : }
562 : }
563 0 : };
564 :
565 : ///////////////////////////////////////////////////////////////////////////////////
566 0 : SideBandSeparatorII::SideBandSeparatorII(const std::vector<std::string>& imagename)
567 0 : : SideBandSeparatorBase (imagename)
568 : {
569 0 : initLocal();
570 0 : checkImage();
571 0 : };
572 :
573 0 : void SideBandSeparatorII::initLocal() {
574 0 : refFreqPix_ = 0.0;
575 0 : refFreqHz_ = -1.0;
576 0 : }
577 :
578 0 : void SideBandSeparatorII::checkImage() {
579 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setImage()", WHERE));
580 :
581 : // check image axes (npix, incr, ref, ndim)
582 0 : os << "Image axes check. Using the first image as the template" << LogIO::POST;
583 0 : CoordinateSystem csys0;
584 0 : IPosition npix0;
585 0 : if (!getImageCoordinate(inputNames_[0], csys0, npix0)) {
586 0 : throw( AipsError("Invalid image "+inputNames_[0]) );
587 : }
588 : // Summary
589 0 : os << "Template image coordinate:" << LogIO::POST;
590 0 : os << "\tndim\t" << csys0.nCoordinates() << LogIO::POST;
591 0 : os << "\tAxes\t" << csys0.worldAxisNames() << LogIO::POST;
592 0 : os << "\tnPix\t" << npix0 << LogIO::POST;
593 0 : os << "\tRefPix\t" << csys0.referencePixel() << LogIO::POST;
594 0 : os << "\tRefPValt" << csys0.referenceValue() << LogIO::POST;
595 0 : os << "\tIncr\t" << csys0.increment() << LogIO::POST;
596 :
597 0 : for (size_t i = 1; i < inputNames_.size(); ++i){
598 0 : if(!compareImageAxes(inputNames_[i], csys0, npix0))
599 0 : throw( AipsError("Image axes mismatch: "+inputNames_[0]) );
600 : }
601 0 : }
602 :
603 0 : bool SideBandSeparatorII::getImageCoordinate(const string& imagename, CoordinateSystem &csys, IPosition &npix) {
604 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setImage()", WHERE));
605 0 : auto ret = ImageFactory::fromFile(imagename);
606 0 : auto imageF = std::get<0>(ret);
607 0 : if (imageF) { //float image
608 0 : os << "Found float image" << LogIO::POST;
609 0 : npix = imageF->shape();
610 0 : ImageMetaData<Float> immd(imageF);
611 0 : vector<Int> myAxes;
612 0 : csys = immd.coordsys(myAxes);
613 0 : return true;
614 0 : } else if (std::get<1>(ret)) { // complex image
615 0 : os << "Found complex image" << LogIO::POST;
616 0 : return false;
617 0 : } else if (std::get<2>(ret)) { // double image
618 0 : os << "Found double image" << LogIO::POST;
619 0 : return false;
620 0 : } else if (std::get<3>(ret)) { // complex double image
621 0 : os << "Found complex double image" << LogIO::POST;
622 0 : return false;
623 : } else {
624 0 : os << LogIO::WARN << "Failed to open " << imagename << LogIO::POST;
625 0 : return false;
626 : }
627 0 : }
628 :
629 0 : bool SideBandSeparatorII::compareImageAxes(const string& imagename, const CoordinateSystem &refcsys, const IPosition &refnpix)
630 : {
631 0 : CoordinateSystem csys;
632 0 : IPosition npix;
633 0 : if (!getImageCoordinate(imagename, csys, npix)) {
634 0 : throw( AipsError("Invalid image "+imagename) );
635 : }
636 0 : uInt ndim = refcsys.nWorldAxes();
637 0 : if (csys.nWorldAxes() != ndim || refnpix.size() != ndim || npix.size() != ndim) {
638 0 : return false;
639 : }
640 0 : bool match = true;
641 0 : constexpr Double frac_tol = 0.1;
642 0 : for (uInt i = 0; i<refcsys.nCoordinates(); ++i) {
643 0 : Double tolerance = frac_tol*abs(refcsys.increment()[i]);
644 0 : match &= (npix[i]==refnpix[i]); // npix
645 0 : match &= (abs(csys.increment()[i]-refcsys.increment()[i]) < tolerance); // incr
646 0 : if (refcsys.type(i) == Coordinate::SPECTRAL) continue; // skip channel comparison
647 0 : match &= (abs(csys.referencePixel()[i]-refcsys.referencePixel()[i]) < frac_tol); // refpix
648 0 : match &= (abs(csys.referenceValue()[i]-refcsys.referenceValue()[i]) < tolerance); // refval
649 : }
650 0 : return match;
651 0 : }
652 :
653 0 : void SideBandSeparatorII::setImageBandFrequency(const double refpix, const Quantity refval)
654 : {
655 0 : LogIO os(LogOrigin("SideBandSeparatorBase","setImageBandFrequency()", WHERE));
656 0 : double freq_hz = refval.getValue("Hz", True);
657 0 : if (freq_hz < 0.0) throw( AipsError("Frequency should be positive") );
658 0 : refFreqPix_ = refpix;
659 0 : refFreqHz_ = freq_hz;
660 0 : os << "Setting frequency axis of image band: " << refFreqHz_*1.e-9
661 0 : << "GHz at channel " << refFreqPix_ << LogIO::POST;
662 0 : }
663 :
664 0 : void SideBandSeparatorII::separate(const string& outfile, const bool overwrite)
665 : {
666 0 : LogIO os(LogOrigin("SideBandSeparatorBase","separate()", WHERE));
667 : /*** Check outputs ***/
668 0 : if (outfile.size() == 0)
669 0 : throw( AipsError("Output file name is undefined."));
670 0 : string const signame = outfile + ".signalband";
671 0 : string const imgname = outfile + ".imageband";
672 0 : if (checkFile(signame, "d") && !overwrite) {
673 0 : throw( AipsError("Image "+signame+" already exists.") );
674 : }
675 0 : if (doboth_ && checkFile(imgname, "d") && !overwrite) {
676 0 : throw( AipsError("Image "+imgname+" already exists.") );
677 : }
678 :
679 : /*** Set up channel shift of image and signal sideband ***/
680 0 : nshift_ = setupShift();
681 0 : if (nshift_ < 2)
682 0 : throw( AipsError("At least 2 IFs are necessary for convolution.") );
683 0 : if (nshift_ != inputNames_.size())
684 0 : throw( AipsError("Internal error: nshift_ and image number differs.") );
685 :
686 : /*** Open input images (data model dependent) ***/
687 0 : vector<SPIIF> images(nshift_);
688 0 : for (size_t i = 0; i < nshift_; ++i) {
689 0 : auto ret = ImageFactory::fromFile(inputNames_[i]);
690 0 : auto imageF = std::get<0>(ret);
691 0 : if (! imageF)
692 0 : throw( AipsError("Float image not found in "+inputNames_[i]) );
693 0 : images[i] = imageF;
694 0 : }
695 : /*** analyze axis of reference image (data model dependent) ***/
696 0 : SPIIF refImage = images[0];
697 0 : IPosition const npix = refImage->shape();
698 0 : uInt const ndim = npix.size();
699 0 : ImageMetaData<Float> const immd(refImage);
700 0 : vector<Int> myAxes;
701 0 : CoordinateSystem const csys = immd.coordsys(myAxes);
702 0 : if (!csys.hasSpectralAxis())
703 0 : throw( AipsError("Could not find spectral axis.") );
704 0 : Int spax_id = csys.spectralAxisNumber();
705 0 : nchan_ = npix[spax_id];
706 0 : IPosition specArrayShape(ndim, 1), specVecShape(1, nchan_);
707 0 : specArrayShape[spax_id] = nchan_;
708 :
709 : /*** prepare output image (data model dependent) ***/
710 0 : SPIIF sigImg, imgImg;
711 0 : sigImg = ImageFactory::createImage<Float>(signame, csys, npix, True, overwrite, nullptr);
712 0 : if (doboth_) {
713 : /*** define frequency coordinate of image sideband ***/
714 0 : CoordinateSystem imgcsys(csys);
715 0 : if (refFreqHz_ < 0.0) {
716 0 : os << LogIO::WARN << "Frequency information of image sideband is not set. Using the value of signal sideband" << LogIO::POST;
717 : } else {
718 0 : Int spax_id = imgcsys.spectralAxisNumber();
719 0 : Vector<Double> refpix = imgcsys.referencePixel();
720 0 : Vector<Double> refval = imgcsys.referenceValue();
721 0 : Vector<Double> incr = imgcsys.increment();
722 0 : refpix[spax_id] = refFreqPix_;
723 0 : refval[spax_id] = refFreqHz_;
724 0 : incr[spax_id] *= -1.0;
725 0 : imgcsys.setReferencePixel(refpix);
726 0 : imgcsys.setReferenceValue(refval);
727 0 : imgcsys.setIncrement(incr);
728 0 : }
729 0 : imgImg = ImageFactory::createImage<Float>(imgname, imgcsys, npix, True, overwrite, nullptr);
730 0 : }
731 :
732 0 : Matrix<float> specMat(nchan_, nshift_);
733 0 : Matrix<bool> maskMat(nchan_, nshift_);
734 0 : Array<float> sigSpec(specVecShape), imgSpec(specVecShape);
735 0 : Array<bool> sigMask(specVecShape), imgMask(specVecShape);
736 0 : vector<uInt> dataIdvec;
737 :
738 : /*** Generate FFTServer ***/
739 0 : fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX);
740 0 : fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL);
741 :
742 : /*** Loop over image pixel and separate sideband (data model dependent) ***/
743 0 : IPosition start(ndim), end(ndim), stride(ndim);
744 0 : for (int ipos = 0; ipos < npix.product()/nchan_; ipos++){
745 0 : dataIdvec.resize(0);
746 : /*** convert 1-D ipos to slicer in image coordinate ***/
747 0 : uInt denominator = 1;
748 0 : for (uInt i = 0; i < ndim; ++i) {
749 0 : if ((Int) i == spax_id) {
750 0 : start(i) = 0;
751 0 : end(i) = nchan_-1;
752 0 : stride(i) = 1;
753 0 : denominator *= nchan_;
754 : } else {
755 0 : uInt pos = ((ipos/denominator) % npix[i]);
756 0 : start(i) = pos;
757 0 : end(i) = pos;
758 0 : stride(i) = npix[i];
759 0 : denominator *= npix[i];
760 : }
761 : }
762 0 : Slicer const specSlicer(start, end, stride, Slicer::endIsLast);
763 : /*** Get a set of spectra to solve (data model dependent) ***/
764 0 : if (!getSpectraToSolve(images, specSlicer, specMat, maskMat, dataIdvec)){
765 0 : sigSpec = 0.0f;
766 0 : imgSpec = 0.0f;
767 0 : sigMask = false;
768 0 : imgMask = false;
769 : } else {
770 : /*** Solve signal sideband ***/
771 0 : sigSpec = solve(specMat, dataIdvec, true);
772 : /*** Aggregate channel mask ***/
773 0 : sigMask = collapseMask(maskMat, dataIdvec, true);
774 0 : if (doboth_) {
775 : /*** Solve image sideband ***/
776 0 : imgSpec = solve(specMat, dataIdvec, false);
777 : /*** Aggregate channel mask ***/
778 0 : imgMask = collapseMask(maskMat, dataIdvec, false);
779 : }
780 : }
781 : /*** Signal side band ***/
782 : /*** now assign spec and mask to output image (data model dependent) ***/
783 0 : sigImg->putSlice(sigSpec.reform(specArrayShape), start, stride);
784 : /*** need to create a pixel mask if not any (data model dependent) ***/
785 0 : if (! sigImg->hasPixelMask()) {
786 0 : casacore::String maskName("");
787 0 : ImageMaskAttacher::makeMask(*sigImg, maskName, true, true, os, true);
788 0 : }
789 0 : sigImg->pixelMask().putSlice(sigMask.reform(specArrayShape), start, stride);
790 : /*** Image side band ***/
791 0 : if (doboth_) {
792 : /*** now assign spec and mask to output image (data model dependent) ***/
793 0 : imgImg->putSlice(imgSpec.reform(specArrayShape), start, stride);
794 : /*** need to create a pixel mask if not any (data model dependent) ***/
795 0 : if (! imgImg->hasPixelMask()) {
796 0 : casacore::String maskName("");
797 0 : ImageMaskAttacher::makeMask(*imgImg, maskName, true, true, os, true);
798 0 : }
799 0 : imgImg->pixelMask().putSlice(imgMask.reform(specArrayShape), start, stride);
800 : }
801 0 : } // end of row loop
802 : /*** Finally, save tables on disk (data model dependent) ***/
803 0 : sigImg.reset();
804 0 : imgImg.reset();
805 :
806 0 : }
807 :
808 0 : bool SideBandSeparatorII::getSpectraToSolve(const vector<SPIIF> &images, const Slicer &slicer,
809 : Matrix<float>& specMat, Matrix<bool>& maskMat, vector<uInt>& imgIdvec){
810 0 : imgIdvec.resize(0);
811 0 : specMat.resize(nchan_, nshift_);
812 0 : maskMat.resize(nchan_, nshift_);
813 0 : Array<float> spec(IPosition(1,nchan_));
814 0 : Array<bool> mask(IPosition(1,nchan_));
815 0 : size_t nspec = 0;
816 0 : for (size_t i = 0; i < images.size(); ++i) {
817 0 : images[i]->getSlice(spec, slicer, False);
818 0 : images[i]->getMaskSlice(mask, slicer, False);
819 : // Do interpolation of masked channels.
820 : // The method return false if there is no valid data (mask is true for valid pixels)
821 : // Skip if no valid channel in this spectrum
822 0 : if (!interpolateMaskedChannels(spec, mask)) continue;
823 : // register this spectrum and mask to matrices to solve
824 0 : specMat.column(nspec) = spec;
825 0 : maskMat.column(nspec) = mask;
826 0 : imgIdvec.push_back((uInt) i);
827 0 : ++nspec;
828 : } // end of image loop
829 0 : if (nspec < nshift_) {
830 0 : specMat.resize(nchan_, nspec, true);
831 0 : maskMat.resize(nchan_, nspec, true);
832 : }
833 0 : return nspec > 0;
834 0 : }
835 :
836 : } //# NAMESPACE CASA - END
|
