LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Image Locators

(Return to Images)

Iterators provide access to an image, pixel by pixel. You often want access to neighbouring pixels (e.g. computing a gradient, or smoothing). Let's consider the problem of smoothing with a

1 2 1
2 4 2
1 2 1

kernel (the code's in image2.cc):

Start by including Image.h defining a namespace for clarity:

#include "lsst/geom.h"
namespace image = lsst::afw::image;
using ImageT = image::Image<int>;
int main() {
afw::table::Key< afw::table::Array< ImagePixelT > > image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.

Declare an Image

Set the image to a ramp

for (int y = 0; y != in.getHeight(); ++y) {
for (ImageT::xy_locator ptr = in.xy_at(0, y), end = in.xy_at(in.getWidth(), y); ptr != end;
++ptr.x()) {
*ptr = y;
}
int end
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:48
}

That didn't gain us much, did it? The code's a little messier than using x_iterator. But now we can add code to calculate the smoothed image. First make an output image, and copy the input pixels:

//
// Convolve with a pseudo-Gaussian kernel ((1, 2, 1), (2, 4, 2), (1, 2, 1))
//
ImageT out(in.getDimensions()); // Make an output image the same size as the input image
out.assign(in);

(we didn't need to copy all of them, just the ones around the edge that we won't smooth, but this is an easy way to do it).

Now do the smoothing:

for (int y = 1; y != in.getHeight() - 1; ++y) {
for (ImageT::xy_locator ptr = in.xy_at(1, y), end = in.xy_at(in.getWidth() - 1, y),
optr = out.xy_at(1, y);
ptr != end; ++ptr.x(), ++optr.x()) {
*optr = ptr(-1, -1) + 2 * ptr(0, -1) + ptr(1, -1) + 2 * ptr(-1, 0) + 4 * ptr(0, 0) +
2 * ptr(1, 0) + ptr(-1, 1) + 2 * ptr(0, 1) + ptr(1, 1);
}
}

(N.b. you don't really want to do this; not only is this kernel separable into 1 2 1 in first the x then the y directions, but lsst::afw::math can do convolutions for you).

Here's a faster way to do the same thing (the use of an Image::Ptr is just for variety)

//
// Do the same thing a faster way, using cached_location_t
//
std::shared_ptr<ImageT> out2(new ImageT(in.getDimensions()));
out2->assign(in);
using xy_loc = ImageT::const_xy_locator;
for (int y = 1; y != in.getHeight() - 1; ++y) {
// "dot" means "cursor location" in emacs
xy_loc dot = in.xy_at(1, y), end = in.xy_at(in.getWidth() - 1, y);
xy_loc::cached_location_t nw = dot.cache_location(-1, -1);
xy_loc::cached_location_t n = dot.cache_location(0, -1);
xy_loc::cached_location_t ne = dot.cache_location(1, -1);
xy_loc::cached_location_t w = dot.cache_location(-1, 0);
xy_loc::cached_location_t c = dot.cache_location(0, 0);
xy_loc::cached_location_t e = dot.cache_location(1, 0);
xy_loc::cached_location_t sw = dot.cache_location(-1, 1);
xy_loc::cached_location_t s = dot.cache_location(0, 1);
xy_loc::cached_location_t se = dot.cache_location(1, 1);
for (ImageT::x_iterator optr = out2->row_begin(y) + 1; dot != end; ++dot.x(), ++optr) {
*optr = dot[nw] + 2 * dot[n] + dot[ne] + 2 * dot[w] + 4 * dot[c] + 2 * dot[e] + dot[sw] +
2 * dot[s] + dot[se];
}
def dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, *args, **kwargs)
Definition: ds9.py:100
double w
Definition: CoaddPsf.cc:69
}

The xy_loc::cached_location_t variables remember relative positions.

We can rewrite this to move setting nw, se etc. out of the loop:

//
// Do the same calculation, but set nw etc. outside the loop
//
xy_loc pix11 = in.xy_at(1, 1);
xy_loc::cached_location_t nw = pix11.cache_location(-1, -1);
xy_loc::cached_location_t n = pix11.cache_location(0, -1);
xy_loc::cached_location_t ne = pix11.cache_location(1, -1);
xy_loc::cached_location_t w = pix11.cache_location(-1, 0);
xy_loc::cached_location_t c = pix11.cache_location(0, 0);
xy_loc::cached_location_t e = pix11.cache_location(1, 0);
xy_loc::cached_location_t sw = pix11.cache_location(-1, 1);
xy_loc::cached_location_t s = pix11.cache_location(0, 1);
xy_loc::cached_location_t se = pix11.cache_location(1, 1);
for (int y = 1; y != in.getHeight() - 1; ++y) {
// "dot" means "cursor location" in emacs
xy_loc dot = in.xy_at(1, y), end = in.xy_at(in.getWidth() - 1, y);
for (ImageT::x_iterator optr = out2->row_begin(y) + 1; dot != end; ++dot.x(), ++optr) {
*optr = dot[nw] + 2 * dot[n] + dot[ne] + 2 * dot[w] + 4 * dot[c] + 2 * dot[e] + dot[sw] +
2 * dot[s] + dot[se];
}
}

You may have noticed that that kernel isn't normalised. We could change the coefficients, but that'd slow things down for integer images (such as the one here); but we can normalise after the fact by making an Image that shares pixels with the central part of out2 and manipulating it via overloaded operator/=

//
// Normalise the kernel. I.e. divide the smoothed parts of image2 by 16
//
{
ImageT center = ImageT(
*out2,
center /= 16;
}
An integer coordinate rectangle.
Definition: Box.h:55

N.b. you can use the iterator embedded in the locator directly if you really want to, e.g.

for (int y = 0; y != in.getHeight(); ++y) {
for (ImageT::xy_x_iterator ptr = in.xy_at(0, y).x(), end = in.xy_at(in.getWidth(), y).x(); ptr != end;
++ptr) {
*ptr = 0;
}
}

we called the iterator xy_x_iterator, not x_iterator, for consistency with MaskedImage.

Finally write some output files and close out main():

//
// Save those images to disk
//
out.writeFits("foo.fits");
out2->writeFits("foo2.fits");
return 0;
}