Line data Source code
1 : //# ComponentImager.cc: this defines ComponentImager which modifies images by ComponentLists
2 : //# Copyright (C) 1999,2000,2001,2002,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: casa-feedback@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: ComponentImager.cc 18855 2005-07-21 08:03:40Z nkilleen $
27 :
28 : #include <imageanalysis/ImageAnalysis/ComponentImager.h>
29 :
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Cube.h>
32 : #include <casacore/casa/Arrays/Vector.h>
33 : #include <casacore/casa/Containers/Block.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 : #include <casacore/casa/Arrays/IPosition.h>
36 : #include <casacore/casa/BasicMath/Math.h>
37 : #include <casacore/casa/BasicSL/Constants.h>
38 : #include <casacore/measures/Measures/MDirection.h>
39 : #include <casacore/measures/Measures/MFrequency.h>
40 : #include <casacore/measures/Measures/MeasRef.h>
41 : #include <casacore/measures/Measures/Stokes.h>
42 : #include <casacore/casa/Logging/LogIO.h>
43 : #include <casacore/casa/Logging/LogOrigin.h>
44 : #include <casacore/casa/Quanta/Unit.h>
45 : #include <casacore/casa/Quanta/UnitMap.h>
46 : #include <casacore/casa/Quanta/UnitVal.h>
47 : #include <casacore/casa/Quanta/MVAngle.h>
48 : #include <casacore/casa/Quanta/MVDirection.h>
49 : #include <casacore/casa/Quanta/MVFrequency.h>
50 : #include <casacore/casa/Quanta/Quantum.h>
51 : #include <casacore/casa/Quanta/QMath.h>
52 : #include <casacore/casa/Utilities/Assert.h>
53 : #include <casacore/casa/BasicSL/String.h>
54 : #include <casacore/casa/Arrays/ArrayLogical.h>
55 :
56 : #include <components/ComponentModels/ComponentList.h>
57 : #include <components/ComponentModels/SpectralModel.h>
58 : #include <components/ComponentModels/Flux.h>
59 : #include <components/ComponentModels/GaussianShape.h>
60 : #include <components/ComponentModels/PointShape.h>
61 : #include <components/ComponentModels/C11Timer.h>
62 : #include <components/ComponentModels/SkyCompRep.h>
63 :
64 : #include <casacore/coordinates/Coordinates/Coordinate.h>
65 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
66 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
67 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
68 : #include <casacore/images/Images/ImageInterface.h>
69 : #include <casacore/images/Images/ImageInfo.h>
70 : #include <casacore/lattices/Lattices/LatticeIterator.h>
71 : #include <casacore/lattices/Lattices/LatticeStepper.h>
72 :
73 : #include <iostream>
74 :
75 : #include <iomanip>
76 :
77 : using namespace casacore;
78 : namespace casa {
79 :
80 19 : ComponentImager::ComponentImager(
81 : const SPIIF image, const Record *const ®ion, const String& mask
82 19 : ) : ImageTask<Float>(image, region, mask, "", false),
83 19 : _image(image) {
84 19 : _construct();
85 19 : }
86 :
87 19 : ComponentImager::~ComponentImager() {}
88 :
89 19 : void ComponentImager::modify(Bool verbose) {
90 19 : *this->_getLog() << LogOrigin(getClass(), __func__);
91 19 : int nelem = _list.nelements();
92 19 : Vector<SkyComponent> mod(nelem);
93 102 : for (int i = 0; i < nelem; ++i) {
94 83 : mod[i] = _list.component(i);
95 : }
96 19 : const auto n = mod.size();
97 19 : ThrowIf(n == 0, "There are no components in the model componentlist");
98 : auto subImage = SubImageFactory<Float>::createSubImageRW(
99 19 : *_image, *this->_getRegion(), this->_getMask(),
100 57 : (verbose ? this->_getLog().get() : nullptr),
101 19 : AxesSpecifier(), this->_getStretch()
102 57 : );
103 : // Allow for subtraction/addition
104 19 : ComponentList cl;
105 102 : for (uInt i = 0; i < n; ++i) {
106 83 : SkyComponent sky = mod(i);
107 83 : if (_subtract) {
108 0 : sky.flux().scaleValue(-1.0);
109 : }
110 83 : cl.add(sky);
111 83 : }
112 : ComponentListImage cli(
113 19 : cl, subImage->coordinates(), subImage->shape(), False
114 19 : );
115 19 : Lattice<Float>* x = &cli;
116 19 : LatticeExpr<Float> expr;
117 19 : const auto& imageUnitName = subImage->units().getName();
118 19 : static const Unit solidAngleUnit("arcsec.arcsec");
119 19 : if (imageUnitName.contains("pixel")) {
120 : // nothing needs to be done
121 : }
122 4 : else if (imageUnitName.contains("beam")) {
123 4 : const auto& imageInfo = subImage->imageInfo();
124 4 : const auto& csys = subImage->coordinates();
125 4 : const auto& dc = csys.directionCoordinate();
126 4 : if (imageInfo.hasSingleBeam()) {
127 4 : auto f = imageInfo.getBeamAreaInPixels(0, 0, dc);
128 4 : expr = cli * f;
129 4 : x = &expr;
130 : }
131 0 : else if (imageInfo.hasMultipleBeams()) {
132 0 : auto nchan = imageInfo.nChannels();
133 0 : auto nstokes = imageInfo.nStokes();
134 0 : auto freqAxis = csys.spectralAxisNumber(False);
135 0 : auto stokesAxis = csys.polarizationAxisNumber(False);
136 0 : auto n = subImage->ndim();
137 0 : IPosition latShape(n, 1);
138 0 : auto hasFreq = freqAxis >= 0;
139 0 : auto hasStokes = stokesAxis >= 0;
140 0 : if (hasFreq) {
141 0 : latShape[freqAxis] = nchan;
142 : }
143 0 : if (hasStokes) {
144 0 : latShape[stokesAxis] = nstokes;
145 : }
146 0 : ArrayLattice<Float> beamAreas(latShape);
147 0 : IPosition pos(n, 0);
148 0 : for (uInt c=0; c<nchan; ++c) {
149 0 : if (hasFreq) {
150 0 : pos[freqAxis] = c;
151 : }
152 0 : for (uInt s=0; s<nstokes; ++s) {
153 0 : if (hasStokes) {
154 0 : pos[stokesAxis] = s;
155 : }
156 0 : beamAreas.putAt(imageInfo.getBeamAreaInPixels(c, s, dc), pos);
157 : }
158 : }
159 0 : expr = cli * beamAreas;
160 0 : x = &expr;
161 0 : }
162 : else {
163 0 : *_getLog() << LogIO::WARN
164 : << "No beam defined even though the image units contain a beam. "
165 0 : << "Will assume the beam area is one pixel" << LogIO::POST;
166 : }
167 : }
168 : else {
169 0 : *_getLog() << LogIO::WARN
170 : << "Image units [" << imageUnitName
171 : << "] are not dimensionally equivalent to Jy/pixel or Jy/beam. "
172 0 : << "Assuming they are equivalent Jy/pixel" << LogIO::POST;
173 : }
174 19 : *subImage += *x;
175 19 : }
176 :
177 0 : void ComponentImager::project(ImageInterface<Float>& image, const ComponentList& list) {
178 0 : const auto& coords = image.coordinates();
179 0 : const auto imageShape = image.shape();
180 0 : LogIO os(LogOrigin("ComponentImager", __func__));
181 :
182 : // I currently REQUIRE that:
183 : // * The list has at least one element.
184 : // * The image has at least one pixel.
185 : // * The image has one direction coordinate (only).
186 : // * The direction coordinate has two pixel and two world axes.
187 : // * Polarization and frequency coordinates are optional, however at most one
188 : // each of these coordinates can exist.
189 : // * If there is a Stokes axis it can only contain Stokes::I,Q,U,V pols.
190 : // * No other coordinate types, like LinearCoordinate, are used.
191 0 : ThrowIf(
192 : ! coords.hasDirectionCoordinate(), "Image does not have a direction coordinate"
193 : );
194 : uInt longAxis, latAxis;
195 : {
196 0 : const Vector<Int> dirAxes = coords.directionAxesNumbers();
197 0 : longAxis = dirAxes(0);
198 0 : latAxis = dirAxes(1);
199 0 : }
200 0 : DirectionCoordinate dirCoord = coords.directionCoordinate();
201 0 : dirCoord.setWorldAxisUnits(Vector<String>(2, "rad"));
202 :
203 : // Make sure get conversion frame, not just the native one
204 :
205 : MDirection::Types dirFrame;
206 0 : dirCoord.getReferenceConversion(dirFrame);
207 0 : const MeasRef<MDirection> dirRef(dirFrame);
208 :
209 0 : MVAngle pixelLatSize, pixelLongSize;
210 : {
211 0 : const Vector<Double> inc = dirCoord.increment();
212 0 : pixelLongSize = MVAngle(abs(inc[0]));
213 0 : pixelLatSize = MVAngle(abs(inc[1]));
214 : ;
215 0 : }
216 :
217 : // Check if there is a Stokes Axes and if so which polarizations. Otherwise
218 : // only grid the I polarisation.
219 : // because the code that puts pixel values needs at least one stokes for this to work.
220 0 : uInt nStokes = 1;
221 0 : if (coords.hasPolarizationCoordinate()) {
222 0 : StokesCoordinate stCoord = coords.stokesCoordinate();
223 0 : Vector<Int> types = stCoord.stokes();
224 0 : nStokes = types.nelements();
225 0 : for (uInt p = 0; p < nStokes; p++) {
226 0 : Stokes::StokesTypes type = Stokes::type(types[p]);
227 0 : ThrowIf(
228 : type != Stokes::I && type != Stokes::Q
229 : && type != Stokes::U && type != Stokes::V,
230 : "Unsupported stokes type " + Stokes::name(type)
231 : );
232 : }
233 0 : }
234 :
235 : // Check if there is a frequency axis and if so get the all the frequencies
236 : // as a Vector<MVFrequency>. Otherwise assume the reference frequency is the
237 : // same as the reference frequency of the first component in the list.
238 :
239 0 : MeasRef<MFrequency> freqRef;
240 0 : uInt nFreqs = 1;
241 0 : Int freqAxis = -1;
242 0 : Vector<MVFrequency> freqValues(nFreqs);
243 0 : if (coords.hasSpectralAxis()) {
244 0 : freqAxis = coords.spectralAxisNumber(false);
245 0 : nFreqs = (uInt)imageShape[freqAxis];
246 0 : freqValues.resize(nFreqs);
247 0 : SpectralCoordinate specCoord = coords.spectralCoordinate();
248 0 : specCoord.setWorldAxisUnits(Vector<String>(1, "Hz"));
249 :
250 : // Create Frequency MeasFrame; this will enable conversions between
251 : // spectral frames (e.g. the CS frame might be TOPO and the CL
252 : // frame LSRK)
253 :
254 : MFrequency::Types specConv;
255 0 : MEpoch epochConv;
256 0 : MPosition posConv;
257 0 : MDirection dirConv;
258 0 : specCoord.getReferenceConversion(specConv, epochConv, posConv, dirConv);
259 0 : MeasFrame measFrame(epochConv, posConv, dirConv);
260 0 : freqRef = MeasRef<MFrequency>(specConv, measFrame);
261 : Double thisFreq;
262 0 : for (uInt f = 0; f < nFreqs; f++) {
263 : // Includes any frame conversion
264 0 : ThrowIf (
265 : !specCoord.toWorld(thisFreq, f),
266 : "cannot convert a frequency value"
267 : );
268 0 : freqValues(f) = MVFrequency(thisFreq);
269 : }
270 0 : }
271 : else {
272 : const MFrequency& defaultFreq =
273 0 : list.component(0).spectrum().refFrequency();
274 0 : freqRef = defaultFreq.getRef();
275 0 : freqValues(0) = defaultFreq.getValue();
276 : }
277 :
278 : // Find out what the units are. Currently allowed units are anything
279 : // dimensionally equivalent to Jy/pixel or Jy/beam. If the former then the
280 : // pixel size at the center of the image is assumed to hold throughout the
281 : // image. If the latter then the beam is fished out of the header and a
282 : // 'beam' unit defined. If the units are not defined or are not one of the
283 : // above they are assumed to be Jy/pixel and a warning message is sent to the
284 : // logger.
285 :
286 0 : Unit fluxUnits;
287 : {
288 0 : Unit imageUnit = image.units();
289 0 : const String& imageUnitName = imageUnit.getName();
290 0 : UnitMap::putUser(
291 0 : "pixel", UnitVal(pixelLatSize.radian() *
292 0 : pixelLongSize.radian(), "rad.rad")
293 : );
294 : // redefine is required to reset Unit Cache
295 0 : const Unit pixel("pixel");
296 0 : if (imageUnitName.contains("pixel")) {
297 : // Get the new definition of the imageUnit which uses the new
298 : // definition of the pixels unit.
299 0 : imageUnit = image.units();
300 0 : fluxUnits.setValue(imageUnit.getValue() * pixel.getValue());
301 0 : fluxUnits.setName(imageUnitName + String(".") + pixel.getName());
302 : }
303 0 : else if (imageUnitName.contains("beam")) {
304 0 : const ImageInfo imageInfo = image.imageInfo();
305 : // FIXME this needs to support multiple beams
306 0 : const GaussianBeam beam = imageInfo.restoringBeam();
307 0 : if (beam.isNull()) {
308 : os << LogIO::WARN
309 : << "No beam defined even though the image units contain a beam"
310 0 : << endl << "Assuming the beam is one pixel" << LogIO::POST;
311 0 : UnitMap::putUser("beam", pixel.getValue());
312 : }
313 : else {
314 0 : const Quantity beamArea = beam.getArea("sr");
315 0 : UnitMap::putUser(
316 0 : "beam", UnitVal(beamArea.getValue(),
317 0 : beamArea.getFullUnit().getName()))
318 : ;
319 0 : }
320 0 : const Unit beamUnit("beam");
321 : const UnitVal fudgeFactor(
322 0 : pixel.getValue().getFac()/
323 0 : beamUnit.getValue().getFac()
324 0 : );
325 :
326 : // Get the new definition of the imageUnit which uses the new
327 : // definition of the beam unit. The re-use of the Unit constructor
328 : // from the String forces the new Unit definitions to take effect
329 0 : imageUnit = Unit(image.units().getName());
330 0 : fluxUnits.setValue(imageUnit.getValue() *
331 0 : beamUnit.getValue() * fudgeFactor);
332 0 : fluxUnits.setName(imageUnitName + String(".") + beamUnit.getName());
333 0 : }
334 : // 20101013 the code above for Jy/pixel doesn't work, since Unit doesn't
335 : // understand that Jy/pixel.pixel == Jy.
336 0 : const Unit jy("Jy");
337 0 : os << "Adding components to image with units [" << fluxUnits.getName() << "]" << LogIO::POST;
338 0 : if (fluxUnits.getName()=="Jy/pixel.pixel") {
339 0 : fluxUnits=jy;
340 : }
341 0 : if (fluxUnits != jy) {
342 : os << LogIO::WARN
343 : << "Image units [" << fluxUnits.getName() << "] are not dimensionally equivalent to "
344 : << "Jy/pixel or Jy/beam " << endl
345 : << "Ignoring the specified units and proceeding on the assumption"
346 0 : << " they are Jy/pixel" << LogIO::POST;
347 0 : fluxUnits = jy;
348 : }
349 0 : }
350 :
351 : // Does the image have a writable mask ? Output pixel values are
352 : // only modified if the mask==T and the coordinate conversions
353 : // succeeded. The mask==F on output if the coordinate conversion
354 : // fails (usually means a pixel is outside of the valid CoordinateSystem)
355 : //
356 0 : Bool doMask = false;
357 0 : if (image.isMasked() && image.hasPixelMask()) {
358 0 : if (image.pixelMask().isWritable()) {
359 0 : doMask = true;
360 : }
361 : else {
362 : os << LogIO::WARN
363 0 : << "The image is masked, but it cannot be written to" << LogIO::POST;
364 : }
365 : }
366 :
367 0 : auto myList = &list;
368 0 : const auto naxis = imageShape.nelements();
369 0 : IPosition pixelPosition(naxis, 0);
370 0 : Int polAxis = coords.polarizationAxisNumber(false);
371 : auto modifiedList = _doPoints(
372 : image, list, longAxis, latAxis, fluxUnits, dirRef,
373 : pixelLatSize, pixelLongSize, freqValues, freqRef,
374 : freqAxis, polAxis, nStokes
375 0 : );
376 :
377 :
378 0 : if (modifiedList) {
379 0 : myList = modifiedList.get();
380 : }
381 :
382 :
383 : // Setup an iterator to step through the image in chunks that can fit into
384 : // memory. Go to a bit of effort to make the chunk size as large as
385 : // possible but still minimize the number of tiles in the cache.
386 0 : auto chunkShape = imageShape;
387 : {
388 0 : const IPosition tileShape = image.niceCursorShape(2048*2048);
389 0 : chunkShape(longAxis) = tileShape(longAxis);
390 0 : chunkShape(latAxis) = tileShape(latAxis);
391 0 : }
392 0 : auto pixelShape = imageShape;
393 0 : pixelShape(longAxis) = pixelShape(latAxis) = 1;
394 0 : LatticeStepper pixelStepper(imageShape, pixelShape, LatticeStepper::RESIZE);
395 0 : LatticeIterator<Float> chunkIter(image, chunkShape);
396 0 : const uInt nDirs = chunkShape(longAxis) * chunkShape(latAxis);
397 0 : Cube<Double> pixelVals(4, nDirs, nFreqs);
398 0 : Vector<MVDirection> dirVals(nDirs);
399 0 : Vector<Bool> coordIsGood(nDirs);
400 0 : Vector<Double> pixelDir(2);
401 : uInt d;
402 :
403 0 : auto doSample = myList->nelements() > 0;
404 0 : Lattice<Bool>* pixelMaskPtr = 0;
405 0 : if (doMask) {
406 0 : pixelMaskPtr = &image.pixelMask();
407 : }
408 0 : std::unique_ptr<Array<Bool> > maskPtr;
409 0 : for (chunkIter.reset(); !chunkIter.atEnd(); chunkIter++) {
410 : // Iterate through sky plane of cursor and do coordinate conversions
411 :
412 0 : const IPosition& blc = chunkIter.position();
413 0 : const IPosition& trc = chunkIter.endPosition();
414 0 : d = 0;
415 0 : pixelDir[1] = blc[latAxis];
416 0 : coordIsGood = true;
417 0 : auto endLat = trc[latAxis];
418 0 : while (pixelDir[1] <= endLat) {
419 0 : pixelDir[0] = blc[longAxis];
420 0 : auto endLong = trc[longAxis];
421 0 : while (pixelDir[0] <= endLong) {
422 0 : if (!dirCoord.toWorld(dirVals[d], pixelDir)) {
423 : // These pixels will be masked
424 0 : coordIsGood[d] = false;
425 : }
426 0 : ++d;
427 0 : ++pixelDir[0];
428 : }
429 0 : ++pixelDir[1];
430 : }
431 0 : if (doSample) {
432 : // Sample model, converting the values in the components
433 : // to the specified direction and spectral frames
434 0 : myList->sample(
435 : pixelVals, fluxUnits, dirVals, dirRef, pixelLatSize,
436 : pixelLongSize, freqValues, freqRef
437 : );
438 : }
439 : else {
440 0 : pixelVals = 0;
441 : }
442 : // Modify data by model for this chunk of data
443 0 : auto& imageChunk = chunkIter.rwCursor();
444 :
445 : // Get input mask values if available
446 0 : if (doMask) {
447 0 : maskPtr.reset(
448 : new Array<Bool>(
449 0 : image.getMaskSlice(chunkIter.position(),
450 0 : chunkIter.cursorShape(), false)
451 0 : )
452 : );
453 : }
454 0 : d = 0;
455 0 : pixelPosition[latAxis] = 0;
456 0 : coordIsGood = true;
457 :
458 0 : while (pixelPosition[latAxis] < chunkShape[latAxis]) {
459 0 : pixelPosition(longAxis) = 0;
460 0 : while (pixelPosition[longAxis] < chunkShape[longAxis]) {
461 0 : if (coordIsGood[d]) {
462 0 : for (uInt f = 0; f < nFreqs; f++) {
463 0 : if (freqAxis >= 0) {
464 0 : pixelPosition(freqAxis) = f;
465 : }
466 0 : for (uInt s = 0; s < nStokes; s++) {
467 0 : if (polAxis >= 0) {
468 0 : pixelPosition(polAxis) = s;
469 : }
470 0 : if (! doMask || (doMask && (*maskPtr)(pixelPosition))) {
471 0 : imageChunk(pixelPosition) += pixelVals(s, d, f);
472 : }
473 : }
474 : }
475 : }
476 0 : else if (doMask) {
477 0 : (*maskPtr)(pixelPosition) = false;
478 : }
479 0 : d++;
480 0 : pixelPosition(longAxis)++;
481 : }
482 0 : pixelPosition(latAxis)++;
483 : }
484 : // Update output mask in appropriate fashion
485 :
486 0 : if (doMask) {
487 0 : pixelMaskPtr->putSlice(*maskPtr, chunkIter.position());
488 : }
489 0 : }
490 0 : }
491 :
492 0 : std::unique_ptr<ComponentList> ComponentImager::_doPoints(
493 : ImageInterface<Float>& image, const ComponentList& list,
494 : int longAxis, int latAxis, const Unit& fluxUnits,
495 : const MeasRef<MDirection>& dirRef, const MVAngle& pixelLatSize,
496 : const MVAngle& pixelLongSize, const Vector<MVFrequency>& freqValues,
497 : const MeasRef<MFrequency>& freqRef, Int freqAxis, Int polAxis, uInt nStokes
498 : ) {
499 : // deal with point sources separately
500 0 : vector<Int> pointSourceIdx;
501 0 : auto n = list.nelements();
502 0 : Vector<Double> pixel;
503 0 : MVDirection imageWorld;
504 0 : auto nFreqs = freqValues.size();
505 0 : Cube<Double> values(4, 1, nFreqs);
506 0 : IPosition pixelPosition(image.ndim(), 0);
507 0 : const auto& dirCoord = image.coordinates().directionCoordinate();
508 0 : const auto imageShape = image.shape();
509 0 : std::unique_ptr<ComponentList> modifiedList;
510 0 : for (uInt i=0; i<n; ++i) {
511 0 : if (list.getShape(i)->type() == ComponentType::POINT) {
512 0 : auto dir = list.getRefDirection(i);
513 0 : dirCoord.toPixel(pixel, dir);
514 0 : pixelPosition[longAxis] = floor(pixel[0] + 0.5);
515 0 : pixelPosition[latAxis] = floor(pixel[1] + 0.5);
516 : // in case the source and the coordinate system have different
517 : // ref frames
518 0 : dirCoord.toWorld(imageWorld, pixel);
519 0 : const auto& point = list.component(i);
520 0 : values = 0;
521 0 : Bool foundPixel = false;
522 0 : if (
523 0 : pixelPosition[longAxis] >= 0 && pixelPosition[latAxis] >= 0
524 0 : && pixelPosition[longAxis] < imageShape[longAxis]
525 0 : && pixelPosition[latAxis] < imageShape[latAxis]
526 : ) {
527 0 : point.sample(
528 0 : values, fluxUnits, Vector<MVDirection>(1, imageWorld),
529 : dirRef, pixelLatSize, pixelLongSize, freqValues, freqRef
530 : );
531 0 : foundPixel = anyNE(values, 0.0);
532 : }
533 0 : if (! foundPixel) {
534 : // look for the pixel in a 3x3 square around the target pixel
535 0 : auto targetPixel = pixelPosition;
536 0 : for (
537 0 : pixelPosition[longAxis]=targetPixel[longAxis]-1;
538 0 : pixelPosition[longAxis]<=targetPixel[longAxis]+1; ++pixelPosition[longAxis]
539 : ) {
540 0 : for (
541 0 : pixelPosition[latAxis]=targetPixel[latAxis]-1;
542 0 : pixelPosition[latAxis]<=targetPixel[latAxis]+1; ++pixelPosition[latAxis]
543 : ) {
544 0 : if (
545 0 : (pixelPosition[longAxis] != targetPixel[longAxis]
546 0 : || pixelPosition[latAxis] != targetPixel[latAxis])
547 0 : && pixelPosition[longAxis] >= 0 && pixelPosition[latAxis] >= 0
548 0 : && pixelPosition[longAxis] < imageShape[longAxis]
549 0 : && pixelPosition[latAxis] < imageShape[latAxis]
550 : ) {
551 0 : dirCoord.toWorld(imageWorld, pixel);
552 0 : point.sample(
553 0 : values, fluxUnits, Vector<MVDirection>(1, imageWorld),
554 : dirRef, pixelLatSize, pixelLongSize, freqValues, freqRef
555 : );
556 0 : foundPixel = anyNE(values, 0.0);
557 0 : if (foundPixel) {
558 0 : break;
559 : }
560 : }
561 0 : if (foundPixel) {
562 0 : break;
563 : }
564 : }
565 : }
566 0 : }
567 0 : if (foundPixel) {
568 0 : pointSourceIdx.push_back(i);
569 0 : for (uInt f = 0; f < nFreqs; f++) {
570 0 : if (freqAxis >= 0) {
571 0 : pixelPosition[freqAxis] = f;
572 : }
573 0 : for (uInt s = 0; s < nStokes; s++) {
574 0 : if (polAxis >= 0) {
575 0 : pixelPosition[polAxis] = s;
576 : }
577 0 : image.putAt(image.getAt(pixelPosition) + values(s, 0, f), pixelPosition);
578 : }
579 : }
580 : }
581 0 : }
582 : }
583 0 : if (! pointSourceIdx.empty()) {
584 0 : modifiedList.reset(new ComponentList(list));
585 0 : modifiedList->remove(Vector<Int>(pointSourceIdx));
586 : }
587 0 : return modifiedList;
588 0 : }
589 :
590 : }
|
