"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Libtmp/Transform/transform.pd" between
PDL-2.077.tar.gz and PDL-2.078.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

transform.pd  (PDL-2.077):transform.pd  (PDL-2.078)
skipping to change at line 1560 skipping to change at line 1560
C<map()> supports bad values in the data array. Bad values in the input C<map()> supports bad values in the data array. Bad values in the input
array are propagated to the output array. The 'g' and 'h' methods will array are propagated to the output array. The 'g' and 'h' methods will
do some smoothing over bad values: if more than 1/3 of the weighted do some smoothing over bad values: if more than 1/3 of the weighted
input-array footprint of a given output pixel is bad, then the output input-array footprint of a given output pixel is bad, then the output
pixel gets marked bad. pixel gets marked bad.
=cut =cut
+==EOD_map_doc== +==EOD_map_doc==
PMCode => <<'+==EOD_map_perlcode==' PMCode => pp_line_numbers(__LINE__, <<'+==EOD_map_perlcode=='),
sub PDL::match { sub PDL::match {
# Set default for rectification to 0 for simple matching... # Set default for rectification to 0 for simple matching...
if( ref($_[$#_]) ne 'HASH' ) { push @_, {} if ref($_[-1]) ne 'HASH';
push(@_,{})
}
my @k = grep(m/^r(e(c(t)?)?)?/,sort keys %{$_[-1]}); my @k = grep(m/^r(e(c(t)?)?)?/,sort keys %{$_[-1]});
unless(@k) { unless(@k) {
$_[-1]->{rectify} = 0; $_[-1]->{rectify} = 0;
} }
t_identity()->map(@_); t_identity()->map(@_);
} }
*PDL::map = \&map; *PDL::map = \&map;
sub map { sub map {
my($me) = shift; my ($me, $in, $tmp, $opt) = @_;
my($in) = shift; ($in, $me) = ($me, $in)
if UNIVERSAL::isa($me,'PDL') && UNIVERSAL::isa($in,'PDL::Transform');
if(UNIVERSAL::isa($me,'PDL') && UNIVERSAL::isa($in,'PDL::Transform')) {
my($x) = $in;
$in = $me;
$me = $x;
}
barf ("PDL::Transform::map: source is not defined or is not a PDL\n") barf ("PDL::Transform::map: source is not defined or is not a PDL\n")
unless(defined $in and UNIVERSAL::isa($in,'PDL')); unless(defined $in and UNIVERSAL::isa($in,'PDL'));
barf ("PDL::Transform::map: source is empty\n") barf ("PDL::Transform::map: source is empty\n")
unless($in->nelem); unless($in->nelem);
my($tmp) = shift;
my($opt) = shift;
# Check for options-but-no-template case # Check for options-but-no-template case
if(ref $tmp eq 'HASH' && !(defined $opt)) { ($opt, $tmp) = ($tmp, undef)
if(!defined($tmp->{NAXIS})) { # FITS headers all have NAXIS. if ref $tmp eq 'HASH' && !defined($opt)
$opt = $tmp; && !defined($tmp->{NAXIS}); # FITS headers all have NAXIS.
$tmp = undef;
}
}
croak("PDL::Transform::map: Option 'p' was ambiguous and has been removed. You probably want 'pix' or 'phot'.") if exists($opt->{'p'}); croak("PDL::Transform::map: Option 'p' was ambiguous and has been removed. You probably want 'pix' or 'phot'.") if exists($opt->{'p'});
$tmp = [$in->dims] unless defined $tmp; # no //= because of Devel::Cover bug
$tmp = [$in->dims] unless(defined($tmp));
# Generate an appropriate output ndarray for values to go in # Generate an appropriate output ndarray for values to go in
my($out); my ($out, @odims, $ohdr);
my(@odims);
my($ohdr);
if(UNIVERSAL::isa($tmp,'PDL')) { if(UNIVERSAL::isa($tmp,'PDL')) {
@odims = $tmp->dims; @odims = $tmp->dims;
my($x); my($x);
if(defined ($x = $tmp->gethdr)) { if(defined ($x = $tmp->gethdr)) {
my(%b) = %{$x}; my(%b) = %{$x};
$ohdr = \%b; $ohdr = \%b;
} }
} elsif(ref $tmp eq 'HASH') { } elsif(ref $tmp eq 'HASH') {
# (must be a fits header -- or would be filtered above) # (must be a fits header -- or would be filtered above)
for my $i(1..$tmp->{NAXIS}){ for my $i(1..$tmp->{NAXIS}){
push(@odims,$tmp->{"NAXIS$i"}); push(@odims,$tmp->{"NAXIS$i"});
} }
skipping to change at line 1872 skipping to change at line 1851
grep /(^CROTA\d*$)|(^(PC)\d+_\d+[A-Z]?$)|(CDELT\d*$)/, keys %{$out ->hdr} grep /(^CROTA\d*$)|(^(PC)\d+_\d+[A-Z]?$)|(CDELT\d*$)/, keys %{$out ->hdr}
}; };
} }
} }
$out->hdrcpy(1); $out->hdrcpy(1);
############################## ##############################
# Sandwich the transform between the input and output plane FITS headers. # Sandwich the transform between the input and output plane FITS headers.
unless($nofits) { unless($nofits) {
$me = !(t_fits($out,{ignore_rgb=>1})) x $me x $f_in; $me = !t_fits($out,{ignore_rgb=>1}) x $me x $f_in;
} }
############################## ##############################
## Figure out the interpND options ## Figure out the interpND options
my $method = _opt($opt,['m','method','Method'],undef); my $method = _opt($opt,['m','method','Method'],undef);
my $bound = _opt($opt,['b','bound','boundary','Boundary'],'t'); my $bound = _opt($opt,['b','bound','boundary','Boundary'],'t');
############################## ##############################
## Rubber meets the road: calculate the inverse transformed points. ## Rubber meets the road: calculate the inverse transformed points.
my $ndc = PDL::Basic::ndcoords(@dd); my $ndc = PDL::Basic::ndcoords(@dd);
skipping to change at line 2054 skipping to change at line 2033
The inverse transform remains connected to the main transform because The inverse transform remains connected to the main transform because
they both point to the original parameters hash. That turns out to be they both point to the original parameters hash. That turns out to be
useful. useful.
=cut =cut
*t_inverse = \&inverse; *t_inverse = \&inverse;
sub inverse { sub inverse {
my($me) = shift; my($me) = shift;
barf("PDL::Transform::inverse: got a transform with no inverse")
unless(defined($me->{inv})) { unless defined $me->{inv};
Carp::cluck("PDL::Transform::inverse: got a transform with no inverse.\n"); bless {
return undef; %$me, # force explicit copy of top-level
} inv => $me->{func},
func => $me->{inv},
my(%out) = %$me; # force explicit copy of top-level (map +("o$_"=>$me->{"i$_"}, "i$_"=>$me->{"o$_"}), qw(dim type unit)),
my($out) = \%out; name => "(inverse ".$me->{name}.")",
is_inverse => !($me->{is_inverse}),
$out->{inv} = $me->{func}; }, ref $me;
$out->{func} = $me->{inv};
$out->{idim} = $me->{odim};
$out->{odim} = $me->{idim};
$out->{otype} = $me->{itype};
$out->{itype} = $me->{otype};
$out->{ounit} = $me->{iunit};
$out->{iunit} = $me->{ounit};
$out->{name} = "(inverse ".$me->{name}.")";
$out->{is_inverse} = !($out->{is_inverse});
bless $out,(ref $me);
return $out;
} }
+======EOD_t_inverse====== +======EOD_t_inverse======
pp_add_exported('t_compose'); pp_add_exported('t_compose');
pp_addpm(<<'+======EOD_t_compose======'); pp_addpm(<<'+======EOD_t_compose======');
=head2 t_compose =head2 t_compose
=for usage =for usage
skipping to change at line 2120 skipping to change at line 2082
No checking is done that the itype/otype and iunit/ounit fields are No checking is done that the itype/otype and iunit/ounit fields are
compatible -- that may happen later, or you can implement it yourself compatible -- that may happen later, or you can implement it yourself
if you like. if you like.
=cut =cut
@PDL::Transform::Composition::ISA = ('PDL::Transform'); @PDL::Transform::Composition::ISA = ('PDL::Transform');
sub PDL::Transform::Composition::stringify { sub PDL::Transform::Composition::stringify {
package PDL::Transform::Composition; package PDL::Transform::Composition;
my($me) = shift; shift->SUPER::stringify;
my($out) = SUPER::stringify $me;
$out;
} }
*t_compose = \&compose; *t_compose = \&compose;
sub compose { sub compose {
local($_); local($_);
my(@funcs) = @_; my(@funcs) = @_;
my($me) = PDL::Transform->new; my($me) = PDL::Transform->new;
# No inputs case: return the identity function # No inputs case: return the identity function
return $me return $me if !@funcs;
if(!@funcs);
$me->{name} = ""; $me->{name} = "";
my($f); my @clist;
my(@clist); for my $f (@funcs) {
for $f(@funcs) {
$me->{idim} = $f->{idim} if($f->{idim} > $me->{idim}); $me->{idim} = $f->{idim} if($f->{idim} > $me->{idim});
$me->{odim} = $f->{odim} if($f->{odim} > $me->{odim}); $me->{odim} = $f->{odim} if($f->{odim} > $me->{odim});
if(UNIVERSAL::isa($f,"PDL::Transform::Composition")) { if(UNIVERSAL::isa($f,"PDL::Transform::Composition")) {
if($f->{is_inverse}) { if($f->{is_inverse}) {
for(reverse(@{$f->{params}->{clist}})) { for(reverse(@{$f->{params}->{clist}})) {
push(@clist,$_->inverse); push(@clist,$_->inverse);
$me->{name} .= " o inverse ( ".$_->{name}." )"; $me->{name} .= " o inverse ( ".$_->{name}." )";
} }
} else { } else {
for(@{$f->{params}->{clist}}) { for(@{$f->{params}->{clist}}) {
push(@clist,$_); push(@clist,$_);
$me->{name} .= " o ".$_->{name}; $me->{name} .= " o ".$_->{name};
} }
} }
} else { # Not a composition -- just push the transform onto the list. } else { # Not a composition -- just push the transform onto the list.
push(@clist,$f); push(@clist,$f);
$me->{name} .= " o ".$f->{name}; $me->{name} .= " o ".$f->{name};
} }
} }
$me->{name}=~ s/^ o //; # Get rid of leading composition mark $me->{name}=~ s/^ o //; # Get rid of leading composition mark
$me->{otype} = $funcs[0]->{otype}; $me->{otype} = $funcs[0]->{otype};
$me->{ounit} = $funcs[0]->{ounit}; $me->{ounit} = $funcs[0]->{ounit};
$me->{itype} = $funcs[-1]->{itype}; $me->{itype} = $funcs[-1]->{itype};
$me->{iunit} = $funcs[-1]->{iunit}; $me->{iunit} = $funcs[-1]->{iunit};
$me->{params}->{clist} = \@clist; $me->{params}->{clist} = \@clist;
$me->{func} = sub { $me->{func} = sub {
my ($data,$p) = @_; my ($data,$p) = @_;
my ($ip) = $data->is_inplace; my ($ip) = $data->is_inplace;
for my $t ( reverse @{$p->{clist}} ) { for my $t ( reverse @{$p->{clist}} ) {
croak("Error: tried to apply a PDL::Transform with no function inside a co mposition!\n offending transform: $t\n") croak("Error: tried to apply a PDL::Transform with no function inside a co mposition!\n offending transform: $t\n")
unless(defined($t->{func}) and ref($t->{func}) eq 'CODE'); unless(defined($t->{func}) and ref($t->{func}) eq 'CODE');
$data = $t->{func}($ip ? $data->inplace : $data, $t->{params}); $data = $t->{func}($ip ? $data->inplace : $data, $t->{params});
} }
$data->is_inplace(0); # clear inplace flag (avoid core bug with inplace) $data->is_inplace(0); # clear inplace flag (avoid core bug with inplace)
$data; $data;
}; };
$me->{inv} = sub { $me->{inv} = sub {
my($data,$p) = @_; my($data,$p) = @_;
my($ip) = $data->is_inplace; my($ip) = $data->is_inplace;
for my $t ( @{$p->{clist}} ) { for my $t ( @{$p->{clist}} ) {
croak("Error: tried to invert a non-invertible PDL::Transform inside a com position!\n offending transform: $t\n") croak("Error: tried to invert a non-invertible PDL::Transform inside a com position!\n offending transform: $t\n")
unless(defined($t->{inv}) and ref($t->{inv}) eq 'CODE'); unless(defined($t->{inv}) and ref($t->{inv}) eq 'CODE');
$data = &{$t->{inv}}($ip ? $data->inplace : $data, $t->{params}); $data = &{$t->{inv}}($ip ? $data->inplace : $data, $t->{params});
} }
$data; $data;
}; };
return bless($me,'PDL::Transform::Composition'); return bless($me,'PDL::Transform::Composition');
} }
+======EOD_t_compose====== +======EOD_t_compose======
pp_add_exported('t_wrap'); pp_add_exported('t_wrap');
pp_addpm(<<'+======EOD_t_wrap======'); pp_addpm(<<'+======EOD_t_wrap======');
=head2 t_wrap =head2 t_wrap
=for usage =for usage
$g1fg = $f->wrap($g); $g1fg = $f->wrap($g);
skipping to change at line 2686 skipping to change at line 2630
#last is undef #last is undef
if(!(ref $o)) { if(!(ref $o)) {
$o = {@_}; $o = {@_};
} }
my($me) = PDL::Transform::new($class); my($me) = PDL::Transform::new($class);
$me->{name} = "linear"; $me->{name} = "linear";
$me->{params}->{pre} = _opt($o,['pre','Pre','preoffset','offset', $me->{params}{pre} = _opt($o,['pre','Pre','preoffset','offset',
'Offset','PreOffset','Preoffset'],0); 'Offset','PreOffset','Preoffset'],0);
$me->{params}->{pre} = pdl($me->{params}->{pre}) $me->{params}{pre} = pdl($me->{params}{pre})
if(defined $me->{params}->{pre}); if(defined $me->{params}{pre});
$me->{params}->{post} = _opt($o,['post','Post','postoffset','PostOffset', $me->{params}{post} = _opt($o,['post','Post','postoffset','PostOffset',
'shift','Shift'],0); 'shift','Shift'],0);
$me->{params}->{post} = pdl($me->{params}->{post}) $me->{params}{post} = pdl($me->{params}{post})
if(defined $me->{params}->{post}); if(defined $me->{params}{post});
$me->{params}->{matrix} = _opt($o,['m','matrix','Matrix','mat','Mat']); $me->{params}{matrix} = _opt($o,['m','matrix','Matrix','mat','Mat']);
$me->{params}->{matrix} = pdl($me->{params}->{matrix}) $me->{params}{matrix} = pdl($me->{params}{matrix})
if(defined $me->{params}->{matrix}); if(defined $me->{params}{matrix});
$me->{params}->{rot} = _opt($o,['r','rot','rota','rotation','Rotation']); $me->{params}{rot} = _opt($o,['r','rot','rota','rotation','Rotation'], 0);
$me->{params}->{rot} = 0 unless defined($me->{params}->{rot}); $me->{params}{rot} = pdl($me->{params}{rot});
$me->{params}->{rot} = pdl($me->{params}->{rot});
my $o_dims = _opt($o,['d','dim','dims','Dims']); my $o_dims = _opt($o,['d','dim','dims','Dims']);
$o_dims = pdl($o_dims) $o_dims = pdl($o_dims) if defined($o_dims);
if defined($o_dims);
my $scale = _opt($o,['s','scale','Scale']); my $scale = _opt($o,['s','scale','Scale']);
$scale = pdl($scale) $scale = pdl($scale) if defined($scale);
if defined($scale);
# Figure out the number of dimensions to transform, and, # Figure out the number of dimensions to transform, and,
# if necessary, generate a new matrix. # if necessary, generate a new matrix.
if(defined($me->{params}->{matrix})) { if(defined($me->{params}{matrix})) {
my $mat = $me->{params}->{matrix} = $me->{params}->{matrix}->slice(":,:"); my $mat = $me->{params}{matrix} = $me->{params}{matrix}->slice(":,:");
$me->{idim} = $mat->dim(0); $me->{idim} = $mat->dim(0);
$me->{odim} = $mat->dim(1); $me->{odim} = $mat->dim(1);
} else { } else {
if(defined($me->{params}->{rot}) && if(defined($me->{params}->{rot}) &&
UNIVERSAL::isa($me->{params}->{rot},'PDL')) { UNIVERSAL::isa($me->{params}->{rot},'PDL')) {
$me->{idim} = $me->{odim} = 2 if($me->{params}->{rot}->nelem == 1); $me->{idim} = $me->{odim} = 2 if($me->{params}{rot}->nelem == 1);
$me->{idim} = $me->{odim} = 3 if($me->{params}->{rot}->nelem == 3); $me->{idim} = $me->{odim} = 3 if($me->{params}{rot}->nelem == 3);
} }
if(defined($scale) && if(defined($scale) &&
UNIVERSAL::isa($scale,'PDL') && UNIVERSAL::isa($scale,'PDL') &&
$scale->getndims > 0) { $scale->getndims > 0) {
$me->{idim} = $me->{odim} = $scale->dim(0); $me->{idim} = $me->{odim} = $scale->dim(0);
$me->{odim} = $scale->dim(0); $me->{odim} = $scale->dim(0);
} elsif(defined($me->{params}->{pre}) && } elsif(defined($me->{params}{pre}) &&
UNIVERSAL::isa($me->{params}->{pre},'PDL') && UNIVERSAL::isa($me->{params}{pre},'PDL') &&
$me->{params}->{pre}->getndims > 0) { $me->{params}{pre}->getndims > 0) {
$me->{idim} = $me->{odim} = $me->{params}->{pre}->dim(0); $me->{idim} = $me->{odim} = $me->{params}{pre}->dim(0);
} elsif(defined($me->{params}->{post}) && } elsif(defined($me->{params}{post}) &&
UNIVERSAL::isa($me->{params}->{post},'PDL') && UNIVERSAL::isa($me->{params}{post},'PDL') &&
$me->{params}->{post}->getndims > 0) { $me->{params}{post}->getndims > 0) {
$me->{idim} = $me->{odim} = $me->{params}->{post}->dim(0); $me->{idim} = $me->{odim} = $me->{params}{post}->dim(0);
} elsif(defined($o_dims)) { } elsif(defined($o_dims)) {
$me->{idim} = $me->{odim} = $o_dims; $me->{idim} = $me->{odim} = $o_dims;
} else { } else {
print "PDL::Transform::Linear: assuming 2-D transform (set dims option to change)\n" if($PDL::Transform::debug); print "PDL::Transform::Linear: assuming 2-D transform (set dims option to change)\n" if($PDL::Transform::debug);
$me->{idim} = $me->{odim} = 2; $me->{idim} = $me->{odim} = 2;
} }
$me->{params}->{matrix} = PDL->zeroes($me->{idim},$me->{odim}); $me->{params}->{matrix} = PDL->zeroes($me->{idim},$me->{odim});
my $tmp; # work around perl -d "feature" my $tmp; # work around perl -d "feature"
($tmp = $me->{params}->{matrix}->diagonal(0,1)) .= 1; ($tmp = $me->{params}->{matrix}->diagonal(0,1)) .= 1;
} }
### Handle rotation option ### Handle rotation option
my $rot = $me->{params}->{rot}; my $rot = $me->{params}{rot};
if(defined($rot)) { if(defined($rot)) {
# Subrotation closure -- rotates from axis $d->(0) --> $d->(1). # Subrotation closure -- rotates from axis $d->(0) --> $d->(1).
my $subrot = sub { my $subrot = sub {
my($d,$angle,$m)=@_; my($d,$angle,$m)=@_;
my($i) = identity($m->dim(0)); my($i) = identity($m->dim(0));
my($subm) = $i->dice($d,$d); my($subm) = $i->dice($d,$d);
$angle = $angle->at(0) $angle = $angle->at(0)
if(UNIVERSAL::isa($angle,'PDL')); if(UNIVERSAL::isa($angle,'PDL'));
my($x) = $angle * $PDL::Transform::DEG2RAD; my($x) = $angle * $PDL::Transform::DEG2RAD;
$subm .= $subm x pdl([cos($x),-sin($x)],[sin($x),cos($x)] ); $subm .= $subm x pdl([cos($x),-sin($x)],[sin($x),cos($x)] );
$m .= $m x $i; $m .= $m x $i;
}; };
if(UNIVERSAL::isa($rot,'PDL') && $rot->nelem > 1) { if(UNIVERSAL::isa($rot,'PDL') && $rot->nelem > 1) {
if($rot->ndims == 2) { if($rot->ndims == 2) {
$me->{params}->{matrix} x= $rot; $me->{params}{matrix} x= $rot;
} elsif($rot->nelem == 3) { } elsif($rot->nelem == 3) {
my $rm = identity(3); my $rm = identity(3);
# Do these in reverse order to make it more like # Do these in reverse order to make it more like
# function composition! # function composition!
&$subrot(pdl(0,1),$rot->at(2),$rm); &$subrot(pdl(0,1),$rot->at(2),$rm);
&$subrot(pdl(2,0),$rot->at(1),$rm); &$subrot(pdl(2,0),$rot->at(1),$rm);
&$subrot(pdl(1,2),$rot->at(0),$rm); &$subrot(pdl(1,2),$rot->at(0),$rm);
$me->{params}->{matrix} .= $me->{params}->{matrix} x $rm; $me->{params}{matrix} .= $me->{params}{matrix} x $rm;
} else { } else {
barf("PDL::Transform::Linear: Got a strange rot option -- giving up.\n") ; barf("PDL::Transform::Linear: Got a strange rot option -- giving up.\n") ;
} }
} else { } else {
if($rot != 0 and $me->{params}->{matrix}->dim(0)>1) { if($rot != 0 and $me->{params}{matrix}->dim(0)>1) {
&$subrot(pdl(0,1),$rot,$me->{params}->{matrix}); &$subrot(pdl(0,1),$rot,$me->{params}{matrix});
} }
} }
} }
#
# Apply scaling # Apply scaling
# $me->{params}{matrix} = $me->{params}{matrix}->slice(":,:");
$me->{params}->{matrix} = $me->{params}->{matrix}->slice(":,:");
my $tmp; # work around perl -d "feature" my $tmp; # work around perl -d "feature"
($tmp = $me->{params}->{matrix}->diagonal(0,1)) *= $scale ($tmp = $me->{params}{matrix}->diagonal(0,1)) *= $scale
if defined($scale); if defined($scale);
#
# Check for an inverse and apply it if possible # Check for an inverse and apply it if possible
#
my($o2); my($o2);
if($me->{params}->{matrix}->det($o2 = {lu=>undef})) { if($me->{params}{matrix}->det($o2 = {lu=>undef})) {
$me->{params}->{inverse} = $me->{params}->{matrix}->inv($o2); $me->{params}{inverse} = $me->{params}{matrix}->inv($o2);
} else { } else {
delete $me->{params}->{inverse}; delete $me->{params}{inverse};
} }
$me->{params}->{idim} = $me->{idim}; $me->{params}{idim} = $me->{idim};
$me->{params}->{odim} = $me->{odim}; $me->{params}{odim} = $me->{odim};
############################## ##############################
# The meat -- just shift, matrix-multiply, and shift again. # The meat -- just shift, matrix-multiply, and shift again.
$me->{func} = sub { $me->{func} = sub {
my($in,$opt) = @_; my($in,$opt) = @_;
my($d) = $opt->{matrix}->dim(0)-1; my($d) = $opt->{matrix}->dim(0)-1;
barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n") barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n")
if($in->dim(0) < $d); if($in->dim(0) < $d);
my $x = $in->slice("0:$d")->copy + $opt->{pre};
my($x) = $in->slice("0:$d")->copy + $opt->{pre}; my $out = $in->is_inplace ? $in : $in->copy;
my($out) = $in->is_inplace ? $in : $in->copy; $out->slice("0:$d") .= $x x $opt->{matrix} + $opt->{post};
my $tmp; # work around perl -d "feature"
($tmp = $out->slice("0:$d")) .= $x x $opt->{matrix} + $opt->{post};
return $out; return $out;
}; };
$me->{inv} = (defined $me->{params}->{inverse}) ? sub { $me->{inv} = (defined $me->{params}{inverse}) ? sub {
my($in,$opt) = @_; my($in,$opt) = @_;
my($d) = $opt->{inverse}->dim(0)-1; my($d) = $opt->{inverse}->dim(0)-1;
barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n") barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n")
if($in->dim(0) < $d); if($in->dim(0) < $d);
my($x) = $in->slice("0:$d")->copy - $opt->{post}; my($x) = $in->slice("0:$d")->copy - $opt->{post};
my($out) = $in->is_inplace ? $in : $in->copy; my($out) = $in->is_inplace ? $in : $in->copy;
my $tmp; # work around perl -d "feature" my $tmp; # work around perl -d "feature"
($tmp = $out->slice("0:$d")) .= $x x $opt->{inverse} - $opt->{pre}; ($tmp = $out->slice("0:$d")) .= $x x $opt->{inverse} - $opt->{pre};
$out; $out;
} : undef; } : undef;
return $me; return $me;
} }
sub PDL::Transform::Linear::stringify { sub PDL::Transform::Linear::stringify {
package PDL::Transform::Linear; package PDL::Transform::Linear;
my($me) = shift; my($out) = SUPER::stringify $me; my($me) = shift; my($out) = SUPER::stringify $me;
my $mp = $me->{params}; my $mp = $me->{params};
skipping to change at line 3016 skipping to change at line 2942
"http://arxiv.org/abs/astro-ph/0207407"). "http://arxiv.org/abs/astro-ph/0207407").
As a special case, you can pass in the boolean option "ignore_rgb" 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 (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 dimension has exactly 3 elements, it will be ignored in the output
transformation. That turns out to be handy for handling rgb images. transformation. That turns out to be handy for handling rgb images.
=cut =cut
sub t_fits { sub t_fits {
my($class) = 'PDL::Transform::Linear';
my($hdr) = shift; my($hdr) = shift;
my($opt) = shift; my($opt) = shift;
if(ref $opt ne 'HASH') { if(ref $opt ne 'HASH') {
$opt = defined $opt ? {$opt,@_} : {} ; $opt = defined $opt ? {$opt,@_} : {} ;
} }
$hdr = $hdr->gethdr $hdr = $hdr->gethdr
if(defined $hdr && UNIVERSAL::isa($hdr,'PDL')); if(defined $hdr && UNIVERSAL::isa($hdr,'PDL'));
croak('PDL::Transform::FITS::new requires a FITS header hash\n') croak('PDL::Transform::FITS::new requires a FITS header hash\n')
if(!defined $hdr || ref $hdr ne 'HASH' || !defined($hdr->{NAXIS})); if(!defined $hdr || ref $hdr ne 'HASH' || !defined($hdr->{NAXIS}));
my($n) = $hdr->{NAXIS}; $n = $n->at(0) if(UNIVERSAL::isa($n,'PDL')); my($n) = $hdr->{NAXIS}; $n = $n->at(0) if(UNIVERSAL::isa($n,'PDL'));
$n = 2 $n = 2
if($opt->{ignore_rgb} && $n==3 && $hdr->{NAXIS3} == 3); if($opt->{ignore_rgb} && $n==3 && $hdr->{NAXIS3} == 3);
my($matrix) = PDL->zeroes($hdr->{NAXIS},$hdr->{NAXIS}); my($matrix) = PDL->zeroes($n,$n);
my($pre) = PDL->zeroes($n); my($pre) = PDL->zeroes($n);
my($post) = PDL->zeroes($n); my($post) = PDL->zeroes($n);
############################## ##############################
# Scaling: Use CDi_j formalism if present; otherwise use the # Scaling: Use CDi_j formalism if present; otherwise use the
# older PCi_j + CDELTi formalism. # older PCi_j + CDELTi formalism.
my(@hgrab); my(@hgrab);
if(@hgrab = grep(m/^CD\d{1,3}_\d{1,3}$/,sort keys %$hdr)) { # assignment if(@hgrab = grep(m/^CD\d{1,3}_\d{1,3}$/,sort keys %$hdr)) { # assignment
skipping to change at line 3118 skipping to change at line 3043
} }
my($i1) = 0; my($i1) = 0;
for my $i(1..$n) { for my $i(1..$n) {
my $tmp; # work around perl -d "feature" my $tmp; # work around perl -d "feature"
($tmp = $pre->slice("($i1)")) .= 1 - $hdr->{"CRPIX$i"}; ($tmp = $pre->slice("($i1)")) .= 1 - $hdr->{"CRPIX$i"};
($tmp = $post->slice("($i1)")) .= $hdr->{"CRVAL$i"}; ($tmp = $post->slice("($i1)")) .= $hdr->{"CRVAL$i"};
$i1++; $i1++;
} }
my($me) = PDL::Transform::Linear::new($class, my $me = PDL::Transform::Linear->new({
{'pre'=>$pre, pre=>$pre, post=>$post, matrix=>$matrix
'post'=>$post, });
'matrix'=>$matrix
});
$me->{name} = 'FITS'; $me->{name} = 'FITS';
my (@otype,@ounit,@itype,@iunit); my (@otype,@ounit,@itype,@iunit);
our (@names) = ('X','Y','Z') unless @names; my @names = qw(X Y Z);
for my $i(1..$hdr->{NAXIS}) { for my $i(1..$hdr->{NAXIS}) {
push(@otype,$hdr->{"CTYPE$i"}); push(@otype,$hdr->{"CTYPE$i"});
push(@ounit,$hdr->{"CUNIT$i"}); push(@ounit,$hdr->{"CUNIT$i"});
push(@itype,"Image ". ( ($i-1<=$#names) ? $names[$i-1] : "${i}th dim" )); push(@itype,"Image ". ( ($i-1<=$#names) ? $names[$i-1] : "${i}th dim" ));
push(@iunit,"Pixels"); push(@iunit,"Pixels");
} }
$me->{otype} = \@otype; $me->{otype} = \@otype;
$me->{itype} = \@itype; $me->{itype} = \@itype;
$me->{ounit} = \@ounit; $me->{ounit} = \@ounit;
$me->{iunit} = \@iunit; $me->{iunit} = \@iunit;
# Check for nonlinear projection... # Check for nonlinear projection...
# if($hdr->{CTYPE1} =~ m/(\w\w\w\w)\-(\w\w\w)/) { # if($hdr->{CTYPE1} =~ m/(\w\w\w\w)\-(\w\w\w)/) {
# print "Nonlinear transformation found... ignoring nonlinear part...\n"; # print "Nonlinear transformation found... ignoring nonlinear part...\n";
# } # }
return $me; return $me;
} }
+======EOD_t_fits====== +======EOD_t_fits======
pp_add_exported('t_code '); pp_add_exported('t_code ');
pp_addpm(<<'+======EOD_t_code ======'); pp_addpm(<<'+======EOD_t_code ======');
=head2 t_code =head2 t_code
=for usage =for usage
skipping to change at line 3351 skipping to change at line 3270
*t_cylindrical = \&t_radial; *t_cylindrical = \&t_radial;
sub t_radial { sub t_radial {
my($class) = 'PDL::Transform'; my($class) = 'PDL::Transform';
my($o) = $_[0]; my($o) = $_[0];
if(ref $o ne 'HASH') { if(ref $o ne 'HASH') {
$o = { @_ }; $o = { @_ };
} }
my($me) = PDL::Transform::new($class); my($me) = PDL::Transform::new($class);
$me->{params}->{origin} = _opt($o,['o','origin','Origin']); $me->{params}{origin} = _opt($o,['o','origin','Origin']);
$me->{params}->{origin} = pdl(0,0) $me->{params}{origin} = pdl(0,0)
unless defined($me->{params}->{origin}); unless defined($me->{params}{origin});
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); $me->{params}{origin} = PDL->pdl($me->{params}{origin});
$me->{params}->{r0} = _opt($o,['r0','R0','c','conformal','Conformal']); $me->{params}{r0} = _opt($o,['r0','R0','c','conformal','Conformal']);
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); $me->{params}{origin} = PDL->pdl($me->{params}{origin});
$me->{params}->{u} = _opt($o,['u','unit','Unit'],'radians'); $me->{params}{u} = _opt($o,['u','unit','Unit'],'radians');
### Replace this kludge with a units call ### Replace this kludge with a units call
$me->{params}->{angunit} = ($me->{params}->{u} =~ m/^d/i) ? $PDL::Transform::R $me->{params}{angunit} = ($me->{params}{u} =~ m/^d/i) ? $PDL::Transform::RAD2D
AD2DEG : 1.0; EG : 1.0;
print "radial: conversion is $me->{params}->{angunit}\n" if($PDL::Transform::d print "radial: conversion is $me->{params}{angunit}\n" if($PDL::Transform::deb
ebug); ug);
$me->{name} = "radial (direct)"; $me->{name} = "radial (direct)";
$me->{idim} = 2; $me->{idim} = 2;
$me->{odim} = 2; $me->{odim} = 2;
if($me->{params}->{r0}) { if($me->{params}{r0}) {
$me->{otype} = ["Azimuth", "Ln radius" . ($me->{params}->{r0} != 1.0 ? "/$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 $me->{ounit} = [$me->{params}->{u},'']; # true-but-null prevents copying
} else { } else {
$me->{otype} = ["Azimuth","Radius"]; $me->{otype} = ["Azimuth","Radius"];
$me->{ounit} = [$me->{params}->{u},'']; # false value copies prev. unit $me->{ounit} = [$me->{params}->{u},'']; # false value copies prev. unit
} }
$me->{func} = sub { $me->{func} = sub {
my($data,$o) = @_; my($data,$o) = @_;
skipping to change at line 3886 skipping to change at line 3805
$me->{odim}=3; $me->{odim}=3;
$me->{params}->{origin} = _opt($o,['o','origin','Origin']); $me->{params}->{origin} = _opt($o,['o','origin','Origin']);
$me->{params}->{origin} = PDL->zeroes(3) $me->{params}->{origin} = PDL->zeroes(3)
unless defined($me->{params}->{origin}); unless defined($me->{params}->{origin});
$me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); $me->{params}->{origin} = PDL->pdl($me->{params}->{origin});
$me->{params}->{deg} = _opt($o,['d','degrees','Degrees']); $me->{params}->{deg} = _opt($o,['d','degrees','Degrees']);
my $unit = _opt($o,['u','unit','Unit']); my $unit = _opt($o,['u','unit','Unit']);
$me->{params}->{angunit} = ($unit =~ m/^d/i) ? $me->{params}{angunit} = ($unit && $unit =~ m/^d/i) ?
$PDL::Transform::DEG2RAD : undef; $PDL::Transform::DEG2RAD : undef;
$me->{name} = "spherical"; $me->{name} = "spherical";
$me->{func} = sub { $me->{func} = sub {
my($data,$o) = @_; my($data,$o) = @_;
my($d) = $data->copy - $o->{origin}; my($d) = $data->copy - $o->{origin};
my($d0,$d1,$d2) = ($d->slice("(0)"),$d->slice("(1)"),$d->slice("(2)")); my($d0,$d1,$d2) = ($d->slice("(0)"),$d->slice("(1)"),$d->slice("(2)"));
my($out) = ($d->is_inplace) ? $data : $data->copy; my($out) = ($d->is_inplace) ? $data : $data->copy;
my $tmp; # work around perl -d "feature" my $tmp; # work around perl -d "feature"
($tmp = $out->slice("(0)")) .= atan2($d1, $d0); ($tmp = $out->slice("(0)")) .= PDL::atan2($d1, $d0);
($tmp = $out->slice("(2)")) .= sqrt($d0*$d0 + $d1*$d1 + $d2*$d2); ($tmp = $out->slice("(2)")) .= PDL::sqrt($d0*$d0 + $d1*$d1 + $d2*$d2);
($tmp = $out->slice("(1)")) .= asin($d2 / $out->slice("(2)")); ($tmp = $out->slice("(1)")) .= PDL::asin($d2 / $out->slice("(2)"));
($tmp = $out->slice("0:1")) /= $o->{angunit} ($tmp = $out->slice("0:1")) /= $o->{angunit}
if(defined $o->{angunit}); if(defined $o->{angunit});
$out; $out;
}; };
$me->{inv} = sub { $me->{inv} = sub {
my($d,$o) = @_; my($d,$o) = @_;
 End of changes. 69 change blocks. 
170 lines changed or deleted 89 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)