Hi,
If I read you correclty,
$rows=sumover sumover ($x*$mask->dummy(2,7))
works. Is it faster?
Ingo
On 08/10/2012 05:31 PM, Tim Haines wrote:
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
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
