File Coverage

blib/lib/PDL/Image2D.pm
Criterion Covered Total %
statement 187 196 95.4
branch 59 90 65.5
condition 30 75 40.0
subroutine 23 24 95.8
pod 0 12 0.0
total 299 397 75.3


line stmt bran cond sub pod time code
1              
2             #
3             # GENERATED WITH PDL::PP! Don't modify!
4             #
5             package PDL::Image2D;
6              
7             @EXPORT_OK = qw( PDL::PP conv2d PDL::PP med2d PDL::PP med2df PDL::PP box2d PDL::PP patch2d PDL::PP patchbad2d PDL::PP max2d_ind PDL::PP centroid2d cc8compt cc4compt PDL::PP ccNcompt polyfill pnpoly polyfillv rotnewsz PDL::PP rot2d PDL::PP bilin2d PDL::PP rescale2d fitwarp2d applywarp2d PDL::PP warp2d warp2d_kernel PDL::PP warp2d_kernel );
8             %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
9              
10 4     4   2351 use PDL::Core;
  4         8  
  4         23  
11 4     4   31 use PDL::Exporter;
  4         9  
  4         19  
12 4     4   22 use DynaLoader;
  4         8  
  4         224  
13              
14              
15              
16            
17             @ISA = ( 'PDL::Exporter','DynaLoader' );
18             push @PDL::Core::PP, __PACKAGE__;
19             bootstrap PDL::Image2D ;
20              
21              
22              
23              
24              
25             =head1 NAME
26              
27             PDL::Image2D - Miscellaneous 2D image processing functions
28              
29             =head1 DESCRIPTION
30              
31             Miscellaneous 2D image processing functions - for want
32             of anywhere else to put them.
33              
34             =head1 SYNOPSIS
35              
36             use PDL::Image2D;
37              
38             =cut
39              
40 4     4   24 use PDL; # ensure qsort routine available
  4         8  
  4         29  
41 4     4   24 use PDL::Math;
  4         8  
  4         21  
42 4     4   26 use Carp;
  4         8  
  4         217  
43              
44 4     4   25 use strict;
  4         8  
  4         3244  
