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