Line data Source code
1 : #include <imageanalysis/Annotations/AnnRegion.h>
2 :
3 : #include <casacore/casa/Exceptions/Error.h>
4 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
5 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
6 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
7 : //#include <imageanalysis/Regions/CasacRegionManager.h>
8 : #include <imageanalysis/IO/ParameterParser.h>
9 : #include <casacore/images/Regions/WCExtension.h>
10 : #include <casacore/images/Regions/WCUnion.h>
11 : #include <casacore/tables/Tables/TableRecord.h>
12 :
13 : #include <iomanip>
14 :
15 : using namespace casacore;
16 : namespace casa {
17 :
18 : const String AnnRegion::_class = "AnnRegion";
19 :
20 139 : AnnRegion::AnnRegion(
21 : const Type shape, const String& dirRefFrameString,
22 : const CoordinateSystem& csys,
23 : const IPosition& imShape,
24 : const Quantity& beginFreq,
25 : const Quantity& endFreq,
26 : const String& freqRefFrame,
27 : const String& dopplerString,
28 : const Quantity& restfreq,
29 : const Vector<Stokes::StokesTypes> stokes,
30 : const Bool annotationOnly,
31 : const Bool requireImageRegion
32 139 : ) : AnnotationBase(
33 : shape, dirRefFrameString, csys, beginFreq, endFreq,
34 : freqRefFrame, dopplerString, restfreq, stokes
35 : ),
36 139 : _requireImageRegion(requireImageRegion),
37 139 : _isAnnotationOnly(annotationOnly),
38 139 : _isDifference(false), _constructing(true), _imShape(imShape),
39 278 : _spectralPixelRange(vector<Double>(0)) {
40 139 : _init();
41 : // just before returning
42 139 : _constructing = false;
43 139 : }
44 :
45 427 : AnnRegion::AnnRegion(
46 : const Type shape,
47 : const CoordinateSystem& csys,
48 : const IPosition& imShape,
49 : const Vector<Stokes::StokesTypes>& stokes,
50 : const Bool requireImageRegion
51 427 : ) : AnnotationBase(shape, csys, stokes),
52 427 : _requireImageRegion(requireImageRegion),
53 427 : _isDifference(false), _constructing(true), _imShape(imShape),
54 854 : _spectralPixelRange(vector<Double>(0)) {
55 427 : _init();
56 : // just before returning
57 427 : _constructing = false;
58 427 : }
59 :
60 0 : AnnRegion::AnnRegion(const AnnRegion& other)
61 : : AnnotationBase(other),
62 0 : _requireImageRegion(other._requireImageRegion),
63 0 : _imageRegion(other._imageRegion),
64 0 : _directionRegion(other._directionRegion),
65 0 : _isAnnotationOnly(other._isAnnotationOnly),
66 0 : _isDifference(other._isDifference),
67 0 : _constructing(other._constructing),
68 0 : _imShape(other._imShape),
69 0 : _spectralPixelRange(other._spectralPixelRange) {}
70 :
71 :
72 566 : AnnRegion::~AnnRegion() {}
73 :
74 0 : AnnRegion& AnnRegion::operator= (const AnnRegion& other) {
75 0 : if (&other != this) {
76 0 : AnnotationBase::operator= (other);
77 0 : _isAnnotationOnly = other._isAnnotationOnly;
78 0 : _requireImageRegion = other._requireImageRegion;
79 0 : _imageRegion = other._imageRegion;
80 0 : _directionRegion = other._directionRegion;
81 0 : _isDifference = other._isDifference;
82 0 : _constructing = other._constructing;
83 0 : _imShape.resize(other._imShape.size());
84 0 : _imShape = other._imShape;
85 0 : _spectralPixelRange = other._spectralPixelRange;
86 : }
87 0 : return *this;
88 : }
89 :
90 0 : Bool AnnRegion::operator== (const AnnRegion& other) const {
91 0 : return &other == this || (
92 0 : _isAnnotationOnly == other._isAnnotationOnly
93 0 : && _requireImageRegion == other._requireImageRegion
94 0 : && _imageRegion == other._imageRegion
95 0 : && _directionRegion == other._directionRegion
96 0 : && _isDifference == other._isDifference
97 0 : && _constructing == other._constructing
98 0 : && _imShape == other._imShape
99 0 : && _spectralPixelRange == other._spectralPixelRange
100 0 : );
101 : }
102 :
103 :
104 0 : void AnnRegion::setAnnotationOnly(const Bool isAnnotationOnly) {
105 0 : _isAnnotationOnly = isAnnotationOnly;
106 0 : }
107 :
108 138 : Bool AnnRegion::isAnnotationOnly() const {
109 138 : return _isAnnotationOnly;
110 : }
111 :
112 566 : Bool AnnRegion::_hasDirectionRegion() {
113 566 : return (_directionRegion.isWCRegion() || _directionRegion.isLCRegion() || _directionRegion.isLCSlicer());
114 : }
115 :
116 565 : Bool AnnRegion::hasImageRegion() const {
117 565 : return (_imageRegion.isWCRegion() || _imageRegion.isLCRegion() || _imageRegion.isLCSlicer());
118 : }
119 :
120 139 : void AnnRegion::setDifference(const Bool difference) {
121 139 : _isDifference = difference;
122 139 : }
123 :
124 276 : Bool AnnRegion::isDifference() const {
125 276 : return _isDifference;
126 : }
127 :
128 0 : TableRecord AnnRegion::asRecord() const {
129 0 : return _imageRegion.toRecord("");
130 : }
131 :
132 0 : ImageRegion AnnRegion::asImageRegion() const {
133 0 : return _imageRegion;
134 : }
135 :
136 427 : CountedPtr<const WCRegion> AnnRegion::getRegion() const {
137 427 : if (hasImageRegion()) {
138 427 : return _imageRegion.asWCRegionPtr()->cloneRegion();
139 : } else {
140 0 : return CountedPtr<const WCRegion>(nullptr);
141 : }
142 : }
143 :
144 138 : std::shared_ptr<const WCRegion> AnnRegion::getRegion2() const {
145 138 : if (hasImageRegion()) {
146 138 : return std::shared_ptr<const WCRegion>(_imageRegion.asWCRegionPtr()->cloneRegion());
147 : } else {
148 0 : return std::shared_ptr<const WCRegion>(nullptr);
149 : }
150 : }
151 :
152 277 : Bool AnnRegion::isRegion() const {
153 277 : return true;
154 : }
155 :
156 566 : void AnnRegion::_init() {
157 566 : if (_imShape.nelements() != getCsys().nPixelAxes()) {
158 0 : ostringstream oss;
159 0 : oss << _class << "::" << __FUNCTION__ << ": Number of coordinate axes ("
160 0 : << getCsys().nPixelAxes() << ") differs from number of dimensions in image shape ("
161 0 : << _imShape.nelements() << ")";
162 0 : throw AipsError(oss.str());
163 0 : }
164 :
165 566 : }
166 :
167 0 : Bool AnnRegion::setFrequencyLimits(
168 : const Quantity& beginFreq,
169 : const Quantity& endFreq,
170 : const String& freqRefFrame,
171 : const String& dopplerString,
172 : const Quantity& restfreq
173 : ) {
174 0 : if (
175 0 : AnnotationBase::setFrequencyLimits(
176 : beginFreq, endFreq, freqRefFrame, dopplerString, restfreq
177 : )
178 : ) {
179 0 : if (! _constructing) {
180 : // have to re-extend the direction region because of new freq range.
181 : // but not during object construction
182 0 : _extend();
183 : }
184 0 : return true;
185 : }
186 : else {
187 0 : return false;
188 : }
189 : }
190 :
191 566 : void AnnRegion::_setDirectionRegion(const ImageRegion& region) {
192 566 : _directionRegion = region;
193 566 : }
194 :
195 :
196 1 : vector<Double> AnnRegion::getSpectralPixelRange() const {
197 1 : return _spectralPixelRange;
198 : }
199 :
200 :
201 566 : void AnnRegion::_extend() {
202 566 : if (!_hasDirectionRegion()) {
203 0 : return; // nothing to extend
204 : }
205 566 : Int stokesAxis = -1;
206 566 : Int spectralAxis = -1;
207 566 : Vector<Quantity> freqRange;
208 566 : uInt nBoxes = 0;
209 566 : Vector<MFrequency> freqLimits = getFrequencyLimits();
210 566 : const CoordinateSystem& csys = getCsys();
211 566 : if (csys.hasSpectralAxis() && freqLimits.size() == 2) {
212 11 : SpectralCoordinate spcoord = csys.spectralCoordinate();
213 11 : spectralAxis = csys.spectralAxisNumber();
214 :
215 11 : if (spectralAxis < 0) {
216 0 : throw AipsError(
217 0 : String(__FUNCTION__) + ": A spectral range was specified "
218 0 : + "but the spectral pixel axis in the supplied coordinate "
219 0 : + "system is not present."
220 0 : );
221 : }
222 :
223 11 : String unit = spcoord.worldAxisUnits()[0];
224 11 : _spectralPixelRange.resize(2);
225 11 : spcoord.toPixel(_spectralPixelRange[0], freqLimits[0]);
226 11 : spcoord.toPixel(_spectralPixelRange[1], freqLimits[1]);
227 11 : freqRange.resize(2);
228 11 : freqRange[0] = freqLimits[0].get(unit);
229 11 : freqRange[1] = freqLimits[1].get(unit);
230 11 : if (_spectralPixelRange[1] <= _spectralPixelRange[0]) {
231 0 : std::swap(_spectralPixelRange[0], _spectralPixelRange[1]);
232 0 : std::swap(freqRange[0], freqRange[1]);
233 : }
234 11 : nBoxes = 1;
235 11 : }
236 566 : vector<Stokes::StokesTypes> stokesRanges;
237 566 : Vector<Stokes::StokesTypes> stokes = getStokes();
238 566 : if (
239 733 : csys.hasPolarizationCoordinate() && stokes.size() > 0
240 733 : && (stokesAxis = csys.polarizationAxisNumber()) >= 0
241 : ) {
242 3 : vector<uInt> stokesNumbers(2*stokes.size());
243 7 : for (uInt i=0; i<stokes.size(); i++) {
244 4 : stokesNumbers[2*i] = (uInt)stokes[i];
245 4 : stokesNumbers[2*i + 1] = stokesNumbers[2*i];
246 : }
247 3 : uInt nSel = 0;
248 : vector<uInt> orderedRanges = ParameterParser::consolidateAndOrderRanges(
249 : nSel, stokesNumbers
250 3 : );
251 11 : for (uInt i=0; i<orderedRanges.size(); i++) {
252 8 : stokesRanges.push_back(Stokes::type(orderedRanges[i]));
253 : }
254 3 : nBoxes = stokesRanges.size()/2;
255 3 : }
256 566 : if (nBoxes == 0) {
257 552 : _imageRegion = _directionRegion;
258 : }
259 : else {
260 14 : uInt nExtendAxes = 0;
261 14 : if (spectralAxis >= 0) {
262 11 : nExtendAxes++;
263 : }
264 14 : if (stokesAxis >= 0) {
265 3 : nExtendAxes++;
266 : }
267 :
268 14 : IPosition pixelAxes(nExtendAxes);
269 14 : uInt n = 0;
270 : // spectral axis must be first to be consistent with _makeExtensionBox()
271 14 : if (spectralAxis > 0) {
272 11 : pixelAxes[n] = spectralAxis;
273 11 : n++;
274 : }
275 14 : if (stokesAxis > 0) {
276 3 : pixelAxes[n] = stokesAxis;
277 3 : n++;
278 : }
279 14 : if (nBoxes == 1) {
280 13 : Vector<Stokes::StokesTypes> const stokesRangesV(stokesRanges);
281 13 : WCBox wbox = _makeExtensionBox(freqRange, stokesRangesV, pixelAxes);
282 13 : _imageRegion = ImageRegion(WCExtension(_directionRegion, wbox));
283 13 : }
284 : else {
285 1 : PtrBlock<const WCRegion*> regions(nBoxes);
286 3 : for (uInt i=0; i<nBoxes; i++) {
287 2 : Vector<Stokes::StokesTypes> stokesRange(2);
288 2 : stokesRange[0] = stokesRanges[2*i];
289 2 : stokesRange[1] = stokesRanges[2*i + 1];
290 2 : WCBox wbox = _makeExtensionBox(freqRange, stokesRange, pixelAxes);
291 2 : regions[i] = new WCExtension(_directionRegion, wbox);
292 2 : }
293 1 : _imageRegion = ImageRegion(WCUnion(true, regions));
294 1 : }
295 14 : }
296 : try {
297 566 : _imageRegion.asWCRegionPtr()->toLCRegion(csys, _imShape);
298 : }
299 0 : catch (const AipsError& x) {
300 0 : throw (ToLCRegionConversionError(x.getMesg()));
301 0 : }
302 566 : }
303 :
304 15 : WCBox AnnRegion::_makeExtensionBox(
305 : const Vector<Quantity>& freqRange,
306 : const Vector<Stokes::StokesTypes>& stokesRange,
307 : const IPosition& pixelAxes
308 : ) const {
309 15 : uInt n = 0;
310 15 : Vector<Quantity> blc(pixelAxes.size());
311 15 : Vector<Quantity> trc(pixelAxes.size());
312 15 : Vector<Int> absRel(pixelAxes.size(), RegionType::Abs);
313 15 : if (freqRange.size() == 2) {
314 11 : blc[n] = freqRange[0];
315 11 : trc[n] = freqRange[1];
316 11 : n++;
317 : }
318 15 : if (stokesRange.size() == 2) {
319 4 : blc[n] = Quantity(stokesRange[0], "");
320 4 : trc[n] = Quantity(stokesRange[1], "");
321 : }
322 15 : WCBox wbox(blc, trc, pixelAxes, getCsys(), absRel);
323 30 : return wbox;
324 15 : }
325 :
326 1061 : Quantity AnnRegion::_lengthToAngle(
327 : const Quantity& quantity, const uInt pixelAxis
328 : ) const {
329 1061 : if(quantity.getUnit() == "pix") {
330 93 : return getCsys().toWorldLength(quantity.getValue(), pixelAxis);
331 : }
332 968 : else if (! quantity.isConform("rad")) {
333 0 : throw AipsError (
334 0 : "Quantity " + String::toString(quantity)
335 0 : + " is not an angular measure nor is it in pixel units."
336 0 : );
337 : }
338 : else {
339 968 : return quantity;
340 : }
341 : }
342 :
343 0 : void AnnRegion::_printPrefix(ostream& os) const {
344 0 : if (isAnnotationOnly()) {
345 0 : os << "ann ";
346 : }
347 0 : else if (isDifference()) {
348 0 : os << "- ";
349 : }
350 0 : }
351 :
352 : }
353 :
|