45              
46              
47              
48              
49              
50              
51              
52             =head1 FUNCTIONS
53              
54              
55              
56             =cut
57              
58              
59              
60              
61              
62              
63              
64              
65              
66              
67              
68              
69              
70              
71              
72              
73              
74              
75              
76              
77              
78             =head2 conv2d
79              
80             =for sig
81              
82             Signature: (a(m,n); kern(p,q); [o]b(m,n); int opt)
83              
84              
85             =for ref
86              
87             2D convolution of an array with a kernel (smoothing)
88              
89             For large kernels, using a FFT routine,
90             such as L in C,
91             will be quicker.
92              
93             =for usage
94              
95             $new = conv2d $old, $kernel, {OPTIONS}
96              
97             =for example
98              
99             $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}
100              
101             =for options
102              
103             Boundary - controls what values are assumed for the image when kernel
104             crosses its edge:
105             => Default - periodic boundary conditions
106             (i.e. wrap around axis)
107             => Reflect - reflect at boundary
108             => Truncate - truncate at boundary
109             => Replicate - repeat boundary pixel values
110              
111              
112              
113             =for bad
114              
115             Unlike the FFT routines, conv2d is able to process bad values.
116              
117             =cut
118              
119              
120              
121              
122              
123              
124             sub PDL::conv2d {
125 8 100   8 0 51 my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  8         33  
126 8 50 33     47 die 'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
127             if $#_<1 || $#_>2;
128 8         22 my($x,$kern) = @_;
129 8 50       38 my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
130             &PDL::_conv2d_int($x,$kern,$c,
131             (!(defined $opt && exists $$opt{Boundary}))?0:
132             (($$opt{Boundary} eq "Reflect") +
133             2*($$opt{Boundary} eq "Truncate") +
134 8 100 66     102216 3*($$opt{Boundary} eq "Replicate")));
135 8         112 return $c;
136             }
137              
138              
139              
140             *conv2d = \&PDL::conv2d;
141              
142              
143              
144              
145              
146             =head2 med2d
147              
148             =for sig
149              
150             Signature: (a(m,n); kern(p,q); [o]b(m,n); int opt)
151              
152              
153             =for ref
154              
155             2D median-convolution of an array with a kernel (smoothing)
156              
157             Note: only points in the kernel E0 are included in the median, other
158             points are weighted by the kernel value (medianing lots of zeroes
159             is rather pointless)
160              
161             =for usage
162              
163             $new = med2d $old, $kernel, {OPTIONS}
164              
165             =for example
166              
167             $smoothed = med2d $image, ones(3,3), {Boundary => Reflect}
168              
169             =for options
170              
171             Boundary - controls what values are assumed for the image when kernel
172             crosses its edge:
173             => Default - periodic boundary conditions (i.e. wrap around axis)
174             => Reflect - reflect at boundary
175             => Truncate - truncate at boundary
176             => Replicate - repeat boundary pixel values
177              
178              
179              
180             =for bad
181              
182             Bad values are ignored in the calculation. If all elements within the
183             kernel are bad, the output is set bad.
184              
185             =cut
186              
187              
188              
189              
190              
191              
192             sub PDL::med2d {
193 2 50   2 0 14 my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  2         10  
194 2 50 33     15 die 'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
195             if $#_<1 || $#_>2;
196 2         6 my($x,$kern) = @_;
197 2 50       74 croak "med2d: kernel must contain some positive elements.\n"
198             if all( $kern <= 0 );
199 2 50       19 my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
200             &PDL::_med2d_int($x,$kern,$c,
201             (!(defined $opt && exists $$opt{Boundary}))?0:
202             (($$opt{Boundary} eq "Reflect") +
203             2*($$opt{Boundary} eq "Truncate") +
204 2 50 33     169 3*($$opt{Boundary} eq "Replicate")));
205 2         27 return $c;
206             }
207              
208              
209              
210             *med2d = \&PDL::med2d;
211              
212              
213              
214              
215              
216             =head2 med2df
217              
218             =for sig
219              
220             Signature: (a(m,n); [o]b(m,n); int __p_size; int __q_size; int opt)
221              
222              
223             =for ref
224              
225             2D median-convolution of an array in a pxq window (smoothing)
226              
227             Note: this routine does the median over all points in a rectangular
228             window and is not quite as flexible as C in this regard
229             but slightly faster instead
230              
231             =for usage
232              
233             $new = med2df $old, $xwidth, $ywidth, {OPTIONS}
234              
235             =for example
236              
237             $smoothed = med2df $image, 3, 3, {Boundary => Reflect}
238              
239             =for options
240              
241             Boundary - controls what values are assumed for the image when kernel
242             crosses its edge:
243             => Default - periodic boundary conditions (i.e. wrap around axis)
244             => Reflect - reflect at boundary
245             => Truncate - truncate at boundary
246             => Replicate - repeat boundary pixel values
247              
248              
249              
250             =for bad
251              
252             med2df does not process bad values.
253             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
254              
255              
256             =cut
257              
258              
259              
260              
261              
262              
263             sub PDL::med2df {
264 1 50   1 0 9 my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  1         5  
265 1 50 33     7 die 'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )'
266             if $#_<2 || $#_>3;
267 1         3 my($x,$p,$q) = @_;
268 1 50 33     4 croak "med2df: kernel must contain some positive elements.\n"
269             if $p == 0 && $q == 0;
270 1 50       5 my $c = $#_ == 3 ? $_[3] : $x->nullcreate;
271             &PDL::_med2df_int($x,$c,$p,$q,
272             (!(defined $opt && exists $$opt{Boundary}))?0:
273             (($$opt{Boundary} eq "Reflect") +
274             2*($$opt{Boundary} eq "Truncate") +
275 1 50 33     62 3*($$opt{Boundary} eq "Replicate")));
276 1         7 return $c;
277             }
278              
279              
280              
281             *med2df = \&PDL::med2df;
282              
283              
284              
285              
286              
287             =head2 box2d
288              
289             =for sig
290              
291             Signature: (a(n,m); [o] b(n,m); int wx; int wy; int edgezero)
292              
293              
294             =for ref
295              
296             fast 2D boxcar average
297              
298             =for example
299              
300             $smoothim = $im->box2d($wx,$wy,$edgezero=1);
301              
302             The edgezero argument controls if edge is set to zero (edgezero=1)
303             or just keeps the original (unfiltered) values.
304              
305             C should be updated to support similar edge options
306             as C and C etc.
307              
308             Boxcar averaging is a pretty crude way of filtering. For serious stuff
309             better filters are around (e.g., use L with the appropriate
310             kernel). On the other hand it is fast and computational cost grows only
311             approximately linearly with window size.
312              
313              
314              
315             =for bad
316              
317             box2d does not process bad values.
318             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
319              
320              
321             =cut
322              
323              
324              
325              
326              
327              
328             *box2d = \&PDL::box2d;
329              
330              
331              
332              
333              
334             =head2 patch2d
335              
336             =for sig
337              
338             Signature: (a(m,n); int bad(m,n); [o]b(m,n))
339              
340              
341             =for ref
342              
343             patch bad pixels out of 2D images using a mask
344              
345             =for usage
346              
347             $patched = patch2d $data, $bad;
348              
349             C<$bad> is a 2D mask array where 1=bad pixel 0=good pixel.
350             Pixels are replaced by the average of their non-bad neighbours;
351             if all neighbours are bad, the original data value is
352             copied across.
353              
354              
355              
356             =for bad
357              
358             This routine does not handle bad values - use L instead
359              
360             =cut
361              
362              
363              
364              
365              
366              
367             *patch2d = \&PDL::patch2d;
368              
369              
370              
371              
372              
373             =head2 patchbad2d
374              
375             =for sig
376              
377             Signature: (a(m,n); [o]b(m,n))
378              
379              
380             =for ref
381              
382             patch bad pixels out of 2D images containing bad values
383              
384             =for usage
385              
386             $patched = patchbad2d $data;
387              
388             Pixels are replaced by the average of their non-bad neighbours;
389             if all neighbours are bad, the output is set bad.
390             If the input piddle contains I bad values, then a straight copy
391             is performed (see L).
392              
393              
394              
395             =for bad
396              
397             patchbad2d handles bad values. The output piddle I contain
398             bad values, depending on the pattern of bad values in the input piddle.
399              
400             =cut
401              
402              
403              
404              
405              
406              
407             *patchbad2d = \&PDL::patchbad2d;
408              
409              
410              
411              
412              
413             =head2 max2d_ind
414              
415             =for sig
416              
417             Signature: (a(m,n); [o]val(); int [o]x(); int[o]y())
418              
419              
420             =for ref
421              
422             Return value/position of maximum value in 2D image
423              
424             Contributed by Tim Jeness
425              
426              
427              
428             =for bad
429              
430             Bad values are excluded from the search. If all pixels
431             are bad then the output is set bad.
432              
433              
434              
435             =cut
436              
437              
438              
439              
440              
441              
442             *max2d_ind = \&PDL::max2d_ind;
443              
444              
445              
446              
447              
448             =head2 centroid2d
449              
450             =for sig
451              
452             Signature: (im(m,n); x(); y(); box(); [o]xcen(); [o]ycen())
453              
454              
455             =for ref
456              
457             Refine a list of object positions in 2D image by centroiding in a box
458              
459             C<$box> is the full-width of the box, i.e. the window
460             is C<+/- $box/2>.
461              
462              
463              
464             =for bad
465              
466             Bad pixels are excluded from the centroid calculation. If all elements are
467             bad (or the pixel sum is 0 - but why would you be centroiding
468             something with negatives in...) then the output values are set bad.
469              
470              
471              
472             =cut
473              
474              
475              
476              
477              
478              
479             *centroid2d = \&PDL::centroid2d;
480              
481              
482              
483              
484             =head2 cc8compt
485              
486             =for ref
487              
488             Connected 8-component labeling of a binary image.
489              
490             Connected 8-component labeling of 0,1 image - i.e. find separate
491             segmented objects and fill object pixels with object number.
492             8-component labeling includes all neighboring pixels.
493             This is just a front-end to ccNcompt. See also L.
494              
495             =for example
496              
497             $segmented = cc8compt( $image > $threshold );
498              
499             =head2 cc4compt
500              
501             =for ref
502              
503             Connected 4-component labeling of a binary image.
504              
505             Connected 4-component labeling of 0,1 image - i.e. find separate
506             segmented objects and fill object pixels with object number.
507             4-component labling does not include the diagonal neighbors.
508             This is just a front-end to ccNcompt. See also L.
509              
510             =for example
511              
512             $segmented = cc4compt( $image > $threshold );
513              
514             =cut
515              
516             sub PDL::cc8compt{
517 1     1 0 48 return ccNcompt(shift,8);
518             }
519             *cc8compt = \&PDL::cc8compt;
520              
521             sub PDL::cc4compt{
522 3     3 0 170 return ccNcompt(shift,4);
523             }
524             *cc4compt = \&PDL::cc4compt;
525              
526              
527              
528              
529              
530             =head2 ccNcompt
531              
532             =for sig
533              
534             Signature: (a(m,n); int+ [o]b(m,n); int con)
535              
536              
537              
538             =for ref
539              
540             Connected component labeling of a binary image.
541              
542             Connected component labeling of 0,1 image - i.e. find separate
543             segmented objects and fill object pixels with object number.
544             See also L and L.
545              
546             The connectivity parameter must be 4 or 8.
547              
548             =for example
549              
550             $segmented = ccNcompt( $image > $threshold, 4);
551              
552             $segmented2 = ccNcompt( $image > $threshold, 8);
553              
554             where the second parameter specifies the connectivity (4 or 8) of the labeling.
555              
556              
557              
558             =for bad
559              
560             ccNcompt ignores the bad-value flag of the input piddles.
561             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
562              
563              
564             =cut
565              
566              
567              
568              
569              
570              
571             *ccNcompt = \&PDL::ccNcompt;
572              
573              
574              
575             =head2 polyfill
576              
577             =for ref
578              
579             fill the area of the given polygon with the given colour.
580              
581             This function works inplace, i.e. modifies C.
582              
583             =for usage
584              
585             polyfill($im,$ps,$colour,[\%options]);
586              
587             The default method of determining which points lie inside of the polygon used
588             is not as strict as the method used in L. Often, it includes vertices
589             and edge points. Set the C option to change this behaviour.
590              
591             =for option
592              
593             Method - Set the method used to determine which points lie in the polygon.
594             => Default - internal PDL algorithm
595             => pnpoly - use the L algorithm
596              
597             =for example
598              
599             # Make a convex 3x3 square of 1s in an image using the pnpoly algorithm
600             $ps = pdl([3,3],[3,6],[6,6],[6,3]);
601             polyfill($im,$ps,1,{'Method' =>'pnpoly'});
602              
603             =cut
604             sub PDL::polyfill {
605 3     3 0 30 my $opt;
606 3 100       43 $opt = pop @_ if ref($_[-1]) eq 'HASH';
607 3 50       11 barf('Usage: polyfill($im,$ps,$colour,[\%options])') unless @_==3;
608 3         7 my ($im,$ps,$colour) = @_;
609              
610 3 100       10 if($opt) {
611 4     4   32 use PDL::Options qw();
  4         16  
  4         1707  
612 1         4 my $parsed = PDL::Options->new({'Method' => undef});
613 1         4 $parsed->options($opt);
614 1 50       2 if( $parsed->current->{'Method'} eq 'pnpoly' ) {
615 1         26 PDL::pnpolyfill_pp($im,$ps,$colour);
616             }
617             }
618             else
619             {
620 2         309 PDL::polyfill_pp($im,$ps,$colour);
621             }
622 2         28 return $im;
623              
624             }
625              
626             *polyfill = \&PDL::polyfill;
627              
628              
629              
630              
631             =head2 pnpoly
632              
633             =for ref
634              
635             'points in a polygon' selection from a 2-D piddle
636              
637             =for usage
638              
639             $mask = $img->pnpoly($ps);
640              
641             # Old style, do not use
642             $mask = pnpoly($x, $y, $px, $py);
643              
644             For a closed polygon determined by the sequence of points in {$px,$py}
645             the output of pnpoly is a mask corresponding to whether or not each
646             coordinate (x,y) in the set of test points, {$x,$y}, is in the interior
647             of the polygon. This is the 'points in a polygon' algorithm from
648             L
649             and vectorized for PDL by Karl Glazebrook.
650              
651             =for example
652              
653             # define a 3-sided polygon (a triangle)
654             $ps = pdl([3, 3], [20, 20], [34, 3]);
655              
656             # $tri is 0 everywhere except for points in polygon interior
657             $tri = $img->pnpoly($ps);
658              
659             With the second form, the x and y coordinates must also be specified.
660             B< I >.
661              
662             $px = pdl( 3, 20, 34 );
663             $py = pdl( 3, 20, 3 );
664             $x = $img->xvals; # get x pixel coords
665             $y = $img->yvals; # get y pixel coords
666              
667             # $tri is 0 everywhere except for points in polygon interior
668             $tri = pnpoly($x,$y,$px,$py);
669              
670             =cut
671              
672             # From: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
673             #
674             # Fixes needed to pnpoly code:
675             #
676             # Use topdl() to ensure piddle args
677             #
678             # Add POD docs for usage
679             #
680             # Calculate first term in & expression to use to fix divide-by-zero when
681             # the test point is in line with a vertical edge of the polygon.
682             # By adding the value of $mask we prevent a divide-by-zero since the &
683             # operator does not "short circuit".
684              
685             sub PDL::pnpoly {
686 2 50 66 2 0 17 barf('Usage: $mask = pnpoly($img, $ps);') unless(@_==2 || @_==4);
687 2         5 my ($tx, $ty, $vertx, $verty) = @_;
688              
689             # if only two inputs, use the pure PP version
690 2 100       5 unless( defined $vertx ) {
691 1 50       6 barf("ps must contain pairwise points.\n") unless $ty->getdim(0) == 2;
692              
693             # Input mapping: $img => $tx, $ps => $ty
694 1         87 return PDL::pnpoly_pp($tx,$ty);
695             }
696              
697 1         4 my $testx = PDL::Core::topdl($tx)->dummy(0);
698 1         4 my $testy = PDL::Core::topdl($ty)->dummy(0);
699 1         4 my $vertxj = PDL::Core::topdl($vertx)->rotate(1);
700 1         10 my $vertyj = PDL::Core::topdl($verty)->rotate(1);
701 1         64 my $mask = ( ($verty>$testy) == ($vertyj>$testy) );
702 1         9 my $c = sumover( ! $mask & ( $testx < ($vertxj-$vertx) * ($testy-$verty)
703             / ($vertyj-$verty+$mask) + $vertx) ) % 2;
704 1         23 return $c;
705             }
706              
707             *pnpoly = \&PDL::pnpoly;
708              
709              
710              
711              
712             =head2 polyfillv
713              
714             =for ref
715              
716             return the (dataflown) area of an image described by a polygon
717              
718             =for usage
719              
720             polyfillv($im,$ps,[\%options]);
721              
722             The default method of determining which points lie inside of the polygon used
723             is not as strict as the method used in L. Often, it includes vertices
724             and edge points. Set the C option to change this behaviour.
725              
726             =for option
727              
728             Method - Set the method used to determine which points lie in the polygon.
729             => Default - internal PDL algorithm
730             => pnpoly - use the L algorithm
731              
732             =for example
733              
734             # increment intensity in area bounded by $poly using the pnpoly algorithm
735             $im->polyfillv($poly,{'Method'=>'pnpoly'})++; # legal in perl >= 5.6
736              
737             # compute average intensity within area bounded by $poly using the default algorithm
738             $av = $im->polyfillv($poly)->avg;
739              
740             =cut
741              
742             sub PDL::polyfillv :lvalue {
743 2     2 0 12 my $opt;
744 2 100       9 $opt = pop @_ if ref($_[-1]) eq 'HASH';
745 2 50       6 barf('Usage: polyfillv($im,$ps,[\%options])') unless @_==2;
746              
747 2         4 my ($im,$ps) = @_;
748 2 50       8 barf("ps must contain pairwise points.\n") unless $ps->getdim(0) == 2;
749              
750 2 100       5 if($opt) {
751 4     4   33 use PDL::Options qw();
  4         11  
  4         6646  
752 1         8 my $parsed = PDL::Options->new({'Method' => undef});
753 1         4 $parsed->options($opt);
754 1 50       3 return $im->where(PDL::pnpoly_pp($im, $ps)) if $parsed->current->{'Method'} eq 'pnpoly';
755             }
756              
757 1         3 my $msk = zeroes(long,$im->dims);
758 1         21 PDL::polyfill_pp($msk, $ps, 1);
759 1         8 return $im->where($msk);
760             }
761             *polyfillv = \&PDL::polyfillv;
762              
763              
764              
765              
766              
767             =head2 rot2d
768              
769             =for sig
770              
771             Signature: (im(m,n); float angle(); bg(); int aa(); [o] om(p,q))
772              
773              
774             =for ref
775              
776             rotate an image by given C
777              
778             =for example
779              
780             # rotate by 10.5 degrees with antialiasing, set missing values to 7
781             $rot = $im->rot2d(10.5,7,1);
782              
783             This function rotates an image through an C between -90 and + 90
784             degrees. Uses/doesn't use antialiasing depending on the C flag.
785             Pixels outside the rotated image are set to C.
786              
787             Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based
788             on "A Fast Algorithm for General Raster Rotation" by Alan Paeth,
789             Graphics Interface '86, pp. 77-81.
790              
791             Use the C function to find out about the dimension of the
792             newly created image
793              
794             ($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;
795              
796             L offers a more general interface to
797             distortions, including rotation, with various types of sampling; but
798             rot2d is faster.
799              
800              
801              
802             =for bad
803              
804             rot2d ignores the bad-value flag of the input piddles.
805             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
806              
807              
808             =cut
809              
810              
811              
812              
813              
814              
815             *rot2d = \&PDL::rot2d;
816              
817              
818              
819              
820              
821             =head2 bilin2d
822              
823             =for sig
824              
825             Signature: (Int(n,m); O(q,p))
826              
827              
828             =for ref
829              
830             Bilinearly maps the first piddle in the second. The
831             interpolated values are actually added to the second
832             piddle which is supposed to be larger than the first one.
833              
834              
835              
836             =for bad
837              
838             bilin2d ignores the bad-value flag of the input piddles.
839             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
840              
841              
842             =cut
843              
844              
845              
846              
847              
848              
849             *bilin2d = \&PDL::bilin2d;
850              
851              
852              
853              
854              
855             =head2 rescale2d
856              
857             =for sig
858              
859             Signature: (Int(m,n); O(p,q))
860              
861              
862             =for ref
863              
864             The first piddle is rescaled to the dimensions of the second
865             (expanding or meaning values as needed) and then added to it in place.
866             Nothing useful is returned.
867              
868             If you want photometric accuracy or automatic FITS header metadata
869             tracking, consider using L
870             instead: it does these things, at some speed penalty compared to
871             rescale2d.
872              
873              
874              
875             =for bad
876              
877             rescale2d ignores the bad-value flag of the input piddles.
878             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
879              
880              
881             =cut
882              
883              
884              
885              
886              
887              
888             *rescale2d = \&PDL::rescale2d;
889              
890              
891              
892              
893              
894             =head2 fitwarp2d
895              
896             =for ref
897              
898             Find the best-fit 2D polynomial to describe
899             a coordinate transformation.
900              
901             =for usage
902              
903             ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf, { options } )
904              
905             Given a set of points in the output plane (C<$u,$v>), find
906             the best-fit (using singular-value decomposition) 2D polynomial
907             to describe the mapping back to the image plane (C<$x,$y>).
908             The order of the fit is controlled by the C<$nf> parameter
909             (the maximum power of the polynomial is C<$nf - 1>), and you
910             can restrict the terms to fit using the C option.
911              
912             C<$px> and C<$py> are C by C element piddles which describe
913             a polynomial mapping (of order C)
914             from the I C<(u,v)> image to the I C<(x,y)> image:
915              
916             x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
917             y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j
918              
919             The transformation is returned for the reverse direction (ie
920             output to input image) since that is what is required by the
921             L routine. The L
922             routine can be used to convert a set of C<$u,$v> points given
923             C<$px> and C<$py>.
924              
925             Options:
926              
927             =for options
928              
929             FIT - which terms to fit? default ones(byte,$nf,$nf)
930              
931             =begin comment
932              
933             old option, caused trouble
934             THRESH - in svd, remove terms smaller than THRESH * max value
935             default is 1.0e-5
936              
937             =end comment
938              
939             =over 4
940              
941             =item FIT
942              
943             C allows you to restrict which terms of the polynomial to fit:
944             only those terms for which the FIT piddle evaluates to true will be
945             evaluated. If a 2D piddle is sent in, then it is
946             used for the x and y polynomials; otherwise
947             C<< $fit->slice(":,:,(0)") >> will be used for C<$px> and
948             C<< $fit->slice(":,:,(1)") >> will be used for C<$py>.
949              
950             =begin comment
951              
952             =item THRESH
953              
954             Remove all singular values whose value is less than C
955             times the largest singular value.
956              
957             =end comment
958              
959             =back
960              
961             The number of points must be at least equal to the number of
962             terms to fit (C<$nf*$nf> points for the default value of C).
963              
964             =for example
965              
966             # points in original image
967             $x = pdl( 0, 0, 100, 100 );
968             $y = pdl( 0, 100, 100, 0 );
969             # get warped to these positions
970             $u = pdl( 10, 10, 90, 90 );
971             $v = pdl( 10, 90, 90, 10 );
972             #
973             # shift of origin + scale x/y axis only
974             $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
975             ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
976             print "px = ${px}py = $py";
977             px =
978             [
979             [-12.5 1.25]
980             [ 0 0]
981             ]
982             py =
983             [
984             [-12.5 0]
985             [ 1.25 0]
986             ]
987             #
988             # Compared to allowing all 4 terms
989             ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
990             print "px = ${px}py = $py";
991             px =
992             [
993             [ -12.5 1.25]
994             [ 1.110223e-16 -1.1275703e-17]
995             ]
996             py =
997             [
998             [ -12.5 1.6653345e-16]
999             [ 1.25 -5.8546917e-18]
1000             ]
1001              
1002             # A higher-degree polynomial should not affect the answer much, but
1003             # will require more control points
1004              
1005             $x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
1006             $y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
1007             $u = $u->glue(0,pdl(73,20,40,20,40));
1008             $v = $v->glue(0,pdl(29,20,40,40,20));
1009             ( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
1010             print "px3 =${px3}py3 =$py3";
1011             px3 =
1012             [
1013             [-6.4981162e+08 71034917 -726498.95]
1014             [ 49902244 -5415096.7 55945.388]
1015             [ -807778.46 88457.191 -902.51612]
1016             ]
1017             py3 =
1018             [
1019             [-6.2732159e+08 68576392 -701354.77]
1020             [ 48175125 -5227679.8 54009.114]
1021             [ -779821.18 85395.681 -871.27997]
1022             ]
1023              
1024             #This illustrates an important point about singular value
1025             #decompositions that are used in fitwarp2d: like all SVDs, the
1026             #rotation matrices are not unique, and so the $px and $py returned
1027             #by fitwarp2d are not guaranteed to be the "simplest" solution.
1028             #They do still work, though:
1029              
1030             ($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
1031             print approx $x3,$x,1e-4;
1032             [1 1 1 1 1 1 1 1 1]
1033             print approx $y3,$y;
1034             [1 1 1 1 1 1 1 1 1]
1035              
1036             =head2 applywarp2d
1037              
1038             =for ref
1039              
1040             Transform a set of points using a 2-D polynomial mapping
1041              
1042             =for usage
1043              
1044             ( $x, $y ) = applywarp2d( $px, $py, $u, $v )
1045              
1046             Convert a set of points (stored in 1D piddles C<$u,$v>)
1047             to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px>
1048             and C<$py>. See L
1049             for more information on the format of C<$px> and C<$py>.
1050              
1051             =cut
1052              
1053             # use SVD to fit data. Assuming no errors.
1054              
1055             =pod
1056              
1057             =begin comment
1058              
1059             Some explanation of the following three subroutines (_svd, _mkbasis,
1060             and fitwarp2d): See Wolberg 1990 (full ref elsewhere in this
1061             documentation), Chapter 3.6 "Polynomial Transformations". The idea is
1062             that, given a set of control points in the input and output images
1063             denoted by coordinates (x,y) and (u,v), we want to create a polynomial
1064             transformation that relates u to linear combinations of powers of x
1065             and y, and another that relates v to powers of x and y.
1066              
1067             The transformations used here and by Wolberg differ slightly, but the
1068             basic idea is something like this. For each u and each v, define a
1069             transform:
1070              
1071             u = (sum over j) (sum over i) a_{ij} x**i * y**j , (eqn 1)
1072             v = (sum over j) (sum over i) b_{ij} x**i * y**j . (eqn 2)
1073              
1074             We can write this in matrix notation. Given that there are M control
1075             points, U is a column vector with M rows. A and B are vectors containing
1076             the N coefficients (related to the degree of the polynomial fit). W
1077             is an MxN matrix of the basis row-vectors (the x**i y**j).
1078              
1079             The matrix equations we are trying to solve are
1080             U = W A , (eqn 3)
1081             V = W B . (eqn 4)
1082              
1083             We need to find the A and B column vectors, those are the coefficients
1084             of the polynomial terms in W. W is not square, so it has no inverse.
1085             But is has a pseudo-inverse W^+ that is NxM. We are going to use that
1086             pseudo-inverse to isolate A and B, like so:
1087              
1088             W^+ U = W^+ W A = A (eqn 5)
1089             W^+ V = W^+ W B = B (eqn 6)
1090              
1091             We are going to get W^+ by performing a singular value decomposition
1092             of W (to use some of the variable names below):
1093              
1094             W = $svd_u x SIGMA x $svd_v->transpose (eqn 7)
1095             W^+ = $svd_v x SIGMA^+ x $svd_u->transpose . (eqn 8)
1096              
1097             Here SIGMA is a square diagonal matrix that contains the singular
1098             values of W that are in the variable $svd_w. SIGMA^+ is the
1099             pseudo-inverse of SIGMA, which is calculated by replacing the non-zero
1100             singular values with their reciprocals, and then taking the transpose
1101             of the matrix (which is a no-op since the matrix is square and
1102             diagonal).
1103              
1104             So the code below does this:
1105              
1106             _mkbasis computes the matrix W, given control coordinates u and v and
1107             the maximum degree of the polynomial (and the terms to use).
1108              
1109             _svd takes the svd of W, computes the pseudo-inverse of W, and then
1110             multiplies that with the U vector (there called $y). The output of
1111             _svd is the A or B vector in eqns 5 & 6 above. Rarely is the matrix
1112             multiplication explicit, unfortunately, so I have added EXPLANATIONs
1113             below.
1114              
1115             =end comment
1116              
1117             =cut
1118              
1119             sub _svd ($$) {
1120 18     18   31 my $basis = shift;
1121 18         21 my $y = shift;
1122             # my $thresh = shift;
1123              
1124             # if we had errors for these points, would normalise the
1125             # basis functions, and the output array, by these errors here
1126              
1127             # perform the SVD
1128 18         2178 my ( $svd_u, $svd_w, $svd_v ) = svd( $basis );
1129              
1130             # DAL, 09/2017: the reason for removing ANY singular values, much less
1131             #the smallest ones, is not clear. I commented the line below since this
1132             #actually removes the LARGEST values in SIGMA^+.
1133             # $svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
1134             # The line below would instead remove the SMALLEST values in SIGMA^+, but I can see no reason to include it either.
1135             # $svd_w *= ( $svd_w <= ($svd_w->min / $thresh ) );
1136              
1137             # perform the back substitution
1138             # EXPLANATION: same thing as $svd_u->transpose x $y->transpose.
1139 18         66 my $tmp = $y x $svd_u;
1140              
1141             #EXPLANATION: the division by (the non-zero elements of) $svd_w (the singular values) is
1142             #equivalent to $sigma_plus x $tmp, so $tmp is now SIGMA^+ x $svd_u x $y
1143 18 50       38 if ( $PDL::Bad::Status ) {
1144 18         236 $tmp /= $svd_w->setvaltobad(0.0);
1145 18         83 $tmp->inplace->setbadtoval(0.0);
1146             } else {
1147             # not checked
1148 0         0 my $mask = ($svd_w == 0.0);
1149 0         0 $tmp /= ( $svd_w + $mask );
1150 0         0 $tmp *= ( 1 - $mask );
1151             }
1152              
1153             #EXPLANATION: $svd_v x SIGMA^+ x $svd_u x $y
1154 18         422 return sumover( $svd_v * $tmp );
1155              
1156             } # sub: _svd()
1157              
1158             #_mkbasis returns a piddle in which the k(=j*n+i)_th column is v**j * u**i
1159             #k=0 j=0 i=0
1160             #k=1 j=0 i=1
1161             #k=2 j=0 i=2
1162             #k=3 j=1 i=0
1163             # ...
1164              
1165             #each row corresponds to a control point
1166             #and the rows for each of the control points look like this, e.g.:
1167             # (1) (u) (u**2) (v) (vu) (v(u**2)) (v**2) ((v**2)u) ((v**2)(u**2))
1168             # and so on for the next control point.
1169              
1170             sub _mkbasis ($$$$) {
1171 14     14   21 my $fit = shift;
1172 14         16 my $npts = shift;
1173 14         19 my $u = shift;
1174 14         35 my $v = shift;
1175              
1176 14         43 my $n = $fit->getdim(0) - 1;
1177 14         31 my $ncoeff = sum( $fit );
1178              
1179 14         43 my $basis = zeroes( $u->type, $ncoeff, $npts );
1180 14         23 my $k = 0;
1181 14         32 foreach my $j ( 0 .. $n ) {
1182 39         546 my $tmp_v = $v**$j;
1183 39         189 foreach my $i ( 0 .. $n ) {
1184 117 100       250 if ( $fit->at($i,$j) ) {
1185 73         245 my $tmp = $basis->slice("($k),:");
1186 73         1456 $tmp .= $tmp_v * $u**$i;
1187 73         557 $k++;
1188             }
1189             }
1190             }
1191 14         30 return $basis;
1192              
1193             } # sub: _mkbasis()
1194              
1195             sub PDL::fitwarp2d {
1196 9 50 66 9 0 903 croak "Usage: (\$px,\$py) = fitwarp2d(x(m);y(m);u(m);v(m);\$nf; { options })"
      33        
1197             if $#_ < 4 or ( $#_ >= 5 and ref($_[5]) ne "HASH" );
1198              
1199 9         17 my $x = shift;
1200 9         16 my $y = shift;
1201 9         12 my $u = shift;
1202 9         12 my $v = shift;
1203 9         16 my $nf = shift;
1204              
1205 9         24 my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf) } ); #, THRESH => 1.0e-5 } );
1206 9 100       49 $opts->options( $_[0] ) if $#_ > -1;
1207 9         23 my $oref = $opts->current();
1208              
1209             # safety checks
1210 9         39 my $npts = $x->nelem;
1211 9 50 33     160 croak "fitwarp2d: x, y, u, and v must be the same size (and 1D)"
      33        
      33        
      33        
      33        
      33        
