Line data Source code
1 : //# RegionManager.cc: framework independent class that provides
2 : //# functionality to tool of same name
3 : //# Copyright (C) 2007
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This program is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU General Public License as published by the Free
8 : //# Software Foundation; either version 2 of the License, or (at your option)
9 : //# any later version.
10 : //#
11 : //# This program is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14 : //# more details.
15 : //#
16 : //# You should have received a copy of the GNU General Public License along
17 : //# with this program; if not, write to the Free Software Foundation, Inc.,
18 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 :
27 : #include <imageanalysis/Regions/CasacRegionManager.h>
28 :
29 : #include <casacore/casa/Containers/Record.h>
30 : #include <casacore/casa/OS/File.h>
31 : #include <casacore/images/Images/TempImage.h>
32 : #include <casacore/images/Images/SubImage.h>
33 : #include <casacore/casa/Utilities/Regex.h>
34 :
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/images/Regions/WCBox.h>
37 : #include <casacore/lattices/LRegions/LCBox.h>
38 : #include <casacore/measures/Measures/Stokes.h>
39 : #include <casacore/tables/Tables/TableRecord.h>
40 :
41 : #include <imageanalysis/Annotations/AnnRegion.h>
42 : #include <imageanalysis/Annotations/RegionTextList.h>
43 : #include <imageanalysis/IO/ParameterParser.h>
44 : #include <imageanalysis/ImageAnalysis/ImageMetaData.h>
45 : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
46 :
47 : #include <casacore/lattices/LRegions/LCSlicer.h>
48 :
49 : #include <casacore/casa/namespace.h>
50 : #include <memory>
51 :
52 : using namespace casacore;
53 : namespace casa { //# name space casa begins
54 :
55 : const String CasacRegionManager::ALL = "ALL";
56 :
57 846 : CasacRegionManager::CasacRegionManager() : RegionManager() {}
58 :
59 29238 : CasacRegionManager::CasacRegionManager(
60 : const CoordinateSystem& csys
61 29238 : ) : RegionManager(csys) {}
62 :
63 30930 : CasacRegionManager::~CasacRegionManager() {}
64 :
65 22342 : vector<uInt> CasacRegionManager::_setPolarizationRanges(
66 : String& specification, const String& firstStokes, const uInt nStokes,
67 : const StokesControl stokesControl
68 : ) const {
69 44684 : LogOrigin origin("CasacRegionManager", __func__);
70 : //const LogIO *myLog = _getLog();
71 22342 : *_getLog() << origin;
72 :
73 22342 : vector<uInt> ranges(0);
74 22342 : CoordinateSystem csys = getcoordsys();
75 22342 : if (! csys.hasPolarizationCoordinate()) {
76 18190 : return ranges;
77 : }
78 4152 : specification.trim();
79 4152 : specification.ltrim('[');
80 4152 : specification.rtrim(']');
81 :
82 4152 : specification.upcase();
83 4152 : if (specification == ALL) {
84 0 : ranges.push_back(0);
85 0 : ranges.push_back(nStokes - 1);
86 0 : return ranges;
87 : }
88 4152 : if (specification.empty()) {
89 : // ranges.resize(2);
90 : // ranges[0] = 0;
91 3999 : ranges.push_back(0);
92 3999 : switch (stokesControl) {
93 85 : case USE_FIRST_STOKES:
94 85 : ranges.push_back(0);
95 85 : specification = firstStokes;
96 85 : break;
97 3914 : case USE_ALL_STOKES:
98 3914 : ranges.push_back(nStokes - 1);
99 3914 : specification = ALL;
100 3914 : break;
101 0 : default:
102 : // bug if we get here
103 0 : ThrowCc("Logic error, unhandled stokes control");
104 : break;
105 : };
106 3999 : return ranges;
107 : }
108 : // First split on commas and semi-colons.
109 : // in the past for polarization specification.
110 :
111 306 : Vector<String> parts = stringToVector(specification, std::regex("[,;]"));
112 306 : Vector<String> polNames = Stokes::allNames(false);
113 153 : uInt nNames = polNames.size();
114 306 : Vector<uInt> nameLengths(nNames);
115 5049 : for (uInt i=0; i<nNames; i++) {
116 4896 : nameLengths[i] = polNames[i].length();
117 : }
118 153 : uInt *lengthData = nameLengths.data();
119 :
120 306 : Vector<uInt> idx(nNames);
121 306 : Sort sorter;
122 153 : sorter.sortKey(lengthData, TpUInt, 0, Sort::Descending);
123 153 : sorter.sort(idx, nNames);
124 :
125 306 : Vector<String> sortedNames(nNames);
126 5049 : for (uInt i=0; i<nNames; i++) {
127 4896 : sortedNames[i] = polNames[idx[i]];
128 4896 : sortedNames[i].upcase();
129 : }
130 301 : for (uInt i=0; i<parts.size(); i++) {
131 153 : String part = parts[i];
132 153 : part.trim();
133 153 : Vector<String>::iterator iter = sortedNames.begin();
134 5092 : while (iter != sortedNames.end() && ! part.empty()) {
135 4939 : if (part.startsWith(*iter)) {
136 151 : Int stokesPix = csys.stokesPixelNumber(*iter);
137 151 : if (stokesPix >= int(nStokes)) {
138 0 : *_getLog() << "Polarization " << *iter << " specified in "
139 0 : << parts[i] << " does not exist in the specified "
140 : << "coordinate system for the specified number of "
141 0 : << "polarization parameters" << LogIO::EXCEPTION;
142 : }
143 151 : ranges.push_back(stokesPix);
144 151 : ranges.push_back(stokesPix);
145 : // consume the string
146 151 : part = part.substr(iter->length());
147 151 : if (! part.empty()) {
148 : // reset the iterator to start over at the beginning of the list for
149 : // the next specified polarization
150 3 : iter = sortedNames.begin();
151 : }
152 : }
153 : else {
154 4788 : iter++;
155 : }
156 : }
157 153 : if (! part.empty()) {
158 5 : *_getLog() << "(Sub)String " << part << " in stokes specification part " << parts[i]
159 5 : << " does not match a known polarization." << LogIO::EXCEPTION;
160 : }
161 158 : }
162 : uInt nSel;
163 148 : return ParameterParser::consolidateAndOrderRanges(nSel, ranges);
164 22352 : }
165 :
166 44647 : Bool CasacRegionManager::_supports2DBox(Bool except) const {
167 44647 : Bool ok = true;
168 44647 : const CoordinateSystem& csys = getcoordsys();
169 44647 : Vector<Int> axes;
170 44647 : if (csys.hasDirectionCoordinate()) {
171 44601 : axes = csys.directionAxesNumbers();
172 : }
173 46 : else if (csys.hasLinearCoordinate()) {
174 14 : axes = csys.linearAxesNumbers();
175 : }
176 : else {
177 32 : ok = false;
178 : }
179 44647 : if (ok) {
180 44615 : uInt nGood = 0;
181 133845 : for (uInt i=0; i<axes.size(); i++) {
182 89230 : if (axes[i] >= 0) {
183 89218 : nGood++;
184 : }
185 : }
186 44615 : if (nGood != 2) {
187 12 : ok = false;
188 : }
189 : }
190 44647 : if (except && ! ok) {
191 0 : *_getLog() << LogOrigin("CasacRegionManager", __func__)
192 0 : << "This image does not have a 2-D direction or linear coordinate";
193 : }
194 44647 : return ok;
195 44647 : }
196 :
197 260 : vector<Double> CasacRegionManager::_setBoxCorners(const String& box) const {
198 260 : if (! box.empty()) {
199 260 : _supports2DBox(true);
200 : }
201 260 : Vector<String> boxParts = stringToVector(box);
202 260 : AlwaysAssert(boxParts.size() % 4 == 0, AipsError);
203 260 : vector<Double> corners(boxParts.size());
204 537 : for(uInt i=0; i<boxParts.size()/4; i++) {
205 277 : uInt index = 4*i;
206 1385 : for (uInt j=0; j<4; j++) {
207 1108 : boxParts[index + j].trim();
208 1108 : if (! boxParts[index + j].matches(RXdouble)) {
209 0 : *_getLog() << "Box spec contains non numeric characters and so is invalid"
210 0 : << LogIO::EXCEPTION;
211 : }
212 1108 : corners[index + j] = String::toDouble(boxParts[index + j]);
213 : }
214 : }
215 520 : return corners;
216 260 : }
217 :
218 29716 : Record CasacRegionManager::fromBCS(
219 : String& diagnostics, uInt& nSelectedChannels, String& stokes,
220 : const Record * const ®ionPtr, const String& regionName,
221 : const String& chans, const StokesControl stokesControl,
222 : const String& box, const IPosition& imShape, const String& imageName,
223 : Bool verbose
224 : ) {
225 59432 : LogOrigin origin("CasacRegionManager", __func__);
226 29716 : Record regionRecord;
227 29716 : if (! box.empty()) {
228 266 : ThrowIf(
229 : regionPtr != nullptr,
230 : "box and regionPtr cannot be simultaneously specified"
231 : );
232 267 : ThrowIf(
233 : box.freq(",") % 4 != 3,
234 : "box not specified correctly"
235 : );
236 265 : if (regionName.empty()) {
237 263 : regionRecord = fromBCS(
238 : diagnostics, nSelectedChannels, stokes,
239 : chans, stokesControl, box, imShape
240 231 : ).toRecord("");
241 231 : if (verbose) {
242 228 : *_getLog() << origin;
243 228 : *_getLog() << LogIO::NORMAL << "Using specified box(es) "
244 228 : << box << LogIO::POST;
245 : }
246 : }
247 : else {
248 2 : auto corners = stringToVector(box, ',');
249 2 : auto iter = corners.begin();
250 2 : auto end = corners.end();
251 2 : String crtfBoxString = "box[[";
252 2 : auto count = 0;
253 10 : for (; iter != end; ++iter, ++count) {
254 8 : iter->trim();
255 8 : if (count > 0 && count % 4 == 0) {
256 0 : crtfBoxString += "\nbox[[";
257 : }
258 8 : crtfBoxString += *iter + "pix";
259 8 : if (count % 2 == 0) {
260 4 : crtfBoxString += ",";
261 : }
262 : else {
263 : // close pixel pair
264 4 : crtfBoxString += "]";
265 4 : if (count - 1 % 4 == 0) {
266 : // first pixel pair done, start second pixel pair
267 2 : crtfBoxString += ",[";
268 : }
269 : else {
270 : // second pixel pair done, close box
271 2 : crtfBoxString += "]";
272 : }
273 : }
274 : }
275 2 : _setRegion(
276 : regionRecord, diagnostics, regionName, imShape,
277 : imageName, crtfBoxString, chans, stokes, verbose
278 : );
279 2 : }
280 : }
281 29450 : else if (regionPtr) {
282 7336 : ThrowIf(
283 : ! (regionName.empty() && chans.empty() && stokes.empty()),
284 : "regionPtr and regionName, chans, and/or stokes cannot "
285 : "be simultaneously specified"
286 : );
287 7336 : _setRegion(regionRecord, diagnostics, regionPtr);
288 7336 : stokes = _stokesFromRecord(regionRecord, stokesControl, imShape);
289 : }
290 22114 : else if (! regionName.empty()) {
291 33 : _setRegion(
292 : regionRecord, diagnostics, regionName,
293 : imShape, imageName, "", chans, stokes, verbose
294 : );
295 11 : if (verbose) {
296 11 : *_getLog() << origin;
297 11 : *_getLog() << LogIO::NORMAL << diagnostics << LogIO::POST;
298 : }
299 11 : stokes = _stokesFromRecord(regionRecord, stokesControl, imShape);
300 : }
301 : else {
302 22092 : vector<uInt> chanEndPts, polEndPts;
303 22092 : regionRecord = fromBCS(
304 : diagnostics, nSelectedChannels, stokes,
305 : chans, stokesControl, box, imShape
306 22069 : ).toRecord("");
307 22069 : if (verbose) {
308 1485 : *_getLog() << origin;
309 1485 : *_getLog() << LogIO::NORMAL << "No directional region specified. Using full positional plane."
310 1485 : << LogIO::POST;
311 : }
312 22069 : const CoordinateSystem& csys = getcoordsys();
313 22069 : if (csys.hasSpectralAxis()) {
314 4126 : if (verbose) {
315 1349 : if (chans.empty()) {
316 1096 : *_getLog() << LogIO::NORMAL << "Using all spectral channels."
317 1096 : << LogIO::POST;
318 : }
319 : else {
320 253 : *_getLog() << LogIO::NORMAL << "Using channel range(s) "
321 253 : << _pairsToString(chanEndPts) << LogIO::POST;
322 : }
323 : }
324 : }
325 22069 : if (csys.hasPolarizationCoordinate() && verbose) {
326 1222 : if (stokes.empty()) {
327 0 : switch (stokesControl) {
328 0 : case USE_ALL_STOKES:
329 0 : *_getLog() << LogIO::NORMAL << "Using all polarizations " << LogIO::POST;
330 0 : break;
331 0 : case USE_FIRST_STOKES:
332 0 : *_getLog() << LogIO::NORMAL << "polarization "
333 0 : << csys.stokesAtPixel(0) << LogIO::POST;
334 0 : break;
335 0 : default:
336 0 : break;
337 : }
338 : }
339 : else {
340 1222 : *_getLog() << LogIO::NORMAL << "Using polarizations " << stokes << LogIO::POST;
341 : }
342 : }
343 22115 : }
344 59298 : return regionRecord;
345 29783 : }
346 :
347 7 : Record CasacRegionManager::regionFromString(
348 : const CoordinateSystem& csys, const String& regionStr,
349 : const String& imageName, const IPosition& imShape,
350 : Bool verbose
351 : ) {
352 7 : CasacRegionManager mgr(csys);
353 7 : Record reg;
354 7 : String diag;
355 10 : mgr._setRegion(
356 : reg, diag, regionStr, imShape, imageName, "", "", "", verbose
357 : );
358 12 : return reg;
359 9 : }
360 :
361 7336 : void CasacRegionManager::_setRegion(
362 : Record& regionRecord, String& diagnostics,
363 : const Record* regionPtr
364 : ) {
365 : // region record pointer provided
366 7336 : regionRecord = *(regionPtr);
367 : // set stokes from the region record
368 7336 : diagnostics = "used provided region record";
369 7336 : }
370 :
371 31 : void CasacRegionManager::_setRegion(
372 : Record& regionRecord, String& diagnostics,
373 : const String& regionName, const IPosition& imShape,
374 : const String& imageName,
375 : const String& prependBox,
376 : const String& globalOverrideChans,
377 : const String& globalStokesOverride,
378 : Bool verbose
379 : ) {
380 31 : if (regionName.empty() && imageName.empty()) {
381 0 : regionRecord = Record();
382 0 : diagnostics = "No region string";
383 4 : return;
384 : }
385 : // region name provided
386 31 : const static Regex image("(.+):(.+)");
387 : // OSX no longer seems to like this one
388 : // const static Regex regionText(
389 : // "^[[:space:]]*[[:alpha:]]+[[:space:]]*\\[(.*)+,(.*)+\\]"
390 : // );
391 : const static Regex regionText(
392 : "^[[:space:]]*[[:alpha:]]+[[:space:]]*\\["
393 31 : );
394 31 : File myFile(regionName);
395 31 : const CoordinateSystem csys = getcoordsys();
396 31 : if (myFile.exists()) {
397 10 : ThrowIf(
398 : ! myFile.isReadable(),
399 : "File " + regionName + " exists but is not readable."
400 : );
401 10 : std::unique_ptr<Record> rec;
402 : try {
403 22 : rec.reset(readImageFile(regionName, ""));
404 : }
405 6 : catch(const AipsError& x) {
406 6 : }
407 10 : if (rec) {
408 4 : ThrowIf(
409 : ! globalOverrideChans.empty() || ! globalStokesOverride.empty()
410 : || ! prependBox.empty(),
411 : "a binary region file and any of box, chans and/or stokes cannot "
412 : "be specified simultaneously"
413 : );
414 4 : regionRecord = *rec;
415 4 : diagnostics = "Region read from binary region file " + regionName;
416 4 : return;
417 : }
418 : try {
419 : // CRTF file attempt
420 : RegionTextList annList(
421 : regionName, csys, imShape, prependBox,
422 : globalOverrideChans, globalStokesOverride
423 11 : );
424 1 : regionRecord = annList.regionAsRecord();
425 1 : diagnostics = "Region read from CRTF file " + regionName;
426 1 : }
427 5 : catch (const AipsError& x) {
428 5 : ThrowCc(
429 : regionName + " is neither a valid binary region file, "
430 : "nor a valid region text file."
431 : );
432 5 : }
433 10 : }
434 21 : else if (regionName.contains(regionText)) {
435 : // region spec is raw CASA region plaintext
436 : try {
437 : RegionTextList annList(
438 : csys, regionName, imShape, prependBox, globalOverrideChans,
439 : globalStokesOverride, verbose
440 13 : );
441 11 : regionRecord = annList.regionAsRecord();
442 11 : diagnostics = "Region read from text string " + regionName;
443 11 : }
444 1 : catch (const AipsError& x) {
445 1 : ThrowCc(x.getMesg());
446 1 : }
447 : }
448 9 : else if (regionName.matches(image) || ! imageName.empty()) {
449 4 : ImageRegion imRegion;
450 4 : String imagename, region;
451 4 : if (regionName.matches(image)) {
452 6 : String res[2];
453 1 : casacore::split(regionName, res, 2, ":");
454 1 : imagename = res[0];
455 1 : region = res[1];
456 3 : }
457 : else {
458 : // imageName is not empty if we get here
459 3 : imagename = imageName;
460 3 : region = regionName;
461 : }
462 : try {
463 4 : Record *myRec = tableToRecord(imagename, region);
464 4 : ThrowIf(
465 : ! globalOverrideChans.empty() || ! globalStokesOverride.empty()
466 : || ! prependBox.empty(),
467 : "a region-in-image and any of box, chans and/or stokes cannot "
468 : "be specified simultaneously"
469 : );
470 4 : if (Table::isReadable(imagename)) {
471 3 : ThrowIf(
472 : myRec == 0,
473 : "Region " + region + " not found in image "
474 : + imagename
475 : );
476 3 : regionRecord = *myRec;
477 6 : diagnostics = "Used region " + region + " from image "
478 3 : + imagename + " table description";
479 : }
480 : else {
481 1 : *_getLog() << "Cannot read image " << imagename
482 1 : << " to get region " << region << LogIO::EXCEPTION;
483 : }
484 : }
485 1 : catch (const AipsError&) {
486 1 : ThrowCc(
487 : "Unable to open region file or region table description "
488 : + region + " in image " + imagename
489 : );
490 1 : }
491 6 : }
492 : else {
493 5 : ostringstream oss;
494 5 : oss << "Unable to open region file or region table description "
495 5 : << regionName << "." << endl
496 5 : << "If it is supposed to be a text string its format is incorrect";
497 5 : ThrowCc(oss.str());
498 5 : }
499 47 : }
500 :
501 22358 : ImageRegion CasacRegionManager::fromBCS(
502 : String& diagnostics, uInt& nSelectedChannels, String& stokes,
503 : const String& chans,
504 : const StokesControl stokesControl, const String& box,
505 : const IPosition& imShape
506 : ) const {
507 22358 : const CoordinateSystem& csys = getcoordsys();
508 : vector<uInt> chanEndPts = setSpectralRanges(
509 : chans, nSelectedChannels, imShape
510 22374 : );
511 22342 : Int polAxisNumber = csys.polarizationAxisNumber();
512 22342 : uInt nTotalPolarizations = polAxisNumber >= 0 ? imShape[polAxisNumber] : 0;
513 22342 : String firstStokes = polAxisNumber >= 0 ? csys.stokesAtPixel(0) : "";
514 : vector<uInt> polEndPts = _setPolarizationRanges(
515 : stokes, firstStokes,
516 : nTotalPolarizations, stokesControl
517 22342 : );
518 22337 : vector<Double> boxCorners;
519 22337 : if (box.empty()) {
520 22077 : if (_supports2DBox(false)) {
521 22055 : if (
522 22055 : csys.hasDirectionCoordinate()
523 22055 : || csys.hasLinearCoordinate()
524 : ) {
525 22055 : Vector<Int> dirAxesNumbers;
526 22055 : if (csys.hasDirectionCoordinate()) {
527 22054 : dirAxesNumbers = csys.directionAxesNumbers();
528 : }
529 : else {
530 1 : dirAxesNumbers = csys.linearAxesNumbers();
531 : }
532 22055 : Vector<Int> dirShape(2);
533 22055 : dirShape[0] = imShape[dirAxesNumbers[0]];
534 22055 : dirShape[1] = imShape[dirAxesNumbers[1]];
535 22055 : boxCorners.resize(4);
536 22055 : boxCorners[0] = 0;
537 22055 : boxCorners[1] = 0;
538 22055 : boxCorners[2] = dirShape[0] - 1;
539 22055 : boxCorners[3] = dirShape[1] - 1;
540 22055 : }
541 : }
542 : }
543 : else {
544 260 : boxCorners = _setBoxCorners(box);
545 : }
546 : return _fromBCS(
547 : diagnostics, boxCorners,
548 : chanEndPts, polEndPts, imShape
549 44674 : );
550 22449 : }
551 :
552 22337 : ImageRegion CasacRegionManager::_fromBCS(
553 : String& diagnostics, const vector<Double>& boxCorners,
554 : const vector<uInt>& chanEndPts, const vector<uInt>& polEndPts,
555 : const IPosition imShape
556 : ) const {
557 44674 : LogOrigin origin("CasacRegionManager", __func__);
558 22337 : *_getLog() << origin;
559 22337 : Vector<Double> blc(imShape.nelements(), 0);
560 22337 : Vector<Double> trc(imShape.nelements(), 0);
561 22337 : const CoordinateSystem csys = getcoordsys();
562 22337 : Vector<Int> directionAxisNumbers = csys.directionAxesNumbers();
563 22337 : vector<Int> linearAxisNumbers = csys.linearAxesNumbers().tovector();
564 : // Stupidly, sometimes the values returned by linearAxesNumbers can be less than 0
565 : // This needs to be fixed in the implementation of that method
566 22337 : vector<Int>::iterator iter = linearAxisNumbers.begin();
567 22337 : vector<Int>::iterator end = linearAxisNumbers.end();
568 22373 : while(iter != end) {
569 36 : if (*iter < 0) {
570 6 : iter = linearAxisNumbers.erase(iter);
571 : }
572 36 : ++iter;
573 : }
574 22337 : Int spectralAxisNumber = csys.spectralAxisNumber();
575 22337 : Int polarizationAxisNumber = csys.polarizationAxisNumber();
576 :
577 22337 : Vector<Double> xCorners(boxCorners.size()/2);
578 22337 : Vector<Double> yCorners(xCorners.size());
579 66958 : for (uInt i=0; i<xCorners.size(); i++) {
580 44648 : Double x = boxCorners[2*i];
581 44648 : Double y = boxCorners[2*i + 1];
582 :
583 44648 : if (x < 0 || y < 0 ) {
584 15 : *_getLog() << "blc in box spec is less than 0" << LogIO::EXCEPTION;
585 : }
586 44633 : if (csys.hasDirectionCoordinate()) {
587 44631 : if (
588 44631 : x >= imShape[directionAxisNumbers[0]]
589 44631 : || y >= imShape[directionAxisNumbers[1]]
590 : ) {
591 12 : *_getLog() << "dAxisNum0=" <<directionAxisNumbers[0] <<" dAxisNum1="<<directionAxisNumbers[1];
592 12 : *_getLog() << "x="<<x<<" imShape[0]="<<imShape[directionAxisNumbers[0]]<< " y="<<y<<" imShape[1]="<<imShape[directionAxisNumbers[1]]<<LogIO::POST;
593 12 : *_getLog() << "trc in box spec is greater than or equal to number "
594 12 : << "of direction coordinate pixels in the image" << LogIO::EXCEPTION;
595 : }
596 : }
597 2 : else if (
598 2 : csys.hasLinearCoordinate()
599 4 : && (
600 2 : x >= imShape[linearAxisNumbers[0]]
601 2 : || y >= imShape[linearAxisNumbers[1]]
602 : )
603 : ) {
604 0 : *_getLog() << "trc in box spec is greater than or equal to number "
605 0 : << "of linear coordinate pixels in the image" << LogIO::EXCEPTION;
606 : }
607 44621 : xCorners[i] = x;
608 44621 : yCorners[i] = y;
609 : }
610 22310 : Vector<Double> polEndPtsDouble(polEndPts.size());
611 30552 : for (uInt i=0; i<polEndPts.size(); ++i) {
612 8242 : polEndPtsDouble[i] = (Double)polEndPts[i];
613 : }
614 :
615 22310 : Bool csysSupports2DBox = _supports2DBox(false);
616 22310 : uInt nRegions = 1;
617 22310 : if (csysSupports2DBox) {
618 22288 : if (csys.hasDirectionCoordinate()) {
619 22287 : nRegions *= boxCorners.size()/4;
620 : }
621 22288 : if (csys.hasLinearCoordinate()) {
622 23 : nRegions *= boxCorners.size()/4;
623 : }
624 : }
625 22310 : if (csys.hasPolarizationCoordinate()) {
626 4120 : nRegions *= polEndPts.size()/2;
627 : }
628 22310 : if (csys.hasSpectralAxis()) {
629 4239 : nRegions *= chanEndPts.size()/2;
630 : }
631 22310 : Vector<Double> extXCorners(2*nRegions, 0);
632 22310 : Vector<Double> extYCorners(2*nRegions, 0);
633 22310 : Vector<Double> extPolEndPts(2*nRegions, 0);
634 22310 : Vector<Double> extChanEndPts(2*nRegions, 0);
635 :
636 22310 : uInt count = 0;
637 :
638 22310 : if (csysSupports2DBox) {
639 44593 : for (uInt i=0; i<max(uInt(1), xCorners.size()/2); i++) {
640 44611 : for (uInt j=0; j<max((uInt)1, polEndPts.size()/2); j++) {
641 44618 : for (uInt k=0; k<max(uInt(1), chanEndPts.size()/2); k++) {
642 22312 : if (
643 22312 : csys.hasDirectionCoordinate()
644 22312 : || csys.hasLinearCoordinate()
645 : ) {
646 22312 : extXCorners[2*count] = xCorners[2*i];
647 22312 : extXCorners[2*count + 1] = xCorners[2*i + 1];
648 22312 : extYCorners[2*count] = yCorners[2*i];
649 22312 : extYCorners[2*count + 1] = yCorners[2*i + 1];
650 : }
651 22312 : if (csys.hasPolarizationCoordinate()) {
652 4133 : extPolEndPts[2*count] = polEndPtsDouble[2*j];
653 4133 : extPolEndPts[2*count + 1] = polEndPtsDouble[2*j + 1];
654 :
655 : }
656 22312 : if (csys.hasSpectralAxis()) {
657 4239 : extChanEndPts[2*count] = chanEndPts[2*k];
658 4239 : extChanEndPts[2*count + 1] = chanEndPts[2*k + 1];
659 : }
660 22312 : count++;
661 : }
662 : }
663 : }
664 : }
665 : else {
666 : // here we have neither a direction nor linear coordinate with two
667 : // pixel axes which are greater than 0
668 44 : for (uInt j=0; j<max((uInt)1, polEndPts.size()/2); j++) {
669 44 : for (uInt k=0; k<max(uInt(1), chanEndPts.size()/2); k++) {
670 22 : if (csys.hasPolarizationCoordinate()) {
671 4 : extPolEndPts[2*count] = polEndPtsDouble[2*j];
672 4 : extPolEndPts[2*count + 1] = polEndPtsDouble[2*j + 1];
673 : }
674 22 : if (csys.hasSpectralAxis()) {
675 22 : extChanEndPts[2*count] = chanEndPts[2*k];
676 22 : extChanEndPts[2*count + 1] = chanEndPts[2*k + 1];
677 : }
678 22 : count++;
679 : }
680 : }
681 : }
682 22310 : map<uInt, Vector<Double> > axisCornerMap;
683 44644 : for (uInt i=0; i<nRegions; i++) {
684 75384 : for (uInt axisNumber=0; axisNumber<csys.nPixelAxes(); axisNumber++) {
685 53050 : if (
686 : (
687 53050 : directionAxisNumbers.size() > 1
688 53016 : && (Int)axisNumber == directionAxisNumbers[0]
689 : )
690 106100 : || (
691 30739 : ! csys.hasDirectionCoordinate()
692 34 : && linearAxisNumbers.size() > 1
693 2 : && (Int)axisNumber == linearAxisNumbers[0]
694 : )
695 : ) {
696 22312 : axisCornerMap[axisNumber] = extXCorners;
697 : }
698 30738 : else if (
699 : (
700 30738 : directionAxisNumbers.size() > 1
701 30705 : && (Int)axisNumber == directionAxisNumbers[1]
702 : )
703 61476 : || (
704 8427 : ! csys.hasDirectionCoordinate()
705 33 : && linearAxisNumbers.size() > 1
706 1 : && (Int)axisNumber == linearAxisNumbers[1]
707 : )
708 : ) {
709 22312 : axisCornerMap[axisNumber] = extYCorners;
710 : }
711 8426 : else if ((Int)axisNumber == spectralAxisNumber) {
712 4261 : axisCornerMap[axisNumber] = extChanEndPts;
713 : }
714 4165 : else if ((Int)axisNumber == polarizationAxisNumber) {
715 4137 : axisCornerMap[axisNumber] = extPolEndPts;
716 : }
717 : else {
718 28 : Vector<Double> range(2, 0);
719 28 : range[1] = imShape[axisNumber] - 1;
720 28 : axisCornerMap[axisNumber] = range;
721 28 : }
722 : }
723 : }
724 22310 : ImageRegion imRegion;
725 44637 : for (uInt i=0; i<nRegions; i++) {
726 75384 : for (uInt axisNumber=0; axisNumber<csys.nPixelAxes(); axisNumber++) {
727 53050 : blc(axisNumber) = axisCornerMap[axisNumber][2*i];
728 53050 : trc(axisNumber) = axisCornerMap[axisNumber][2*i + 1];
729 : }
730 22334 : LCBox lcBox(blc, trc, imShape);
731 22327 : WCBox wcBox(lcBox, csys);
732 22327 : ImageRegion thisRegion(wcBox);
733 : imRegion = (i == 0)
734 : ? thisRegion
735 22327 : : imRegion = *(doUnion(imRegion, thisRegion));
736 22327 : }
737 22303 : ostringstream os;
738 22303 : os << "Used image region from " << endl;
739 22303 : if (csys.hasDirectionCoordinate()) {
740 22280 : os << " position box corners: ";
741 44577 : for (uInt i=0; i<boxCorners.size()/4; i++) {
742 22297 : os << boxCorners[4*i] << ", " << boxCorners[4*i + 1]
743 22297 : << ", " << boxCorners[4*i + 2] << ", " << boxCorners[4*i + 3];
744 22297 : if (i < boxCorners.size()/4 - 1) {
745 17 : os << "; ";
746 : }
747 : }
748 : }
749 22303 : if (getcoordsys().hasSpectralAxis()) {
750 4232 : os << " spectral channel ranges: " << _pairsToString(chanEndPts);
751 : }
752 22303 : if (getcoordsys().hasPolarizationCoordinate()) {
753 4113 : os << " polarization pixel ranges: " << _pairsToString(polEndPts);
754 : }
755 22303 : diagnostics = os.str();
756 44606 : return imRegion;
757 22624 : }
758 :
759 8598 : String CasacRegionManager::_pairsToString(const vector<uInt>& pairs) const {
760 8598 : ostringstream os;
761 8598 : uInt nPairs = pairs.size()/2;
762 16947 : for (uInt i=0; i<nPairs; i++) {
763 8349 : os << pairs[2*i] << " - " << pairs[2*i + 1];
764 8349 : if (i < nPairs - 1) {
765 4 : os << "; ";
766 : }
767 : }
768 25794 : return os.str();
769 8598 : }
770 :
771 7347 : String CasacRegionManager::_stokesFromRecord(
772 : const Record& region, const StokesControl stokesControl, const IPosition& shape
773 : ) const {
774 : // FIXME This implementation is incorrect for complex, recursive records
775 7347 : String stokes = "";
776 7347 : CoordinateSystem csys = getcoordsys();
777 7347 : if(! csys.hasPolarizationCoordinate()) {
778 1143 : return stokes;
779 : }
780 6204 : Int polAxis = csys.polarizationAxisNumber();
781 6204 : if (shape[polAxis] == 1) {
782 : // degenerate stokes axis
783 6040 : return csys.stokesAtPixel(0);
784 : }
785 164 : uInt stokesBegin = 0;
786 164 : uInt stokesEnd = 0;
787 164 : if (region.empty()) {
788 24 : stokesEnd = stokesControl == USE_FIRST_STOKES
789 : ? 0
790 24 : : csys.stokesCoordinate().stokes().size() - 1;
791 : }
792 : else {
793 280 : std::unique_ptr<ImageRegion> imreg(ImageRegion::fromRecord(region, ""));
794 140 : Array<Float> blc, trc;
795 140 : Bool oneRelAccountedFor = false;
796 140 : if (imreg->isLCSlicer()) {
797 74 : blc = imreg->asLCSlicer().blc();
798 74 : if ((Int)blc.size() <= polAxis) {
799 0 : *_getLog() << LogIO::WARN << "Cannot determine stokes. "
800 : << "blc of input region does not include the polarization coordinate."
801 0 : << LogIO::POST;
802 0 : return stokes;
803 : }
804 74 : trc = imreg->asLCSlicer().trc();
805 74 : if ((Int)trc.size() <= polAxis) {
806 0 : *_getLog() << LogIO::WARN << "Cannot determine stokes. "
807 : << "trc of input region does not include the polarization coordinate."
808 0 : << LogIO::POST;
809 0 : return stokes;
810 : }
811 74 : stokesBegin = (uInt)((Vector<Float>)blc)[polAxis];
812 74 : stokesEnd = (uInt)((Vector<Float>)trc)[polAxis];
813 74 : oneRelAccountedFor = true;
814 : }
815 66 : else if (
816 132 : RegionManager::isPixelRegion(
817 66 : *(ImageRegion::fromRecord(region, ""))
818 : )
819 : ) {
820 0 : region.toArray("blc", blc);
821 0 : region.toArray("trc", trc);
822 0 : stokesBegin = (uInt)((Vector<Float>)blc)[polAxis];
823 0 : stokesEnd = (uInt)((Vector<Float>)trc)[polAxis];
824 : }
825 66 : else if (region.fieldNumber("x") >= 0 && region.fieldNumber("y") >= 0) {
826 : // world polygon
827 1 : oneRelAccountedFor = true;
828 1 : stokesBegin = 0;
829 2 : stokesEnd = stokesControl == USE_FIRST_STOKES
830 1 : ? 0 : shape[polAxis];
831 : }
832 65 : else if (region.fieldNumber("blc") >= 0 && region.fieldNumber("blc") >= 0) {
833 : // world box
834 59 : Record blcRec = region.asRecord("blc");
835 59 : Record trcRec = region.asRecord("trc");
836 59 : String polField = "*" + String::toString(polAxis + 1);
837 59 : stokesBegin = blcRec.isDefined(polField)
838 177 : ? (Int)blcRec.asRecord(
839 118 : String(polField)
840 59 : ).asDouble("value")
841 : : 0;
842 59 : stokesEnd = trcRec.isDefined(polField)
843 118 : ? (Int)blcRec.asRecord(
844 118 : String(polField)
845 59 : ).asDouble("value")
846 : : stokesControl == USE_FIRST_STOKES
847 : ? 0
848 0 : : shape[polAxis];
849 59 : if (! blcRec.isDefined(polField)) {
850 0 : oneRelAccountedFor = true;
851 : }
852 59 : }
853 : else {
854 : // FIXME not very nice, but until all can be implemented this will have to do
855 6 : *_getLog() << LogIO::WARN << "Stokes cannot be determined because this "
856 : << "region type is not handled yet. But chances are very good this is no need to be alarmed."
857 6 : << LogIO::POST;
858 6 : return stokes;
859 : }
860 134 : if (
861 134 : ! oneRelAccountedFor
862 193 : && region.isDefined("oneRel")
863 193 : && region.asBool("oneRel")
864 : ) {
865 59 : stokesBegin--;
866 59 : stokesEnd--;
867 : }
868 152 : }
869 488 : for (uInt i=stokesBegin; i<=stokesEnd; i++) {
870 330 : stokes += csys.stokesAtPixel(i);
871 : }
872 158 : return stokes;
873 7347 : }
874 :
875 22214 : vector<uInt> CasacRegionManager::_initSpectralRanges(
876 : uInt& nSelectedChannels, const IPosition& imShape
877 : ) const {
878 22214 : vector<uInt> ranges(0);
879 22214 : if (! getcoordsys().hasSpectralAxis()) {
880 18248 : nSelectedChannels = 0;
881 18248 : return ranges;
882 : }
883 3966 : uInt specNum = getcoordsys().spectralAxisNumber();
884 3966 : ThrowIf(
885 : specNum >= imShape.size(),
886 : "Spectral axis number " + String::toString(specNum)
887 : + " must be less than number of shape dimensions "
888 : + String::toString(imShape.size())
889 : );
890 3966 : uInt nChannels = imShape[getcoordsys().spectralAxisNumber()];
891 3966 : ranges.push_back(0);
892 3966 : ranges.push_back(nChannels - 1);
893 3966 : nSelectedChannels = nChannels;
894 3966 : return ranges;
895 0 : }
896 :
897 326 : vector<uInt> CasacRegionManager::setSpectralRanges(
898 : uInt& nSelectedChannels,
899 : const Record *const regionRec, const IPosition& imShape
900 : ) const {
901 326 : if (regionRec == 0 || ! getcoordsys().hasSpectralAxis()) {
902 177 : return _initSpectralRanges(nSelectedChannels, imShape);
903 : }
904 : else {
905 149 : return _spectralRangeFromRegionRecord(nSelectedChannels, regionRec, imShape);
906 : }
907 : }
908 :
909 22388 : vector<uInt> CasacRegionManager::setSpectralRanges(
910 : String specification, uInt& nSelectedChannels,
911 : /*const String& globalChannelOverride,
912 : const String& globalStokesOverride, */ const IPosition& imShape
913 : ) const {
914 44776 : LogOrigin origin("CasacRegionManager", __func__);
915 22388 : *_getLog() << origin;
916 22388 : specification.trim();
917 22388 : String x = specification;
918 22388 : x.upcase();
919 22388 : if (x.empty() || x == ALL) {
920 22037 : return _initSpectralRanges(nSelectedChannels, imShape);
921 : }
922 351 : else if (! getcoordsys().hasSpectralAxis()) {
923 0 : *_getLog() << LogIO::WARN << "Channel specification is "
924 : << "not empty but the coordinate system has no spectral axis."
925 0 : << "Channel specification will be ignored" << LogIO::POST;
926 0 : nSelectedChannels = 0;
927 0 : return vector<uInt>(0);
928 : }
929 351 : else if (specification.contains("range")) {
930 : // this is a specification in the "new" ASCII region format
931 1 : return _spectralRangeFromRangeFormat(nSelectedChannels, specification, imShape);
932 : }
933 : else {
934 350 : uInt nChannels = imShape[getcoordsys().spectralAxisNumber()];
935 : return ParameterParser::spectralRangesFromChans(
936 : nSelectedChannels, specification, nChannels
937 350 : );
938 : }
939 22409 : }
940 :
941 149 : vector<uInt> CasacRegionManager::_spectralRangeFromRegionRecord(
942 : uInt& nSelectedChannels, const Record *const regionRec,
943 : const IPosition& imShape
944 : ) const {
945 149 : const CoordinateSystem& csys = getcoordsys();
946 149 : TempImage<Float> x(imShape, csys);
947 149 : x.set(0);
948 298 : SPCIIF subimage = SubImageFactory<Float>::createSubImageRO(
949 298 : x, *regionRec, "", _getLog(), AxesSpecifier(), false, true
950 149 : );
951 149 : uInt nChan = 0;
952 : {
953 149 : ImageMetaData<Float> md(subimage);
954 149 : nChan = md.nChannels();
955 149 : }
956 149 : const SpectralCoordinate& subsp = subimage->coordinates().spectralCoordinate();
957 : Double subworld;
958 149 : subsp.toWorld(subworld, 0);
959 149 : const SpectralCoordinate& imsp = csys.spectralCoordinate();
960 : Double pixOff;
961 149 : imsp.toPixel(pixOff, subworld);
962 149 : Int specAxisNumber = csys.spectralAxisNumber();
963 149 : IPosition start(subimage->ndim(), 0);
964 149 : ArrayLattice<Bool> mask(subimage->getMask());
965 149 : IPosition planeShape = subimage->shape();
966 149 : vector<uInt> myList;
967 149 : planeShape[specAxisNumber] = 1;
968 313 : for (uInt i=0; i<nChan; i++) {
969 164 : start[specAxisNumber] = i;
970 164 : Array<Bool> maskSlice;
971 164 : mask.getSlice(maskSlice, start, planeShape);
972 164 : if (anyTrue(maskSlice)) {
973 164 : uInt real = i + (uInt)rint(pixOff);
974 164 : myList.push_back(real);
975 164 : myList.push_back(real);
976 : }
977 164 : }
978 298 : return ParameterParser::consolidateAndOrderRanges(nSelectedChannels, myList);
979 149 : }
980 :
981 1 : vector<uInt> CasacRegionManager::_spectralRangeFromRangeFormat(
982 : uInt& nSelectedChannels, const String& specification,
983 : const IPosition& imShape /*, const String& globalChannelOverride,
984 : const String& globalStokesOverride */
985 : ) const {
986 : Bool spectralParmsUpdated;
987 : // check and make sure there are no disallowed parameters
988 1 : const CoordinateSystem csys = getcoordsys();
989 : RegionTextParser::ParamSet parms = RegionTextParser::getParamSet(
990 1 : spectralParmsUpdated, *_getLog(),
991 2 : specification, "", csys, std::shared_ptr<std::pair<MFrequency, MFrequency> >(nullptr),
992 1 : std::shared_ptr<Vector<Stokes::StokesTypes> >(nullptr)
993 3 : );
994 1 : RegionTextParser::ParamSet::const_iterator end = parms.end();
995 1 : for (
996 1 : RegionTextParser::ParamSet::const_iterator iter=parms.begin();
997 2 : iter!=end; iter++
998 : ) {
999 1 : AnnotationBase::Keyword key = iter->first;
1000 1 : ThrowIf(
1001 : key != AnnotationBase::FRAME && key != AnnotationBase::RANGE
1002 : && key != AnnotationBase::VELTYPE && key != AnnotationBase::RESTFREQ,
1003 : "Non-frequency related keyword '"
1004 : + AnnotationBase::keywordToString(key)
1005 : + "' found."
1006 : );
1007 : }
1008 : // Parameters OK. We need to modify the input string so we can construct an AnnRegion
1009 : // from which to get the spectral range information
1010 1 : String regSpec = "box[[0pix, 0pix], [1pix, 1pix]] " + specification;
1011 2 : RegionTextParser parser(csys, imShape, regSpec);
1012 1 : vector<uInt> range(2);
1013 1 : ThrowIf(
1014 : parser.getLines().empty(),
1015 : "The specified spectral range " + specification
1016 : + " does not intersect the image spectral range."
1017 : );
1018 1 : auto ann = parser.getLines()[0].getAnnotationBase();
1019 1 : /*const AnnRegion*/ auto *reg = dynamic_cast<const AnnRegion*>(ann.get());
1020 1 : ThrowIf(! reg, "Dynamic cast failed");
1021 1 : /*vector<Double>*/ const auto drange = reg->getSpectralPixelRange();
1022 1 : range[0] = uInt(max(0.0, floor(drange[0] + 0.5)));
1023 1 : range[1] = uInt(floor(drange[1] + 0.5));
1024 1 : nSelectedChannels = range[1] - range[0] + 1;
1025 2 : return range;
1026 1 : }
1027 :
1028 : }
|
