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 0 : CasacRegionManager::CasacRegionManager() : RegionManager() {}
58 :
59 0 : CasacRegionManager::CasacRegionManager(
60 : const CoordinateSystem& csys
61 0 : ) : RegionManager(csys) {}
62 :
63 0 : CasacRegionManager::~CasacRegionManager() {}
64 :
65 0 : vector<uInt> CasacRegionManager::_setPolarizationRanges(
66 : String& specification, const String& firstStokes, const uInt nStokes,
67 : const StokesControl stokesControl
68 : ) const {
69 0 : LogOrigin origin("CasacRegionManager", __func__);
70 : //const LogIO *myLog = _getLog();
71 0 : *_getLog() << origin;
72 :
73 0 : vector<uInt> ranges(0);
74 0 : CoordinateSystem csys = getcoordsys();
75 0 : if (! csys.hasPolarizationCoordinate()) {
76 0 : return ranges;
77 : }
78 0 : specification.trim();
79 0 : specification.ltrim('[');
80 0 : specification.rtrim(']');
81 :
82 0 : specification.upcase();
83 0 : if (specification == ALL) {
84 0 : ranges.push_back(0);
85 0 : ranges.push_back(nStokes - 1);
86 0 : return ranges;
87 : }
88 0 : if (specification.empty()) {
89 : // ranges.resize(2);
90 : // ranges[0] = 0;
91 0 : ranges.push_back(0);
92 0 : switch (stokesControl) {
93 0 : case USE_FIRST_STOKES:
94 0 : ranges.push_back(0);
95 0 : specification = firstStokes;
96 0 : break;
97 0 : case USE_ALL_STOKES:
98 0 : ranges.push_back(nStokes - 1);
99 0 : specification = ALL;
100 0 : break;
101 0 : default:
102 : // bug if we get here
103 0 : ThrowCc("Logic error, unhandled stokes control");
104 : break;
105 : };
106 0 : return ranges;
107 : }
108 : // First split on commas and semi-colons.
109 : // in the past for polarization specification.
110 :
111 0 : Vector<String> parts = stringToVector(specification, std::regex("[,;]"));
112 0 : Vector<String> polNames = Stokes::allNames(false);
113 0 : uInt nNames = polNames.size();
114 0 : Vector<uInt> nameLengths(nNames);
115 0 : for (uInt i=0; i<nNames; i++) {
116 0 : nameLengths[i] = polNames[i].length();
117 : }
118 0 : uInt *lengthData = nameLengths.data();
119 :
120 0 : Vector<uInt> idx(nNames);
121 0 : Sort sorter;
122 0 : sorter.sortKey(lengthData, TpUInt, 0, Sort::Descending);
123 0 : sorter.sort(idx, nNames);
124 :
125 0 : Vector<String> sortedNames(nNames);
126 0 : for (uInt i=0; i<nNames; i++) {
127 0 : sortedNames[i] = polNames[idx[i]];
128 0 : sortedNames[i].upcase();
129 : }
130 0 : for (uInt i=0; i<parts.size(); i++) {
131 0 : String part = parts[i];
132 0 : part.trim();
133 0 : Vector<String>::iterator iter = sortedNames.begin();
134 0 : while (iter != sortedNames.end() && ! part.empty()) {
135 0 : if (part.startsWith(*iter)) {
136 0 : Int stokesPix = csys.stokesPixelNumber(*iter);
137 0 : 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 0 : ranges.push_back(stokesPix);
144 0 : ranges.push_back(stokesPix);
145 : // consume the string
146 0 : part = part.substr(iter->length());
147 0 : if (! part.empty()) {
148 : // reset the iterator to start over at the beginning of the list for
149 : // the next specified polarization
150 0 : iter = sortedNames.begin();
151 : }
152 : }
153 : else {
154 0 : iter++;
155 : }
156 : }
157 0 : if (! part.empty()) {
158 0 : *_getLog() << "(Sub)String " << part << " in stokes specification part " << parts[i]
159 0 : << " does not match a known polarization." << LogIO::EXCEPTION;
160 : }
161 0 : }
162 : uInt nSel;
163 0 : return ParameterParser::consolidateAndOrderRanges(nSel, ranges);
164 0 : }
165 :
166 0 : Bool CasacRegionManager::_supports2DBox(Bool except) const {
167 0 : Bool ok = true;
168 0 : const CoordinateSystem& csys = getcoordsys();
169 0 : Vector<Int> axes;
170 0 : if (csys.hasDirectionCoordinate()) {
171 0 : axes = csys.directionAxesNumbers();
172 : }
173 0 : else if (csys.hasLinearCoordinate()) {
174 0 : axes = csys.linearAxesNumbers();
175 : }
176 : else {
177 0 : ok = false;
178 : }
179 0 : if (ok) {
180 0 : uInt nGood = 0;
181 0 : for (uInt i=0; i<axes.size(); i++) {
182 0 : if (axes[i] >= 0) {
183 0 : nGood++;
184 : }
185 : }
186 0 : if (nGood != 2) {
187 0 : ok = false;
188 : }
189 : }
190 0 : 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 0 : return ok;
195 0 : }
196 :
197 0 : vector<Double> CasacRegionManager::_setBoxCorners(const String& box) const {
198 0 : if (! box.empty()) {
199 0 : _supports2DBox(true);
200 : }
201 0 : Vector<String> boxParts = stringToVector(box);
202 0 : AlwaysAssert(boxParts.size() % 4 == 0, AipsError);
203 0 : vector<Double> corners(boxParts.size());
204 0 : for(uInt i=0; i<boxParts.size()/4; i++) {
205 0 : uInt index = 4*i;
206 0 : for (uInt j=0; j<4; j++) {
207 0 : boxParts[index + j].trim();
208 0 : 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 0 : corners[index + j] = String::toDouble(boxParts[index + j]);
213 : }
214 : }
215 0 : return corners;
216 0 : }
217 :
218 0 : 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 0 : LogOrigin origin("CasacRegionManager", __func__);
226 0 : Record regionRecord;
227 0 : if (! box.empty()) {
228 0 : ThrowIf(
229 : regionPtr != nullptr,
230 : "box and regionPtr cannot be simultaneously specified"
231 : );
232 0 : ThrowIf(
233 : box.freq(",") % 4 != 3,
234 : "box not specified correctly"
235 : );
236 0 : if (regionName.empty()) {
237 0 : regionRecord = fromBCS(
238 : diagnostics, nSelectedChannels, stokes,
239 : chans, stokesControl, box, imShape
240 0 : ).toRecord("");
241 0 : if (verbose) {
242 0 : *_getLog() << origin;
243 0 : *_getLog() << LogIO::NORMAL << "Using specified box(es) "
244 0 : << box << LogIO::POST;
245 : }
246 : }
247 : else {
248 0 : auto corners = stringToVector(box, ',');
249 0 : auto iter = corners.begin();
250 0 : auto end = corners.end();
251 0 : String crtfBoxString = "box[[";
252 0 : auto count = 0;
253 0 : for (; iter != end; ++iter, ++count) {
254 0 : iter->trim();
255 0 : if (count > 0 && count % 4 == 0) {
256 0 : crtfBoxString += "\nbox[[";
257 : }
258 0 : crtfBoxString += *iter + "pix";
259 0 : if (count % 2 == 0) {
260 0 : crtfBoxString += ",";
261 : }
262 : else {
263 : // close pixel pair
264 0 : crtfBoxString += "]";
265 0 : if (count - 1 % 4 == 0) {
266 : // first pixel pair done, start second pixel pair
267 0 : crtfBoxString += ",[";
268 : }
269 : else {
270 : // second pixel pair done, close box
271 0 : crtfBoxString += "]";
272 : }
273 : }
274 : }
275 0 : _setRegion(
276 : regionRecord, diagnostics, regionName, imShape,
277 : imageName, crtfBoxString, chans, stokes, verbose
278 : );
279 0 : }
280 : }
281 0 : else if (regionPtr) {
282 0 : ThrowIf(
283 : ! (regionName.empty() && chans.empty() && stokes.empty()),
284 : "regionPtr and regionName, chans, and/or stokes cannot "
285 : "be simultaneously specified"
286 : );
287 0 : _setRegion(regionRecord, diagnostics, regionPtr);
288 0 : stokes = _stokesFromRecord(regionRecord, stokesControl, imShape);
289 : }
290 0 : else if (! regionName.empty()) {
291 0 : _setRegion(
292 : regionRecord, diagnostics, regionName,
293 : imShape, imageName, "", chans, stokes, verbose
294 : );
295 0 : if (verbose) {
296 0 : *_getLog() << origin;
297 0 : *_getLog() << LogIO::NORMAL << diagnostics << LogIO::POST;
298 : }
299 0 : stokes = _stokesFromRecord(regionRecord, stokesControl, imShape);
300 : }
301 : else {
302 0 : vector<uInt> chanEndPts, polEndPts;
303 0 : regionRecord = fromBCS(
304 : diagnostics, nSelectedChannels, stokes,
305 : chans, stokesControl, box, imShape
306 0 : ).toRecord("");
307 0 : if (verbose) {
308 0 : *_getLog() << origin;
309 0 : *_getLog() << LogIO::NORMAL << "No directional region specified. Using full positional plane."
310 0 : << LogIO::POST;
311 : }
312 0 : const CoordinateSystem& csys = getcoordsys();
313 0 : if (csys.hasSpectralAxis()) {
314 0 : if (verbose) {
315 0 : if (chans.empty()) {
316 0 : *_getLog() << LogIO::NORMAL << "Using all spectral channels."
317 0 : << LogIO::POST;
318 : }
319 : else {
320 0 : *_getLog() << LogIO::NORMAL << "Using channel range(s) "
321 0 : << _pairsToString(chanEndPts) << LogIO::POST;
322 : }
323 : }
324 : }
325 0 : if (csys.hasPolarizationCoordinate() && verbose) {
326 0 : 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 0 : *_getLog() << LogIO::NORMAL << "Using polarizations " << stokes << LogIO::POST;
341 : }
342 : }
343 0 : }
344 0 : return regionRecord;
345 0 : }
346 :
347 0 : Record CasacRegionManager::regionFromString(
348 : const CoordinateSystem& csys, const String& regionStr,
349 : const String& imageName, const IPosition& imShape,
350 : Bool verbose
351 : ) {
352 0 : CasacRegionManager mgr(csys);
353 0 : Record reg;
354 0 : String diag;
355 0 : mgr._setRegion(
356 : reg, diag, regionStr, imShape, imageName, "", "", "", verbose
357 : );
358 0 : return reg;
359 0 : }
360 :
361 0 : void CasacRegionManager::_setRegion(
362 : Record& regionRecord, String& diagnostics,
363 : const Record* regionPtr
364 : ) {
365 : // region record pointer provided
366 0 : regionRecord = *(regionPtr);
367 : // set stokes from the region record
368 0 : diagnostics = "used provided region record";
369 0 : }
370 :
371 0 : 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 0 : if (regionName.empty() && imageName.empty()) {
381 0 : regionRecord = Record();
382 0 : diagnostics = "No region string";
383 0 : return;
384 : }
385 : // region name provided
386 0 : 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 0 : );
394 0 : File myFile(regionName);
395 0 : const CoordinateSystem csys = getcoordsys();
396 0 : if (myFile.exists()) {
397 0 : ThrowIf(
398 : ! myFile.isReadable(),
399 : "File " + regionName + " exists but is not readable."
400 : );
401 0 : std::unique_ptr<Record> rec;
402 : try {
403 0 : rec.reset(readImageFile(regionName, ""));
404 : }
405 0 : catch(const AipsError& x) {
406 0 : }
407 0 : if (rec) {
408 0 : 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 0 : regionRecord = *rec;
415 0 : diagnostics = "Region read from binary region file " + regionName;
416 0 : return;
417 : }
418 : try {
419 : // CRTF file attempt
420 : RegionTextList annList(
421 : regionName, csys, imShape, prependBox,
422 : globalOverrideChans, globalStokesOverride
423 0 : );
424 0 : regionRecord = annList.regionAsRecord();
425 0 : diagnostics = "Region read from CRTF file " + regionName;
426 0 : }
427 0 : catch (const AipsError& x) {
428 0 : ThrowCc(
429 : regionName + " is neither a valid binary region file, "
430 : "nor a valid region text file."
431 : );
432 0 : }
433 0 : }
434 0 : 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 0 : );
441 0 : regionRecord = annList.regionAsRecord();
442 0 : diagnostics = "Region read from text string " + regionName;
443 0 : }
444 0 : catch (const AipsError& x) {
445 0 : ThrowCc(x.getMesg());
446 0 : }
447 : }
448 0 : else if (regionName.matches(image) || ! imageName.empty()) {
449 0 : ImageRegion imRegion;
450 0 : String imagename, region;
451 0 : if (regionName.matches(image)) {
452 0 : String res[2];
453 0 : casacore::split(regionName, res, 2, ":");
454 0 : imagename = res[0];
455 0 : region = res[1];
456 0 : }
457 : else {
458 : // imageName is not empty if we get here
459 0 : imagename = imageName;
460 0 : region = regionName;
461 : }
462 : try {
463 0 : Record *myRec = tableToRecord(imagename, region);
464 0 : 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 0 : if (Table::isReadable(imagename)) {
471 0 : ThrowIf(
472 : myRec == 0,
473 : "Region " + region + " not found in image "
474 : + imagename
475 : );
476 0 : regionRecord = *myRec;
477 0 : diagnostics = "Used region " + region + " from image "
478 0 : + imagename + " table description";
479 : }
480 : else {
481 0 : *_getLog() << "Cannot read image " << imagename
482 0 : << " to get region " << region << LogIO::EXCEPTION;
483 : }
484 : }
485 0 : catch (const AipsError&) {
486 0 : ThrowCc(
487 : "Unable to open region file or region table description "
488 : + region + " in image " + imagename
489 : );
490 0 : }
491 0 : }
492 : else {
493 0 : ostringstream oss;
494 0 : oss << "Unable to open region file or region table description "
495 0 : << regionName << "." << endl
496 0 : << "If it is supposed to be a text string its format is incorrect";
497 0 : ThrowCc(oss.str());
498 0 : }
499 0 : }
500 :
501 0 : 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 0 : const CoordinateSystem& csys = getcoordsys();
508 : vector<uInt> chanEndPts = setSpectralRanges(
509 : chans, nSelectedChannels, imShape
510 0 : );
511 0 : Int polAxisNumber = csys.polarizationAxisNumber();
512 0 : uInt nTotalPolarizations = polAxisNumber >= 0 ? imShape[polAxisNumber] : 0;
513 0 : String firstStokes = polAxisNumber >= 0 ? csys.stokesAtPixel(0) : "";
514 : vector<uInt> polEndPts = _setPolarizationRanges(
515 : stokes, firstStokes,
516 : nTotalPolarizations, stokesControl
517 0 : );
518 0 : vector<Double> boxCorners;
519 0 : if (box.empty()) {
520 0 : if (_supports2DBox(false)) {
521 0 : if (
522 0 : csys.hasDirectionCoordinate()
523 0 : || csys.hasLinearCoordinate()
524 : ) {
525 0 : Vector<Int> dirAxesNumbers;
526 0 : if (csys.hasDirectionCoordinate()) {
527 0 : dirAxesNumbers = csys.directionAxesNumbers();
528 : }
529 : else {
530 0 : dirAxesNumbers = csys.linearAxesNumbers();
531 : }
532 0 : Vector<Int> dirShape(2);
533 0 : dirShape[0] = imShape[dirAxesNumbers[0]];
534 0 : dirShape[1] = imShape[dirAxesNumbers[1]];
535 0 : boxCorners.resize(4);
536 0 : boxCorners[0] = 0;
537 0 : boxCorners[1] = 0;
538 0 : boxCorners[2] = dirShape[0] - 1;
539 0 : boxCorners[3] = dirShape[1] - 1;
540 0 : }
541 : }
542 : }
543 : else {
544 0 : boxCorners = _setBoxCorners(box);
545 : }
546 : return _fromBCS(
547 : diagnostics, boxCorners,
548 : chanEndPts, polEndPts, imShape
549 0 : );
550 0 : }
551 :
552 0 : 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 0 : LogOrigin origin("CasacRegionManager", __func__);
558 0 : *_getLog() << origin;
559 0 : Vector<Double> blc(imShape.nelements(), 0);
560 0 : Vector<Double> trc(imShape.nelements(), 0);
561 0 : const CoordinateSystem csys = getcoordsys();
562 0 : Vector<Int> directionAxisNumbers = csys.directionAxesNumbers();
563 0 : 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 0 : vector<Int>::iterator iter = linearAxisNumbers.begin();
567 0 : vector<Int>::iterator end = linearAxisNumbers.end();
568 0 : while(iter != end) {
569 0 : if (*iter < 0) {
570 0 : iter = linearAxisNumbers.erase(iter);
571 : }
572 0 : ++iter;
573 : }
574 0 : Int spectralAxisNumber = csys.spectralAxisNumber();
575 0 : Int polarizationAxisNumber = csys.polarizationAxisNumber();
576 :
577 0 : Vector<Double> xCorners(boxCorners.size()/2);
578 0 : Vector<Double> yCorners(xCorners.size());
579 0 : for (uInt i=0; i<xCorners.size(); i++) {
580 0 : Double x = boxCorners[2*i];
581 0 : Double y = boxCorners[2*i + 1];
582 :
583 0 : if (x < 0 || y < 0 ) {
584 0 : *_getLog() << "blc in box spec is less than 0" << LogIO::EXCEPTION;
585 : }
586 0 : if (csys.hasDirectionCoordinate()) {
587 0 : if (
588 0 : x >= imShape[directionAxisNumbers[0]]
589 0 : || y >= imShape[directionAxisNumbers[1]]
590 : ) {
591 0 : *_getLog() << "dAxisNum0=" <<directionAxisNumbers[0] <<" dAxisNum1="<<directionAxisNumbers[1];
592 0 : *_getLog() << "x="<<x<<" imShape[0]="<<imShape[directionAxisNumbers[0]]<< " y="<<y<<" imShape[1]="<<imShape[directionAxisNumbers[1]]<<LogIO::POST;
593 0 : *_getLog() << "trc in box spec is greater than or equal to number "
594 0 : << "of direction coordinate pixels in the image" << LogIO::EXCEPTION;
595 : }
596 : }
597 0 : else if (
598 0 : csys.hasLinearCoordinate()
599 0 : && (
600 0 : x >= imShape[linearAxisNumbers[0]]
601 0 : || 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 0 : xCorners[i] = x;
608 0 : yCorners[i] = y;
609 : }
610 0 : Vector<Double> polEndPtsDouble(polEndPts.size());
611 0 : for (uInt i=0; i<polEndPts.size(); ++i) {
612 0 : polEndPtsDouble[i] = (Double)polEndPts[i];
613 : }
614 :
615 0 : Bool csysSupports2DBox = _supports2DBox(false);
616 0 : uInt nRegions = 1;
617 0 : if (csysSupports2DBox) {
618 0 : if (csys.hasDirectionCoordinate()) {
619 0 : nRegions *= boxCorners.size()/4;
620 : }
621 0 : if (csys.hasLinearCoordinate()) {
622 0 : nRegions *= boxCorners.size()/4;
623 : }
624 : }
625 0 : if (csys.hasPolarizationCoordinate()) {
626 0 : nRegions *= polEndPts.size()/2;
627 : }
628 0 : if (csys.hasSpectralAxis()) {
629 0 : nRegions *= chanEndPts.size()/2;
630 : }
631 0 : Vector<Double> extXCorners(2*nRegions, 0);
632 0 : Vector<Double> extYCorners(2*nRegions, 0);
633 0 : Vector<Double> extPolEndPts(2*nRegions, 0);
634 0 : Vector<Double> extChanEndPts(2*nRegions, 0);
635 :
636 0 : uInt count = 0;
637 :
638 0 : if (csysSupports2DBox) {
639 0 : for (uInt i=0; i<max(uInt(1), xCorners.size()/2); i++) {
640 0 : for (uInt j=0; j<max((uInt)1, polEndPts.size()/2); j++) {
641 0 : for (uInt k=0; k<max(uInt(1), chanEndPts.size()/2); k++) {
642 0 : if (
643 0 : csys.hasDirectionCoordinate()
644 0 : || csys.hasLinearCoordinate()
645 : ) {
646 0 : extXCorners[2*count] = xCorners[2*i];
647 0 : extXCorners[2*count + 1] = xCorners[2*i + 1];
648 0 : extYCorners[2*count] = yCorners[2*i];
649 0 : extYCorners[2*count + 1] = yCorners[2*i + 1];
650 : }
651 0 : if (csys.hasPolarizationCoordinate()) {
652 0 : extPolEndPts[2*count] = polEndPtsDouble[2*j];
653 0 : extPolEndPts[2*count + 1] = polEndPtsDouble[2*j + 1];
654 :
655 : }
656 0 : if (csys.hasSpectralAxis()) {
657 0 : extChanEndPts[2*count] = chanEndPts[2*k];
658 0 : extChanEndPts[2*count + 1] = chanEndPts[2*k + 1];
659 : }
660 0 : 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 0 : for (uInt j=0; j<max((uInt)1, polEndPts.size()/2); j++) {
669 0 : for (uInt k=0; k<max(uInt(1), chanEndPts.size()/2); k++) {
670 0 : if (csys.hasPolarizationCoordinate()) {
671 0 : extPolEndPts[2*count] = polEndPtsDouble[2*j];
672 0 : extPolEndPts[2*count + 1] = polEndPtsDouble[2*j + 1];
673 : }
674 0 : if (csys.hasSpectralAxis()) {
675 0 : extChanEndPts[2*count] = chanEndPts[2*k];
676 0 : extChanEndPts[2*count + 1] = chanEndPts[2*k + 1];
677 : }
678 0 : count++;
679 : }
680 : }
681 : }
682 0 : map<uInt, Vector<Double> > axisCornerMap;
683 0 : for (uInt i=0; i<nRegions; i++) {
684 0 : for (uInt axisNumber=0; axisNumber<csys.nPixelAxes(); axisNumber++) {
685 0 : if (
686 : (
687 0 : directionAxisNumbers.size() > 1
688 0 : && (Int)axisNumber == directionAxisNumbers[0]
689 : )
690 0 : || (
691 0 : ! csys.hasDirectionCoordinate()
692 0 : && linearAxisNumbers.size() > 1
693 0 : && (Int)axisNumber == linearAxisNumbers[0]
694 : )
695 : ) {
696 0 : axisCornerMap[axisNumber] = extXCorners;
697 : }
698 0 : else if (
699 : (
700 0 : directionAxisNumbers.size() > 1
701 0 : && (Int)axisNumber == directionAxisNumbers[1]
702 : )
703 0 : || (
704 0 : ! csys.hasDirectionCoordinate()
705 0 : && linearAxisNumbers.size() > 1
706 0 : && (Int)axisNumber == linearAxisNumbers[1]
707 : )
708 : ) {
709 0 : axisCornerMap[axisNumber] = extYCorners;
710 : }
711 0 : else if ((Int)axisNumber == spectralAxisNumber) {
712 0 : axisCornerMap[axisNumber] = extChanEndPts;
713 : }
714 0 : else if ((Int)axisNumber == polarizationAxisNumber) {
715 0 : axisCornerMap[axisNumber] = extPolEndPts;
716 : }
717 : else {
718 0 : Vector<Double> range(2, 0);
719 0 : range[1] = imShape[axisNumber] - 1;
720 0 : axisCornerMap[axisNumber] = range;
721 0 : }
722 : }
723 : }
724 0 : ImageRegion imRegion;
725 0 : for (uInt i=0; i<nRegions; i++) {
726 0 : for (uInt axisNumber=0; axisNumber<csys.nPixelAxes(); axisNumber++) {
727 0 : blc(axisNumber) = axisCornerMap[axisNumber][2*i];
728 0 : trc(axisNumber) = axisCornerMap[axisNumber][2*i + 1];
729 : }
730 0 : LCBox lcBox(blc, trc, imShape);
731 0 : WCBox wcBox(lcBox, csys);
732 0 : ImageRegion thisRegion(wcBox);
733 : imRegion = (i == 0)
734 : ? thisRegion
735 0 : : imRegion = *(doUnion(imRegion, thisRegion));
736 0 : }
737 0 : ostringstream os;
738 0 : os << "Used image region from " << endl;
739 0 : if (csys.hasDirectionCoordinate()) {
740 0 : os << " position box corners: ";
741 0 : for (uInt i=0; i<boxCorners.size()/4; i++) {
742 0 : os << boxCorners[4*i] << ", " << boxCorners[4*i + 1]
743 0 : << ", " << boxCorners[4*i + 2] << ", " << boxCorners[4*i + 3];
744 0 : if (i < boxCorners.size()/4 - 1) {
745 0 : os << "; ";
746 : }
747 : }
748 : }
749 0 : if (getcoordsys().hasSpectralAxis()) {
750 0 : os << " spectral channel ranges: " << _pairsToString(chanEndPts);
751 : }
752 0 : if (getcoordsys().hasPolarizationCoordinate()) {
753 0 : os << " polarization pixel ranges: " << _pairsToString(polEndPts);
754 : }
755 0 : diagnostics = os.str();
756 0 : return imRegion;
757 0 : }
758 :
759 0 : String CasacRegionManager::_pairsToString(const vector<uInt>& pairs) const {
760 0 : ostringstream os;
761 0 : uInt nPairs = pairs.size()/2;
762 0 : for (uInt i=0; i<nPairs; i++) {
763 0 : os << pairs[2*i] << " - " << pairs[2*i + 1];
764 0 : if (i < nPairs - 1) {
765 0 : os << "; ";
766 : }
767 : }
768 0 : return os.str();
769 0 : }
770 :
771 0 : 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 0 : String stokes = "";
776 0 : CoordinateSystem csys = getcoordsys();
777 0 : if(! csys.hasPolarizationCoordinate()) {
778 0 : return stokes;
779 : }
780 0 : Int polAxis = csys.polarizationAxisNumber();
781 0 : if (shape[polAxis] == 1) {
782 : // degenerate stokes axis
783 0 : return csys.stokesAtPixel(0);
784 : }
785 0 : uInt stokesBegin = 0;
786 0 : uInt stokesEnd = 0;
787 0 : if (region.empty()) {
788 0 : stokesEnd = stokesControl == USE_FIRST_STOKES
789 : ? 0
790 0 : : csys.stokesCoordinate().stokes().size() - 1;
791 : }
792 : else {
793 0 : std::unique_ptr<ImageRegion> imreg(ImageRegion::fromRecord(region, ""));
794 0 : Array<Float> blc, trc;
795 0 : Bool oneRelAccountedFor = false;
796 0 : if (imreg->isLCSlicer()) {
797 0 : blc = imreg->asLCSlicer().blc();
798 0 : 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 0 : trc = imreg->asLCSlicer().trc();
805 0 : 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 0 : stokesBegin = (uInt)((Vector<Float>)blc)[polAxis];
812 0 : stokesEnd = (uInt)((Vector<Float>)trc)[polAxis];
813 0 : oneRelAccountedFor = true;
814 : }
815 0 : else if (
816 0 : RegionManager::isPixelRegion(
817 0 : *(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 0 : else if (region.fieldNumber("x") >= 0 && region.fieldNumber("y") >= 0) {
826 : // world polygon
827 0 : oneRelAccountedFor = true;
828 0 : stokesBegin = 0;
829 0 : stokesEnd = stokesControl == USE_FIRST_STOKES
830 0 : ? 0 : shape[polAxis];
831 : }
832 0 : else if (region.fieldNumber("blc") >= 0 && region.fieldNumber("blc") >= 0) {
833 : // world box
834 0 : Record blcRec = region.asRecord("blc");
835 0 : Record trcRec = region.asRecord("trc");
836 0 : String polField = "*" + String::toString(polAxis + 1);
837 0 : stokesBegin = blcRec.isDefined(polField)
838 0 : ? (Int)blcRec.asRecord(
839 0 : String(polField)
840 0 : ).asDouble("value")
841 : : 0;
842 0 : stokesEnd = trcRec.isDefined(polField)
843 0 : ? (Int)blcRec.asRecord(
844 0 : String(polField)
845 0 : ).asDouble("value")
846 : : stokesControl == USE_FIRST_STOKES
847 : ? 0
848 0 : : shape[polAxis];
849 0 : if (! blcRec.isDefined(polField)) {
850 0 : oneRelAccountedFor = true;
851 : }
852 0 : }
853 : else {
854 : // FIXME not very nice, but until all can be implemented this will have to do
855 0 : *_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 0 : << LogIO::POST;
858 0 : return stokes;
859 : }
860 0 : if (
861 0 : ! oneRelAccountedFor
862 0 : && region.isDefined("oneRel")
863 0 : && region.asBool("oneRel")
864 : ) {
865 0 : stokesBegin--;
866 0 : stokesEnd--;
867 : }
868 0 : }
869 0 : for (uInt i=stokesBegin; i<=stokesEnd; i++) {
870 0 : stokes += csys.stokesAtPixel(i);
871 : }
872 0 : return stokes;
873 0 : }
874 :
875 0 : vector<uInt> CasacRegionManager::_initSpectralRanges(
876 : uInt& nSelectedChannels, const IPosition& imShape
877 : ) const {
878 0 : vector<uInt> ranges(0);
879 0 : if (! getcoordsys().hasSpectralAxis()) {
880 0 : nSelectedChannels = 0;
881 0 : return ranges;
882 : }
883 0 : uInt specNum = getcoordsys().spectralAxisNumber();
884 0 : 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 0 : uInt nChannels = imShape[getcoordsys().spectralAxisNumber()];
891 0 : ranges.push_back(0);
892 0 : ranges.push_back(nChannels - 1);
893 0 : nSelectedChannels = nChannels;
894 0 : return ranges;
895 0 : }
896 :
897 0 : vector<uInt> CasacRegionManager::setSpectralRanges(
898 : uInt& nSelectedChannels,
899 : const Record *const regionRec, const IPosition& imShape
900 : ) const {
901 0 : if (regionRec == 0 || ! getcoordsys().hasSpectralAxis()) {
902 0 : return _initSpectralRanges(nSelectedChannels, imShape);
903 : }
904 : else {
905 0 : return _spectralRangeFromRegionRecord(nSelectedChannels, regionRec, imShape);
906 : }
907 : }
908 :
909 0 : vector<uInt> CasacRegionManager::setSpectralRanges(
910 : String specification, uInt& nSelectedChannels,
911 : /*const String& globalChannelOverride,
912 : const String& globalStokesOverride, */ const IPosition& imShape
913 : ) const {
914 0 : LogOrigin origin("CasacRegionManager", __func__);
915 0 : *_getLog() << origin;
916 0 : specification.trim();
917 0 : String x = specification;
918 0 : x.upcase();
919 0 : if (x.empty() || x == ALL) {
920 0 : return _initSpectralRanges(nSelectedChannels, imShape);
921 : }
922 0 : 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 0 : else if (specification.contains("range")) {
930 : // this is a specification in the "new" ASCII region format
931 0 : return _spectralRangeFromRangeFormat(nSelectedChannels, specification, imShape);
932 : }
933 : else {
934 0 : uInt nChannels = imShape[getcoordsys().spectralAxisNumber()];
935 : return ParameterParser::spectralRangesFromChans(
936 : nSelectedChannels, specification, nChannels
937 0 : );
938 : }
939 0 : }
940 :
941 0 : vector<uInt> CasacRegionManager::_spectralRangeFromRegionRecord(
942 : uInt& nSelectedChannels, const Record *const regionRec,
943 : const IPosition& imShape
944 : ) const {
945 0 : const CoordinateSystem& csys = getcoordsys();
946 0 : TempImage<Float> x(imShape, csys);
947 0 : x.set(0);
948 0 : SPCIIF subimage = SubImageFactory<Float>::createSubImageRO(
949 0 : x, *regionRec, "", _getLog(), AxesSpecifier(), false, true
950 0 : );
951 0 : uInt nChan = 0;
952 : {
953 0 : ImageMetaData<Float> md(subimage);
954 0 : nChan = md.nChannels();
955 0 : }
956 0 : const SpectralCoordinate& subsp = subimage->coordinates().spectralCoordinate();
957 : Double subworld;
958 0 : subsp.toWorld(subworld, 0);
959 0 : const SpectralCoordinate& imsp = csys.spectralCoordinate();
960 : Double pixOff;
961 0 : imsp.toPixel(pixOff, subworld);
962 0 : Int specAxisNumber = csys.spectralAxisNumber();
963 0 : IPosition start(subimage->ndim(), 0);
964 0 : ArrayLattice<Bool> mask(subimage->getMask());
965 0 : IPosition planeShape = subimage->shape();
966 0 : vector<uInt> myList;
967 0 : planeShape[specAxisNumber] = 1;
968 0 : for (uInt i=0; i<nChan; i++) {
969 0 : start[specAxisNumber] = i;
970 0 : Array<Bool> maskSlice;
971 0 : mask.getSlice(maskSlice, start, planeShape);
972 0 : if (anyTrue(maskSlice)) {
973 0 : uInt real = i + (uInt)rint(pixOff);
974 0 : myList.push_back(real);
975 0 : myList.push_back(real);
976 : }
977 0 : }
978 0 : return ParameterParser::consolidateAndOrderRanges(nSelectedChannels, myList);
979 0 : }
980 :
981 0 : 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 0 : const CoordinateSystem csys = getcoordsys();
989 : RegionTextParser::ParamSet parms = RegionTextParser::getParamSet(
990 0 : spectralParmsUpdated, *_getLog(),
991 0 : specification, "", csys, std::shared_ptr<std::pair<MFrequency, MFrequency> >(nullptr),
992 0 : std::shared_ptr<Vector<Stokes::StokesTypes> >(nullptr)
993 0 : );
994 0 : RegionTextParser::ParamSet::const_iterator end = parms.end();
995 0 : for (
996 0 : RegionTextParser::ParamSet::const_iterator iter=parms.begin();
997 0 : iter!=end; iter++
998 : ) {
999 0 : AnnotationBase::Keyword key = iter->first;
1000 0 : 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 0 : String regSpec = "box[[0pix, 0pix], [1pix, 1pix]] " + specification;
1011 0 : RegionTextParser parser(csys, imShape, regSpec);
1012 0 : vector<uInt> range(2);
1013 0 : ThrowIf(
1014 : parser.getLines().empty(),
1015 : "The specified spectral range " + specification
1016 : + " does not intersect the image spectral range."
1017 : );
1018 0 : auto ann = parser.getLines()[0].getAnnotationBase();
1019 0 : /*const AnnRegion*/ auto *reg = dynamic_cast<const AnnRegion*>(ann.get());
1020 0 : ThrowIf(! reg, "Dynamic cast failed");
1021 0 : /*vector<Double>*/ const auto drange = reg->getSpectralPixelRange();
1022 0 : range[0] = uInt(max(0.0, floor(drange[0] + 0.5)));
1023 0 : range[1] = uInt(floor(drange[1] + 0.5));
1024 0 : nSelectedChannels = range[1] - range[0] + 1;
1025 0 : return range;
1026 0 : }
1027 :
1028 : }
|
