1 2# 3# GENERATED WITH PDL::PP! Don't modify! 4# 5package PDL::Transform; 6 7@EXPORT_OK = qw( apply invert map PDL::PP map unmap t_inverse t_compose t_wrap t_identity t_lookup t_linear t_scale t_offset t_rot t_fits t_code t_cylindrical t_radial t_quadratic t_cubic t_quadratic t_spherical t_projective ); 8%EXPORT_TAGS = (Func=>[@EXPORT_OK]); 9 10use PDL::Core; 11use PDL::Exporter; 12use DynaLoader; 13 14 15 16 17 @ISA = ( 'PDL::Exporter','DynaLoader' ); 18 push @PDL::Core::PP, __PACKAGE__; 19 bootstrap PDL::Transform ; 20 21 22 23 24 25=head1 NAME 26 27PDL::Transform - Coordinate transforms, image warping, and N-D functions 28 29=head1 SYNOPSIS 30 31use PDL::Transform; 32 33 my $t = new PDL::Transform::<type>(<opt>) 34 35 $out = $t->apply($in) # Apply transform to some N-vectors (Transform method) 36 $out = $in->apply($t) # Apply transform to some N-vectors (PDL method) 37 38 $im1 = $t->map($im); # Transform image coordinates (Transform method) 39 $im1 = $im->map($t); # Transform image coordinates (PDL method) 40 41 $t2 = $t->compose($t1); # compose two transforms 42 $t2 = $t x $t1; # compose two transforms (by analogy to matrix mult.) 43 44 $t3 = $t2->inverse(); # invert a transform 45 $t3 = !$t2; # invert a transform (by analogy to logical "not") 46 47=head1 DESCRIPTION 48 49PDL::Transform is a convenient way to represent coordinate 50transformations and resample images. It embodies functions mapping 51R^N -> R^M, both with and without inverses. Provision exists for 52parametrizing functions, and for composing them. You can use this 53part of the Transform object to keep track of arbitrary functions 54mapping R^N -> R^M with or without inverses. 55 56The simplest way to use a Transform object is to transform vector 57data between coordinate systems. The L<apply|/apply> method 58accepts a PDL whose 0th dimension is coordinate index (all other 59dimensions are threaded over) and transforms the vectors into the new 60coordinate system. 61 62Transform also includes image resampling, via the L<map|/map> method. 63You define a coordinate transform using a Transform object, then use 64it to remap an image PDL. The output is a remapped, resampled image. 65 66You can define and compose several transformations, then apply them 67all at once to an image. The image is interpolated only once, when 68all the composed transformations are applied. 69 70In keeping with standard practice, but somewhat counterintuitively, 71the L<map|/map> engine uses the inverse transform to map coordinates 72FROM the destination dataspace (or image plane) TO the source dataspace; 73hence PDL::Transform keeps track of both the forward and inverse transform. 74 75For terseness and convenience, most of the constructors are exported 76into the current package with the name C<< t_<transform> >>, so the following 77(for example) are synonyms: 78 79 $t = new PDL::Transform::Radial(); # Long way 80 81 $t = t_radial(); # Short way 82 83Several math operators are overloaded, so that you can compose and 84invert functions with expression syntax instead of method syntax (see below). 85 86=head1 EXAMPLE 87 88Coordinate transformations and mappings are a little counterintuitive 89at first. Here are some examples of transforms in action: 90 91 use PDL::Transform; 92 $a = rfits('m51.fits'); # Substitute path if necessary! 93 $ts = t_linear(Scale=>3); # Scaling transform 94 95 $w = pgwin(xs); 96 $w->imag($a); 97 98 ## Grow m51 by a factor of 3; origin is at lower left. 99 $b = $ts->map($a,{pix=>1}); # pix option uses direct pixel coord system 100 $w->imag($b); 101 102 ## Shrink m51 by a factor of 3; origin still at lower left. 103 $c = $ts->unmap($a, {pix=>1}); 104 $w->imag($c); 105 106 ## Grow m51 by a factor of 3; origin is at scientific origin. 107 $d = $ts->map($a,$a->hdr); # FITS hdr template prevents autoscaling 108 $w->imag($d); 109 110 ## Shrink m51 by a factor of 3; origin is still at sci. origin. 111 $e = $ts->unmap($a,$a->hdr); 112 $w->imag($e); 113 114 ## A no-op: shrink m51 by a factor of 3, then autoscale back to size 115 $f = $ts->map($a); # No template causes autoscaling of output 116 117=head1 OPERATOR OVERLOADS 118 119=over 3 120 121=item '!' 122 123The bang is a unary inversion operator. It binds exactly as 124tightly as the normal bang operator. 125 126=item 'x' 127 128By analogy to matrix multiplication, 'x' is the compose operator, so these 129two expressions are equivalent: 130 131 $f->inverse()->compose($g)->compose($f) # long way 132 !$f x $g x $f # short way 133 134Both of those expressions are equivalent to the mathematical expression 135f^-1 o g o f, or f^-1(g(f(x))). 136 137=item '**' 138 139By analogy to numeric powers, you can apply an operator a positive 140integer number of times with the ** operator: 141 142 $f->compose($f)->compose($f) # long way 143 $f**3 # short way 144 145=back 146 147=head1 INTERNALS 148 149Transforms are perl hashes. Here's a list of the meaning of each key: 150 151=over 3 152 153=item func 154 155Ref to a subroutine that evaluates the transformed coordinates. It's 156called with the input coordinate, and the "params" hash. This 157springboarding is done via explicit ref rather than by subclassing, 158for convenience both in coding new transforms (just add the 159appropriate sub to the module) and in adding custom transforms at 160run-time. Note that, if possible, new C<func>s should support 161L<inplace|PDL::Core/inplace> operation to save memory when the data are flagged 162inplace. But C<func> should always return its result even when 163flagged to compute in-place. 164 165C<func> should treat the 0th dimension of its input as a dimensional 166index (running 0..N-1 for R^N operation) and thread over all other input 167dimensions. 168 169=item inv 170 171Ref to an inverse method that reverses the transformation. It must 172accept the same "params" hash that the forward method accepts. This 173key can be left undefined in cases where there is no inverse. 174 175=item idim, odim 176 177Number of useful dimensions for indexing on the input and output sides 178(ie the order of the 0th dimension of the coordinates to be fed in or 179that come out). If this is set to 0, then as many are allocated as needed. 180 181=item name 182 183A shorthand name for the transformation (convenient for debugging). 184You should plan on using UNIVERAL::isa to identify classes of 185transformation, e.g. all linear transformations should be subclasses 186of PDL::Transform::Linear. That makes it easier to add smarts to, 187e.g., the compose() method. 188 189=item itype 190 191An array containing the name of the quantity that is expected from the 192input piddle for the transform, for each dimension. This field is advisory, 193and can be left blank if there's no obvious quantity associated with 194the transform. This is analogous to the CTYPEn field used in FITS headers. 195 196=item oname 197 198Same as itype, but reporting what quantity is delivered for each 199dimension. 200 201=item iunit 202 203The units expected on input, if a specific unit (e.g. degrees) is expected. 204This field is advisory, and can be left blank if there's no obvious 205unit associated with the transform. 206 207=item ounit 208 209Same as iunit, but reporting what quantity is delivered for each dimension. 210 211=item params 212 213Hash ref containing relevant parameters or anything else the func needs to 214work right. 215 216=item is_inverse 217 218Bit indicating whether the transform has been inverted. That is useful 219for some stringifications (see the PDL::Transform::Linear 220stringifier), and may be useful for other things. 221 222=back 223 224Transforms should be inplace-aware where possible, to prevent excessive 225memory usage. 226 227If you define a new type of transform, consider generating a new stringify 228method for it. Just define the sub "stringify" in the subclass package. 229It should call SUPER::stringify to generate the first line (though 230the PDL::Transform::Composition bends this rule by tweaking the 231top-level line), then output (indented) additional lines as necessary to 232fully describe the transformation. 233 234=head1 NOTES 235 236Transforms have a mechanism for labeling the units and type of each 237coordinate, but it is just advisory. A routine to identify and, if 238necessary, modify units by scaling would be a good idea. Currently, 239it just assumes that the coordinates are correct for (e.g.) FITS 240scientific-to-pixel transformations. 241 242Composition works OK but should probably be done in a more 243sophisticated way so that, for example, linear transformations are 244combined at the matrix level instead of just strung together 245pixel-to-pixel. 246 247=head1 MODULE INTERFACE 248 249There are both operators and constructors. The constructors are all 250exported, all begin with "t_", and all return objects that are subclasses 251of PDL::Transform. 252 253The L<apply|/apply>, L<invert|/invert>, L<map|/map>, 254and L<unmap|/unmap> methods are also exported to the C<PDL> package: they 255are both Transform methods and PDL methods. 256 257=cut 258 259 260 261 262 263 264 265=head1 FUNCTIONS 266 267 268 269=cut 270 271 272 273 274 275=head2 apply 276 277=for sig 278 279 Signature: (data(); PDL::Transform t) 280 281=for usage 282 283 $out = $data->apply($t); 284 $out = $t->apply($data); 285 286=for ref 287 288Apply a transformation to some input coordinates. 289 290In the example, C<$t> is a PDL::Transform and C<$data> is a PDL to 291be interpreted as a collection of N-vectors (with index in the 0th 292dimension). The output is a similar but transformed PDL. 293 294For convenience, this is both a PDL method and a Transform method. 295 296=cut 297 298use Carp; 299 300*PDL::apply = \&apply; 301sub apply { 302 my($me) = shift; 303 my($from) = shift; 304 305 if(UNIVERSAL::isa($me,'PDL')){ 306 my($a) = $from; 307 $from = $me; 308 $me = $a; 309 } 310 311 if(UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($from,'PDL')){ 312 croak "Applying a PDL::Transform with no func! Oops.\n" unless(defined($me->{func}) and ref($me->{func}) eq 'CODE'); 313 my $result = &{$me->{func}}($from,$me->{params}); 314 $result->is_inplace(0); # clear inplace flag, just in case. 315 return $result; 316 } else { 317 croak "apply requires both a PDL and a PDL::Transform.\n"; 318 } 319 320} 321 322 323 324 325=head2 invert 326 327=for sig 328 329 Signature: (data(); PDL::Transform t) 330 331=for usage 332 333 $out = $t->invert($data); 334 $out = $data->invert($t); 335 336=for ref 337 338Apply an inverse transformation to some input coordinates. 339 340In the example, C<$t> is a PDL::Transform and C<$data> is a piddle to 341be interpreted as a collection of N-vectors (with index in the 0th 342dimension). The output is a similar piddle. 343 344For convenience this is both a PDL method and a PDL::Transform method. 345 346=cut 347 348*PDL::invert = \&invert; 349sub invert { 350 my($me) = shift; 351 my($data) = shift; 352 353 if(UNIVERSAL::isa($me,'PDL')){ 354 my($a) = $data; 355 $data = $me; 356 $me = $a; 357 } 358 359 if(UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($data,'PDL')){ 360 croak "Inverting a PDL::Transform with no inverse! Oops.\n" unless(defined($me->{inv}) and ref($me->{inv}) eq 'CODE'); 361 my $result = &{$me->{inv}}($data, $me->{params}); 362 $result->is_inplace(0); # make sure inplace flag is clear. 363 return $result; 364 } else { 365 croak("invert requires a PDL and a PDL::Transform (did you want 'inverse' instead?)\n"); 366 } 367} 368 369 370 371 372 373=head2 map 374 375=for sig 376 377 Signature: (k0(); SV *in; SV *out; SV *map; SV *boundary; SV *method; 378 SV *big; SV *blur; SV *sv_min; SV *flux; SV *bv) 379 380 381=head2 match 382 383=for usage 384 385 $b = $a->match($c); # Match $c's header and size 386 $b = $a->match([100,200]); # Rescale to 100x200 pixels 387 $b = $a->match([100,200],{rect=>1}); # Rescale and remove rotation/skew. 388 389=for ref 390 391Resample a scientific image to the same coordinate system as another. 392 393The example above is syntactic sugar for 394 395 $b = $a->map(t_identity, $c, ...); 396 397it resamples the input PDL with the identity transformation in 398scientific coordinates, and matches the pixel coordinate system to 399$c's FITS header. 400 401There is one difference between match and map: match makes the 402C<rectify> option to C<map> default to 0, not 1. This only affects 403matching where autoscaling is required (i.e. the array ref example 404above). By default, that example simply scales $a to the new size and 405maintains any rotation or skew in its scientiic-to-pixel coordinate 406transform. 407 408=head2 map 409 410=for usage 411 412 $b = $a->map($xform,[<template>],[\%opt]); # Distort $a with transform $xform 413 $b = $a->map(t_identity,[$pdl],[\%opt]); # rescale $a to match $pdl's dims. 414 415=for ref 416 417Resample an image or N-D dataset using a coordinate transform. 418 419The data are resampled so that the new pixel indices are proportional 420to the transformed coordinates rather than the original ones. 421 422The operation uses the inverse transform: each output pixel location 423is inverse-transformed back to a location in the original dataset, and 424the value is interpolated or sampled appropriately and copied into the 425output domain. A variety of sampling options are available, trading 426off speed and mathematical correctness. 427 428For convenience, this is both a PDL method and a PDL::Transform method. 429 430C<map> is FITS-aware: if there is a FITS header in the input data, 431then the coordinate transform acts on the scientific coordinate system 432rather than the pixel coordinate system. 433 434By default, the output coordinates are separated from pixel coordinates 435by a single layer of indirection. You can specify the mapping between 436output transform (scientific) coordinates to pixel coordinates using 437the C<orange> and C<irange> options (see below), or by supplying a 438FITS header in the template. 439 440If you don't specify an output transform, then the output is 441autoscaled: C<map> transforms a few vectors in the forward direction 442to generate a mapping that will put most of the data on the image 443plane, for most transformations. The calculated mapping gets stuck in the 444output's FITS header. 445 446Autoscaling is especially useful for rescaling images -- if you specify 447the identity transform and allow autoscaling, you duplicate the 448functionality of L<rescale2d|PDL::Image2D/rescale2d>, but with more 449options for interpolation. 450 451You can operate in pixel space, and avoid autoscaling of the output, 452by setting the C<nofits> option (see below). 453 454The output has the same data type as the input. This is a feature, 455but it can lead to strange-looking banding behaviors if you use 456interpolation on an integer input variable. 457 458The C<template> can be one of: 459 460=over 3 461 462=item * a PDL 463 464The PDL and its header are copied to the output array, which is then 465populated with data. If the PDL has a FITS header, then the FITS 466transform is automatically applied so that $t applies to the output 467scientific coordinates and not to the output pixel coordinates. In 468this case the NAXIS fields of the FITS header are ignored. 469 470=item * a FITS header stored as a hash ref 471 472The FITS NAXIS fields are used to define the output array, and the 473FITS transformation is applied to the coordinates so that $t applies to the 474output scientific coordinates. 475 476=item * a list ref 477 478This is a list of dimensions for the output array. The code estimates 479appropriate pixel scaling factors to fill the available space. The 480scaling factors are placed in the output FITS header. 481 482=item * nothing 483 484In this case, the input image size is used as a template, and scaling 485is done as with the list ref case (above). 486 487=back 488 489OPTIONS: 490 491The following options are interpreted: 492 493=over 3 494 495=item b, bound, boundary, Boundary (default = 'truncate') 496 497This is the boundary condition to be applied to the input image; it is 498passed verbatim to L<range|PDL::Slices/range> or 499L<interpND|PDL::Primitive/interpND> in the sampling or interpolating 500stage. Other values are 'forbid','extend', and 'periodic'. You can 501abbreviate this to a single letter. The default 'truncate' causes the 502entire notional space outside the original image to be filled with 0. 503 504=item pix, Pixel, nf, nofits, NoFITS (default = 0) 505 506If you set this to a true value, then FITS headers and interpretation 507are ignored; the transformation is treated as being in raw pixel coordinates. 508 509=item j, J, just, justify, Justify (default = 0) 510 511If you set this to 1, then output pixels are autoscaled to have unit 512aspect ratio in the output coordinates. If you set it to a non-1 513value, then it is the aspect ratio between the first dimension and all 514subsequent dimensions -- or, for a 2-D transformation, the scientific 515pixel aspect ratio. Values less than 1 shrink the scale in the first 516dimension compared to the other dimensions; values greater than 1 517enlarge it compared to the other dimensions. (This is the same sense 518as in the L<PGPLOT|PDL::Graphics::PGPLOT>interface.) 519 520=item ir, irange, input_range, Input_Range 521 522This is a way to modify the autoscaling. It specifies the range of 523input scientific (not necessarily pixel) coordinates that you want to be 524mapped to the output image. It can be either a nested array ref or 525a piddle. The 0th dim (outside coordinate in the array ref) is 526dimension index in the data; the 1st dim should have order 2. 527For example, passing in either [[-1,2],[3,4]] or pdl([[-1,2],[3,4]]) 528limites the map to the quadrilateral in input space defined by the 529four points (-1,3), (-1,4), (2,4), and (2,3). 530 531As with plain autoscaling, the quadrilateral gets sparsely sampled by 532the autoranger, so pathological transformations can give you strange 533results. 534 535This parameter is overridden by C<orange>, below. 536 537=item or, orange, output_range, Output_Range 538 539This sets the window of output space that is to be sampled onto the 540output array. It works exactly like C<irange>, except that it specifies 541a quadrilateral in output space. Since the output pixel array is itself 542a quadrilateral, you get pretty much exactly what you asked for. 543 544This parameter overrides C<irange>, if both are specified. It forces 545rectification of the output (so that scientific axes align with the pixel 546grid). 547 548=item r, rect, rectify 549 550This option defaults TRUE and controls how autoscaling is performed. If 551it is true or undefined, then autoscaling adjusts so that pixel coordinates 552in the output image are proportional to individual scientific coordinates. 553If it is false, then autoscaling automatically applies the inverse of any 554input FITS transformation *before* autoscaling the pixels. In the special 555case of linear transformations, this preserves the rectangular shape of the 556original pixel grid and makes output pixel coordinate proportional to input 557coordinate. 558 559=item m, method, Method 560 561This option controls the interpolation method to be used. 562Interpolation greatly affects both speed and quality of output. For 563most cases the option is directly passed to 564L<interpND|PDL::Primitive/interpnd> for interpolation. Possible 565options, in order from fastest to slowest, are: 566 567=over 3 568 569 570=item * s, sample (default for ints) 571 572Pixel values in the output plane are sampled from the closest data value 573in the input plane. This is very fast but not very accurate for either 574magnification or decimation (shrinking). It is the default for templates 575of integer type. 576 577=item * l, linear (default for floats) 578 579Pixel values are linearly interpolated from the closest data value in the 580input plane. This is reasonably fast but only accurate for magnification. 581Decimation (shrinking) of the image causes aliasing and loss of photometry 582as features fall between the samples. It is the default for floating-point 583templates. 584 585=item * c, cubic 586 587Pixel values are interpolated using an N-cubic scheme from a 4-pixel 588N-cube around each coordinate value. As with linear interpolation, 589this is only accurate for magnification. 590 591=item * f, fft 592 593Pixel values are interpolated using the term coefficients of the 594Fourier transform of the original data. This is the most appropriate 595technique for some kinds of data, but can yield undesired "ringing" for 596expansion of normal images. Best suited to studying images with 597repetitive or wavelike features. 598 599=item * h, hanning 600 601Pixel values are filtered through a spatially-variable filter tuned to 602the computed Jacobian of the transformation, with hanning-window 603(cosine) pixel rolloff in each dimension. This prevents aliasing in the 604case where the image is distorted or shrunk, but allows small amounts 605of aliasing at pixel edges wherever the image is enlarged. 606 607=item * g, gaussian, j, jacobian 608 609Pixel values are filtered through a spatially-variable filter tuned to 610the computed Jacobian of the transformation, with radial Gaussian 611rolloff. This is the most accurate resampling method, in the sense of 612introducing the fewest artifacts into a properly sampled data set. 613This method uses a lookup table to speed up calculation of the Gaussian 614weighting. 615 616=item * G 617 618This method works exactly like 'g' (above), except that the Gaussian 619values are computed explicitly for every sample -- which is considerably 620slower. 621 622=back 623 624=item blur, Blur (default = 1.0) 625 626This value scales the input-space footprint of each output pixel in 627the gaussian and hanning methods. It's retained for historical 628reasons. Larger values yield blurrier images; values significantly 629smaller than unity cause aliasing. The parameter has slightly 630different meanings for Gaussian and Hanning interpolation. For 631Hanning interpolation, numbers smaller than unity control the 632sharpness of the border at the edge of each pixel (so that blur=>0 is 633equivalent to sampling for non-decimating transforms). For 634Gaussian interpolation, the blur factor parameter controls the 635width parameter of the Gaussian around the center of each pixel. 636 637=item sv, SV (default = 1.0) 638 639This value lets you set the lower limit of the transformation's 640singular values in the hanning and gaussian methods, limiting the 641minimum radius of influence associated with each output pixel. Large 642numbers yield smoother interpolation in magnified parts of the image 643but don't affect reduced parts of the image. 644 645=item big, Big (default = 0.5) 646 647This is the largest allowable input spot size which may be mapped to a 648single output pixel by the hanning and gaussian methods, in units of 649the largest non-thread input dimension. (i.e. the default won't let 650you reduce the original image to less than 5 pixels across). This places 651a limit on how long the processing can take for pathological transformations. 652Smaller numbers keep the code from hanging for a long time; larger numbers 653provide for photometric accuracy in more pathological cases. Numbers larer 654than 1.0 are silly, because they allow the entire input array to be compressed 655into a region smaller than a single pixel. 656 657Wherever an output pixel would require averaging over an area that is too 658big in input space, it instead gets NaN or the bad value. 659 660=item phot, photometry, Photometry 661 662This lets you set the style of photometric conversion to be used in the 663hanning or gaussian methods. You may choose: 664 665=over 3 666 667=item * 0, s, surf, surface, Surface (default) 668 669(this is the default): surface brightness is preserved over the transformation, 670so features maintain their original intensity. This is what the sampling 671and interpolation methods do. 672 673=item * 1, f, flux, Flux 674 675Total flux is preserved over the transformation, so that the brightness 676integral over image regions is preserved. Parts of the image that are 677shrunk wind up brighter; parts that are enlarged end up fainter. 678 679=back 680 681=back 682 683VARIABLE FILTERING: 684 685The 'hanning' and 'gaussian' methods of interpolation give 686photometrically accurate resampling of the input data for arbitrary 687transformations. At each pixel, the code generates a linear 688approximation to the input transformation, and uses that linearization 689to estimate the "footprint" of the output pixel in the input space. 690The output value is a weighted average of the appropriate input spaces. 691 692A caveat about these methods is that they assume the transformation is 693continuous. Transformations that contain discontinuities will give 694incorrect results near the discontinuity. In particular, the 180th 695meridian isn't handled well in lat/lon mapping transformations (see 696L<PDL::Transform::Cartography>) -- pixels along the 180th meridian get 697the average value of everything along the parallel occupied by the 698pixel. This flaw is inherent in the assumptions that underly creating 699a Jacobian matrix. Maybe someone will write code to work around it. 700Maybe that someone is you. 701 702BAD VALUES: 703 704If your PDL was compiled with bad value support, C<map()> supports 705bad values in the data array. Bad values in the input array are 706propagated to the output array. The 'g' and 'h' methods will do some 707smoothing over bad values: if more than 1/3 of the weighted input-array 708footprint of a given output pixel is bad, then the output pixel gets marked 709bad. 710 711 712 713=for bad 714 715map does not process bad values. 716It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 717 718 719=cut 720 721 722 723 724 725sub PDL::match { 726 # Set default for rectification to 0 for simple matching... 727 if( ref($_[$#_]) ne 'HASH' ) { 728 push(@_,{}) 729 } 730 my @k = grep(m/^r(e(c(t)?)?)?/,keys %{$_[$#_]}); 731 unless(@k) { 732 $_[$#_]->{rectify} = 0; 733 } 734 t_identity()->map(@_); 735} 736 737 738*PDL::map = \↦ 739sub map { 740 my($me) = shift; 741 my($in) = shift; 742 743 if(UNIVERSAL::isa($me,'PDL') && UNIVERSAL::isa($in,'PDL::Transform')) { 744 my($a) = $in; 745 $in = $me; 746 $me = $a; 747 } 748 749 barf ("PDL::Transform::map: source is not defined or is not a PDL\n") 750 unless(defined $in and UNIVERSAL::isa($in,'PDL')); 751 752 barf ("PDL::Transform::map: source is empty\n") 753 unless($in->nelem); 754 755 my($tmp) = shift; 756 my($opt) = shift; 757 758 # Check for options-but-no-template case 759 if(ref $tmp eq 'HASH' && !(defined $opt)) { 760 if(!defined($tmp->{NAXIS})) { # FITS headers all have NAXIS. 761 $opt = $tmp; 762 $tmp = undef; 763 } 764 } 765 766 croak("PDL::Transform::map: Option 'p' was ambiguous and has been removed. You probably want 'pix' or 'phot'.") if exists($opt->{'p'}); 767 768 $tmp = [$in->dims] unless(defined($tmp)); 769 770 # Generate an appropriate output piddle for values to go in 771 my($out); 772 my(@odims); 773 my($ohdr); 774 if(UNIVERSAL::isa($tmp,'PDL')) { 775 @odims = $tmp->dims; 776 777 my($a); 778 if(defined ($a = $tmp->gethdr)) { 779 my(%b) = %{$a}; 780 $ohdr = \%b; 781 } 782 } elsif(ref $tmp eq 'HASH') { 783 # (must be a fits header -- or would be filtered above) 784 for my $i(1..$tmp->{NAXIS}){ 785 push(@odims,$tmp->{"NAXIS$i"}); 786 } 787 # deep-copy fits header into output 788 my %foo = %{$tmp}; 789 $ohdr = \%foo; 790 } elsif(ref $tmp eq 'ARRAY') { 791 @odims = @$tmp; 792 } else { 793 barf("map: confused about dimensions of the output array...\n"); 794 } 795 796 if(scalar(@odims) < scalar($in->dims)) { 797 my @idims = $in->dims; 798 push(@odims, splice(@idims,scalar(@odims))); 799 } 800 801 $out = PDL::new_from_specification('PDL',$in->type,@odims); 802 $out->sethdr($ohdr) if defined($ohdr); 803 804 if($PDL::Bad::Status) { 805 # set badflag on output all the time if possible, to account for boundary violations 806 $out->badflag(1); 807 } 808 809 ############################## 810 ## Figure out the dimensionality of the 811 ## transform itself (extra dimensions come along for the ride) 812 my $nd = $me->{odim} || $me->{idim} || 2; 813 my @sizes = $out->dims; 814 my @dd = @sizes; 815 816 splice @dd,$nd; # Cut out dimensions after the end 817 818 # Check that there are elements in the output fields... 819 barf "map: output has no dims!\n" 820 unless(@dd); 821 my $ddtotal = 1; 822 map {$ddtotal *= $_} @dd; 823 barf "map: output has no elements (at least one dim is 0)!\n" 824 unless($ddtotal); 825 826 827 ############################## 828 # If necessary, generate an appropriate FITS header for the output. 829 830 my $nofits = _opt($opt, ['nf','nofits','NoFITS','pix','pixel','Pixel']); 831 832 ############################## 833 # Autoscale by transforming a subset of the input points' coordinates 834 # to the output range, and pick a FITS header that fits the output 835 # coordinates into the given template. 836 # 837 # Autoscaling always produces a simple, linear mapping in the FITS header. 838 # We support more complex mappings (via t_fits) but only to match a pre-existing 839 # FITS header (which doesn't use autoscaling). 840 # 841 # If the rectify option is set (the default) then the image is rectified 842 # in scientific coordinates; if it is clear, then the existing matrix 843 # is used, preserving any shear or rotation in the coordinate system. 844 # Since we eschew CROTA whenever possible, the CDi_j formalism is used instead. 845 my $f_in = (defined($in->hdr->{NAXIS}) ? t_fits($in,{ignore_rgb=>1}) : t_identity()); 846 847 unless((defined $out->gethdr && $out->hdr->{NAXIS}) or $nofits) { 848 print "generating output FITS header..." if($PDL::Transform::debug); 849 850 $out->sethdr($in->hdr_copy) # Copy extraneous fields... 851 if(defined $in->hdr); 852 853 my $samp_ratio = 300; 854 855 my $orange = _opt($opt, ['or','orange','output_range','Output_Range'], 856 undef); 857 858 my $omin; 859 my $omax; 860 my $osize; 861 862 863 my $rectify = _opt($opt,['r','rect','rectify','Rectify'],1); 864 865 866 if (defined $orange) { 867 # orange always rectifies the coordinates -- the output scientific 868 # coordinates *must* align with the axes, or orange wouldn't make 869 # sense. 870 print "using user's orange..." if($PDL::Transform::debug); 871 $orange = pdl($orange) unless(UNIVERSAL::isa($orange,'PDL')); 872 barf "map: orange must be 2xN for an N-D transform" 873 unless ( (($orange->dim(1)) == $nd ) 874 && $orange->ndims == 2); 875 876 $omin = $orange->slice("(0)"); 877 $omax = $orange->slice("(1)"); 878 $osize = $omax - $omin; 879 880 $rectify = 1; 881 882 } else { 883 884 ############################## 885 # Real autoscaling happens here. 886 887 if(!$rectify and ref( $f_in ) !~ /Linear/i) { 888 if( $f_in->{name} ne 'identity' ) { 889 print STDERR "Warning: map can't preserve nonlinear FITS distortions while autoscaling.\n"; 890 } 891 $rectify=1; 892 } 893 894 my $f_tr = ( $rectify ? 895 $me x $f_in : 896 ( ($me->{name} eq 'identity') ? # Simple optimization for match() 897 $me : # identity -- just matching 898 !$f_in x $me x $f_in # common case 899 ) 900 ); 901 902 my $samps = (pdl(($in->dims)[0..$nd-1]))->clip(0,$samp_ratio); 903 904 my $coords = ndcoords(($samps + 1)->list); 905 906 my $t; 907 my $irange = _opt($opt, ['ir','irange','input_range','Input_Range'], 908 undef); 909 910 # If input range is defined, sample that quadrilateral -- else 911 # sample the quad defined by the boundaries of the input image. 912 if(defined $irange) { 913 print "using user's irange..." if($PDL::Transform::debug); 914 $irange = pdl($irange) unless(UNIVERSAL::isa($irange,'PDL')); 915 barf "map: irange must be 2xN for an N-D transform" 916 unless ( (($irange->dim(1)) == $nd ) 917 && $irange->ndims == 2); 918 919 $coords *= ($irange->slice("(1)") - $irange->slice("(0)")) / $samps; 920 $coords += $irange->slice("(0)"); 921 $coords -= 0.5; # offset to pixel corners... 922 $t = $me; 923 } else { 924 $coords *= pdl(($in->dims)[0..$nd-1]) / $samps; 925 $coords -= 0.5; # offset to pixel corners... 926 $t = $f_tr; 927 } 928 my $ocoords = $t->apply($coords)->mv(0,-1)->clump($nd); 929 930 # discard non-finite entries 931 my $oc2 = $ocoords->range( 932 which( 933 $ocoords-> 934 xchg(0,1)-> 935 sumover-> 936 isfinite 937 ) 938 ->dummy(0,1) 939 ); 940 941 $omin = $oc2->minimum; 942 $omax = $oc2->maximum; 943 944 $osize = $omax - $omin; 945 my $tosize; 946 ($tosize = $osize->where($osize == 0)) .= 1.0; 947 } 948 949 my ($scale) = $osize / pdl(($out->dims)[0..$nd-1]); 950 951 my $justify = _opt($opt,['j','J','just','justify','Justify'],0); 952 if($justify) { 953 my $tmp; # work around perl -d "feature" 954 ($tmp = $scale->slice("0")) *= $justify; 955 $scale .= $scale->max; 956 $scale->slice("0") /= $justify; 957 } 958 959 print "done with autoscale. Making fits header....\n" if($PDL::Transform::debug); 960 if( $rectify ) { 961 # Rectified header generation -- make a simple coordinate header with no 962 # rotation or skew. 963 print "rectify\n" if($PDL::Transform::debug); 964 for my $d(1..$nd) { 965 $out->hdr->{"CRPIX$d"} = 1 + ($out->dim($d-1)-1)/2 ; 966 $out->hdr->{"CDELT$d"} = $scale->at($d-1); 967 $out->hdr->{"CRVAL$d"} = ( $omin->at($d-1) + $omax->at($d-1) ) /2 ; 968 $out->hdr->{"NAXIS$d"} = $out->dim($d-1); 969 $out->hdr->{"CTYPE$d"} = ( (defined($me->{otype}) ? 970 $me->{otype}->[$d-1] : "") 971 || $in->hdr->{"CTYPE$d"} 972 || ""); 973 $out->hdr->{"CUNIT$d"} = ( (defined($me->{ounit}) ? 974 $me->{ounit}->[$d-1] : "") 975 || $in->hdr->{"CUNIT$d"} 976 || $in->hdr->{"CTYPE$d"} 977 || ""); 978 } 979 $out->hdr->{"NAXIS"} = $nd; 980 981 $out->hdr->{"SIMPLE"} = 'T'; 982 $out->hdr->{"HISTORY"} .= "Header written by PDL::Transform::Cartography::map"; 983 984 ### Eliminate fancy newfangled output header pointing tags if they exist 985 ### These are the CROTA<n>, PCi_j, and CDi_j. 986 for $k(keys %{$out->hdr}) { 987 if( $k=~m/(^CROTA\d*$)|(^(CD|PC)\d+_\d+[A-Z]?$)/ ){ 988 delete $out->hdr->{$k}; 989 } 990 } 991 } else { 992 # Non-rectified output -- generate a CDi_j matrix instead of the simple formalism. 993 # We have to deal with a linear transformation: we've got: (scaling) x !input x (t x input), 994 # where input is a linear transformation with offset and scaling is a simple scaling. We have 995 # the scaling parameters and the matrix for !input -- we have only to merge them and then we 996 # can write out the FITS header in CDi_j form. 997 print "non-rectify\n" if($PDL::Transform::debug); 998 my $midpoint_val = (pdl(($out->dims)[0..$nd-1])/2 * $scale)->apply( $f_in ); 999 print "midpoint_val is $midpoint_val\n" if($PDL::Transform::debug); 1000 # linear transformation 1001 unless(ref($f_in) =~ m/Linear/) { 1002 croak("Whups -- got a nonlinear t_fits transformation. Can't deal with it."); 1003 } 1004 1005 my $inv_sc_mat = zeroes($nd,$nd); 1006 $inv_sc_mat->diagonal(0,1) .= $scale; 1007 my $mat = $f_in->{params}->{matrix} x $inv_sc_mat; 1008 print "scale is $scale; mat is $mat\n" if($PDL::Transform::debug); 1009 1010 print "looping dims 1..$nd: " if($PDL::Transform::debug); 1011 for my $d(1..$nd) { 1012 print "$d..." if($PDL::Transform::debug); 1013 $out->hdr->{"CRPIX$d"} = 1 + ($out->dim($d-1)-1)/2; 1014 $out->hdr->{"CRVAL$d"} = $midpoint_val->at($d-1); 1015 $out->hdr->{"NAXIS$d"} = $out->dim($d-1); 1016 $out->hdr->{"CTYPE$d"} = ( (defined($me->{otype}) ? 1017 $me->{otype}->[$d-1] : "") 1018 || $in->hdr->{"CTYPE$d"} 1019 || ""); 1020 $out->hdr->{"CUNIT$d"} = ( (defined($me->{ounit}) ? 1021 $me->{ounit}->[$d-1] : "") 1022 || $in->hdr->{"CUNIT$d"} 1023 || $in->hdr->{"CTYPE$d"} 1024 || ""); 1025 for my $e(1..$nd) { 1026 $out->hdr->{"CD${d}_${e}"} = $mat->at($d-1,$e-1); 1027 print "setting CD${d}_${e} to ".$mat->at($d-1,$e-1)."\n" if($PDL::Transform::debug); 1028 } 1029 } 1030 1031 ## Eliminate competing header pointing tags if they exist 1032 for $k(keys %{$out->hdr}) { 1033 if( $k =~ m/(^CROTA\d*$)|(^(PC)\d+_\d+[A-Z]?$)|(CDELT\d*$)/ ) { 1034 delete $out->hdr->{$k}; 1035 } 1036 } 1037 } 1038 1039 1040 1041 } 1042 1043 $out->hdrcpy(1); 1044 1045 ############################## 1046 # Sandwich the transform between the input and output plane FITS headers. 1047 unless($nofits) { 1048 $me = !(t_fits($out,{ignore_rgb=>1})) x $me x $f_in; 1049 } 1050 1051 ############################## 1052 ## Figure out the interpND options 1053 my $method = _opt($opt,['m','method','Method'],undef); 1054 my $bound = _opt($opt,['b','bound','boundary','Boundary'],'t'); 1055 1056 1057 ############################## 1058 ## Rubber meets the road: calculate the inverse transformed points. 1059 my $ndc = PDL::Basic::ndcoords(@dd); 1060 my $idx = $me->invert($ndc->double); 1061 1062 barf "map: Transformation had no inverse\n" unless defined($idx); 1063 1064 ############################## 1065 ## Integrate ? (Jacobian, Gaussian, Hanning) 1066 my $integrate = ($method =~ m/^[jghJGH]/) if defined($method); 1067 1068 ############################## 1069 ## Sampling code: 1070 ## just transform and interpolate. 1071 ## ( Kind of an anticlimax after all that, eh? ) 1072 if(!$integrate) { 1073 my $a = $in->interpND($idx,{method=>$method, bound=>$bound}); 1074 my $tmp; # work around perl -d "feature" 1075 ($tmp = $out->slice(":")) .= $a; # trivial slice prevents header overwrite... 1076 return $out; 1077 } 1078 1079 ############################## 1080 ## Anti-aliasing code: 1081 ## Condition the input and call the pixelwise C interpolator. 1082 ## 1083 1084 barf("PDL::Transform::map: Too many dims in transformation\n") 1085 if($in->ndims < $idx->ndims-1); 1086 1087 #################### 1088 ## Condition the threading -- pixelwise interpolator only threads 1089 ## in 1 dimension, so squish all thread dimensions into 1, if necessary 1090 my @iddims = $idx->dims; 1091 if($in->ndims == $#iddims) { 1092 $in2 = $in->dummy(-1,1); 1093 } else { 1094 $in2 = ( $in 1095 ->reorder($nd..$in->ndims-1, 0..$nd-1) 1096 ->clump($in->ndims - $nd) 1097 ->mv(0,-1) 1098 ); 1099 } 1100 1101 #################### 1102 # Allocate the output array 1103 my $o2 = PDL->new_from_specification($in2->type, 1104 @iddims[1..$#iddims], 1105 $in2->dim(-1) 1106 ); 1107 1108 #################### 1109 # Pack boundary string if necessary 1110 if(defined $bound) { 1111 if(ref $bound eq 'ARRAY') { 1112 my ($s,$el); 1113 foreach $el(@$bound) { 1114 barf "Illegal boundary value '$el' in range" 1115 unless( $el =~ m/^([0123fFtTeEpPmM])/ ); 1116 $s .= $1; 1117 } 1118 $bound = $s; 1119 } 1120 elsif($bound !~ m/^[0123ftepx]+$/ && $bound =~ m/^([0123ftepx])/i ) { 1121 $bound = $1; 1122 } 1123 } 1124 1125 #################### 1126 # Get the blur and minimum-sv values 1127 my $blur = _opt($opt,['blur','Blur'],1.0); 1128 my $svmin = _opt($opt,['sv','SV'],1.0); 1129 my $big = _opt($opt,['big','Big'],1.0); 1130 my $flux = _opt($opt,['phot','photometry'],0); 1131 my @idims = $in->dims; 1132 1133 $flux = ($flux =~ m/^[1fF]/); 1134 $big = $big * max(pdl(@idims[0..$nd])); 1135 $blur = $blur->at(0) if(ref $blur); 1136 $svmin = $svmin->at(0) if(ref $svmin); 1137 1138 my $bv; 1139 if($PDL::Bad::Status and $in->badflag){ 1140 $bv = $in->badvalue; 1141 } else { 1142 $bv = 0; 1143 } 1144 1145 ### The first argument is a dummy to set $GENERIC. 1146 $idx = double($idx) unless($idx->type == double); 1147 print "Calling _map_int...\n" if($PDL::Transform::debug); 1148 &PDL::_map_int( $in2->flat->index(0), 1149 $in2, $o2, $idx, 1150 $bound, $method, $big, $blur, $svmin, $flux, $bv); 1151 1152 my @rdims = (@iddims[1..$#iddims], @idims[$#iddims..$#idims]); 1153 { 1154 my $tmp; # work around perl -d "feature" 1155 ($tmp = $out->slice(":")) .= $o2->reshape(@rdims); 1156 } 1157 return $out; 1158} 1159 1160 1161 1162*map = \&PDL::map; 1163 1164 1165 1166 1167###################################################################### 1168 1169=head2 unmap 1170 1171=for sig 1172 1173 Signature: (data(); PDL::Transform a; template(); \%opt) 1174 1175=for usage 1176 1177 $out_image = $in_image->unmap($t,[<options>],[<template>]); 1178 $out_image = $t->unmap($in_image,[<options>],[<template>]); 1179 1180=for ref 1181 1182Map an image or N-D dataset using the inverse as a coordinate transform. 1183 1184This convenience function just inverts $t and calls L<map|/map> on 1185the inverse; everything works the same otherwise. For convenience, it 1186is both a PDL method and a PDL::Transform method. 1187 1188=cut 1189 1190*PDL::unmap = \&unmap; 1191sub unmap { 1192 my($me) = shift; 1193 my($data) = shift; 1194 my(@params) = @_; 1195 1196 if(UNIVERSAL::isa($data,'PDL::Transform') && UNIVERSAL::isa($me,'PDL')) { 1197 my $a = $data; 1198 $data = $me; 1199 $me = $a; 1200 } 1201 1202 return $me->inverse->map($data,@params); 1203} 1204 1205 1206 1207 1208=head2 t_inverse 1209 1210=for usage 1211 1212 $t2 = t_inverse($t); 1213 $t2 = $t->inverse; 1214 $t2 = $t ** -1; 1215 $t2 = !$t; 1216 1217=for ref 1218 1219Return the inverse of a PDL::Transform. This just reverses the 1220func/inv, idim/odim, itype/otype, and iunit/ounit pairs. Note that 1221sometimes you end up with a transform that cannot be applied or 1222mapped, because either the mathematical inverse doesn't exist or the 1223inverse func isn't implemented. 1224 1225You can invert a transform by raising it to a negative power, or by 1226negating it with '!'. 1227 1228The inverse transform remains connected to the main transform because 1229they both point to the original parameters hash. That turns out to be 1230useful. 1231 1232=cut 1233 1234*t_inverse = \&inverse; 1235 1236sub inverse { 1237 my($me) = shift; 1238 1239 unless(defined($me->{inv})) { 1240 Carp::cluck("PDL::Transform::inverse: got a transform with no inverse.\n"); 1241 return undef; 1242 } 1243 1244 my(%out) = %$me; # force explicit copy of top-level 1245 my($out) = \%out; 1246 1247 $out->{inv} = $me->{func}; 1248 $out->{func} = $me->{inv}; 1249 1250 $out->{idim} = $me->{odim}; 1251 $out->{odim} = $me->{idim}; 1252 1253 $out->{otype} = $me->{itype}; 1254 $out->{itype} = $me->{otype}; 1255 1256 $out->{ounit} = $me->{iunit}; 1257 $out->{iunit} = $me->{ounit}; 1258 1259 $out->{name} = "(inverse ".$me->{name}.")"; 1260 1261 $out->{is_inverse} = !($out->{is_inverse}); 1262 1263 bless $out,(ref $me); 1264 return $out; 1265} 1266 1267 1268 1269 1270=head2 t_compose 1271 1272=for usage 1273 1274 $f2 = t_compose($f, $g,[...]); 1275 $f2 = $f->compose($g[,$h,$i,...]); 1276 $f2 = $f x $g x ...; 1277 1278=for ref 1279 1280Function composition: f(g(x)), f(g(h(x))), ... 1281 1282You can also compose transforms using the overloaded matrix-multiplication 1283(nee repeat) operator 'x'. 1284 1285This is accomplished by inserting a splicing code ref into the C<func> 1286and C<inv> slots. It combines multiple compositions into a single 1287list of transforms to be executed in order, fram last to first (in 1288keeping with standard mathematical notation). If one of the functions is 1289itself a composition, it is interpolated into the list rather than left 1290separate. Ultimately, linear transformations may also be combined within 1291the list. 1292 1293No checking is done that the itype/otype and iunit/ounit fields are 1294compatible -- that may happen later, or you can implement it yourself 1295if you like. 1296 1297=cut 1298 1299@PDL::Transform::Composition::ISA = ('PDL::Transform'); 1300sub PDL::Transform::Composition::stringify { 1301 package PDL::Transform::Composition; 1302 my($me) = shift; 1303 my($out) = SUPER::stringify $me; 1304 $out; 1305} 1306 1307*t_compose = \&compose; 1308 1309sub compose { 1310 local($_); 1311 my(@funcs) = @_; 1312 my($me) = PDL::Transform->new; 1313 1314 # No inputs case: return the identity function 1315 return $me 1316 if(!@funcs); 1317 1318 $me->{name} = ""; 1319 my($f); 1320 my(@clist); 1321 1322 for $f(@funcs) { 1323 1324 $me->{idim} = $f->{idim} if($f->{idim} > $me->{idim}); 1325 $me->{odim} = $f->{odim} if($f->{odim} > $me->{odim}); 1326 1327 if(UNIVERSAL::isa($f,"PDL::Transform::Composition")) { 1328 if($f->{is_inverse}) { 1329 for(reverse(@{$f->{params}->{clist}})) { 1330 push(@clist,$_->inverse); 1331 $me->{name} .= " o inverse ( ".$_->{name}." )"; 1332 } 1333 } else { 1334 for(@{$f->{params}->{clist}}) { 1335 push(@clist,$_); 1336 $me->{name} .= " o ".$_->{name}; 1337 } 1338 } 1339 } else { # Not a composition -- just push the transform onto the list. 1340 push(@clist,$f); 1341 $me->{name} .= " o ".$f->{name}; 1342 } 1343 } 1344 1345 $me->{name}=~ s/^ o //; # Get rid of leading composition mark 1346 1347 $me->{otype} = $funcs[0]->{otype}; 1348 $me->{ounit} = $funcs[0]->{ounit}; 1349 1350 $me->{itype} = $funcs[-1]->{itype}; 1351 $me->{iunit} = $funcs[-1]->{iunit}; 1352 1353 $me->{params}->{clist} = \@clist; 1354 1355 $me->{func} = sub { 1356 my ($data,$p) = @_; 1357 my ($ip) = $data->is_inplace; 1358 for my $t ( reverse @{$p->{clist}} ) { 1359 croak("Error: tried to apply a PDL::Transform with no function inside a composition!\n offending transform: $t\n") 1360 unless(defined($t->{func}) and ref($t->{func}) eq 'CODE'); 1361 $data = $t->{func}($ip ? $data->inplace : $data, $t->{params}); 1362 } 1363 $data->is_inplace(0); # clear inplace flag (avoid core bug with inplace) 1364 $data; 1365 }; 1366 1367 $me->{inv} = sub { 1368 my($data,$p) = @_; 1369 my($ip) = $data->is_inplace; 1370 for my $t ( @{$p->{clist}} ) { 1371 croak("Error: tried to invert a non-invertible PDL::Transform inside a composition!\n offending transform: $t\n") 1372 unless(defined($t->{inv}) and ref($t->{inv}) eq 'CODE'); 1373 $data = &{$t->{inv}}($ip ? $data->inplace : $data, $t->{params}); 1374 } 1375 $data; 1376 }; 1377 1378 return bless($me,'PDL::Transform::Composition'); 1379} 1380 1381 1382 1383 1384=head2 t_wrap 1385 1386=for usage 1387 1388 $g1fg = $f->wrap($g); 1389 $g1fg = t_wrap($f,$g); 1390 1391=for ref 1392 1393Shift a transform into a different space by 'wrapping' it with a second. 1394 1395This is just a convenience function for two 1396L<t_compose|/t_compose> calls. C<< $a->wrap($b) >> is the same as 1397C<(!$b) x $a x $b>: the resulting transform first hits the data with 1398$b, then with $a, then with the inverse of $b. 1399 1400For example, to shift the origin of rotation, do this: 1401 1402 $im = rfits('m51.fits'); 1403 $tf = t_fits($im); 1404 $tr = t_linear({rot=>30}); 1405 $im1 = $tr->map($tr); # Rotate around pixel origin 1406 $im2 = $tr->map($tr->wrap($tf)); # Rotate round FITS scientific origin 1407 1408=cut 1409 1410*t_wrap = \&wrap; 1411 1412sub wrap { 1413 my($f) = shift; 1414 my($g) = shift; 1415 1416 return $g->inverse->compose($f,$g); 1417} 1418 1419 1420 1421###################################################################### 1422 1423# Composition operator -- handles 'x'. 1424sub _compose_op { 1425 my($a,$b,$c) = @_; 1426 $c ? compose($b,$a) : compose($a,$b); 1427} 1428 1429# Raise-to-power operator -- handles '**'. 1430 1431sub _pow_op { 1432 my($a,$b,$c) = @_; 1433 1434 barf("%s", "Can't raise anything to the power of a transform") 1435 if($c || UNIVERSAL::isa($b,'PDL::Transform')) ; 1436 1437 $a = $a->inverse 1438 if($b < 0); 1439 1440 return $a if(abs($b) == 1); 1441 return new PDL::Transform if(abs($b) == 0); 1442 1443 my(@l); 1444 for my $i(1..abs($b)) { 1445 push(@l,$a); 1446 } 1447 1448 t_compose(@l); 1449} 1450 1451 1452 1453 1454=head2 t_identity 1455 1456=for usage 1457 1458 my $xform = t_identity 1459 my $xform = new PDL::Transform; 1460 1461=for ref 1462 1463Generic constructor generates the identity transform. 1464 1465This constructor really is trivial -- it is mainly used by the other transform 1466constructors. It takes no parameters and returns the identity transform. 1467 1468=cut 1469 1470sub _identity { return shift; } 1471sub t_identity { new PDL::Transform(@_) }; 1472 1473sub new { 1474 my($class) = shift; 1475 my $me = {name=>'identity', 1476 idim => 0, 1477 odim => 0, 1478 func=>\&PDL::Transform::_identity, 1479 inv=>\&PDL::Transform::_identity, 1480 params=>{} 1481 }; 1482 1483 return bless $me,$class; 1484} 1485 1486 1487 1488 1489=head2 t_lookup 1490 1491=for usage 1492 1493 $f = t_lookup($lookup, {<options>}); 1494 1495=for ref 1496 1497Transform by lookup into an explicit table. 1498 1499You specify an N+1-D PDL that is interpreted as an N-D lookup table of 1500column vectors (vector index comes last). The last dimension has 1501order equal to the output dimensionality of the transform. 1502 1503For added flexibility in data space, You can specify pre-lookup linear 1504scaling and offset of the data. Of course you can specify the 1505interpolation method to be used. The linear scaling stuff is a little 1506primitive; if you want more, try composing the linear transform with 1507this one. 1508 1509The prescribed values in the lookup table are treated as 1510pixel-centered: that is, if your input array has N elements per row 1511then valid data exist between the locations (-0.5) and (N-0.5) in 1512lookup pixel space, because the pixels (which are numbered from 0 to 1513N-1) are centered on their locations. 1514 1515Lookup is done using L<interpND|PDL::Primitive/interpnd>, so the boundary conditions 1516and threading behaviour follow from that. 1517 1518The indexed-over dimensions come first in the table, followed by a 1519single dimension containing the column vector to be output for each 1520set of other dimensions -- ie to output 2-vectors from 2 input 1521parameters, each of which can range from 0 to 49, you want an index 1522that has dimension list (50,50,2). For the identity lookup table 1523you could use C<cat(xvals(50,50),yvals(50,50))>. 1524 1525If you want to output a single value per input vector, you still need 1526that last index threading dimension -- if necessary, use C<dummy(-1,1)>. 1527 1528The lookup index scaling is: out = lookup[ (scale * data) + offset ]. 1529 1530A simplistic table inversion routine is included. This means that 1531you can (for example) use the C<map> method with C<t_lookup> transformations. 1532But the table inversion is exceedingly slow, and not practical for tables 1533larger than about 100x100. The inversion table is calculated in its 1534entirety the first time it is needed, and then cached until the object is 1535destroyed. 1536 1537Options are listed below; there are several synonyms for each. 1538 1539=over 3 1540 1541=item s, scale, Scale 1542 1543(default 1.0) Specifies the linear amount of scaling to be done before 1544lookup. You can feed in a scalar or an N-vector; other values may cause 1545trouble. If you want to save space in your table, then specify smaller 1546scale numbers. 1547 1548=item o, offset, Offset 1549 1550(default 0.0) Specifies the linear amount of offset before lookup. 1551This is only a scalar, because it is intended to let you switch to 1552corner-centered coordinates if you want to (just feed in o=-0.25). 1553 1554=item b, bound, boundary, Boundary 1555 1556Boundary condition to be fed to L<interpND|PDL::Primitive/interpND> 1557 1558=item m, method, Method 1559 1560Interpolation method to be fed to L<interpND|PDL::Primitive/interpND> 1561 1562=back 1563 1564EXAMPLE 1565 1566To scale logarithmically the Y axis of m51, try: 1567 1568 $a = float rfits('m51.fits'); 1569 $lookup = xvals(128,128) -> cat( 10**(yvals(128,128)/50) * 256/10**2.55 ); 1570 $t = t_lookup($lookup); 1571 $b = $t->map($a); 1572 1573To do the same thing but with a smaller lookup table, try: 1574 1575 $lookup = 16 * xvals(17,17)->cat(10**(yvals(17,17)/(100/16)) * 16/10**2.55); 1576 $t = t_lookup($lookup,{scale=>1/16.0}); 1577 $b = $t->map($a); 1578 1579(Notice that, although the lookup table coordinates are is divided by 16, 1580it is a 17x17 -- so linear interpolation works right to the edge of the original 1581domain.) 1582 1583NOTES 1584 1585Inverses are not yet implemented -- the best way to do it might be by 1586judicious use of map() on the forward transformation. 1587 1588the type/unit fields are ignored. 1589 1590=cut 1591 1592sub t_lookup { 1593 my($class) = 'PDL::Transform'; 1594 my($source)= shift; 1595 my($o) = shift; 1596 1597 if(!defined($o) && ((ref $source) eq 'HASH')) { 1598 Carp::cluck("lookup transform called as sub not method; using 'PDL::Transform' as class...\n"); 1599 $o = $source; 1600 $source = $class; 1601 $class = "PDL::Transform"; 1602 } 1603 1604 $o = {} unless(ref $o eq 'HASH'); 1605 1606 my($me) = PDL::Transform::new($class); 1607 1608 my($bound) = _opt($o,['b','bound','boundary','Boundary']); 1609 my($method)= _opt($o,['m','meth','method','Method']); 1610 1611 $me->{idim} = $source->ndims - 1; 1612 $me->{odim} = $source->dim($source->ndims-1); 1613 1614 $me->{params} = { 1615 table => $source, 1616 scale => _opt($o,['s','scale','Scale'],1.0), 1617 offset => _opt($o,['o','off','offset','Offset'],0.0), 1618 interpND_opt => { 1619 method => $method, 1620 bound => $bound, 1621 bad => _opt($o,['bad'],0) 1622 } 1623 }; 1624 1625 1626 my $lookup_func = sub { 1627 my($data,$p,$table,$scale,$offset) = @_; 1628 1629 $data = pdl($data) 1630 unless ((ref $data) && (UNIVERSAL::isa($data,'PDL'))); 1631 $main::foo = $data; 1632 1633 if($data->dim(0) > $me->{idim}) { 1634 barf("Too many dims (".$data->dim(0).") for your table (".$me->{idim}.")\n"); 1635 }; 1636 1637 my($a)= ($table 1638 ->interpND(float($data) * $scale + $offset, 1639 $p->{interpND_opt} 1640 ) 1641 ); 1642 1643 1644 # Put the index dimension (and threaded indices) back at the front of 1645 # the dimension list. 1646 my($dnd) = $data->ndims - 1; 1647 return ($a -> ndims > $data->ndims - 1) ? 1648 ($a->reorder( $dnd..($dnd + $table->ndims - $data->dim(0)-1) 1649 , 0..$data->ndims-2 1650 ) 1651 ) : $a; 1652 }; 1653 1654 $me->{func} = sub {my($data,$p) = @_; &$lookup_func($data,$p,$p->{table},$p->{scale},$p->{offset})}; 1655 1656 ####### 1657 ## Lazy inverse -- find it if and only if we need it... 1658 $me->{inv} = sub { 1659 my $data = shift; 1660 my $p = shift; 1661 if(!defined($p->{'itable'})) { 1662 if($me->{idim} != $me->{odim}) { 1663 barf("t_lookup: can't calculate an inverse of a projection operation! (idim != odim)"); 1664 } 1665 print "t_lookup: Warning, table inversion is only weakly supported. \n I'll try to invert it using a pretty boneheaded algorithm...\n (If it takes too long, consider writing a faster algorithm!)\n Calculating inverse table (will be cached)...\n" if($PDL::verbose || $PDL::debug || $PDL::Transform::debug); 1666 my $itable = zeroes($p->{table}); 1667 my $minvals = $p->{table}->clump($me->{idim})->minimum; 1668 my $maxvals = $p->{table}->clump($me->{idim})->maximum; 1669 1670 # Scale so that the range runs from 0 through the top pixel in the table 1671 my $scale = ( pdl( $itable->dims )->slice("0:-2")-1 ) / 1672 (($maxvals - $minvals)+ (($maxvals-$minvals) == 0)); 1673 my $offset = - ($minvals * $scale); 1674 1675 $p->{iscale} = $scale; 1676 $p->{ioffset} = $offset; 1677 1678 my $enl_scale = $p->{'enl_scale'} || 10; 1679 my $smallcoords = ndcoords((pdl($enl_scale * 2 + 1)->at(0)) x $me->{idim})/ $enl_scale - 1; 1680 1681 # $oloop runs over (point, index) for all points in the output table, in 1682 # $p->{table} output space 1683 $oloop = ndcoords($itable->mv(-1,0)->slice("(0)"))-> 1684 double-> 1685 mv(0,-1)-> 1686 clump($itable->ndims-1); # oloop: (pixel, index) 1687 { 1688 my $tmp; # work around perl -d "feature" 1689 ($tmp = $oloop->mv(-1,0)) -= $offset; 1690 ($tmp = $oloop->mv(-1,0)) /= $scale; 1691 } 1692 1693 # The Right Thing to do here is to take the outer product of $itable and $otable, then collapse 1694 # to find minimum distance. But memory demands for that would be HUGE. Instead we do an 1695 # elementwise calculation. 1696 1697 print "t_lookup: inverting ".$oloop->dim(0)." points...\n" if($PDL::verbose || $PDL::debug || $PDL::Transform::debug); 1698 my $pt = $p->{table}->mv(-1,0); # pt runs (index, x,y,...) 1699 1700 my $itable_flattened = zeroes($oloop); 1701 1702 for $i(0..$oloop->dim(0)-1) { 1703 1704 my $olp = $oloop->slice("($i)"); # olp runs (index) 1705 my $diff = ($pt - $olp); # diff runs (index, x, y, ...) 1706 my $r2 = ($diff * $diff)->sumover; # r2 runs (x,y,...) 1707 my $c = whichND($r2==$r2->min)->slice(":,(0)"); # c runs (index) 1708 1709 # Now zero in on the neighborhood around the point of closest approach. 1710 my $neighborhood = $p->{table}->interpND($c + $smallcoords,{method=>'linear',bound=>'t'})-> 1711 mv(-1,0); # neighborhood runs (index, dx, dy,...) 1712 $diff = $neighborhood - $olp; # diff runs (index, dx, dy, ...) 1713 $r2 = ($diff * $diff)->sumover; # r2 runs (dx,dy,...) 1714 my $dc = $smallcoords->mv(0,-1)->indexND(0+whichND($r2==$r2->min)->slice(":,(0)")); 1715 1716 1717 my $coord = $c + $dc; 1718 # At last, we've found the best-fit to an enl_scale'th of an input-table pixel. 1719 # Back it out to input-science coordinates, and stuff it in the inverse table. 1720 my $tmp; # work around perl -d "feature" 1721 ($tmp = $itable_flattened->slice("($i)")) .= $coord; 1722 1723 print " $i..." if( ($i%1000==0) && ( $PDL::verbose || $PDL::debug || $PDL::Transform::debug)); 1724 } 1725 1726 { 1727 my $tmp; # work around perl -d "feature" 1728 ($tmp = $itable->clump($itable->ndims-1)) .= $itable_flattened; 1729 ($tmp = $itable->mv(-1,0)) -= $p->{offset}; 1730 ($tmp = $itable->mv(-1,0)) /= $p->{scale}; 1731 } 1732 1733 $p->{itable} = $itable; 1734 } 1735 &$lookup_func($data,$p, $p->{itable},$p->{iscale},$p->{ioffset}) ; 1736 }; 1737 1738 1739 $me->{name} = 'Lookup'; 1740 1741 return $me; 1742} 1743 1744 1745 1746 1747=head2 t_linear 1748 1749=for usage 1750 1751$f = t_linear({options}); 1752 1753=for ref 1754 1755Linear (affine) transformations with optional offset 1756 1757t_linear implements simple matrix multiplication with offset, 1758also known as the affine transformations. 1759 1760You specify the linear transformation with pre-offset, a mixing 1761matrix, and a post-offset. That overspecifies the transformation, so 1762you can choose your favorite method to specify the transform you want. 1763The inverse transform is automagically generated, provided that it 1764actually exists (the transform matrix is invertible). Otherwise, the 1765inverse transform just croaks. 1766 1767Extra dimensions in the input vector are ignored, so if you pass a 17683xN vector into a 3-D linear transformation, the final dimension is passed 1769through unchanged. 1770 1771The options you can usefully pass in are: 1772 1773=over 3 1774 1775=item s, scale, Scale 1776 1777A scaling scalar (heh), vector, or matrix. If you specify a vector 1778it is treated as a diagonal matrix (for convenience). It gets 1779left-multiplied with the transformation matrix you specify (or the 1780identity), so that if you specify both a scale and a matrix the 1781scaling is done after the rotation or skewing or whatever. 1782 1783=item r, rot, rota, rotation, Rotation 1784 1785A rotation angle in degrees -- useful for 2-D and 3-D data only. If 1786you pass in a scalar, it specifies a rotation from the 0th axis toward 1787the 1st axis. If you pass in a 3-vector as either a PDL or an array 1788ref (as in "rot=>[3,4,5]"), then it is treated as a set of Euler 1789angles in three dimensions, and a rotation matrix is generated that 1790does the following, in order: 1791 1792=over 3 1793 1794=item * Rotate by rot->(2) degrees from 0th to 1st axis 1795 1796=item * Rotate by rot->(1) degrees from the 2nd to the 0th axis 1797 1798=item * Rotate by rot->(0) degrees from the 1st to the 2nd axis 1799 1800=back 1801 1802The rotation matrix is left-multiplied with the transformation matrix 1803you specify, so that if you specify both rotation and a general matrix 1804the rotation happens after the more general operation -- though that is 1805deprecated. 1806 1807Of course, you can duplicate this functionality -- and get more 1808general -- by generating your own rotation matrix and feeding it in 1809with the C<matrix> option. 1810 1811=item m, matrix, Matrix 1812 1813The transformation matrix. It does not even have to be square, if you want 1814to change the dimensionality of your input. If it is invertible (note: 1815must be square for that), then you automagically get an inverse transform too. 1816 1817=item pre, preoffset, offset, Offset 1818 1819The vector to be added to the data before they get multiplied by the matrix 1820(equivalent of CRVAL in FITS, if you are converting from scientific to 1821pixel units). 1822 1823=item post, postoffset, shift, Shift 1824 1825The vector to be added to the data after it gets multiplied by the matrix 1826(equivalent of CRPIX-1 in FITS, if youre converting from scientific to pixel 1827units). 1828 1829=item d, dim, dims, Dims 1830 1831Most of the time it is obvious how many dimensions you want to deal 1832with: if you supply a matrix, it defines the transformation; if you 1833input offset vectors in the C<pre> and C<post> options, those define 1834the number of dimensions. But if you only supply scalars, there is no way 1835to tell and the default number of dimensions is 2. This provides a way 1836to do, e.g., 3-D scaling: just set C<{s=><scale-factor>, dims=>3}> and 1837you are on your way. 1838 1839=back 1840 1841NOTES 1842 1843the type/unit fields are currently ignored by t_linear. 1844 1845=cut 1846 1847@PDL::Transform::Linear::ISA = ('PDL::Transform'); 1848 1849sub t_linear { new PDL::Transform::Linear(@_); } 1850 1851sub PDL::Transform::Linear::new { 1852 my($class) = shift; 1853 my($o) = $_[0]; 1854 pop @_ if (($#_ % 2 ==0) && !defined($_[-1])); 1855 #suppresses a warning if @_ has an odd number of elements and the 1856 #last is undef 1857 1858 if(!(ref $o)) { 1859 $o = {@_}; 1860 } 1861 1862 my($me) = PDL::Transform::new($class); 1863 1864 $me->{name} = "linear"; 1865 1866 $me->{params}->{pre} = _opt($o,['pre','Pre','preoffset','offset', 1867 'Offset','PreOffset','Preoffset'],0); 1868 $me->{params}->{pre} = pdl($me->{params}->{pre}) 1869 if(defined $me->{params}->{pre}); 1870 1871 $me->{params}->{post} = _opt($o,['post','Post','postoffset','PostOffset', 1872 'shift','Shift'],0); 1873 $me->{params}->{post} = pdl($me->{params}->{post}) 1874 if(defined $me->{params}->{post}); 1875 1876 $me->{params}->{matrix} = _opt($o,['m','matrix','Matrix','mat','Mat']); 1877 $me->{params}->{matrix} = pdl($me->{params}->{matrix}) 1878 if(defined $me->{params}->{matrix}); 1879 1880 $me->{params}->{rot} = _opt($o,['r','rot','rota','rotation','Rotation']); 1881 $me->{params}->{rot} = 0 unless defined($me->{params}->{rot}); 1882 $me->{params}->{rot} = pdl($me->{params}->{rot}); 1883 1884 my $o_dims = _opt($o,['d','dim','dims','Dims']); 1885 $o_dims = pdl($o_dims) 1886 if defined($o_dims); 1887 1888 my $scale = _opt($o,['s','scale','Scale']); 1889 $scale = pdl($scale) 1890 if defined($scale); 1891 1892 # Figure out the number of dimensions to transform, and, 1893 # if necessary, generate a new matrix. 1894 1895 if(defined($me->{params}->{matrix})) { 1896 my $mat = $me->{params}->{matrix} = $me->{params}->{matrix}->slice(":,:"); 1897 $me->{idim} = $mat->dim(0); 1898 $me->{odim} = $mat->dim(1); 1899 1900 } else { 1901 if(defined($me->{params}->{rot}) && 1902 UNIVERSAL::isa($me->{params}->{rot},'PDL')) { 1903 $me->{idim} = $me->{odim} = 2 if($me->{params}->{rot}->nelem == 1); 1904 $me->{idim} = $me->{odim} = 3 if($me->{params}->{rot}->nelem == 3); 1905 } 1906 1907 if(defined($scale) && 1908 UNIVERSAL::isa($scale,'PDL') && 1909 $scale->getndims > 0) { 1910 $me->{idim} = $me->{odim} = $scale->dim(0); 1911 $me->{odim} = $scale->dim(0); 1912 1913 } elsif(defined($me->{params}->{pre}) && 1914 UNIVERSAL::isa($me->{params}->{pre},'PDL') && 1915 $me->{params}->{pre}->getndims > 0) { 1916 $me->{idim} = $me->{odim} = $me->{params}->{pre}->dim(0); 1917 1918 } elsif(defined($me->{params}->{post}) && 1919 UNIVERSAL::isa($me->{params}->{post},'PDL') && 1920 $me->{params}->{post}->getndims > 0) { 1921 $me->{idim} = $me->{odim} = $me->{params}->{post}->dim(0); 1922 } elsif(defined($o_dims)) { 1923 $me->{idim} = $me->{odim} = $o_dims; 1924 } else { 1925 print "PDL::Transform::Linear: assuming 2-D transform (set dims option to change)\n" if($PDL::Transform::debug); 1926 $me->{idim} = $me->{odim} = 2; 1927 } 1928 1929 $me->{params}->{matrix} = PDL->zeroes($me->{idim},$me->{odim}); 1930 my $tmp; # work around perl -d "feature" 1931 ($tmp = $me->{params}->{matrix}->diagonal(0,1)) .= 1; 1932 1933 } 1934 1935 ### Handle rotation option 1936 my $rot = $me->{params}->{rot}; 1937 if(defined($rot)) { 1938 # Subrotation closure -- rotates from axis $d->(0) --> $d->(1). 1939 my $subrot = sub { 1940 my($d,$angle,$m)=@_; 1941 my($i) = identity($m->dim(0)); 1942 my($subm) = $i->dice($d,$d); 1943 1944 $angle = $angle->at(0) 1945 if(UNIVERSAL::isa($angle,'PDL')); 1946 1947 my($a) = $angle * $DEG2RAD; 1948 $subm .= $subm x pdl([cos($a),-sin($a)],[sin($a),cos($a)]); 1949 $m .= $m x $i; 1950 }; 1951 1952 if(UNIVERSAL::isa($rot,'PDL') && $rot->nelem > 1) { 1953 if($rot->ndims == 2) { 1954 $me->{params}->{matrix} x= $rot; 1955 } elsif($rot->nelem == 3) { 1956 my $rm = identity(3); 1957 1958 # Do these in reverse order to make it more like 1959 # function composition! 1960 &$subrot(pdl(0,1),$rot->at(2),$rm); 1961 &$subrot(pdl(2,0),$rot->at(1),$rm); 1962 &$subrot(pdl(1,2),$rot->at(0),$rm); 1963 1964 $me->{params}->{matrix} .= $me->{params}->{matrix} x $rm; 1965 } else { 1966 barf("PDL::Transform::Linear: Got a strange rot option -- giving up.\n"); 1967 } 1968 } else { 1969 if($rot != 0 and $me->{params}->{matrix}->dim(0)>1) { 1970 &$subrot(pdl(0,1),$rot,$me->{params}->{matrix}); 1971 } 1972 } 1973 } 1974 1975 1976 # 1977 # Apply scaling 1978 # 1979 $me->{params}->{matrix} = $me->{params}->{matrix}->slice(":,:"); 1980 my $tmp; # work around perl -d "feature" 1981 ($tmp = $me->{params}->{matrix}->diagonal(0,1)) *= $scale 1982 if defined($scale); 1983 1984 # 1985 # Check for an inverse and apply it if possible 1986 # 1987 my($o2); 1988 if($me->{params}->{matrix}->det($o2 = {lu=>undef})) { 1989 $me->{params}->{inverse} = $me->{params}->{matrix}->inv($o2); 1990 } else { 1991 delete $me->{params}->{inverse}; 1992 } 1993 1994 $me->{params}->{idim} = $me->{idim}; 1995 $me->{params}->{odim} = $me->{odim}; 1996 1997 1998 ############################## 1999 # The meat -- just shift, matrix-multiply, and shift again. 2000 $me->{func} = sub { 2001 my($in,$opt) = @_; 2002 2003 my($d) = $opt->{matrix}->dim(0)-1; 2004 2005 barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n") 2006 if($in->dim(0) < $d); 2007 2008 my($a) = $in->slice("0:$d")->copy + $opt->{pre}; 2009 my($out) = $in->is_inplace ? $in : $in->copy; 2010 2011 my $tmp; # work around perl -d "feature" 2012 ($tmp = $out->slice("0:$d")) .= $a x $opt->{matrix} + $opt->{post}; 2013 2014 return $out; 2015 }; 2016 2017 2018 $me->{inv} = (defined $me->{params}->{inverse}) ? sub { 2019 my($in,$opt) = @_; 2020 2021 my($d) = $opt->{inverse}->dim(0)-1; 2022 barf("Linear transform: transform is $d-D; data only ".($in->dim(0))."\n") 2023 if($in->dim(0) < $d); 2024 2025 my($a) = $in->slice("0:$d")->copy - $opt->{post}; 2026 my($out) = $in->is_inplace ? $in : $in->copy; 2027 2028 my $tmp; # work around perl -d "feature" 2029 ($tmp = $out->slice("0:$d")) .= $a x $opt->{inverse} - $opt->{pre}; 2030 2031 $out; 2032 } : undef; 2033 2034 return $me; 2035} 2036 2037sub PDL::Transform::Linear::stringify { 2038 package PDL::Transform::Linear; 2039 my($me) = shift; my($out) = SUPER::stringify $me; 2040 my $mp = $me->{params}; 2041 2042 if(!($me->{is_inverse})){ 2043 $out .= "Pre-add: ".($mp->{pre})."\n" 2044 if(defined $mp->{pre}); 2045 2046 $out .= "Post-add: ".($mp->{post})."\n" 2047 if(defined $mp->{post}); 2048 2049 $out .= "Forward matrix:".($mp->{matrix}) 2050 if(defined $mp->{matrix}); 2051 2052 $out .= "Inverse matrix:".($mp->{inverse}) 2053 if(defined $mp->{inverse}); 2054 } else { 2055 $out .= "Pre-add: ".(-$mp->{post})."\n" 2056 if(defined $mp->{post}); 2057 2058 $out .= "Post-add: ".(-$mp->{pre})."\n" 2059 if(defined $mp->{pre}); 2060 2061 $out .= "Forward matrix:".($mp->{inverse}) 2062 if(defined $mp->{inverse}); 2063 2064 $out .= "Inverse matrix:".($mp->{matrix}) 2065 if(defined $mp->{matrix}); 2066 } 2067 2068 $out =~ s/\n/\n /go; 2069 $out; 2070} 2071 2072 2073 2074 2075=head2 t_scale 2076 2077=for usage 2078 2079 $f = t_scale(<scale>) 2080 2081=for ref 2082 2083Convenience interface to L<t_linear|/t_linear>. 2084 2085t_scale produces a transform that scales around the origin by a fixed 2086amount. It acts exactly the same as C<t_linear(Scale=>\<scale\>)>. 2087 2088=cut 2089 2090sub t_scale { 2091 my($scale) = shift; 2092 my($b) = shift; 2093 return t_linear(scale=>$scale,%{$b}) 2094 if(ref $b eq 'HASH'); 2095 t_linear(Scale=>$scale,$b,@_); 2096} 2097 2098 2099 2100 2101=head2 t_offset 2102 2103=for usage 2104 2105 $f = t_offset(<shift>) 2106 2107=for ref 2108 2109Convenience interface to L<t_linear|/t_linear>. 2110 2111t_offset produces a transform that shifts the origin to a new location. 2112It acts exactly the same as C<t_linear(Pre=>\<shift\>)>. 2113 2114=cut 2115 2116sub t_offset { 2117 my($pre) = shift; 2118 my($b) = shift; 2119 return t_linear(pre=>$pre,%{$b}) 2120 if(ref $b eq 'HASH'); 2121 2122 t_linear(pre=>$pre,$b,@_); 2123} 2124 2125 2126 2127 2128=head2 t_rot 2129 2130=for usage 2131 2132 $f = t_rot(<rotation-in-degrees>) 2133 2134=for ref 2135 2136Convenience interface to L<t_linear|/t_linear>. 2137 2138t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or 2139N-D (matrix). It acts exactly the same as C<t_linear(Rot=>\<shift\>)>. 2140 2141=cut 2142 2143*t_rot = \&t_rotate; 2144sub t_rotate { 2145 my $rot = shift; 2146 my($b) = shift; 2147 return t_linear(rot=>$rot,%{$b}) 2148 if(ref $b eq 'HASH'); 2149 2150 t_linear(rot=>$rot,$b,@_); 2151} 2152 2153 2154 2155 2156=head2 t_fits 2157 2158=for usage 2159 2160 $f = t_fits($fits,[option]); 2161 2162=for ref 2163 2164FITS pixel-to-scientific transformation with inverse 2165 2166You feed in a hash ref or a PDL with one of those as a header, and you 2167get back a transform that converts 0-originated, pixel-centered 2168coordinates into scientific coordinates via the transformation in the 2169FITS header. For most FITS headers, the transform is reversible, so 2170applying the inverse goes the other way. This is just a convenience 2171subclass of PDL::Transform::Linear, but with unit/type support 2172using the FITS header you supply. 2173 2174For now, this transform is rather limited -- it really ought to 2175accept units differences and stuff like that, but they are just 2176ignored for now. Probably that would require putting units into 2177the whole transform framework. 2178 2179This transform implements the linear transform part of the WCS FITS 2180standard outlined in Greisen & Calabata 2002 (A&A in press; find it at 2181"http://arxiv.org/abs/astro-ph/0207407"). 2182 2183As a special case, you can pass in the boolean option "ignore_rgb" 2184(default 0), and if you pass in a 3-D FITS header in which the last 2185dimension has exactly 3 elements, it will be ignored in the output 2186transformation. That turns out to be handy for handling rgb images. 2187 2188=cut 2189 2190sub t_fits { 2191 my($class) = 'PDL::Transform::Linear'; 2192 my($hdr) = shift; 2193 my($opt) = shift; 2194 2195 if(ref $opt ne 'HASH') { 2196 $opt = defined $opt ? {$opt,@_} : {} ; 2197 } 2198 2199 $hdr = $hdr->gethdr 2200 if(defined $hdr && UNIVERSAL::isa($hdr,'PDL')); 2201 2202 croak('PDL::Transform::FITS::new requires a FITS header hash\n') 2203 if(!defined $hdr || ref $hdr ne 'HASH' || !defined($hdr->{NAXIS})); 2204 2205 my($n) = $hdr->{NAXIS}; $n = $n->at(0) if(UNIVERSAL::isa($n,'PDL')); 2206 2207 $n = 2 2208 if($opt->{ignore_rgb} && $n==3 && $hdr->{NAXIS3} == 3); 2209 2210 my($matrix) = PDL->zeroes($hdr->{NAXIS},$hdr->{NAXIS}); 2211 my($pre) = PDL->zeroes($n); 2212 my($post) = PDL->zeroes($n); 2213 2214 ############################## 2215 # Scaling: Use CDi_j formalism if present; otherwise use the 2216 # older PCi_j + CDELTi formalism. 2217 2218 my(@hgrab); 2219 2220 if(@hgrab = grep(m/^CD\d{1,3}_\d{1,3}$/,keys %$hdr)) { # assignment 2221 # 2222 # CDi_j formalism 2223 # 2224 for my $h(@hgrab) { 2225 $h =~ m/CD(\d{1,3})_(\d{1,3})/; # Should always match 2226 my $tmp; # work around perl -d "feature" 2227 ($tmp = $matrix->slice("(".($1-1)."),(".($2-1).")")) .= $hdr->{$h}; 2228 } 2229 print "PDL::Transform::FITS: Detected CDi_j matrix: \n",$matrix,"\n" 2230 if($PDL::Transform::debug); 2231 2232 } else { 2233 2234 # 2235 # PCi_j + CDELTi formalism 2236 # If PCi_j aren't present, and N=2, then try using CROTA or 2237 # CROTA2 to generate a rotation matrix instea. 2238 # 2239 2240 my($cdm) = PDL->zeroes($n,$n); 2241 my($cd) = $cdm->diagonal(0,1); 2242 2243 my($cpm) = PDL->zeroes($n,$n); 2244 my $tmp; # work around perl -d "feature" 2245 ($tmp = $cpm->diagonal(0,1)) .= 1; # PC: diagonal defaults to unity 2246 $cd .= 1; 2247 2248 2249 if( @hgrab = grep(m/^PC\d{1,3}_\d{1,3}$/,keys %$hdr) ) { # assignment 2250 2251 for my $h(@hgrab) { 2252 $h =~ m/PC(\d{1,3})_(\d{1,3})$/ || die "t_fits - match failed! This should never happen!"; 2253 my $tmp; # work around perl -d "feature" 2254 ($tmp = $cpm->slice("(".($1-1)."),(".($2-1).")")) .= $hdr->{$h}; 2255 } 2256 print "PDL::Transform::FITS: Detected PCi_j matrix: \n",$cpm,"\n" 2257 if($PDL::Transform::debug && @hgrab); 2258 2259 } elsif($n==2 && ( defined $hdr->{CROTA} || defined $hdr->{CROTA1} || defined $hdr->{CROTA2}) ) { 2260 2261 ## CROTA is deprecated; CROTA1 was used for a while but is unofficial; 2262 ## CROTA2 is encouraged instead. 2263 my $cr; 2264 $cr = $hdr->{CROTA2} unless defined $cr; 2265 $cr = $hdr->{CROTA} unless defined $cr; 2266 $cr = $hdr->{CROTA1} unless defined $cr; 2267 2268 $cr *= $DEG2RAD; 2269 # Rotation matrix rotates counterclockwise to get from sci to pixel coords 2270 # (detector has been rotated ccw, according to FITS standard) 2271 $cpm .= pdl( [cos($cr), sin($cr)],[-sin($cr),cos($cr)] ); 2272 2273 } 2274 2275 for my $i(1..$n) { 2276 my $tmp; # work around perl -d "feature" 2277 ($tmp = $cd->slice("(".($i-1).")")) .= $hdr->{"CDELT$i"}; 2278 } 2279#If there are no CDELTs, then we assume they are all 1.0, 2280#as in PDL::Graphics::PGPLOT::Window::_FITS_tr. 2281 $cd+=1 if (all($cd==0)); 2282 2283 $matrix = $cdm x $cpm; 2284 } 2285 2286 my($i1) = 0; 2287 for my $i(1..$n) { 2288 my $tmp; # work around perl -d "feature" 2289 ($tmp = $pre->slice("($i1)")) .= 1 - $hdr->{"CRPIX$i"}; 2290 ($tmp = $post->slice("($i1)")) .= $hdr->{"CRVAL$i"}; 2291 $i1++; 2292 } 2293 2294 my($me) = PDL::Transform::Linear::new($class, 2295 {'pre'=>$pre, 2296 'post'=>$post, 2297 'matrix'=>$matrix 2298 }); 2299 2300 $me->{name} = 'FITS'; 2301 2302 my (@otype,@ounit,@itype,@iunit); 2303 our (@names) = ('X','Y','Z') unless @names; 2304 2305 for my $i(1..$hdr->{NAXIS}) { 2306 push(@otype,$hdr->{"CTYPE$i"}); 2307 push(@ounit,$hdr->{"CUNIT$i"}); 2308 push(@itype,"Image ". ( ($i-1<=$#names) ? $names[$i-1] : "${i}th dim" )); 2309 push(@iunit,"Pixels"); 2310 } 2311 2312 $me->{otype} = \@otype; 2313 $me->{itype} = \@itype; 2314 $me->{ounit} = \@ounit; 2315 $me->{iunit} = \@iunit; 2316 2317 # Check for nonlinear projection... 2318# if($hdr->{CTYPE1} =~ m/(\w\w\w\w)\-(\w\w\w)/) { 2319# print "Nonlinear transformation found... ignoring nonlinear part...\n"; 2320# } 2321 2322 return $me; 2323 2324 2325} 2326 2327 2328 2329 2330=head2 t_code 2331 2332=for usage 2333 2334 $f = t_code(<func>,[<inv>],[options]); 2335 2336=for ref 2337 2338Transform implementing arbitrary perl code. 2339 2340This is a way of getting quick-and-dirty new transforms. You pass in 2341anonymous (or otherwise) code refs pointing to subroutines that 2342implement the forward and, optionally, inverse transforms. The 2343subroutines should accept a data PDL followed by a parameter hash ref, 2344and return the transformed data PDL. The parameter hash ref can be 2345set via the options, if you want to. 2346 2347Options that are accepted are: 2348 2349=over 3 2350 2351=item p,params 2352 2353The parameter hash that will be passed back to your code (defaults to the 2354empty hash). 2355 2356=item n,name 2357 2358The name of the transform (defaults to "code"). 2359 2360=item i, idim (default 2) 2361 2362The number of input dimensions (additional ones should be passed through 2363unchanged) 2364 2365=item o, odim (default 2) 2366 2367The number of output dimensions 2368 2369=item itype 2370 2371The type of the input dimensions, in an array ref (optional and advisiory) 2372 2373=item otype 2374 2375The type of the output dimension, in an array ref (optional and advisory) 2376 2377=item iunit 2378 2379The units that are expected for the input dimensions (optional and advisory) 2380 2381=item ounit 2382 2383The units that are returned in the output (optional and advisory). 2384 2385=back 2386 2387The code variables are executable perl code, either as a code ref or 2388as a string that will be eval'ed to produce code refs. If you pass in 2389a string, it gets eval'ed at call time to get a code ref. If it compiles 2390OK but does not return a code ref, then it gets re-evaluated with "sub { 2391... }" wrapped around it, to get a code ref. 2392 2393Note that code callbacks like this can be used to do really weird 2394things and generate equally weird results -- caveat scriptor! 2395 2396=cut 2397 2398sub t_code { 2399 my($class) = 'PDL::Transform'; 2400 my($func, $inv, $o) = @_; 2401 if(ref $inv eq 'HASH') { 2402 $o = $inv; 2403 $inv = undef; 2404 } 2405 2406 my($me) = PDL::Transform::new($class); 2407 $me->{name} = _opt($o,['n','name','Name']) || "code"; 2408 $me->{func} = $func; 2409 $me->{inv} = $inv; 2410 $me->{params} = _opt($o,['p','params','Params']) || {}; 2411 $me->{idim} = _opt($o,['i','idim']) || 2; 2412 $me->{odim} = _opt($o,['o','odim']) || 2; 2413 $me->{itype} = _opt($o,['itype']) || []; 2414 $me->{otype} = _opt($o,['otype']) || []; 2415 $me->{iunit} = _opt($o,['iunit']) || []; 2416 $me->{ounit} = _opt($o,['ounit']) || []; 2417 2418 $me; 2419} 2420 2421 2422 2423 2424=head2 t_cylindrical 2425 2426C<t_cylindrical> is an alias for C<t_radial> 2427 2428=head2 t_radial 2429 2430=for usage 2431 2432 $f = t_radial(<options>); 2433 2434=for ref 2435 2436Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with inverse) 2437 2438Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose 2439direct or conformal conversion. Direct conversion preserves radial 2440distance from the origin; conformal conversion preserves local angles, 2441so that each small-enough part of the image only appears to be scaled 2442and rotated, not stretched. Conformal conversion puts the radius on a 2443logarithmic scale, so that scaling of the original image plane is 2444equivalent to a simple offset of the transformed image plane. 2445 2446If you use three or more dimensions, the higher dimensions are ignored, 2447yielding a conversion from Cartesian to cylindrical coordinates, which 2448is why there are two aliases for the same transform. If you use higher 2449dimensionality than 2, you must manually specify the origin or you will 2450get dimension mismatch errors when you apply the transform. 2451 2452Theta runs B<clockwise> instead of the more usual counterclockwise; that is 2453to preserve the mirror sense of small structures. 2454 2455OPTIONS: 2456 2457=over 3 2458 2459=item d, direct, Direct 2460 2461Generate (theta,r) coordinates out (this is the default); incompatible 2462with Conformal. Theta is in radians, and the radial coordinate is 2463in the units of distance in the input plane. 2464 2465=item r0, c, conformal, Conformal 2466 2467If defined, this floating-point value causes t_radial to generate 2468(theta, ln(r/r0)) coordinates out. Theta is in radians, and the 2469radial coordinate varies by 1 for each e-folding of the r0-scaled 2470distance from the input origin. The logarithmic scaling is useful for 2471viewing both large and small things at the same time, and for keeping 2472shapes of small things preserved in the image. 2473 2474=item o, origin, Origin [default (0,0,0)] 2475 2476This is the origin of the expansion. Pass in a PDL or an array ref. 2477 2478=item u, unit, Unit [default 'radians'] 2479 2480This is the angular unit to be used for the azimuth. 2481 2482=back 2483 2484EXAMPLES 2485 2486These examples do transformations back into the same size image as they 2487started from; by suitable use of the "transform" option to 2488L<unmap|/unmap> you can send them to any size array you like. 2489 2490Examine radial structure in M51: 2491Here, we scale the output to stretch 2*pi radians out to the 2492full image width in the horizontal direction, and to stretch 1 radius out 2493to a diameter in the vertical direction. 2494 2495 $a = rfits('m51.fits'); 2496 $ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image 2497 $tu = t_radial(o => [130,130]); # Expand around galactic core 2498 $b = $a->map($ts x $tu); 2499 2500Examine radial structure in M51 (conformal): 2501Here, we scale the output to stretch 2*pi radians out to the full image width 2502in the horizontal direction, and scale the vertical direction by the exact 2503same amount to preserve conformality of the operation. Notice that 2504each piece of the image looks "natural" -- only scaled and not stretched. 2505 2506 $a = rfits('m51.fits') 2507 $ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale. 2508 $tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image 2509 $b = $ts->compose($tu)->unmap($a); 2510 2511 2512=cut 2513 2514*t_cylindrical = \&t_radial; 2515sub t_radial { 2516 my($class) = 'PDL::Transform'; 2517 my($o) = $_[0]; 2518 if(ref $o ne 'HASH') { 2519 $o = { @_ }; 2520 } 2521 2522 my($me) = PDL::Transform::new($class); 2523 2524 $me->{params}->{origin} = _opt($o,['o','origin','Origin']); 2525 $me->{params}->{origin} = pdl(0,0) 2526 unless defined($me->{params}->{origin}); 2527 $me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); 2528 2529 2530 $me->{params}->{r0} = _opt($o,['r0','R0','c','conformal','Conformal']); 2531 $me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); 2532 2533 $me->{params}->{u} = _opt($o,['u','unit','Unit'],'radians'); 2534 ### Replace this kludge with a units call 2535 $me->{params}->{angunit} = ($me->{params}->{u} =~ m/^d/i) ? $RAD2DEG : 1.0; 2536 print "radial: conversion is $me->{params}->{angunit}\n" if($PDL::Transform::debug); 2537 2538 $me->{name} = "radial (direct)"; 2539 2540 $me->{idim} = 2; 2541 $me->{odim} = 2; 2542 2543 if($me->{params}->{r0}) { 2544 $me->{otype} = ["Azimuth", "Ln radius" . ($me->{params}->{r0} != 1.0 ? "/$me->{params}->{r0}" : "")]; 2545 $me->{ounit} = [$me->{params}->{u},'']; # true-but-null prevents copying 2546 } else { 2547 $me->{otype} = ["Azimuth","Radius"]; 2548 $me->{ounit} = [$me->{params}->{u},'']; # false value copies prev. unit 2549 } 2550 2551 $me->{func} = sub { 2552 2553 my($data,$o) = @_; 2554 2555 my($out) = ($data->new_or_inplace); 2556 2557 my($d) = $data->copy; 2558 my $tmp; # work around perl -d "feature" 2559 ($tmp = $d->slice("0:1")) -= $o->{origin}; 2560 2561 my($d0) = $d->slice("(0)"); 2562 my($d1) = $d->slice("(1)"); 2563 2564 # (mod operator on atan2 puts everything in the interval [0,2*PI).) 2565 ($tmp = $out->slice("(0)")) .= (atan2(-$d1,$d0) % (2*$PI)) * $me->{params}->{angunit}; 2566 2567 ($tmp = $out->slice("(1)")) .= (defined $o->{r0}) ? 2568 0.5 * log( ($d1*$d1 + $d0 * $d0) / ($o->{r0} * $o->{r0}) ) : 2569 sqrt($d1*$d1 + $d0*$d0); 2570 2571 $out; 2572 }; 2573 2574 $me->{inv} = sub { 2575 2576 my($d,$o) = @_; 2577 my($d0,$d1,$out)= 2578 ( ($d->is_inplace) ? 2579 ($d->slice("(0)")->copy, $d->slice("(1)")->copy->dummy(0,2), $d) : 2580 ($d->slice("(0)"), $d->slice("(1)")->dummy(0,2), $d->copy) 2581 ); 2582 2583 $d0 /= $me->{params}->{angunit}; 2584 2585 my($os) = $out->slice("0:1"); 2586 $os .= append(cos($d0)->dummy(0,1),-sin($d0)->dummy(0,1)); 2587 $os *= defined $o->{r0} ? ($o->{r0} * exp($d1)) : $d1; 2588 $os += $o->{origin}; 2589 2590 $out; 2591 }; 2592 2593 2594 $me; 2595} 2596 2597 2598 2599 2600=head2 t_quadratic 2601 2602=for usage 2603 2604 $t = t_quadratic(<options>); 2605 2606=for ref 2607 2608Quadratic scaling -- cylindrical pincushion (n-d; with inverse) 2609 2610Quadratic scaling emulates pincushion in a cylindrical optical system: 2611separate quadratic scaling is applied to each axis. You can apply 2612separate distortion along any of the principal axes. If you want 2613different axes, use L<t_wrap|/t_wrap> and L<t_linear|/t_linear> to rotate 2614them to the correct angle. The scaling options may be scalars or 2615vectors; if they are scalars then the expansion is isotropic. 2616 2617The formula for the expansion is: 2618 2619 f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1) 2620 2621where <strength> is a scaling coefficient and <L_0> is a fundamental 2622length scale. Negative values of <strength> result in a pincushion 2623contraction. 2624 2625Note that, because quadratic scaling does not have a strict inverse for 2626coordinate systems that cross the origin, we cheat slightly and use 2627$x * abs($x) rather than $x**2. This does the Right thing for pincushion 2628and barrel distortion, but means that t_quadratic does not behave exactly 2629like t_cubic with a null cubic strength coefficient. 2630 2631OPTIONS 2632 2633=over 3 2634 2635=item o,origin,Origin 2636 2637The origin of the pincushion. (default is the, er, origin). 2638 2639=item l,l0,length,Length,r0 2640 2641The fundamental scale of the transformation -- the radius that remains 2642unchanged. (default=1) 2643 2644=item s,str,strength,Strength 2645 2646The relative strength of the pincushion. (default = 0.1) 2647 2648=item honest (default=0) 2649 2650Sets whether this is a true quadratic coordinate transform. The more 2651common form is pincushion or cylindrical distortion, which switches 2652branches of the square root at the origin (for symmetric expansion). 2653Setting honest to a non-false value forces true quadratic behavior, 2654which is not mirror-symmetric about the origin. 2655 2656=item d, dim, dims, Dims 2657 2658The number of dimensions to quadratically scale (default is the 2659dimensionality of your input vectors) 2660 2661 2662=back 2663 2664=cut 2665 2666sub t_quadratic { 2667 my($class) = 'PDL::Transform'; 2668 my($o) = $_[0]; 2669 if(ref $o ne 'HASH') { 2670 $o = {@_}; 2671 } 2672 my($me) = PDL::Transform::new($class); 2673 2674 $me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0)); 2675 $me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1)); 2676 $me->{params}->{str} = _opt($o,['s','str','strength','Strength'],pdl(0.1)); 2677 $me->{params}->{dim} = _opt($o,['d','dim','dims','Dims']); 2678 $me->{name} = "quadratic"; 2679 2680 $me->{func} = sub { 2681 my($data,$o) = @_; 2682 my($dd) = $data->copy - $o->{origin}; 2683 my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd; 2684 $d += $o->{str} * ($d * abs($d)) / $o->{l0}; 2685 $d /= (abs($o->{str}) + 1); 2686 $d += $o->{origin}; 2687 if($data->is_inplace) { 2688 $data .= $dd; 2689 return $data; 2690 } 2691 $dd; 2692 }; 2693 2694 $me->{inv} = sub { 2695 my($data,$opt) = @_; 2696 my($dd) = $data->copy ; 2697 my($d) = (defined $opt->{dim}) ? $dd->slice("0:".($opt->{dim}-1)) : $dd; 2698 my($o) = $opt->{origin}; 2699 my($s) = $opt->{str}; 2700 my($l) = $opt->{l0}; 2701 2702 $d .= ((-1 + sqrt(1 + 4 * $s/$l * abs($d-$o) * (1+abs($s)))) 2703 / 2 / $s * $l) * (1 - 2*($d < $o)); 2704 $d += $o; 2705 if($data->is_inplace) { 2706 $data .= $dd; 2707 return $data; 2708 } 2709 $dd; 2710 }; 2711 $me; 2712} 2713 2714 2715 2716 2717=head2 t_cubic 2718 2719=for usage 2720 2721 $t = t_cubic(<options>); 2722 2723=for ref 2724 2725Cubic scaling - cubic pincushion (n-d; with inverse) 2726 2727Cubic scaling is a generalization of t_quadratic to a purely 2728cubic expansion. 2729 2730The formula for the expansion is: 2731 2732 f(a) = ( a' + st * a'^3/L_0^2 ) / (1 + abs(st)) + origin 2733 2734where a'=(a-origin). That is a simple pincushion 2735expansion/contraction that is fixed at a distance of L_0 from the 2736origin. 2737 2738Because there is no quadratic term the result is always invertible with 2739one real root, and there is no mucking about with complex numbers or 2740multivalued solutions. 2741 2742 2743OPTIONS 2744 2745=over 3 2746 2747=item o,origin,Origin 2748 2749The origin of the pincushion. (default is the, er, origin). 2750 2751=item l,l0,length,Length,r0 2752 2753The fundamental scale of the transformation -- the radius that remains 2754unchanged. (default=1) 2755 2756 2757=item d, dim, dims, Dims 2758 2759The number of dimensions to treat (default is the 2760dimensionality of your input vectors) 2761 2762=back 2763 2764=cut 2765 2766sub t_cubic { 2767 my ($class) = 'PDL::Transform'; 2768 my($o) = $_[0]; 2769 if(ref $o ne 'HASH') { 2770 $o = {@_}; 2771 } 2772 my($me) = PDL::Transform::new($class); 2773 2774 $me->{params}->{dim} = _opt($o,['d','dim','dims','Dims'],undef); 2775 $me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0)); 2776 $me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1)); 2777 $me->{params}->{st} = _opt($o,['s','st','str'],pdl(0)); 2778 $me->{name} = "cubic"; 2779 2780 $me->{params}->{cuberoot} = sub { 2781 my $a = shift; 2782 my $as = 1 - 2*($a<0); 2783 return $as * ( abs($a) ** (1/3) ); 2784 }; 2785 2786 $me->{func} = sub { 2787 my($data,$o) = @_; 2788 my($dd) = $data->copy; 2789 my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd; 2790 2791 $d -= $o->{origin}; 2792 2793 my $dl0 = $d / $o->{l0}; 2794 2795 # f(x) = x + x**3 * ($o->{st} / $o->{l0}**2): 2796 # A = ($o->{st}/$dl0**2) 2797 # B = 0 2798 # C = 1 2799 # D = -f(x) 2800 $d += $o->{st} * $d * $dl0 * $dl0; 2801 $d /= ($o->{st}**2 + 1); 2802 2803 $d += $o->{origin}; 2804 2805 if($data->is_inplace) { 2806 $data .= $dd; 2807 return $data; 2808 } 2809 return $dd; 2810 }; 2811 2812 $me->{inv} = sub { 2813 my($data,$o) = @_; 2814 my($l) = $o->{l0}; 2815 2816 my($dd) = $data->copy; 2817 my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd; 2818 2819 $d -= $o->{origin}; 2820 $d *= ($o->{st}+1); 2821 2822 # Now we have to solve: 2823 # A x^3 + B X^2 + C x + D dd = 0 2824 # with the assignments above for A,B,C,D. 2825 # Since B is zero, this is greatly simplified - the discriminant is always negative, 2826 # so there is always exactly one real root. 2827 # 2828 # The only real hassle is creating a symmetric cube root; for convenience 2829 # is stashed in the params hash. 2830 2831 # First: create coefficients for mnemonics. 2832 my ($A, $C, $D) = ( $o->{st} / $l / $l, 1, - $d ); 2833 my $alpha = 27 * $A * $A * $D; 2834 my $beta = 3 * $A * $C; 2835 2836 my $inner_root = sqrt( $alpha * $alpha + 4 * $beta * $beta * $beta ); 2837 2838 $d .= (-1 / (3 * $A)) * 2839 ( 2840 + &{$o->{cuberoot}}( 0.5 * ( $alpha + $inner_root ) ) 2841 + &{$o->{cuberoot}}( 0.5 * ( $alpha - $inner_root ) ) 2842 ); 2843 2844 $d += $origin; 2845 2846 if($data->is_inplace) { 2847 $data .= $dd; 2848 return $data; 2849 } else { 2850 return $dd; 2851 } 2852 }; 2853 2854 $me; 2855} 2856 2857 2858 2859 2860=head2 t_quartic 2861 2862=for usage 2863 2864 $t = t_quartic(<options>); 2865 2866=for ref 2867 2868Quartic scaling -- cylindrical pincushion (n-d; with inverse) 2869 2870Quartic scaling is a generalization of t_quadratic to a quartic 2871expansion. Only even powers of the input coordinates are retained, 2872and (as with t_quadratic) sign is preserved, making it an odd function 2873although a true quartic transformation would be an even function. 2874 2875You can apply separate distortion along any of the principal axes. If 2876you want different axes, use L<t_wrap|/t_wrap> and 2877L<t_linear|/t_linear> to rotate them to the correct angle. The 2878scaling options may be scalars or vectors; if they are scalars then 2879the expansion is isotropic. 2880 2881The formula for the expansion is: 2882 2883 f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1) 2884 2885where <strength> is a scaling coefficient and <L_0> is a fundamental 2886length scale. Negative values of <strength> result in a pincushion 2887contraction. 2888 2889Note that, because quadratic scaling does not have a strict inverse for 2890coordinate systems that cross the origin, we cheat slightly and use 2891$x * abs($x) rather than $x**2. This does the Right thing for pincushion 2892and barrel distortion, but means that t_quadratic does not behave exactly 2893like t_cubic with a null cubic strength coefficient. 2894 2895OPTIONS 2896 2897=over 3 2898 2899=item o,origin,Origin 2900 2901The origin of the pincushion. (default is the, er, origin). 2902 2903=item l,l0,length,Length,r0 2904 2905The fundamental scale of the transformation -- the radius that remains 2906unchanged. (default=1) 2907 2908=item s,str,strength,Strength 2909 2910The relative strength of the pincushion. (default = 0.1) 2911 2912=item honest (default=0) 2913 2914Sets whether this is a true quadratic coordinate transform. The more 2915common form is pincushion or cylindrical distortion, which switches 2916branches of the square root at the origin (for symmetric expansion). 2917Setting honest to a non-false value forces true quadratic behavior, 2918which is not mirror-symmetric about the origin. 2919 2920=item d, dim, dims, Dims 2921 2922The number of dimensions to quadratically scale (default is the 2923dimensionality of your input vectors) 2924 2925 2926=back 2927 2928=cut 2929 2930sub t_quartic { 2931 my($class) = 'PDL::Transform'; 2932 my($o) = $_[0]; 2933 if(ref $o ne 'HASH') { 2934 $o = {@_}; 2935 } 2936 my($me) = PDL::Transform::new($class); 2937 2938 $me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0)); 2939 $me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1)); 2940 $me->{params}->{str} = _opt($o,['s','str','strength','Strength'],pdl(0.1)); 2941 $me->{params}->{dim} = _opt($o,['d','dim','dims','Dims']); 2942 $me->{name} = "quadratic"; 2943 2944 $me->{func} = sub { 2945 my($data,$o) = @_; 2946 my($dd) = $data->copy - $o->{origin}; 2947 my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd; 2948 $d += $o->{str} * ($d * abs($d)) / $o->{l0}; 2949 $d /= (abs($o->{str}) + 1); 2950 $d += $o->{origin}; 2951 if($data->is_inplace) { 2952 $data .= $dd; 2953 return $data; 2954 } 2955 $dd; 2956 }; 2957 2958 $me->{inv} = sub { 2959 my($data,$opt) = @_; 2960 my($dd) = $data->copy ; 2961 my($d) = (defined $opt->{dim}) ? $dd->slice("0:".($opt->{dim}-1)) : $dd; 2962 my($o) = $opt->{origin}; 2963 my($s) = $opt->{str}; 2964 my($l) = $opt->{l0}; 2965 2966 $d .= ((-1 + sqrt(1 + 4 * $s/$l * abs($d-$o) * (1+abs($s)))) 2967 / 2 / $s * $l) * (1 - 2*($d < $o)); 2968 $d += $o; 2969 if($data->is_inplace) { 2970 $data .= $dd; 2971 return $data; 2972 } 2973 $dd; 2974 }; 2975 $me; 2976} 2977 2978 2979 2980 2981=head2 t_spherical 2982 2983=for usage 2984 2985 $t = t_spherical(<options>); 2986 2987=for ref 2988 2989Convert Cartesian to spherical coordinates. (3-D; with inverse) 2990 2991Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta 2992is longitude, centered on 0, and phi is latitude, also centered on 0. 2993Unless you specify Euler angles, the pole points in the +Z direction 2994and the prime meridian is in the +X direction. The default is for 2995theta and phi to be in radians; you can select degrees if you want 2996them. 2997 2998Just as the L<t_radial|/t_radial> 2-D transform acts like a 3-D 2999cylindrical transform by ignoring third and higher dimensions, 3000Spherical acts like a hypercylindrical transform in four (or higher) 3001dimensions. Also as with L<t_radial|/t_radial>, you must manually specify 3002the origin if you want to use more dimensions than 3. 3003 3004To deal with latitude & longitude on the surface of a sphere (rather 3005than full 3-D coordinates), see 3006L<t_unit_sphere|PDL::Transform::Cartography/t_unit_sphere>. 3007 3008OPTIONS: 3009 3010=over 3 3011 3012=item o, origin, Origin [default (0,0,0)] 3013 3014This is the Cartesian origin of the spherical expansion. Pass in a PDL 3015or an array ref. 3016 3017=item e, euler, Euler [default (0,0,0)] 3018 3019This is a 3-vector containing Euler angles to change the angle of the 3020pole and ordinate. The first two numbers are the (theta, phi) angles 3021of the pole in a (+Z,+X) spherical expansion, and the last is the 3022angle that the new prime meridian makes with the meridian of a simply 3023tilted sphere. This is implemented by composing the output transform 3024with a PDL::Transform::Linear object. 3025 3026=item u, unit, Unit (default radians) 3027 3028This option sets the angular unit to be used. Acceptable values are 3029"degrees","radians", or reasonable substrings thereof (e.g. "deg", and 3030"rad", but "d" and "r" are deprecated). Once genuine unit processing 3031comes online (a la Math::Units) any angular unit should be OK. 3032 3033=back 3034 3035=cut 3036 3037sub t_spherical { 3038 my($class) = 'PDL::Transform'; 3039 my($o) = $_[0]; 3040 if(ref $o ne 'HASH') { 3041 $o = { @_ } ; 3042 } 3043 3044 my($me) = PDL::Transform::new($class); 3045 3046 $me->{idim}=3; 3047 $me->{odim}=3; 3048 3049 $me->{params}->{origin} = _opt($o,['o','origin','Origin']); 3050 $me->{params}->{origin} = PDL->zeroes(3) 3051 unless defined($me->{params}->{origin}); 3052 $me->{params}->{origin} = PDL->pdl($me->{params}->{origin}); 3053 3054 $me->{params}->{deg} = _opt($o,['d','degrees','Degrees']); 3055 3056 my $unit = _opt($o,['u','unit','Unit']); 3057 $me->{params}->{angunit} = ($unit =~ m/^d/i) ? 3058 $DEG2RAD : undef; 3059 3060 $me->{name} = "spherical"; 3061 3062 $me->{func} = sub { 3063 my($data,$o) = @_; 3064 my($d) = $data->copy - $o->{origin}; 3065 3066 my($d0,$d1,$d2) = ($d->slice("(0)"),$d->slice("(1)"),$d->slice("(2)")); 3067 my($out) = ($d->is_inplace) ? $data : $data->copy; 3068 3069 my $tmp; # work around perl -d "feature" 3070 ($tmp = $out->slice("(0)")) .= atan2($d1, $d0); 3071 ($tmp = $out->slice("(2)")) .= sqrt($d0*$d0 + $d1*$d1 + $d2*$d2); 3072 ($tmp = $out->slice("(1)")) .= asin($d2 / $out->slice("(2)")); 3073 3074 ($tmp = $out->slice("0:1")) /= $o->{angunit} 3075 if(defined $o->{angunit}); 3076 3077 $out; 3078 }; 3079 3080 $me->{inv} = sub { 3081 my($d,$o) = @_; 3082 3083 my($theta,$phi,$r,$out) = 3084 ( ($d->is_inplace) ? 3085 ($d->slice("(0)")->copy, $d->slice("(1)")->copy, $d->slice("(2)")->copy, $d) : 3086 ($d->slice("(0)"), $d->slice("(1)"), $d->slice("(2)"), $d->copy) 3087 ); 3088 3089 3090 my($x,$y,$z) = 3091 ($out->slice("(0)"),$out->slice("(1)"),$out->slice("(2)")); 3092 3093 my($ph,$th); 3094 if(defined $o->{angunit}){ 3095 $ph = $o->{angunit} * $phi; 3096 $th = $o->{angunit} * $theta; 3097 } else { 3098 $ph = $phi; 3099 $th = $theta; 3100 } 3101 3102 $z .= $r * sin($ph); 3103 $x .= $r * cos($ph); 3104 $y .= $x * sin($th); 3105 $x *= cos($th); 3106 $out += $o->{origin}; 3107 3108 $out; 3109 }; 3110 3111 $me; 3112 } 3113 3114 3115 3116 3117=head2 t_projective 3118 3119=for usage 3120 3121 $t = t_projective(<options>); 3122 3123=for ref 3124 3125Projective transformation 3126 3127Projective transforms are simple quadratic, quasi-linear 3128transformations. They are the simplest transformation that 3129can continuously warp an image plane so that four arbitrarily chosen 3130points exactly map to four other arbitrarily chosen points. They 3131have the property that straight lines remain straight after transformation. 3132 3133You can specify your projective transformation directly in homogeneous 3134coordinates, or (in 2 dimensions only) as a set of four unique points that 3135are mapped one to the other by the transformation. 3136 3137Projective transforms are quasi-linear because they are most easily 3138described as a linear transformation in homogeneous coordinates 3139(e.g. (x',y',w) where w is a normalization factor: x = x'/w, etc.). 3140In those coordinates, an N-D projective transformation is represented 3141as simple multiplication of an N+1-vector by an N+1 x N+1 matrix, 3142whose lower-right corner value is 1. 3143 3144If the bottom row of the matrix consists of all zeroes, then the 3145transformation reduces to a linear affine transformation (as in 3146L<t_linear|/t_linear>). 3147 3148If the bottom row of the matrix contains nonzero elements, then the 3149transformed x,y,z,etc. coordinates are related to the original coordinates 3150by a quadratic polynomial, because the normalization factor 'w' allows 3151a second factor of x,y, and/or z to enter the equations. 3152 3153OPTIONS: 3154 3155=over 3 3156 3157=item m, mat, matrix, Matrix 3158 3159If specified, this is the homogeneous-coordinate matrix to use. It must 3160be N+1 x N+1, for an N-dimensional transformation. 3161 3162=item p, point, points, Points 3163 3164If specified, this is the set of four points that should be mapped one to the other. 3165The homogeneous-coordinate matrix is calculated from them. You should feed in a 31662x2x4 PDL, where the 0 dimension runs over coordinate, the 1 dimension runs between input 3167and output, and the 2 dimension runs over point. For example, specifying 3168 3169 p=> pdl([ [[0,1],[0,1]], [[5,9],[5,8]], [[9,4],[9,3]], [[0,0],[0,0]] ]) 3170 3171maps the origin and the point (0,1) to themselves, the point (5,9) to (5,8), and 3172the point (9,4) to (9,3). 3173 3174This is similar to the behavior of fitwarp2d with a quadratic polynomial. 3175 3176=back 3177 3178=cut 3179 3180sub t_projective { 3181 my($class) = 'PDL::Transform'; 3182 my($o) = $_[0]; 3183 if(ref $o ne 'HASH') { 3184 $o = { @_ }; 3185 } 3186 3187 my($me) = PDL::Transform::new($class); 3188 3189 $me->{name} = "projective"; 3190 3191 ### Set options... 3192 3193 3194 $me->{params}->{idim} = $me->{idim} = _opt($o,['d','dim','Dim']); 3195 3196 my $matrix; 3197 if(defined ($matrix=_opt($o,['m','matrix','Matrix']))) { 3198 $matrix = pdl($matrix); 3199 die "t_projective: needs a square matrix" 3200 if($matrix->dims != 2 || $matrix->dim(0) != $matrix->dim(1)); 3201 3202 $me->{params}->{idim} = $matrix->dim(0)-1 3203 unless(defined($me->{params}->{idim})); 3204 3205 $me->{idim} = $me->{params}->{idim}; 3206 3207 die "t_projective: matrix not compatible with given dimension (should be N+1xN+1)\n" 3208 unless($me->{params}->{idim}==$matrix->dim(0)-1); 3209 3210 my $inv = $matrix->inv; 3211 print STDERR "t_projective: warning - received non-invertible matrix\n" 3212 unless(all($inv*0 == 0)); 3213 3214 $me->{params}->{matrix} = $matrix; 3215 $me->{params}->{matinv} = $inv; 3216 3217 } elsif(defined ($p=pdl(_opt($o,['p','point','points','Point','Points'])))) { 3218 die "t_projective: points array should be 2(x,y) x 2(in,out) x 4(point)\n\t(only 2-D points spec is available just now, sorry)\n" 3219 unless($p->dims==3 && all(pdl($p->dims)==pdl(2,2,4))); 3220 3221 # Solve the system of N equations to find the projective transform 3222 my ($p0,$p1,$p2,$p3) = ( $p->slice(":,(0),(0)"), $p->slice(":,(0),(1)"), $p->slice(":,(0),(2)"), $p->slice(":,(0),(3)") ); 3223 my ($P0,$P1,$P2,$P3) = ( $p->slice(":,(1),(0)"), $p->slice(":,(1),(1)"), $p->slice(":,(1),(2)"), $p->slice(":,(1),(3)") ); 3224#print "declaring PDL...\n" ; 3225 my $M = pdl( [ [$p0->at(0), $p0->at(1), 1, 0, 0, 0, -$P0->at(0)*$p0->at(0), -$P0->at(0)*$p0->at(1)], 3226 [ 0, 0, 0, $p0->at(0), $p0->at(1), 1, -$P0->at(1)*$p0->at(0), -$P0->at(1)*$p0->at(1)], 3227 [$p1->at(0), $p1->at(1), 1, 0, 0, 0, -$P1->at(0)*$p1->at(0), -$P1->at(0)*$p1->at(1)], 3228 [ 0, 0, 0, $p1->at(0), $p1->at(1), 1, -$P1->at(1)*$p1->at(0), -$P1->at(1)*$p1->at(1)], 3229 [$p2->at(0), $p2->at(1), 1, 0, 0, 0, -$P2->at(0)*$p2->at(0), -$P2->at(0)*$p2->at(1)], 3230 [ 0, 0, 0, $p2->at(0), $p2->at(1), 1, -$P2->at(1)*$p2->at(0), -$P2->at(1)*$p2->at(1)], 3231 [$p3->at(0), $p3->at(1), 1, 0, 0, 0, -$P3->at(0)*$p3->at(0), -$P3->at(0)*$p3->at(1)], 3232 [ 0, 0, 0, $p3->at(0), $p3->at(1), 1, -$P3->at(1)*$p3->at(0), -$P3->at(1)*$p3->at(1)] 3233 ] ); 3234#print "ok. Finding inverse...\n"; 3235 my $h = ($M->inv x $p->slice(":,(1),:")->flat->slice("*1"))->slice("(0)"); 3236# print "ok\n"; 3237 my $matrix = ones(3,3); 3238 my $tmp; # work around perl -d "feature" 3239 ($tmp = $matrix->flat->slice("0:7")) .= $h; 3240 $me->{params}->{matrix} = $matrix; 3241 3242 $me->{params}->{matinv} = $matrix->inv; 3243 } 3244 3245 3246 $me->{params}->{idim} = 2 unless defined $me->{params}->{idim}; 3247 $me->{params}->{odim} = $me->{params}->{idim}; 3248 $me->{idim} = $me->{params}->{idim}; 3249 $me->{odim} = $me->{params}->{odim}; 3250 3251 $me->{func} = sub { 3252 my($data,$o) = @_; 3253 my($id) = $data->dim(0); 3254 my($d) = $data->glue(0,ones($data->slice("0"))); 3255 my($out) = ($o->{matrix} x $d->slice("*1"))->slice("(0)"); 3256 return ($out->slice("0:".($id-1))/$out->slice("$id")); 3257 }; 3258 3259 $me->{inv} = sub { 3260 my($data,$o) = @_; 3261 my($id) = $data->dim(0); 3262 my($d) = $data->glue(0,ones($data->slice("0"))); 3263 my($out) = ($o->{matinv} x $d->slice("*1"))->slice("(0)"); 3264 return ($out->slice("0:".($id-1))/$out->slice("$id")); 3265 }; 3266 3267 $me; 3268} 3269 3270 3271 3272; 3273 3274 3275=head1 AUTHOR 3276 3277Copyright 2002, 2003 Craig DeForest. There is no warranty. You are allowed 3278to redistribute this software under certain conditions. For details, 3279see the file COPYING in the PDL distribution. If this file is 3280separated from the PDL distribution, the copyright notice should be 3281included in the file. 3282 3283=cut 3284 3285package PDL::Transform; 3286use Carp; 3287use overload '""' => \&_strval; 3288use overload 'x' => \&_compose_op; 3289use overload '**' => \&_pow_op; 3290use overload '!' => \&t_inverse; 3291 3292use PDL; 3293use PDL::MatrixOps; 3294 3295our $PI = 3.1415926535897932384626; 3296our $DEG2RAD = $PI / 180; 3297our $RAD2DEG = 180/$PI; 3298our $E = exp(1); 3299 3300 3301#### little helper kludge parses a list of synonyms 3302sub _opt { 3303 my($hash) = shift; 3304 my($synonyms) = shift; 3305 my($alt) = shift; # default is undef -- ok. 3306 local($_); 3307 foreach $_(@$synonyms){ 3308 return (UNIVERSAL::isa($alt,'PDL')) ? PDL->pdl($hash->{$_}) : $hash->{$_} 3309 if defined($hash->{$_}) ; 3310 } 3311 return $alt; 3312} 3313 3314###################################################################### 3315# 3316# Stringification hack. _strval just does a method search on stringify 3317# for the object itself. This gets around the fact that stringification 3318# overload is a subroutine call, not a method search. 3319# 3320 3321sub _strval { 3322 my($me) = shift; 3323 $me->stringify(); 3324} 3325 3326###################################################################### 3327# 3328# PDL::Transform overall stringifier. Subclassed stringifiers should 3329# call this routine first then append auxiliary information. 3330# 3331sub stringify { 3332 my($me) = shift; 3333 my($mestr) = (ref $me); 3334 $mestr =~ s/PDL::Transform:://; 3335 my $out = $mestr . " (" . $me->{name} . "): "; 3336 $out .= "fwd ". ((defined ($me->{func})) ? ( (ref($me->{func}) eq 'CODE') ? "ok" : "non-CODE(!!)" ): "missing")."; "; 3337 $out .= "inv ". ((defined ($me->{inv})) ? ( (ref($me->{inv}) eq 'CODE') ? "ok" : "non-CODE(!!)" ):"missing").".\n"; 3338} 3339 3340 3341 3342 3343 3344# Exit with OK status 3345 33461; 3347 3348