Hmmm... I wouldn't put much stock in the flux preservation around singularities
-- it works by dividing by the determinant of the transformation's Jacobian,
which is fine in smooth areas but really awful for that last pixel. It's
probably better to use direct calculation:
$out = $x->map(...)->average * (xvals($out_r_dim)+0.5)
where the "..." is all the parameters we discussed earlier but without the
"phot=>flux" option). You might also want to look into either using "pix=>1"
to prevent autoscaling (in which case you need to compose t_radial with an
appropriate linear transformation to handle the output pixel grid:
$t = t_linear( scale=>[$output_dim_0 / 2 / PI, 1] ) x t_radial(...);
and use $t instead of the t_radial() expression in your map() call), or
"orange=>[...]" to set the output coordinate range for the autoscaler.
Cheers,
Craig
On Aug 10, 2012, at 11:50 AM, Tim Haines wrote:
> 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
