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