Line data Source code
1 : #include <imageanalysis/ImageAnalysis/StatImageCreator.h>
2 :
3 : #include <imageanalysis/Annotations/AnnCenterBox.h>
4 : #include <imageanalysis/Annotations/AnnCircle.h>
5 : #include <casacore/lattices/LEL/LatticeExpr.h>
6 :
7 : namespace casa {
8 :
9 0 : StatImageCreator::StatImageCreator(
10 : const SPCIIF image, const Record *const region,
11 : const String& mask, const String& outname, Bool overwrite
12 0 : ) : ImageStatsBase(image, region, mask, outname, overwrite) {
13 0 : this->_construct();
14 0 : auto da = _getImage()->coordinates().directionAxesNumbers();
15 0 : _dirAxes[0] = da[0];
16 0 : _dirAxes[1] = da[1];
17 0 : setAnchorPosition(0, 0);
18 0 : }
19 :
20 0 : void StatImageCreator::setRadius(const Quantity& radius) {
21 0 : ThrowIf(
22 : ! (radius.getUnit().startsWith("pix") || radius.isConform("rad")),
23 : "radius units " + radius.getUnit() + " must be pixel or angular"
24 : );
25 0 : _xlen = radius;
26 0 : _ylen = Quantity(0,"");
27 0 : }
28 :
29 0 : void StatImageCreator::setRectangle(
30 : const Quantity& xLength, const Quantity& yLength
31 : ) {
32 0 : ThrowIf(
33 : ! (xLength.getUnit().startsWith("pix") || xLength.isConform("rad")),
34 : "xLength units " + xLength.getUnit() + " must be pixel or angular"
35 : );
36 0 : ThrowIf(
37 : ! (yLength.getUnit().startsWith("pix") || yLength.isConform("rad")),
38 : "xLength units " + yLength.getUnit() + " must be pixel or angular"
39 : );
40 0 : _xlen = xLength;
41 0 : _ylen = yLength;
42 0 : }
43 :
44 0 : void StatImageCreator::setAnchorPosition(Int x, Int y) {
45 0 : Vector<Double> anchorPixel(_getImage()->ndim(), 0);
46 0 : anchorPixel[_dirAxes[0]] = x;
47 0 : anchorPixel[_dirAxes[1]] = y;
48 0 : _anchor = _getImage()->coordinates().toWorld(anchorPixel);
49 0 : }
50 :
51 0 : void StatImageCreator::useReferencePixelAsAnchor() {
52 0 : const auto refPix = _getImage()->coordinates().referencePixel();
53 0 : Int x = round(refPix[_dirAxes[0]]);
54 0 : Int y = round(refPix[_dirAxes[1]]);
55 0 : *_getLog() << LogIO::NORMAL << LogOrigin("StatImageCreator", __func__)
56 : << "Anchor being set at pixel [" << x << "," << y
57 0 : << "], at/near image reference pixel." << LogIO::POST;
58 0 : setAnchorPosition(x, y);
59 0 : }
60 :
61 0 : void StatImageCreator::setGridSpacing(uInt x, uInt y) {
62 0 : _grid.first = x;
63 0 : _grid.second = y;
64 0 : }
65 :
66 0 : SPIIF StatImageCreator::compute() {
67 0 : static const AxesSpecifier dummyAxesSpec;
68 0 : *_getLog() << LogOrigin(getClass(), __func__);
69 0 : auto mylog = _getLog().get();
70 : auto subImageRO = SubImageFactory<Float>::createSubImageRO(
71 0 : *_getImage(), *_getRegion(), _getMask(),
72 0 : mylog, dummyAxesSpec, _getStretch()
73 0 : );
74 : auto subImageRW = SubImageFactory<Float>::createImage(
75 0 : *_getImage(), "", *_getRegion(), _getMask(),
76 0 : dummyAxesSpec, False, False, _getStretch()
77 0 : );
78 0 : _doMask = ! ImageMask::isAllMaskTrue(*subImageRO);
79 0 : const auto imshape = subImageRO->shape();
80 0 : const auto xshape = imshape[_dirAxes[0]];
81 0 : const auto yshape = imshape[_dirAxes[1]];
82 0 : const auto& csys = subImageRO->coordinates();
83 0 : auto anchorPixel = csys.toPixel(_anchor);
84 0 : Int xanchor = rint(anchorPixel[_dirAxes[0]]);
85 0 : Int yanchor = rint(anchorPixel[_dirAxes[1]]);
86 : // ensure xanchor and yanchor are positive
87 0 : if (xanchor < 0) {
88 : // ugh, mod of a negative number in C++ doesn't do what I want it to
89 : // integer division
90 0 : xanchor += (abs(xanchor)/_grid.first + 1)*_grid.first;
91 : }
92 0 : if (yanchor < 0) {
93 0 : yanchor += (abs(yanchor)/_grid.second + 1)*_grid.second;
94 : }
95 : // xstart and ystart are the pixel location in the
96 : // subimage of the lower left corner of the grid,
97 : // ie they are the pixel location of the grid point
98 : // with the smallest non-negative x and y values
99 0 : Int xstart = xanchor % _grid.first;
100 0 : Int ystart = yanchor % _grid.second;
101 0 : if (xstart < 0) {
102 0 : xstart += _grid.first;
103 : }
104 0 : if (ystart < 0) {
105 0 : ystart += _grid.second;
106 : }
107 0 : Array<Float> tmp;
108 : // integer division
109 0 : auto nxpts = (xshape - xstart)/_grid.first;
110 0 : if ((xshape - xstart) % _grid.first != 0) {
111 0 : ++nxpts;
112 : }
113 0 : auto nypts = (yshape - ystart)/ _grid.second;
114 0 : if ((yshape - ystart) % _grid.second != 0) {
115 0 : ++nypts;
116 : }
117 0 : IPosition storeShape = imshape;
118 0 : storeShape[_dirAxes[0]] = nxpts;
119 0 : storeShape[_dirAxes[1]] = nypts;
120 0 : auto interpolate = _grid.first > 1 || _grid.second > 1;
121 0 : TempImage<Float>* writeTo = dynamic_cast<TempImage<Float> *>(subImageRW.get());
122 0 : std::unique_ptr<TempImage<Float>> store;
123 0 : if (interpolate) {
124 0 : store.reset(new TempImage<Float>(storeShape, csys));
125 0 : store->set(0.0);
126 0 : if (_doMask) {
127 0 : store->attachMask(ArrayLattice<Bool>(storeShape));
128 0 : store->pixelMask().set(True);
129 : }
130 0 : writeTo = store.get();
131 : }
132 0 : _computeStat(*writeTo, subImageRO, nxpts, nypts, xstart, ystart);
133 0 : if (_doProbit) {
134 0 : writeTo->copyData((LatticeExpr<Float>)(*writeTo * C::probit_3_4));
135 : }
136 0 : if (interpolate) {
137 0 : _doInterpolation(
138 0 : subImageRW, *store, subImageRO,
139 : nxpts, nypts, xstart, ystart
140 : );
141 : }
142 0 : auto res = _prepareOutputImage(*subImageRW);
143 0 : return res;
144 0 : }
145 :
146 0 : void StatImageCreator::_computeStat(
147 : TempImage<Float>& writeTo, SPCIIF subImage,
148 : uInt nxpts, uInt nypts, Int xstart, Int ystart
149 : ) {
150 0 : std::shared_ptr<Array<Bool>> regionMask;
151 0 : uInt xBlcOff = 0;
152 0 : uInt yBlcOff = 0;
153 0 : uInt xChunkSize = 0;
154 0 : uInt yChunkSize = 0;
155 0 : _nominalChunkInfo(
156 : regionMask, xBlcOff, yBlcOff,
157 : xChunkSize, yChunkSize, subImage
158 : );
159 0 : Array<Bool> regMaskCopy;
160 0 : if (regionMask) {
161 0 : regMaskCopy = regionMask->copy();
162 0 : if (! writeTo.hasPixelMask()) {
163 0 : ArrayLattice<Bool> mymask(writeTo.shape());
164 0 : mymask.set(True);
165 0 : writeTo.attachMask(mymask);
166 0 : }
167 : }
168 0 : auto imshape = subImage->shape();
169 0 : auto ndim = imshape.size();
170 0 : IPosition chunkShape(ndim, 1);
171 0 : chunkShape[_dirAxes[0]] = xChunkSize;
172 0 : chunkShape[_dirAxes[1]] = yChunkSize;
173 0 : auto xsize = imshape[_dirAxes[0]];
174 0 : auto ysize = imshape[_dirAxes[1]];
175 0 : IPosition planeShape(imshape.size(), 1);
176 0 : planeShape[_dirAxes[0]] = xsize;
177 0 : planeShape[_dirAxes[1]] = ysize;
178 : RO_MaskedLatticeIterator<Float> lattIter (
179 0 : *subImage, planeShape, True
180 0 : );
181 0 : String algName;
182 0 : auto alg = _getStatsAlgorithm(algName);
183 0 : std::set<StatisticsData::STATS> statToCompute;
184 0 : statToCompute.insert(_statType);
185 0 : alg->setStatsToCalculate(statToCompute);
186 0 : auto ngrid = nxpts*nypts;
187 0 : auto nPts = ngrid;
188 0 : IPosition loopAxes;
189 0 : for (uInt i=0; i<ndim; ++i) {
190 0 : if (i != _dirAxes[0] && i != _dirAxes[1]) {
191 0 : loopAxes.append(IPosition(1, i));
192 0 : nPts *= imshape[i];
193 : }
194 : }
195 0 : auto doCircle = _ylen.getValue() <= 0;
196 0 : *_getLog() << LogOrigin(getClass(), __func__) << LogIO::NORMAL
197 0 : << "Using ";
198 0 : if (doCircle) {
199 0 : *_getLog() << "circular region of radius " << _xlen;
200 : }
201 : else {
202 0 : *_getLog() << "rectangular region of specified dimensions " << _xlen
203 0 : << " x " << _ylen;
204 : }
205 0 : *_getLog() << " (because of centering and "
206 : << "rounding to use whole pixels, actual dimensions of bounding box are "
207 0 : << xChunkSize << " pix x " << yChunkSize << " pix";
208 0 : if (doCircle) {
209 0 : *_getLog() << " and there are " << ntrue(*regionMask)
210 0 : << " good pixels in the circle that are being used";
211 : }
212 0 : *_getLog() << ") to choose pixels for computing " << _statName
213 : << " using the " << algName << " algorithm around each of "
214 0 : << ngrid << " grid points in " << (nPts/ngrid) << " planes." << LogIO::POST;
215 0 : auto imageChunkShape = chunkShape;
216 0 : _doStatsLoop(
217 : writeTo, lattIter, nxpts, nypts, xstart, ystart,
218 : xBlcOff, yBlcOff, xChunkSize, yChunkSize, imshape,
219 : chunkShape, regionMask, alg, regMaskCopy, loopAxes, nPts
220 : );
221 0 : }
222 :
223 0 : void StatImageCreator::_doStatsLoop(
224 : TempImage<Float>& writeTo, RO_MaskedLatticeIterator<Float>& lattIter,
225 : uInt nxpts, uInt nypts, Int xstart, Int ystart, uInt xBlcOff,
226 : uInt yBlcOff, uInt xChunkSize, uInt yChunkSize, const IPosition& imshape,
227 : const IPosition& chunkShape, std::shared_ptr<Array<Bool>> regionMask,
228 : std::shared_ptr<
229 : StatisticsAlgorithm<
230 : Double, Array<Float>::const_iterator, Array<Bool>::const_iterator,
231 : Array<Float>::const_iterator
232 : >
233 : >& alg, const Array<Bool>& regMaskCopy, const IPosition& loopAxes, uInt nPts
234 : ) {
235 0 : auto ximshape = imshape[_dirAxes[0]];
236 0 : auto yimshape = imshape[_dirAxes[1]];
237 0 : auto ndim = imshape.size();
238 0 : IPosition planeBlc(ndim, 0);
239 0 : auto& xPlaneBlc = planeBlc[_dirAxes[0]];
240 0 : auto& yPlaneBlc = planeBlc[_dirAxes[1]];
241 0 : auto planeTrc = planeBlc + chunkShape - 1;
242 0 : auto& xPlaneTrc = planeTrc[_dirAxes[0]];
243 0 : auto& yPlaneTrc = planeTrc[_dirAxes[1]];
244 0 : uInt nCount = 0;
245 0 : auto hasLoopAxes = ! loopAxes.empty();
246 0 : ProgressMeter pm(0, nPts, "Processing stats at grid points");
247 : // Vector rather than IPosition because blc values can be negative
248 0 : Vector<Int> blc(ndim, 0);
249 0 : blc[_dirAxes[0]] = xstart - xBlcOff;
250 0 : blc[_dirAxes[1]] = ystart - yBlcOff;
251 0 : for (lattIter.atStart(); ! lattIter.atEnd(); ++lattIter) {
252 0 : auto plane = lattIter.cursor();
253 0 : auto lattMask = lattIter.getMask();
254 0 : auto checkGrid = ! allTrue(lattMask);
255 0 : auto outPos = lattIter.position();
256 0 : auto& xOutPos = outPos[_dirAxes[0]];
257 0 : auto& yOutPos = outPos[_dirAxes[1]];
258 : // inPos is the position of the current grid point
259 : // This is the position in the current chunk, not the entire lattice
260 0 : auto inPos = plane.shape() - 1;
261 0 : inPos[_dirAxes[0]] = xstart;
262 0 : inPos[_dirAxes[1]] = ystart;
263 0 : auto& xInPos = inPos[_dirAxes[0]];
264 0 : auto& yInPos = inPos[_dirAxes[1]];
265 0 : Int yblc = ystart - yBlcOff;
266 0 : for (uInt yCount=0; yCount<nypts; ++yCount, yblc+=_grid.second, yInPos+=_grid.second) {
267 0 : yOutPos = yCount;
268 0 : yPlaneBlc = max(0, yblc);
269 0 : yPlaneTrc = min(
270 0 : yblc + (Int)yChunkSize - 1, (Int)yimshape - 1
271 : );
272 0 : IPosition regMaskStart(ndim, 0);
273 0 : auto& xRegMaskStart = regMaskStart[_dirAxes[0]];
274 0 : auto& yRegMaskStart = regMaskStart[_dirAxes[1]];
275 0 : auto regMaskLength = regMaskStart;
276 0 : auto& xRegMaskLength = regMaskLength[_dirAxes[0]];
277 0 : auto& yRegMaskLength = regMaskLength[_dirAxes[1]];
278 0 : Bool yDoMaskSlice = False;
279 0 : if (regionMask) {
280 0 : regMaskLength = regionMask->shape();
281 0 : if (yblc < 0) {
282 0 : yRegMaskStart = -yblc;
283 0 : yRegMaskLength += yblc;
284 0 : yDoMaskSlice = True;
285 : }
286 0 : else if (yblc + yChunkSize > yimshape) {
287 0 : yRegMaskLength = yimshape - yblc;
288 0 : yDoMaskSlice = True;
289 : }
290 : }
291 0 : xInPos = xstart;
292 0 : Int xblc = xstart - xBlcOff;
293 0 : for (uInt xCount=0; xCount<nxpts; ++xCount, xblc+=_grid.first, xInPos+=_grid.first) {
294 0 : Float res = 0;
295 0 : xOutPos = xCount;
296 0 : if (checkGrid && ! lattMask(inPos)) {
297 : // grid point is masked, no computation should be done
298 0 : writeTo.putAt(res, outPos);
299 0 : writeTo.pixelMask().putAt(False, outPos);
300 : }
301 : else {
302 0 : xPlaneBlc = max(0, xblc);
303 0 : xPlaneTrc = min(
304 0 : xblc + (Int)xChunkSize - 1, (Int)ximshape - 1
305 : );
306 0 : std::shared_ptr<Array<Bool>> subRegionMask;
307 0 : if (regionMask) {
308 0 : auto doMaskSlice = yDoMaskSlice;
309 0 : xRegMaskStart = 0;
310 0 : if (xblc < 0) {
311 0 : xRegMaskStart = -xblc;
312 0 : xRegMaskLength = regionMask->shape()[_dirAxes[0]] + xblc;
313 0 : doMaskSlice = True;
314 : }
315 0 : else if (xblc + xChunkSize > ximshape) {
316 0 : regMaskLength[_dirAxes[0]] = ximshape - xblc;
317 0 : doMaskSlice = True;
318 : }
319 : else {
320 0 : xRegMaskLength = xChunkSize;
321 : }
322 0 : if (doMaskSlice) {
323 0 : Slicer sl(regMaskStart, regMaskLength);
324 0 : subRegionMask.reset(new Array<Bool>(regMaskCopy(sl)));
325 0 : }
326 : else {
327 0 : subRegionMask.reset(new Array<Bool>(regMaskCopy));
328 : }
329 : }
330 0 : auto maskChunk = lattMask(planeBlc, planeTrc).copy();
331 0 : if (subRegionMask) {
332 0 : maskChunk = maskChunk && *subRegionMask;
333 : }
334 0 : if (anyTrue(maskChunk)) {
335 0 : auto chunk = plane(planeBlc, planeTrc);
336 0 : if (allTrue(maskChunk)) {
337 0 : alg->setData(chunk.begin(), chunk.size());
338 : }
339 : else {
340 0 : alg->setData(chunk.begin(), maskChunk.begin(), chunk.size());
341 : }
342 0 : res = alg->getStatistic(_statType);
343 0 : }
344 : else {
345 : // no pixels in region on which to do stats, mask grid point
346 0 : writeTo.pixelMask().putAt(False, outPos);
347 : }
348 0 : writeTo.putAt(res, outPos);
349 0 : }
350 : }
351 0 : nCount += nxpts;
352 0 : pm.update(nCount);
353 0 : }
354 0 : if (hasLoopAxes) {
355 0 : Bool done = False;
356 0 : uInt idx = 0;
357 0 : while (! done) {
358 0 : auto targetAxis = loopAxes[idx];
359 0 : ++blc[targetAxis];
360 0 : if (blc[targetAxis] == imshape[targetAxis]) {
361 0 : blc[targetAxis] = 0;
362 0 : ++idx;
363 0 : done = idx == loopAxes.size();
364 : }
365 : else {
366 0 : done = True;
367 : }
368 : }
369 : }
370 0 : }
371 0 : }
372 :
373 : std::shared_ptr<
374 : StatisticsAlgorithm<Double, Array<Float>::const_iterator,
375 : Array<Bool>::const_iterator>
376 : >
377 0 : StatImageCreator::_getStatsAlgorithm(String& algName) const {
378 0 : auto ac = _getAlgConf();
379 0 : switch(ac.algorithm) {
380 0 : case StatisticsData::CLASSICAL:
381 0 : algName = "classical";
382 0 : return std::shared_ptr<
383 : ClassicalStatistics<Double, Array<Float>::const_iterator,
384 : Array<Bool>::const_iterator>
385 : >(
386 : new ClassicalStatistics<
387 : Double, Array<Float>::const_iterator,
388 : Array<Bool>::const_iterator
389 0 : >()
390 0 : );
391 0 : case StatisticsData::CHAUVENETCRITERION:
392 0 : algName = "Chauvenet criterion/z-score";
393 0 : return std::shared_ptr<ChauvenetCriterionStatistics<
394 : Double, Array<Float>::const_iterator,
395 : Array<Bool>::const_iterator>
396 : >(
397 : new ChauvenetCriterionStatistics<
398 : Double, Array<Float>::const_iterator,
399 : Array<Bool>::const_iterator
400 : >(
401 : ac.zs, ac.mi
402 0 : )
403 0 : );
404 0 : default:
405 0 : ThrowCc("Unsupported statistics algorithm");
406 : }
407 : }
408 :
409 0 : void StatImageCreator::_nominalChunkInfo(
410 : std::shared_ptr<Array<Bool>>& templateMask, uInt& xBlcOff, uInt& yBlcOff,
411 : uInt& xChunkSize, uInt& yChunkSize, SPCIIF subimage
412 : ) const {
413 0 : Double xPixLength = 0;
414 0 : const auto& csys = subimage->coordinates();
415 0 : if (_xlen.getUnit().startsWith("pix")) {
416 0 : xPixLength = _xlen.getValue();
417 : }
418 : else {
419 0 : const auto& dc = csys.directionCoordinate();
420 0 : auto units = dc.worldAxisUnits();
421 0 : auto inc = dc.increment();
422 0 : Quantity worldPixSize(abs(inc[0]), units[0]);
423 0 : xPixLength = _xlen.getValue("rad")/worldPixSize.getValue("rad");
424 0 : }
425 0 : const auto shape = subimage->shape();
426 0 : ThrowIf(
427 : xPixLength >= shape[_dirAxes[0]] - 4,
428 : "x region length is nearly as large as or larger than the subimage extent of "
429 : + String::toString(shape[_dirAxes[0]]) + " pixels. Such a configuration is not supported"
430 : );
431 0 : ThrowIf(
432 : xPixLength <= 0.5,
433 : "x region length is only " + String::toString(xPixLength)
434 : + " pixels. Such a configuration is not supported"
435 : );
436 0 : Double yPixLength = xPixLength;
437 0 : if (_ylen.getValue() > 0) {
438 0 : if (_ylen.getUnit().startsWith("pix")) {
439 0 : yPixLength = _ylen.getValue();
440 : }
441 : else {
442 0 : const auto& dc = csys.directionCoordinate();
443 0 : auto units = dc.worldAxisUnits();
444 0 : auto inc = dc.increment();
445 0 : Quantity worldPixSize(abs(inc[1]), units[1]);
446 0 : yPixLength = _ylen.getValue("rad")/worldPixSize.getValue("rad");
447 0 : }
448 : }
449 0 : ThrowIf(
450 : yPixLength >= shape[_dirAxes[1]] - 4,
451 : "y region length is nearly as large as or larger than the subimage extent of "
452 : + String::toString(_dirAxes[1]) + " pixels. Such a configuration is not supported"
453 : );
454 0 : ThrowIf(
455 : yPixLength <= 0.5,
456 : "y region length is only " + String::toString(yPixLength)
457 : + " pixels. Such a configuration is not supported"
458 : );
459 0 : IPosition templateShape(shape.size(), 1);
460 0 : templateShape[_dirAxes[0]] = shape[_dirAxes[0]];
461 0 : templateShape[_dirAxes[1]] = shape[_dirAxes[1]];
462 0 : TempImage<Float> templateImage(templateShape, csys);
463 : // integer division, xq, yq must have integral values
464 0 : uInt xcenter = templateShape[_dirAxes[0]]/2;
465 0 : uInt ycenter = templateShape[_dirAxes[1]]/2;
466 0 : Quantity xq(xcenter, "pix");
467 0 : Quantity yq(ycenter, "pix");
468 0 : IPosition centerPix(templateShape.size(), 0);
469 0 : centerPix[_dirAxes[0]] = xcenter;
470 0 : centerPix[_dirAxes[1]] = ycenter;
471 0 : auto world = csys.toWorld(centerPix);
472 0 : static const Vector<Stokes::StokesTypes> dummyStokes;
473 0 : Record reg;
474 0 : if (_ylen.getValue() > 0) {
475 : AnnCenterBox box(
476 0 : xq, yq, _xlen, _ylen, csys, templateShape, dummyStokes
477 0 : );
478 0 : reg = box.getRegion()->toRecord("");
479 0 : }
480 : else {
481 : AnnCircle circle(
482 0 : xq, yq, _xlen, csys, templateShape, dummyStokes
483 0 : );
484 0 : reg = circle.getRegion()->toRecord("");
485 0 : }
486 0 : static const AxesSpecifier dummyAxesSpec;
487 : auto chunkImage = SubImageFactory<Float>::createSubImageRO(
488 : templateImage, reg, "", nullptr, dummyAxesSpec, False
489 0 : );
490 0 : auto blcOff = chunkImage->coordinates().toPixel(world);
491 0 : xBlcOff = rint(blcOff[_dirAxes[0]]);
492 0 : yBlcOff = rint(blcOff[_dirAxes[1]]);
493 0 : auto chunkShape = chunkImage->shape();
494 0 : xChunkSize = chunkShape[_dirAxes[0]];
495 0 : yChunkSize = chunkShape[_dirAxes[1]];
496 0 : templateMask.reset();
497 0 : if (chunkImage->isMasked()) {
498 0 : auto mask = chunkImage->getMask();
499 0 : if (! allTrue(mask)) {
500 0 : templateMask.reset(new Array<Bool>(mask));
501 : }
502 0 : }
503 0 : }
504 :
505 0 : void StatImageCreator::_doInterpolation(
506 : SPIIF output, TempImage<Float>& store,
507 : SPCIIF subImage, uInt nxpts, uInt nypts,
508 : Int xstart, Int ystart
509 : ) const {
510 0 : const auto imshape = subImage->shape();
511 0 : const auto xshape = imshape[_dirAxes[0]];
512 0 : const auto yshape = imshape[_dirAxes[1]];
513 0 : const auto ndim = subImage->ndim();
514 0 : *_getLog() << LogOrigin(getClass(), __func__)
515 : << LogIO::NORMAL << "Interpolate using "
516 0 : << _interpName << " algorithm." << LogIO::POST;
517 0 : Matrix<Float> result(xshape, yshape);
518 0 : Matrix<Bool> resultMask;
519 0 : if (_doMask) {
520 0 : resultMask.resize(IPosition(2, xshape, yshape));
521 0 : resultMask.set(True);
522 : }
523 0 : Matrix<Float> storage(nxpts, nypts);
524 0 : Matrix<Bool> storeMask(nxpts, nypts);
525 0 : std::pair<uInt, uInt> start;
526 0 : start.first = xstart;
527 0 : start.second = ystart;
528 0 : if (ndim == 2) {
529 0 : store.get(storage);
530 0 : store.getMask(storeMask);
531 0 : _interpolate(result, resultMask, storage, storeMask, start);
532 0 : output->put(result);
533 0 : if (_doMask) {
534 0 : output->pixelMask().put(resultMask && subImage->pixelMask().get());
535 : }
536 : }
537 : else {
538 : // get each direction plane in the storage lattice at each chan/stokes
539 0 : IPosition cursorShape(ndim, 1);
540 0 : cursorShape[_dirAxes[0]] = nxpts;
541 0 : cursorShape[_dirAxes[1]] = nypts;
542 0 : auto axisPath = _dirAxes;
543 0 : axisPath.append((IPosition::otherAxes(ndim, _dirAxes)));
544 0 : LatticeStepper stepper(store.shape(), cursorShape, axisPath);
545 0 : Slicer slicer(stepper.position(), stepper.endPosition(), Slicer::endIsLast);
546 0 : IPosition myshape(ndim, 1);
547 0 : myshape[_dirAxes[0]] = xshape;
548 0 : myshape[_dirAxes[1]] = yshape;
549 0 : for (stepper.reset(); ! stepper.atEnd(); stepper++) {
550 0 : auto pos = stepper.position();
551 0 : slicer.setStart(pos);
552 0 : slicer.setEnd(stepper.endPosition());
553 0 : storage = store.getSlice(slicer, True);
554 0 : storeMask = store.getMaskSlice(slicer, True);
555 0 : _interpolate(result, resultMask, storage, storeMask, start);
556 0 : output->putSlice(result, pos);
557 0 : if (_doMask) {
558 0 : output->pixelMask().putSlice(
559 0 : resultMask && subImage->pixelMask().getSlice(pos, myshape, True),
560 : pos
561 : );
562 : }
563 0 : }
564 0 : }
565 0 : }
566 :
567 0 : void StatImageCreator::setInterpAlgorithm(Interpolate2D::Method alg) {
568 0 : switch (alg) {
569 0 : case Interpolate2D::CUBIC:
570 0 : _interpName = "CUBIC";
571 0 : break;
572 0 : case Interpolate2D::LANCZOS:
573 0 : _interpName = "LANCZOS";
574 0 : break;
575 0 : case Interpolate2D::LINEAR:
576 0 : _interpName = "LINEAR";
577 0 : break;
578 0 : case Interpolate2D::NEAREST:
579 0 : _interpName = "NEAREST";
580 0 : break;
581 0 : default:
582 0 : ThrowCc("Unhandled interpolation method " + String::toString(alg));
583 : }
584 0 : _interpolater = Interpolate2D(alg);
585 0 : }
586 :
587 0 : void StatImageCreator::setStatType(StatisticsData::STATS s) {
588 0 : _statType = s;
589 0 : _statName = StatisticsData::toString(s);
590 0 : _statName.upcase();
591 0 : }
592 :
593 0 : void StatImageCreator::setStatType(const String& s) {
594 0 : String m = s;
595 : StatisticsData::STATS stat;
596 0 : m.downcase();
597 0 : if (m.startsWith("i")) {
598 0 : stat = StatisticsData::INNER_QUARTILE_RANGE;
599 : }
600 0 : else if (m.startsWith("max")) {
601 0 : stat = StatisticsData::MAX;
602 : }
603 0 : else if (m.startsWith("mea")) {
604 0 : stat = StatisticsData::MEAN;
605 : }
606 0 : else if (m.startsWith("meda") || m.startsWith("mad") || m.startsWith("x")) {
607 0 : stat = StatisticsData::MEDABSDEVMED;
608 : }
609 0 : else if (m.startsWith("medi")) {
610 0 : stat = StatisticsData::MEDIAN;
611 : }
612 0 : else if (m.startsWith("mi")) {
613 0 : stat = StatisticsData::MIN;
614 : }
615 0 : else if (m.startsWith("n")) {
616 0 : stat = StatisticsData::NPTS;
617 : }
618 0 : else if (m.startsWith("q1")) {
619 0 : stat = StatisticsData::FIRST_QUARTILE;
620 : }
621 0 : else if (m.startsWith("q3")) {
622 0 : stat = StatisticsData::THIRD_QUARTILE;
623 : }
624 0 : else if (m.startsWith("r")) {
625 0 : stat = StatisticsData::RMS;
626 : }
627 0 : else if (m.startsWith("si") || m.startsWith("st")) {
628 0 : stat = StatisticsData::STDDEV;
629 : }
630 0 : else if (m.startsWith("sums")) {
631 0 : stat = StatisticsData::SUMSQ;
632 : }
633 0 : else if (m.startsWith("sum")) {
634 0 : stat = StatisticsData::SUM;
635 : }
636 0 : else if (m.startsWith("v")) {
637 0 : stat = StatisticsData::VARIANCE;
638 : }
639 : else {
640 0 : ThrowCc("Statistic " + s + " not supported.");
641 : }
642 0 : _doProbit = m.startsWith("x");
643 0 : setStatType(stat);
644 0 : }
645 :
646 0 : void StatImageCreator::_interpolate(
647 : Matrix<Float>& result, Matrix<Bool>& resultMask,
648 : const Matrix<Float>& storage,
649 : const Matrix<Bool>& storeMask,
650 : const std::pair<uInt, uInt>& start
651 : ) const {
652 0 : auto shape = result.shape();
653 0 : auto imax = shape[0];
654 0 : auto jmax = shape[1];
655 : // x values change fastest in the iterator
656 0 : auto iter = result.begin();
657 0 : auto miter = resultMask.begin();
658 : // xcell and ycell are the current positions within the current
659 : // grid cell represented by an integer modulo pointsPerCell
660 0 : auto xCell = start.first == 0 ? 0 : _grid.first - start.first;
661 0 : auto yCell = start.second == 0 ? 0 : _grid.second - start.second;
662 0 : Int xStoreInt = start.first == 0 ? 0 : -1;
663 0 : Int yStoreInt = start.second == 0 ? 0 : -1;
664 0 : Double xStoreFrac = (Double)xCell/(Double)_grid.first;
665 0 : Double yStoreFrac = (Double)yCell/(Double)_grid.second;
666 0 : Vector<Double> storeLoc(2);
667 0 : storeLoc[0] = xStoreInt + xStoreFrac;
668 0 : storeLoc[1] = yStoreInt + yStoreFrac;
669 0 : for (uInt j=0; j<jmax; ++j, ++yCell) {
670 0 : Bool onYGrid = yCell == 0;
671 0 : if (yCell == _grid.second) {
672 0 : yCell = 0;
673 0 : ++yStoreInt;
674 0 : yStoreFrac = 0;
675 0 : onYGrid = True;
676 : }
677 : else {
678 0 : yStoreFrac = (Double)yCell/(Double)_grid.second;
679 : }
680 0 : storeLoc[1] = yStoreInt + yStoreFrac;
681 0 : xCell = start.first == 0 ? 0 : _grid.first - start.first;
682 0 : xStoreInt = start.first == 0 ? 0 : -1;
683 0 : xStoreFrac = (Double)xCell/(Double)_grid.first;
684 0 : for (uInt i=0; i<imax; ++i, ++xCell) {
685 0 : Bool onXGrid = xCell == 0;
686 0 : if (xCell == _grid.first) {
687 0 : xCell = 0;
688 0 : ++xStoreInt;
689 0 : onXGrid = True;
690 : }
691 0 : if (onXGrid && onYGrid) {
692 : // exactly on a grid point, no interpolation needed
693 : // just copy value directly from storage matrix
694 0 : *iter = storage(xStoreInt, yStoreInt);
695 0 : if (_doMask) {
696 0 : *miter = storeMask(xStoreInt, yStoreInt);
697 : }
698 : }
699 : else {
700 0 : xStoreFrac = (Double)xCell/(Double)_grid.first;
701 0 : storeLoc[0] = xStoreInt + xStoreFrac;
702 0 : if (! _interpolater.interp(*iter, storeLoc, storage, storeMask)) {
703 0 : ThrowIf(
704 : ! _doMask,
705 : "Logic error: bad interpolation but there is no mask to set"
706 : );
707 0 : *miter = False;
708 : }
709 : }
710 0 : ++iter;
711 0 : if (_doMask) {
712 0 : ++miter;
713 : }
714 : }
715 : }
716 0 : }
717 :
718 : }
|
