(Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).
| image2d.pd (PDL-2.081) | : | image2d.pd (PDL-2.082) |
|---|
| | |
| skipping to change at line 28 | | skipping to change at line 28 |
|---|
| =head1 SYNOPSIS | | =head1 SYNOPSIS |
| | |
| use PDL::Image2D; | | use PDL::Image2D; |
| | |
| =cut | | =cut |
| | |
| use PDL; # ensure qsort routine available | | use PDL; # ensure qsort routine available |
| use PDL::Math; | | use PDL::Math; |
| use Carp; | | use Carp; |
| | |
|
| | my %boundary2value = (Reflect=>1, Truncate=>2, Replicate=>3); |
| EOD | | EOD |
| | |
| pp_addpm({At=>'Bot'},<<'EOD'); | | pp_addpm({At=>'Bot'},<<'EOD'); |
| | |
| =head1 AUTHORS | | =head1 AUTHORS |
| | |
| Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams | | Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams |
|
| (rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu), | | (rjrw@ast.leeds.ac.uk), Tim Jenness (timj@jach.hawaii.edu), |
| and Doug Burke (burke@ifa.hawaii.edu). | | and Doug Burke (burke@ifa.hawaii.edu). |
| | |
| All rights reserved. There is no warranty. You are allowed | | All rights reserved. There is no warranty. You are allowed |
| to redistribute this software / documentation under certain | | to redistribute this software / documentation under certain |
| conditions. For details, see the file COPYING in the PDL | | conditions. For details, see the file COPYING in the PDL |
| distribution. If this file is separated from the PDL distribution, | | distribution. If this file is separated from the PDL distribution, |
| the copyright notice should be included in the file. | | the copyright notice should be included in the file. |
| | |
| =cut | | =cut |
| | |
| | |
| skipping to change at line 188 | | skipping to change at line 189 |
|---|
| break; | | break; |
| default: | | default: |
| REALMOD($loop2,$size); | | REALMOD($loop2,$size); |
| } | | } |
| map${var}\[$loop] = $loop2; | | map${var}\[$loop] = $loop2; |
| }\n"; | | }\n"; |
| | |
| } # sub: init_map() | | } # sub: init_map() |
| | |
| sub init_vars { | | sub init_vars { |
|
| my $href = shift || { }; | | |
| my $sizevars = shift || [qw(m n p q)]; | | my $sizevars = shift || [qw(m n p q)]; |
| my $compvars = shift || []; | | my $compvars = shift || []; |
|
| $href->{vars} = '' unless defined $href->{vars}; | | my $str = "int i,j, i1,j1, i2,j2, poff, qoff;"; |
| $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 .= | | $str .= |
| 'int opt = $COMP(opt); | | 'int opt = $COMP(opt); |
| ' . join("\n", map "int ${_}_size = \$SIZE(${_});", @$sizevars) . ' | | ' . join("\n", map "int ${_}_size = \$SIZE(${_});", @$sizevars) . ' |
|
| ' . join("\n", map "int ${_}_size = \$COMP(__${_}_size);", @$compvars) . | | ' . join("\n", map "int ${_}_size = \$COMP(${_}_size);", @$compvars) . ' |
| ' | | PDL_Indx *mapi = $P(mapi), *mapj = $P(mapj); |
| int *mapi, *mapj; | | if ((mapi==NULL) || (mapj==NULL)) $CROAK("Out of Memory"); |
| | |
| 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 .= ' $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<PDL::FFT/fftconvolve>, |
| will be quicker. | | will be quicker. |
| | |
| =for usage | | =for usage |
| | |
| $new = conv2d $old, $kernel, {OPTIONS} | | $new = conv2d $old, $kernel, {OPTIONS} |
| | |
| =for example | | =for example |
| | |
| $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect} | | $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect} |
| | |
| | |
| skipping to change at line 249 | | skipping to change at line 238 |
|---|
| (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); indx [t]mapi(isize); indx [t]mapj |
| | (jsize)', |
| | RedoDimsCode => ' |
| | $SIZE(isize) = $SIZE(p) + $SIZE(m); |
| | $SIZE(jsize) = $SIZE(q) + $SIZE(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: | | (!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}} |
| (($$opt{Boundary} eq "Reflect") + | | ); |
| 2*($$opt{Boundary} eq "Truncate") + | | |
| 3*($$opt{Boundary} eq "Replicate"))); | | |
| return $c; | | return $c; |
| } | | } |
| ', | | ', |
| GenericTypes => $A, | | GenericTypes => $A, |
| Code => | | Code => |
|
| init_vars( { vars => 'PDL_CLDouble tmp; int flag;' } ) . | | init_vars() . |
| init_map("i") . | | init_map("i") . |
| init_map("j") . | | init_map("j") . |
| ' | | ' |
|
| | PDL_CLDouble tmp; int flag; |
| 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) { |
| for(i1=0; i1<p_size; i1++) { | | for(i1=0; i1<p_size; i1++) { |
| i2 = mapi[m-i1]; | | i2 = mapi[m-i1]; |
| if (i2 >= 0) { | | if (i2 >= 0) { |
| | |
| skipping to change at line 295 | | skipping to change at line 287 |
|---|
| } /* if: good */ | | } /* if: good */ |
| } /* if: i2 >= 0 */ | | } /* if: i2 >= 0 */ |
| } /* for: i1 */ | | } /* for: i1 */ |
| } /* if: j2 >= 0 */ | | } /* if: j2 >= 0 */ |
| } /* for: j1 */ | | } /* for: j1 */ |
| PDL_IF_BAD(if ( !flag ) { $SETBAD(b()); } | | PDL_IF_BAD(if ( !flag ) { $SETBAD(b()); } |
| else,) { $b() = tmp; } | | else,) { $b() = tmp; } |
| %} | | %} |
| %} | | %} |
| %} | | %} |
|
| 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) |
| | |
| skipping to change at line 334 | | skipping to change at line 326 |
|---|
| => 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); indx [t]mapi(
isize); indx [t]mapj(jsize)', |
| 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: | | (!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}} |
| (($$opt{Boundary} eq "Reflect") + | | ); |
| 2*($$opt{Boundary} eq "Truncate") + | | |
| 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); |
| | $SIZE(isize) = $SIZE(p) + $SIZE(m); |
| | $SIZE(jsize) = $SIZE(q) + $SIZE(n); |
| | ', |
| Code => | | Code => |
|
| init_vars( { vars => 'int flag;' } ) . | | init_vars() . |
| init_map("i") . | | init_map("i") . |
| init_map("j") . | | init_map("j") . |
| ' | | ' |
|
| | int flag; |
| PDL_LDouble kk, aa; | | 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++) { |
| | |
| skipping to change at line 385 | | skipping to change at line 380 |
|---|
| if ( kk > 0 ) { | | if ( kk > 0 ) { |
| $tmp(pq=>count++) = aa * kk; | | $tmp(pq=>count++) = aa * kk; |
| } | | } |
| } | | } |
| } /* if: i2 >= 0 */ | | } /* if: i2 >= 0 */ |
| } /* for: i1 */ | | } /* for: i1 */ |
| } /* for: j1 */ | | } /* for: j1 */ |
| PDL_IF_BAD(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_$PPSYM(tmp)( $P(tmp), 0, count-1 ); |
| $b(m=>i,n=>j) = $tmp(pq=>(count-1)/2); | | $b(m=>i,n=>j) = $tmp(pq=>(count-1)/2); |
| } | | } |
| } /* for: i */ | | } /* for: i */ |
| } /* for: j */ | | } /* for: j */ |
| %} | | %} |
|
| 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', |
| | |
| =for ref | | =for ref |
| | |
| 2D median-convolution of an array in a pxq window (smoothing) | | 2D median-convolution of an array in a pxq window (smoothing) |
| | |
| | |
| skipping to change at line 426 | | skipping to change at line 420 |
|---|
| 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 |
|
| Pars => 'a(m,n); [o]b(m,n);', | | Pars => 'a(m,n); [o]b(m,n); indx [t]mapi(isize); indx [t]mapj(jsize)', |
| # funny parameter names to avoid special case in 'init_vars' | | OtherPars => 'int p_size=>p; int q_size=>q; int opt;', |
| OtherPars => 'int __p_size; int __q_size; int opt;', | | RedoDimsCode => ' |
| | $SIZE(isize) = $SIZE(p) + $SIZE(m); |
| | $SIZE(jsize) = $SIZE(q) + $SIZE(n); |
| | ', |
| PMCode => <<'EOF', | | 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: | | (!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}} |
| (($$opt{Boundary} eq "Reflect") + | | ); |
| 2*($$opt{Boundary} eq "Truncate") + | | |
| 3*($$opt{Boundary} eq "Replicate"))); | | |
| return $c; | | return $c; |
| } | | } |
| EOF | | EOF |
| Code => | | Code => |
|
| init_vars( { vars => '$GENERIC() kk; int count;' }, [qw(m n)], [qw(p
q)] ) . | | init_vars() . |
| init_map("i") . | | init_map("i") . |
| init_map("j") . | | init_map("j") . |
| ' | | ' |
|
| $GENERIC() tmp[p_size*q_size]; | | int count; |
| | $GENERIC() kk, 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; |
| 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 474 | | skipping to change at line 470 |
|---|
| } /* if: i2 >= 0 */ | | } /* if: i2 >= 0 */ |
| } /* for: i1 */ | | } /* for: i1 */ |
| } /* for: j1 */ | | } /* for: j1 */ |
| | |
| $b(m=>i,n=>j) = | | $b(m=>i,n=>j) = |
| quick_select_$PPSYM() (tmp, count ); | | quick_select_$PPSYM() (tmp, count ); |
| | |
| } /* for: i */ | | } /* for: i */ |
| } /* for: j */ | | } /* for: j */ |
| %} | | %} |
|
| free(mapj+1-q_size); free(mapi+1-p_size); | | |
| ', | | ', |
| | |
| ); # 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', |
| | |
| skipping to change at line 710 | | skipping to change at line 705 |
|---|
| ', | | ', |
| ); | | ); |
| | |
| 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 Jenness |
| | |
| =cut | | =cut |
| | |
| EOD | | EOD |
| | |
| 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. |
| | |
| | |
| skipping to change at line 1486 | | skipping to change at line 1481 |
|---|
| | |
| C<$px> and C<$py> are C<np> by C<np> element ndarrays which describe | | C<$px> and C<$py> are C<np> by C<np> element ndarrays which describe |
| a polynomial mapping (of order C<np-1>) | | a polynomial mapping (of order C<np-1>) |
| from the I<output> C<(u,v)> image to the I<input> C<(x,y)> image: | | from the I<output> C<(u,v)> image to the I<input> C<(x,y)> image: |
| | |
| x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j | | 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 | | 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 | | The transformation is returned for the reverse direction (ie |
| output to input image) since that is what is required by the | | output to input image) since that is what is required by the |
|
| L<warp2d()|/warp2d> routine. The L<applywarp2d()|/applywarp2d> | | L</warp2d> routine. The L</applywarp2d> |
| routine can be used to convert a set of C<$u,$v> points given | | routine can be used to convert a set of C<$u,$v> points given |
| C<$px> and C<$py>. | | C<$px> and C<$py>. |
| | |
| Options: | | Options: |
| | |
| =for options | | =for options |
| | |
| FIT - which terms to fit? default ones(byte,$nf,$nf) | | FIT - which terms to fit? default ones(byte,$nf,$nf) |
| | |
| =begin comment | | =begin comment |
| | |
| skipping to change at line 1613 | | skipping to change at line 1608 |
|---|
| =for ref | | =for ref |
| | |
| Transform a set of points using a 2-D polynomial mapping | | Transform a set of points using a 2-D polynomial mapping |
| | |
| =for usage | | =for usage |
| | |
| ( $x, $y ) = applywarp2d( $px, $py, $u, $v ) | | ( $x, $y ) = applywarp2d( $px, $py, $u, $v ) |
| | |
| Convert a set of points (stored in 1D ndarrays C<$u,$v>) | | Convert a set of points (stored in 1D ndarrays C<$u,$v>) |
| to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px> | | to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px> |
|
| and C<$py>. See L<fitwarp2d()|/fitwarp2d> | | and C<$py>. See L</fitwarp2d> |
| for more information on the format of C<$px> and C<$py>. | | for more information on the format of C<$px> and C<$py>. |
| | |
| =cut | | =cut |
| | |
| # use SVD to fit data. Assuming no errors. | | # use SVD to fit data. Assuming no errors. |
| | |
| =pod | | =pod |
| | |
| =begin comment | | =begin comment |
| | |
| | |
| skipping to change at line 1889 | | skipping to change at line 1884 |
|---|
| *applywarp2d = \&PDL::applywarp2d; | | *applywarp2d = \&PDL::applywarp2d; |
| | |
| ' ); | | ' ); |
| | |
| ## resampling routines taken from v3.6-0 of the Eclipse package | | ## resampling routines taken from v3.6-0 of the Eclipse package |
| ## http://www.eso.org/eclipse by Nicolas Devillard | | ## http://www.eso.org/eclipse by Nicolas Devillard |
| ## | | ## |
| | |
| pp_addhdr( '#include "resample.h"' . "\n" ); | | pp_addhdr( '#include "resample.h"' . "\n" ); |
| | |
|
| # pod for warp2d | | pp_def( |
| # and support routine | | 'warp2d', |
| # | | Doc=> <<'EOD', |
| pp_addpm( <<'EOD'); | | |
| | |
| =head2 warp2d | | |
| | |
| =for sig | | |
| | |
| Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { opt | | |
| ions }) | | |
| | | |
| =for ref | | =for ref |
| | |
| Warp a 2D image given a polynomial describing the I<reverse> mapping. | | Warp a 2D image given a polynomial describing the I<reverse> mapping. |
| | |
| =for usage | | =for usage |
| | |
| $out = warp2d( $img, $px, $py, { options } ); | | $out = warp2d( $img, $px, $py, { options } ); |
| | |
| Apply the polynomial transformation encoded in the C<$px> and | | Apply the polynomial transformation encoded in the C<$px> and |
| C<$py> ndarrays to warp the input image C<$img> into the output | | C<$py> ndarrays to warp the input image C<$img> into the output |
| image C<$out>. | | image C<$out>. |
| | |
| The format for the polynomial transformation is described in | | The format for the polynomial transformation is described in |
|
| the documentation for the L<fitwarp2d()|/fitwarp2d> routine. | | the documentation for the L</fitwarp2d> routine. |
| | |
| At each point C<x,y>, the closest 16 pixel values are combined | | At each point C<x,y>, the closest 16 pixel values are combined |
| with an interpolation kernel to calculate the value at C<u,v>. | | with an interpolation kernel to calculate the value at C<u,v>. |
| The interpolation is therefore done in the image, rather than | | The interpolation is therefore done in the image, rather than |
| Fourier, domain. | | Fourier, domain. |
| By default, a C<tanh> kernel is used, but this can be changed | | By default, a C<tanh> kernel is used, but this can be changed |
| using the C<KERNEL> option discussed below | | using the C<KERNEL> option discussed below |
| (the choice of kernel depends on the frequency content of the input image). | | (the choice of kernel depends on the frequency content of the input image). |
| | |
| The routine is based on the C<warping> command from | | The routine is based on the C<warping> command from |
| | |
| skipping to change at line 1958 | | skipping to change at line 1945 |
|---|
| | |
| The options are: | | The options are: |
| | |
| =for options | | =for options |
| | |
| KERNEL - default value is tanh | | KERNEL - default value is tanh |
| NOVAL - default value is 0 | | NOVAL - default value is 0 |
| | |
| C<KERNEL> is used to specify which interpolation kernel to use | | C<KERNEL> is used to specify which interpolation kernel to use |
| (to see what these kernels look like, use the | | (to see what these kernels look like, use the |
|
| L<warp2d_kernel()|/warp2d_kernel> routine). | | L</warp2d_kernel> routine). |
| The options are: | | The options are: |
| | |
| =over 4 | | =over 4 |
| | |
| =item tanh | | =item tanh |
| | |
| Hyperbolic tangent: the approximation of an ideal box filter by the | | Hyperbolic tangent: the approximation of an ideal box filter by the |
| product of symmetric tanh functions. | | product of symmetric tanh functions. |
| | |
| =item sinc | | =item sinc |
| | |
| skipping to change at line 2007 | | skipping to change at line 1994 |
|---|
| | |
| =item hamming | | =item hamming |
| | |
| This kernel uses the same C<H(x)> as the Hann filter, but with | | This kernel uses the same C<H(x)> as the Hann filter, but with |
| C<a = 0.54>. | | C<a = 0.54>. |
| | |
| =back | | =back |
| | |
| C<NOVAL> gives the value used to indicate that a pixel in the | | C<NOVAL> gives the value used to indicate that a pixel in the |
| output image does not map onto one in the input image. | | output image does not map onto one in the input image. |
|
| | | EOD |
| =cut | | HandleBad => 0, |
| | | Pars => 'img(m,n); ldouble px(np,np); ldouble py(np,np); [o] warp(m,n); ldou |
| | ble [t] poly(np); ldouble [t] kernel(ns)', |
| | OtherPars => 'char *kernel_type; double noval; int nsamples => ns', |
| | GenericTypes => [ qw(F D E) ], |
| | PMCode => <<'EOF', |
| # support routine | | # support routine |
| { | | { |
| my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann ); | | my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann ); |
|
| | | |
| # note: convert to lower case | | # note: convert to lower case |
| sub _check_kernel ($$) { | | sub _check_kernel ($$) { |
| my $kernel = lc shift; | | my $kernel = lc shift; |
| my $code = shift; | | my $code = shift; |
| barf "Unknown kernel $kernel sent to $code\n" . | | barf "Unknown kernel $kernel sent to $code\n" . |
| "\tmust be one of [" . join(',',keys %warp2d) . "]\n" | | "\tmust be one of [" . join(',',keys %warp2d) . "]\n" |
| unless exists $warp2d{$kernel}; | | unless exists $warp2d{$kernel}; |
| return $kernel; | | return $kernel; |
| } | | } |
| } | | } |
| | |
|
| EOD | | |
| | |
| pp_def( | | |
| 'warp2d', | | |
| Doc=> undef, | | |
| HandleBad => 0, | | |
| Pars => 'img(m,n); ldouble px(np,np); ldouble py(np,np); [o] warp(m,n);', | | |
| OtherPars => 'char *kernel_type; double noval;', | | |
| GenericTypes => [ qw(F D E) ], | | |
| 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; |
| my $py = shift; | | my $py = shift; |
| my $out = $#_ == -1 ? PDL->null() : shift; | | my $out = $#_ == -1 ? PDL->null() : shift; |
|
| | | |
| # safety checks | | # safety checks |
| 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}, _get_kernel |
| &PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} ); | | _size() ); |
| return $out; | | return $out; |
| } | | } |
|
| | | EOF |
| ', | | Code => <<'EOF', |
| Code => ' | | |
| | | |
| 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, i, j, k, lx_3, ly_3, px, py, tabx, taby; |
| PDL_Indx da[16], db[16] ; | | PDL_Indx da[16], db[16] ; |
|
| PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel, *poly ; | | PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel = $P(kernel), *
poly = $P(poly); |
| | |
|
| /* Generate default interpolation kernel */ | | /* Generate interpolation kernel */ |
| kernel = generate_interpolation_kernel( $COMP(kernel_type) ) ; | | if (!generate_interpolation_kernel($COMP(kernel_type), $SIZE(ns), kernel)) |
| if (kernel == NULL) { | | $CROAK("Invalid kernel type '%s'", $COMP(kernel_type)); |
| $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? */ |
| ly_out = $SIZE(n); | | ly_out = $SIZE(n); |
| lx_3 = lx_out - 3; | | lx_3 = lx_out - 3; |
| ly_3 = ly_out - 3; | | ly_3 = ly_out - 3; |
| | |
| /* Pre compute leaps for 16 closest neighbors positions */ | | /* Pre compute leaps for 16 closest neighbors positions */ |
| da[0] = -1; db[0] = -1; | | da[0] = -1; db[0] = -1; |
| | |
| skipping to change at line 2096 | | skipping to change at line 2067 |
|---|
| da[8] = -1; db[8] = 1; | | da[8] = -1; db[8] = 1; |
| 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 */ | | |
| poly = malloc( ncoeff * sizeof(PDL_LDouble) ); | | |
| if ( poly == NULL ) { | | |
| $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] = (PDL_LDouble) n * poly[k-1]; | | poly[k] = (PDL_LDouble) n * poly[k-1]; |
| } | | } |
| | |
| skipping to change at line 2177 | | skipping to change at line 2143 |
|---|
| rsc[3]*neighbors[15] ) ; | | rsc[3]*neighbors[15] ) ; |
| | |
| /* Copy the value to the output image */ | | /* Copy the value to the output image */ |
| $warp() = ($GENERIC()) (cur/sumrs); | | $warp() = ($GENERIC()) (cur/sumrs); |
| | |
| } /* if: edge or interior */ | | } /* if: edge or interior */ |
| | |
| %} /* loop(m) */ | | %} /* loop(m) */ |
| %} /* loop(n) */ | | %} /* loop(n) */ |
| %} /* broadcastloop */ | | %} /* broadcastloop */ |
|
| | | EOF |
| free(poly); | | |
| kernel_free(kernel) ; | | |
| | |
| ', | | |
| ); # pp_def: warp2d | | ); # pp_def: warp2d |
| | |
| pp_addxs( ' | | pp_addxs( ' |
| | |
| int | | int |
| _get_kernel_size() | | _get_kernel_size() |
| PROTOTYPE: | | PROTOTYPE: |
| CODE: | | CODE: |
| RETVAL = KERNEL_SAMPLES; | | RETVAL = KERNEL_SAMPLES; |
| OUTPUT: | | OUTPUT: |
| RETVAL | | RETVAL |
| | |
| '); | | '); |
| | |
|
| pp_add_exported('', 'warp2d_kernel'); | | pp_def( 'warp2d_kernel', |
| pp_addpm( ' | | Doc => <<'EOF', |
| | |
| =head2 warp2d_kernel | | |
| | | |
| =for ref | | =for ref |
| | |
| Return the specified kernel, as used by L</warp2d> | | Return the specified kernel, as used by L</warp2d> |
| | |
| =for usage | | =for usage |
| | |
| ( $x, $k ) = warp2d_kernel( $name ) | | ( $x, $k ) = warp2d_kernel( $name ) |
| | |
| The valid values for C<$name> are the same as the C<KERNEL> option | | The valid values for C<$name> are the same as the C<KERNEL> option |
|
| of L<warp2d()|/warp2d>. | | of L</warp2d>. |
| | |
| =for example | | =for example |
| | |
| line warp2d_kernel( "hamming" ); | | line warp2d_kernel( "hamming" ); |
|
| | | EOF |
| =cut | | |
| | |
| '); # pp_addpm | | |
| | |
| # this is not very clever, but it's a pain to create a valid | | |
| # ndarray in XS code | | |
| # | | |
| pp_def( 'warp2d_kernel', | | |
| Doc => undef, | | |
| HandleBad => 0, | | HandleBad => 0, |
| PMCode => ' | | PMCode => ' |
|
| | | |
| sub PDL::warp2d_kernel ($) { | | sub PDL::warp2d_kernel ($) { |
| my $kernel = _check_kernel( shift, "warp2d_kernel" ); | | my $kernel = _check_kernel( shift, "warp2d_kernel" ); |
|
| | | &PDL::_warp2d_kernel_int( my $x=PDL->null, my $k=PDL->null, $kernel, _get_ke |
| my $nelem = _get_kernel_size(); | | rnel_size() ); |
| my $x = zeroes( $nelem ); | | |
| my $k = zeroes( $nelem ); | | |
| | |
| &PDL::_warp2d_kernel_int( $x, $k, $kernel ); | | |
| return ( $x, $k ); | | return ( $x, $k ); |
|
| | | |
| # return _get_kernel( $kernel ); | | |
| } | | } |
|
| *warp2d_kernel = \&PDL::warp2d_kernel; | | |
| | | |
| ', | | ', |
|
| Pars => '[o] x(n); [o] k(n);', | | Pars => '[o] x(n); [o] k(n); ldouble [t] kernel(n)', |
| OtherPars => 'char *name;', | | OtherPars => 'char *name; PDL_Indx nsize => n', |
| GenericTypes => [ qw(F D E) ], | | GenericTypes => [ qw(F D E) ], |
|
| Code => ' | | Code => <<'EOF', |
| | | PDL_LDouble *kernel = $P(kernel); |
| PDL_LDouble *kernel, xx; | | if (!generate_interpolation_kernel($COMP(name), $SIZE(n), kernel)) |
| | $CROAK("Invalid kernel type '%s'", $COMP(name)); |
| 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 ndarrays */ | | /* fill in ndarrays */ |
|
| xx = 0.0; | | PDL_LDouble xx = 0.0; |
| broadcastloop %{ | | broadcastloop %{ |
| loop (n) %{ | | loop (n) %{ |
| $x() = xx; | | $x() = xx; |
| $k() = kernel[n]; | | $k() = kernel[n]; |
| xx += 1.0 / (double) TABSPERPIX; | | xx += 1.0 / (double) TABSPERPIX; |
| %} | | %} |
| %} | | %} |
|
| | | EOF |
| /* free the kernel */ | | ); |
| kernel_free( kernel ); | | |
| | | |
| '); | | |
| | |
| pp_done(); | | pp_done(); |
| | | |
| End of changes. 51 change blocks. |
|---|
| 149 lines changed or deleted | | 79 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 