Hello, all.
I have a datacube of spatially-resolved spectra where each pixel in one of
the NxM planes is one of the K flux values along the wavelength axis. I am
trying to sum the spectra that lie on concentric annuli in the planes as
determined by rvals. Here is what I have thus far.
use PDL;
my $x = sequence(double,5,6,7);
my $annuli = $x->slice(':,:,(0)')->sever->long->inplace->rvals;
my $mask = $annuli == $annuli->uniq->dummy(0)->dummy(0);
# Letting C = $annuli->nelem, we thread the $mask cube planes
# over $x by inserting a single dimension at the third dim index
# of $mask.
#
# NxMxK * NxMx1xC -> (clump/sumover) -> KxC image
my $rows = ($x * $mask->dummy(2,1))->clump(2)->sumover;
Where $annuli is
[
[3 3 3 3 3]
[2 2 2 2 2]
[2 1 1 1 2]
[2 1 0 1 2]
[2 1 1 1 2]
[2 2 2 2 2]
]
and (for example) the mask for the third annulus ($mask->slice(':,:,(2)'))
is
[
[0 0 0 0 0]
[1 1 1 1 1]
[1 0 0 0 1]
[1 0 0 0 1]
[1 0 0 0 1]
[1 1 1 1 1]
]
as expected.
The last line produces the desired result, but at great cost in memory
footprint. Threading each plane in the $mask cube to select only the $x
cube spectra lying along the current annulus creates a 4D hypercube. Given
that some of my cubes have ~20 annuli and are ~20 MB in size, this gives me
a memory footprint of ~400 MB for the intermediate hypercube. This is
roughly a factor of 1000 greater than the size of the result in $rows.
Clearly, I need to find a way of looking up which spectra I want, rather
than generating the whole result and then clumping to sum. I tried using
whichND, but, while it gives me the correct set of indices for each
annulus...
whichND($radii == $radii->uniq->dummy(0)->dummy(0))->xchg(0,1) =
[
[2 1 2 3 1 3 1 2 3 0 1 2 3 4 0 4 0 4 0 4 0 1 2 3 4 0 1 2 3 4]
[3 2 2 2 3 3 4 4 4 1 1 1 1 1 2 2 3 3 4 4 5 5 5 5 5 0 0 0 0 0]
[0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3] <-- want
to thread over these
]
, I couldn't figure out how to make it group the indices by annulus number
so that I can thread over them. I considered range, but it seems to be more
for selecting contiguous regions. indexND looked plausible, but I couldn't
figure out how to get it to thread over the Z dimension of my cube. If I
could find a way of making whichND group the annuli, I would be able to use
dice, but no such luck.
Any thoughts are greatly appreciated.
Thanks.
- Tim
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
