#################################################################### Using lsst.afw.geom.SpanSet to represent irregular regions of pixels #################################################################### A ``SpanSet`` is a class that represents regions of pixels on an image. As the name implies, a ``SpanSet`` is a collection of ``Span`` objects which behave as a mathematical set. The class provides methods for doing set operations, as well as convenience method for operating on data products at the locations defined by the ``SpanSet``. Creating SpanSets ================= The following example shows a few different ways to create a ``SpanSet`` in Python (C++ is similar but with differences in syntax). .. code-block:: python from lsst.afw.geom import Span, SpanSet, Stencil # Construct a SpanSet from a list of Spans # Spans are constructed as y, begin x, end x coordinates spanList = [Span(1, 0, 2), Span(2, 0, 2), Span(3, 0, 2)] spanSetFromList = SpanSet(spanList) # Construct the same SpanSet from a shape radius = 1 spanSetFromShape = SpanSet.fromShape(radius, Stencil.BOX, offset=(1, 2)) # Verify that the two are the same assert(spanSetFromList == spanSetFromShape) fromShape --------- The last example introduced a useful factory method for creating ``SpanSets`` that are worth exploring further. ``fromShape`` can be called in one of two ways, specifying a radius, stencil (defaults to a circle), and an offset (defaults to 0, 0), or passing in an ``afw.geom.ellipse.Ellipse`` object. There are three options when creating a ``SpanSet`` from a ``afw.geom.Stencil``; ``BOX``, ``CIRCLE``, and ``MANHATTAN`` (a diamond). Each shape is centered on the 0, 0 pixel and extends the given radius in all directions according to the shape selected. If an origin other than 0, 0 is desired an x, y offset to the origin can be provided. SpanSets are immutable ====================== When working with ``SpanSet`` objects it is important to note that they are immutable once they are created. As such, all method calls on a ``SpanSet`` object which operate on a ``SpanSet`` create, modify and return a new ``SpanSet`` object. This design is intended to offer a safety guarantee to programmers. If two different bits of code each have a shared pointer to the same ``SpanSet`` one can not modify the region and possible corrupt what the other bit of code is trying to do. SpanSets do not need to be Spatially contiguous =============================================== Another important note is that regions defined by a ``SpanSet`` are not guaranteed to be contiguous. There is however a method which will return if a given ``SpanSet`` is contiguous, and a method which will split a non-contiguous ``SpanSet`` into a standard vector (or list in Python) of contiguous ``SpanSet`` objects. The following example demonstrates this behavior alongside set operations. .. code-block:: python from lsst.afw.geom import SpanSet, Stencil # Construct some SpanSets radius = 5 ss1 = SpanSet.fromShape(radius) ss2 = SpanSet.fromShape(radius, offset=(6, 6)) # Create a union of the immutable sets ss3 = ss1.union(ss2) # Erode the combine SpanSet with a circular kernel newRadius = 4 ssEroded = ss3.eroded(newRadius, Stencil.CIRCLE) # Result is now non-contiguous assert(ssEroded.isContiguous() == False) # Split the non contiguous SpanSet spanSetSplitList = ssEroded.split() assert(len(spanSetSplitList) == 2) SpanSets, Masks, Images, and Ndarray ==================================== Masks ----- Masks and ``SpanSet``\s can often be used in similar fashion to represent areas of interest in an image. As such ``SpanSet``\s provide many utility functions to make transporting this information between these two data types similar. As Masks are similar to ``SpanSet``\s they can even participate in ``SpanSet`` mathematical operations. Examples for some of the ways masks and ``SpanSet``\s work together can be found in the example below. .. code-block:: python from lsst.afw.geom import SpanSet from lsst.afw.image import MaskU # Create a mask to be populated size = 10 mask = MaskU(size, size) # Create a SpanSet which represents the pixels to be set in the mask, and # set bit two radius = 4 ss = SpanSet.fromShape(radius, offset=(4, 4)) bitMask = 2 ss.setMask(mask, bitMask) # Intersect not (~) the SpanSet with the mask, the result should be a null # SpanSet ssIntersectNot = ss.intersectNot(mask) # Convert the mask into a SpanSet and verify it evaluates equal to the # original newSS = SpanSet.fromMask(mask) assert(newSS == ss) Images ------ As mentioned above the ``SpanSet`` class is used to encode sets of x, y locations on an image. These locations can be used to interact with additional images through a series of convenience methods demonstrated in the subsequent example. .. code-block:: python from lsst.afw.geom import SpanSet from lsst.afw.image import ImageI # Define two different spans sets of differing sized centered at different # positions radius1 = 3 radius2 = 2 spanSet1 = SpanSet.fromShape(radius1, offset=(3, 3)) spanSet2 = SpanSet.fromShape(radius2, offset=(7, 7)) # Create two different Images, of the same size imageSize = 10 image1 = ImageI(imageSize, imageSize) image2 = ImageI(imageSize, imageSize) # Use the SpanSets to set pixels in each image to an arbitrary value spanSet1.setImage(image1, 10) spanSet2.setImage(image2, 15) # Use the second SpanSet to copy the values from image2 into image1 at the # positions defined in spanSet2 spanSet2.copyImage(image2, image1) # Show the results print(image1.getArray()) # Output: #[[ 0 0 0 10 0 0 0 0 0 0] # [ 0 10 10 10 10 10 0 0 0 0] # [ 0 10 10 10 10 10 0 0 0 0] # [10 10 10 10 10 10 10 0 0 0] # [ 0 10 10 10 10 10 0 0 0 0] # [ 0 10 10 10 10 10 0 15 0 0] # [ 0 0 0 10 0 0 15 15 15 0] # [ 0 0 0 0 0 15 15 15 15 15] # [ 0 0 0 0 0 0 15 15 15 0] # [ 0 0 0 0 0 0 0 15 0 0]] Ndarrays -------- A ``SpanSet`` can also be used to extract or insert values from / into ndarrays while expanding or reducing dimensionality. The ``flatten`` method extracts data from an array at locations defined by the ``SpanSet`` and returns (or inserts into) an ``ndarray`` with one less dimension. The ``unflatten`` method does the opposite. The ``flatten`` method takes the first two dimensions of the ndarray as the dimensions to flattened, indexed at the locations of the ``SpanSet``. If a ``SpanSet`` is defined to cover a 5x5 area, and is used to flatten a 5x5x4x10 array, the resulting array will be 25x4x10. Below is small example. .. code-block:: python import numpy as np from lsst.afw.geom import SpanSet, Stencil # Create a 2D array with ascending values and a SpanSet of a sub region dims = 5 array = np.arange(dims * dims).reshape(dims, dims) radius = 1 ss = SpanSet.fromShape(radius, Stencil.BOX, offset=(1, 1)) # Extract the sub region into a flattened array flat = ss.flatten(array) # Show the flattned values print(flat.shape) # Output: # (9,) print(flat) # Output: # [ 0 1 2 5 6 7 10 11 12] Using indices (Python only) =========================== A ``SpanSet`` is a representation of coordinates that is very efficient in terms of memory usage. This however does not always lend itself to the Python / Numpy style of programming, owing to the need to do a double loop to access the actual coordinates. In order to support a more natural way of programming with Python / Numpy the ``SpanSet`` class provides an ``indices`` method. This method, when called, returns a tuple of two lists. The first list contains the y coordinate for each point in the ``SpanSet``, and the second provides the corresponding x coordinates. Note this is different that the x, y ordering common through other parts of the LSST code base, but was chosen to be similar to numpy.indices, and the ordering of Numpy arrays. This representation is less memory efficient and should be used thoughtfully, but enables coding styles similar to the following example. .. code-block:: python import numpy as np from lsst.afw.geom import SpanSet, Stencil # Create a numpy array to work with arrayDim = 5 dataArray = np.zeros((arrayDim, arrayDim)) # Create a SpanSet which indexes all the x, y locations in the data array radius = 2 ss = SpanSet.fromShape(2, Stencil.BOX, offset=(2, 2)) # Get the indices corresponding to the SpanSet and use it to set values in # the data array yind, xind = ss.indices() dataArray[yind, xind] = 9 # Show the modified data array print(dataArray) # Output: # [[ 9. 9. 9. 9. 9.] # [ 9. 9. 9. 9. 9.] # [ 9. 9. 9. 9. 9.] # [ 9. 9. 9. 9. 9.] # [ 9. 9. 9. 9. 9.]] Using applyFunctor (C++ only) ============================= When a ``SpanSet`` class is used in C++ there is a useful convenience function that is unavailable from the Python interface called ``applyFunctor``. This method is meant to simplify the complexities of doing operations at the locations defined by a ``SpanSet`` on data of mixed types and shapes. The key point to this functionality is specifying some function-like object that contains the operation to be done as if each pixel only represented a single value. The ``applyFunctor`` method then takes this operation, and the data to be operated on, and iterates over the data types in such a way that the operator is supplied only one set of values at a time. The method can handle data types of ``Image``, ``MaskedImage``, ``ndarrays``, numeric values (i.e. a float), and iterators. Any number of data types may be supplied (constrained by the number of arguments the supplied function operator takes). The following contains a snippet of C++ code as a demonstration, a full working example can be found in the C++ unit test, and more detail on syntax can be found in the applyFunctor doxygen. .. code-block:: cpp afwImage::Image image(5, 5, 1); afwImage::Image outputImage(5, 5, 0); std::vector vec(5*5, 2); ndarray::Array ndAr = ndarray::allocate(ndarray::makeVector(5,5)); ndAr.deep() = 1; int constant = 2; auto ss = afwGeom::SpanSet::fromShape(2, afwGeom::Stencil::BOX, offset=afwGeom::Box2I(2,2)) // The Point2I argument says where in the SpanSet the operator is being applied // but is unused in this example ss.applyFunctor([]( afGeom::Point2I const &, int & out, int const & inIm, int const & inVec, int const & ndAr, int number){ out = inIm * inVec * ndAr / number; }, outputImage, image, vec, ndarray::ndImage(ndAr), constant);