"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Libtmp/Image2D/image2d.pd" between
PDL-2.080.tar.gz and PDL-2.081.tar.gz

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

image2d.pd  (PDL-2.080)	:	image2d.pd  (PDL-2.081)
use strict;		use strict;
use warnings;		use warnings;
		use PDL::Types qw(ppdefs_all);
		my $A = [ppdefs_all];
pp_addpm({At=>'Top'},<<'EOD');		pp_addpm({At=>'Top'},<<'EOD');
use strict;		use strict;
use warnings;		use warnings;
=head1 NAME		=head1 NAME
PDL::Image2D - Miscellaneous 2D image processing functions		PDL::Image2D - Miscellaneous 2D image processing functions
=head1 DESCRIPTION		=head1 DESCRIPTION
skipping to change at line 265		skipping to change at line 266
my($x,$kern) = @_;		my($x,$kern) = @_;
my $c = $#_ == 2 ? $_[2] : $x->nullcreate;		my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
PDL::_conv2d_int($x,$kern,$c,		PDL::_conv2d_int($x,$kern,$c,
(!(defined $opt && exists $$opt{Boundary}))?0:		(!(defined $opt && exists $$opt{Boundary}))?0:
(($$opt{Boundary} eq "Reflect") +		(($$opt{Boundary} eq "Reflect") +
2*($$opt{Boundary} eq "Truncate") +		2*($$opt{Boundary} eq "Truncate") +
3*($$opt{Boundary} eq "Replicate")));		3*($$opt{Boundary} eq "Replicate")));
return $c;		return $c;
}		}
',		',
		GenericTypes => $A,
