For applications like that I use
$a->map(t_radial(o=>$ctr),[@newdims],{m=>'g',b=>'pt'})->average
or something like that. More to follow when I am not out hiking.
(Mobile)
On Aug 10, 2012, at 9:31 AM, Tim Haines <[email protected]> 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
