image2d.pd (PDL-2.077) | : | image2d.pd (PDL-2.078) | ||
---|---|---|---|---|
skipping to change at line 217 | skipping to change at line 217 | |||
$str .= ' $CROAK("Out of Memory"); | $str .= ' $CROAK("Out of Memory"); | |||
poff = p_size/2; mapi += p_size-1; | poff = p_size/2; mapi += p_size-1; | |||
qoff = q_size/2; mapj += q_size-1; | qoff = q_size/2; mapj += q_size-1; | |||
'; | '; | |||
return $str; | return $str; | |||
} # sub: init_vars() | } # sub: init_vars() | |||
pp_def('conv2d', Doc=><<'EOD', | pp_def('conv2d', Doc=><<'EOD', | |||
=for ref | =for ref | |||
2D convolution of an array with a kernel (smoothing) | 2D convolution of an array with a kernel (smoothing) | |||
For large kernels, using a FFT routine, | For large kernels, using a FFT routine, | |||
such as L<fftconvolve()|PDL::FFT/fftconvolve()> in C<PDL::FFT>, | such as L<fftconvolve()|PDL::FFT/fftconvolve()> in C<PDL::FFT>, | |||
will be quicker. | will be quicker. | |||
=for usage | =for usage | |||
skipping to change at line 245 | skipping to change at line 244 | |||
Boundary - controls what values are assumed for the image when kernel | Boundary - controls what values are assumed for the image when kernel | |||
crosses its edge: | crosses its edge: | |||
=> Default - periodic boundary conditions | => Default - periodic boundary conditions | |||
(i.e. wrap around axis) | (i.e. wrap around axis) | |||
=> 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 => | |||
'Unlike the FFT routines, conv2d is able to process bad values.', | 'Unlike the FFT routines, conv2d is able to process bad values.', | |||
HandleBad => 1, | HandleBad => 1, | |||
Pars => 'a(m,n); kern(p,q); [o]b(m,n);', | Pars => 'a(m,n); kern(p,q); [o]b(m,n);', | |||
OtherPars => 'int opt;', | OtherPars => 'int opt;', | |||
PMCode => ' | PMCode => ' | |||
sub PDL::conv2d { | sub PDL::conv2d { | |||
my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; | my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; | |||
die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' | die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' | |||
if $#_<1 || $#_>2; | if $#_<1 || $#_>2; | |||
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; | |||
} | } | |||
', | ', | |||
Code => | Code => | |||
init_vars( { vars => 'PDL_Double tmp;' } ) . | ||||
init_map("i") . | ||||
init_map("j") . | ||||
' | ||||
broadcastloop %{ | ||||
for(j=0; j<n_size; j++) { | ||||
for(i=0; i<m_size; i++) { | ||||
tmp = 0; | ||||
for(j1=0; j1<q_size; j1++) { | ||||
j2 = mapj[j-j1]; | ||||
if (j2 >= 0) { | ||||
for(i1=0; i1<p_size; i1++) { | ||||
i2 = mapi[i-i1]; | ||||
if (i2 >= 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_vars( { vars => 'PDL_Double tmp; int flag;' } ) . | |||
init_map("i") . | init_map("i") . | |||
init_map("j") . | init_map("j") . | |||
' | ' | |||
broadcastloop %{ | broadcastloop %{ | |||
for(j=0; j<n_size; j++) { | loop(n) %{ | |||
for(i=0; i<m_size; i++) { | loop(m) %{ | |||
tmp = 0; | tmp = 0; | |||
for(j1=0; j1<q_size; j1++) { | for(j1=0; j1<q_size; j1++) { | |||
j2 = mapj[j-j1]; | j2 = mapj[n-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[m-i1]; | |||
if (i2 >= 0) { | if (i2 >= 0) { | |||
if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(kern(p=>i1,q =>j1)) ) { | PDL_IF_BAD(if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(k ern(p=>i1,q=>j1)) ),) { | |||
tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1); | tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1); | |||
flag = 1; | flag = 1; | |||
} /* if: good */ | } /* if: good */ | |||
} /* if: i2 >= 0 */ | } /* if: i2 >= 0 */ | |||
} /* for: i1 */ | } /* for: i1 */ | |||
} /* if: j2 >= 0 */ | } /* if: j2 >= 0 */ | |||
} /* for: j1 */ | } /* for: j1 */ | |||
if ( flag ) { $b(m=>i,n=>j) = tmp; } | PDL_IF_BAD(if ( !flag ) { $SETBAD(b()); } | |||
else { $SETBAD(b(m=>i,n=>j)); } | else,) { $b() = tmp; } | |||
} /* for: i */ | %} | |||
} /* for: j */ | %} | |||
%} | %} | |||
free(mapj+1-q_size); free(mapi+1-p_size);', | free(mapj+1-q_size); free(mapi+1-p_size);', | |||
); # pp_def: conv2d | ); # pp_def: conv2d | |||
pp_addhdr(<<'EOF'); | pp_addhdr(<<'EOF'); | |||
PDL_TYPELIST_REAL(PDL_QSORT) | PDL_TYPELIST_REAL(PDL_QSORT) | |||
EOF | EOF | |||
pp_def('med2d', Doc=> <<'EOD', | pp_def('med2d', Doc=> <<'EOD', | |||
=for ref | =for ref | |||
2D median-convolution of an array with a kernel (smoothing) | 2D median-convolution of an array with a kernel (smoothing) | |||
Note: only points in the kernel E<gt>0 are included in the median, other | Note: only points in the kernel E<gt>0 are included in the median, other | |||
points are weighted by the kernel value (medianing lots of zeroes | points are weighted by the kernel value (medianing lots of zeroes | |||
is rather pointless) | is rather pointless) | |||
=for usage | =for usage | |||
skipping to change at line 359 | skipping to change at line 327 | |||
=for options | =for options | |||
Boundary - controls what values are assumed for the image when kernel | Boundary - controls what values are assumed for the image when kernel | |||
crosses its edge: | crosses its edge: | |||
=> Default - periodic boundary conditions (i.e. wrap around axis) | => Default - periodic boundary conditions (i.e. wrap around axis) | |||
=> 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);', | Pars => 'a(m,n); kern(p,q); [o]b(m,n); double [t]tmp(pq);', | |||
OtherPars => 'int opt;', | OtherPars => 'int opt;', | |||
PMCode => ' | 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 | ||||
', | RedoDimsCode => '$SIZE(pq) = $SIZE(p)*$SIZE(q);', | |||
Code => | Code => | |||
init_vars( { vars => 'PDL_Double *tmp, kk; int count;', | init_vars( { vars => 'PDL_Double kk, aa; int flag;' } ) . | |||
malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double) | ||||
);', | ||||
check => '(tmp==NULL) || ' } ) . | ||||
init_map("i") . | ||||
init_map("j") . | ||||
' | ||||
broadcastloop %{ | ||||
for(j=0; j<n_size; j++) { | ||||
for(i=0; i<m_size; i++) { | ||||
count = 0; | ||||
for(j1=0; j1<q_size; j1++) { | ||||
j2 = mapj[j-j1]; | ||||
if (j2 >= 0) | ||||
for(i1=0; i1<p_size; i1++) { | ||||
i2 = mapi[i-i1]; | ||||
if (i2 >= 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 */ | ||||
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 kk, aa; int count, flag;' } ) . | ||||
init_map("i") . | init_map("i") . | |||
init_map("j") . | init_map("j") . | |||
' | ' | |||
PDL_Double tmp[p_size*q_size]; | ||||
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++) { | |||
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]; | |||
if (i2 >= 0) { | if (i2 >= 0) { | |||
kk = $kern(p=>i1,q=>j1); | kk = $kern(p=>i1,q=>j1); | |||
aa = $a(m=>i2,n=>j2); | aa = $a(m=>i2,n=>j2); | |||
if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(aa,a) ) { | PDL_IF_BAD(if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(a a,a) ),) { | |||
flag = 1; | flag = 1; | |||
if ( kk > 0 ) { | if ( kk > 0 ) { | |||
tmp[count++] = aa * kk; | $tmp(pq=>count++) = aa * kk; | |||
} | } | |||
} | } | |||
} /* if: i2 >= 0 */ | } /* if: i2 >= 0 */ | |||
} /* for: i1 */ | } /* for: i1 */ | |||
} /* for: j1 */ | } /* for: j1 */ | |||
if ( flag == 0 ) { | PDL_IF_BAD(if ( flag == 0 ) { | |||
$SETBAD(b(m=>i,n=>j)); | $SETBAD(b(m=>i,n=>j)); | |||
} else { | } else,) { | |||
qsort_D( $P(tmp), 0, count-1 ); | ||||
qsort_D( tmp, 0, count-1 ); | $b(m=>i,n=>j) = $tmp(pq=>(count-1)/2); | |||
$b(m=>i,n=>j) = tmp[(count-1)/2]; | ||||
} | } | |||
} /* for: i */ | } /* for: i */ | |||
} /* for: j */ | } /* for: j */ | |||
%} | %} | |||
free(mapj+1-q_size); free(mapi+1-p_size); | free(mapj+1-q_size); free(mapi+1-p_size); | |||
' | ' | |||
); # pp_def: med2d | ); # pp_def: med2d | |||
pp_def('med2df', Doc=> <<'EOD', | pp_def('med2df', Doc=> <<'EOD', | |||
skipping to change at line 498 | skipping to change at line 426 | |||
=> 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 | |||
Pars => 'a(m,n); [o]b(m,n);', | Pars => 'a(m,n); [o]b(m,n);', | |||
# funny parameter names to avoid special case in 'init_vars' | # funny parameter names to avoid special case in 'init_vars' | |||
OtherPars => 'int __p_size; int __q_size; int opt;', | OtherPars => 'int __p_size; int __q_size; int opt;', | |||
PMCode => ' | PMCode => <<'EOF', | |||
sub PDL::med2df { | sub PDL::med2df { | |||
my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; | my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH'; | |||
die \'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )\' | die 'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )' | |||
if $#_<2 || $#_>3; | if $#_<2 || $#_>3; | |||
my($x,$p,$q) = @_; | my($x,$p,$q) = @_; | |||
croak "med2df: kernel must contain some positive elements.\n" | croak "med2df: kernel must contain some positive elements.\n" | |||
if $p == 0 && $q == 0; | if $p == 0 && $q == 0; | |||
my $c = $#_ == 3 ? $_[3] : $x->nullcreate; | my $c = $#_ == 3 ? $_[3] : $x->nullcreate; | |||
&PDL::_med2df_int($x,$c,$p,$q, | &PDL::_med2df_int($x,$c,$p,$q, | |||
(!(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 | ||||
', | ||||
Code => | Code => | |||
init_vars( { vars => '$GENERIC() kk; int count;' }, [qw(m n)], [qw(p q)] ) . | init_vars( { vars => '$GENERIC() kk; int count;' }, [qw(m n)], [qw(p q)] ) . | |||
init_map("i") . | init_map("i") . | |||
init_map("j") . | init_map("j") . | |||
' | ' | |||
$GENERIC() tmp[p_size*q_size]; | $GENERIC() tmp[p_size*q_size]; | |||
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++) { | |||
count = 0; | count = 0; | |||
skipping to change at line 731 | skipping to change at line 657 | |||
is performed (see L</patch2d>). | is performed (see L</patch2d>). | |||
=cut | =cut | |||
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);', | |||
Code => 'loop(n,m) %{ $b() = $a(); %}', # just copy | ||||
CopyBadStatusCode => '', # handled by BadCode | CopyBadStatusCode => '', # handled by BadCode | |||
BadCode => | Code => | |||
'int m_size, n_size, i,j, i1,j1, i2,j2, norm, flag; | 'int i1,j1, flag = 0; | |||
double tmp; | PDL_Indx m_size = $SIZE(m), n_size = $SIZE(n); | |||
$GENERIC(a) a_val; | ||||
flag = 0; | ||||
m_size = $SIZE(m); n_size = $SIZE(n); | ||||
broadcastloop %{ | broadcastloop %{ | |||
loop(n) %{ | ||||
for(j=0; j<n_size; j++) { | loop(m) %{ | |||
for(i=0; i<m_size; i++) { | $GENERIC(a) a_val = $a(); | |||
PDL_IF_BAD(if ( !$ISGOODVAR(a_val,a) ) { | ||||
a_val = $a(m=>i,n=>j); | double tmp = 0; PDL_Indx norm=0; | |||
if ( $ISGOODVAR(a_val,a) ) { | ||||
$b(m=>i,n=>j) = a_val; | ||||
} else { | ||||
tmp = 0; norm=0; | ||||
for(j1=-1; j1<=1; j1++) { | for(j1=-1; j1<=1; j1++) { | |||
j2 = j+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 ) { | |||
i2 = i+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) ) { | |||
tmp += a_val; | tmp += a_val; | |||
norm++; | norm++; | |||
} | } | |||
} | } | |||
} /* if: i1!=0 || j1!=0 */ | } /* if: i1!=0 || j1!=0 */ | |||
} /* for: i1 */ | } /* for: i1 */ | |||
} | } | |||
} /* for: j1 */ | } /* for: j1 */ | |||
/* Patch */ | /* Patch */ | |||
if (norm>0) { | if (norm>0) { | |||
$b(m=>i,n=>j) = tmp/norm; | $b() = tmp/norm; | |||
} else { | } else { | |||
$SETBAD(b(m=>i,n=>j)); | $SETBAD(b()); | |||
flag = 1; | flag = 1; | |||
} | } | |||
} else,) { | ||||
$b() = a_val; | ||||
} /* if: ISGOODVAR() */ | } /* if: ISGOODVAR() */ | |||
%} | ||||
} /* for: i */ | %} | |||
} /* for: j */ | ||||
%} /* broadcastloop */ | %} /* broadcastloop */ | |||
/* handle bad flag */ | /* handle bad flag */ | |||
if ( flag ) $PDLSTATESETBAD(b); | PDL_IF_BAD(if ( flag ) $PDLSTATESETBAD(b);,) | |||
', # BadCode | ', | |||
); | ); | |||
pp_def('max2d_ind', | pp_def('max2d_ind', | |||
Doc=><<'EOD', | Doc=><<'EOD', | |||
=for ref | =for ref | |||
Return value/position of maximum value in 2D image | Return value/position of maximum value in 2D image | |||
Contributed by Tim Jeness | Contributed by Tim Jeness | |||
skipping to change at line 814 | 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; int curind1; int curind2; | double cur = 0; | |||
curind1=0; | PDL_Indx curind1 = PDL_IF_BAD(-1,0), curind2 = PDL_IF_BAD(-1,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(m) %{ | |||
loop(n) %{ | loop(n) %{ | |||
if( $ISGOOD(a()) && ( (!n && !m) || ($a() > 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; | |||
} | } | |||
%} | %} | |||
%} | %} | |||
if ( curind1 < 0 ) { | PDL_IF_BAD(if ( curind1 < 0 ) { | |||
$SETBAD(val()); | $SETBAD(val()); | |||
$SETBAD(x()); | $SETBAD(x()); | |||
$SETBAD(y()); | $SETBAD(y()); | |||
} else { | } else,) { | |||
$val() = cur; | $val() = cur; | |||
$x() = curind1; | $x() = curind1; | |||
$y() = curind2; | $y() = curind2; | |||
} | } | |||
'); | '); | |||
pp_def('centroid2d', | pp_def('centroid2d', | |||
Doc=><<'EOD', | Doc=><<'EOD', | |||
=for ref | =for ref | |||
Refine a list of object positions in 2D image by centroiding in a box | 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 | C<$box> is the full-width of the box, i.e. the window | |||
is C<+/- $box/2>. | is C<+/- $box/2>. | |||
=cut | =cut | |||
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 => ' | |||
int i,j,i1,i2,j1,j2,m_size,n_size; | PDL_Indx i,j,i1,i2,j1,j2,m_size = $SIZE(m),n_size = $SIZE(n); | |||
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; | double sum,data,sumx,sumy; | |||
m_size = $SIZE(m); n_size = $SIZE(n); | ||||
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); | |||
if ( $ISGOODVAR(data,im) ) { | PDL_IF_BAD(if ( $ISGOODVAR(data,im) ),) { | |||
sum += data; | sum += data; | |||
sumx += data*i; | sumx += data*i; | |||
sumy += data*j; | sumy += data*j; | |||
} | } | |||
} | } | |||
} | } | |||
/* | /* | |||
* if sum == 0 then we will flag as bad -- although it could just mean that | * if sum == 0 then we will flag as bad -- although it could just mean that | |||
* there is negative values in the dataset. | * there is negative values in the dataset. | |||
* - should use a better check than != 0.0 ... | * - should use a better check than != 0.0 ... | |||
*/ | */ | |||
if ( sum != 0.0 ) { | PDL_IF_BAD(if ( sum == 0.0 ) { | |||
$xcen() = sumx/sum; | ||||
$ycen() = sumy/sum; | ||||
} else { | ||||
$SETBAD(xcen()); | $SETBAD(xcen()); | |||
$SETBAD(ycen()); | $SETBAD(ycen()); | |||
} else,) { | ||||
$xcen() = sumx/sum; | ||||
$ycen() = sumy/sum; | ||||
} | } | |||
' | ' | |||
); | ); | |||
pp_addhdr(' | pp_addhdr(' | |||
void AddEquiv ( PDL_Long* equiv, PDL_Long i, PDL_Long j); | void AddEquiv ( PDL_Long* equiv, PDL_Long i, PDL_Long j); | |||
'); | '); | |||
pp_add_exported('', 'cc8compt','cc4compt'); | pp_add_exported('', 'cc8compt','cc4compt'); | |||
pp_addpm(<<'EOPM'); | pp_addpm(<<'EOPM'); | |||
skipping to change at line 2163 | skipping to change at line 2026 | |||
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 ; | int i, j, k ; | |||
int ncoeff, lx_out, ly_out ; | PDL_Indx ncoeff, lx_out, ly_out ; | |||
int lx_3, ly_3 ; | PDL_Indx lx_3, ly_3 ; | |||
double cur ; | double cur ; | |||
double neighbors[16] ; | double neighbors[16] ; | |||
double rsc[8], | double rsc[8], | |||
sumrs ; | sumrs ; | |||
double x, y ; | double x, y ; | |||
int px, py ; | int px, py ; | |||
int tabx, taby ; | int tabx, taby ; | |||
double *kernel, *poly ; | double *kernel, *poly ; | |||
int da[16], db[16] ; | int da[16], db[16] ; | |||
End of changes. 60 change blocks. | ||||
198 lines changed or deleted | 60 lines changed or added |