4
pp_addpm({At=>'Top'},<<'EOD');
8
PDL::Image2D - Miscellaneous 2D image processing functions
12
Miscellaneous 2D image processing functions - for want
13
of anywhere else to put them
21
use PDL; # ensure qsort routine available
29
pp_addpm({At=>'Bot'},<<'EOD');
33
Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams
34
(rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu),
35
and Doug Burke (burke@ifa.hawaii.edu).
37
All rights reserved. There is no warranty. You are allowed
38
to redistribute this software / documentation under certain
39
conditions. For details, see the file COPYING in the PDL
40
distribution. If this file is separated from the PDL distribution,
41
the copyright notice should be included in the file.
50
#define IsNaN(x) (x != x)
52
/* Fast Modulus with proper negative behaviour */
54
#define REALMOD(a,b) {while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);}
56
/* rint is missing on some platforms (eg Win32) */
59
#define rint(X) floor( X + 0.5 )
64
for (keys %PDL::Types::typehash) {
65
my $ctype = $PDL::Types::typehash{$_}{ctype};
66
my $ppsym = $PDL::Types::typehash{$_}{ppsym};
71
* this routine is based on code referenced from
72
* http://www.eso.org/~ndevilla/median/
73
* the original algorithm is described in Numerical Recipes
76
#define ELEM_SWAP(a,b) { register $ctype t=(a);(a)=(b);(b)=t; }
78
$ctype quick_select_$ppsym($ctype arr[], int n)
84
low = 0 ; high = n-1 ; median = (low + high) / 2;
86
if (high <= low) /* One element only */
89
if (high == low + 1) { /* Two elements only */
90
if (arr[low] > arr[high])
91
ELEM_SWAP(arr[low], arr[high]) ;
95
/* Find median of low, middle and high items; swap into position low */
96
middle = (low + high) / 2;
97
if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
98
if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
99
if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
101
/* Swap low item (now in position middle) into position (low+1) */
102
ELEM_SWAP(arr[middle], arr[low+1]) ;
104
/* Nibble from each end towards middle, swapping items when stuck */
108
do ll++; while (arr[low] > arr[ll]) ;
109
do hh--; while (arr[hh] > arr[low]) ;
114
ELEM_SWAP(arr[ll], arr[hh]) ;
117
/* Swap middle item (in position low) back into correct position */
118
ELEM_SWAP(arr[low], arr[hh]) ;
120
/* Re-set active partition */
135
i => { size => 'm_size', off => 'poff', init => '1-p_size' },
136
j => { size => 'n_size', off => 'qoff', init => '1-q_size' },
139
# requires 'int $var, ${var}2' to have been declared in the c code
140
# (along with [pq]off and [pq]_size)
146
my $loop2 = "${var}2";
148
my $href = $init{$var} ||
149
die "ERROR: unknown variable sent to init_map()\n";
150
my $size = $href->{size} ||
151
die "ERROR: unable to find size for $var\n";
152
my $off = $href->{off} ||
153
die "ERROR: unable to find off for $var\n";
154
my $init = $href->{init} ||
155
die "ERROR: unable to find init for $var\n";
158
"for ( $loop = $init; $loop< $size; ${loop}++) {
159
$loop2 = $loop + $off;
161
case 1: /* REFLECT */
163
$loop2 = -${loop2}-1;
164
else if ($loop2 >= $size)
165
$loop2 = 2*${size}-(${loop2}+1);
167
case 2: /* TRUNCATE */
168
if (${loop2}<0 || ${loop2} >= $size)
172
REALMOD($loop2,$size);
174
map${var}\[$loop] = $loop2;
180
my $href = shift || { };
181
$href->{vars} = '' unless defined $href->{vars};
182
$href->{malloc} = '' unless defined $href->{malloc};
183
$href->{check} = '' unless defined $href->{check};
185
my $str = $href->{vars};
186
$str .= "int i,j, i1,j1, i2,j2, poff, qoff;";
188
'int opt = $COMP(opt);
189
int m_size = $COMP(__m_size);
190
int n_size = $COMP(__n_size);
191
int p_size = $COMP(__p_size);
192
int q_size = $COMP(__q_size);
195
mapi = (int *) malloc((p_size+m_size)*sizeof(int));
196
mapj = (int *) malloc((q_size+n_size)*sizeof(int));
198
$str .= $href->{malloc} . "\n";
199
$str .= "if ($href->{check} (mapi==NULL) || (mapj==NULL))\n";
200
$str .= ' barf("Out of Memory");
202
poff = p_size/2; mapi += p_size-1;
203
qoff = q_size/2; mapj += q_size-1;
209
pp_def('conv2d', Doc=><<'EOD',
212
2D convolution of an array with a kernel (smoothing)
214
For large kernels, using a FFT routine,
215
such as L<fftconvolve()|PDL::FFT/fftconvolve> in C<PDL::FFT>,
220
$new = conv2d $old, $kernel, {OPTIONS}
224
$smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}
228
Boundary - controls what values are assumed for the image when kernel
230
=> Default - periodic boundary conditions
231
(i.e. wrap around axis)
232
=> Reflect - reflect at boundary
233
=> Truncate - truncate at boundary
237
'Unlike the FFT routines, conv2d is able to process bad values.',
239
Pars => 'a(m,n); kern(p,q); [o]b(m,n);',
240
OtherPars => 'int opt;',
244
my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\';
245
die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\'
248
my $c = $#_ == 2 ? $_[2] : $a->nullcreate;
249
&PDL::_conv2d_int($a,$kern,$c,
250
(!(defined $opt && exists $$opt{Boundary}))?0:
251
(($$opt{Boundary} eq "Reflect") +
252
2*($$opt{Boundary} eq "Truncate")));
258
init_vars( { vars => 'PDL_Double tmp;' } ) .
263
for(j=0; j<n_size; j++) {
264
for(i=0; i<m_size; i++) {
266
for(j1=0; j1<q_size; j1++) {
270
for(i1=0; i1<p_size; i1++) {
273
tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1);
281
free(mapj+1-q_size); free(mapi+1-p_size);',
283
init_vars( { vars => 'PDL_Double tmp; int flag;' } ) .
288
for(j=0; j<n_size; j++) {
289
for(i=0; i<m_size; i++) {
291
for(j1=0; j1<q_size; j1++) {
295
for(i1=0; i1<p_size; i1++) {
298
if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(kern(p=>i1,q=>j1)) ) {
299
tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1);
306
if ( flag ) { $b(m=>i,n=>j) = tmp; }
307
else { $SETBAD(b(m=>i,n=>j)); }
311
free(mapj+1-q_size); free(mapi+1-p_size);',
315
pp_def('med2d', Doc=> <<'EOD',
318
2D median-convolution of an array with a kernel (smoothing)
320
Note: only points in the kernel E<gt>0 are included in the median, other
321
points are weighted by the kernel value (medianing lots of zeroes
326
$new = med2d $old, $kernel, {OPTIONS}
330
$smoothed = med2d $image, ones(3,3), {Boundary => Reflect}
334
Boundary - controls what values are assumed for the image when kernel
336
=> Default - periodic boundary conditions (i.e. wrap around axis)
337
=> Reflect - reflect at boundary
338
=> Truncate - truncate at boundary
342
'Bad values are ignored in the calculation. If all elements within the
343
kernel are bad, the output is set bad.',
345
Pars => 'a(m,n); kern(p,q); [o]b(m,n);',
346
OtherPars => 'int opt;',
350
my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\';
351
die \'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\'
354
croak "med2d: kernel must contain some positive elements.\n"
355
if all( $kern <= 0 );
356
my $c = $#_ == 2 ? $_[2] : $a->nullcreate;
357
&PDL::_med2d_int($a,$kern,$c,
358
(!(defined $opt && exists $opt->{Boundary}))?0:
359
(($$opt{Boundary} eq "Reflect") +
360
2*($$opt{Boundary} eq "Truncate")));
366
init_vars( { vars => 'PDL_Double *tmp, kk; int count;',
367
malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));',
368
check => '(tmp==NULL) || ' } ) .
373
for(j=0; j<n_size; j++) {
374
for(i=0; i<m_size; i++) {
376
for(j1=0; j1<q_size; j1++) {
380
for(i1=0; i1<p_size; i1++) {
383
kk = $kern(p=>i1,q=>j1);
385
tmp[count++] = $a(m=>i2,n=>j2) * kk;
391
PDL->qsort_D( tmp, 0, count-1 );
392
$b(m=>i,n=>j) = tmp[(count-1)/2];
397
free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
400
init_vars( { vars => 'PDL_Double *tmp, kk, aa; int count, flag;',
401
malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));',
402
check => '(tmp==NULL) || ' } ) .
407
for(j=0; j<n_size; j++) {
408
for(i=0; i<m_size; i++) {
411
for(j1=0; j1<q_size; j1++) {
415
for(i1=0; i1<p_size; i1++) {
418
kk = $kern(p=>i1,q=>j1);
419
aa = $a(m=>i2,n=>j2);
420
if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(aa,a) ) {
423
tmp[count++] = aa * kk;
430
$SETBAD(b(m=>i,n=>j));
433
PDL->qsort_D( tmp, 0, count-1 );
434
$b(m=>i,n=>j) = tmp[(count-1)/2];
440
free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
445
pp_def('med2df', Doc=> <<'EOD',
448
2D median-convolution of an array in a pxq window (smoothing)
450
Note: this routine does the median over all points in a rectangular
451
window and is not quite as flexible as C<med2d> in this regard
452
but slightly faster instead
456
$new = med2df $old, $xwidth, $ywidth, {OPTIONS}
460
$smoothed = med2df $image, 3, 3, {Boundary => Reflect}
464
Boundary - controls what values are assumed for the image when kernel
466
=> Default - periodic boundary conditions (i.e. wrap around axis)
467
=> Reflect - reflect at boundary
468
=> Truncate - truncate at boundary
471
Pars => 'a(m,n); [o]b(m,n);',
472
# funny parameter names to avoid special case in 'init_vars'
473
OtherPars => 'int __p_size; int __q_size; int opt;',
477
my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\';
478
die \'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )\'
481
croak "med2df: kernel must contain some positive elements.\n"
482
if $p == 0 && $q == 0;
483
my $c = $#_ == 3 ? $_[3] : $a->nullcreate;
484
&PDL::_med2df_int($a,$c,$p,$q,
485
(!(defined $opt && exists $opt->{Boundary}))?0:
486
(($$opt{Boundary} eq "Reflect") +
487
2*($$opt{Boundary} eq "Truncate")));
493
init_vars( { vars => '$GENERIC() *tmp, kk; int count;',
494
malloc => 'tmp = malloc(p_size*q_size*sizeof($GENERIC()));',
495
check => '(tmp==NULL) || ' } ) .
500
for(j=0; j<n_size; j++) {
501
for(i=0; i<m_size; i++) {
503
for(j1=0; j1<q_size; j1++) {
507
for(i1=0; i1<p_size; i1++) {
510
tmp[count++] = $a(m=>i2,n=>j2);
516
quick_select_$TBSULFD(B,S,U,L,F,D) (tmp, count );
521
free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
530
patch bad pixels out of 2D images using a mask
534
$patched = patch2d $data, $bad;
536
C<$bad> is a 2D mask array where 1=bad pixel 0=good pixel.
537
Pixels are replaced by the average of their non-bad neighbours;
538
if all neighbours are bad, the original data value is
543
'This routine does not handle bad values - use L<patchbad2d|/patchbad2d> instead',
545
Pars => 'a(m,n); int bad(m,n); [o]b(m,n);',
547
'int m_size, n_size, i,j, i1,j1, i2,j2, norm;
550
m_size = $COMP(__m_size); n_size = $COMP(__n_size);
554
for(j=0; j<n_size; j++) {
555
for(i=0; i<m_size; i++) {
557
$b(m=>i,n=>j) = $a(m=>i,n=>j);
559
if ( $bad(m=>i,n=>j)==1 ) {
561
for(j1=-1; j1<=1; j1++) {
563
if ( j2>=0 && j2<n_size ) {
564
for(i1=-1; i1<=1; i1++) {
565
/* ignore central pixel, which we know is bad */
566
if ( i1!=0 || j1!=0 ) {
568
if ( i2>=0 && i2<m_size && $bad(m=>i2,n=>j2)!=1 ) {
569
tmp += $a(m=>i2,n=>j2);
572
} /* if: i1!=0 || j1!=0 */
577
if (norm>0) { /* Patch */
578
$b(m=>i,n=>j) = tmp/norm;
595
patch bad pixels out of 2D images containing bad values
599
$patched = patchbad2d $data;
601
Pixels are replaced by the average of their non-bad neighbours;
602
if all neighbours are bad, the output is set bad.
603
If the input piddle contains I<no> bad values, then a straight copy
604
is performed (see L<patch2d|/patch2d>).
608
'patchbad2d handles bad values. The output piddle I<may> contain
609
bad values, depending on the pattern of bad values in the input piddle.',
611
Pars => 'a(m,n); [o]b(m,n);',
612
Code => 'loop(n,m) %{ $b() = $a(); %}', # just copy
613
CopyBadStatusCode => '', # handled by BadCode
615
'int m_size, n_size, i,j, i1,j1, i2,j2, norm, flag;
620
m_size = $COMP(__m_size); n_size = $COMP(__n_size);
624
for(j=0; j<n_size; j++) {
625
for(i=0; i<m_size; i++) {
627
a_val = $a(m=>i,n=>j);
628
if ( $ISGOODVAR(a_val,a) ) {
629
$b(m=>i,n=>j) = a_val;
633
for(j1=-1; j1<=1; j1++) {
635
if ( j2>=0 && j2<n_size ) {
636
for(i1=-1; i1<=1; i1++) {
637
/* ignore central pixel, which we know is bad */
638
if ( i1!=0 || j1!=0 ) {
640
if ( i2>=0 && i2<m_size ) {
641
a_val = $a(m=>i2,n=>j2);
642
if ( $ISGOODVAR(a_val,a) ) {
647
} /* if: i1!=0 || j1!=0 */
654
$b(m=>i,n=>j) = tmp/norm;
656
$SETBAD(b(m=>i,n=>j));
660
} /* if: ISGOODVAR() */
667
/* handle bad flag */
668
if ( flag ) $PDLSTATESETBAD(b);
676
Return value/position of maximum value in 2D image
678
Contributed by Tim Jeness
684
Bad values are excluded from the search. If all pixels
685
are bad then the output is set bad.
690
Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();',
692
double cur; int curind1; int curind2;
697
if((!m && !n) || $a() > cur || IsNaN(cur)) {
698
cur = $a(); curind1 = m; curind2 = n;
707
double cur; int curind1; int curind2;
712
if( $ISGOOD(a()) && ( (!n && !m) || ($a() > cur) ) ) {
713
cur = $a(); curind1 = m; curind2 = n;
733
Refine a list of object positions in 2D image by centroiding in a box
735
C<$box> is the full-width of the box, i.e. the window
741
Bad pixels are excluded from the centroid calculation. If all elements are
742
bad (or the pixel sum is 0 - but why would you be centroiding
743
something with negatives in...) then the output values are set bad.
748
Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();',
751
int i,j,i1,i2,j1,j2,m_size,n_size;
752
double sum,data,sumx,sumy;
754
m_size = $SIZE(m); n_size = $SIZE(n);
756
i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;
757
i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;
758
j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
759
j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;
761
sum = sumx = sumy = 0;
762
for(j=j1; j<=j2; j++) { for(i=i1; i<=i2; i++) {
763
data = $im(m=>i,n=>j);
773
int i,j,i1,i2,j1,j2,m_size,n_size;
774
double sum,data,sumx,sumy;
776
m_size = $SIZE(m); n_size = $SIZE(n);
778
i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;
779
i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;
780
j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
781
j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;
783
sum = sumx = sumy = 0;
784
for(j=j1; j<=j2; j++) {
785
for(i=i1; i<=i2; i++) {
786
data = $im(m=>i,n=>j);
787
if ( $ISGOODVAR(data,im) ) {
795
* if sum == 0 then we will flag as bad -- although it could just mean that
796
* there is negative values in the dataset.
797
* - should use a better check than != 0.0 ...
811
/* Add an equivalence to a list - used by pdl_cc8compt */
813
void AddEquiv ( PDL_Long* equiv, PDL_Long i, PDL_Long j) {
823
} while ( k != j && k != i );
835
pp_def('cc8compt',Doc=>'
838
Connected 8-component labeling of a binary image.
840
Connected 8-component labeling of 0,1 image - i.e. find seperate
841
segmented objects and fill object pixels with object number
845
$segmented = cc8compt( $image > $threshold );
848
HandleBad => 0, # a marker
849
Pars => 'a(m,n); [o]b(m,n);',
854
PDL_Long neighbour[4];
856
PDL_Long pass,count,next,this;
859
PDL_Long nx = $SIZE(m);
860
PDL_Long ny = $SIZE(n);
862
loop(n) %{ loop(m) %{ /* Copy */
866
/* 1st pass counts max possible compts, 2nd records equivalences */
868
for (pass = 0; pass<2; pass++) {
871
equiv = (PDL_Long*) malloc((newlabel+1)*sizeof(PDL_Long));
872
if (equiv==(PDL_Long*)0)
873
barf("Out of memory");
874
for(i=0;i<=newlabel;i++)
878
newlabel = 1; /* Running label */
880
for(j=0; j<ny; j++) { for(i=0; i<nx; i++) { /* Loop over image pixels */
882
nfound = 0; /* Number of neighbour >0 */
884
i1 = i-1; j1 = j-1; i2 = i+1;
886
if ($b(m=>i, n=>j) > 0) { /* Check 4 neighbour already seen */
888
if (i>0 && $b(m=>i1, n=>j)>0)
889
neighbour[nfound++] = $b(m=>i1, n=>j); /* Store label of it */
890
if (j>0 && $b(m=>i, n=>j1)>0)
891
neighbour[nfound++] = $b(m=>i, n=>j1);
892
if (j>0 && i>0 && $b(m=>i1, n=>j1)>0)
893
neighbour[nfound++] = $b(m=>i1, n=>j1);
894
if (j>0 && i<(nx-1) && $b(m=>i2, n=>j1)>0)
895
neighbour[nfound++] = $b(m=>i2, n=>j1);
897
if (nfound==0) { /* Assign new label */
898
$b(m=>i, n=>j) = newlabel++;
901
$b(m=>i, n=>j) = neighbour[0];
902
if (nfound>1 && pass == 1) { /* Assign equivalents */
903
for(k=1; k<nfound; k++)
904
AddEquiv( equiv, (PDL_Long)$b(m=>i, n=>j),
910
else { /* No label */
915
}} /* End of image loop */
919
/* Replace each cycle by single label */
922
for (i = 1; i <= newlabel; i++)
923
if ( i <= equiv[i] ) {
926
while ( equiv[this] != i ) {
935
/* Now remove equivalences */
937
for(j=0; j<ny; j++) { for(i=0; i<nx; i++) { /* Loop over image pixels */
938
$b(m=>i, n=>j) = equiv[ (PDL_Long) $b(m=>i, n=>j) ] ;
941
free(equiv); /* Tidy */
946
#define line(x1, x2, y) for (k=x1;k<=x2;k++) \
947
{ /* printf("line from %d to %d\n",x1,x2); */ \
948
image[k+wx*y] = col; }
949
#define PX(n) ps[2*n]
950
#define PY(n) ps[2*n+1]
952
void polyfill(PDL_Long *image, int wx, int wy, float *ps, int n,
953
PDL_Long col, int *ierr)
955
int ymin, ymax, xmin, xmax, fwrd = 1, i, j, k, nsect;
956
int x[MAXSEC], temp, l;
957
float s1, s2, t1, t2;
959
ymin = PY(0); ymax = PY(0);
960
xmin = PX(0); xmax = PX(0);
963
for (i=1; i<n; i++) {
964
ymin = ymin > PY(i) ? PY(i) : ymin;
965
ymax = ymax < PY(i) ? PY(i) : ymax;
966
xmin = xmin > PX(i) ? PX(i) : xmin;
967
xmax = xmax < PX(i) ? PX(i) : xmax;
969
if (xmin < 0 || xmax >= wx || ymin < 0 || ymax >= wy) {
970
*ierr = 1; /* clipping */
975
for (l=ymin; l<= ymax; l++) {
978
for (i=0; i<n; i++) {
981
if ((t1 < l && l <= t2) || (t1 >= l && l > t2)) {
982
if (nsect > MAXSEC) {
983
*ierr = 2; /* too complex */
986
x[nsect] = (s1+(s2-s1)*((l-t1)/(t2-t1)));
992
/* sort the intersections */
993
for (i=1; i<nsect; i++)
1001
for (i=0; i<nsect-1; i += 2)
1002
line(x[i],x[i+1],l);
1005
for (i=nsect-1; i>0; i -= 2)
1006
line(x[i-1],x[i],l);
1015
HandleBad => 0, # a marker
1016
Pars => 'int [o,nc] im(m,n); float ps(two=2,np); int col()',
1017
Code => 'int ierr = 0, nerr;
1019
polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr);
1020
ierr = ierr < nerr ? nerr : ierr;
1022
if (ierr) warn("errors during polygonfilling");
1027
fill the area inside the given polygon with a given colour
1029
This function works inplace, i.e. modifies C<im>.
1034
pp_add_exported('', 'polyfillv');
1041
return the (dataflown) area of an image within a polygon
1045
# increment intensity in area bounded by $poly
1046
$im->polyfillv($pol)++; # legal in perl >= 5.6
1047
# compute average intensity within area bounded by $poly
1048
$av = $im->polyfillv($poly)->avg;
1052
sub PDL::polyfillv {
1054
my $msk = zeroes(long,$im->dims);
1055
polyfill($msk, $ps, 1);
1056
return $im->where($msk == 1);
1058
*polyfillv = \&PDL::polyfillv;
1062
pp_addhdr('#include "rotate.c"'."\n\n");
1063
pp_add_exported('','rotnewsz');
1072
int newcols, newrows;
1074
if (getnewsize(m,n,angle,&newcols,&newrows) != 0)
1075
croak("wrong angle (should be between -90 and +90)");
1077
PUSHs(sv_2mortal(newSVnv(newcols)));
1078
PUSHs(sv_2mortal(newSVnv(newrows)));
1083
Pars => 'im(m,n); float angle(); bg(); int aa(); [o] om(p,q)',
1085
if ((ierr = rotate($P(im),$P(om),$SIZE(m),$SIZE(n),$SIZE(p),
1086
$SIZE(q),$angle(),$bg(),$aa())) != 0)
1088
croak("error during rotate, wrong angle");
1090
croak("wrong output dims, did you set them?");
1092
# ugly workaround since $SIZE(m) and $SIZE(n) are not initialized
1093
# when the redodimscode is called
1095
RedoDimsCode => 'int ncols, nrows;
1096
if ($PDL(im)->ndims < 2)
1097
croak("need > 2d piddle");
1098
if (getnewsize($PDL(im)->dims[0],$PDL(im)->dims[1],
1101
croak("error during rotate, wrong angle");
1102
/* printf("o: %d, p: %d\n",ncols,nrows); */
1105
GenericTypes => ['B'],
1109
rotate an image by given C<angle>
1113
# rotate by 10.5 degrees with antialiasing, set missing values to 7
1114
$rot = $im->rot2d(10.5,7,1);
1116
This function rotates an image through an C<angle> between -90 and + 90
1117
degrees. Uses/doesn't use antialiasing depending on the C<aa> flag.
1118
Pixels outside the rotated image are set to C<bg>.
1120
Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based
1121
on "A Fast Algorithm for General Raster Rotation" by Alan Paeth,
1122
Graphics Interface '86, pp. 77-81.
1124
Use the C<rotnewsz> function to find out about the dimension of the
1127
($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;
1134
Pars => 'I(n,m); O(q,p)',
1138
Bilineary maps the first piddle in the second. The
1139
interpolated values are actually added to the second
1140
piddle which is supposed to be larger than the first one.
1145
int i,j,ii,jj,ii1,jj1,num;
1146
double x,y,dx,dy,y1,y2,y3,y4,t,u,sum;
1148
if ($SIZE(q)>=$SIZE(n) && $SIZE(p)>=$SIZE(m)) {
1150
dx = ((double) ($SIZE(n)-1)) / ($SIZE(q)-1);
1151
dy = ((double) ($SIZE(m)-1)) / ($SIZE(p)-1);
1152
for(i=0,x=0;i<$SIZE(q);i++,x+=dx) {
1153
for(j=0,y=0;j<$SIZE(p);j++,y+=dy) {
1154
ii = (int) floor(x);
1155
if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2;
1156
jj = (int) floor(y);
1157
if (jj>=($SIZE(m)-1)) jj = $SIZE(m)-2;
1160
y1 = $I(n=>ii,m=>jj);
1161
y2 = $I(n=>ii1,m=>jj);
1162
y3 = $I(n=>ii1,m=>jj1);
1163
y4 = $I(n=>ii,m=>jj1);
1166
$O(q=>i,p=>j) += (1-t)*(1-u)*y1 + t*(1-u)*y2 + t*u*y3 + (1-t)*u*y4;
1172
barf("the second matrix must be greater than first! (bilin2d)");
1178
Pars => 'I(n,m); O(q,p)',
1182
The first piddle is rescaled to the dimensions of the second
1183
(expandind or meaning values as needed) and then added to it.
1188
int ix,iy,ox,oy,i,j,lx,ly,cx,cy,xx,yy,num;
1196
if(ox >= ix && oy >= iy) {
1198
kx = ((double) (ox)) / (ix);
1199
ky = ((double) (oy)) / (iy);
1204
cx = rint((i+1)*kx)-1;
1205
cy = rint((j+1)*ky)-1;
1206
for(xx=lx;xx<=cx;xx++)
1207
for(yy=ly;yy<=cy;yy++) {
1208
/* fprintf(stderr,"i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1209
$O(p=>xx,q=>yy) += $I(n=>j,m=>i);
1217
else if(ox < ix && oy < iy) {
1219
kx = ((double) (ix)) / (ox);
1220
ky = ((double) (iy)) / (oy);
1225
cx = rint((i+1)*kx)-1;
1226
cy = rint((j+1)*ky)-1;
1229
for(xx=lx;xx<=cx;xx++)
1230
for(yy=ly;yy<=cy;yy++) {
1231
/* fprintf(stderr,"i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1232
temp += $I(n=>yy,m=>xx);
1235
$O(p=>i,q=>j) += temp/num;
1242
else if(ox >= ix && oy < iy) {
1244
kx = ((double) (ox)) / (ix);
1245
ky = ((double) (iy)) / (oy);
1250
cx = rint((i+1)*kx)-1;
1251
cy = rint((j+1)*ky)-1;
1254
for(yy=ly;yy<=cy;yy++) {
1255
/* fprintf(stderr,"1 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1256
temp += $I(n=>yy,m=>i);
1259
for(xx=lx;xx<=cx;xx++) {
1260
/* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1261
$O(p=>xx,q=>j) += temp/num;
1269
else if(ox < ix && oy >= iy) {
1271
kx = ((double) (ix)) / (ox);
1272
ky = ((double) (oy)) / (iy);
1277
cx = rint((i+1)*kx)-1;
1278
cy = rint((j+1)*ky)-1;
1281
for(xx=lx;xx<=cx;xx++) {
1282
/* fprintf(stderr,"1 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1283
temp += $I(n=>j,m=>xx);
1286
for(yy=ly;yy<=cy;yy++) {
1287
/* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */
1288
$O(p=>i,q=>yy) += temp/num;
1296
else barf("I am not supposed to be here, please report the bug to <chri@infis.univ.ts.it>");
1299
# functions to make handling 2D polynomial mappings a bit easier
1301
pp_add_exported('', 'fitwarp2d applywarp2d');
1309
Find the best-fit 2D polynomial to describe
1310
a coordinate transformation.
1314
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf. { options } )
1316
Given a set of points in the output plane (C<$u,$v>), find
1317
the best-fit (using singular-value decomposition) 2D polynomial
1318
to describe the mapping back to the image plane (C<$x,$y>).
1319
The order of the fit is controlled by the C<$nf> parameter
1320
(the maximum power of the polynomial is C<$nf - 1>), and you
1321
can restrict the terms to fit using the C<FIT> option.
1323
C<$px> and C<$py> are C<np> by C<np> element piddles which describe
1324
a polynomial mapping (of order C<np-1>)
1325
from the I<output> C<(u,v)> image to the I<input> C<(x,y)> image:
1327
x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
1328
y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j
1330
The transformation is returned for the reverse direction (ie
1331
output to input image) since that is what is required by the
1332
L<warp2d()|/warp2d> routine. The L<applywarp2d()|/applywarp2d>
1333
routine can be used to convert a set of C<$u,$v> points given
1340
FIT - which terms to fit? default ones(byte,$nf,$nf)
1341
THRESH - in svd, remove terms smaller than THRESH * max value
1348
C<FIT> allows you to restrict which terms of the polynomial to fit:
1349
only those terms for which the FIT piddle evaluates to true will be
1350
evaluated. If a 2D piddle is sent in, then it is
1351
used for the x and y polynomials; otherwise
1352
C<$fit-E<gt>slice(":,:,(0)")> will be used for C<$px> and
1353
C<$fit-E<gt>slice(":,:,(1)")> will be used for C<$py>.
1357
Remove all singular values whose valus is less than C<THRESH>
1358
times the largest singular value.
1362
The number of points must be at least equal to the number of
1363
terms to fit (C<$nf*$nf> points for the default value of C<FIT>).
1367
# points in original image
1368
$x = pdl( 0, 0, 100, 100 );
1369
$y = pdl( 0, 100, 100, 0 );
1370
# get warped to these positions
1371
$u = pdl( 10, 10, 90, 90 );
1372
$v = pdl( 10, 90, 90, 10 );
1374
# shift of origin + scale x/y axis only
1375
$fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
1376
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
1377
print "px = ${px}py = $py";
1389
# Compared to allowing all 4 terms
1390
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
1391
print "px = ${px}py = $py";
1395
[ 1.110223e-16 -1.1275703e-17]
1399
[ -12.5 1.6653345e-16]
1400
[ 1.25 -5.8546917e-18]
1407
Transform a set of points using a 2-D polynomial mapping
1411
( $x, $y ) = applywarp2d( $px, $py, $u, $v )
1413
Convert a set of points (stored in 1D piddles C<$u,$v>)
1414
to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px>
1415
and C<$py>. See L<fitwarp2d()|/fitwarp2d>
1416
for more information on the format of C<$px> and C<$py>.
1420
# use SVD to fit data. Assuming no errors.
1426
# if we had errors for these points, would normalise the
1427
# basis functions, and the output array, by these errors here
1430
my ( $svd_u, $svd_w, $svd_v ) = svd( $basis );
1432
# remove any singular values
1433
$svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
1435
# perform the back substitution
1437
my $tmp = $y x $svd_u;
1438
if ( $PDL::Bad::Status ) {
1439
$tmp /= $svd_w->setvaltobad(0.0);
1440
$tmp->inplace->setbadtoval(0.0);
1443
my $mask = ($svd_w == 0.0);
1444
$tmp /= ( $svd_w + $mask );
1445
$tmp *= ( 1 - $mask );
1448
my $ans = sumover( $svd_v * $tmp );
1454
sub _mkbasis ($$$$) {
1460
my $n = $fit->getdim(0) - 1;
1461
my $ncoeff = sum( $fit );
1463
my $basis = zeroes( $u->type, $ncoeff, $npts );
1465
foreach my $j ( 0 .. $n ) {
1467
foreach my $i ( 0 .. $n ) {
1468
if ( $fit->at($i,$j) ) {
1469
my $tmp = $basis->slice("($k),:");
1470
$tmp .= $tmp_v * $u**$i;
1479
sub PDL::fitwarp2d {
1480
croak "Usage: (\$px,\$py) = fitwarp2d(x(m);y(m);u(m);v(m);\$nf; { options })"
1481
if $#_ < 4 or ( $#_ >= 5 and ref($_[5]) ne "HASH" );
1489
my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf), THRESH => 1.0e-5 } );
1490
$opts->options( $_[0] ) if $#_ > -1;
1491
my $oref = $opts->current();
1494
my $npts = $x->nelem;
1495
croak "fitwarp2d: x, y, u, and v must be the same size (and 1D)"
1496
unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem
1497
and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1;
1499
my $svd_thresh = $$oref{THRESH};
1500
croak "fitwarp2d: THRESH option must be >= 0."
1503
my $fit = $$oref{FIT};
1504
my $fit_ndim = $fit->getndims();
1505
croak "fitwarp2d: FIT option must be sent a (\$nf,\$nf[,2]) element piddle"
1506
unless UNIVERSAL::isa($fit,"PDL") and
1507
($fit_ndim == 2 or ($fit_ndim == 3 and $fit->getdim(2) == 2)) and
1508
$fit->getdim(0) == $nf and $fit->getdim(1) == $nf;
1510
# how many coeffs to fit (first we ensure $fit is either 0 or 1)
1511
$fit = convert( $fit != 0, byte );
1513
my ( $fitx, $fity, $ncoeffx, $ncoeffy, $ncoeff );
1514
if ( $fit_ndim == 2 ) {
1517
$ncoeff = $ncoeffx = $ncoeffy = sum( $fit );
1519
$fitx = $fit->slice(",,(0)");
1520
$fity = $fit->slice(",,(1)");
1521
$ncoeffx = sum($fitx);
1522
$ncoeffy = sum($fity);
1523
$ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy;
1526
croak "fitwarp2d: number of points must be >= \$ncoeff"
1527
unless $npts >= $ncoeff;
1529
# create the basis functions for the SVD fitting
1530
my ( $basisx, $basisy );
1531
$basisx = _mkbasis( $fitx, $npts, $u, $v );
1532
if ( $fit_ndim == 2 ) {
1535
$basisy = _mkbasis( $fity, $npts, $u, $v );
1538
my $px = _svd( $basisx, $x, $svd_thresh );
1539
my $py = _svd( $basisy, $y, $svd_thresh );
1541
# convert into $nf x $nf element piddles, if necessary
1542
my $nf2 = $nf * $nf;
1544
return ( $px->reshape($nf,$nf), $py->reshape($nf,$nf) )
1545
if $ncoeff == $nf2 and $ncoeffx == $ncoeffy;
1547
# re-create the matrix
1548
my $xtmp = zeroes( $nf, $nf );
1549
my $ytmp = zeroes( $nf, $nf );
1553
foreach my $i ( 0 .. ($nf - 1) ) {
1554
foreach my $j ( 0 .. ($nf - 1) ) {
1555
if ( $fitx->at($i,$j) ) {
1556
$xtmp->set($i,$j, $px->at($kx) );
1559
if ( $fity->at($i,$j) ) {
1560
$ytmp->set($i,$j, $py->at($ky) );
1566
return ( $xtmp, $ytmp )
1570
*fitwarp2d = \&PDL::fitwarp2d;
1572
sub PDL::applywarp2d {
1574
croak "Usage: (\$x,\$y) = applywarp2d(px(nf,nf);py(nf,nf);u(m);v(m);)"
1581
my $npts = $u->nelem;
1584
croak "applywarp2d: u and v must be the same size (and 1D)"
1585
unless $npts == $u->nelem and $npts == $v->nelem
1586
and $u->getndims == 1 and $v->getndims == 1;
1588
my $nf = $px->getdim(0);
1589
my $nf2 = $nf * $nf;
1591
# could remove terms with 0 coeff here
1592
# (would also have to remove them from px/py for
1593
# the matrix multiplication below)
1595
my $mat = _mkbasis( ones(byte,$nf,$nf), $npts, $u, $v );
1597
my $x = reshape( $mat x $px->clump(-1)->transpose(), $npts );
1598
my $y = reshape( $mat x $py->clump(-1)->transpose(), $npts );
1601
} # sub: applywarp2d
1603
*applywarp2d = \&PDL::applywarp2d;
1607
## resampling routines taken from v3.6-0 of the Eclipse package
1608
## http://www.eso.org/eclipse by Nicolas Devillard
1611
pp_addhdr( '#include "resample.h"' . "\n" );
1614
# and support routine
1622
Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { options })
1626
Warp a 2D image given a polynomial describing the I<reverse> mapping.
1630
$out = warp2d( $img, $px, $py, { options } );
1632
Apply the polynomial transformation encoded in the C<$px> and
1633
C<$py> piddles to warp the input image C<$img> into the output
1636
The format for the polynomial transformation is described in
1637
the documentation for the L<fitwarp2d()|/fitwarp2d> routine.
1639
At each point C<x,y>, the closest 16 pixel values are combined
1640
with an interpolation kernel to calculate the value at C<u,v>.
1641
The interpolation is therefore done in the image, rather than
1643
By default, a C<tanh> kernel is used, but this can be changed
1644
using the C<KERNEL> option discussed below
1645
(the choice of kernel depends on the frequency content of the input image).
1647
The routine is based on the C<warping> command from
1648
the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and
1649
for further details on image resampling see
1650
Wolberg, G., "Digital Image Warping", 1990, IEEE Computer
1651
Society Press ISBN 0-8186-8944-7).
1653
Currently the output image is the same size as the input one,
1654
which means data will be lost if the transformation reduces
1655
the pixel scale. This will (hopefully) be changed soon.
1659
$img = rvals(byte,501,501);
1660
imag $img, { JUSTIFY => 1 };
1662
# use a not-particularly-obvious transformation:
1663
# x = -10 + 0.5 * $u - 0.1 * $v
1664
# y = -20 + $v - 0.002 * $u * $v
1666
$px = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
1667
$py = pdl( [ -20, 0 ], [ 1, 0.002 ] );
1668
$wrp = warp2d( $img, $px, $py );
1670
# see the warped image
1671
imag $warp, { JUSTIFY => 1 };
1677
KERNEL - default value is tanh
1678
NOVAL - default value is 0
1680
C<KERNEL> is used to specify which interpolation kernel to use
1681
(to see what these kernels look like, use the
1682
L<warp2d_kernel()|/warp2d_kernel> routine).
1689
Hyperbolic tangent: the approximation of an ideal box filter by the
1690
product of symmetric tanh functions.
1694
For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle,
1695
which produces a C<sinc> interpolation kernel in the spatial domain:
1697
sinc(x) = sin(pi * x) / (pi * x)
1699
However, it is not ideal for the C<4x4> pixel region used here.
1703
This is the square of the sinc function.
1707
Although defined differently to the C<tanh> kernel, the result is very
1708
similar in the spatial domain. The Lanczos function is defined as
1710
L(x) = sinc(x) * sinc(x/2) if abs(x) < 2
1715
This kernel is derived from the following function:
1717
H(x) = a + (1-a) * cos(2*pi*x/(N-1)) if abs(x) < 0.5*(N-1)
1720
with C<a = 0.5> and N currently equal to 2001.
1724
This kernel uses the same C<H(x)> as the Hann filter, but with
1729
C<NOVAL> gives the value used to indicate that a pixel in the
1730
output image does not map onto one in the input image.
1736
my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann );
1738
# note: convert to lower case
1739
sub _check_kernel ($$) {
1740
my $kernel = lc shift;
1742
barf "Unknown kernel $kernel sent to $code\n" .
1743
"\tmust be one of [" . join(',',keys %warp2d) . "]\n"
1744
unless exists $warp2d{$kernel};
1755
Pars => 'img(m,n); double px(np,np); double py(np,np); [o] warp(m,n);',
1756
OtherPars => 'char *kernel_type; double noval;',
1757
GenericTypes => [ 'F', 'D' ],
1761
my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );
1762
$opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";
1764
die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} )"
1769
my $out = $#_ == -1 ? PDL->null() : shift;
1772
my $copt = $opts->current();
1773
my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );
1775
&PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} );
1783
int ncoeff, lx_out, ly_out ;
1786
double neighbors[16] ;
1792
double *kernel, *poly ;
1793
int da[16], db[16] ;
1795
/* Generate default interpolation kernel */
1796
kernel = generate_interpolation_kernel( $COMP(kernel_type) ) ;
1797
if (kernel == NULL) {
1798
croak( "Ran out of memory building kernel\n" );
1803
lx_out = $SIZE(m); /* is this right? */
1808
/* Pre compute leaps for 16 closest neighbors positions */
1809
da[0] = -1; db[0] = -1;
1810
da[1] = 0; db[1] = -1;
1811
da[2] = 1; db[2] = -1;
1812
da[3] = 2; db[3] = -1;
1814
da[4] = -1; db[4] = 0;
1815
da[5] = 0; db[5] = 0;
1816
da[6] = 1; db[6] = 0;
1817
da[7] = 2; db[7] = 0;
1819
da[8] = -1; db[8] = 1;
1820
da[9] = 0; db[9] = 1;
1821
da[10] = 1; db[10] = 1;
1822
da[11] = 2; db[11] = 1;
1824
da[12] = -1; db[12] = 2;
1825
da[13] = 0; db[13] = 2;
1826
da[14] = 1; db[14] = 2;
1827
da[15] = 2; db[15] = 2;
1829
/* allocate memory for polynomial */
1830
poly = malloc( ncoeff * sizeof(double) );
1831
if ( poly == NULL ) {
1832
croak( "Ran out of memory\n" );
1836
/* Loop over the output image */
1840
/* fill in poly array */
1841
for ( k = 1; k < ncoeff; k++ ) {
1842
poly[k] = (double) n * poly[k-1];
1847
/* Compute the original source for this pixel */
1848
x = poly2d_compute( ncoeff, $P(px), (double) m, poly );
1849
y = poly2d_compute( ncoeff, $P(py), (double) m, poly );
1851
/* Which is the closest integer positioned neighbor? */
1855
if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
1856
$warp() = ($GENERIC()) $COMP(noval);
1859
/* Now feed the positions for the closest 16 neighbors */
1860
for (k=0 ; k<16 ; k++) {
1861
i = px + da[k]; j = py + db[k];
1862
neighbors[k] = (double) $img( m => i, n => j );
1865
/* Which tabulated value index shall we use? */
1866
tabx = (x - (double)px) * (double)(TABSPERPIX) ;
1867
taby = (y - (double)py) * (double)(TABSPERPIX) ;
1869
/* Compute resampling coefficients */
1870
/* rsc[0..3] in x, rsc[4..7] in y */
1872
rsc[0] = kernel[TABSPERPIX + tabx] ;
1873
rsc[1] = kernel[tabx] ;
1874
rsc[2] = kernel[TABSPERPIX - tabx] ;
1875
rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
1876
rsc[4] = kernel[TABSPERPIX + taby] ;
1877
rsc[5] = kernel[taby] ;
1878
rsc[6] = kernel[TABSPERPIX - taby] ;
1879
rsc[7] = kernel[2 * TABSPERPIX - taby] ;
1881
sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
1882
(rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
1884
/* Compute interpolated pixel now */
1885
cur = rsc[4] * ( rsc[0]*neighbors[0] +
1886
rsc[1]*neighbors[1] +
1887
rsc[2]*neighbors[2] +
1888
rsc[3]*neighbors[3] ) +
1889
rsc[5] * ( rsc[0]*neighbors[4] +
1890
rsc[1]*neighbors[5] +
1891
rsc[2]*neighbors[6] +
1892
rsc[3]*neighbors[7] ) +
1893
rsc[6] * ( rsc[0]*neighbors[8] +
1894
rsc[1]*neighbors[9] +
1895
rsc[2]*neighbors[10] +
1896
rsc[3]*neighbors[11] ) +
1897
rsc[7] * ( rsc[0]*neighbors[12] +
1898
rsc[1]*neighbors[13] +
1899
rsc[2]*neighbors[14] +
1900
rsc[3]*neighbors[15] ) ;
1902
/* Copy the value to the output image */
1903
$warp() = ($GENERIC()) (cur/sumrs);
1905
} /* if: edge or interior */
1923
RETVAL = KERNEL_SAMPLES;
1929
pp_add_exported('', 'warp2d_kernel');
1932
=head2 warp2d_kernel
1936
Return the specified kernel, as used by L<warp2d|/warp2d>
1940
( $x, $k ) = warp2d_kernel( $name )
1942
The valid values for C<$name> are the same as the C<KERNEL> option
1943
of L<warp2d()|/warp2d>.
1947
line warp2d_kernel( "hamming" );
1953
# this is not very clever, but it's a pain to create a valid
1956
pp_def( 'warp2d_kernel',
1961
sub PDL::warp2d_kernel ($) {
1962
my $kernel = _check_kernel( shift, "warp2d_kernel" );
1964
my $nelem = _get_kernel_size();
1965
my $x = zeroes( $nelem );
1966
my $k = zeroes( $nelem );
1968
&PDL::_warp2d_kernel_int( $x, $k, $kernel );
1971
# return _get_kernel( $kernel );
1973
*warp2d_kernel = \&PDL::warp2d_kernel;
1976
Pars => '[o] x(n); [o] k(n);',
1977
OtherPars => 'char *name;',
1978
GenericTypes => [ 'D' ],
1983
if ( $SIZE(n) != KERNEL_SAMPLES ) {
1984
croak( "Internal error in warp2d_kernel - mismatch in kernel size\n" );
1987
kernel = generate_interpolation_kernel($COMP(name));
1988
if ( kernel == NULL ) {
1989
croak( "unable to allocate memory for kernel" );
1992
/* fill in piddles */
1998
xx += 1.0 / (double) TABSPERPIX;
2002
/* free the kernel */
added
removed
Lib/Image2D/image2d.pd
