Line data Source code
1 : //# ComponentListImage.cc: defines the PagedImage class
2 : //# Copyright (C) 1994,1995,1996,1997,1998,1999,2000,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <imageanalysis/Images/ComponentListImage.h>
29 :
30 : #include <casacore/casa/OS/Path.h>
31 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
32 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
33 : #include <casacore/images/Images/ImageOpener.h>
34 : #include <casacore/images/Regions/RegionHandlerTable.h>
35 : #include <components/ComponentModels/ComponentShape.h>
36 : #include <components/ComponentModels/SpectralModel.h>
37 : #include <casacore/casa/Quanta/UnitMap.h>
38 : #ifdef _OPENMP
39 : #include <omp.h>
40 : #endif
41 : #include <set>
42 :
43 : using namespace casacore;
44 :
45 : namespace casa {
46 :
47 : const String ComponentListImage::IMAGE_TYPE = "ComponentListImage";
48 :
49 0 : ComponentListImage::ComponentListImage(
50 : const ComponentList& compList, const CoordinateSystem& csys,
51 : const IPosition& shape, const String& tableName, Bool doCache
52 0 : ) : ImageInterface<Float>(RegionHandlerTable(getTable, this)),
53 : _cl(compList.copy()), _modifiedCL(compList.copy()),
54 0 : _cache(doCache) {
55 0 : ThrowIf(tableName.empty(), "Table name cannot be empty");
56 0 : _cl.rename(Path(tableName));
57 : // resizing must precede setting of coordinate system
58 0 : _resize(shape);
59 0 : setCoordinateInfo(csys);
60 0 : setUnits("Jy/pixel");
61 0 : }
62 :
63 0 : ComponentListImage::ComponentListImage(
64 : const ComponentList& compList, const CoordinateSystem& csys,
65 : const IPosition& shape, Bool doCache
66 0 : ) : ImageInterface<Float>(), _cl(compList.copy()),
67 0 : _modifiedCL(compList.copy()), _cache(doCache) {
68 : // resizing must precede setting of coordinate system
69 0 : _resize(shape);
70 0 : setCoordinateInfo(csys);
71 0 : setUnits("Jy/pixel");
72 0 : }
73 :
74 0 : ComponentListImage::ComponentListImage(
75 : const String& filename, Bool doCache, MaskSpecifier maskSpec
76 0 : ) : ImageInterface<Float>(RegionHandlerTable(getTable, this)), _cl(filename, False, False),
77 : _modifiedCL(_cl.copy()),
78 0 : _cache(doCache) {
79 0 : _openLogTable();
80 0 : _restoreAll(_cl.getTable().keywordSet());
81 0 : _applyMaskSpecifier(maskSpec);
82 0 : }
83 :
84 0 : ComponentListImage::ComponentListImage(const ComponentListImage& image)
85 0 : : ImageInterface<Float>(image), _cl(image._cl),
86 0 : _modifiedCL(image._modifiedCL), _shape(image._shape),
87 0 : _latAxis(image._latAxis),
88 0 : _longAxis(image._longAxis), _freqAxis(image._freqAxis),
89 0 : _stokesAxis(image._stokesAxis), _pixelLongSize(image._pixelLongSize),
90 0 : _pixelLatSize(image._pixelLatSize), _dirRef(image._dirRef),
91 0 : _dirVals(image._dirVals), _freqRef(image._freqRef),
92 0 : _freqVals(image._freqVals), _ptSourcePixelVals(image._ptSourcePixelVals),
93 0 : _pixelToIQUV(image._pixelToIQUV), _cache(image._cache),
94 0 : _hasFreq(image._hasFreq), _hasStokes(image._hasStokes),
95 0 : _computedPtSources(image._computedPtSources),
96 0 : _defaultFreq(image._defaultFreq) {
97 0 : if (image._mask) {
98 0 : _mask.reset(new LatticeRegion(*image._mask));
99 : }
100 0 : if (image._valueCache) {
101 0 : _valueCache.reset(
102 0 : dynamic_cast<TempImage<Float>*>(image._valueCache->cloneII())
103 : );
104 : }
105 0 : }
106 :
107 0 : ComponentListImage::~ComponentListImage() {}
108 :
109 0 : ComponentListImage& ComponentListImage::operator=(const ComponentListImage& other) {
110 0 : if (this != &other) {
111 0 : ImageInterface<Float>::operator= (other);
112 0 : _cl = other._cl;
113 0 : _modifiedCL = other._modifiedCL;
114 0 : _shape = other._shape;
115 0 : if (other._mask) {
116 0 : _mask.reset(new LatticeRegion(*other._mask));
117 : }
118 0 : _latAxis = other._latAxis;
119 0 : _longAxis = other._longAxis;
120 0 : _freqAxis = other._freqAxis;
121 0 : _stokesAxis = other._stokesAxis;
122 0 : _dirRef = other._dirRef;
123 0 : _dirVals = other._dirVals;
124 0 : _freqRef = other._freqRef;;
125 0 : _freqVals = other._freqVals;
126 0 : _pixelLongSize = other._pixelLongSize;
127 0 : _pixelLatSize = other._pixelLatSize;
128 0 : _ptSourcePixelVals = other._ptSourcePixelVals;
129 0 : _pixelToIQUV = other._pixelToIQUV;
130 0 : _cache = other._cache;
131 0 : _hasFreq = other._hasFreq;
132 0 : _hasStokes = other._hasStokes;
133 0 : _computedPtSources = other._computedPtSources;
134 0 : _defaultFreq = other._defaultFreq;
135 0 : if (other._valueCache) {
136 0 : _valueCache.reset(
137 0 : dynamic_cast<TempImage<Float>*>(other._valueCache->cloneII())
138 : );
139 : }
140 : }
141 0 : return *this;
142 : }
143 0 : void ComponentListImage::apply(casacore::Float (*)(casacore::Float)) {
144 0 : ThrowCc("There is no general way to run " + String(__func__) + " on a ComponentList");
145 : }
146 :
147 0 : void ComponentListImage::apply(casacore::Float (*)(const casacore::Float&)) {
148 0 : ThrowCc("There is no general way to run " + String(__func__) + " on a ComponentList");
149 : }
150 :
151 0 : void ComponentListImage::apply(
152 : const casacore::Functional<casacore::Float, casacore::Float>&
153 : ) {
154 0 : ThrowCc("There is no general way to run " + String(__func__) + " on a ComponentList");
155 : }
156 :
157 0 : ImageInterface<Float>* ComponentListImage::cloneII() const {
158 0 : return new ComponentListImage(*this);
159 : }
160 :
161 0 : void ComponentListImage::copyData (const casacore::Lattice<casacore::Float>&) {
162 0 : ThrowCc("There is no general way to run " + String(__func__) + " on a ComponentList");
163 : }
164 :
165 0 : const ComponentList& ComponentListImage::componentList() const {
166 0 : return _cl;
167 : }
168 :
169 0 : Bool ComponentListImage::doGetMaskSlice(Array<Bool>& buffer, const Slicer& section) {
170 0 : if (_mask) {
171 0 : return _mask->doGetSlice(buffer, section);
172 : }
173 : else {
174 : // If no mask, base implementation returns a True mask.
175 0 : return MaskedLattice<Float>::doGetMaskSlice (buffer, section);
176 : }
177 : }
178 :
179 0 : Bool ComponentListImage::doGetSlice(
180 : Array<Float>& buffer, const Slicer& section
181 : ) {
182 0 : if (! *_computedPtSources) {
183 0 : _computePointSourcePixelValues();
184 : }
185 0 : const auto& chunkShape = section.length();
186 0 : if (_cache) {
187 : // if the values have already been computed and cached, just use them
188 0 : Array<Bool> mask(chunkShape);
189 0 : _valueCache->getMaskSlice(mask, section);
190 0 : if (allTrue(mask)) {
191 0 : _valueCache->getSlice(buffer, section);
192 0 : return True;
193 : }
194 : }
195 0 : auto secStart = section.start();
196 0 : const auto nDirs = chunkShape(_longAxis) * chunkShape(_latAxis);
197 0 : const auto nFreqs = _hasFreq ? chunkShape[_freqAxis] : 1;
198 0 : Cube<Double> pixelVals(4, nDirs, nFreqs);
199 0 : if (_modifiedCL.nelements() == 0) {
200 0 : pixelVals = 0;
201 : }
202 : else {
203 0 : Vector<MVDirection> dirVals(nDirs);
204 0 : auto secEnd = section.end();
205 0 : auto endLong = secEnd[_longAxis];
206 0 : auto endLat = secEnd[_latAxis];
207 0 : const auto& dirCoord = coordinates().directionCoordinate();
208 0 : if (_cache) {
209 0 : _getDirValsDoCache(dirVals, secStart, endLong, endLat, dirCoord);
210 : }
211 : else {
212 0 : _getDirValsNoCache(dirVals, secStart, endLong, endLat, dirCoord);
213 : }
214 0 : Vector<MVFrequency> freqVals(nFreqs);
215 0 : if (_hasFreq) {
216 0 : const auto& specCoord = coordinates().spectralCoordinate();
217 0 : if (_cache) {
218 0 : _getFreqValsDoCache(freqVals, secStart, nFreqs, specCoord);
219 : }
220 : else {
221 0 : _getFreqValsNoCache(freqVals, secStart, nFreqs, specCoord);
222 : }
223 : }
224 : else {
225 0 : freqVals[0] = _defaultFreq.getValue();
226 : }
227 0 : _modifiedCL.sample(
228 0 : pixelVals, Unit("Jy"), dirVals, _dirRef, _pixelLatSize,
229 0 : _pixelLongSize, freqVals, _freqRef
230 : );
231 : }
232 0 : _fillBuffer(
233 : buffer, chunkShape, secStart, nFreqs, pixelVals
234 : );
235 0 : if (_cache) {
236 0 : _valueCache->putSlice(buffer, secStart);
237 0 : Array<Bool> trueChunk(chunkShape, True);
238 0 : _valueCache->pixelMask().putSlice(trueChunk, secStart);
239 : }
240 0 : return True;
241 : }
242 :
243 0 : void ComponentListImage::doPutSlice (
244 : const Array<Float>&, const IPosition&, const IPosition&
245 : ) {
246 0 : ThrowCc(
247 : "ComponentListImage::" + String(__func__)
248 : + ": pixel values cannot be modified in a ComponentListImage"
249 : );
250 : }
251 :
252 0 : Table& ComponentListImage::getTable (void* imagePtr, Bool writable) {
253 0 : ComponentListImage* im = static_cast<ComponentListImage*>(imagePtr);
254 0 : if (writable) {
255 0 : im->_reopenRW();
256 : }
257 0 : return im->_cl._getTable();
258 : }
259 :
260 0 : const LatticeRegion* ComponentListImage::getRegionPtr() const {
261 0 : return _mask.get();
262 : }
263 :
264 0 : Bool ComponentListImage::hasPixelMask() const {
265 0 : return Bool(_mask);
266 : }
267 :
268 0 : String ComponentListImage::imageType() const {
269 0 : return (_cl.getTable().isNull() ? "Temporary " : "Persistent ") + IMAGE_TYPE;
270 : }
271 :
272 0 : Bool ComponentListImage::isMasked() const {
273 0 : return Bool(_mask);
274 : }
275 :
276 0 : Bool ComponentListImage::isPaged() const {
277 0 : return ! _cl.getTable().isNull();
278 : }
279 :
280 0 : Bool ComponentListImage::isPersistent() const {
281 0 : return ! _cl.getTable().isNull();
282 : }
283 :
284 0 : Bool ComponentListImage::isWritable() const {
285 0 : return False;
286 : }
287 :
288 0 : String ComponentListImage::name(bool stripPath) const {
289 0 : const static String tempIndicator = "Temporary ComponentListImage";
290 0 : const auto& table = _cl.getTable();
291 0 : if (table.isNull()) {
292 0 : return tempIndicator;
293 : }
294 : else {
295 0 : Path path(table.tableName());
296 0 : if (!stripPath) {
297 0 : return path.absoluteName();
298 : }
299 0 : return path.baseName();
300 : }
301 : }
302 :
303 0 : bool ComponentListImage::ok() const {
304 0 : auto n = _shape.size();
305 0 : auto x = coordinates().nPixelAxes();
306 0 : return n > 1 && n < 5 && x > 1 && x < 5;
307 : }
308 :
309 0 : LatticeBase* ComponentListImage::openCLImage(
310 : const String& name, const MaskSpecifier&
311 : ) {
312 0 : return new ComponentListImage(name);
313 : }
314 :
315 0 : const Lattice<Bool>& ComponentListImage::pixelMask() const {
316 0 : ThrowIf(
317 : ! _mask, "ComponentListImage::" + String(__func__)
318 : + " - no mask attached"
319 : );
320 0 : return *_mask;
321 : }
322 :
323 0 : void ComponentListImage::registerOpenFunction() {
324 0 : ImageOpener::registerOpenImageFunction(
325 : ImageOpener::COMPLISTIMAGE, &openCLImage
326 : );
327 0 : }
328 :
329 0 : void ComponentListImage::removeRegion(
330 : const String& name, RegionHandler::GroupType type, Bool throwIfUnknown
331 : ) {
332 : // Remove the default mask if it is the region to be removed.
333 0 : if (name == getDefaultMask()) {
334 0 : setDefaultMask("");
335 : }
336 0 : ImageInterface<Float>::removeRegion(name, type, throwIfUnknown);
337 0 : }
338 :
339 0 : void ComponentListImage::resize(const TiledShape& newShape) {
340 0 : _resize(newShape);
341 0 : _deleteCache();
342 0 : if (_cache) {
343 0 : _initCache();
344 : }
345 0 : }
346 :
347 0 : void ComponentListImage::set(const casacore::Float&) {
348 0 : ThrowCc(
349 : "There is no general way to run " + String(__func__)
350 : + " on a ComponentList"
351 : );
352 : }
353 :
354 0 : void ComponentListImage::setCache(casacore::Bool doCache) {
355 0 : if (doCache == _cache) {
356 : // already set to what is wanted
357 0 : return;
358 : }
359 0 : _cache = doCache;
360 0 : _initCache();
361 : }
362 :
363 0 : Bool ComponentListImage::setCoordinateInfo (const CoordinateSystem& csys) {
364 0 : auto x = csys.nPixelAxes();
365 0 : ThrowIf(
366 : x != _shape.size(),
367 : "Coordinate system must have the same dimensionality as the shape"
368 : );
369 0 : ThrowIf(
370 : ! csys.hasDirectionCoordinate(),
371 : "coordinate system must have a direction coordinate"
372 : );
373 0 : auto polAxisNum = csys.polarizationAxisNumber(False);
374 0 : _hasStokes = polAxisNum > 0;
375 0 : ThrowIf(
376 : _hasStokes && _shape[polAxisNum] > 4,
377 : "Polarization axis can have no more than four pixels"
378 : );
379 0 : if (_hasStokes) {
380 0 : const auto& stCoord = csys.stokesCoordinate();
381 0 : auto idx = stCoord.stokes();
382 0 : for (auto stokesIndex: idx) {
383 0 : Stokes::StokesTypes type = Stokes::type(stokesIndex);
384 0 : ThrowIf(
385 : type != Stokes::I && type != Stokes::Q
386 : && type != Stokes::U && type != Stokes::V,
387 : "Unsupported stokes type " + Stokes::name(type)
388 : );
389 : }
390 : }
391 : // implementation copied from PagedImage and tweaked.
392 0 : auto& table = _cl._getTable();
393 0 : ThrowIf(
394 : ! table.isNull() && ! table.isWritable(),
395 : "Image is not writable, cannot save coordinate system"
396 : );
397 0 : ThrowIf(
398 : ! ImageInterface<Float>::setCoordinateInfo(csys),
399 : "Could not set coordinate system"
400 : );
401 0 : if (! table.isNull()) {
402 : // we've tested for writability above, so if here it's writable
403 0 : _reopenRW();
404 : // Update the coordinates
405 0 : if (table.keywordSet().isDefined("coords")) {
406 0 : table.rwKeywordSet().removeField("coords");
407 : }
408 0 : ThrowIf(
409 : ! csys.save(table.rwKeywordSet(), "coords"),
410 : "Error writing coordinate system to image"
411 : );
412 : }
413 0 : _cacheCoordinateInfo(csys);
414 0 : return True;
415 : }
416 :
417 0 : void ComponentListImage::setDefaultMask(const String& regionName) {
418 : // Use the new region as the image's mask.
419 0 : _applyMask(regionName);
420 : // Store the new default name.
421 0 : ImageInterface<Float>::setDefaultMask(regionName);
422 0 : }
423 :
424 0 : Bool ComponentListImage::setImageInfo(const ImageInfo& info) {
425 0 : ThrowIf(info.hasBeam(), "A ComponentListImage may not have a beam(s)");
426 : // Set imageinfo in base class.
427 0 : Bool ok = ImageInterface<Float>::setImageInfo(info);
428 0 : auto& tab = _cl._getTable();
429 0 : if (ok && ! tab.isNull()) {
430 : // Make persistent in table keywords.
431 0 : _reopenRW();
432 0 : if (tab.isWritable()) {
433 : // Delete existing one if there.
434 0 : if (tab.keywordSet().isDefined("imageinfo")) {
435 0 : tab.rwKeywordSet().removeField("imageinfo");
436 : }
437 : // Convert info to a record and save as keyword.
438 0 : TableRecord rec;
439 0 : String error;
440 0 : if (imageInfo().toRecord(error, rec)) {
441 0 : tab.rwKeywordSet().defineRecord("imageinfo", rec);
442 : }
443 : else {
444 : // Could not convert to record.
445 0 : LogIO os;
446 0 : os << LogIO::SEVERE << "Error saving ImageInfo in image " << name()
447 0 : << "; " << error << LogIO::POST;
448 0 : ok = False;
449 : }
450 : }
451 : else {
452 : // Table not writable.
453 0 : LogIO os;
454 : os << LogIO::SEVERE
455 0 : << "Image " << name() << " is not writable; not saving ImageInfo"
456 0 : << LogIO::POST;
457 : }
458 : }
459 0 : return ok;
460 : }
461 :
462 0 : Bool ComponentListImage::setMiscInfo(const RecordInterface& newInfo) {
463 0 : setMiscInfoMember(newInfo);
464 0 : auto& tab = _cl._getTable();
465 0 : if (! tab.isNull()) {
466 0 : _reopenRW();
467 0 : if (! tab.isWritable()) {
468 0 : return False;
469 : }
470 0 : if (tab.keywordSet().isDefined("miscinfo")) {
471 0 : tab.rwKeywordSet().removeField("miscinfo");
472 : }
473 0 : tab.rwKeywordSet().defineRecord("miscinfo", newInfo);
474 : }
475 0 : return True;
476 : }
477 :
478 0 : Bool ComponentListImage::setUnits(const Unit& newUnits) {
479 0 : ThrowIf(
480 : newUnits.getName() != "Jy/pixel", "units must be Jy/pixel"
481 : );
482 0 : setUnitMember (newUnits);
483 0 : if (! _cl.getTable().isNull()) {
484 0 : _reopenRW();
485 0 : auto& tab = _cl._getTable();
486 0 : if (! tab.isWritable()) {
487 0 : return False;
488 : }
489 0 : if (tab.keywordSet().isDefined("units")) {
490 0 : tab.rwKeywordSet().removeField("units");
491 : }
492 0 : tab.rwKeywordSet().define("units", newUnits.getName());
493 : }
494 0 : return True;
495 : }
496 :
497 0 : IPosition ComponentListImage::shape() const {
498 0 : return _shape;
499 : }
500 :
501 0 : void ComponentListImage::useMask(MaskSpecifier spec) {
502 0 : _applyMaskSpecifier(spec);
503 0 : }
504 :
505 0 : void ComponentListImage::handleMath(const Lattice<Float>&, int) {
506 0 : ThrowCc(
507 : "There is no general way to run " + String(__func__)
508 : + " on a ComponentList"
509 : );
510 : }
511 :
512 0 : void ComponentListImage::_applyMask(const String& maskName) {
513 : // No region if no mask name is given.
514 0 : if (maskName.empty()) {
515 0 : _mask.reset();
516 0 : return;
517 : }
518 : // Reconstruct the ImageRegion object.
519 : // Turn the region into lattice coordinates.
520 : std::unique_ptr<ImageRegion> regPtr(
521 0 : getImageRegionPtr(maskName, RegionHandler::Masks)
522 0 : );
523 : std::shared_ptr<LatticeRegion> latReg(
524 0 : new LatticeRegion(regPtr->toLatticeRegion(coordinates(), shape()))
525 0 : );
526 : // The mask has to cover the entire image.
527 0 : ThrowIf(
528 : latReg->shape() != shape(),
529 : "region " + maskName + " does not cover the full image"
530 : );
531 0 : _mask = latReg;
532 : }
533 :
534 0 : void ComponentListImage::_applyMaskSpecifier (const MaskSpecifier& spec) {
535 : // Use default mask if told to do so.
536 : // If it does not exist, use no mask.
537 0 : auto name = spec.name();
538 0 : if (spec.useDefault()) {
539 0 : name = getDefaultMask();
540 0 : if (! hasRegion(name, RegionHandler::Masks)) {
541 0 : name = "";
542 : }
543 : }
544 0 : _applyMask(name);
545 0 : }
546 :
547 0 : void ComponentListImage::_cacheCoordinateInfo(const CoordinateSystem& csys) {
548 : // cache often used info
549 0 : const auto dirAxes = csys.directionAxesNumbers();
550 0 : _longAxis = dirAxes[0];
551 0 : _latAxis = dirAxes[1];
552 0 : const auto& dirCoord = csys.directionCoordinate();
553 : // use the conversion frame, not the native one
554 : MDirection::Types dirFrame;
555 0 : dirCoord.getReferenceConversion(dirFrame);
556 0 : _dirRef = MeasRef<MDirection>(dirFrame);
557 0 : const auto inc = dirCoord.increment();
558 0 : const auto units = dirCoord.worldAxisUnits();
559 0 : _pixelLongSize = MVAngle(Quantity(abs(inc[_longAxis]), units[_longAxis]));
560 0 : _pixelLatSize = MVAngle(Quantity(abs(inc[_latAxis]), units[_latAxis]));
561 0 : _freqAxis = csys.spectralAxisNumber(false);
562 0 : _hasFreq = _freqAxis > 0;
563 0 : if (_hasFreq) {
564 0 : const auto& specCoord = csys.spectralCoordinate();
565 :
566 : // Create Frequency MeasFrame; this will enable conversions between
567 : // spectral frames (e.g. the CS frame might be TOPO and the CL
568 : // frame LSRK)
569 :
570 : MFrequency::Types specConv;
571 0 : MEpoch epochConv;
572 0 : MPosition posConv;
573 0 : MDirection dirConv;
574 0 : specCoord.getReferenceConversion(
575 : specConv, epochConv, posConv, dirConv
576 : );
577 0 : MeasFrame measFrame(epochConv, posConv, dirConv);
578 0 : _freqRef = MeasRef<MFrequency>(specConv, measFrame);
579 : }
580 : else {
581 0 : _defaultFreq = _cl.component(0).spectrum().refFrequency();
582 0 : _freqRef = _defaultFreq.getRef();
583 : }
584 0 : _stokesAxis = csys.polarizationAxisNumber(False);
585 0 : _hasStokes = _stokesAxis > 0;
586 0 : if (_hasStokes) {
587 0 : auto nstokes = _shape[_stokesAxis];
588 0 : _pixelToIQUV.resize(nstokes);
589 0 : for (uInt s=0; s<nstokes; ++s) {
590 0 : auto mystokes = csys.stokesAtPixel(s);
591 0 : if (mystokes == "I") {
592 0 : _pixelToIQUV[s] = 0;
593 : }
594 0 : else if (mystokes == "Q") {
595 0 : _pixelToIQUV[s] = 1;
596 : }
597 0 : else if (mystokes == "U") {
598 0 : _pixelToIQUV[s] = 2;
599 : }
600 0 : else if (mystokes == "V") {
601 0 : _pixelToIQUV[s] = 3;
602 : }
603 : else {
604 0 : ThrowCc(
605 : "Unsupported stokes value " + mystokes + " at pixel "
606 : + String::toString(s)
607 : );
608 : }
609 : }
610 : }
611 : else {
612 0 : _pixelToIQUV.resize(1);
613 0 : _pixelToIQUV[0] = 0;
614 : }
615 0 : _deleteCache();
616 0 : if (_cache) {
617 0 : _initCache();
618 : }
619 0 : }
620 :
621 0 : void ComponentListImage::_computePointSourcePixelValues() {
622 0 : if (*_computedPtSources) {
623 0 : return;
624 : }
625 0 : _ptSourcePixelVals->clear();
626 0 : auto n = _cl.nelements();
627 0 : std::set<uInt> pointSourceIdx;
628 0 : for (uInt i=0; i<n; ++i) {
629 0 : if (_cl.getShape(i)->type() == ComponentType::POINT) {
630 0 : pointSourceIdx.insert(i);
631 : }
632 : }
633 0 : if (pointSourceIdx.empty()) {
634 0 : *_computedPtSources = True;
635 0 : return;
636 : }
637 0 : std::vector<uInt> foundPtSourceIdx;
638 0 : uInt nInside = 0;
639 0 : uInt nOutside = 0;
640 0 : _findPointSoures(foundPtSourceIdx, nInside, nOutside, pointSourceIdx);
641 0 : LogIO log(LogOrigin(IMAGE_TYPE, __func__));
642 0 : auto npts = pointSourceIdx.size();
643 0 : if (nInside > 0) {
644 : log << LogIO::NORMAL << "Found " << nInside << " of " << npts
645 : << " point sources located within the image and cached their "
646 0 : << "pixel coordinates." << LogIO::POST;
647 : }
648 0 : if (nOutside > 0) {
649 : log << LogIO::NORMAL << "Found " << nOutside << " of " << npts
650 : << " point sources to be located outside the image, so will ignore "
651 0 : << "those." << LogIO::POST;
652 : }
653 0 : if (! foundPtSourceIdx.empty()) {
654 0 : _modifiedCL.remove(Vector<Int>(foundPtSourceIdx));
655 : }
656 0 : *_computedPtSources = True;
657 : }
658 :
659 0 : void ComponentListImage::_deleteCache() {
660 0 : _freqVals.resize();
661 0 : _dirVals.resize();
662 0 : if (_valueCache) {
663 0 : _valueCache->resize(TiledShape());
664 : }
665 0 : _ptSourcePixelVals->clear();
666 0 : *_computedPtSources = False;
667 0 : }
668 :
669 0 : void ComponentListImage::_fillBuffer(
670 : Array<Float>& buffer, const IPosition& chunkShape,
671 : const IPosition& secStart, uInt nFreqs,
672 : const Cube<Double>& pixelVals
673 : ) const {
674 0 : if (buffer.size() != chunkShape) {
675 0 : buffer.resize(chunkShape, False);
676 : }
677 0 : const auto lookForPtSources = ! _ptSourcePixelVals->empty();
678 0 : const auto nLong = chunkShape[_longAxis];
679 0 : const auto nLat = chunkShape[_latAxis];
680 0 : const auto nStokes = _hasStokes ? chunkShape[_stokesAxis] : 1;
681 0 : const uInt ilatStart = secStart[_latAxis];
682 0 : const uInt startLong = secStart[_longAxis];
683 0 : const uInt startFreq = _hasFreq ? secStart[_freqAxis] : 0;
684 0 : const uInt startPol = _hasStokes ? secStart[_stokesAxis] : 0;
685 0 : const uInt ndim = chunkShape.size();
686 : #ifdef _OPENMP
687 0 : #pragma omp parallel for collapse(4)
688 : #endif
689 : for(uInt blat=0; blat<nLat; ++blat) {
690 : for(uInt blong=0; blong<nLong; ++blong) {
691 : for (uInt bfreq=0; bfreq<nFreqs; ++bfreq) {
692 : for (uInt bpol=0; bpol<nStokes; ++bpol) {
693 : uInt ilat = ilatStart + blat;
694 : IPosition posInArray(ndim);
695 : posInArray[_latAxis] = blat;
696 : std::pair<uInt, uInt> dirPos;
697 : uInt ilong = startLong + blong;
698 : dirPos.first = ilong;
699 : uInt d = blat * nLong + blong;
700 : posInArray[_longAxis] = blong;
701 : dirPos.second = ilat;
702 : uInt ifreq = startFreq + bfreq;
703 : if (_hasFreq) {
704 : posInArray[_freqAxis] = bfreq;
705 : }
706 : uInt ipol = startPol + bpol;
707 : if (_hasStokes) {
708 : posInArray[_stokesAxis] = bpol;
709 : }
710 : buffer(posInArray) = pixelVals(_pixelToIQUV[ipol], d, bfreq);
711 : if (lookForPtSources) {
712 : auto ptSourceContrib = _ptSourcePixelVals->find(dirPos);
713 : if (ptSourceContrib != _ptSourcePixelVals->end()) {
714 : buffer(posInArray) += ptSourceContrib->second(ifreq, ipol);
715 : }
716 : }
717 : }
718 : }
719 : }
720 : }
721 0 : }
722 :
723 0 : ComponentListImage::PtFound ComponentListImage::_findPixel(
724 : Cube<Double>& values, const IPosition& pixelPosition,
725 : const DirectionCoordinate& dirCoord, const SkyComponent& point,
726 : const Vector<MVFrequency>& freqValues
727 : ) const {
728 0 : MVDirection imageWorld;
729 0 : static const Unit jy("Jy");
730 0 : Vector<Double> pixel(2);
731 0 : pixel[0] = pixelPosition[0];
732 0 : pixel[1] = pixelPosition[1];
733 0 : dirCoord.toWorld(imageWorld, pixel);
734 0 : point.sample(
735 0 : values, jy, Vector<MVDirection>(1, imageWorld),
736 0 : _dirRef, _pixelLatSize, _pixelLongSize, freqValues, _freqRef
737 0 : );
738 0 : if (anyNE(values, 0.0)) {
739 0 : if (
740 0 : pixelPosition[_longAxis] >= 0 && pixelPosition[_latAxis] >= 0
741 0 : && pixelPosition[_longAxis] < _shape[_longAxis]
742 0 : && pixelPosition[_latAxis] < _shape[_latAxis]
743 : ) {
744 0 : return FOUND_INSIDE;
745 : }
746 : else {
747 0 : return FOUND_OUTSIDE;
748 : }
749 : }
750 0 : return NOT_FOUND;
751 : }
752 :
753 0 : ComponentListImage::PtFound ComponentListImage::_findPixelIn3x3Box(
754 : IPosition& pixelPosition, Cube<Double>& values,
755 : const DirectionCoordinate& dirCoord, const SkyComponent& point,
756 : const Vector<MVFrequency>& freqValues
757 : ) const {
758 : // look for the pixel in a 3x3 square around the target pixel
759 0 : auto targetPixel = pixelPosition;
760 0 : for (
761 0 : pixelPosition[_longAxis]=targetPixel[_longAxis]-1;
762 0 : pixelPosition[_longAxis]<=targetPixel[_longAxis]+1; ++pixelPosition[_longAxis]
763 : ) {
764 0 : for (
765 0 : pixelPosition[_latAxis]=targetPixel[_latAxis]-1;
766 0 : pixelPosition[_latAxis]<=targetPixel[_latAxis]+1; ++pixelPosition[_latAxis]
767 : ) {
768 : // we've already looked at the target pixel position before this method
769 : // was called and didn't find the point source there, so don't waste
770 : // time doing it again
771 0 : if (
772 : (
773 0 : pixelPosition[_longAxis] != targetPixel[_longAxis]
774 0 : || pixelPosition[_latAxis] != targetPixel[_latAxis]
775 : )
776 : ) {
777 0 : auto res = _findPixel(values, pixelPosition, dirCoord, point, freqValues);
778 0 : if (res != NOT_FOUND) {
779 0 : return res;
780 : }
781 : }
782 : }
783 : }
784 0 : return NOT_FOUND;
785 : }
786 :
787 0 : void ComponentListImage::_findPointSoures(
788 : std::vector<uInt>& foundPtSourceIdx, uInt& nInside, uInt& nOutside,
789 : const std::set<uInt>& pointSourceIdx
790 : ) {
791 0 : Vector<Double> pixel;
792 0 : auto nFreqs = _hasFreq ? _shape[_freqAxis] : 1;
793 0 : auto freqValues = _getAllFreqValues(nFreqs);
794 0 : Cube<Double> values(4, 1, nFreqs);
795 0 : IPosition pixelPosition(_shape.size(), 0);
796 0 : const auto& dirCoord = coordinates().directionCoordinate();
797 0 : uInt nStokes = _hasStokes ? _shape[_stokesAxis] : 1;
798 0 : std::pair<uInt, uInt> dirPos;
799 0 : for (const auto i: pointSourceIdx) {
800 0 : dirCoord.toPixel(pixel, _cl.getRefDirection(i));
801 0 : pixelPosition[_longAxis] = floor(pixel[0] + 0.5);
802 0 : pixelPosition[_latAxis] = floor(pixel[1] + 0.5);
803 0 : const auto& point = _cl.component(i);
804 0 : values = 0;
805 : // allow some slop at the image boundaries on the first pass
806 0 : if (
807 0 : pixelPosition[_longAxis] < -1
808 0 : || pixelPosition[_latAxis] < -1
809 0 : || pixelPosition[_longAxis] > _shape[_longAxis]
810 0 : || pixelPosition[_latAxis] > _shape[_latAxis]
811 : ) {
812 0 : ++nOutside;
813 0 : foundPtSourceIdx.push_back(i);
814 : }
815 : else {
816 0 : auto foundPixel = _findPixel(
817 : values, pixelPosition, dirCoord, point, freqValues
818 : );
819 0 : if (foundPixel == NOT_FOUND) {
820 0 : foundPixel = _findPixelIn3x3Box(
821 : pixelPosition, values, dirCoord, point, freqValues
822 : );
823 : }
824 0 : if (foundPixel == FOUND_INSIDE) {
825 0 : ++nInside;
826 0 : foundPtSourceIdx.push_back(i);
827 0 : dirPos.first = pixelPosition[_longAxis];
828 0 : dirPos.second = pixelPosition[_latAxis];
829 0 : auto myiter = _ptSourcePixelVals->find(dirPos);
830 0 : if (myiter == _ptSourcePixelVals->end()) {
831 0 : (*_ptSourcePixelVals)[dirPos]
832 0 : = Matrix<Float>(nFreqs, nStokes, 0);
833 : }
834 0 : for (uInt f = 0; f < nFreqs; ++f) {
835 0 : for (uInt s = 0; s < nStokes; ++s) {
836 0 : auto v = values(_pixelToIQUV[s], 0, f);
837 0 : (*_ptSourcePixelVals)[dirPos](f, s) += v;
838 : }
839 : }
840 : }
841 0 : else if (foundPixel == FOUND_OUTSIDE) {
842 0 : ++nOutside;
843 0 : foundPtSourceIdx.push_back(i);
844 : }
845 : }
846 : }
847 0 : }
848 :
849 0 : Vector<MVFrequency> ComponentListImage::_getAllFreqValues(uInt nFreqs) {
850 0 : Vector<MVFrequency> freqValues(nFreqs);
851 0 : if (! _hasFreq) {
852 0 : freqValues[0] = _defaultFreq.get("Hz");
853 0 : if (_cache) {
854 0 : _freqVals[0].reset(new MVFrequency(freqValues[0]));
855 : }
856 0 : return freqValues;
857 : }
858 : Double thisFreq;
859 0 : const auto& specCoord = coordinates().spectralCoordinate();
860 0 : auto freqUnit = specCoord.worldAxisUnits()[0];
861 0 : Quantity freq(0, freqUnit);
862 0 : for (Double pixelFreq=0; pixelFreq<nFreqs; ++pixelFreq) {
863 : // Includes any frame conversion
864 0 : ThrowIf (
865 : ! specCoord.toWorld(thisFreq, pixelFreq),
866 : "cannot convert a frequency value"
867 : );
868 0 : freq.setValue(thisFreq);
869 0 : freqValues[pixelFreq] = MVFrequency(freq);
870 0 : if (_cache) {
871 0 : _freqVals[pixelFreq].reset(new MVFrequency(freq));
872 : }
873 : }
874 0 : return freqValues;
875 : }
876 :
877 0 : void ComponentListImage::_getDirValsDoCache(
878 : Vector<MVDirection>& dirVals, const IPosition& secStart, uInt endLong,
879 : uInt endLat, const DirectionCoordinate& dirCoord
880 : ) {
881 0 : Vector<Double> pixelDir(2);
882 0 : IPosition iPixelDir(2);
883 0 : uInt d = 0;
884 0 : for(
885 0 : pixelDir[1]=secStart[_latAxis], iPixelDir[1]=secStart[_latAxis];
886 0 : pixelDir[1]<=endLat; ++pixelDir[1], ++iPixelDir[1]
887 : ) {
888 0 : for (
889 0 : pixelDir[0]=secStart[_longAxis], iPixelDir[0]=secStart[_longAxis];
890 0 : pixelDir[0]<=endLong; ++pixelDir[0], ++iPixelDir[0], ++d
891 : ) {
892 0 : auto dirVal = _dirVals(iPixelDir);
893 0 : if (dirVal) {
894 : // cached value exists, use it
895 0 : dirVals[d] = *dirVal;
896 : }
897 : else {
898 0 : std::shared_ptr<MVDirection> newDirVal(new MVDirection);
899 0 : if (! dirCoord.toWorld(*newDirVal, pixelDir)) {
900 0 : ostringstream oss;
901 0 : oss << "Unable to convert to world direction at pixel "
902 0 : << pixelDir;
903 0 : ThrowCc(oss.str());
904 : }
905 : // cache it
906 0 : _dirVals(iPixelDir) = newDirVal;
907 0 : dirVals[d] = *newDirVal;
908 : }
909 : }
910 : }
911 0 : }
912 :
913 0 : void ComponentListImage::_getDirValsNoCache(
914 : Vector<MVDirection>& dirVals, const IPosition& secStart, uInt endLong,
915 : uInt endLat, const DirectionCoordinate& dirCoord
916 : ) const {
917 0 : Vector<Double> pixelDir(2);
918 0 : IPosition iPixelDir(2);
919 0 : uInt d = 0;
920 0 : for(
921 0 : pixelDir[1]=secStart[_latAxis], iPixelDir[1]=secStart[_latAxis];
922 0 : pixelDir[1]<=endLat; ++pixelDir[1], ++iPixelDir[1]
923 : ) {
924 0 : for (
925 0 : pixelDir[0]=secStart[_longAxis], iPixelDir[0]=secStart[_longAxis];
926 0 : pixelDir[0]<=endLong; ++pixelDir[0], ++iPixelDir[0], ++d
927 : ) {
928 0 : if (! dirCoord.toWorld(dirVals[d], pixelDir)) {
929 0 : ostringstream oss;
930 0 : oss << "Unable to convert to world direction at pixel "
931 0 : << pixelDir;
932 0 : ThrowCc(oss.str());
933 : }
934 : }
935 : }
936 0 : }
937 :
938 0 : void ComponentListImage::_getFreqValsDoCache(
939 : Vector<MVFrequency>& freqVals, const IPosition& secStart,
940 : uInt nFreqs, const SpectralCoordinate& specCoord
941 : ) {
942 0 : uInt f = 0;
943 : Double thisFreq;
944 0 : auto freqUnit = specCoord.worldAxisUnits()[0];
945 0 : Quantity freq(0, freqUnit);
946 0 : for (Double pixelFreq=secStart[_freqAxis]; f<nFreqs; ++pixelFreq, ++f) {
947 0 : auto freqVal = _freqVals[pixelFreq];
948 0 : if (freqVal) {
949 : // already cached
950 0 : freqVals[f] = *freqVal;
951 : }
952 : else {
953 : // Includes any frame conversion
954 0 : ThrowIf (
955 : ! specCoord.toWorld(thisFreq, pixelFreq),
956 : "cannot convert a frequency value"
957 : );
958 0 : freq.setValue(thisFreq);
959 0 : freqVals[f] = freq;
960 0 : _freqVals[pixelFreq].reset(new MVFrequency(freq));
961 : }
962 : }
963 0 : }
964 :
965 0 : void ComponentListImage::_getFreqValsNoCache(
966 : Vector<MVFrequency>& freqVals, const IPosition& secStart,
967 : uInt nFreqs, const SpectralCoordinate& specCoord
968 : ) const {
969 0 : uInt f = 0;
970 : Double thisFreq;
971 0 : auto freqUnit = specCoord.worldAxisUnits()[0];
972 0 : Quantity freq(0, freqUnit);
973 0 : for (Double pixelFreq=secStart[_freqAxis]; f<nFreqs; ++pixelFreq, ++f) {
974 : // Includes any frame conversion
975 0 : ThrowIf (
976 : ! specCoord.toWorld(thisFreq, pixelFreq),
977 : "cannot convert a frequency value"
978 : );
979 0 : freq.setValue(thisFreq);
980 0 : freqVals[f] = freq;
981 : }
982 0 : }
983 :
984 0 : void ComponentListImage::_initCache() {
985 0 : if (_cache) {
986 0 : auto nLong = _shape[_longAxis];
987 0 : auto nLat = _shape[_latAxis];
988 0 : _dirVals = Matrix<std::shared_ptr<MVDirection>>(nLong, nLat);
989 0 : _dirVals.set(std::shared_ptr<MVDirection>(nullptr));
990 0 : auto nFreqs = _hasFreq ? _shape[_freqAxis] : 1;
991 0 : _freqVals.resize(nFreqs);
992 0 : _freqVals.set(std::shared_ptr<MVFrequency>(nullptr));
993 0 : _valueCache.reset(new TempImage<Float>(_shape, coordinates()));
994 0 : _valueCache->attachMask(TempLattice<Bool>(_shape, False));
995 0 : *_computedPtSources = False;
996 0 : _ptSourcePixelVals->clear();
997 : }
998 : else {
999 0 : _deleteCache();
1000 : }
1001 0 : }
1002 :
1003 0 : void ComponentListImage::_openLogTable() {
1004 : // Open logtable as readonly if main table is not writable.
1005 0 : auto& tab = _cl._getTable();
1006 0 : setLogMember (LoggerHolder (name() + "/logtable", tab.isWritable()));
1007 : // Insert the keyword if possible and if it does not exist yet.
1008 0 : if (tab.isWritable() && ! tab.keywordSet().isDefined ("logtable")) {
1009 0 : tab.rwKeywordSet().defineTable("logtable", Table(name() + "/logtable"));
1010 : }
1011 0 : }
1012 :
1013 0 : void ComponentListImage::_reopenRW() {
1014 : // implementation copied from PagedImage and tweaked
1015 0 : auto& table = _cl._getTable();
1016 : //# Open for write if not done yet and if writable.
1017 0 : if (!table.isWritable()) {
1018 0 : table.reopenRW();
1019 : }
1020 0 : }
1021 :
1022 0 : void ComponentListImage::_resize(const TiledShape& newShape) {
1023 0 : auto shape = newShape.shape();
1024 0 : ThrowIf(
1025 : shape.size() < 2 || shape.size() > 4,
1026 : "ComponentListImages must have 2, 3, or 4 dimensions"
1027 : );
1028 0 : ThrowIf(anyLE(shape.asVector(), 0), "All shape elements must be positive");
1029 0 : if (! _shape.conform(shape)) {
1030 0 : _shape.resize(shape.size(), False);
1031 : }
1032 0 : _shape = shape;
1033 0 : auto& table = _cl._getTable();
1034 0 : if (! table.isNull()) {
1035 0 : _reopenRW();
1036 0 : if (table.isWritable()) {
1037 : // Update the shape
1038 0 : if (table.keywordSet().isDefined("shape")) {
1039 0 : table.rwKeywordSet().removeField("shape");
1040 : }
1041 0 : table.rwKeywordSet().define("shape", _shape.asVector());
1042 : }
1043 : else {
1044 0 : LogIO os;
1045 0 : os << LogIO::SEVERE << "Image " << name()
1046 : << " is not writable; not saving shape"
1047 0 : << LogIO::POST;
1048 : }
1049 : }
1050 0 : }
1051 :
1052 0 : void ComponentListImage::_restoreAll(const TableRecord& rec) {
1053 : // must do shape before coordinates
1054 0 : ThrowIf(! rec.isDefined("shape"), "shape is not present in " + name());
1055 0 : ThrowIf(
1056 : rec.type(rec.fieldNumber("shape")) != TpArrayInt,
1057 : "shape found in " + name() + " is not stored as an integer array"
1058 : );
1059 0 : _resize(IPosition(rec.asArrayInt("shape")));
1060 : // Restore the coordinates.
1061 : std::unique_ptr<CoordinateSystem> restoredCoords(
1062 : CoordinateSystem::restore(rec, "coords")
1063 0 : );
1064 0 : ThrowIf(! restoredCoords, "Error restoring coordinate system");
1065 0 : setCoordinateInfo(*restoredCoords);
1066 : // Restore the image info.
1067 0 : _restoreImageInfo(rec);
1068 : // Restore the units.
1069 0 : _restoreUnits(rec);
1070 : // Restore the miscinfo.
1071 0 : _restoreMiscInfo(rec);
1072 0 : }
1073 :
1074 0 : void ComponentListImage::_restoreImageInfo (const TableRecord& rec) {
1075 0 : if (rec.isDefined("imageinfo")) {
1076 0 : String error;
1077 0 : ImageInfo info;
1078 0 : Bool ok = info.fromRecord (error, rec.asRecord("imageinfo"));
1079 0 : if (ok) {
1080 0 : setImageInfoMember(info);
1081 : }
1082 : else {
1083 0 : LogIO os;
1084 0 : os << LogIO::WARN << "Failed to restore the ImageInfo in image " << name()
1085 0 : << "; " << error << LogIO::POST;
1086 : }
1087 : }
1088 0 : }
1089 :
1090 0 : void ComponentListImage::_restoreMiscInfo (const TableRecord& rec) {
1091 0 : if (rec.isDefined("miscinfo") && rec.dataType("miscinfo") == TpRecord) {
1092 0 : setMiscInfoMember (rec.asRecord ("miscinfo"));
1093 : }
1094 0 : }
1095 :
1096 0 : void ComponentListImage::_restoreUnits (const TableRecord& rec) {
1097 0 : Unit retval;
1098 0 : String unitName;
1099 0 : if (rec.isDefined("units")) {
1100 0 : if (rec.dataType("units") != TpString) {
1101 0 : LogIO os;
1102 0 : os << LogOrigin("PagedImage<T>", "units()", WHERE)
1103 : << "'units' keyword in image table is not a string! Units not restored."
1104 0 : << LogIO::SEVERE << LogIO::POST;
1105 : }
1106 : else {
1107 0 : rec.get("units", unitName);
1108 : }
1109 : }
1110 0 : if (! unitName.empty()) {
1111 : // OK, non-empty unit, see if it's valid, if not try some known things to
1112 : // make a valid unit out of it.
1113 0 : if (! UnitVal::check(unitName)) {
1114 : // Beam and Pixel are the most common undefined units
1115 0 : UnitMap::putUser("Pixel",UnitVal(1.0),"Pixel unit");
1116 0 : UnitMap::putUser("Beam",UnitVal(1.0),"Beam area");
1117 : }
1118 0 : if (! UnitVal::check(unitName)) {
1119 : // OK, maybe we need FITS
1120 0 : UnitMap::addFITS();
1121 : }
1122 0 : if (!UnitVal::check(unitName)) {
1123 : // I give up!
1124 0 : LogIO os;
1125 0 : UnitMap::putUser(unitName, UnitVal(1.0, UnitDim::Dnon), unitName);
1126 : os << LogIO::WARN << "FITS unit \"" << unitName
1127 : << "\" unknown to CASA - will treat it as non-dimensional."
1128 0 : << LogIO::POST;
1129 0 : retval.setName(unitName);
1130 0 : retval.setValue(UnitVal(1.0, UnitDim::Dnon));
1131 : }
1132 : else {
1133 0 : retval = Unit(unitName);
1134 : }
1135 : }
1136 0 : setUnitMember(retval);
1137 0 : }
1138 :
1139 : }
|
