Line data Source code
1 : //# ImageDecomposer.cc: defines the ImageDecomposer class
2 : //# Copyright (C) 1994,1995,1996,1997,1998,1999,2000,2001,2002
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: ImageDecomposer.tcc 20739 2009-09-29 01:15:15Z Malte.Marquarding $
27 :
28 : #include <imageanalysis/ImageAnalysis/ImageDecomposer.h>
29 :
30 : #include <casacore/casa/Arrays/Slicer.h>
31 : #include <casacore/casa/Arrays/Cube.h>
32 : #include <casacore/casa/IO/ArrayIO.h>
33 : #include <casacore/lattices/Lattices/TiledShape.h>
34 : #include <casacore/scimath/Fitting/FitGaussian.h>
35 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
36 : #include <casacore/lattices/LRegions/LatticeRegion.h>
37 : #include <casacore/lattices/LRegions/LCMask.h>
38 : #include <casacore/lattices/LEL/LatticeExpr.h>
39 : #include <casacore/lattices/LEL/LatticeExprNode.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <casacore/casa/BasicMath/Math.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 : #include <casacore/images/Images/TempImage.h>
44 : #include <casacore/images/Images/SubImage.h>
45 :
46 :
47 : namespace casa {
48 :
49 : template <class T>
50 0 : ImageDecomposer<T>::ImageDecomposer(const casacore::ImageInterface<T>& image)
51 0 : : itsImagePtr(image.cloneII()),
52 0 : itsMapPtr(0),
53 0 : itsShape(itsImagePtr->shape()),
54 0 : itsDim(itsShape.nelements()),
55 0 : itsNRegions(0),
56 0 : itsNComponents(0),
57 0 : itsDeblendIt(true),
58 0 : itsThresholdVal(0.1),
59 0 : itsNContour(11),
60 0 : itsMinRange(2),
61 0 : itsNAxis(2),
62 0 : itsFitIt(true),
63 0 : itsMaximumRMS(0.1),
64 0 : itsMaxRetries(-1),
65 0 : itsMaxIter(256),
66 0 : itsConvCriteria(0.001)
67 : {
68 0 : itsMapPtr = new casacore::TempLattice<casacore::Int>(casacore::TiledShape(itsShape), 1);
69 0 : if (!itsMapPtr) {
70 0 : delete itsImagePtr;
71 0 : throw(casacore::AipsError("Failed to create internal TempLattice"));
72 : }
73 : //
74 0 : itsMapPtr->set(0);
75 0 : }
76 :
77 :
78 : template <class T>
79 0 : ImageDecomposer<T>::ImageDecomposer(const ImageDecomposer<T>& other)
80 :
81 0 : : itsImagePtr(other.itsImagePtr->cloneII()),
82 0 : itsMapPtr(0),
83 0 : itsShape(other.itsShape),
84 0 : itsDim(other.itsDim),
85 0 : itsNRegions(0),
86 0 : itsNComponents(0)
87 : {
88 0 : itsMapPtr = new casacore::TempLattice<casacore::Int>(casacore::TiledShape(itsShape), 1);
89 0 : if (!itsMapPtr) {
90 0 : delete itsImagePtr;
91 0 : throw(casacore::AipsError("Failed to create internal TempLattice"));
92 : }
93 : //
94 0 : itsNRegions = other.itsNRegions;
95 0 : itsNComponents = other.itsNComponents;
96 0 : itsList = other.itsList.copy();
97 :
98 0 : copyOptions(other);
99 :
100 : // Copy values from other.itsMapPtr
101 :
102 0 : itsMapPtr->copyData(*(other.itsMapPtr));
103 0 : }
104 : /*
105 : template <class T>
106 : ImageDecomposer<T> &ImageDecomposer<T>::operator=(const ImageDecomposer<T> &other)
107 : {
108 : if (this != &other) {
109 : delete itsImagePtr;
110 : delete itsMapPtr;
111 :
112 : itsImagePtr = other.itsImagePtr->cloneII();
113 : itsShape = other.itsShape;
114 : itsDim = other.itsDim;
115 : itsNRegions = 0;
116 : itsNComponents = 0;
117 :
118 : itsMapPtr = new casacore::TempLattice<casacore::Int>(casacore::TiledShape(itsShape), 1);
119 : itsMapPtr->copyData(*(other.itsMapPtr));
120 : itsList = other.itsList.copy();
121 :
122 : copyOptions(other);
123 : }
124 : return *this;
125 : }
126 : */
127 :
128 : template <class T>
129 0 : ImageDecomposer<T>::~ImageDecomposer()
130 : {
131 0 : if (itsImagePtr) {
132 0 : delete itsImagePtr;
133 0 : itsImagePtr = 0;
134 : }
135 0 : if (itsMapPtr) {
136 0 : delete itsMapPtr;
137 0 : itsMapPtr = 0;
138 : }
139 0 : }
140 :
141 : // Basic set/get functions
142 : // -----------------------
143 :
144 : /*
145 : template <class T>
146 : void ImageDecomposer<T>::setImage(casacore::ImageInterface<T>& image)
147 : {
148 : if (itsImagePtr) {
149 : delete itsImagePtr;
150 : itsImagePtr = 0;
151 : }
152 : if (itsMapPtr) {
153 : delete itsMapPtr;
154 : itsMapPtr = 0;
155 : }
156 : //
157 : itsImagePtr = image.cloneII();
158 : itsShape.resize(0); //necessary for dimension change
159 : itsShape = itsImagePtr->shape();
160 : itsDim = itsShape.nelements();
161 : itsNRegions = 0;
162 : itsNComponents = 0;
163 : //
164 : itsMapPtr = new casacore::TempLattice<casacore::Int>(casacore::TiledShape(itsShape), 1);
165 : if (!itsMapPtr) {
166 : delete itsImagePtr;
167 : throw(casacore::AipsError("Failed to create internal TempLattice"));
168 : }
169 : itsMapPtr->set(0);
170 : }
171 : */
172 :
173 : template <class T>
174 0 : void ImageDecomposer<T>::setDeblend(casacore::Bool deblendIt)
175 : {
176 0 : itsDeblendIt = deblendIt;
177 0 : return;
178 : }
179 :
180 : template <class T>
181 0 : void ImageDecomposer<T>::setDeblendOptions(T thresholdVal,casacore::uInt nContour,
182 : casacore::Int minRange, casacore::Int nAxis)
183 : {
184 0 : itsThresholdVal = thresholdVal;
185 0 : itsNContour = nContour;
186 0 : itsMinRange = minRange;
187 0 : itsNAxis = nAxis;
188 0 : return;
189 : }
190 :
191 : template <class T>
192 0 : void ImageDecomposer<T>::setFit(casacore::Bool fitIt)
193 : {
194 0 : itsFitIt = fitIt;
195 0 : return;
196 : }
197 :
198 : template <class T>
199 : void ImageDecomposer<T>::setFitOptions(T maximumRMS, casacore::Int maxRetries,
200 : casacore::uInt maxIter, T convCriteria)
201 : {
202 : itsMaximumRMS = maximumRMS;
203 : itsMaxRetries = maxRetries;
204 : itsMaxIter = maxIter;
205 : itsConvCriteria = convCriteria;
206 : return;
207 : }
208 :
209 :
210 : template <class T>
211 0 : void ImageDecomposer<T>::copyOptions(const ImageDecomposer<T> &other)
212 : {
213 0 : itsDeblendIt = other.itsDeblendIt;
214 0 : itsThresholdVal = other.itsThresholdVal;
215 0 : itsNContour = other.itsNContour;
216 0 : itsMinRange = other.itsMinRange;
217 0 : itsNAxis = other.itsNAxis;
218 0 : itsFitIt = other.itsFitIt;
219 0 : itsMaximumRMS = other.itsMaximumRMS;
220 0 : itsMaxRetries = other.itsMaxRetries;
221 0 : itsMaxIter = other.itsMaxIter;
222 0 : itsConvCriteria = other.itsConvCriteria;
223 0 : return;
224 : }
225 :
226 :
227 : template <class T>
228 : int ImageDecomposer<T>::getCell(casacore::Int x, casacore::Int y) const
229 : {
230 : return itsMapPtr->getAt(casacore::IPosition(2,x,y));
231 : }
232 :
233 : template <class T>
234 : int ImageDecomposer<T>::getCell(casacore::Int x, casacore::Int y, casacore::Int z) const
235 : {
236 : return itsMapPtr->getAt(casacore::IPosition(3,x,y,z));
237 : }
238 : template <class T>
239 0 : int ImageDecomposer<T>::getCell(const casacore::IPosition& coord) const
240 : {
241 0 : return itsMapPtr->getAt(coord); //note: 3D casacore::IPosition works on 2D image
242 : }
243 :
244 :
245 : template <class T>
246 : void ImageDecomposer<T>::setCell(casacore::Int x, casacore::Int y, casacore::Int sval)
247 : {
248 : itsMapPtr->putAt(sval, casacore::IPosition(2,x,y));
249 : return;
250 : }
251 : template <class T>
252 : void ImageDecomposer<T>::setCell(casacore::Int x, casacore::Int y, casacore::Int z, casacore::Int sval)
253 : {
254 : itsMapPtr->putAt(sval, casacore::IPosition(3,x,y,z));
255 : return;
256 : }
257 : template <class T>
258 0 : void ImageDecomposer<T>::setCell(const casacore::IPosition& coord, casacore::Int sval)
259 : {
260 0 : itsMapPtr->putAt(sval, coord);
261 0 : return;
262 : }
263 :
264 : template <class T>
265 0 : casacore::IPosition ImageDecomposer<T>::shape() const
266 : {
267 0 : return itsShape;
268 : }
269 :
270 : template <class T>
271 0 : int ImageDecomposer<T>::shape(casacore::uInt axis) const
272 : {
273 0 : if (itsDim > axis) return itsShape(axis);
274 0 : return 1;
275 : }
276 :
277 : //The numRegions and numComponents functions are both getters for their
278 : //respective variables and flags telling whether the image has been derived
279 : //and/or decomposed.
280 :
281 : template <class T>
282 0 : casacore::uInt ImageDecomposer<T>::numRegions() const
283 : {
284 0 : return itsNRegions;
285 : }
286 :
287 : template <class T>
288 0 : casacore::uInt ImageDecomposer<T>::numComponents() const
289 : {
290 0 : return itsNComponents;
291 : }
292 :
293 : template <class T>
294 0 : bool ImageDecomposer<T>::isDerived() const
295 : {
296 0 : return itsNRegions>0;
297 : }
298 :
299 : template <class T>
300 0 : bool ImageDecomposer<T>::isDecomposed() const
301 : {
302 0 : return itsNComponents>0;
303 : }
304 :
305 : template <class T>
306 0 : T ImageDecomposer<T>::getImageVal(casacore::Int x, casacore::Int y) const
307 : {
308 0 : return getImageVal(casacore::IPosition(2,x,y));
309 : }
310 :
311 : template <class T>
312 0 : T ImageDecomposer<T>::getImageVal(casacore::Int x, casacore::Int y, casacore::Int z) const
313 : {
314 0 : return getImageVal(casacore::IPosition(3,x,y,z));
315 : }
316 :
317 : template <class T>
318 0 : T ImageDecomposer<T>::getImageVal(casacore::IPosition coord) const
319 : {
320 0 : return itsImagePtr->getAt(coord);
321 : }
322 :
323 :
324 : template <class T>
325 : int ImageDecomposer<T>::getContourVal(casacore::Int x, casacore::Int y, const casacore::Vector<T>& clevels) const
326 : {
327 : return getContourVal(casacore::IPosition(2,x,y), clevels);
328 : }
329 :
330 : template <class T>
331 : int ImageDecomposer<T>::getContourVal(casacore::Int x, casacore::Int y, casacore::Int z,
332 : const casacore::Vector<T>& clevels) const
333 : {
334 : return getContourVal(casacore::IPosition(3,x,y,z), clevels);
335 : }
336 :
337 : template <class T>
338 : int ImageDecomposer<T>::getContourVal(casacore::IPosition coord,
339 : const casacore::Vector<T>& clevels) const
340 : {
341 : T val = itsImagePtr->getAt(coord);
342 : for (casacore::uInt c = 0; c < clevels.nelements(); c++) {
343 : if (val < clevels(c)) return c - 1;
344 : }
345 : return clevels.nelements()-1;
346 : }
347 :
348 : template <class T>
349 0 : int ImageDecomposer<T>::getContourVal(T val, const casacore::Vector<T>& clevels) const
350 : {
351 0 : for (casacore::uInt c = 0; c < clevels.nelements(); c++) {
352 0 : if (val < clevels(c)) return c - 1;
353 : }
354 0 : return clevels.nelements()-1;
355 : }
356 :
357 : template <class T>
358 : casacore::Vector<T> ImageDecomposer<T>::autoContour(T mincon, T maxcon, T inc) const
359 : {
360 : if (inc == T(0)) {
361 : throw(casacore::AipsError("Vector<T> ImageDecomposer<T>::autocontour"
362 : "T mincon, T maxcon, T inc) - inc cannot be zero"));
363 : }
364 : if ((maxcon - mincon) * inc < 0) inc = -inc;
365 :
366 : casacore::Int c = 0;
367 : for (T cl = mincon; cl <= maxcon; cl += inc) c++;
368 :
369 : casacore::Vector<T> contours(c);
370 : c = 0;
371 : for (T cl = mincon; cl <= maxcon; cl += inc) {
372 : contours(c++) = cl;
373 : }
374 : return contours;
375 : }
376 :
377 : template <class T>
378 0 : casacore::Vector<T> ImageDecomposer<T>::autoContour(casacore::Int nContours, T minValue) const
379 : {
380 : // IMPR: a noise estimate to determine default value of lowest contour
381 : // would be useful.
382 :
383 0 : casacore::Vector<T> contours(nContours);
384 : T maxValue;
385 : //
386 0 : maxValue = findAreaGlobalMax(casacore::IPosition(itsDim,0), shape());
387 0 : maxValue -= (maxValue-minValue)/((nContours-1)*3);
388 : //cerr << "Autocontour: minvalue, maxvalue = " << minValue << ", " << maxValue << endl;
389 :
390 : // Make maximum contour ~1/3 contour increment less than max value of image
391 :
392 0 : for (casacore::Int i=0; i<nContours; i++) {
393 0 : contours(i) = minValue + (maxValue-minValue)*i/(nContours-1);
394 : }
395 : //
396 0 : return contours;
397 0 : }
398 :
399 : template <class T>
400 : casacore::Vector<T> ImageDecomposer<T>::autoContour(const casacore::Function1D<T>& fn,
401 : casacore::Int ncontours, T minvalue) const
402 : {
403 : // NOTE: This function has not been recently tested.
404 :
405 : casacore::Vector<T> contours(ncontours);
406 : T maxvalue;
407 : T calibzero, calibmax;
408 : //
409 : for (casacore::Int i=1; i<ncontours; i++) {
410 : if (fn(T(i-1))>fn(T(i))) {
411 : throw(casacore::AipsError("ImageDecomposer<T>::autoContour-"
412 : " fn must be nondecreasing in domain"));
413 : }
414 : }
415 : //
416 : maxvalue = findAreaGlobalMax(casacore::IPosition(itsDim,0), shape());
417 : maxvalue -= (maxvalue-minvalue)/((ncontours-1)*10); //much closer to top
418 : calibzero = minvalue - fn(T(0));
419 : calibmax = (maxvalue - calibzero) / fn(ncontours - 1);
420 : //
421 : for (casacore::Int i=0; i<ncontours; i++) {
422 : contours(i) = calibzero + calibmax*fn(i);
423 : }
424 : //
425 : return contours;
426 : }
427 :
428 : template <class T>
429 0 : casacore::Matrix<T> ImageDecomposer<T>::componentList() const
430 : {
431 : //IMPR: the pixel->world conversion shouldn't have to be done every time.
432 :
433 0 : casacore::Matrix<T> worldList;
434 0 : worldList = itsList;
435 :
436 0 : for (casacore::uInt g = 0; g < itsNComponents; g++)
437 : {
438 0 : casacore::Vector<casacore::Double> centercoords(itsDim); //casacore::CoordinateSystem uses casacore::Double only
439 0 : casacore::Vector<casacore::Double> compwidth(itsDim);
440 0 : for (casacore::uInt d = 0; d < itsDim; d++)
441 : {
442 0 : centercoords(d) = itsList(g,1+d);
443 : }
444 0 : for (casacore::uInt d = 0; d < itsDim; d++)
445 : {
446 0 : compwidth(d) = itsList(g,1+itsDim+d);
447 : }
448 :
449 0 : itsImagePtr->coordinates().toWorld(centercoords, centercoords);
450 0 : itsImagePtr->coordinates().toWorld(compwidth, compwidth);
451 0 : itsImagePtr->coordinates().makeWorldRelative(compwidth);
452 :
453 0 : for (casacore::uInt d = 0; d < itsDim; d++)
454 : {
455 0 : worldList(g,1+d) = centercoords(d);
456 : }
457 0 : for (casacore::uInt d = 0; d < itsDim; d++)
458 : {
459 0 : if (itsDim == 2 && d == 1) continue; // 2d: axis ratio, not x-width
460 0 : worldList(g,1+itsDim+d) = compwidth(d);
461 : }
462 :
463 : }
464 :
465 0 : return worldList;
466 0 : }
467 :
468 : template <class T>
469 : void ImageDecomposer<T>::componentMap() const
470 : {
471 : // does nothing yet.
472 : return;
473 : }
474 :
475 :
476 :
477 :
478 : // Maxima functions
479 : // ----------------
480 :
481 : template <class T>
482 0 : T ImageDecomposer<T>::findAreaGlobalMax(casacore::IPosition blc, casacore::IPosition trc) const
483 : {
484 : T val;
485 0 : T maxval = 0.0;
486 0 : correctBlcTrc(blc,trc);
487 : //
488 : {
489 0 : casacore::IPosition pos(blc);
490 0 : decrement(pos);
491 0 : while (increment(pos,trc)) {
492 0 : val = getImageVal(pos);
493 0 : if (val > maxval) maxval = val;
494 : }
495 0 : }
496 : //
497 0 : return maxval;
498 : }
499 :
500 :
501 : template <class T>
502 : void ImageDecomposer<T>::findAreaGlobalMax(T& maxval, casacore::IPosition& maxvalpos,
503 : casacore::IPosition blc, casacore::IPosition trc) const
504 : {
505 : T val;
506 :
507 : maxvalpos = casacore::IPosition(itsDim,0);
508 : maxval = 0.0;
509 : correctBlcTrc (blc,trc);
510 : //
511 : {
512 : casacore::IPosition pos(blc); decrement(pos);
513 : while (increment(pos,trc)) {
514 : val = getImageVal(pos);
515 : if (val > maxval) {maxval = val; maxvalpos = pos;}
516 : }
517 : }
518 :
519 : }
520 :
521 : template <class T>
522 : casacore::Vector<T> ImageDecomposer<T>::findAreaLocalMax(casacore::IPosition blc, casacore::IPosition trc,
523 : casacore::Int naxis) const
524 : {
525 : casacore::uInt const blocksize = 10;
526 : casacore::uInt maxn = 0;
527 : casacore::Vector<T> maxvals;
528 : correctBlcTrc (blc, trc);
529 : //
530 : {
531 : casacore::IPosition pos(blc);
532 : decrement(pos);
533 : while (increment(pos,trc)) {
534 : if (isLocalMax(pos,naxis)) {
535 : if (maxn % blocksize == 0) {
536 : maxvals.resize(maxn+blocksize, true);
537 : }
538 : maxvals(maxn) = getImageVal(pos);
539 : maxn++;
540 : }
541 : }
542 : }
543 : maxvals.resize(maxn, true);
544 : return maxvals;
545 : }
546 :
547 :
548 : template <class T>
549 : void ImageDecomposer<T>::findAreaLocalMax(casacore::Vector<T>& maxvals,
550 : casacore::Block<casacore::IPosition>& maxvalpos,
551 : casacore::IPosition blc, casacore::IPosition trc,
552 : casacore::Int naxis) const
553 : {
554 : casacore::uInt const blocksize = 10;
555 : casacore::uInt maxn = 0;
556 : maxvals.resize();
557 : maxvalpos.resize(0);
558 : correctBlcTrc(blc, trc);
559 : //
560 : {
561 : casacore::IPosition pos(blc);
562 : decrement(pos);
563 : while (increment(pos,trc)) {
564 : if (isLocalMax(pos,naxis)) {
565 : if (maxn % blocksize == 0) {
566 : maxvals.resize(maxn+blocksize, true);
567 : maxvalpos.resize(maxn+blocksize, false, true);
568 : }
569 : maxvals(maxn) = getImageVal(pos);
570 : maxvalpos[maxn] = pos;
571 : maxn++;
572 : }
573 : }
574 : }
575 : maxvals.resize(maxn, true);
576 : maxvalpos.resize(maxn, true, true);
577 : return;
578 : }
579 :
580 :
581 :
582 : template <class T>
583 : casacore::Vector<T> ImageDecomposer<T>::findRegionLocalMax(casacore::Int regionID, casacore::Int naxis) const
584 : {
585 : casacore::uInt const blocksize = 10;
586 : casacore::uInt maxn = 0;
587 : casacore::Vector<T> maxvals;
588 : {
589 : casacore::IPosition pos(itsDim,0);
590 : decrement(pos);
591 : while (increment(pos,shape())) {
592 : if ((getCell(pos) == regionID) && isLocalMax(pos,naxis)) {
593 : if (maxn % blocksize == 0) {
594 : maxvals.resize(maxn+blocksize, true);
595 : }
596 : maxvals(maxn) = getImageVal(pos);
597 : maxn++;
598 : }
599 : }
600 : }
601 : maxvals.resize(maxn, true);
602 : return maxvals;
603 : }
604 :
605 :
606 : template <class T>
607 0 : void ImageDecomposer<T>::findRegionLocalMax(casacore::Vector<T>& maxvals,
608 : casacore::Block<casacore::IPosition>& maxvalpos,
609 : casacore::Int regionID, casacore::Int naxis) const
610 : {
611 0 : casacore::uInt const blocksize = 10;
612 0 : casacore::uInt maxn = 0;
613 0 : maxvals.resize();
614 0 : maxvalpos.resize(0);
615 : //
616 : {
617 0 : casacore::IPosition pos(itsDim,0);
618 0 : decrement(pos);
619 0 : while (increment(pos,shape())) {
620 0 : if ((getCell(pos) == regionID) && isLocalMax(pos,naxis)) {
621 0 : cout << "Local max at " << pos << endl;
622 0 : if (maxn % blocksize == 0) {
623 0 : maxvals.resize(maxn+blocksize, true);
624 0 : maxvalpos.resize(maxn+blocksize, false, true);
625 : }
626 0 : maxvals(maxn) = getImageVal(pos);
627 0 : maxvalpos[maxn] = pos;
628 0 : maxn++;
629 : }
630 : }
631 0 : }
632 : //
633 0 : maxvals.resize(maxn, true);
634 0 : maxvalpos.resize(maxn, true, true);
635 : //
636 0 : return;
637 : }
638 :
639 : template <class T>
640 0 : casacore::Vector<T> ImageDecomposer<T>::findAllRegionGlobalMax() const
641 : {
642 : //NOTE: while the regions are identified in the itsMapPtr with #s starting at
643 : //one, the array returned by this function begin with zero, so there is
644 : //an offset of one between itsMapPtr IDs and those used by this function.
645 :
646 : casacore::Int r;
647 : T val;
648 0 : casacore::Vector<T> maxval(itsNRegions);
649 0 : maxval = 0.0;
650 : //
651 : {
652 0 : casacore::IPosition pos(itsDim,0);
653 0 : decrement(pos);
654 0 : while (increment(pos,shape())) {
655 0 : r = getCell(pos);
656 0 : if (r > 0) {
657 0 : val = getImageVal(pos);
658 0 : if (val > maxval(r-1)) maxval(r-1) = val;
659 : }
660 : }
661 0 : }
662 0 : return maxval;
663 0 : }
664 :
665 : template <class T>
666 0 : void ImageDecomposer<T>::findAllRegionGlobalMax(casacore::Vector<T>& maxvals,
667 : casacore::Block<casacore::IPosition>& maxvalpos) const
668 : {
669 : //NOTE: while the regions are identified in the itsMapPtr with #s starting at
670 : //one, the arrays returned by this function begin with zero, so there is
671 : //an offset of one between itsMapPtr IDs and those used by this function.
672 :
673 : casacore::Int r;
674 : T val;
675 :
676 0 : maxvals.resize(itsNRegions);
677 0 : maxvalpos.resize(itsNRegions);
678 0 : maxvals = 0; //note: wholly negative images still return 0
679 :
680 : {
681 0 : casacore::IPosition pos(itsDim,0);
682 0 : decrement(pos);
683 0 : while (increment(pos,shape())) {
684 0 : r = getCell(pos);
685 0 : if (r > 0) {
686 0 : val = getImageVal(pos);
687 0 : if (val > maxvals(r-1)) {
688 0 : maxvals(r-1) = val;
689 0 : maxvalpos[r-1] = pos;
690 : }
691 : }
692 : }
693 0 : }
694 :
695 0 : return;
696 : }
697 :
698 :
699 : template <class T>
700 0 : bool ImageDecomposer<T>::isLocalMax(const casacore::IPosition& pos, casacore::Int naxis) const
701 : {
702 0 : if (pos.nelements()==2) {
703 0 : return isLocalMax(pos(0), pos(1), naxis);
704 0 : } else if (pos.nelements()==3) {
705 0 : return isLocalMax(pos(0), pos(1), pos(2),naxis);
706 : } else {
707 0 : throw(casacore::AipsError("ImageDecomposer<T>::localmax(IPosition pos, Int naxis)"
708 : " - pos has wrong number of dimensions"));
709 : }
710 : return false;
711 : }
712 :
713 : template <class T>
714 0 : bool ImageDecomposer<T>::isLocalMax(casacore::Int x, casacore::Int y, casacore::Int naxis) const
715 : {
716 0 : T val = getImageVal(x,y);
717 0 : casacore::Int ximin = (x>0)? -1:0;
718 0 : casacore::Int yimin = (y>0)? -1:0;
719 0 : casacore::Int ximax = (x+1<shape(0))? 1:0;
720 0 : casacore::Int yimax = (y+1<shape(1))? 1:0;
721 0 : for (casacore::Int xi=ximin; xi<=ximax; xi++) {
722 0 : for (casacore::Int yi=yimin; yi<=yimax; yi++) {
723 0 : if ( ((naxis > 0) || !(xi || yi))
724 0 : && ((naxis > 1) || !(xi && yi))
725 0 : && (getImageVal(x+xi,y+yi) > val)) {
726 0 : return false;
727 : }
728 : }
729 : }
730 : //
731 0 : return true;
732 : }
733 :
734 : template <class T>
735 0 : bool ImageDecomposer<T>::isLocalMax(casacore::Int x, casacore::Int y, casacore::Int z, casacore::Int naxis) const
736 : {
737 0 : T maxval = getImageVal(x,y,z);
738 0 : casacore::Int ximin = (x>0)? -1:0;
739 0 : casacore::Int yimin = (y>0)? -1:0;
740 0 : casacore::Int zimin = (z>0)? -1:0;
741 0 : casacore::Int ximax = (x+1<shape(0))? 1:0;
742 0 : casacore::Int yimax = (y+1<shape(1))? 1:0;
743 0 : casacore::Int zimax = (z+1<shape(2))? 1:0;
744 0 : for (casacore::Int xi=ximin; xi<=ximax; xi++) {
745 0 : for (casacore::Int yi=yimin; yi<=yimax; yi++) {
746 0 : for (casacore::Int zi=zimin; zi<=zimax; zi++) {
747 0 : if ( ((naxis > 0) || !(xi || yi || zi))
748 0 : && ((naxis > 1) || !((xi && yi) || (xi && zi) || (yi && zi) ))
749 0 : && ((naxis > 2) || !(xi && yi && zi))
750 0 : && (getImageVal(x+xi,y+yi,z+zi) > maxval)) {
751 0 : return false;
752 : }
753 : }
754 : }
755 : }
756 : //
757 0 : return true;
758 : }
759 :
760 :
761 : // Estimation Functions
762 : // --------------------
763 :
764 :
765 : template <class T>
766 0 : void ImageDecomposer<T>::estimateComponentWidths(casacore::Matrix<T>& width,
767 : const casacore::Block<casacore::IPosition>& maxvalpos)
768 : const
769 : {
770 : // Finds a rough estimate of the width of each component.
771 : // Requires the location of each component.
772 :
773 : // This function is now mostly obsolete; its functionality replaced by
774 : // calculateMoments() except on non-deblended images.
775 :
776 0 : width.resize(maxvalpos.nelements(), itsDim);
777 : casacore::Bool dblflag;
778 : //
779 0 : for (casacore::uInt r = 0; r < maxvalpos.nelements(); r++) {
780 0 : casacore::IPosition lpos(itsDim);
781 0 : casacore::IPosition rpos(itsDim);
782 0 : casacore::IPosition maxpos(itsDim);
783 0 : maxpos = maxvalpos[r];
784 0 : T maxvalr = getImageVal(maxpos);
785 0 : T thrval = maxvalr*0.25;
786 : T val, prevval;
787 0 : for (casacore::uInt a = 0; a < itsDim; a++) {
788 0 : dblflag = 0;
789 0 : lpos = maxpos;
790 0 : val = maxvalr;
791 0 : prevval = val;
792 0 : while ((lpos(a) > 0) && (val >= thrval) && (val <= prevval)) {
793 0 : prevval = val;
794 0 : lpos(a) --;
795 0 : val = getImageVal(lpos);
796 : }
797 0 : if (val < thrval) {
798 0 : width(r,a) = T(maxpos(a)-lpos(a)) - (thrval-val) / (prevval-val);
799 0 : } else if (val > prevval) {
800 0 : width(r,a) = T(maxpos(a)-lpos(a));
801 : } else { //lpos == 0
802 0 : {width(r,a) = 0; dblflag = 1;} //can't find left limit;
803 : //use 2xright limit instead
804 : }
805 :
806 0 : rpos = maxpos;
807 0 : val = maxvalr;
808 0 : prevval = val;
809 0 : while ((rpos(a)<shape(a)-1) && (val >= thrval) && (val <= prevval)) {
810 0 : prevval = val;
811 0 : rpos(a) ++;
812 0 : val = getImageVal(rpos);
813 : }
814 0 : if (val < thrval) {
815 0 : width(r,a) += T(rpos(a)-maxpos(a)) - (thrval-val) / (prevval-val);
816 0 : } else if (val > prevval) {
817 0 : width(r,a) += T(rpos(a)-maxpos(a));
818 : } else {
819 0 : if (!dblflag) {
820 0 : dblflag = 1; //use 2x left side instead
821 : } else {
822 0 : width(r,a) += T(rpos(a)-maxpos(a)) * (maxvalr-thrval)/(maxvalr-val);
823 0 : dblflag = 1;
824 : }
825 : }
826 : //
827 0 : if (width(r,a) <= 0.0) width(r,a) = shape(a);//gaussian bigger than image
828 0 : if (!dblflag) width(r,a) /= 2.0;
829 0 : if (casacore::isNaN(width(r,a)))
830 : {
831 0 : width(r,a) = 1.0;
832 0 : cerr << "WARNING: Nonphysical estimate, setting width to 1.0" << endl;
833 : }
834 :
835 : }
836 : }
837 : //
838 0 : return;
839 : }
840 :
841 : template <class T>
842 0 : casacore::Array<T> ImageDecomposer<T>::calculateMoments(casacore::Int region) const
843 : {
844 : // Calculates the moments of an image region up to second order.
845 :
846 : // The current implementation is inefficient because it must scan the entire
847 : // image for each region. It would be better to return an array of Arrays,
848 : // or a casacore::Matrix with the M array collapsed to 1D and the region # along the
849 : // other axis.
850 :
851 0 : casacore::IPosition pos(itsDim);
852 0 : casacore::IPosition start(itsDim,0);
853 0 : decrement(start);
854 : T I;
855 :
856 0 : if (itsDim == 2)
857 : {
858 0 : casacore::Matrix<T> M(3,3, 0.0);
859 0 : pos = start;
860 0 : while (increment(pos,shape())) {
861 0 : if (getCell(pos) == region)
862 : {
863 0 : I = getImageVal(pos);
864 0 : M(0,0) += I;
865 0 : M(1,0) += pos(0) * I;
866 0 : M(0,1) += pos(1) * I;
867 : }
868 : }
869 :
870 :
871 0 : T xc = M(1,0) / M(0,0);
872 0 : T yc = M(0,1) / M(0,0);
873 :
874 0 : pos = start;
875 0 : while (increment(pos,shape())) {
876 0 : if (getCell(pos) == region)
877 : {
878 0 : I = getImageVal(pos);
879 0 : T xd = pos(0) - xc;
880 0 : T yd = pos(1) - yc;
881 0 : M(1,1) += xd * yd * I;
882 0 : M(2,0) += xd * xd * I;
883 0 : M(0,2) += yd * yd * I;
884 : }
885 : }
886 :
887 0 : return M;
888 : //Does not calculate higher level components (2,1; 1,2; 2,2)
889 0 : }
890 :
891 0 : if (itsDim == 3)
892 : {
893 0 : casacore::Cube<T> M(3,3,3, 0.0);
894 :
895 0 : pos = start;
896 0 : while (increment(pos,shape())) {
897 0 : if (getCell(pos) == region)
898 : {
899 0 : I = getImageVal(pos);
900 0 : M(0,0,0) += I;
901 0 : M(1,0,0) += pos(0) * I;
902 0 : M(0,1,0) += pos(1) * I;
903 0 : M(0,0,1) += pos(2) * I;
904 : }
905 : }
906 :
907 :
908 0 : T xc = M(1,0,0) / M(0,0,0);
909 0 : T yc = M(0,1,0) / M(0,0,0);
910 0 : T zc = M(0,0,1) / M(0,0,0);
911 0 : pos = start;
912 0 : while (increment(pos,shape())) {
913 0 : if (getCell(pos) == region)
914 : {
915 0 : I = getImageVal(pos);
916 0 : T xd = pos(0) - xc;
917 0 : T yd = pos(1) - yc;
918 0 : T zd = pos(2) - zc;
919 0 : M(1,1,0) += xd * yd * I;
920 0 : M(1,0,1) += xd * zd * I;
921 0 : M(0,1,1) += yd * zd * I;
922 0 : M(2,0,0) += xd * xd * I;
923 0 : M(0,2,0) += yd * yd * I;
924 0 : M(0,0,2) += zd * zd * I;
925 : }
926 : }
927 :
928 0 : return M;
929 0 : }
930 :
931 0 : return casacore::Array<T>();
932 0 : }
933 :
934 :
935 : // Region editing functions
936 : // ------------------------
937 :
938 : template <class T>
939 0 : void ImageDecomposer<T>::boundRegions(casacore::Block<casacore::IPosition>& blc,
940 : casacore::Block<casacore::IPosition>& trc)
941 : {
942 : // Boxes each region in the componentmap:
943 : // blc is set to the lowest coordinate value in each region;
944 : // trc is set to one above the highest coordinate value in each region.
945 :
946 0 : DebugAssert(blc.nelements() == itsNRegions, casacore::AipsError);
947 0 : DebugAssert(trc.nelements() == itsNRegions, casacore::AipsError);
948 :
949 0 : for (casacore::uInt r=0; r<itsNRegions; r++) {
950 0 : blc[r] = itsShape;
951 0 : trc[r] = casacore::IPosition(itsDim,0);
952 : }
953 :
954 : {
955 0 : casacore::IPosition pos(itsDim,0);
956 0 : decrement(pos);
957 : casacore::Int r;
958 0 : while (increment(pos,shape())) {
959 0 : r = getCell(pos);
960 0 : if (r > 0) {
961 0 : for (casacore::uInt i = 0; i < itsDim; i++) {
962 0 : if (blc[r-1](i) > pos(i)) blc[r-1](i) = pos(i);
963 0 : if (trc[r-1](i) <= pos(i)) trc[r-1](i) = pos(i)+1;
964 : }
965 : }
966 : }
967 0 : }
968 :
969 0 : return;
970 : }
971 :
972 : template <class T>
973 0 : void ImageDecomposer<T>::zero()
974 : {
975 : // Set all elements in the component map to zero and clear the component list.
976 :
977 0 : itsMapPtr->set(0);
978 0 : itsNRegions = 0;
979 0 : itsNComponents = 0;
980 0 : itsList.resize();
981 0 : return;
982 : }
983 :
984 : template <class T>
985 0 : void ImageDecomposer<T>::clear()
986 : {
987 : //Clear the component map
988 0 : casacore::LatticeIterator<casacore::Int> iter(*itsMapPtr);
989 : casacore::Bool deleteIt;
990 0 : casacore::Int* p = 0;
991 0 : for (iter.reset(); !iter.atEnd(); iter++) {
992 0 : casacore::Array<casacore::Int>& tmp = iter.rwCursor();
993 0 : p = tmp.getStorage(deleteIt);
994 0 : for (casacore::uInt i=0; i<tmp.nelements(); i++) if (p[i] != MASKED) p[i] = 0;
995 0 : tmp.putStorage(p, deleteIt);
996 : }
997 0 : itsNRegions = 0;
998 :
999 : //Clear the component list
1000 0 : itsNComponents = 0;
1001 0 : itsList.resize();
1002 0 : return;
1003 0 : }
1004 :
1005 :
1006 : template <class T>
1007 0 : void ImageDecomposer<T>::destroyRegions(const casacore::Vector<casacore::Bool>& killRegion)
1008 : {
1009 : //Wipes out any regions whose corresponding values in killRegion are true
1010 : // by setting all pixel values in the componentmap set to that region to
1011 : // zero. Zero-oriented; there is an offset of one between the index in
1012 : // killRegion and the actual region in the componentmap.
1013 :
1014 : {
1015 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1016 0 : while (increment(pos,shape()))
1017 : {
1018 0 : casacore::Int reg = getCell(pos);
1019 0 : if (reg > 0 && killRegion(reg-1)) setCell(pos,0);
1020 : }
1021 0 : }
1022 :
1023 0 : renumberRegions();
1024 0 : return;
1025 : }
1026 :
1027 :
1028 : template <class T>
1029 0 : void ImageDecomposer<T>::renumberRegions()
1030 : {
1031 : //Eliminates redundant regions (components with no representative cells in
1032 : //the component map) by renumbering higher-numbered regions to fill in
1033 : //the gaps. For example..
1034 : // 011 011
1035 : // 113 becomes 112
1036 : // 113 112
1037 :
1038 0 : casacore::Vector<casacore::Bool> regpresent(itsNRegions+1, 0);
1039 0 : casacore::Vector<casacore::Int> renumregs(itsNRegions+1);
1040 0 : casacore::uInt const ngpar = itsDim * 3;
1041 :
1042 : //any region that has any pixel members is flagged as 1, others left 0
1043 :
1044 : {
1045 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1046 0 : while (increment(pos,shape()))
1047 : {
1048 0 : casacore::Int reg = getCell(pos);
1049 0 : if (reg >= 0) regpresent(reg) = 1;
1050 : }
1051 0 : }
1052 :
1053 : //determine new # of regions and which regions will be renumbered to what
1054 0 : casacore::uInt newnr = 0;
1055 0 : for (casacore::uInt r = 1; r <= itsNRegions; r++)
1056 0 : if (regpresent(r)) renumregs(r) = ++newnr;
1057 :
1058 0 : if (newnr < itsNRegions)
1059 : {
1060 0 : itsNRegions = newnr;
1061 :
1062 : //actually renumber the regions in the pmap
1063 : {
1064 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1065 0 : while (increment(pos,shape()))
1066 : {
1067 0 : casacore::Int reg = getCell(pos);
1068 0 : if (reg >= 0) setCell(pos, renumregs(reg));
1069 : }
1070 0 : }
1071 :
1072 0 : if (isDecomposed())
1073 : {
1074 : //eliminate componentlist entries of lost components
1075 0 : casacore::Matrix<T> oldcomponentlist(itsList);
1076 :
1077 0 : itsList.resize(newnr, ngpar);
1078 0 : for (casacore::Int c = 0; c < casacore::Int(itsNComponents); c++)
1079 0 : if (regpresent(c+1) && (c+1 != renumregs(c+1)))
1080 0 : for (casacore::uInt p = 0; p < 9; p++)
1081 0 : itsList(renumregs(c+1)-1,p) = oldcomponentlist(c+1,p);
1082 :
1083 0 : itsNComponents = newnr;
1084 0 : }
1085 : }
1086 :
1087 0 : return;
1088 0 : }
1089 :
1090 : template <class T>
1091 0 : void ImageDecomposer<T>::synthesize(const ImageDecomposer<T>& subdecomposer,
1092 : casacore::IPosition blc)
1093 : {
1094 : // Overlays a smaller map onto an empty region of a larger map,
1095 : // and adds submap component list to main component list.
1096 :
1097 : // The user should exercise caution with this function and synthesize submaps
1098 : // only into regions of the main map that are truly empty (0), because the
1099 : // program does not perform any blending between components.
1100 : // Otherwise, false detections are likely.
1101 :
1102 0 : casacore::uInt ngpar = 0;
1103 0 : if (itsDim == 2) ngpar = 6;
1104 0 : if (itsDim == 3) ngpar = 9;
1105 :
1106 : // Scan to the edge of the boundary of the host map or the submap, whichever
1107 : // is closer.
1108 :
1109 0 : casacore::IPosition scanlimit(itsDim); //submap-indexed
1110 0 : for (casacore::uInt i=0; i<itsDim; i++) {
1111 0 : if (subdecomposer.shape(i) > shape(i) - blc(i)) {
1112 0 : scanlimit(i) = shape(i) - blc(i);
1113 : } else {
1114 0 : scanlimit(i) = subdecomposer.shape(i);
1115 : }
1116 : }
1117 :
1118 : // Write pixels in sub- component map to main component map.
1119 :
1120 : {
1121 0 : casacore::IPosition pos(itsDim,0); //submap-indexed
1122 0 : decrement(pos);
1123 0 : while (increment(pos,scanlimit)) {
1124 0 : if (subdecomposer.getCell(pos) > 0) {
1125 0 : setCell(pos + blc, itsNRegions + subdecomposer.getCell(pos));
1126 : }
1127 : }
1128 0 : }
1129 0 : itsNRegions += subdecomposer.numRegions();
1130 :
1131 : // Add components in subdecomposer to components in main decomposer.
1132 :
1133 0 : if (subdecomposer.isDecomposed()) {
1134 :
1135 0 : casacore::Matrix<T> oldList;
1136 0 : oldList = itsList;
1137 0 : itsList.resize(itsNComponents+subdecomposer.numComponents(),ngpar);
1138 0 : for (casacore::uInt c = 0; c < itsNComponents; c++) {
1139 0 : for (casacore::uInt p = 0; p < ngpar; p++) {
1140 0 : itsList(c,p) = oldList(c,p); //copy after resize
1141 : }
1142 : }
1143 :
1144 0 : for (casacore::uInt subc = 0; subc < subdecomposer.numComponents(); subc++) {
1145 0 : for (casacore::uInt p = 0; p < ngpar; p++) {
1146 0 : itsList(itsNComponents+subc,p)=subdecomposer.itsList(subc,p);
1147 : }
1148 : // make adjustments to center values due to offset
1149 0 : if (itsDim == 2) {
1150 0 : itsList(itsNComponents+subc,1) += blc(0);
1151 0 : itsList(itsNComponents+subc,2) += blc(1);
1152 0 : } else if (itsDim == 3) {
1153 0 : itsList(itsNComponents+subc,1) += blc(0);
1154 0 : itsList(itsNComponents+subc,2) += blc(1);
1155 0 : itsList(itsNComponents+subc,3) += blc(2);
1156 : }
1157 : }
1158 0 : itsNComponents += subdecomposer.numComponents();
1159 :
1160 0 : }
1161 : //
1162 : //renumberRegions(); //this should have no effect unless this function
1163 : //was used to overwrite an entire region, which
1164 : //should be avoided.
1165 0 : }
1166 :
1167 :
1168 :
1169 : // Analysis functions
1170 : // ------------------
1171 :
1172 :
1173 : template <class T>
1174 0 : casacore::uInt ImageDecomposer<T>::identifyRegions(T thrval, casacore::Int naxis)
1175 : {
1176 : // Performs a single threshold scan on the image. In other words,
1177 : // identifies all contigous blocks of pixels in the target image above the
1178 : // threshold value thrval, assigning each unique block to an integer,
1179 : // starting at one. All pixels with target image values below thrval are set
1180 : // to zero.
1181 :
1182 : // NOTE: Formerly a specialization existed for 2D which may have been slightly
1183 : // more efficient. However, it complicated the code excessively (this
1184 : // program is already far too long.) It could be resurrected if necessary.
1185 :
1186 0 : casacore::Int const blocksize = 1024; //increment to grow size of anchor array
1187 0 : casacore::Int const pageexpsize = 128;
1188 0 : casacore::Int cnum = 0; //region number
1189 0 : if (naxis > casacore::Int(itsDim)) naxis = itsDim;
1190 :
1191 : // The program first scans through the image until it finds any pixel in
1192 : // any region. Once there, it immediately scans all 6 adjacent pixels,
1193 : // registering any that fall within the region and setting them as anchors.
1194 : // It then moves to the first of these anchors and repeats the process,
1195 : // until no more anchors can be found and every existing anchor has already
1196 : // been explored. It then scans until it locates a new region, and repeats
1197 : // until every region has been similarly scanned.
1198 :
1199 : // This has the equivalent effect of a 'surface' of active anchors
1200 : // sweeping across each region starting at a seed point until the
1201 : // region is fully explored.
1202 :
1203 : // The naxis parameter must be 2 or greater to avoid spurious detections
1204 : // along horizontal ridges. However, this slows down performance by a factor
1205 : // of roughly two, so in instances where objects are clearly defined and
1206 : // not closely blended and the image is very large naxis=1 may be better.
1207 :
1208 0 : casacore::IPosition scanpos(itsDim,0); //decrement(scanpos);
1209 0 : while (true) {
1210 : //First find any pixel in next region.
1211 : //Stop scanning when an unassigned, unmasked pixel is found exceeding
1212 : //the threshold value.
1213 :
1214 0 : while (getImageVal(scanpos) < thrval || getCell(scanpos)) {
1215 0 : if (!increment(scanpos,shape())) {
1216 0 : return itsNRegions = cnum;
1217 : }
1218 : //scanned w/out finding new region
1219 : }
1220 : //
1221 0 : casacore::IPosition pos(scanpos);
1222 0 : cnum++;
1223 :
1224 : // As many anchors will be required as pixels in the region (the volume) -
1225 : // but only a small fraction of that at any one time (the active surface).
1226 : // So the program allocates 'pages' of anchors as they become necessary,
1227 : // but continually deletes pages of anchors that have already been used.
1228 :
1229 0 : casacore::Int seta = -1; //index of highest established anchor
1230 0 : casacore::Int geta = -1; //index of highest analyzed anchor
1231 : //the active surface is all anchors geta < a < seta
1232 0 : casacore::PtrBlock<casacore::Matrix<casacore::Int> *> aindex(pageexpsize); //anchor structure
1233 0 : casacore::Int setblock = -1; //index of page containing highest established anchor
1234 0 : casacore::Int getblock = -1; //index of page containing highest analyzed anchor
1235 0 : setCell(pos,cnum);
1236 :
1237 : do { //(geta < seta)
1238 : //cout << geta << " / " << seta << ", " << pos <<
1239 : // " = " << getCell(pos) << endl;
1240 :
1241 : //Analyze the cell -
1242 : //Scan the cells around the active cell as follows:
1243 : //naxis = 1: scan 6 adjacent cells (axes only)
1244 : //naxis = 2: scan 18 adjacent cells (axes and 2-axis diagonals)
1245 : //naxis = 3: scan 26 adjacent cells (axes and 2/3-axis diagonals)
1246 :
1247 0 : casacore::Int ximin = (pos(0)>0)? -1:0;
1248 0 : casacore::Int yimin = (pos(1)>0)? -1:0;
1249 0 : casacore::Int zimin = ((itsDim>=3)&&(pos(2)>0))? -1:0;
1250 0 : casacore::Int ximax = (pos(0)+1<shape(0))? 1:0;
1251 0 : casacore::Int yimax = (pos(1)+1<shape(1))? 1:0;
1252 0 : casacore::Int zimax = ((itsDim>=3)&&(pos(2)+1<shape(2)))? 1:0; //safe for 2D
1253 : //
1254 0 : for (casacore::Int xi=ximin; xi<=ximax; xi++) {
1255 0 : for (casacore::Int yi=yimin; yi<=yimax; yi++) {
1256 0 : for (casacore::Int zi=zimin; zi<=zimax; zi++) {
1257 0 : if ( (xi || yi || zi) &&
1258 0 : ((naxis > 1) || !((xi && yi) || (xi && zi) || (yi && zi) )) &&
1259 0 : ((naxis > 2) || !(xi && yi && zi))) {
1260 0 : casacore::IPosition ipos(pos);
1261 0 : ipos(0) += xi; ipos(1) += yi; if (itsDim == 3) ipos(2) += zi;
1262 :
1263 : //if any contiguous pixel is
1264 0 : if ((getImageVal(ipos) >= thrval) // above threshold and
1265 0 : && getCell(ipos) != cnum) { // not yet scanned...
1266 : //record its location as an anchor and come back later.
1267 0 : seta++;
1268 :
1269 0 : if ((seta % blocksize) == 0) {
1270 :
1271 : //current block out of memory: allocate new memory block
1272 :
1273 0 : setblock++;
1274 0 : if ((setblock % pageexpsize == 0) && setblock) {
1275 0 : aindex.resize(((setblock/pageexpsize)+1)*pageexpsize);
1276 : }
1277 0 : aindex[setblock] = new casacore::Matrix<casacore::Int>(blocksize,itsDim);
1278 : }
1279 :
1280 : //set new anchor
1281 0 : for (casacore::uInt axis = 0; axis < itsDim; axis ++) {
1282 0 : (*(aindex[setblock]))(seta%blocksize,axis) = ipos(axis);
1283 : }
1284 :
1285 : //cout<<'A'<<seta<<'['<<x+xi<<','<<y+yi<<','<<z+zi<<']'<<endl;
1286 :
1287 0 : setCell(ipos, cnum);
1288 : }
1289 0 : }
1290 : }
1291 : }
1292 : }
1293 0 : geta++;
1294 :
1295 :
1296 0 : if (((geta % blocksize) == 0) || (geta>seta)) {
1297 0 : if (getblock>=0) {
1298 : //memory block data obsolete: delete old memory block
1299 0 : delete aindex[getblock];
1300 : }
1301 0 : getblock++;
1302 : }
1303 :
1304 0 : if (geta <= seta) {
1305 : //go to lowest anchor
1306 0 : for (casacore::uInt axis = 0; axis < itsDim; axis++) {
1307 0 : pos(axis) = (*(aindex[getblock]))(geta%blocksize,axis);
1308 : }
1309 : //cout << ">>A" << geta << " " << endl;
1310 : }
1311 :
1312 0 : } while (geta <= seta);
1313 : }
1314 0 : }
1315 :
1316 : template <class T>
1317 0 : void ImageDecomposer<T>::deblendRegions(const casacore::Vector<T>& contours,
1318 : casacore::Int minRange, casacore::Int naxis)
1319 : {
1320 :
1321 : // Performs the contour decomposition on a blended image to generate a
1322 : // component map that can detect components blended above any threshold(s),
1323 : // by performing threshold scans at each contour level and recognizing
1324 : // as individual any components that are distinct above any such level.
1325 :
1326 0 : casacore::Int const printIntermediates = 0;
1327 0 : casacore::Int const blocksize = 3;
1328 0 : casacore::Int ncontours = contours.nelements();
1329 :
1330 0 : ImageDecomposer<T> contourMap(*this);//Component map thresholded at current
1331 : //contour level; "lower" contour
1332 0 : casacore::Block<casacore::IPosition> regcenter; //Coordinates of first pixel found in each
1333 : //component (a rough estimate of the center)
1334 : //Indexing for this array is offset by one.
1335 0 : casacore::Vector<casacore::Int> originlevel; //first contour in which region was detected
1336 : //Indexing for this array is offset by one.
1337 : //If set to -1 the region is defunct.
1338 :
1339 :
1340 0 : if (isDerived()) zero();
1341 :
1342 : // Component decomposition:
1343 : // This is the main deblending algorithm. The program starts by performing
1344 : // a threshold scan at the top contour value, separating any regions
1345 : // that are visibly distinct at that threshold. It continues downward though
1346 : // the contour vector, forming a component map from each contour and comparing
1347 : // that with the contour map that is gradually forming to distinguish new
1348 : // regions, correspondent regions, and blending regions, and assigns the
1349 : // pixels in the main itsMapPtr based on this information.
1350 :
1351 0 : for (casacore::Int c = ncontours-1; c >= 0; c--) {
1352 : if (printIntermediates == 2) cout << endl << "CONTOUR " << c << endl;
1353 :
1354 : casacore::Int lowreg, highreg; //number of regions in lower contour, current pmap
1355 0 : contourMap.clear(); //only necessary if region grows between contours
1356 :
1357 0 : lowreg = contourMap.identifyRegions(contours(c),naxis);
1358 0 : highreg = itsNRegions;
1359 :
1360 : if (printIntermediates) {
1361 : cout << "Now comparing current pmap to contour " << c << '.' << endl;
1362 : display();
1363 : contourMap.display();
1364 : }
1365 :
1366 0 : casacore::Vector<casacore::Int> root(highreg+1, 0); //Indicates corresponding object
1367 : // in lower contour
1368 0 : casacore::Vector<casacore::Int> nOffspring(lowreg+1, 0); //Total number of different objects
1369 : // above this region
1370 0 : casacore::Block<casacore::Vector<casacore::Int>*> offspring(lowreg+1); //Region IDs of objects above
1371 : // this region
1372 :
1373 : // Can't finish allocation until nOffspring is known.
1374 :
1375 : // Scan through current pmap ("higher contour") and find the root of all
1376 : // regions as defined in the lower contour. Simultaneously mark all
1377 : // regions in the lower contour that have "offspring" above.
1378 :
1379 0 : for (casacore::Int r = 1; r <= highreg; r++) {
1380 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1381 0 : while (increment(pos,shape()) && !root(r)){
1382 0 : if (getCell(pos) == r) {
1383 0 : root(r) = contourMap.getCell(pos);
1384 0 : nOffspring(contourMap.getCell(pos)) += 1;
1385 : }
1386 : }
1387 : }
1388 :
1389 : // Set up offspring array
1390 :
1391 0 : for (casacore::Int r = 1; r <= lowreg; r++) {
1392 0 : offspring[r] = new casacore::Vector<casacore::Int>(nOffspring(r));
1393 : }
1394 :
1395 0 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1396 0 : casacore::Int f = 0;
1397 0 : for (casacore::Int hr = 1; hr <= highreg; hr++) {
1398 0 : if (root(hr) == lr) (*offspring[lr])(f++) = hr;
1399 : }
1400 : }
1401 :
1402 : if (printIntermediates == 2) {
1403 : cout << "Contour Level " << c << endl;
1404 : cout << highreg << ' ' << lowreg << endl;
1405 :
1406 : for (casacore::Int hr = 1; hr <= highreg; hr++) {
1407 : cout << "root of " << hr << ':' << root(hr) << endl;
1408 : }
1409 :
1410 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1411 : for (casacore::Int f = 0; f < nOffspring(lr); f++) {
1412 : cout << "offspring of " << lr << ':' << (*offspring[lr])(f) << endl;
1413 : }
1414 : }
1415 : }
1416 :
1417 :
1418 :
1419 : // If a newly discovered component merges with another too quickly
1420 : // (within minRange contours) rescind its status and treat it as part
1421 : // of the other contour.
1422 : // "f" refers somewhat cryptically to offspring index and "fr" to its
1423 : // region number in the pmap
1424 :
1425 0 : if (minRange > 1) {
1426 0 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1427 0 : if (nOffspring(lr) > 1) {
1428 0 : for (casacore::Int f = 0; f < nOffspring(lr); f++) {
1429 0 : casacore::Int fr = (*offspring[lr])(f);
1430 :
1431 0 : if (originlevel(fr-1) - c < minRange)
1432 : {
1433 : //Find the closest offspring
1434 0 : casacore::Int mindistsq = 1073741823; //maximum casacore::Int value
1435 0 : casacore::Int frabs = 0; //closest region
1436 0 : for (casacore::Int f2 = 0; f2 < nOffspring(lr); f2++) {
1437 0 : if (f2 == f) continue;
1438 0 : casacore::Int fr2 = (*offspring[lr])(f2);
1439 0 : casacore::Int distsq = 0;
1440 0 : if (originlevel(fr2-1) == -1) continue;
1441 0 : for (casacore::uInt a = 0; a < itsDim; a++) {
1442 0 : distsq += casacore::square(casacore::Int(regcenter[fr2-1](a)-regcenter[fr-1](a)));
1443 : }
1444 0 : if (distsq < mindistsq) {
1445 0 : frabs = fr2;
1446 0 : mindistsq = distsq;
1447 : }
1448 : }
1449 0 : if (frabs == 0)
1450 : {
1451 : //No valid closest offspring - find biggest offspring
1452 0 : for (casacore::Int f2 = 0; f2 < nOffspring(lr); f2++) {
1453 0 : if (f2 == f) continue;
1454 0 : casacore::Int fr2 = (*offspring[lr])(f2);
1455 0 : casacore::Int maxlevel = 0;
1456 0 : if (originlevel(fr2-1) == -1) continue;
1457 0 : if (originlevel(fr2-1) > maxlevel) {
1458 0 : frabs = fr2;
1459 0 : maxlevel = originlevel(fr2-1);
1460 : }
1461 : }
1462 : }
1463 0 : if (frabs == 0)
1464 : {
1465 : //must be the only offspring left! Don't absorb.
1466 0 : break;
1467 : }
1468 :
1469 :
1470 : if (printIntermediates == 2) {
1471 : cout << "Absorbing region " << fr << " (origin at "
1472 : << originlevel(fr-1) << ") into region " << frabs << endl;
1473 : }
1474 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1475 0 : while (increment(pos,shape())){
1476 0 : if (getCell(pos) == fr) setCell(pos,frabs);
1477 : }
1478 0 : originlevel(fr-1) = -1;
1479 0 : }
1480 : }
1481 : }
1482 : }
1483 : }
1484 :
1485 : // Need to check if this works properly...
1486 :
1487 :
1488 :
1489 : // The object/region/component labeling is done in three steps to make the
1490 : // order most closely match the order in which they were found, and the
1491 : // order of the peak intensities.
1492 :
1493 : // First (and most complex) step is to deblend the large consolidated
1494 : // regions (nOffspring >= 2). The algorithm scans all pixels that
1495 : // exist in the new contour but not in the current pmap and belong to
1496 : // a region with multiple offspring, then scans the surrounding vicinity
1497 : // in the current pmap until it finds a pixel belonging to a region there,
1498 : // to which the new pixel is assigned.
1499 :
1500 0 : casacore::Int cgindex = 1; //component index base, global for this contour
1501 0 : ImageDecomposer<T> copyMap(*this); //The original data must be temporarily
1502 : //preserved while the componentmap is
1503 : //overwritten for comparisons between
1504 : //nearby pixels
1505 :
1506 0 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1507 0 : if (nOffspring(lr) >= 2) {
1508 0 : casacore::IPosition pos(itsDim,0); decrement(pos);
1509 0 : while (increment(pos,shape())) {
1510 :
1511 : // Renumber pixels that already exist in the pmap
1512 0 : if ((contourMap.getCell(pos)==lr)&&(copyMap.getCell(pos))) {
1513 0 : for (casacore::Int f = 0; f < nOffspring(lr); f++) {
1514 :
1515 : // Translate old pmap id to new pmap id
1516 0 : if ((*offspring[lr])(f) == copyMap.getCell(pos))
1517 0 : setCell(pos, cgindex+f);
1518 : }
1519 : }
1520 :
1521 : // Number new pixels.
1522 :
1523 0 : if ((contourMap.getCell(pos)==lr)&&(!copyMap.getCell(pos))) {
1524 0 : casacore::Int pinh = 0; //region as defined in pmap pixel is set to
1525 0 : casacore::Int srad = 1; //search radius
1526 0 : casacore::uInt naxis = 1;
1527 :
1528 : // cout << "Searching for nearest cell to " << pos << endl;
1529 :
1530 0 : while(!pinh && srad < 250) { //search increasing naxis, srad
1531 : // IMPR: an N-dimensional structure would be better here.
1532 : casacore::Int xi, yi, zi;
1533 0 : casacore::Int ximin = (pos(0)-srad < 0)? -pos(0) : -srad;
1534 0 : casacore::Int ximax = (pos(0)+srad >= shape(0))? shape(0)-pos(0)-1 : srad;
1535 0 : casacore::Int yimin = (pos(1)-srad < 0)? -pos(1) : -srad;
1536 0 : casacore::Int yimax = (pos(1)+srad >= shape(1))? shape(1)-pos(1)-1 : srad;
1537 0 : casacore::Int zimin = 0, zimax = 0;
1538 0 : if (itsDim == 2) {
1539 0 : zimin = 0; zimax = 0;
1540 : }
1541 0 : if (itsDim >= 3) {
1542 0 : zimin = (pos(2)-srad < 0)? -pos(2) : -srad;
1543 0 : zimax = (pos(2)+srad >= shape(2))? shape(2)-pos(2)-1 : srad;
1544 : }
1545 0 : while (!pinh && naxis <= itsDim) {
1546 0 : for (xi = ximin; xi <= ximax; xi++) {
1547 0 : for (yi = yimin; yi <= yimax; yi++) {
1548 0 : for (zi = zimin; zi <= zimax; zi++) {
1549 0 : casacore::IPosition ipos(pos);
1550 0 : ipos(0) += xi; ipos(1) += yi;
1551 0 : if (itsDim==3) ipos(2) += zi;
1552 :
1553 0 : if (abs(xi)<srad && abs(yi)<srad && abs(zi)<srad) {
1554 0 : continue; //border of radius only
1555 : }
1556 :
1557 0 : if ( ((naxis==1) && ((xi&&yi) || (xi&&zi) || (yi&&zi)))
1558 0 : || ((naxis==2) && (xi && yi && zi))) {
1559 0 : continue;
1560 : }
1561 :
1562 0 : casacore::Int inh = copyMap.getCell(ipos);
1563 0 : if (inh<=0) continue;
1564 0 : if (!pinh) {
1565 0 : pinh = inh;
1566 0 : for (casacore::Int f = 0; f < nOffspring(lr); f++) {
1567 :
1568 : // Translate old pmap id to new pmap id
1569 :
1570 0 : if ((*offspring[lr])(f) == inh) {
1571 0 : setCell(pos, cgindex+f);
1572 : }
1573 : //reassign pixel to new component
1574 : }
1575 0 : } else if (pinh!=inh) {
1576 : //equidistant to different objects:
1577 : // temporarily flag as nonexistant object
1578 0 : setCell(pos, INDETERMINATE);
1579 : }
1580 : }
1581 : }
1582 : }
1583 0 : naxis++;
1584 : }
1585 0 : naxis = 1; srad++;
1586 : }
1587 : }
1588 : }
1589 0 : cgindex += nOffspring(lr);
1590 0 : }
1591 : }
1592 :
1593 : // Now scan nonforked regions that exist in both contours.
1594 : // This is as simple as just renumbering the region.
1595 :
1596 0 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1597 0 : if (nOffspring(lr) == 1) {
1598 0 : casacore::IPosition pos(itsDim,0);
1599 0 : decrement(pos);
1600 0 : while (increment(pos,shape())) {
1601 0 : if (contourMap.getCell(pos) == lr) setCell(pos, cgindex);
1602 : }
1603 0 : cgindex++;
1604 0 : }
1605 : } //IMPR: probably can make this work with single scan
1606 : //similar to above algorithm
1607 :
1608 : // Finally, scan regions that only exist in lower contour. Same as above,
1609 : // but since these are new regions, add their initial positions to the seed
1610 : // arrays and increment the region count.
1611 :
1612 0 : for (casacore::Int lr = 1; lr <= lowreg; lr++) {
1613 0 : if (nOffspring(lr) == 0) {
1614 :
1615 0 : casacore::IPosition newregioncenter(itsDim,0);
1616 0 : casacore::uInt topcells = 0;
1617 : {
1618 0 : casacore::IPosition pos(itsDim,0);
1619 0 : decrement(pos);
1620 0 : while (increment(pos,shape())) {
1621 : //cout << pos << endl;
1622 0 : if (contourMap.getCell(pos) == lr) {
1623 0 : setCell(pos, cgindex);
1624 0 : newregioncenter += pos;
1625 0 : topcells++;
1626 : }
1627 : }
1628 0 : }
1629 0 : newregioncenter /= topcells; //note: integer division. may or may not
1630 : //want to keep this
1631 :
1632 0 : cgindex++;
1633 :
1634 0 : if ((itsNRegions % blocksize) == 0) {
1635 0 : regcenter.resize(itsNRegions+blocksize, true);
1636 0 : originlevel.resize(itsNRegions+blocksize, true);
1637 : }
1638 :
1639 : // Add to region center array
1640 0 : regcenter[itsNRegions] = newregioncenter;
1641 0 : originlevel(itsNRegions) = c;
1642 0 : itsNRegions++;
1643 0 : }
1644 : }
1645 : }
1646 :
1647 :
1648 : // At end of scan, assign all flagged pixels to region containing the
1649 : // nearest "seed" - the first point identified in each region.
1650 :
1651 : if (printIntermediates == 2) {
1652 : cout << "Located the following seeds:" << endl;
1653 : for (casacore::uInt s = 0; s < itsNRegions; s++)
1654 : cout << s << " at " << regcenter[s]
1655 : << " in component #" << getCell(regcenter[s]) << endl;
1656 : }
1657 :
1658 : {
1659 0 : casacore::IPosition pos(itsDim,0);
1660 0 : decrement(pos);
1661 0 : while (increment(pos,shape())) {
1662 0 : if (getCell(pos) == INDETERMINATE) {
1663 0 : casacore::Int mindistsq = 1073741823; //maximum casacore::Int value
1664 0 : for (casacore::uInt s = 0; s < itsNRegions; s++) {
1665 0 : casacore::Int distsq = 0;
1666 0 : if (originlevel(s) == -1) continue; //defunct region
1667 0 : for (casacore::uInt a = 0; a < itsDim; a++) {
1668 0 : distsq += (pos(a) - regcenter[s](a)) * (pos(a) - regcenter[s](a));
1669 : }
1670 0 : if (distsq < mindistsq) {
1671 0 : setCell(pos, getCell(regcenter[s]) );
1672 0 : mindistsq = distsq;
1673 : }
1674 : }
1675 : }
1676 : }
1677 0 : }
1678 :
1679 0 : if (minRange > 1) renumberRegions();
1680 :
1681 0 : return;
1682 0 : }
1683 :
1684 :
1685 :
1686 : template <class T>
1687 0 : void ImageDecomposer<T>::decomposeImage()
1688 : {
1689 0 : if (!itsDeblendIt)
1690 : {
1691 : // Find contiguous regions via thresholding
1692 :
1693 0 : casacore::uInt nRegions = identifyRegions(itsThresholdVal,itsNAxis);
1694 0 : cout << "Found " << nRegions << " regions" << endl;
1695 :
1696 : // Fit each region. A further pass is done through each region
1697 : // to ascertain whether there are multiple components to be
1698 : // fit within it.
1699 :
1700 0 : if (itsFitIt) fitRegions();
1701 :
1702 : }
1703 : else
1704 : {
1705 0 : const casacore::Bool showProcess = false;
1706 0 : casacore::Bool varyContours = false; //this is always false now...
1707 :
1708 : // Make a local decomposer
1709 :
1710 0 : ImageDecomposer<T> thresholdMap(*itsImagePtr);
1711 :
1712 :
1713 : // Generate global contours
1714 :
1715 0 : casacore::Vector<T> mainContours(itsNContour);
1716 0 : if (!varyContours) {
1717 0 : mainContours = autoContour(itsNContour, itsThresholdVal);
1718 : //cerr << "Main contours = " << mainContours << endl;
1719 : if (showProcess) displayContourMap(mainContours);
1720 : }
1721 :
1722 : // Find contiguous regions via thresholding
1723 :
1724 0 : casacore::uInt nRegions = thresholdMap.identifyRegions(itsThresholdVal, itsNAxis);
1725 0 : casacore::uInt nBadRegions = 0;
1726 0 : if (itsMinRange > 1)
1727 : {
1728 : // Eliminate weak regions
1729 0 : casacore::Vector<T> maxvals;
1730 0 : maxvals = thresholdMap.findAllRegionGlobalMax();
1731 0 : casacore::Vector<casacore::Bool> killRegion(nRegions,0);
1732 0 : for (casacore::uInt r = 0; r < nRegions; r++) {
1733 0 : if (thresholdMap.getContourVal(maxvals(r),mainContours)+1 <itsMinRange)
1734 : {
1735 0 : killRegion(r) = 1;
1736 0 : nBadRegions++;
1737 : }
1738 : }
1739 0 : thresholdMap.destroyRegions(killRegion);
1740 0 : }
1741 :
1742 0 : nRegions -= nBadRegions;
1743 :
1744 : if (showProcess) {
1745 : cout << "Located " << nRegions << " regions above threshold "
1746 : << itsThresholdVal << "." << endl;
1747 : thresholdMap.display();
1748 : }
1749 :
1750 :
1751 : // Find a blc and a trc for each region
1752 :
1753 0 : casacore::Block<casacore::IPosition> blc(nRegions);
1754 0 : casacore::Block<casacore::IPosition> trc(nRegions);
1755 0 : thresholdMap.boundRegions(blc, trc);
1756 0 : if (isDerived()) zero();
1757 :
1758 : if (showProcess) {
1759 : cout << "Bounded " << nRegions<<" regions for decomposition and fitting:"
1760 : << endl;
1761 : for (casacore::uInt r = 0; r < nRegions; r++) {
1762 : cout << r+1 <<": " << blc[r] << trc[r] << endl;
1763 : }
1764 : }
1765 :
1766 : // For each distinct region above the threshold, form a componentmap
1767 : // (subcomponentmap) and perform the fitting. Regions are treated as
1768 : // independent entities by the fitter - even if one region contains a
1769 : // very high-sigma component that may extend into another region, this
1770 : // is not taken into account by the other region
1771 :
1772 0 : for (casacore::uInt r=0; r<nRegions; r++) {
1773 :
1774 : // Make a decomposer for this region
1775 :
1776 0 : casacore::Slicer sl(blc[r], trc[r]-blc[r], casacore::Slicer::endIsLength);
1777 0 : casacore::SubImage<T> subIm(*itsImagePtr, sl);
1778 0 : ImageDecomposer<T> subpmap(subIm);
1779 0 : subpmap.copyOptions(*this);
1780 :
1781 : // Flag pixels outside the target region (this makes sure that other
1782 : // regions that happen to overlap part of the target region's bounding
1783 : // rectangle are not counted twice, and that only the target region
1784 : // pixels are used in fitting.)
1785 : {
1786 0 : casacore::IPosition pos(subpmap.itsDim,0); decrement(pos);
1787 0 : while (increment(pos,subpmap.shape())) {
1788 0 : if (thresholdMap.getCell(blc[r] + pos) != casacore::Int(r+1)) {
1789 0 : subpmap.setCell(pos, MASKED);
1790 : }
1791 : }
1792 0 : }
1793 :
1794 0 : casacore::Vector<T> subContours(itsNContour);
1795 : casacore::Vector<T> *contourPtr;
1796 :
1797 : // Generate contours for this region or use global
1798 : // ones for entire image
1799 :
1800 0 : if (varyContours) {
1801 0 : subContours = subpmap.autoContour(itsNContour, itsThresholdVal);
1802 0 : contourPtr = &subContours;
1803 : } else {
1804 0 : contourPtr = &mainContours;
1805 : }
1806 :
1807 : if (showProcess) {
1808 : cout << "-----------------------------------------" << endl;
1809 : cout << "Subimage " << r << endl;
1810 : cout << "Contour Map:" << endl;
1811 : subpmap.displayContourMap(*contourPtr);
1812 : }
1813 :
1814 : // Deblend components in this region
1815 :
1816 0 : subpmap.deblendRegions(*contourPtr, itsMinRange, itsNAxis);
1817 : if (showProcess) {
1818 : cout << "Component Map:" << endl;
1819 : subpmap.display();
1820 : }
1821 :
1822 0 : if (itsFitIt) {
1823 : // Fit gaussians to region
1824 0 : subpmap.fitComponents();
1825 : }
1826 : else {
1827 : // Just estimate
1828 0 : subpmap.itsNComponents = subpmap.itsNRegions;
1829 0 : subpmap.itsList.resize();
1830 0 : subpmap.itsList = subpmap.estimateComponents();
1831 : }
1832 :
1833 : if (showProcess) {
1834 : cout << "Object " << r+1 << " subcomponents: " << endl;
1835 : subpmap.printComponents();
1836 : cout << endl;
1837 : }
1838 :
1839 : // Add this region back into the main component map
1840 :
1841 0 : synthesize(subpmap, blc[r]);
1842 : }
1843 0 : }
1844 0 : return;
1845 : }
1846 :
1847 :
1848 :
1849 : // Fitting functions
1850 : // -----------------
1851 :
1852 : template <class T>
1853 0 : casacore::Matrix<T> ImageDecomposer<T>::fitRegion(casacore::Int nregion)
1854 : {
1855 0 : cout << "Fit Region " << nregion << endl;
1856 :
1857 : // Fits multiple gaussians to a single region. First performs a local
1858 : // maximum scan to estimate the number of components in the region.
1859 :
1860 0 : casacore::uInt nGaussians = 0;
1861 0 : casacore::uInt npoints = 0;
1862 0 : casacore::uInt ngpar = 0;
1863 0 : if (itsDim == 2) ngpar = 6;
1864 0 : if (itsDim == 3) ngpar = 9;
1865 0 : if (!isDerived()) nregion = 0; //fit to all data.
1866 :
1867 : // Determine number of data points in the region
1868 :
1869 : {
1870 0 : casacore::IPosition pos(itsDim,0);
1871 0 : decrement(pos);
1872 0 : while (increment(pos,shape())) {
1873 0 : if (getCell(pos) == nregion) npoints++;
1874 : }
1875 0 : }
1876 :
1877 : // Fill data and positions arrays
1878 :
1879 0 : casacore::Matrix<T> positions(npoints,itsDim);
1880 0 : casacore::Vector<T> dataValues(npoints);
1881 0 : casacore::Vector<T> sigma(npoints);
1882 0 : sigma = 1.0;
1883 0 : casacore::uInt k = 0;
1884 : {
1885 0 : casacore::IPosition pos(itsDim,0);
1886 0 : decrement(pos);
1887 0 : while (increment(pos,shape())) {
1888 0 : if (getCell(pos) == nregion) {
1889 0 : for (casacore::uInt i = 0; i<itsDim; i++) {
1890 0 : positions(k,i) = T(pos(i));
1891 : }
1892 0 : dataValues(k) = getImageVal(pos);
1893 0 : k++;
1894 : }
1895 : }
1896 0 : }
1897 :
1898 : // Estimate the initial parameters.
1899 :
1900 0 : casacore::Matrix<T> width;
1901 0 : casacore::Vector<T> maxval;
1902 0 : casacore::Block<casacore::IPosition> maxvalpos;
1903 :
1904 : // This estimates whether there are multiple components
1905 : // in the current region or not.
1906 :
1907 0 : findRegionLocalMax(maxval, maxvalpos, nregion, 2);
1908 0 : estimateComponentWidths(width, maxvalpos);
1909 0 : nGaussians = maxval.nelements();
1910 0 : cout << "Found " << nGaussians << " components" << endl;
1911 :
1912 0 : casacore::Matrix<T> initestimate(nGaussians, ngpar);
1913 0 : casacore::Matrix<T> solution(nGaussians, ngpar);
1914 :
1915 0 : if (itsDim == 2) {
1916 0 : for (casacore::uInt r = 0; r < nGaussians; r++) {
1917 0 : initestimate(r,0) = maxval(r);
1918 0 : initestimate(r,1) = maxvalpos[r](0);
1919 0 : initestimate(r,2) = maxvalpos[r](1);
1920 0 : initestimate(r,3) = width(r,1);
1921 0 : initestimate(r,4) = width(r,0)/width(r,1);
1922 0 : initestimate(r,5) = 0;
1923 : }
1924 : }
1925 0 : if (itsDim == 3) {
1926 0 : for (casacore::uInt r = 0; r < nGaussians; r++) {
1927 0 : initestimate(r,0) = maxval(r);
1928 0 : initestimate(r,1) = maxvalpos[r](0);
1929 0 : initestimate(r,2) = maxvalpos[r](1);
1930 0 : initestimate(r,3) = maxvalpos[r](2);
1931 0 : initestimate(r,4) = width(r,0);
1932 0 : initestimate(r,5) = width(r,1);
1933 0 : initestimate(r,6) = width(r,2);
1934 0 : initestimate(r,7) = (0.0);
1935 0 : initestimate(r,8) = (0.0);
1936 : }
1937 : }
1938 :
1939 : // Fit for nGaussians simultaneously
1940 :
1941 0 : solution = fitGauss(positions, dataValues, initestimate);
1942 0 : return solution;
1943 :
1944 0 : }
1945 :
1946 : template <class T>
1947 0 : void ImageDecomposer<T>::fitRegions()
1948 : {
1949 : // Fits gaussians to an image; multiple gaussians per region in the pmap.
1950 : // The regions are fit sequentially and independently, so this function
1951 : // can be used on the main image.
1952 : // If the map is not yet thresholded, will fit to the entire image as if it
1953 : // were a single composite object, which will be very slow.
1954 :
1955 0 : casacore::uInt ngpar = 0;
1956 0 : if (itsDim == 2) ngpar = 6;
1957 0 : if (itsDim == 3) ngpar = 9;
1958 :
1959 0 : if (itsNRegions == 0) { //not deblended.
1960 0 : itsList = fitRegion(0);
1961 0 : return;
1962 : }
1963 : //
1964 0 : for (casacore::uInt r = 1; r <= itsNRegions; r++) {
1965 0 : casacore::Matrix<T> subitsList;
1966 0 : casacore::Matrix<T> olditsList;
1967 0 : subitsList = fitRegion(r);
1968 0 : olditsList = itsList;
1969 0 : itsList.resize(itsNComponents + subitsList.nrow(), ngpar);
1970 : //
1971 0 : for (casacore::uInt c = 0; c < itsNComponents; c++) {
1972 0 : for (casacore::uInt p = 0; p < ngpar; p++) {
1973 0 : itsList(c,p) = olditsList(c,p);
1974 : }
1975 : }
1976 :
1977 0 : for (casacore::uInt subc = 0; subc < subitsList.nrow(); subc++) {
1978 0 : for (casacore::uInt p = 0; p < ngpar; p++) {
1979 0 : itsList(itsNComponents+subc, p) = subitsList(subc, p);
1980 : }
1981 : }
1982 0 : itsNComponents += subitsList.nrow();
1983 : }
1984 :
1985 0 : return;
1986 : }
1987 :
1988 :
1989 : template <class T>
1990 0 : void ImageDecomposer<T>::fitComponents()
1991 : {
1992 : // Fits gaussians to an image; one gaussian per region in the pmap.
1993 : // This function is intended to be used only by ImageDecomposer on its
1994 : // intermediary subimages; using it at higher level will execute a full
1995 : // gaussian fit on the main image and will be extremely slow. Every
1996 : // nonflagged object pixel in the image is used in fitting.
1997 :
1998 : // If the deblended flag is true, the function will treat each region as
1999 : // an individual component and will fit that many gaussians to the image
2000 :
2001 0 : casacore::uInt ngpar = itsDim * 3;
2002 0 : casacore::uInt npoints = 0;
2003 :
2004 0 : if (!isDerived()) {
2005 0 : throw(casacore::AipsError("Cannot fit until components are deblended"
2006 : " - use identifyRegions() or deblendRegions()"));
2007 : }
2008 :
2009 :
2010 : // Determine number of data points in all the regions
2011 : // and get data and position vectors
2012 :
2013 : {
2014 0 : casacore::IPosition pos(itsDim,0);
2015 0 : decrement(pos);
2016 0 : while (increment(pos,shape())) {
2017 0 : if (getCell(pos) > 0) npoints++;
2018 : }
2019 0 : }
2020 :
2021 0 : casacore::Matrix<T> positions(npoints,itsDim);
2022 0 : casacore::Vector<T> dataValues(npoints);
2023 :
2024 : {
2025 0 : casacore::IPosition pos(itsDim,0);
2026 0 : decrement(pos);
2027 0 : casacore::uInt p = 0;
2028 0 : while (increment(pos,shape())) {
2029 0 : if (getCell(pos) > 0) {
2030 0 : for (casacore::uInt i=0; i<itsDim; i++) {
2031 0 : positions(p,i) = T(pos(i));
2032 : }
2033 0 : dataValues(p) = getImageVal(pos);
2034 0 : p++;
2035 : }
2036 : }
2037 0 : }
2038 :
2039 : // Estimate the initial parameters.
2040 :
2041 0 : casacore::Matrix<T> initestimate(itsNRegions, ngpar);
2042 0 : casacore::Matrix<T> solution(itsNRegions, ngpar);
2043 :
2044 0 : initestimate = estimateComponents();
2045 :
2046 0 : solution = fitGauss(positions, dataValues, initestimate);
2047 :
2048 0 : itsNComponents = itsNRegions;
2049 0 : itsList.resize(solution.shape());
2050 0 : itsList = solution;
2051 :
2052 0 : return;
2053 0 : }
2054 :
2055 :
2056 : template <class T>
2057 0 : casacore::Matrix<T> ImageDecomposer<T>::estimateComponents()
2058 : {
2059 0 : casacore::uInt ngaussians = itsNRegions;
2060 0 : casacore::uInt ngpar = 0;
2061 0 : if (itsDim == 2) ngpar = 6;
2062 0 : if (itsDim == 3) ngpar = 9;
2063 0 : if (!isDerived()) {
2064 0 : throw(casacore::AipsError("Cannot estimate until components are deblended"
2065 : " - use identifyRegions() or deblendRegions()"));
2066 : }
2067 :
2068 :
2069 0 : casacore::Matrix<T> estimate(ngaussians, ngpar);
2070 0 : casacore::Matrix<T> width;
2071 0 : casacore::Vector<T> maxval;
2072 0 : casacore::Block<casacore::IPosition> maxvalpos;
2073 :
2074 0 : ngaussians = itsNRegions;
2075 0 : findAllRegionGlobalMax(maxval, maxvalpos);
2076 :
2077 : // This is based on the moment analysis methods given in Jarvis & Tyson 1981,
2078 : // though they have been generalized to 3D. The formula to estimate phi
2079 : // (3D parameter 8) was not derived rigorously and may not be correct, though
2080 : // it works very well overall.
2081 :
2082 0 : for (casacore::uInt r = 0; r < ngaussians; r++)
2083 : {
2084 0 : if (itsDim == 2) {
2085 0 : casacore::Matrix<T> M;
2086 0 : M = calculateMoments(r+1);
2087 0 : estimate(r,0) = maxval[r];
2088 0 : estimate(r,1) = M(1,0) / M(0,0);
2089 0 : estimate(r,2) = M(0,1) / M(0,0);
2090 0 : estimate(r,3) = 2.84 * sqrt(M(0,2)/M(0,0));
2091 0 : estimate(r,4) = sqrt(M(2,0)/M(0,2));
2092 0 : estimate(r,5) = 0.5 * atan(2 * M(1,1) / (M(2,0) - M(0,2)));
2093 0 : if (estimate(r,3) < 1.0) estimate(r,3) = 1.0;
2094 0 : if (estimate(r,4)*estimate(r,3) < 1.0) estimate(r,4) = 1.0/estimate(r,3);
2095 0 : if (casacore::isNaN(estimate(r,4))) estimate(r,4) = 1.0/estimate(r,3);
2096 0 : if (casacore::isNaN(estimate(r,5))) estimate(r,5) = 0;
2097 0 : } else if (itsDim == 3) {
2098 0 : casacore::Cube<T> M;
2099 0 : M = calculateMoments(r+1);
2100 0 : estimate(r,0) = maxval[r];
2101 0 : estimate(r,1) = M(1,0,0) / M(0,0,0);
2102 0 : estimate(r,2) = M(0,1,0) / M(0,0,0);
2103 0 : estimate(r,3) = M(0,0,1) / M(0,0,0);
2104 0 : estimate(r,4) = 2.84 * sqrt(M(2,0,0) / M(0,0,0));
2105 0 : estimate(r,5) = 2.84 * sqrt(M(0,2,0) / M(0,0,0));
2106 0 : estimate(r,6) = 2.84 * sqrt(M(0,0,2) / M(0,0,0));
2107 0 : estimate(r,7) = 0.5 * atan(2 * M(1,1,0) / (M(2,0,0) - M(0,2,0)));
2108 0 : estimate(r,8) = -0.5 * atan(2 * M(1,0,1) /
2109 0 : ((M(2,0,0) - M(0,0,2))*cos(estimate(r,7)) +
2110 0 : (M(0,2,0) - M(0,0,2))*sin(estimate(r,7))));
2111 0 : if (estimate(r,4) < 1.0) estimate(r,4) = 1.0;
2112 0 : if (estimate(r,5) < 1.0) estimate(r,5) = 1.0;
2113 0 : if (estimate(r,6) < 1.0) estimate(r,6) = 1.0;
2114 0 : if (casacore::isNaN(estimate(r,8))) estimate(r,8) = 0;
2115 0 : }
2116 : }
2117 :
2118 0 : return estimate;
2119 0 : }
2120 :
2121 :
2122 :
2123 :
2124 :
2125 :
2126 :
2127 : template <class T>
2128 0 : casacore::Matrix<T> ImageDecomposer<T>::fitGauss(const casacore::Matrix<T>& positions,
2129 : const casacore::Vector<T>& dataValues,
2130 : const casacore::Matrix<T>& initestimate) const
2131 : {
2132 : // Fits the specified number of 3D gaussians to the data, and returns
2133 : // solution in image (world) coordinates.
2134 :
2135 : //casacore::uInt ngpar = 0;
2136 0 : casacore::uInt ngaussians = initestimate.nrow();
2137 : //if (itsDim == 2) ngpar = 6;
2138 : //if (itsDim == 3) ngpar = 9;
2139 :
2140 0 : casacore::Matrix<T> solution;
2141 : //T chisquare;
2142 :
2143 : // Might be useful to send to screen in AIPS++
2144 : //cout << "Primary estimation matrix:" << endl;
2145 : //for (casacore::uInt r = 0; r < ngaussians; r++)
2146 : // cout << initestimate.row(r) << endl;
2147 :
2148 0 : casacore::FitGaussian<T> fitter(itsDim,ngaussians);
2149 0 : fitter.setFirstEstimate(initestimate);
2150 0 : if (itsMaxRetries < 0) {
2151 0 : fitter.setMaxRetries(itsDim * ngaussians);
2152 : }
2153 : else {
2154 0 : fitter.setMaxRetries(itsMaxRetries);
2155 : }
2156 :
2157 : try{
2158 0 : solution = fitter.fit(positions, dataValues, itsMaximumRMS, itsMaxIter,
2159 0 : itsConvCriteria);
2160 0 : } catch (casacore::AipsError fiterr) {
2161 0 : cout << fiterr.getMesg() << endl;
2162 0 : cout << "Fitting failed." << endl;
2163 0 : solution = 0;
2164 : //chisquare = -1.0;
2165 0 : return solution;
2166 : }
2167 :
2168 0 : if (fitter.converged()) {
2169 : //chisquare = fitter.chisquared();
2170 : } else {
2171 0 : cout << "Fitting did not converge to a reasonable result - using estimate."
2172 0 : << endl;
2173 0 : solution = initestimate;
2174 : //chisquare = -1.0;
2175 0 : return solution;
2176 : }
2177 :
2178 0 : return solution;
2179 0 : }
2180 :
2181 :
2182 :
2183 :
2184 :
2185 : template <class T>
2186 : void ImageDecomposer<T>::display() const
2187 : {
2188 : // Displays the componentmap in a terminal environment as an array of
2189 : // characters on screen.
2190 :
2191 : casacore::Int windowwidth = 80;
2192 : casacore::Int const spacing = 4;
2193 : casacore::Int const cposinc = shape(0) * 2 + spacing;
2194 : if (cposinc > windowwidth) windowwidth = cposinc;
2195 : casacore::Int cpos = 0;
2196 : //
2197 : casacore::Int z = 0;
2198 : casacore::Int benchz = 0;
2199 : //
2200 : //cerr << "shape = " << shape() << endl;
2201 : while (z < shape(2)) {
2202 : for (casacore::Int y = 0; y < shape(1); y++) {
2203 : z = benchz;
2204 : while ((cpos += cposinc) < windowwidth && z < shape(2)) {
2205 : for (casacore::Int x = 0; x < itsShape(0); x++) {
2206 : if (getCell(x,y,z) >= 0) {
2207 : cout << getCell(x,y,z);
2208 : } else if (getCell(x,y,z) == INDETERMINATE) {
2209 : cout << '*';
2210 : } else if (getCell(x,y,z) == MASKED) {
2211 : cout << '-';
2212 : }
2213 : if (getCell(x,y,z) < 10) {
2214 : cout << ' ';
2215 : }
2216 : }
2217 : z++;
2218 : cout << " ";
2219 : }
2220 : cpos = 0;
2221 : cout << endl;
2222 : }
2223 : benchz = z;
2224 : cout << endl;
2225 : }
2226 : return;
2227 : }
2228 :
2229 :
2230 : // Console display functions
2231 : // -------------------------
2232 :
2233 :
2234 : template <class T>
2235 : void ImageDecomposer<T>::displayContourMap(const casacore::Vector<T>& clevels) const
2236 : {
2237 : // Displays the target image as a contourmap in a terminal environment as
2238 : // an array of characters on screen.
2239 :
2240 : casacore::Int windowwidth = 80;
2241 : casacore::Int const spacing = 4;
2242 : casacore::Int const cposinc = shape(0) * 2 + spacing;
2243 : if (cposinc > windowwidth) windowwidth = cposinc;
2244 : casacore::Int cpos = 0;
2245 :
2246 : casacore::Int z = 0;
2247 : casacore::Int benchz = 0;
2248 : cout << "Contour levels:" << clevels << endl;
2249 :
2250 : if (itsDim == 2) {
2251 : for (casacore::Int y = 0; y < shape(1); y++) {
2252 : for (casacore::Int x = 0; x < shape(0); x++) {
2253 : if (getContourVal(x,y,clevels) >= 0) {
2254 : cout << getContourVal(x,y,clevels);
2255 : }
2256 : else if (getContourVal(x,y,clevels) <= -1) {
2257 : cout << '-';
2258 : }
2259 : if (getContourVal(x,y,clevels) < 10) {
2260 : cout << ' ';
2261 : }
2262 : }
2263 : cout << endl;
2264 : }
2265 : cout << endl;
2266 : }
2267 :
2268 : if (itsDim == 3){
2269 : //this actually works in 2 dimensions on a casacore::TempImage, but not on a
2270 : //casacore::SubImage, where there is a failure inside the getImageVal command
2271 : //on a failed assertion involving a casacore::LatticeRegion object. As a result
2272 : //the above specialization was written, but it would be nice if 3-D
2273 : //IPositions worked on 2-D images in casacore::SubImage as well as TempImage.
2274 : while (z < shape(2)) {
2275 : for (casacore::Int y = 0; y < shape(1); y++) {
2276 : z = benchz;
2277 : while ((cpos += cposinc) < windowwidth && (z < shape(2))) {
2278 : for (casacore::Int x = 0; x < shape(0); x++) {
2279 : if (getContourVal(x,y,z,clevels) >= 0) {
2280 : cout << getContourVal(x,y,z,clevels);
2281 : }
2282 : else if (getContourVal(x,y,z,clevels) <= -1) {
2283 : cout << '-';
2284 : }
2285 : if (getContourVal(x,y,z,clevels) < 10) {
2286 : cout << ' ';
2287 : }
2288 : }
2289 : z++;
2290 : cout << " ";
2291 : }
2292 : cpos = 0;
2293 : cout << endl;
2294 : }
2295 : benchz = z;
2296 : cout << endl;
2297 : }
2298 : cout << endl;
2299 : }
2300 :
2301 : return;
2302 : }
2303 :
2304 :
2305 : template <class T>
2306 0 : void ImageDecomposer<T>::printComponents() const
2307 : {
2308 : //Prints the components as formatted output on screen.
2309 : //IMPR: Probably could be modified as an ostream output function.
2310 :
2311 0 : casacore::Matrix<T> clist;
2312 0 : clist = componentList();
2313 :
2314 0 : for (casacore::uInt g = 0; g < clist.nrow(); g++)
2315 : {
2316 0 : cout << g+1 << ": ";
2317 0 : if (itsList(g,0) == T(0)) {
2318 0 : cout << "Failed";
2319 : } else {
2320 0 : cout << "Peak: " << clist(g,0) << " ";
2321 :
2322 0 : if (itsDim == 2) {
2323 0 : cout << "Mu: [" << clist(g,1)
2324 0 : << ", " << clist(g,2) << "] ";
2325 0 : cout << "Axes: [" << clist(g,3)
2326 0 : << ", " << clist(g,3) * clist(g,4) << "] ";
2327 0 : cout << "Rotation: " << clist(g,5) /* * 180 / C::pi */;
2328 : }
2329 0 : if (itsDim == 3) {
2330 0 : cout << "Mu: [" << clist(g,1)
2331 0 : << ", " << clist(g,2)
2332 0 : << ", " << clist(g,3) << "] ";
2333 0 : cout << "Axes: [" << clist(g,4)
2334 0 : << ", " << clist(g,5)
2335 0 : << ", " << clist(g,6) << "] ";
2336 0 : cout << "Rotation: [" << clist(g,7)/* *180/C::pi */
2337 0 : << ", " << clist(g,8) /* *180/C::pi */ << "]";
2338 : }
2339 : }
2340 0 : cout << endl;
2341 : }
2342 :
2343 0 : return;
2344 0 : }
2345 :
2346 :
2347 : // Functions associated only with other classes
2348 : // --------------------------------------------
2349 :
2350 : template <class T>
2351 0 : void ImageDecomposer<T>::correctBlcTrc(casacore::IPosition& blc, casacore::IPosition& trc) const
2352 : {
2353 : //Ensures blc and trc correctly describe the limits of a rectangular block.
2354 :
2355 : casacore::Int t;
2356 0 : for (casacore::uInt i = 0; i<itsDim; i++) {
2357 0 : if (blc(i)<0) blc(i) = 0;
2358 0 : if (trc(i)>shape(i)) trc(i) = shape(i);
2359 0 : if (trc(i)<0) trc(i) = 0;
2360 : //
2361 0 : if (blc(i)>shape(i)) blc(i) = shape(i);
2362 0 : if (blc(i)>trc(i)) { // nebk - ERROR should be trc(a) not trc(0) ??
2363 0 : t = blc(i);
2364 0 : blc(i) = trc(i);
2365 0 : trc(i) = t;
2366 : }
2367 : }
2368 0 : return;
2369 : }
2370 :
2371 : template <class T>
2372 0 : inline casacore::Bool ImageDecomposer<T>::increment(casacore::IPosition& pos,
2373 : const casacore::IPosition& limit) const
2374 : {
2375 : // N-Dimensional looping function: use in place of nested for loops
2376 : // Returns false when pos reaches limit.
2377 : // Use as follows: while(increment(pos,limit))
2378 : // IMPR: this function probably should be global or in IPosition. Or even
2379 : // better, omitted completely by using LatticeIterators.
2380 :
2381 0 : pos(itsDim-1)++;
2382 0 : for (casacore::uInt i = itsDim-1; i>0; i--) {
2383 0 : if (pos(i) == limit(i)) {
2384 0 : pos(i) = 0;
2385 0 : pos(i-1)++;
2386 : } else {
2387 0 : return true;
2388 : }
2389 : }
2390 : //
2391 0 : if (pos(0) == limit(0)) {
2392 0 : return false;
2393 : } else {
2394 0 : return true;
2395 : }
2396 : }
2397 :
2398 : template <class T>
2399 0 : inline void ImageDecomposer<T>::decrement(casacore::IPosition& pos) const
2400 : {
2401 : //To ensure while loop starts at 0,0,0..., decrement before first increment
2402 0 : pos(itsDim-1)--;
2403 0 : }
2404 :
2405 :
2406 :
2407 : /*
2408 : casacore::RO_LatticeIterator<casacore::Int> otherIter(*(other.itsMapPtr));
2409 : casacore::Bool deleteOther, deleteThis;
2410 : const casacore::Int* pOther = 0;
2411 : casacore::Int* pThis = 0;
2412 : for (otherIter.reset(); !otherIter.atEnd(); otherIter++) {
2413 : const casacore::Array<casacore::Int>& otherArray = otherIter.cursor();
2414 : casacore::Array<casacore::Int> thisArray = itsMapPtr->getSlice(otherIter.position(), otherIter.cursorShape());
2415 : //
2416 : pOther = otherArray.getStorage(deleteOther);
2417 : pThis = thisArray.getStorage(deleteThis);
2418 : //
2419 : for (casacore::uInt i=0; i<otherIter.cursor().nelements(); i++) {
2420 : if (pOther[i] != regionID) pThis[i] = MASKED;
2421 : }
2422 : //
2423 : otherArray.freeStorage(pOther, deleteOther);
2424 : thisArray.putStorage(pThis, deleteThis);
2425 : }
2426 : */
2427 :
2428 :
2429 : } //# NAMESPACE CASA - END
2430 :
|