Code =>		Code =>
init_vars( { vars => 'PDL_Double tmp; int flag;' } ) .		init_vars( { vars => 'PDL_CLDouble tmp; int flag;' } ) .
init_map("i") .		init_map("i") .
init_map("j") .		init_map("j") .
'		'
broadcastloop %{		broadcastloop %{
loop(n) %{		loop(n) %{
loop(m) %{		loop(m) %{
tmp = 0;		tmp = 0;
for(j1=0; j1<q_size; j1++) {		for(j1=0; j1<q_size; j1++) {
j2 = mapj[n-j1];		j2 = mapj[n-j1];
if (j2 >= 0) {		if (j2 >= 0) {
skipping to change at line 332		skipping to change at line 334
=> Reflect - reflect at boundary		=> Reflect - reflect at boundary
=> Truncate - truncate at boundary		=> Truncate - truncate at boundary
=> Replicate - repeat boundary pixel values		=> Replicate - repeat boundary pixel values
=cut		=cut
EOD		EOD
BadDoc =>		BadDoc =>
'Bad values are ignored in the calculation. If all elements within the		'Bad values are ignored in the calculation. If all elements within the
kernel are bad, the output is set bad.',		kernel are bad, the output is set bad.',
HandleBad => 1,		HandleBad => 1,
Pars => 'a(m,n); kern(p,q); [o]b(m,n); double [t]tmp(pq);',		Pars => 'a(m,n); kern(p,q); [o]b(m,n); double+ [t]tmp(pq);',
OtherPars => 'int opt;',		OtherPars => 'int opt;',
PMCode => <<'EOF',		PMCode => <<'EOF',
sub PDL::med2d {		sub PDL::med2d {
my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';		my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
die 'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'		die 'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
if $#_<1 || $#_>2;		if $#_<1 || $#_>2;
my($x,$kern) = @_;		my($x,$kern) = @_;
croak "med2d: kernel must contain some positive elements.\n"		croak "med2d: kernel must contain some positive elements.\n"
if all( $kern <= 0 );		if all( $kern <= 0 );
my $c = $#_ == 2 ? $_[2] : $x->nullcreate;		my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
PDL::_med2d_int($x,$kern,$c,		PDL::_med2d_int($x,$kern,$c,
(!(defined $opt && exists $$opt{Boundary}))?0:		(!(defined $opt && exists $$opt{Boundary}))?0:
(($$opt{Boundary} eq "Reflect") +		(($$opt{Boundary} eq "Reflect") +
2*($$opt{Boundary} eq "Truncate") +		2*($$opt{Boundary} eq "Truncate") +
3*($$opt{Boundary} eq "Replicate")));		3*($$opt{Boundary} eq "Replicate")));
return $c;		return $c;
}		}
EOF		EOF
RedoDimsCode => '$SIZE(pq) = $SIZE(p)*$SIZE(q);',		RedoDimsCode => '$SIZE(pq) = $SIZE(p)*$SIZE(q);',
Code =>		Code =>
init_vars( { vars => 'PDL_Double kk, aa; int flag;' } ) .		init_vars( { vars => 'int flag;' } ) .
init_map("i") .		init_map("i") .
init_map("j") .		init_map("j") .
'		'
		PDL_LDouble kk, aa;
broadcastloop %{		broadcastloop %{
for(j=0; j<n_size; j++) {		for(j=0; j<n_size; j++) {
for(i=0; i<m_size; i++) {		for(i=0; i<m_size; i++) {
PDL_Indx count = 0;		PDL_Indx count = 0;
flag = 0;		flag = 0;
for(j1=0; j1<q_size; j1++) {		for(j1=0; j1<q_size; j1++) {
j2 = mapj[j-j1];		j2 = mapj[j-j1];
if (j2 >= 0)		if (j2 >= 0)
for(i1=0; i1<p_size; i1++) {		for(i1=0; i1<p_size; i1++) {
i2 = mapi[i-i1];		i2 = mapi[i-i1];
skipping to change at line 482		skipping to change at line 485
',		',
); # pp_def: med2df		); # pp_def: med2df
pp_addhdr(<<'EOH');		pp_addhdr(<<'EOH');
#define EZ(x) ez ? 0 : (x)		#define EZ(x) ez ? 0 : (x)
EOH		EOH
pp_def('box2d',		pp_def('box2d',
Pars => 'a(n,m); [o] b(n,m)',		Pars => 'a(n,m); [o] b(n,m)',
OtherPars => 'int wx; int wy; int edgezero',		OtherPars => 'int wx; int wy; int edgezero',
		GenericTypes => $A,
Code => '		Code => '
register int nx = 0.5*$COMP(wx);		register int nx = 0.5*$COMP(wx);
register int ny = 0.5*$COMP(wy);		register int ny = 0.5*$COMP(wy);
register int xs = $SIZE(n);		register int xs = $SIZE(n);
register int ys = $SIZE(m);		register int ys = $SIZE(m);
register int ez = $COMP(edgezero);		register int ez = $COMP(edgezero);
double div, sum, lsum;		PDL_CLDouble div, sum, lsum;
int xx,yy,y,ind1,ind2,first;		int xx,yy,y,ind1,ind2,first;
div = 1/((2.0*nx+1)*(2.0*ny+1));		div = 1/((2.0*nx+1)*(2.0*ny+1));
broadcastloop %{		broadcastloop %{
first = 1;		first = 1;
for (y=0;y<ys;y++)		for (y=0;y<ys;y++)
for (xx=0; xx<nx; xx++) {		for (xx=0; xx<nx; xx++) {
ind1 = xs-1-xx; /* rightmost strip */		ind1 = xs-1-xx; /* rightmost strip */
$b(n=>xx,m=>y) = EZ($a(n=>xx,m=>y));		$b(n=>xx,m=>y) = EZ($a(n=>xx,m=>y));
skipping to change at line 588		skipping to change at line 592
if all neighbours are bad, the original data value is		if all neighbours are bad, the original data value is
copied across.		copied across.
=cut		=cut
EOD		EOD
BadDoc =>		BadDoc =>
'This routine does not handle bad values - use L</patchbad2d> instead',		'This routine does not handle bad values - use L</patchbad2d> instead',
HandleBad => 0,		HandleBad => 0,
Pars => 'a(m,n); int bad(m,n); [o]b(m,n);',		Pars => 'a(m,n); int bad(m,n); [o]b(m,n);',
		GenericTypes => $A,
Code =>		Code =>
'int m_size, n_size, i,j, i1,j1, i2,j2, norm;		'PDL_Indx m_size, n_size, i,j, i1,j1, i2,j2, norm;
double tmp;		PDL_CLDouble tmp;
m_size = $SIZE(m); n_size = $SIZE(n);		m_size = $SIZE(m); n_size = $SIZE(n);
broadcastloop %{		broadcastloop %{
for(j=0; j<n_size; j++) {		for(j=0; j<n_size; j++) {
for(i=0; i<m_size; i++) {		for(i=0; i<m_size; i++) {
$b(m=>i,n=>j) = $a(m=>i,n=>j);		$b(m=>i,n=>j) = $a(m=>i,n=>j);
skipping to change at line 652		skipping to change at line 657
Pixels are replaced by the average of their non-bad neighbours;		Pixels are replaced by the average of their non-bad neighbours;
if all neighbours are bad, the output is set bad.		if all neighbours are bad, the output is set bad.
If the input ndarray contains I<no> bad values, then a straight copy		If the input ndarray contains I<no> bad values, then a straight copy
is performed (see L</patch2d>).		is performed (see L</patch2d>).
EOD		EOD
BadDoc =>		BadDoc =>
'patchbad2d handles bad values. The output ndarray I<may> contain		'patchbad2d handles bad values. The output ndarray I<may> contain
bad values, depending on the pattern of bad values in the input ndarray.',		bad values, depending on the pattern of bad values in the input ndarray.',
HandleBad => 1,		HandleBad => 1,
Pars => 'a(m,n); [o]b(m,n);',		Pars => 'a(m,n); [o]b(m,n);',
		GenericTypes => $A,
Code =>		Code =>
'int i1,j1, flag = 0;		'int i1,j1, flag = 0;
PDL_Indx m_size = $SIZE(m), n_size = $SIZE(n);		PDL_Indx m_size = $SIZE(m), n_size = $SIZE(n);
broadcastloop %{		broadcastloop %{
loop(n) %{		loop(n) %{
loop(m) %{		loop(m) %{
$GENERIC(a) a_val = $a();		$GENERIC(a) a_val = $a();
PDL_IF_BAD(if ( !$ISGOODVAR(a_val,a) ) {		PDL_IF_BAD(if ( !$ISGOODVAR(a_val,a) ) {
double tmp = 0; PDL_Indx norm=0;		PDL_CLDouble tmp = 0; PDL_Indx norm=0;
for(j1=-1; j1<=1; j1++) {		for(j1=-1; j1<=1; j1++) {
PDL_Indx j2 = n+j1;		PDL_Indx j2 = n+j1;
if ( j2>=0 && j2<n_size ) {		if ( j2>=0 && j2<n_size ) {
for(i1=-1; i1<=1; i1++) {		for(i1=-1; i1<=1; i1++) {
/* ignore central pixel, which we know is bad */		/* ignore central pixel, which we know is bad */
if ( i1!=0 || j1!=0 ) {		if ( i1!=0 || j1!=0 ) {
PDL_Indx i2 = m+i1;		PDL_Indx i2 = m+i1;
if ( i2>=0 && i2<m_size ) {		if ( i2>=0 && i2<m_size ) {
a_val = $a(m=>i2,n=>j2);		a_val = $a(m=>i2,n=>j2);
if ( $ISGOODVAR(a_val,a) ) {		if ( $ISGOODVAR(a_val,a) ) {
skipping to change at line 720		skipping to change at line 726
BadDoc=><<'EOD',		BadDoc=><<'EOD',
Bad values are excluded from the search. If all pixels		Bad values are excluded from the search. If all pixels
are bad then the output is set bad.		are bad then the output is set bad.
EOD		EOD
HandleBad => 1,		HandleBad => 1,
Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();',		Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();',
Code => '		Code => '
double cur = 0;		PDL_LDouble cur = 0;
PDL_Indx curind1 = PDL_IF_BAD(-1,0), curind2 = PDL_IF_BAD(-1,0);		PDL_Indx curind1 = PDL_IF_BAD(-1,0), curind2 = PDL_IF_BAD(-1,0);
loop(m) %{		loop(m) %{
loop(n) %{		loop(n) %{
if( PDL_IF_BAD($ISGOOD(a()) &&,) ( (!n && !m) || $a() > cur || IsNaN( cur)) ) {		if( PDL_IF_BAD($ISGOOD(a()) &&,) ( (!n && !m) || $a() > cur || IsNaN( cur)) ) {
cur = $a(); curind1 = m; curind2 = n;		cur = $a(); curind1 = m; curind2 = n;
}		}
%}		%}
%}		%}
PDL_IF_BAD(if ( curind1 < 0 ) {		PDL_IF_BAD(if ( curind1 < 0 ) {
$SETBAD(val());		$SETBAD(val());
skipping to change at line 760		skipping to change at line 766
EOD		EOD
BadDoc=><<'EOD',		BadDoc=><<'EOD',
Bad pixels are excluded from the centroid calculation. If all elements are		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		bad (or the pixel sum is 0 - but why would you be centroiding
something with negatives in...) then the output values are set bad.		something with negatives in...) then the output values are set bad.
EOD		EOD
HandleBad => 1,		HandleBad => 1,
Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();',		Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();',
Code => '		Code => '
PDL_Indx i,j,i1,i2,j1,j2,m_size = $SIZE(m),n_size = $SIZE(n);		PDL_Indx i,j,i1,i2,j1,j2,m_size = $SIZE(m),n_size = $SIZE(n);
double sum,data,sumx,sumy;		PDL_LDouble sum,data,sumx,sumy;
i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;		i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;
i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;		i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;
j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;		j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;		j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;
sum = sumx = sumy = 0;		sum = sumx = sumy = 0;
for(j=j1; j<=j2; j++) {		for(j=j1; j<=j2; j++) {
for(i=i1; i<=i2; i++) {		for(i=i1; i<=i2; i++) {
data = $im(m=>i,n=>j);		data = $im(m=>i,n=>j);
PDL_IF_BAD(if ( $ISGOODVAR(data,im) ),) {		PDL_IF_BAD(if ( $ISGOODVAR(data,im) ),) {
sum += data;		sum += data;
skipping to change at line 1216		skipping to change at line 1222
if (getnewsize(m,n,angle,&newcols,&newrows) != 0)		if (getnewsize(m,n,angle,&newcols,&newrows) != 0)
croak("wrong angle (should be between -90 and +90)");		croak("wrong angle (should be between -90 and +90)");
EXTEND(sp,2);		EXTEND(sp,2);
PUSHs(sv_2mortal(newSVnv(newcols)));		PUSHs(sv_2mortal(newSVnv(newcols)));
PUSHs(sv_2mortal(newSVnv(newrows)));		PUSHs(sv_2mortal(newSVnv(newrows)));
');		');
pp_def('rot2d',		pp_def('rot2d',
HandleBad => 0,		HandleBad => 0,
Pars => 'im(m,n); float angle(); bg(); int aa(); [o] om(p,q)',		Pars => 'im(m,n); float angle(); bg(); int aa(); [o] om(p,q)',
		GenericTypes => $A,
Code => 'int ierr;		Code => 'int ierr;
if ((ierr = rotate($P(im),$P(om),$SIZE(m),$SIZE(n),$SIZE(p),		if ((ierr = rotate($P(im),$P(om),$SIZE(m),$SIZE(n),$SIZE(p),
$SIZE(q),$angle(),$bg(),$aa())) != 0) {		$SIZE(q),$angle(),$bg(),$aa())) != 0) {
if (ierr == -1)		if (ierr == -1)
$CROAK("error during rotate, wrong angle");		$CROAK("error during rotate, wrong angle");
else		else
$CROAK("wrong output dims, did you set them?");		$CROAK("wrong output dims, did you set them?");
}',		}',
RedoDimsCode => '		RedoDimsCode => '
if (getnewsize($SIZE(m),$SIZE(n),		if (getnewsize($SIZE(m),$SIZE(n),
skipping to change at line 1280		skipping to change at line 1287
=for ref		=for ref
Bilinearly maps the first ndarray in the second. The		Bilinearly maps the first ndarray in the second. The
interpolated values are actually added to the second		interpolated values are actually added to the second
ndarray which is supposed to be larger than the first one.		ndarray which is supposed to be larger than the first one.
=cut		=cut
EOD		EOD
,		,
		GenericTypes => $A,
Code =>'		Code =>'
int i,j,ii,jj,ii1,jj1,num;		PDL_Indx i,j,ii,jj,ii1,jj1,num;
double x,y,dx,dy,y1,y2,y3,y4,t,u,sum;		PDL_CLDouble x,y,dx,dy,y1,y2,y3,y4,t,u,sum;
if ($SIZE(q)>=$SIZE(n) && $SIZE(p)>=$SIZE(m)) {		if ($SIZE(q)>=$SIZE(n) && $SIZE(p)>=$SIZE(m)) {
broadcastloop %{		broadcastloop %{
dx = ((double) ($SIZE(n)-1)) / ($SIZE(q)-1);		dx = ((double) ($SIZE(n)-1)) / ($SIZE(q)-1);
dy = ((double) ($SIZE(m)-1)) / ($SIZE(p)-1);		dy = ((double) ($SIZE(m)-1)) / ($SIZE(p)-1);
for(i=0,x=0;i<$SIZE(q);i++,x+=dx) {		for(i=0,x=0;i<$SIZE(q);i++,x+=dx) {
for(j=0,y=0;j<$SIZE(p);j++,y+=dy) {		for(j=0,y=0;j<$SIZE(p);j++,y+=dy) {
ii = (int) floor(x);		ii = (int) floor(x);
if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2;		if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2;
jj = (int) floor(y);		jj = (int) floor(y);
skipping to change at line 1332		skipping to change at line 1340
If you want photometric accuracy or automatic FITS header metadata		If you want photometric accuracy or automatic FITS header metadata
tracking, consider using L<PDL::Transform::map|PDL::Transform/map>		tracking, consider using L<PDL::Transform::map|PDL::Transform/map>
instead: it does these things, at some speed penalty compared to		instead: it does these things, at some speed penalty compared to
rescale2d.		rescale2d.
=cut		=cut
EOD		EOD
,		,
		GenericTypes => $A,
Code =>'		Code =>'
int ix,iy,ox,oy,i,j,lx,ly,cx,cy,xx,yy,num;		PDL_Indx ix,iy,ox,oy,i,j,lx,ly,cx,cy,xx,yy,num;
double kx,ky,temp;		PDL_CLDouble kx,ky,temp;
ix = $SIZE(m);		ix = $SIZE(m);
iy = $SIZE(n);		iy = $SIZE(n);
ox = $SIZE(p);		ox = $SIZE(p);
oy = $SIZE(q);		oy = $SIZE(q);
if(ox >= ix && oy >= iy) {		if(ox >= ix && oy >= iy) {
broadcastloop %{		broadcastloop %{
kx = ((double) (ox)) / (ix);		kx = ((double) (ox)) / (ix);
ky = ((double) (oy)) / (iy);		ky = ((double) (oy)) / (iy);
skipping to change at line 2022		skipping to change at line 2031
return $kernel;		return $kernel;
}		}
}		}
EOD		EOD
pp_def(		pp_def(
'warp2d',		'warp2d',
Doc=> undef,		Doc=> undef,
HandleBad => 0,		HandleBad => 0,
Pars => 'img(m,n); double px(np,np); double py(np,np); [o] warp(m,n);',		Pars => 'img(m,n); ldouble px(np,np); ldouble py(np,np); [o] warp(m,n);',
OtherPars => 'char *kernel_type; double noval;',		OtherPars => 'char *kernel_type; double noval;',
GenericTypes => [ 'F', 'D' ],		GenericTypes => [ qw(F D E) ],
PMCode => '		PMCode => '
sub PDL::warp2d {		sub PDL::warp2d {
my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );		my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );
$opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";		$opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";
die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} ) "		die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} ) "
if $#_<2 || $#_>3;		if $#_<2 || $#_>3;
my $img = shift;		my $img = shift;
my $px = shift;		my $px = shift;
skipping to change at line 2049		skipping to change at line 2058
my $copt = $opts->current();		my $copt = $opts->current();
my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );		my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );
&PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} );		&PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} );
return $out;		return $out;
}		}
',		',
Code => '		Code => '
int i, j, k ;		PDL_Indx ncoeff, lx_out, ly_out, i, j, k, lx_3, ly_3, px, py, tabx, taby;
PDL_Indx ncoeff, lx_out, ly_out ;		PDL_Indx da[16], db[16] ;
PDL_Indx lx_3, ly_3 ;		PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel, *poly ;
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 */		/* Generate default interpolation kernel */
kernel = generate_interpolation_kernel( $COMP(kernel_type) ) ;		kernel = generate_interpolation_kernel( $COMP(kernel_type) ) ;
if (kernel == NULL) {		if (kernel == NULL) {
$CROAK( "Ran out of memory building kernel\n" );		$CROAK( "Ran out of memory building kernel\n" );
}		}
/* Compute sizes */		/* Compute sizes */
ncoeff = $SIZE(np);		ncoeff = $SIZE(np);
lx_out = $SIZE(m); /* is this right? */		lx_out = $SIZE(m); /* is this right? */
skipping to change at line 2097		skipping to change at line 2097
da[9] = 0; db[9] = 1;		da[9] = 0; db[9] = 1;
da[10] = 1; db[10] = 1;		da[10] = 1; db[10] = 1;
da[11] = 2; db[11] = 1;		da[11] = 2; db[11] = 1;
da[12] = -1; db[12] = 2;		da[12] = -1; db[12] = 2;
da[13] = 0; db[13] = 2;		da[13] = 0; db[13] = 2;
da[14] = 1; db[14] = 2;		da[14] = 1; db[14] = 2;
da[15] = 2; db[15] = 2;		da[15] = 2; db[15] = 2;
/* allocate memory for polynomial */		/* allocate memory for polynomial */
poly = malloc( ncoeff * sizeof(double) );		poly = malloc( ncoeff * sizeof(PDL_LDouble) );
if ( poly == NULL ) {		if ( poly == NULL ) {
$CROAK( "Ran out of memory\n" );		$CROAK( "Ran out of memory\n" );
}		}
poly[0] = 1.0;		poly[0] = 1.0;
/* Loop over the output image */		/* Loop over the output image */
broadcastloop %{		broadcastloop %{
loop(n) %{		loop(n) %{
/* fill in poly array */		/* fill in poly array */
for ( k = 1; k < ncoeff; k++ ) {		for ( k = 1; k < ncoeff; k++ ) {
poly[k] = (double) n * poly[k-1];		poly[k] = (PDL_LDouble) n * poly[k-1];
}		}
loop(m) %{		loop(m) %{
/* Compute the original source for this pixel */		/* Compute the original source for this pixel */
x = poly2d_compute( ncoeff, $P(px), (double) m, poly );		x = poly2d_compute( ncoeff, $P(px), (PDL_LDouble) m, poly );
y = poly2d_compute( ncoeff, $P(py), (double) m, poly );		y = poly2d_compute( ncoeff, $P(py), (PDL_LDouble) m, poly );
/* Which is the closest integer positioned neighbor? */		/* Which is the closest integer positioned neighbor? */
px = (int)x ;		px = (int)x ;
py = (int)y ;		py = (int)y ;
if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))		if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
$warp() = ($GENERIC()) $COMP(noval);		$warp() = ($GENERIC()) $COMP(noval);
else {		else {
/* Now feed the positions for the closest 16 neighbors */		/* Now feed the positions for the closest 16 neighbors */
for (k=0 ; k<16 ; k++) {		for (k=0 ; k<16 ; k++) {
i = px + da[k]; j = py + db[k];		i = px + da[k]; j = py + db[k];
neighbors[k] = (double) $img( m => i, n => j );		neighbors[k] = (PDL_LDouble) $img( m => i, n => j );
}		}
/* Which tabulated value index shall we use? */		/* Which tabulated value index shall we use? */
tabx = (x - (double)px) * (double)(TABSPERPIX) ;		tabx = (x - (PDL_LDouble)px) * (PDL_LDouble)(TABSPERPIX) ;
taby = (y - (double)py) * (double)(TABSPERPIX) ;		taby = (y - (PDL_LDouble)py) * (PDL_LDouble)(TABSPERPIX) ;
/* Compute resampling coefficients */		/* Compute resampling coefficients */
/* rsc[0..3] in x, rsc[4..7] in y */		/* rsc[0..3] in x, rsc[4..7] in y */
rsc[0] = kernel[TABSPERPIX + tabx] ;		rsc[0] = kernel[TABSPERPIX + tabx] ;
rsc[1] = kernel[tabx] ;		rsc[1] = kernel[tabx] ;
rsc[2] = kernel[TABSPERPIX - tabx] ;		rsc[2] = kernel[TABSPERPIX - tabx] ;
rsc[3] = kernel[2 * TABSPERPIX - tabx] ;		rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
rsc[4] = kernel[TABSPERPIX + taby] ;		rsc[4] = kernel[TABSPERPIX + taby] ;
rsc[5] = kernel[taby] ;		rsc[5] = kernel[taby] ;
skipping to change at line 2245		skipping to change at line 2245
&PDL::_warp2d_kernel_int( $x, $k, $kernel );		&PDL::_warp2d_kernel_int( $x, $k, $kernel );
return ( $x, $k );		return ( $x, $k );
# return _get_kernel( $kernel );		# return _get_kernel( $kernel );
}		}
*warp2d_kernel = \&PDL::warp2d_kernel;		*warp2d_kernel = \&PDL::warp2d_kernel;
',		',
Pars => '[o] x(n); [o] k(n);',		Pars => '[o] x(n); [o] k(n);',
OtherPars => 'char *name;',		OtherPars => 'char *name;',
GenericTypes => [ 'D' ],		GenericTypes => [ qw(F D E) ],
Code => '		Code => '
double *kernel, xx;		PDL_LDouble *kernel, xx;
if ( $SIZE(n) != KERNEL_SAMPLES ) {		if ( $SIZE(n) != KERNEL_SAMPLES ) {
$CROAK( "Internal error in warp2d_kernel - mismatch in kernel size\n" );		$CROAK( "Internal error in warp2d_kernel - mismatch in kernel size\n" );
}		}
kernel = generate_interpolation_kernel($COMP(name));		kernel = generate_interpolation_kernel($COMP(name));
if ( kernel == NULL ) {		if ( kernel == NULL ) {
$CROAK( "unable to allocate memory for kernel" );		$CROAK( "unable to allocate memory for kernel" );
}		}
skipping to change at line 2272		skipping to change at line 2272
loop (n) %{		loop (n) %{
$x() = xx;		$x() = xx;
$k() = kernel[n];		$k() = kernel[n];
xx += 1.0 / (double) TABSPERPIX;		xx += 1.0 / (double) TABSPERPIX;
%}		%}
%}		%}
/* free the kernel */		/* free the kernel */
kernel_free( kernel );		kernel_free( kernel );
'); # pp_addpm		');
pp_done();		pp_done();
End of changes. 30 change blocks.
38 lines changed or deleted		38 lines changed or added

---