Shapelet Image Decomposition
Very brief notes about the "Shapelet Decomposition" query
VisIt provides a query that decomposes an image (or rectilinear grid) by a 2d combination of cartesian shapelet basis functions. The full decomposition result & coefficients are accessible via the GetQueryOutputObject() method in VisIt's cli. GetQueryOutputObject()returns a python representation of a "MapNode" (basically a nested python dictionary tree) that contains the following values:
- nmax - Basis combination limiter (set by used when query is called)
- beta - Basis Max (set by used when query is called)
- width - Input image width
- height - Input image height
- coeffs - Vector of packed basis coefficients
- error - Reconstruction error
Decomposition coefficients are packed in row-major order (the x component basis index is fast varying). With this packing some higher frequency (and lower power) coefficients are inserted before some lower frequency (and higher power) coefficients causing a saw tooth pattern. For example - the coefficient for basis (i,j) = (63,0), shows up before (i,j) = (0,1).
Here is a simple python method to compute the proper flattened index given a basis combination (i,j):
def shapelet_basis_coeff_index(i,j,nmax): bsize = (nmax * nmax +nmax) / 2 idx = nmax -j idx = bsize - (idx * idx + idx) / 2 + i return idx
This ordering is simple to produce with a nested for loop, but in retrospect was not a great choice because nmax influences the packing.
Here is a simple python mapping for a better ordering not influenced by nmax:
def shapelet_basis_index(i,j): ij = i+j idx = i + (ij * ij + ij) / 2 return idx
With this scheme the basis coefficients from a decomposition using a lower nmax would be directly comparable to a subset of the coefficients from a decomposition using a higher nmax. (This is the one the IDL Shapelets package uses) Unpacking & repacking the basis coefficients into a better order is pretty cheap - for example with nmax=64 you have 2080 basis coefficients. The numerics of evaluation of the high order gauss-hermite polynomials already limit the foreseeable values nmax - so for the time being we don't have any plans to switch the output to the better packing scheme.