option.
=item m, matrix, Matrix
The transformation matrix. It does not even have to be square, if you want
to change the dimensionality of your input. If it is invertible (note:
must be square for that), then you automagically get an inverse transform too.
=item pre, preoffset, offset, Offset
The vector to be added to the data before they get multiplied by the matrix
(equivalent of CRVAL in FITS, if you are converting from scientific to
pixel units).
=item post, postoffset, shift, Shift
The vector to be added to the data after it gets multiplied by the matrix
(equivalent of CRPIX-1 in FITS, if youre converting from scientific to pixel
units).
=item d, dim, dims, Dims
Most of the time it is obvious how many dimensions you want to deal
with: if you supply a matrix, it defines the transformation; if you
input offset vectors in the C and C options, those define
the number of dimensions. But if you only supply scalars, there is no way
to tell and the default number of dimensions is 2. This provides a way
to do, e.g., 3-D scaling: just set C<{s=>, dims=>3}> and
you are on your way.
=back
NOTES
the type/unit fields are currently ignored by t_linear.
=cut
@PDL::Transform::Linear::ISA = ('PDL::Transform');
sub t_linear { new PDL::Transform::Linear(@_); }
sub PDL::Transform::Linear::new {
my($class) = shift;
my($o) = $_[0];
if(!(ref $o)) {
$o = {@_};
}
my($me) = PDL::Transform::new($class);
$me->{name} = "linear";
$me->{params}->{pre} = _opt($o,['pre','Pre','preoffset','offset',
'Offset','PreOffset','Preoffset'],0);
$me->{params}->{pre} = pdl($me->{params}->{pre})
if(defined $me->{params}->{pre});
$me->{params}->{post} = _opt($o,['post','Post','postoffset','PostOffset',
'shift','Shift'],0);
$me->{params}->{post} = pdl($me->{params}->{post})
if(defined $me->{params}->{post});
$me->{params}->{matrix} = _opt($o,['m','matrix','Matrix','mat','Mat']);
$me->{params}->{matrix} = pdl($me->{params}->{matrix})
if(defined $me->{params}->{matrix});
$me->{params}->{rot} = _opt($o,['r','rot','rota','rotation','Rotation']);
$me->{params}->{rot} = 0 unless defined($me->{params}->{rot});
$me->{params}->{rot} = pdl($me->{params}->{rot});
my $o_dims = _opt($o,['d','dim','dims','Dims']);
$o_dims = pdl($o_dims)
if defined($o_dims);
my $scale = _opt($o,['s','scale','Scale']);
$scale = pdl($scale)
if defined($scale);
# Figure out the number of dimensions to transform, and,
# if necessary, generate a new matrix.
if(defined($me->{params}->{matrix})) {
$me->{idim} = $me->{params}->{matrix}->dim(0);
$me->{odim} = $me->{params}->{matrix}->dim(1);
} else {
if(defined($scale) &&
UNIVERSAL::isa($scale,'PDL') &&
$scale->getndims > 0) {
$me->{idim} = $me->{odim} = $scale->dim(0);
$me->{odim} = $scale->dim(0);
} elsif(defined($me->{params}->{pre}) &&
UNIVERSAL::isa($me->{params}->{pre},'PDL') &&
$me->{params}->{pre}->getndims > 0) {
$me->{idim} = $me->{odim} = $me->{params}->{pre}->dim(0);
} elsif(defined($me->{params}->{post}) &&
UNIVERSAL::isa($me->{params}->{post},'PDL') &&
$me->{params}->{post}->getndims > 0) {
$me->{idim} = $me->{odim} = $me->{params}->{post}->dim(0);
} elsif(defined($o_dims)) {
$me->{idim} = $me->{odim} = $o_dims;
} else {
print "PDL::Transform::Linear: assuming 2-D transform (set dims option to change)\n" if($PDL::Transform::debug);
$me->{idim} = $me->{odim} = 2;
}
$me->{params}->{matrix} = PDL->zeroes($me->{idim},$me->{odim});
$me->{params}->{matrix}->diagonal(0,1) .= 1;
}
### Handle rotation option
my $rot = $me->{params}->{rot};
if(defined($rot)) {
# Subrotation closure -- rotates from axis $d->(0) --> $d->(1).
my $subrot = sub {
my($d,$angle,$m)=@_;
my($i) = identity($m->dim(0));
my($subm) = $i->dice($d,$d);
$angle = $angle->at(0)
if(UNIVERSAL::isa($angle,'PDL'));
my($a) = $angle * $DEG2RAD;
$subm .= $subm x pdl([cos($a),-sin($a)],[sin($a),cos($a)]);
$m .= $m x $i;
};
if(UNIVERSAL::isa($rot,'PDL') && $rot->nelem > 1) {
if($rot->ndims == 2) {
$me->{params}->{matrix} x= $rot;
} elsif($rot->nelem == 3) {
my $rm = identity(3);
# Do these in reverse order to make it more like
# function composition!
&$subrot(pdl(0,1),$rot->at(2),$rm);
&$subrot(pdl(2,0),$rot->at(1),$rm);
&$subrot(pdl(1,2),$rot->at(0),$rm);
$me->{params}->{matrix} .= $me->{params}->{matrix} x $rm;
} else {
barf("PDL::Transform::Linear: Got a strange rot option -- giving up.\n");
}
} else {
&$subrot(pdl(0,1),$rot,$me->{params}->{matrix});
}
}
#
# Apply scaling
#
$me->{params}->{matrix}->diagonal(0,1) *= $scale
if defined($scale);
#
# Check for an inverse and apply it if possible
#
my($o2);
if($me->{params}->{matrix}->det($o2 = {lu=>undef})) {
$me->{params}->{inverse} = $me->{params}->{matrix}->inv($o2);
} else {
delete $me->{params}->{inverse};
}
$me->{params}->{idim} = $me->{idim};
$me->{params}->{odim} = $me->{odim};
##############################
# The meat -- just shift, matrix-multiply, and shift again.
$me->{func} = sub {
my($in,$opt) = @_;
my($d) = $opt->{matrix}->dim(0)-1;
barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n")
if($in->dim(0) < $d);
my($a) = $in->(0:$d)->copy + $opt->{pre};
my($out) = $in->is_inplace ? $in : $in->copy;
$out->(0:$d) .= $a x $opt->{matrix} + $opt->{post};
return $out;
};
$me->{inv} = (defined $me->{params}->{inverse}) ? sub {
my($in,$opt) = @_;
my($d) = $opt->{inverse}->dim(0)-1;
barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n")
if($in->dim(0) < $d);
my($a) = $in->(0:$d)->copy - $opt->{post};
my($out) = $in->is_inplace ? $in : $in->copy;
$out->(0:$d) .= $a x $opt->{inverse} - $opt->{pre};
$out;
} : undef;
return $me;
}
sub PDL::Transform::Linear::stringify {
package PDL::Transform::Linear;
my($me) = shift; my($out) = SUPER::stringify $me;
my $mp = $me->{params};
if(!($me->{is_inverse})){
$out .= "Pre-add: ".($mp->{pre})."\n"
if(defined $mp->{pre});
$out .= "Post-add: ".($mp->{post})."\n"
if(defined $mp->{post});
$out .= "Forward matrix:".($mp->{matrix})
if(defined $mp->{matrix});
$out .= "Inverse matrix:".($mp->{inverse})
if(defined $mp->{inverse});
} else {
$out .= "Pre-add: ".(-$mp->{post})."\n"
if(defined $mp->{post});
$out .= "Post-add: ".(-$mp->{pre})."\n"
if(defined $mp->{pre});
$out .= "Forward matrix:".($mp->{inverse})
if(defined $mp->{inverse});
$out .= "Inverse matrix:".($mp->{matrix})
if(defined $mp->{matrix});
}
$out =~ s/\n/\n /go;
$out;
}
======EOD_t_linear======
pp_add_exported('t_scale');
pp_addpm(<<'======EOD_t_scale======');
=head2 t_scale
=for usage
$f = t_scale()
=for ref
Convenience interface to L.
t_scale produces a tranform that scales around the origin by a fixed
amount. It acts exactly the same as C\)>.
=cut
sub t_scale {
my($scale) = shift;
my($b) = shift;
return t_linear(scale=>$scale,%{$b})
if(ref $b eq 'HASH');
t_linear(Scale=>$scale,$b,@_);
}
======EOD_t_scale======
pp_add_exported('t_offset ');
pp_addpm(<<'======EOD_t_offset ======');
=head2 t_offset
=for usage
$f = t_offset()
=for ref
Convenience interface to L.
t_offset produces a transform that shifts the origin to a new location.
It acts exactly the same as C\)>.
=cut
sub t_offset {
my($pre) = shift;
my($b) = shift;
return t_linear(pre=>$pre,%{$b})
if(ref $b eq 'HASH');
t_linear(pre=>$pre,$b,@_);
}
======EOD_t_offset ======
pp_add_exported('t_rot');
pp_addpm(<<'======EOD_t_rot======');
=head2 t_rot
=for usage
$f = t_rot()
=for ref
Convenience interface to L.
t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or
N-D (matrix). It acts exactly the same as C\)>.
=cut
*t_rot = \&t_rotate;
sub t_rotate {
my $rot = shift;
my($b) = shift;
return t_linear(rot=>$rot,%{$b})
if(ref $b eq 'HASH');
t_linear(rot=>$rot,$b,@_);
}
======EOD_t_rot======
pp_add_exported('t_fits');
pp_addpm(<<'======EOD_t_fits======');
=head2 t_fits
=for usage
$f = t_fits($fits,[option]);
=for ref
FITS pixel-to-scientific transformation with inverse
You feed in a hash ref or a PDL with one of those as a header, and you
get back a transform that converts 0-originated, pixel-centered
coordinates into scientific coordinates via the transformation in the
FITS header. For most FITS headers, the transform is reversible, so
applying the inverse goes the other way. This is just a convenience
subclass of PDL::Transform::Linear, but with unit/type support
using the FITS header you supply.
For now, this transform is rather limited -- it really ought to
accept units differences and stuff like that, but they are just
ignored for now. Probably that would require putting units into
the whole transform framework.
This transform implements the linear transform part of the WCS FITS
standard outlined in Greisen & Calabata 2002 (A&A in press; find it at
"http://arxiv.org/abs/astro-ph/0207407").
As a special case, you can pass in the boolean option "ignore_rgb"
(default 0), and if you pass in a 3-D FITS header in which the last
dimension has exactly 3 elements, it will be ignored in the output
transformation. That turns out to be handy for handling rgb images.
=cut
sub t_fits {
my($class) = 'PDL::Transform::Linear';
my($hdr) = shift;
my($opt) = shift;
if(ref $opt ne 'HASH') {
$opt = defined $opt ? {$opt,@_} : {} ;
}
$hdr = $hdr->gethdr
if(defined $hdr && UNIVERSAL::isa($hdr,'PDL'));
croak('PDL::Transform::FITS::new requires a FITS header hash\n')
if(!defined $hdr || ref $hdr ne 'HASH' || !defined($hdr->{NAXIS}));
my($n) = $hdr->{NAXIS}; $n = $n->at(0) if(UNIVERSAL::isa($n,'PDL'));
$n = 2
if($opt->{ignore_rgb} && $n==3 && $hdr->{NAXIS3} == 3);
my($matrix) = PDL->zeroes($hdr->{NAXIS});
my($pre) = PDL->zeroes($n);
my($post) = PDL->zeroes($n);
##############################
# Scaling: Use CDi_j formalism if present (mostly in Hubble
# datasets); otherwise use CPi_j + CDELTi formalism.
my(@hgrab);
if(@hgrab = grep(m/^CD\d{1,3}_\d{1,3}/,keys %$hdr)) { # assignment
#
# CDi_j formalism
#
for my $h(@hgrab) {
$h =~ m/CD(\d{1,3})_(\d{1,3})/; # Should always match
$matrix->(($1),($2)) .= $hdr->{h};
}
print "PDL::Transform::FITS: Detected CDi_j matrix: \n",$matrix,"\n"
if($PDL::Transform::debug);
} else {
#
# CPi_j + CDELTi formalism
# If CPi_j arent present, and N=2, then try using CROTA or
# CROTA1 to generate a rotation matrix instea.
#
my($cdm) = PDL->zeroes($n,$n);
my($cd) = $cdm->diagonal(0,1);
my($cpm) = PDL->zeroes($n,$n);
$cpm->diagonal(0,1) .= 1; # CP: diagonal defaults to unity
$cd .= 1;
if( @hgrab = grep(m/^CP\d{1,3}_\d{1,3}/,keys %$hdr) ) { # assignment
for my $h(@hgrab) {
$h =~ m/CP(\d{1,3})_(\d{1,3})/; # Should always match
$cpm->(($1),($2)) .= $hdr->{h};
}
print "PDL::Transform::FITS: Detected CPi_j matrix: \n",$cpm,"\n"
if($PDL::Transform::debug && @hgrab);
} elsif($n==2 && ( defined $hdr->{CROTA} || defined $hdr->{CROTA1} ) ) {
my $cr = $hdr->{CROTA};
$cr = $hdr->{CROTA1} unless defined $cr;
$cr *= $DEG2RAD;
$cpm .= pdl( [cos($cr), sin($cr)],[-sin($cr),cos($cr)] );
}
for my $i(1..$n) {
$cd->(($i-1)) .= $hdr->{"CDELT$i"};
}
$matrix = $cdm x $cpm
}
my($i1) = 0;
for my $i(1..$n) {
$pre->(($i1)) .= 1 - $hdr->{"CRPIX$i"};
$post->(($i1)).= $hdr->{"CRVAL$i"};
$i1++;
}
my($me) = PDL::Transform::Linear::new($class,
{'pre'=>$pre,
'post'=>$post,
'matrix'=>$matrix
});
$me->{name} = 'FITS';
my (@otype,@ounit,@itype,@iunit);
our (@names) = ('X','Y','Z') unless defined(@names);
for my $i(1..$hdr->{NAXIS}) {
push(@otype,$hdr->{"CTYPE$i"});
push(@ounit,$hdr->{"CUNIT$i"});
push(@itype,"Image ". ( ($i<$#names) ? $names[$i] : "${i}th dim" ));
push(@iunit,"Pixels");
}
$me->{otype} = \@otype;
$me->{itype} = \@itype;
$me->{ounit} = \@ounit;
$me->{iunit} = \@iunit;
return $me;
}
======EOD_t_fits======
pp_add_exported('t_code ');
pp_addpm(<<'======EOD_t_code ======');
=head2 t_code
=for usage
$f = t_code(,[],[options]);
=for ref
Transform implementing arbitrary perl code.
This is a way of getting quick-and-dirty new transforms. You pass in
anonymous (or otherwise) code refs pointing to subroutines that
implement the forward and, optionally, inverse transforms. The
subroutines should accept a data PDL followed by a parameter hash ref,
and return the transformed data PDL. The parameter hash ref can be
set via the options, if you want to.
Options that are accepted are:
=over 3
=item p,params
The parameter hash that will be passed back to your code (defaults to the
empty hash).
=item n,name
The name of the transform (defaults to "code").
=item i, idim (default 2)
The number of input dimensions (additional ones should be passed through
unchanged)
=item o, odim (default 2)
The number of output dimensions
=item itype
The type of the input dimensions, in an array ref (optional and advisiory)
=item otype
The type of the output dimension, in an array ref (optional and advisory)
=item iunit
The units that are expected for the input dimensions (optional and advisory)
=item ounit
The units that are returned in the output (optional and advisory).
=back
The code variables are executable perl code, either as a code ref or
as a string that will be eval'ed to produce code refs. If you pass in
a string, it gets eval'ed at call time to get a code ref. If it compiles
OK but does not return a code ref, then it gets re-evaluated with "sub {
... }" wrapped around it, to get a code ref.
Note that code callbacks like this can be used to do really weird
things and generate equally weird results -- caveat scriptor!
=cut
sub t_code {
my($class) = 'PDL::Transform';
my($func, $inv, $o) = @_;
if(ref $inv eq 'HASH') {
$o = $inv;
$inv = undef;
}
my($me) = PDL::Transform::new($class);
$me->{name} = _opt($o,['n','name','Name']) || "code";
$me->{func} = $func;
$me->{inv} = $inv;
$me->{params} = _opt($o,['p','params','Params']) || {};
$me->{idim} = _opt($o,['i','idim']) || 2;
$me->{odim} = _opt($o,['o','odim']) || 2;
$me->{itype} = _opt($o,['itype']) || [];
$me->{otype} = _opt($o,['otype']) || [];
$me->{iunit} = _opt($o,['iunit']) || [];
$me->{ounit} = _opt($o,['ounit']) || [];
$me;
}
======EOD_t_code ======
pp_add_exported('t_cylindrical');
pp_add_exported('t_radial');
pp_addpm(<<'======EOD_t_cylindrical======');
=head2 t_cylindrical
=head2 t_radial
=for usage
$f = t_radial();
=for ref
Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with inverse)
Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose
direct or conformal conversion. Direct conversion preserves radial
distance from the origin; conformal conversion preserves local angles,
so that each small-enough part of the image only appears to be scaled
and rotated, not stretched. Conformal conversion puts the radius on a
logarithmic scale, so that scaling of the original image plane is
equivalent to a simple offset of the transformed image plane.
If you use three or more dimensions, the higher dimensions are ignored,
yielding a conversion from Cartesian to cylindrical coordinates, which
is why there are two aliases for the same transform. If you use higher
dimensionality than 2, you must manually specify the origin or you will
get dimension mismatch errors when you apply the transform.
Theta runs B instead of the more usual counterclockwise; that is
to preserve the mirror sense of small structures.
OPTIONS:
=over 3
=item d, direct, Direct
Generate (theta,r) coordinates out (this is the default); incompatible
with Conformal. Theta is in radians, and the radial coordinate is
in the units of distance in the input plane.
=item r0, c, conformal, Conformal
If defined, this floating-point value causes t_radial to generate
(theta, ln(r/r0)) coordinates out. Theta is in radians, and the
radial coordinate varies by 1 for each e-folding of the r0-scaled
distance from the input origin. The logarithmic scaling is useful for
viewing both large and small things at the same time, and for keeping
shapes of small things preserved in the image.
=item o, origin, Origin [default (0,0,0)]
This is the origin of the expansion. Pass in a PDL or an array ref.
=item u, unit, Unit [default 'radians']
This is the angular unit to be used for the azimuth.
=back
EXAMPLES
These examples do transformations back into the same size image as they
started from; by suitable use of the "transform" option to
L you can send them to any size array you like.
Examine radial structure in M51:
Here, we scale the output to stretch 2*pi radians out to the
full image width in the horizontal direction, and to stretch 1 radius out
to a diameter in the vertical direction.
$a = rfits('m51.fits');
$ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image
$tu = t_radial(o => [130,130]); # Expand around galactic core
$b = $a->map($ts x $tu);
Examine radial structure in M51 (conformal):
Here, we scale the output to stretch 2*pi radians out to the full image width
in the horizontal direction, and scale the vertical direction by the exact
same amount to preserve conformality of the operation. Notice that
each piece of the image looks "natural" -- only scaled and not stretched.
$a = rfits('m51.fits')
$ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale.
$tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image
$b = $ts->compose($tu)->unmap($a);
=cut
*t_cylindrical = \&t_radial;
sub t_radial {
my($class) = 'PDL::Transform';
my($o) = $_[0];
if(ref $o ne 'HASH') {
$o = { @_ };
}
my($me) = PDL::Transform::new($class);
$me->{params}->{origin} = _opt($o,['o','origin','Origin']);
$me->{params}->{origin} = pdl(0,0)
unless defined($me->{params}->{origin});
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin});
$me->{params}->{r0} = _opt($o,['r0','R0','c','conformal','Conformal']);
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin});
$me->{params}->{u} = _opt($o,['u','unit','Unit'],'radians');
### Replace this kludge with a units call
$me->{params}->{angunit} = ($me->{params}->{u} =~ m/^d/i) ? $RAD2DEG : 1.0;
print "radial: conversion is $me->{params}->{angunit}\n" if($PDL::Transform::debug);
$me->{name} = "radial (direct)";
$me->{idim} = 2;
$me->{odim} = 2;
if($me->{params}->{r0}) {
$me->{otype} = ["Azimuth", "Ln radius" . ($me->{params}->{r0} != 1.0 ? "/$me->{params}->{r0}" : "")];
$me->{ounit} = [$me->{params}->{u},'']; # true-but-null prevents copying
} else {
$me->{otype} = ["Azimuth","Radius"];
$me->{ounit} = [$me->{params}->{u},'']; # false value copies prev. unit
}
$me->{func} = sub {
my($data,$o) = @_;
my($out) = ($data->new_or_inplace);
my($d) = $data->copy;
$d->(0:1) -= $o->{origin};
my($d0) = $d->((0));
my($d1) = $d->((1));
# (mod operator on atan2 puts everything in the interval [0,2*PI).)
$out->((0)) .= (atan2(-$d1,$d0) % (2*$PI)) * $me->{params}->{angunit};
$out->((1)) .= (defined $o->{r0}) ?
0.5 * log( ($d1*$d1 + $d0 * $d0) / ($o->{r0} * $o->{r0}) ) :
sqrt($d1*$d1 + $d0*$d0);
$out;
};
$me->{inv} = sub {
my($d,$o) = @_;
my($d0,$d1,$out)=
( ($d->is_inplace) ?
($d->((0))->copy, $d->((1))->copy->dummy(0,2), $d) :
($d->((0)), $d->((1))->dummy(0,2), $d->copy)
);
$d0 /= $me->{params}->{angunit};
my($os) = $out->(0:1);
$os .= append(cos($d0)->dummy(0,1),-sin($d0)->dummy(0,1));
$os *= defined $o->{r0} ? ($o->{r0} * exp($d1)) : $d1;
$os += $o->{origin};
$out;
};
$me;
}
======EOD_t_cylindrical======
pp_add_exported('t_quadratic');
pp_addpm(<<'======EOD_t_quadratic======');
=head2 t_quadratic
=for usage
$t = t_quadratic();
=for ref
Quadratic scaling -- cylindrical pincushion (n-d; with inverse)
Quadratic scaling emulates pincushion in a cylindrical optical system:
separate quadratic scaling is applied to each axis. You can apply
separate distortion along any of the principal axes. If you want
different axes, use L and L to rotate
them to the correct angle. The scaling options may be scalars or
vectors; if they are scalars then the expansion is isotropic.
The formula for the expansion is:
f(a) = ( + * a^2/ ) / (abs() + 1)
where is a scaling coefficient and is a fundamental
length scale. Negative values of result in a pincushion
contraction.
OPTIONS
=over 3
=item o,origin,Origin
The origin of the pincushion. (default is the, er, origin).
=item l,l0,length,Length,r0
The fundamental scale of the transformation -- the radius that remains
unchanged. (default=1)
=item s,str,strength,Strength
The relative strength of the pincushion. (default = 0.1)
=item d, dim, dims, Dims
The number of dimensions to quadratically scale (default is the
dimensionality of your input vectors)
=back
=cut
sub t_quadratic {
my($class) = 'PDL::Transform';
my($o) = $_[0];
if(ref $o ne 'HASH') {
$o = {@_};
}
my($me) = PDL::Transform::new($class);
$me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0));
$me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1));
$me->{params}->{str} = _opt($o,['s','str','strength','Strength'],pdl(0.1));
$me->{params}->{dim} = _opt($o,['d','dim','dims','Dims']);
$me->{name} = "quadratic";
$me->{func} = sub {
my($data,$o) = @_;
my($dd) = $data->copy - $o->{origin};
my($d) = (defined $o->{dim}) ? $dd->(0:($o->{dim}-1)) : $dd;
$d += $o->{str} * ($d * abs($d)) / $o->{l0};
$d /= (abs($o->{str}) + 1);
$d += $o->{origin};
if($data->is_inplace) {
$data .= $dd;
return $data;
}
$dd;
};
$me->{inv} = sub {
my($data,$opt) = @_;
my($dd) = $data->copy ;
my($d) = (defined $opt->{dim}) ? $dd->(0:($opt->{dim}-1)) : $dd;
my($o) = $opt->{origin};
my($s) = $opt->{str};
my($l) = $opt->{l0};
$d .= ((-1 + sqrt(1 + 4 * $s/$l * abs($d-$o) * (1+abs($s))))
/ 2 / $s * $l) * (1 - 2*($d < $o));
$d += $o;
if($data->is_inplace) {
$data .= $dd;
return $data;
}
$dd;
};
$me;
}
======EOD_t_quadratic======
pp_add_exported('t_spherical');
pp_addpm(<<'======EOD_t_spherical======');
=head2 t_spherical
=for usage
$t = t_spherical();
=for ref
Convert Cartesian to spherical coordinates. (3-D; with inverse)
Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta
is longitude, centered on 0, and phi is latitude, also centered on 0.
Unless you specify Euler angles, the pole points in the +Z direction
and the prime meridian is in the +X direction. The default is for
theta and phi to be in radians; you can select degrees if you want
them.
Just as the L 2-D transform acts like a 3-D
cylindrical transform by ignoring third and higher dimensions,
Spherical acts like a hypercylindrical transform in four (or higher)
dimensions. Also as with L, you must manually specify
the origin if you want to use more dimensions than 3.
To deal with latitude & longitude on the surface of a sphere (rather than
full 3-D coordinates), see L.
OPTIONS:
=over 3
=item o, origin, Origin [default (0,0,0)]
This is the Cartesian origin of the spherical expansion. Pass in a PDL
or an array ref.
=item e, euler, Euler [default (0,0,0)]
This is a 3-vector containing Euler angles to change the angle of the
pole and ordinate. The first two numbers are the (theta, phi) angles
of the pole in a (+Z,+X) spherical expansion, and the last is the
angle that the new prime meridian makes with the meridian of a simply
tilted sphere. This is implemented by composing the output transform
with a PDL::Transform::Linear object.
=item u, unit, Unit (default radians)
This option sets the angular unit to be used. Acceptable values are
"degrees","radians", or reasonable substrings thereof (e.g. "deg", and
"rad", but "d" and "r" are deprecated). Once genuine unit processing
comes online (a la Math::Units) any angular unit should be OK.
=back
=cut
sub t_spherical {
my($class) = 'PDL::Transform';
my($o) = $_[0];
if(ref $o ne 'HASH') {
$o = { @_ } ;
}
my($me) = PDL::Transform::new($class);
$me->{params}->{idim} = 3;
$me->{params}->{odim} = 3;
$me->{params}->{origin} = _opt($o,['o','origin','Origin']);
$me->{params}->{origin} = PDL->zeroes(3)
unless defined($me->{params}->{origin});
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin});
$me->{params}->{deg} = _opt($o,['d','degrees','Degrees']);
my $unit = _opt($o,['u','unit','Unit']);
$me->{params}->{angunit} = ($unit =~ m/^d/i) ?
$DEG2RAD : undef;
$me->{name} = "spherical";
$me->{func} = sub {
my($data,$o) = @_;
my($d) = $data->copy - $o->{origin};
my($d0,$d1,$d2) = ($d->((0)),$d->((1)),$d->((2)));
my($out) = ($d->is_inplace) ? $data : $data->copy;
$out->((0)) .= sqrt($d0*$d0 + $d1*$d1 + $d2*$d2);
$out->((1)) .= atan2($d1, $d0);
$out->((2)) .= asin($d2 / $out->((0)));
$out->(1:2) /= $o->{angunit}
if(defined $o->{angunit});
$out;
};
$me->{inv} = sub {
my($d,$o) = @_;
my($theta,$phi,$r,$out) =
( ($d->is_inplace) ?
($d->((0))->copy, $d->((1))->copy, $d->((2))->copy, $d) :
($d->((0)), $d->((1)), $d->((2)), $d->copy)
);
my($x,$y,$z) =
($out->((0)),$out->((1)),$out->((2)));
my($ph,$th);
if(defined $o->{angunit}){
$ph = $o->{angunit} * $phi;
$th = $o->{angunit} * $theta;
} else {
$ph = $phi;
$th = $theta;
}
$z .= $r * sin($ph);
$x .= $r * cos($ph);
$y .= $x * sin($th);
$x *= cos($th);
$out += $o->{origin};
$out;
};
$me;
}
======EOD_t_spherical======
pp_done();