1212             unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem
1213             and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1;
1214              
1215             # my $svd_thresh = $$oref{THRESH};
1216             # croak "fitwarp2d: THRESH option must be >= 0."
1217             # if $svd_thresh < 0;
1218              
1219 9         19 my $fit = $$oref{FIT};
1220 9         18 my $fit_ndim = $fit->getndims();
1221 9 50 66     120 croak "fitwarp2d: FIT option must be sent a (\$nf,\$nf[,2]) element piddle"
      33        
      33        
      33        
1222             unless UNIVERSAL::isa($fit,"PDL") and
1223             ($fit_ndim == 2 or ($fit_ndim == 3 and $fit->getdim(2) == 2)) and
1224             $fit->getdim(0) == $nf and $fit->getdim(1) == $nf;
1225              
1226             # how many coeffs to fit (first we ensure $fit is either 0 or 1)
1227 9         188 $fit = convert( $fit != 0, byte );
1228              
1229 9         37 my ( $fitx, $fity, $ncoeffx, $ncoeffy, $ncoeff );
1230 9 100       22 if ( $fit_ndim == 2 ) {
1231 5         11 $fitx = $fit;
1232 5         5 $fity = $fit;
1233 5         16 $ncoeff = $ncoeffx = $ncoeffy = sum( $fit );
1234             } else {
1235 4         14 $fitx = $fit->slice(",,(0)");
1236 4         13 $fity = $fit->slice(",,(1)");
1237 4         17 $ncoeffx = sum($fitx);
1238 4         11 $ncoeffy = sum($fity);
1239 4 50       12 $ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy;
1240             }
1241              
1242 9 50       24 croak "fitwarp2d: number of points ($npts) must be >= \$ncoeff ($ncoeff)"
1243             unless $npts >= $ncoeff;
1244              
1245             # create the basis functions for the SVD fitting
1246 9         17 my ( $basisx, $basisy );
1247              
1248 9         21 $basisx = _mkbasis( $fitx, $npts, $u, $v );
1249              
1250 9 100       18 if ( $fit_ndim == 2 ) {
1251 5         13 $basisy = $basisx;
1252             } else {
1253 4         9 $basisy = _mkbasis( $fity, $npts, $u, $v );
1254             }
1255              
1256 9         28 my $px = _svd( $basisx, $x ); # $svd_thresh);
1257 9         35 my $py = _svd( $basisy, $y ); # $svd_thresh);
1258              
1259             # convert into $nf x $nf element piddles, if necessary
1260 9         31 my $nf2 = $nf * $nf;
1261              
1262 9 100 66     47 return ( $px->reshape($nf,$nf), $py->reshape($nf,$nf) )
1263             if $ncoeff == $nf2 and $ncoeffx == $ncoeffy;
1264              
1265             # re-create the matrix
1266 4         12 my $xtmp = zeroes( $nf, $nf );
1267 4         9 my $ytmp = zeroes( $nf, $nf );
1268              
1269 4         8 my $kx = 0;
1270 4         7 my $ky = 0;
1271 4         11 foreach my $i ( 0 .. ($nf - 1) ) {
1272 11         20 foreach my $j ( 0 .. ($nf - 1) ) {
1273 33 100       60 if ( $fitx->at($i,$j) ) {
1274 11         25 $xtmp->set($i,$j, $px->at($kx) );
1275 11         14 $kx++;
1276             }
1277 33 100       55 if ( $fity->at($i,$j) ) {
1278 11         22 $ytmp->set($i,$j, $py->at($ky) );
1279 11         16 $ky++;
1280             }
1281             }
1282             }
1283              
1284 4         62 return ( $xtmp, $ytmp )
1285              
1286             } # sub: fitwarp2d
1287              
1288             *fitwarp2d = \&PDL::fitwarp2d;
1289              
1290             sub PDL::applywarp2d {
1291             # checks
1292 1 50   1 0 8 croak "Usage: (\$x,\$y) = applywarp2d(px(nf,nf);py(nf,nf);u(m);v(m);)"
1293             if $#_ != 3;
1294              
1295 1         3 my $px = shift;
1296 1         3 my $py = shift;
1297 1         2 my $u = shift;
1298 1         3 my $v = shift;
1299 1         4 my $npts = $u->nelem;
1300              
1301             # safety check
1302 1 50 33     20 croak "applywarp2d: u and v must be the same size (and 1D)"
      33        
      33        
1303             unless $npts == $u->nelem and $npts == $v->nelem
1304             and $u->getndims == 1 and $v->getndims == 1;
1305              
1306 1         3 my $nf = $px->getdim(0);
1307 1         2 my $nf2 = $nf * $nf;
1308              
1309             # could remove terms with 0 coeff here
1310             # (would also have to remove them from px/py for
1311             # the matrix multiplication below)
1312             #
1313 1         5 my $mat = _mkbasis( ones(byte,$nf,$nf), $npts, $u, $v );
1314              
1315 1         7 my $x = reshape( $mat x $px->clump(-1)->transpose(), $npts );
1316 1         8 my $y = reshape( $mat x $py->clump(-1)->transpose(), $npts );
1317 1         8 return ( $x, $y );
1318              
1319             } # sub: applywarp2d
1320              
1321             *applywarp2d = \&PDL::applywarp2d;
1322              
1323              
1324              
1325              
1326             =head2 warp2d
1327              
1328             =for sig
1329              
1330             Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { options })
1331              
1332             =for ref
1333              
1334             Warp a 2D image given a polynomial describing the I mapping.
1335              
1336             =for usage
1337              
1338             $out = warp2d( $img, $px, $py, { options } );
1339              
1340             Apply the polynomial transformation encoded in the C<$px> and
1341             C<$py> piddles to warp the input image C<$img> into the output
1342             image C<$out>.
1343              
1344             The format for the polynomial transformation is described in
1345             the documentation for the L routine.
1346              
1347             At each point C, the closest 16 pixel values are combined
1348             with an interpolation kernel to calculate the value at C.
1349             The interpolation is therefore done in the image, rather than
1350             Fourier, domain.
1351             By default, a C kernel is used, but this can be changed
1352             using the C option discussed below
1353             (the choice of kernel depends on the frequency content of the input image).
1354              
1355             The routine is based on the C command from
1356             the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and
1357             for further details on image resampling see
1358             Wolberg, G., "Digital Image Warping", 1990, IEEE Computer
1359             Society Press ISBN 0-8186-8944-7).
1360              
1361             Currently the output image is the same size as the input one,
1362             which means data will be lost if the transformation reduces
1363             the pixel scale. This will (hopefully) be changed soon.
1364              
1365             =for example
1366              
1367             $img = rvals(byte,501,501);
1368             imag $img, { JUSTIFY => 1 };
1369             #
1370             # use a not-particularly-obvious transformation:
1371             # x = -10 + 0.5 * $u - 0.1 * $v
1372             # y = -20 + $v - 0.002 * $u * $v
1373             #
1374             $px = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
1375             $py = pdl( [ -20, 0 ], [ 1, 0.002 ] );
1376             $wrp = warp2d( $img, $px, $py );
1377             #
1378             # see the warped image
1379             imag $warp, { JUSTIFY => 1 };
1380              
1381             The options are:
1382              
1383             =for options
1384              
1385             KERNEL - default value is tanh
1386             NOVAL - default value is 0
1387              
1388             C is used to specify which interpolation kernel to use
1389             (to see what these kernels look like, use the
1390             L routine).
1391             The options are:
1392              
1393             =over 4
1394              
1395             =item tanh
1396              
1397             Hyperbolic tangent: the approximation of an ideal box filter by the
1398             product of symmetric tanh functions.
1399              
1400             =item sinc
1401              
1402             For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle,
1403             which produces a C interpolation kernel in the spatial domain:
1404              
1405             sinc(x) = sin(pi * x) / (pi * x)
1406              
1407             However, it is not ideal for the C<4x4> pixel region used here.
1408              
1409             =item sinc2
1410              
1411             This is the square of the sinc function.
1412              
1413             =item lanczos
1414              
1415             Although defined differently to the C kernel, the result is very
1416             similar in the spatial domain. The Lanczos function is defined as
1417              
1418             L(x) = sinc(x) * sinc(x/2) if abs(x) < 2
1419             = 0 otherwise
1420              
1421             =item hann
1422              
1423             This kernel is derived from the following function:
1424              
1425             H(x) = a + (1-a) * cos(2*pi*x/(N-1)) if abs(x) < 0.5*(N-1)
1426             = 0 otherwise
1427              
1428             with C and N currently equal to 2001.
1429              
1430             =item hamming
1431              
1432             This kernel uses the same C as the Hann filter, but with
1433             C.
1434              
1435             =back
1436              
1437             C gives the value used to indicate that a pixel in the
1438             output image does not map onto one in the input image.
1439              
1440             =cut
1441              
1442             # support routine
1443             {
1444             my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann );
1445              
1446             # note: convert to lower case
1447             sub _check_kernel ($$) {
1448 6     6   16 my $kernel = lc shift;
1449 6         10 my $code = shift;
1450             barf "Unknown kernel $kernel sent to $code\n" .
1451             "\tmust be one of [" . join(',',keys %warp2d) . "]\n"
1452 6 50       17 unless exists $warp2d{$kernel};
1453 6         19 return $kernel;
1454             }
1455             }
1456              
1457              
1458              
1459              
1460              
1461             sub PDL::warp2d {
1462 6     6 0 52 my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );
1463 6 50       20 $opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";
1464              
1465 6 50 33     28 die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} )"
1466             if $#_<2 || $#_>3;
1467 6         9 my $img = shift;
1468 6         8 my $px = shift;
1469 6         8 my $py = shift;
1470 6 50       22 my $out = $#_ == -1 ? PDL->null() : shift;
1471              
1472             # safety checks
1473 6         19 my $copt = $opts->current();
1474 6         16 my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );
1475              
1476 6         40476 &PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} );
1477 6         85 return $out;
1478             }
1479              
1480              
1481              
1482             *warp2d = \&PDL::warp2d;
1483              
1484              
1485              
1486              
1487              
1488             =head2 warp2d_kernel
1489              
1490             =for ref
1491              
1492             Return the specified kernel, as used by L
1493              
1494             =for usage
1495              
1496             ( $x, $k ) = warp2d_kernel( $name )
1497              
1498             The valid values for C<$name> are the same as the C option
1499             of L.
1500              
1501             =for example
1502              
1503             line warp2d_kernel( "hamming" );
1504              
1505             =cut
1506              
1507              
1508              
1509              
1510              
1511             sub PDL::warp2d_kernel ($) {
1512 0     0 0   my $kernel = _check_kernel( shift, "warp2d_kernel" );
1513              
1514 0           my $nelem = _get_kernel_size();
1515 0           my $x = zeroes( $nelem );
1516 0           my $k = zeroes( $nelem );
1517              
1518 0           &PDL::_warp2d_kernel_int( $x, $k, $kernel );
1519 0           return ( $x, $k );
1520              
1521             # return _get_kernel( $kernel );
1522             }
1523             *warp2d_kernel = \&PDL::warp2d_kernel;
1524              
1525              
1526              
1527             *warp2d_kernel = \&PDL::warp2d_kernel;
1528              
1529              
1530              
1531             ;
1532              
1533              
1534             =head1 AUTHORS
1535              
1536             Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams
1537             (rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu),
1538             and Doug Burke (burke@ifa.hawaii.edu).
1539              
1540             All rights reserved. There is no warranty. You are allowed
1541             to redistribute this software / documentation under certain
1542             conditions. For details, see the file COPYING in the PDL
1543             distribution. If this file is separated from the PDL distribution,
1544             the copyright notice should be included in the file.
1545              
1546             =cut
1547              
1548              
1549              
1550              
1551              
1552             # Exit with OK status
1553              
1554             1;
1555              
1556