use strict; use PDL::Types; pp_addpm({At=>'Top'},<<'EOD'); =head1 NAME PDL::Image2D - Miscellaneous 2D image processing functions =head1 DESCRIPTION Miscellaneous 2D image processing functions - for want of anywhere else to put them. =head1 SYNOPSIS use PDL::Image2D; =cut use PDL; # ensure qsort routine available use PDL::Math; use Carp; use strict; EOD pp_addpm({At=>'Bot'},<<'EOD'); =head1 AUTHORS Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams (rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu), and Doug Burke (burke@ifa.hawaii.edu). All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file. =cut EOD pp_addhdr(' #define IsNaN(x) (x != x) /* Fast Modulus with proper negative behaviour */ #define REALMOD(a,b) {while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);} /* rint is missing on some platforms (eg Win32) */ #ifdef NEEDS_RINT #define rint(X) floor( X + 0.5 ) #endif '); for (keys %PDL::Types::typehash) { my $ctype = $PDL::Types::typehash{$_}{ctype}; my $ppsym = $PDL::Types::typehash{$_}{ppsym}; pp_addhdr << "EOH"; /* * this routine is based on code referenced from * http://www.eso.org/~ndevilla/median/ * the original algorithm is described in Numerical Recipes */ #define ELEM_SWAP(a,b) { register $ctype t=(a);(a)=(b);(b)=t; } $ctype quick_select_$ppsym($ctype arr[], int n) { int low, high ; int median; int middle, ll, hh; low = 0 ; high = n-1 ; median = (low + high) / 2; for (;;) { if (high <= low) /* One element only */ return arr[median] ; if (high == low + 1) { /* Two elements only */ if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; return arr[median] ; } /* Find median of low, middle and high items; swap into position low */ middle = (low + high) / 2; if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; /* Swap low item (now in position middle) into position (low+1) */ ELEM_SWAP(arr[middle], arr[low+1]) ; /* Nibble from each end towards middle, swapping items when stuck */ ll = low + 1; hh = high; for (;;) { do ll++; while (arr[low] > arr[ll]) ; do hh--; while (arr[hh] > arr[low]) ; if (hh < ll) break; ELEM_SWAP(arr[ll], arr[hh]) ; } /* Swap middle item (in position low) back into correct position */ ELEM_SWAP(arr[low], arr[hh]) ; /* Re-set active partition */ if (hh <= median) low = ll; if (hh >= median) high = hh - 1; } } #undef ELEM_SWAP EOH } my %init = ( i => { size => 'm_size', off => 'poff', init => '1-p_size' }, j => { size => 'n_size', off => 'qoff', init => '1-q_size' }, ); # requires 'int $var, ${var}2' to have been declared in the c code # (along with [pq]off and [pq]_size) # sub init_map { my $var = shift; my $loop = $var; my $loop2 = "${var}2"; my $href = $init{$var} || die "ERROR: unknown variable sent to init_map()\n"; my $size = $href->{size} || die "ERROR: unable to find size for $var\n"; my $off = $href->{off} || die "ERROR: unable to find off for $var\n"; my $init = $href->{init} || die "ERROR: unable to find init for $var\n"; return "for ( $loop = $init; $loop< $size; ${loop}++) { $loop2 = $loop + $off; switch (opt) { case 1: /* REFLECT */ if (${loop2}<0) $loop2 = -${loop2}-1; else if ($loop2 >= $size) $loop2 = 2*${size}-(${loop2}+1); break; case 2: /* TRUNCATE */ if (${loop2}<0 || ${loop2} >= $size) $loop2 = -1; break; default: REALMOD($loop2,$size); } map${var}\[$loop] = $loop2; }\n"; } # sub: init_map() sub init_vars { my $href = shift || { }; $href->{vars} = '' unless defined $href->{vars}; $href->{malloc} = '' unless defined $href->{malloc}; $href->{check} = '' unless defined $href->{check}; my $str = $href->{vars}; $str .= "int i,j, i1,j1, i2,j2, poff, qoff;"; $str .= 'int opt = $COMP(opt); int m_size = $COMP(__m_size); int n_size = $COMP(__n_size); int p_size = $COMP(__p_size); int q_size = $COMP(__q_size); int *mapi, *mapj; mapi = (int *) malloc((p_size+m_size)*sizeof(int)); mapj = (int *) malloc((q_size+n_size)*sizeof(int)); '; $str .= $href->{malloc} . "\n"; $str .= "if ($href->{check} (mapi==NULL) || (mapj==NULL))\n"; $str .= ' barf("Out of Memory"); poff = p_size/2; mapi += p_size-1; qoff = q_size/2; mapj += q_size-1; '; return $str; } # sub: init_vars() pp_def('conv2d', Doc=><<'EOD', =for ref 2D convolution of an array with a kernel (smoothing) For large kernels, using a FFT routine, such as L in C, will be quicker. =for usage $new = conv2d $old, $kernel, {OPTIONS} =for example $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect} =for options Boundary - controls what values are assumed for the image when kernel crosses its edge: => Default - periodic boundary conditions (i.e. wrap around axis) => Reflect - reflect at boundary => Truncate - truncate at boundary EOD BadDoc => 'Unlike the FFT routines, conv2d is able to process bad values.', HandleBad => 1, Pars => 'a(m,n); kern(p,q); [o]b(m,n);', OtherPars => 'int opt;', PMCode => ' sub PDL::conv2d { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' if $#_<1 || $#_>2; my($a,$kern) = @_; my $c = $#_ == 2 ? $_[2] : $a->nullcreate; &PDL::_conv2d_int($a,$kern,$c, (!(defined $opt && exists $$opt{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; } ', Code => init_vars( { vars => 'PDL_Double tmp;' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j= 0) { for(i1=0; i1= 0) tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1); } /* for: i1 */ } /* if: j2 >= 0 */ } /* for: j1 */ $b(m=>i,n=>j) = tmp; } /* for: i */ } /* for: j */ %} free(mapj+1-q_size); free(mapi+1-p_size);', BadCode => init_vars( { vars => 'PDL_Double tmp; int flag;' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j= 0) { for(i1=0; i1= 0) { if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(kern(p=>i1,q=>j1)) ) { tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1); flag = 1; } /* if: good */ } /* if: i2 >= 0 */ } /* for: i1 */ } /* if: j2 >= 0 */ } /* for: j1 */ if ( flag ) { $b(m=>i,n=>j) = tmp; } else { $SETBAD(b(m=>i,n=>j)); } } /* for: i */ } /* for: j */ %} free(mapj+1-q_size); free(mapi+1-p_size);', ); # pp_def: conv2d pp_def('med2d', Doc=> <<'EOD', =for ref 2D median-convolution of an array with a kernel (smoothing) Note: only points in the kernel E0 are included in the median, other points are weighted by the kernel value (medianing lots of zeroes is rather pointless) =for usage $new = med2d $old, $kernel, {OPTIONS} =for example $smoothed = med2d $image, ones(3,3), {Boundary => Reflect} =for options Boundary - controls what values are assumed for the image when kernel crosses its edge: => Default - periodic boundary conditions (i.e. wrap around axis) => Reflect - reflect at boundary => Truncate - truncate at boundary EOD BadDoc => 'Bad values are ignored in the calculation. If all elements within the kernel are bad, the output is set bad.', HandleBad => 1, Pars => 'a(m,n); kern(p,q); [o]b(m,n);', OtherPars => 'int opt;', PMCode => ' sub PDL::med2d { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' if $#_<1 || $#_>2; my($a,$kern) = @_; croak "med2d: kernel must contain some positive elements.\n" if all( $kern <= 0 ); my $c = $#_ == 2 ? $_[2] : $a->nullcreate; &PDL::_med2d_int($a,$kern,$c, (!(defined $opt && exists $opt->{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; } ', Code => init_vars( { vars => 'PDL_Double *tmp, kk; int count;', malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));', check => '(tmp==NULL) || ' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j= 0) for(i1=0; i1= 0) { kk = $kern(p=>i1,q=>j1); if (kk>0) { tmp[count++] = $a(m=>i2,n=>j2) * kk; } } /* if: i2 >= 0 */ } /* for: i1 */ } /* for: j1 */ PDL->qsort_D( tmp, 0, count-1 ); $b(m=>i,n=>j) = tmp[(count-1)/2]; } /* for: i */ } /* for: j */ %} free(mapj+1-q_size); free(mapi+1-p_size); free(tmp); ', BadCode => init_vars( { vars => 'PDL_Double *tmp, kk, aa; int count, flag;', malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));', check => '(tmp==NULL) || ' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j= 0) for(i1=0; i1= 0) { kk = $kern(p=>i1,q=>j1); aa = $a(m=>i2,n=>j2); if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(aa,a) ) { flag = 1; if ( kk > 0 ) { tmp[count++] = aa * kk; } } } /* if: i2 >= 0 */ } /* for: i1 */ } /* for: j1 */ if ( flag == 0 ) { $SETBAD(b(m=>i,n=>j)); } else { PDL->qsort_D( tmp, 0, count-1 ); $b(m=>i,n=>j) = tmp[(count-1)/2]; } } /* for: i */ } /* for: j */ %} free(mapj+1-q_size); free(mapi+1-p_size); free(tmp); ' ); # pp_def: med2d pp_def('med2df', Doc=> <<'EOD', =for ref 2D median-convolution of an array in a pxq window (smoothing) Note: this routine does the median over all points in a rectangular window and is not quite as flexible as C in this regard but slightly faster instead =for usage $new = med2df $old, $xwidth, $ywidth, {OPTIONS} =for example $smoothed = med2df $image, 3, 3, {Boundary => Reflect} =for options Boundary - controls what values are assumed for the image when kernel crosses its edge: => Default - periodic boundary conditions (i.e. wrap around axis) => Reflect - reflect at boundary => Truncate - truncate at boundary EOD Pars => 'a(m,n); [o]b(m,n);', # funny parameter names to avoid special case in 'init_vars' OtherPars => 'int __p_size; int __q_size; int opt;', PMCode => ' sub PDL::med2df { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )\' if $#_<2 || $#_>3; my($a,$p,$q) = @_; croak "med2df: kernel must contain some positive elements.\n" if $p == 0 && $q == 0; my $c = $#_ == 3 ? $_[3] : $a->nullcreate; &PDL::_med2df_int($a,$c,$p,$q, (!(defined $opt && exists $opt->{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; } ', Code => init_vars( { vars => '$GENERIC() *tmp, kk; int count;', malloc => 'tmp = malloc(p_size*q_size*sizeof($GENERIC()));', check => '(tmp==NULL) || ' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j= 0) for(i1=0; i1= 0) { tmp[count++] = $a(m=>i2,n=>j2); } /* if: i2 >= 0 */ } /* for: i1 */ } /* for: j1 */ $b(m=>i,n=>j) = quick_select_$TBSULQFD(B,S,U,L,Q,F,D) (tmp, count ); } /* for: i */ } /* for: j */ %} free(mapj+1-q_size); free(mapi+1-p_size); free(tmp); ', ); # pp_def: med2df pp_addhdr(<<'EOH'); #define EZ(x) ez ? 0 : (x) EOH pp_def('box2d', Pars => 'a(n,m); [o] b(n,m)', OtherPars => 'int wx; int wy; int edgezero', Code => ' register int nx = 0.5*$COMP(wx); register int ny = 0.5*$COMP(wy); register int xs = $SIZE(n); register int ys = $SIZE(m); register int ez = $COMP(edgezero); double div, sum, lsum; int xx,yy,y,ind1,ind2,first; div = 1/((2.0*nx+1)*(2.0*ny+1)); threadloop %{ first = 1; for (y=0;yxx,m=>y) = EZ($a(n=>xx,m=>y)); $b(n=>ind1,m=>y) = EZ($a(n=>ind1,m=>y)); } for (xx=0;xxxx,m=>y) = EZ($a(n=>xx,m=>y)); $b(n=>xx,m=>ind1) = EZ($a(n=>xx,m=>ind1)); } for (y=ny;yxx,m=>yy); } else { ind1 = y-ny-1; ind2 = y+ny; for (xx=0;xx<=2*nx;xx++) { lsum -= $a(n=>xx,m=>ind1); /* remove top pixels */ lsum += $a(n=>xx,m=>ind2); /* add bottom pixels */ } } sum = lsum; $b(n=>nx,m=>y) = div*sum; /* and assign */ for (xx=nx+1;xxind1,m=>yy); /* remove leftmost data */ sum += $a(n=>ind2,m=>yy); /* and add rightmost */ } $b(n=>xx,m=>y) = div*sum; /* and assign */ } } %}', Doc => << 'EOD', =for ref fast 2D boxcar average =for example $smoothim = $im->box2d($wx,$wy,$edgezero=1); The edgezero argument controls if edge is set to zero (edgezero=1) or just keeps the original (unfiltered) values. C should be updated to support similar edge options as C and C etc. Boxcar averaging is a pretty crude way of filtering. For serious stuff better filters are around (e.g., use L with the appropriate kernel). On the other hand it is fast and computational cost grows only approximately linearly with window size. EOD ); # pp_def box2d pp_def('patch2d', Doc=><<'EOD', =for ref patch bad pixels out of 2D images using a mask =for usage $patched = patch2d $data, $bad; C<$bad> is a 2D mask array where 1=bad pixel 0=good pixel. Pixels are replaced by the average of their non-bad neighbours; if all neighbours are bad, the original data value is copied across. EOD BadDoc => 'This routine does not handle bad values - use L instead', HandleBad => 0, Pars => 'a(m,n); int bad(m,n); [o]b(m,n);', Code => 'int m_size, n_size, i,j, i1,j1, i2,j2, norm; double tmp; m_size = $COMP(__m_size); n_size = $COMP(__n_size); threadloop %{ for(j=0; ji,n=>j) = $a(m=>i,n=>j); if ( $bad(m=>i,n=>j)==1 ) { tmp = 0; norm=0; for(j1=-1; j1<=1; j1++) { j2 = j+j1; if ( j2>=0 && j2=0 && i2i2,n=>j2)!=1 ) { tmp += $a(m=>i2,n=>j2); norm++; } } /* if: i1!=0 || j1!=0 */ } /* for: i1 */ } } /* for: j1 */ if (norm>0) { /* Patch */ $b(m=>i,n=>j) = tmp/norm; } } /* if: bad() */ } /* for: i */ } /* for: j */ %} /* threadloop */ ', # Code ); pp_def('patchbad2d', Doc=><<'EOD', =for ref patch bad pixels out of 2D images containing bad values =for usage $patched = patchbad2d $data; Pixels are replaced by the average of their non-bad neighbours; if all neighbours are bad, the output is set bad. If the input piddle contains I bad values, then a straight copy is performed (see L). EOD BadDoc => 'patchbad2d handles bad values. The output piddle I contain bad values, depending on the pattern of bad values in the input piddle.', HandleBad => 1, Pars => 'a(m,n); [o]b(m,n);', Code => 'loop(n,m) %{ $b() = $a(); %}', # just copy CopyBadStatusCode => '', # handled by BadCode BadCode => 'int m_size, n_size, i,j, i1,j1, i2,j2, norm, flag; double tmp; $GENERIC(a) a_val; flag = 0; m_size = $COMP(__m_size); n_size = $COMP(__n_size); threadloop %{ for(j=0; ji,n=>j); if ( $ISGOODVAR(a_val,a) ) { $b(m=>i,n=>j) = a_val; } else { tmp = 0; norm=0; for(j1=-1; j1<=1; j1++) { j2 = j+j1; if ( j2>=0 && j2=0 && i2i2,n=>j2); if ( $ISGOODVAR(a_val,a) ) { tmp += a_val; norm++; } } } /* if: i1!=0 || j1!=0 */ } /* for: i1 */ } } /* for: j1 */ /* Patch */ if (norm>0) { $b(m=>i,n=>j) = tmp/norm; } else { $SETBAD(b(m=>i,n=>j)); flag = 1; } } /* if: ISGOODVAR() */ } /* for: i */ } /* for: j */ %} /* threadloop */ /* handle bad flag */ if ( flag ) $PDLSTATESETBAD(b); ', # BadCode ); pp_def('max2d_ind', Doc=><<'EOD', =for ref Return value/position of maximum value in 2D image Contributed by Tim Jeness EOD BadDoc=><<'EOD', Bad values are excluded from the search. If all pixels are bad then the output is set bad. EOD HandleBad => 1, Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();', Code => ' double cur; int curind1; int curind2; curind1=0; curind2=0; loop(m) %{ loop(n) %{ if((!m && !n) || $a() > cur || IsNaN(cur)) { cur = $a(); curind1 = m; curind2 = n; } %} %} $val() = cur; $x() = curind1; $y() = curind2; ', BadCode => ' double cur; int curind1; int curind2; curind1 = -1; curind2 = -1; loop(m) %{ loop(n) %{ if( $ISGOOD(a()) && ( (!n && !m) || ($a() > cur) ) ) { cur = $a(); curind1 = m; curind2 = n; } %} %} if ( curind1 < 0 ) { $SETBAD(val()); $SETBAD(x()); $SETBAD(y()); } else { $val() = cur; $x() = curind1; $y() = curind2; } '); pp_def('centroid2d', Doc=><<'EOD', =for ref Refine a list of object positions in 2D image by centroiding in a box C<$box> is the full-width of the box, i.e. the window is C<+/- $box/2>. EOD BadDoc=><<'EOD', Bad pixels are excluded from the centroid calculation. If all elements are bad (or the pixel sum is 0 - but why would you be centroiding something with negatives in...) then the output values are set bad. EOD HandleBad => 1, Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();', Code => ' int i,j,i1,i2,j1,j2,m_size,n_size; double sum,data,sumx,sumy; m_size = $SIZE(m); n_size = $SIZE(n); i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1; i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2; j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1; j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2; sum = sumx = sumy = 0; for(j=j1; j<=j2; j++) { for(i=i1; i<=i2; i++) { data = $im(m=>i,n=>j); sum += data; sumx += data*i; sumy += data*j; }} $xcen() = sumx/sum; $ycen() = sumy/sum; ', BadCode => ' int i,j,i1,i2,j1,j2,m_size,n_size; double sum,data,sumx,sumy; m_size = $SIZE(m); n_size = $SIZE(n); i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1; i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2; j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1; j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2; sum = sumx = sumy = 0; for(j=j1; j<=j2; j++) { for(i=i1; i<=i2; i++) { data = $im(m=>i,n=>j); if ( $ISGOODVAR(data,im) ) { sum += data; sumx += data*i; sumy += data*j; } } } /* * if sum == 0 then we will flag as bad -- although it could just mean that * there is negative values in the dataset. * - should use a better check than != 0.0 ... */ if ( sum != 0.0 ) { $xcen() = sumx/sum; $ycen() = sumy/sum; } else { $SETBAD(xcen()); $SETBAD(ycen()); } ' ); pp_addhdr(' /* Add an equivalence to a list - used by pdl_cc8compt */ void AddEquiv ( PDL_Long* equiv, PDL_Long i, PDL_Long j) { PDL_Long k, tmp; if (i==j) return; k = j; do { k = equiv[k]; } while ( k != j && k != i ); if ( k == j ) { tmp = equiv[i]; equiv[i] = equiv[j]; equiv[j] = tmp; } } '); pp_def('cc8compt',Doc=>' =for ref Connected 8-component labeling of a binary image. Connected 8-component labeling of 0,1 image - i.e. find seperate segmented objects and fill object pixels with object number =for example $segmented = cc8compt( $image > $threshold ); ', HandleBad => 0, # a marker Pars => 'a(m,n); [o]b(m,n);', Code => ' PDL_Long i,j,k; PDL_Long newlabel; PDL_Long neighbour[4]; PDL_Long nfound; PDL_Long pass,count,next,this; PDL_Long *equiv; PDL_Long i1,j1,i2; PDL_Long nx = $SIZE(m); PDL_Long ny = $SIZE(n); loop(n) %{ loop(m) %{ /* Copy */ $b() = $a(); %} %} /* 1st pass counts max possible compts, 2nd records equivalences */ for (pass = 0; pass<2; pass++) { if (pass==1) { equiv = (PDL_Long*) malloc((newlabel+1)*sizeof(PDL_Long)); if (equiv==(PDL_Long*)0) barf("Out of memory"); for(i=0;i<=newlabel;i++) equiv[i]=i; } newlabel = 1; /* Running label */ for(j=0; j0 */ i1 = i-1; j1 = j-1; i2 = i+1; if ($b(m=>i, n=>j) > 0) { /* Check 4 neighbour already seen */ if (i>0 && $b(m=>i1, n=>j)>0) neighbour[nfound++] = $b(m=>i1, n=>j); /* Store label of it */ if (j>0 && $b(m=>i, n=>j1)>0) neighbour[nfound++] = $b(m=>i, n=>j1); if (j>0 && i>0 && $b(m=>i1, n=>j1)>0) neighbour[nfound++] = $b(m=>i1, n=>j1); if (j>0 && i<(nx-1) && $b(m=>i2, n=>j1)>0) neighbour[nfound++] = $b(m=>i2, n=>j1); if (nfound==0) { /* Assign new label */ $b(m=>i, n=>j) = newlabel++; } else { $b(m=>i, n=>j) = neighbour[0]; if (nfound>1 && pass == 1) { /* Assign equivalents */ for(k=1; ki, n=>j), neighbour[k] ); } } } else { /* No label */ $b(m=>i, n=>j) = 0; } }} /* End of image loop */ } /* Passes */ /* Replace each cycle by single label */ count = 0; for (i = 1; i <= newlabel; i++) if ( i <= equiv[i] ) { count++; this = i; while ( equiv[this] != i ) { next = equiv[this]; equiv[this] = count; this = next; } equiv[this] = count; } /* Now remove equivalences */ for(j=0; ji, n=>j) = equiv[ (PDL_Long) $b(m=>i, n=>j) ] ; }} free(equiv); /* Tidy */ '); pp_addhdr(' #define MAXSEC 32 #define line(x1, x2, y) for (k=x1;k<=x2;k++) \ { /* printf("line from %d to %d\n",x1,x2); */ \ image[k+wx*y] = col; } #define PX(n) ps[2*n] #define PY(n) ps[2*n+1] void polyfill(PDL_Long *image, int wx, int wy, float *ps, int n, PDL_Long col, int *ierr) { int ymin, ymax, xmin, xmax, fwrd = 1, i, j, k, nsect; int x[MAXSEC], temp, l; float s1, s2, t1, t2; ymin = PY(0); ymax = PY(0); xmin = PX(0); xmax = PX(0); *ierr = 0; for (i=1; i PY(i) ? PY(i) : ymin; ymax = ymax < PY(i) ? PY(i) : ymax; xmin = xmin > PX(i) ? PX(i) : xmin; xmax = xmax < PX(i) ? PX(i) : xmax; } if (xmin < 0 || xmax >= wx || ymin < 0 || ymax >= wy) { *ierr = 1; /* clipping */ return; } s1 = PX(n-1); t1 = PY(n-1); for (l=ymin; l<= ymax; l++) { nsect = 0; fwrd = 1; for (i=0; i= l && l > t2)) { if (nsect > MAXSEC) { *ierr = 2; /* too complex */ return; } x[nsect] = (s1+(s2-s1)*((l-t1)/(t2-t1))); nsect += 1; } s1 = s2; t1 = t2; } /* sort the intersections */ for (i=1; i x[i]) { temp = x[j]; x[j] = x[i]; x[i] = temp; } if (fwrd) { for (i=0; i0; i -= 2) line(x[i-1],x[i],l); fwrd = 1; } } } '); pp_def('polyfill', HandleBad => 0, # a marker Pars => 'int [o,nc] im(m,n); float ps(two=2,np); int col()', Code => 'int ierr = 0, nerr; threadloop %{ polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr); ierr = ierr < nerr ? nerr : ierr; %} if (ierr) warn("errors during polygonfilling"); ', Doc => << 'EOD', =for ref fill the area inside the given polygon with a given colour This function works inplace, i.e. modifies C. EOD ); pp_add_exported('', 'polyfillv'); pp_addpm(<<'EOPM'); =head2 polyfillv =for ref return the (dataflown) area of an image within a polygon =for example # increment intensity in area bounded by $poly $im->polyfillv($pol)++; # legal in perl >= 5.6 # compute average intensity within area bounded by $poly $av = $im->polyfillv($poly)->avg; =cut sub PDL::polyfillv { my ($im, $ps) = @_; my $msk = zeroes(long,$im->dims); polyfill($msk, $ps, 1); return $im->where($msk == 1); } *polyfillv = \&PDL::polyfillv; EOPM pp_addhdr('#include "rotate.c"'."\n\n"); pp_add_exported('','rotnewsz'); pp_addxs(' void rotnewsz(m,n,angle) int m int n float angle PPCODE: int newcols, newrows; if (getnewsize(m,n,angle,&newcols,&newrows) != 0) croak("wrong angle (should be between -90 and +90)"); EXTEND(sp,2); PUSHs(sv_2mortal(newSVnv(newcols))); PUSHs(sv_2mortal(newSVnv(newrows))); '); pp_def('rot2d', HandleBad => 0, Pars => 'im(m,n); float angle(); bg(); int aa(); [o] om(p,q)', Code => 'int ierr; if ((ierr = rotate($P(im),$P(om),$SIZE(m),$SIZE(n),$SIZE(p), $SIZE(q),$angle(),$bg(),$aa())) != 0) if (ierr == -1) croak("error during rotate, wrong angle"); else croak("wrong output dims, did you set them?"); ', # ugly workaround since $SIZE(m) and $SIZE(n) are not initialized # when the redodimscode is called # need to fix this! RedoDimsCode => 'int ncols, nrows; if ($PDL(im)->ndims < 2) croak("need > 2d piddle"); if (getnewsize($PDL(im)->dims[0],$PDL(im)->dims[1], $angle(), &ncols, &nrows) != 0) croak("error during rotate, wrong angle"); /* printf("o: %d, p: %d\n",ncols,nrows); */ $SIZE(p) = ncols; $SIZE(q) = nrows;', GenericTypes => ['B'], Doc => << 'EOD', =for ref rotate an image by given C =for example # rotate by 10.5 degrees with antialiasing, set missing values to 7 $rot = $im->rot2d(10.5,7,1); This function rotates an image through an C between -90 and + 90 degrees. Uses/doesn't use antialiasing depending on the C flag. Pixels outside the rotated image are set to C. Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based on "A Fast Algorithm for General Raster Rotation" by Alan Paeth, Graphics Interface '86, pp. 77-81. Use the C function to find out about the dimension of the newly created image ($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle; L offers a more general interface to distortions, including rotation, with various types of sampling; but rot2d is faster. EOD ); pp_def('bilin2d', HandleBad => 0, Pars => 'I(n,m); O(q,p)', Doc=><<'EOD', =for ref Bilinearly maps the first piddle in the second. The interpolated values are actually added to the second piddle which is supposed to be larger than the first one. EOD , Code =>' int i,j,ii,jj,ii1,jj1,num; double x,y,dx,dy,y1,y2,y3,y4,t,u,sum; if ($SIZE(q)>=$SIZE(n) && $SIZE(p)>=$SIZE(m)) { threadloop %{ dx = ((double) ($SIZE(n)-1)) / ($SIZE(q)-1); dy = ((double) ($SIZE(m)-1)) / ($SIZE(p)-1); for(i=0,x=0;i<$SIZE(q);i++,x+=dx) { for(j=0,y=0;j<$SIZE(p);j++,y+=dy) { ii = (int) floor(x); if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2; jj = (int) floor(y); if (jj>=($SIZE(m)-1)) jj = $SIZE(m)-2; ii1 = ii+1; jj1 = jj+1; y1 = $I(n=>ii,m=>jj); y2 = $I(n=>ii1,m=>jj); y3 = $I(n=>ii1,m=>jj1); y4 = $I(n=>ii,m=>jj1); t = x-ii; u = y-jj; $O(q=>i,p=>j) += (1-t)*(1-u)*y1 + t*(1-u)*y2 + t*u*y3 + (1-t)*u*y4; } } %} } else { barf("the second matrix must be greater than first! (bilin2d)"); } '); pp_def('rescale2d', HandleBad => 0, Pars => 'I(m,n); O(p,q)', Doc=><<'EOD', =for ref The first piddle is rescaled to the dimensions of the second (expanding or meaning values as needed) and then added to it in place. Nothing useful is returned. If you want photometric accuracy or automatic FITS header metadata tracking, consider using L instead: it does these things, at some speed penalty compared to rescale2d. EOD , Code =>' int ix,iy,ox,oy,i,j,lx,ly,cx,cy,xx,yy,num; double kx,ky,temp; ix = $SIZE(m); iy = $SIZE(n); ox = $SIZE(p); oy = $SIZE(q); if(ox >= ix && oy >= iy) { threadloop %{ kx = ((double) (ox)) / (ix); ky = ((double) (oy)) / (iy); lx = 0; for(i=0;ixx,q=>yy) += $I(m=>i,n=>j); } ly = cy + 1; } lx = cx + 1; } %} } else if(ox < ix && oy < iy) { threadloop %{ kx = ((double) (ix)) / (ox); ky = ((double) (iy)) / (oy); lx = 0; for(i=0;iyy,m=>xx); num++; } $O(p=>i,q=>j) += temp/num; ly = cy + 1; } lx = cx + 1; } %} } else if(ox >= ix && oy < iy) { threadloop %{ kx = ((double) (ox)) / (ix); ky = ((double) (iy)) / (oy); lx = 0; for(i=0;iyy,m=>i); num++; } for(xx=lx;xx<=cx;xx++) { /* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ $O(p=>xx,q=>j) += temp/num; } ly = cy + 1; } lx = cx + 1; } %} } else if(ox < ix && oy >= iy) { threadloop %{ kx = ((double) (ix)) / (ox); ky = ((double) (oy)) / (iy); lx = 0; for(i=0;ij,m=>xx); num++; } for(yy=ly;yy<=cy;yy++) { /* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ $O(p=>i,q=>yy) += temp/num; } ly = cy + 1; } lx = cx + 1; } %} } else barf("I am not supposed to be here, please report the bug to "); '); # functions to make handling 2D polynomial mappings a bit easier # pp_add_exported('', 'fitwarp2d applywarp2d'); pp_addpm( ' =head2 fitwarp2d =for ref Find the best-fit 2D polynomial to describe a coordinate transformation. =for usage ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf. { options } ) Given a set of points in the output plane (C<$u,$v>), find the best-fit (using singular-value decomposition) 2D polynomial to describe the mapping back to the image plane (C<$x,$y>). The order of the fit is controlled by the C<$nf> parameter (the maximum power of the polynomial is C<$nf - 1>), and you can restrict the terms to fit using the C option. C<$px> and C<$py> are C by C element piddles which describe a polynomial mapping (of order C) from the I C<(u,v)> image to the I C<(x,y)> image: x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j The transformation is returned for the reverse direction (ie output to input image) since that is what is required by the L routine. The L routine can be used to convert a set of C<$u,$v> points given C<$px> and C<$py>. Options: =for options FIT - which terms to fit? default ones(byte,$nf,$nf) THRESH - in svd, remove terms smaller than THRESH * max value default is 1.0e-5 =over 4 =item FIT C allows you to restrict which terms of the polynomial to fit: only those terms for which the FIT piddle evaluates to true will be evaluated. If a 2D piddle is sent in, then it is used for the x and y polynomials; otherwise C<$fit-Eslice(":,:,(0)")> will be used for C<$px> and C<$fit-Eslice(":,:,(1)")> will be used for C<$py>. =item THRESH Remove all singular values whose valus is less than C times the largest singular value. =back The number of points must be at least equal to the number of terms to fit (C<$nf*$nf> points for the default value of C). =for example # points in original image $x = pdl( 0, 0, 100, 100 ); $y = pdl( 0, 100, 100, 0 ); # get warped to these positions $u = pdl( 10, 10, 90, 90 ); $v = pdl( 10, 90, 90, 10 ); # # shift of origin + scale x/y axis only $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] ); ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } ); print "px = ${px}py = $py"; px = [ [-12.5 1.25] [ 0 0] ] py = [ [-12.5 0] [ 1.25 0] ] # # Compared to allowing all 4 terms ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 ); print "px = ${px}py = $py"; px = [ [ -12.5 1.25] [ 1.110223e-16 -1.1275703e-17] ] py = [ [ -12.5 1.6653345e-16] [ 1.25 -5.8546917e-18] ] =head2 applywarp2d =for ref Transform a set of points using a 2-D polynomial mapping =for usage ( $x, $y ) = applywarp2d( $px, $py, $u, $v ) Convert a set of points (stored in 1D piddles C<$u,$v>) to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px> and C<$py>. See L for more information on the format of C<$px> and C<$py>. =cut # use SVD to fit data. Assuming no errors. sub _svd ($$$) { my $basis = shift; my $y = shift; my $thresh = shift; # if we had errors for these points, would normalise the # basis functions, and the output array, by these errors here # perform the SVD my ( $svd_u, $svd_w, $svd_v ) = svd( $basis ); # remove any singular values $svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) ); # perform the back substitution # my $tmp = $y x $svd_u; if ( $PDL::Bad::Status ) { $tmp /= $svd_w->setvaltobad(0.0); $tmp->inplace->setbadtoval(0.0); } else { # not checked my $mask = ($svd_w == 0.0); $tmp /= ( $svd_w + $mask ); $tmp *= ( 1 - $mask ); } my $ans = sumover( $svd_v * $tmp ); return $ans; } # sub: _svd() sub _mkbasis ($$$$) { my $fit = shift; my $npts = shift; my $u = shift; my $v = shift; my $n = $fit->getdim(0) - 1; my $ncoeff = sum( $fit ); my $basis = zeroes( $u->type, $ncoeff, $npts ); my $k = 0; foreach my $j ( 0 .. $n ) { my $tmp_v = $v**$j; foreach my $i ( 0 .. $n ) { if ( $fit->at($i,$j) ) { my $tmp = $basis->slice("($k),:"); $tmp .= $tmp_v * $u**$i; $k++; } } } return $basis; } # sub: _mkbasis() sub PDL::fitwarp2d { croak "Usage: (\$px,\$py) = fitwarp2d(x(m);y(m);u(m);v(m);\$nf; { options })" if $#_ < 4 or ( $#_ >= 5 and ref($_[5]) ne "HASH" ); my $x = shift; my $y = shift; my $u = shift; my $v = shift; my $nf = shift; my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf), THRESH => 1.0e-5 } ); $opts->options( $_[0] ) if $#_ > -1; my $oref = $opts->current(); # safety checks my $npts = $x->nelem; croak "fitwarp2d: x, y, u, and v must be the same size (and 1D)" unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1; my $svd_thresh = $$oref{THRESH}; croak "fitwarp2d: THRESH option must be >= 0." if $svd_thresh < 0; my $fit = $$oref{FIT}; my $fit_ndim = $fit->getndims(); croak "fitwarp2d: FIT option must be sent a (\$nf,\$nf[,2]) element piddle" unless UNIVERSAL::isa($fit,"PDL") and ($fit_ndim == 2 or ($fit_ndim == 3 and $fit->getdim(2) == 2)) and $fit->getdim(0) == $nf and $fit->getdim(1) == $nf; # how many coeffs to fit (first we ensure $fit is either 0 or 1) $fit = convert( $fit != 0, byte ); my ( $fitx, $fity, $ncoeffx, $ncoeffy, $ncoeff ); if ( $fit_ndim == 2 ) { $fitx = $fit; $fity = $fit; $ncoeff = $ncoeffx = $ncoeffy = sum( $fit ); } else { $fitx = $fit->slice(",,(0)"); $fity = $fit->slice(",,(1)"); $ncoeffx = sum($fitx); $ncoeffy = sum($fity); $ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy; } croak "fitwarp2d: number of points must be >= \$ncoeff" unless $npts >= $ncoeff; # create the basis functions for the SVD fitting my ( $basisx, $basisy ); $basisx = _mkbasis( $fitx, $npts, $u, $v ); if ( $fit_ndim == 2 ) { $basisy = $basisx; } else { $basisy = _mkbasis( $fity, $npts, $u, $v ); } my $px = _svd( $basisx, $x, $svd_thresh ); my $py = _svd( $basisy, $y, $svd_thresh ); # convert into $nf x $nf element piddles, if necessary my $nf2 = $nf * $nf; return ( $px->reshape($nf,$nf), $py->reshape($nf,$nf) ) if $ncoeff == $nf2 and $ncoeffx == $ncoeffy; # re-create the matrix my $xtmp = zeroes( $nf, $nf ); my $ytmp = zeroes( $nf, $nf ); my $kx = 0; my $ky = 0; foreach my $i ( 0 .. ($nf - 1) ) { foreach my $j ( 0 .. ($nf - 1) ) { if ( $fitx->at($i,$j) ) { $xtmp->set($i,$j, $px->at($kx) ); $kx++; } if ( $fity->at($i,$j) ) { $ytmp->set($i,$j, $py->at($ky) ); $ky++; } } } return ( $xtmp, $ytmp ) } # sub: fitwarp2d *fitwarp2d = \&PDL::fitwarp2d; sub PDL::applywarp2d { # checks croak "Usage: (\$x,\$y) = applywarp2d(px(nf,nf);py(nf,nf);u(m);v(m);)" if $#_ != 3; my $px = shift; my $py = shift; my $u = shift; my $v = shift; my $npts = $u->nelem; # safety check croak "applywarp2d: u and v must be the same size (and 1D)" unless $npts == $u->nelem and $npts == $v->nelem and $u->getndims == 1 and $v->getndims == 1; my $nf = $px->getdim(0); my $nf2 = $nf * $nf; # could remove terms with 0 coeff here # (would also have to remove them from px/py for # the matrix multiplication below) # my $mat = _mkbasis( ones(byte,$nf,$nf), $npts, $u, $v ); my $x = reshape( $mat x $px->clump(-1)->transpose(), $npts ); my $y = reshape( $mat x $py->clump(-1)->transpose(), $npts ); return ( $x, $y ); } # sub: applywarp2d *applywarp2d = \&PDL::applywarp2d; ' ); ## resampling routines taken from v3.6-0 of the Eclipse package ## http://www.eso.org/eclipse by Nicolas Devillard ## pp_addhdr( '#include "resample.h"' . "\n" ); # pod for warp2d # and support routine # pp_addpm( <<'EOD'); =head2 warp2d =for sig Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { options }) =for ref Warp a 2D image given a polynomial describing the I mapping. =for usage $out = warp2d( $img, $px, $py, { options } ); Apply the polynomial transformation encoded in the C<$px> and C<$py> piddles to warp the input image C<$img> into the output image C<$out>. The format for the polynomial transformation is described in the documentation for the L routine. At each point C, the closest 16 pixel values are combined with an interpolation kernel to calculate the value at C. The interpolation is therefore done in the image, rather than Fourier, domain. By default, a C kernel is used, but this can be changed using the C option discussed below (the choice of kernel depends on the frequency content of the input image). The routine is based on the C command from the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and for further details on image resampling see Wolberg, G., "Digital Image Warping", 1990, IEEE Computer Society Press ISBN 0-8186-8944-7). Currently the output image is the same size as the input one, which means data will be lost if the transformation reduces the pixel scale. This will (hopefully) be changed soon. =for example $img = rvals(byte,501,501); imag $img, { JUSTIFY => 1 }; # # use a not-particularly-obvious transformation: # x = -10 + 0.5 * $u - 0.1 * $v # y = -20 + $v - 0.002 * $u * $v # $px = pdl( [ -10, 0.5 ], [ -0.1, 0 ] ); $py = pdl( [ -20, 0 ], [ 1, 0.002 ] ); $wrp = warp2d( $img, $px, $py ); # # see the warped image imag $warp, { JUSTIFY => 1 }; The options are: =for options KERNEL - default value is tanh NOVAL - default value is 0 C is used to specify which interpolation kernel to use (to see what these kernels look like, use the L routine). The options are: =over 4 =item tanh Hyperbolic tangent: the approximation of an ideal box filter by the product of symmetric tanh functions. =item sinc For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle, which produces a C interpolation kernel in the spatial domain: sinc(x) = sin(pi * x) / (pi * x) However, it is not ideal for the C<4x4> pixel region used here. =item sinc2 This is the square of the sinc function. =item lanczos Although defined differently to the C kernel, the result is very similar in the spatial domain. The Lanczos function is defined as L(x) = sinc(x) * sinc(x/2) if abs(x) < 2 = 0 otherwise =item hann This kernel is derived from the following function: H(x) = a + (1-a) * cos(2*pi*x/(N-1)) if abs(x) < 0.5*(N-1) = 0 otherwise with C and N currently equal to 2001. =item hamming This kernel uses the same C as the Hann filter, but with C. =back C gives the value used to indicate that a pixel in the output image does not map onto one in the input image. =cut # support routine { my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann ); # note: convert to lower case sub _check_kernel ($$) { my $kernel = lc shift; my $code = shift; barf "Unknown kernel $kernel sent to $code\n" . "\tmust be one of [" . join(',',keys %warp2d) . "]\n" unless exists $warp2d{$kernel}; return $kernel; } } EOD pp_def( 'warp2d', Doc=> undef, HandleBad => 0, Pars => 'img(m,n); double px(np,np); double py(np,np); [o] warp(m,n);', OtherPars => 'char *kernel_type; double noval;', GenericTypes => [ 'F', 'D' ], PMCode => ' sub PDL::warp2d { my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } ); $opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH"; die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} )" if $#_<2 || $#_>3; my $img = shift; my $px = shift; my $py = shift; my $out = $#_ == -1 ? PDL->null() : shift; # safety checks my $copt = $opts->current(); my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" ); &PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} ); return $out; } ', Code => ' int i, j, k ; int ncoeff, lx_out, ly_out ; int lx_3, ly_3 ; double cur ; double neighbors[16] ; double rsc[8], sumrs ; double x, y ; int px, py ; int tabx, taby ; double *kernel, *poly ; int da[16], db[16] ; /* Generate default interpolation kernel */ kernel = generate_interpolation_kernel( $COMP(kernel_type) ) ; if (kernel == NULL) { croak( "Ran out of memory building kernel\n" ); } /* Compute sizes */ ncoeff = $SIZE(np); lx_out = $SIZE(m); /* is this right? */ ly_out = $SIZE(n); lx_3 = lx_out - 3; ly_3 = ly_out - 3; /* Pre compute leaps for 16 closest neighbors positions */ da[0] = -1; db[0] = -1; da[1] = 0; db[1] = -1; da[2] = 1; db[2] = -1; da[3] = 2; db[3] = -1; da[4] = -1; db[4] = 0; da[5] = 0; db[5] = 0; da[6] = 1; db[6] = 0; da[7] = 2; db[7] = 0; da[8] = -1; db[8] = 1; da[9] = 0; db[9] = 1; da[10] = 1; db[10] = 1; da[11] = 2; db[11] = 1; da[12] = -1; db[12] = 2; da[13] = 0; db[13] = 2; da[14] = 1; db[14] = 2; da[15] = 2; db[15] = 2; /* allocate memory for polynomial */ poly = malloc( ncoeff * sizeof(double) ); if ( poly == NULL ) { croak( "Ran out of memory\n" ); } poly[0] = 1.0; /* Loop over the output image */ threadloop %{ loop(n) %{ /* fill in poly array */ for ( k = 1; k < ncoeff; k++ ) { poly[k] = (double) n * poly[k-1]; } loop(m) %{ /* Compute the original source for this pixel */ x = poly2d_compute( ncoeff, $P(px), (double) m, poly ); y = poly2d_compute( ncoeff, $P(py), (double) m, poly ); /* Which is the closest integer positioned neighbor? */ px = (int)x ; py = (int)y ; if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3)) $warp() = ($GENERIC()) $COMP(noval); else { /* Now feed the positions for the closest 16 neighbors */ for (k=0 ; k<16 ; k++) { i = px + da[k]; j = py + db[k]; neighbors[k] = (double) $img( m => i, n => j ); } /* Which tabulated value index shall we use? */ tabx = (x - (double)px) * (double)(TABSPERPIX) ; taby = (y - (double)py) * (double)(TABSPERPIX) ; /* Compute resampling coefficients */ /* rsc[0..3] in x, rsc[4..7] in y */ rsc[0] = kernel[TABSPERPIX + tabx] ; rsc[1] = kernel[tabx] ; rsc[2] = kernel[TABSPERPIX - tabx] ; rsc[3] = kernel[2 * TABSPERPIX - tabx] ; rsc[4] = kernel[TABSPERPIX + taby] ; rsc[5] = kernel[taby] ; rsc[6] = kernel[TABSPERPIX - taby] ; rsc[7] = kernel[2 * TABSPERPIX - taby] ; sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) * (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ; /* Compute interpolated pixel now */ cur = rsc[4] * ( rsc[0]*neighbors[0] + rsc[1]*neighbors[1] + rsc[2]*neighbors[2] + rsc[3]*neighbors[3] ) + rsc[5] * ( rsc[0]*neighbors[4] + rsc[1]*neighbors[5] + rsc[2]*neighbors[6] + rsc[3]*neighbors[7] ) + rsc[6] * ( rsc[0]*neighbors[8] + rsc[1]*neighbors[9] + rsc[2]*neighbors[10] + rsc[3]*neighbors[11] ) + rsc[7] * ( rsc[0]*neighbors[12] + rsc[1]*neighbors[13] + rsc[2]*neighbors[14] + rsc[3]*neighbors[15] ) ; /* Copy the value to the output image */ $warp() = ($GENERIC()) (cur/sumrs); } /* if: edge or interior */ %} /* loop(m) */ %} /* loop(n) */ %} /* threadloop */ free(poly); free(kernel) ; ', ); # pp_def: warp2d pp_addxs( ' int _get_kernel_size() PROTOTYPE: CODE: RETVAL = KERNEL_SAMPLES; OUTPUT: RETVAL '); pp_add_exported('', 'warp2d_kernel'); pp_addpm( ' =head2 warp2d_kernel =for ref Return the specified kernel, as used by L =for usage ( $x, $k ) = warp2d_kernel( $name ) The valid values for C<$name> are the same as the C option of L. =for example line warp2d_kernel( "hamming" ); =cut '); # pp_addpm # this is not very clever, but it's a pain to create a valid # piddle in XS code # pp_def( 'warp2d_kernel', Doc => undef, HandleBad => 0, PMCode => ' sub PDL::warp2d_kernel ($) { my $kernel = _check_kernel( shift, "warp2d_kernel" ); my $nelem = _get_kernel_size(); my $x = zeroes( $nelem ); my $k = zeroes( $nelem ); &PDL::_warp2d_kernel_int( $x, $k, $kernel ); return ( $x, $k ); # return _get_kernel( $kernel ); } *warp2d_kernel = \&PDL::warp2d_kernel; ', Pars => '[o] x(n); [o] k(n);', OtherPars => 'char *name;', GenericTypes => [ 'D' ], Code => ' double *kernel, xx; if ( $SIZE(n) != KERNEL_SAMPLES ) { croak( "Internal error in warp2d_kernel - mismatch in kernel size\n" ); } kernel = generate_interpolation_kernel($COMP(name)); if ( kernel == NULL ) { croak( "unable to allocate memory for kernel" ); } /* fill in piddles */ xx = 0.0; threadloop %{ loop (n) %{ $x() = xx; $k() = kernel[n]; xx += 1.0 / (double) TABSPERPIX; %} %} /* free the kernel */ free( kernel ); '); # pp_addpm pp_done();