Sorry for the brief answer earlier. There are a couple of ways to do it. For
selection (as you're contemplating) you can use whichND and indexND.
(Unfortunately, the most obvious technique, which I won't even describe here,
doesn't work because of the way non-affine dataflow works under the hood). For
selection, you can avoid the memory footprint with:
$rv = rvals($x->dim(0),$x->dim(1));
$out = zeroes($rmax, $x->dim(2));
for $r(0..$rmax) {
$dex = whichND( ($rv->floor < $r) & ($rv->ceil >= $r) );
$out->indexND($dex) .= $x->indexND($dex)->sumover;
}
That works because indexND will thread over the spectral dimension - so you
select just the (x,y) pixels you want and accumulate them, while preserving the
threaded spectral dimension. At each iteration, $dex gets a (2xP) set of
2-coordinates to accumulate, and $out->indexND($dex) has dim PxK.
The earlier example I used - map - may even do what you want even better:
use PDL::Transform
$out = $x->map(t_radial(o=>$ctr), # Radially expand around the coordinate
($ctr)
[360,$x->dim(1)/2], # Map into a 360xr output
{method=>'gaussian', # Use spatially-variable antialias
filtering
bound=>'trunc', # Truncate original data set
phot=>'flux'}) # Preserve flux rather than value
->sumover; # Sum over theta
If you have a FITS header in your input with a scientific coordinate system
defined (via CRPIXn/CRVALn/CDELTn/CROTA), then $ctr (which defaults to (0,0) if
you leave it out) is taken to be in scientific coordinates. If you have no
FITS header, then it'll use pixel coordinates.
The advantages of using map are (a) that it will actually interpolate smoothly
and, arguably, correctly around the original irregular pixel grid and (b) that
separating coordinate transformations from resampling turns out to be really
useful for a lot of things. The disadvantage is that it has more moving parts
than a simple index and sum.
On Aug 10, 2012, at 10:10 AM, Craig DeForest wrote:
> 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
