Hi, Craig. This looks very interesting. I have never played with PDL::Transform before. Indeed, one of the concerns I had using my rvals sampling technique was the non-roundness of my annuli. My immediate concern with using the radial transform is that it does not provide information about the center pixel (spectrum), and I need to retain the it for my data analysis as it has the highest S/N.
# My indexing method [ [ 17 47 77 107 137 167 197] [ 136 376 616 856 1096 1336 1576] [ 272 752 1232 1712 2192 2672 3152] [ 10 160 310 460 610 760 910] ] # t_radial [ [ 154.48246 429.07182 703.66119 978.25055 1252.8399 1527.4293 1802.0186] [ 195.70768 586.3652 977.02273 1367.6802 1758.3378 2148.9953 2539.6528] [ 33.919814 176.84179 319.76376 462.68574 605.60771 748.52968 891.45166] ] I would just cat the central spectrum to the output of map, but I am concerned about flux scaling. From P::Transform::map for phot=>flux, Total flux is preserved over the transformation, so that the brightness integral over image regions is preserved. Parts of the image that are shrunk wind up brighter; parts that are enlarged end up fainter. Since the central pixel is most likely neither shrunk nor enlarged, is it safe to assume that it is unaffected by flux scaling? Many thanks for your help. - Tim On Fri, Aug 10, 2012 at 12:05 PM, Craig DeForest <[email protected]>wrote: > 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
