"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/GENERATED/PDL/ImageND.pm" (28 May 2022, 16397 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ImageND.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

    1 #
    2 # GENERATED WITH PDL::PP! Don't modify!
    3 #
    4 package PDL::ImageND;
    5 
    6 our @EXPORT_OK = qw(kernctr convolve ninterpol rebin circ_mean circ_mean_p convolveND );
    7 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
    8 
    9 use PDL::Core;
   10 use PDL::Exporter;
   11 use DynaLoader;
   12 
   13 
   14    
   15    our @ISA = ( 'PDL::Exporter','DynaLoader' );
   16    push @PDL::Core::PP, __PACKAGE__;
   17    bootstrap PDL::ImageND ;
   18 
   19 
   20 
   21 
   22 
   23 
   24 #line 4 "imagend.pd"
   25 
   26 
   27 =head1 NAME
   28 
   29 PDL::ImageND - useful image processing in N dimensions
   30 
   31 =head1 DESCRIPTION
   32 
   33 These routines act on PDLs as N-dimensional objects, not as broadcasted
   34 sets of 0-D or 1-D objects.  The file is sort of a catch-all for 
   35 broadly functional routines, most of which could legitimately
   36 be filed elsewhere (and probably will, one day).  
   37 
   38 ImageND is not a part of the PDL core (v2.4) and hence must be explicitly
   39 loaded.
   40 
   41 =head1 SYNOPSIS
   42 
   43  use PDL::ImageND;
   44 
   45  $y = $x->convolveND($kernel,{bound=>'periodic'});
   46  $y = $x->rebin(50,30,10);
   47  
   48 =cut
   49 
   50 use strict;
   51 use warnings;
   52 #line 53 "ImageND.pm"
   53 
   54 
   55 
   56 
   57 
   58 
   59 =head1 FUNCTIONS
   60 
   61 =cut
   62 
   63 
   64 
   65 
   66 #line 95 "imagend.pd"
   67 
   68 
   69 use Carp;
   70 #line 71 "ImageND.pm"
   71 
   72 
   73 
   74 #line 948 "../../blib/lib/PDL/PP.pm"
   75 
   76 
   77 
   78 =head2 convolve
   79 
   80 =for sig
   81 
   82   Signature: (a(m); b(n); indx adims(p); indx bdims(q); [o]c(m))
   83 
   84 =for ref
   85 
   86 N-dimensional convolution (Deprecated; use convolveND)
   87 
   88 =for usage
   89 
   90   $new = convolve $x, $kernel
   91 
   92 Convolve an array with a kernel, both of which are N-dimensional.  This 
   93 routine does direct convolution (by copying) but uses quasi-periodic
   94 boundary conditions: each dim "wraps around" to the next higher row in
   95 the next dim.  
   96 
   97 This routine is kept for backwards compatibility with earlier scripts; 
   98 for most purposes you want L<convolveND|PDL::ImageND/convolveND> instead:
   99 it runs faster and handles a variety of boundary conditions.
  100 
  101 
  102 
  103 =for bad
  104 
  105 convolve does not process bad values.
  106 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  107 
  108 
  109 =cut
  110 #line 111 "ImageND.pm"
  111 
  112 
  113 
  114 #line 949 "../../blib/lib/PDL/PP.pm"
  115 
  116 
  117 
  118 # Custom Perl wrapper
  119 
  120 sub PDL::convolve{
  121     my($x,$y,$c) = @_;
  122     barf("Usage: convolve(a(*), b(*), [o]c(*)") if $#_<1 || $#_>2;
  123     $c = PDL->null if $#_<2;
  124     &PDL::_convolve_int( $x->clump(-1), $y->clump(-1),
  125        long([$x->dims]), long([$y->dims]),
  126        ($c->getndims>1? $c->clump(-1) : $c)
  127      );
  128      $c->setdims([$x->dims]);
  129 
  130     if($x->is_inplace) {
  131       $x .= $c;
  132       $x->set_inplace(0);
  133       return $x;
  134     }
  135     return $c;
  136 }
  137 #line 138 "ImageND.pm"
  138 
  139 
  140 
  141 #line 950 "../../blib/lib/PDL/PP.pm"
  142 
  143 *convolve = \&PDL::convolve;
  144 #line 145 "ImageND.pm"
  145 
  146 
  147 
  148 #line 225 "imagend.pd"
  149 
  150 
  151 =head2 ninterpol()
  152 
  153 =for ref
  154 
  155 N-dimensional interpolation routine
  156 
  157 =for sig
  158 
  159  Signature: ninterpol(point(),data(n),[o]value())
  160 
  161 =for usage
  162 
  163       $value = ninterpol($point, $data);
  164 
  165 C<ninterpol> uses C<interpol> to find a linearly interpolated value in
  166 N dimensions, assuming the data is spread on a uniform grid.  To use
  167 an arbitrary grid distribution, need to find the grid-space point from
  168 the indexing scheme, then call C<ninterpol> -- this is far from
  169 trivial (and ill-defined in general).
  170 
  171 See also L<interpND|PDL::Primitive/interpND>, which includes boundary 
  172 conditions and allows you to switch the method of interpolation, but
  173 which runs somewhat slower.
  174 
  175 =cut
  176 
  177 
  178 *ninterpol = \&PDL::ninterpol;
  179 
  180 sub PDL::ninterpol {
  181     use PDL::Math 'floor';
  182     use PDL::Primitive 'interpol';
  183     print 'Usage: $x = ninterpolate($point(s), $data);' if $#_ != 1;
  184     my ($p, $y) = @_;
  185     my ($ip) = floor($p);
  186     # isolate relevant N-cube
  187     $y = $y->slice(join (',',map($_.':'.($_+1),list $ip)));
  188     for (list ($p-$ip)) { $y = interpol($_,$y->xvals,$y); }
  189     $y;
  190 }
  191 #line 192 "ImageND.pm"
  192 
  193 
  194 
  195 #line 948 "../../blib/lib/PDL/PP.pm"
  196 
  197 
  198 
  199 =head2 rebin
  200 
  201 =for sig
  202 
  203   Signature: (a(m); [o]b(n); int ns => n)
  204 
  205 =for ref
  206 
  207 N-dimensional rebinning algorithm
  208 
  209 =for usage
  210 
  211   $new = rebin $x, $dim1, $dim2,..;.
  212   $new = rebin $x, $template;
  213   $new = rebin $x, $template, {Norm => 1};
  214 
  215 Rebin an N-dimensional array to newly specified dimensions.
  216 Specifying `Norm' keeps the sum constant, otherwise the intensities
  217 are kept constant.  If more template dimensions are given than for the
  218 input pdl, these dimensions are created; if less, the final dimensions
  219 are maintained as they were.
  220 
  221 So if C<$x> is a 10 x 10 pdl, then C<rebin($x,15)> is a 15 x 10 pdl,
  222 while C<rebin($x,15,16,17)> is a 15 x 16 x 17 pdl (where the values
  223 along the final dimension are all identical).
  224 
  225 Expansion is performed by sampling; reduction is performed by averaging.
  226 If you want different behavior, use L<PDL::Transform::map|PDL::Transform/map>
  227 instead.  PDL::Transform::map runs slower but is more flexible.
  228 
  229 
  230 
  231 =for bad
  232 
  233 rebin does not process bad values.
  234 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  235 
  236 
  237 =cut
  238 #line 239 "ImageND.pm"
  239 
  240 
  241 
  242 #line 949 "../../blib/lib/PDL/PP.pm"
  243 
  244 
  245 
  246 # Custom Perl wrapper
  247 
  248 sub PDL::rebin {
  249     my($x) = shift;
  250     my($opts) = ref $_[-1] eq "HASH" ? pop : {};
  251     my(@idims) = $x->dims;
  252     my(@odims) = ref $_[0] ? $_[0]->dims : @_;
  253     my($i,$y);
  254     foreach $i (0..$#odims) {
  255       if ($i > $#idims) {  # Just dummy extra dimensions
  256           $x = $x->dummy($i,$odims[$i]);
  257           next;
  258       # rebin_int can cope with all cases, but code
  259       # 1->n and n->1 separately for speed
  260       } elsif ($odims[$i] != $idims[$i]) {       # If something changes
  261          if (!($odims[$i] % $idims[$i])) {      # Cells map 1 -> n
  262                my ($r) = $odims[$i]/$idims[$i];
  263                $y = $x->mv($i,0)->dupN($r);
  264          } elsif (!($idims[$i] % $odims[$i])) { # Cells map n -> 1
  265                my ($r) = $idims[$i]/$odims[$i];
  266                $x = $x->mv($i,0);
  267                # -> copy so won't corrupt input PDL
  268                $y = $x->slice("0:-1:$r")->copy;
  269                foreach (1..$r-1) {
  270                   $y += $x->slice("$_:-1:$r");
  271                }
  272                $y /= $r;
  273          } else {                               # Cells map n -> m
  274              &PDL::_rebin_int($x->mv($i,0), $y = null, $odims[$i]);
  275          }
  276          $x = $y->mv(0,$i);
  277       }
  278     }
  279     if (exists $opts->{Norm} and $opts->{Norm}) {
  280       my ($norm) = 1;
  281       for $i (0..$#odims) {
  282          if ($i > $#idims) {
  283               $norm /= $odims[$i];
  284          } else {
  285               $norm *= $idims[$i]/$odims[$i];
  286          }
  287       }
  288       return $x * $norm;
  289     } else {
  290       # Explicit copy so i) can't corrupt input PDL through this link
  291       #                 ii) don't waste space on invisible elements
  292       return $x -> copy;
  293     }
  294 }
  295 #line 296 "ImageND.pm"
  296 
  297 
  298 
  299 #line 950 "../../blib/lib/PDL/PP.pm"
  300 
  301 *rebin = \&PDL::rebin;
  302 #line 303 "ImageND.pm"
  303 
  304 
  305 
  306 #line 378 "imagend.pd"
  307 
  308 
  309 =head2 circ_mean_p
  310 
  311 =for ref
  312 
  313 Calculates the circular mean of an n-dim image and returns
  314 the projection. Optionally takes the center to be used.
  315 
  316 =for usage
  317 
  318    $cmean=circ_mean_p($im);
  319    $cmean=circ_mean_p($im,{Center => [10,10]});
  320 
  321 =cut
  322 
  323 
  324 sub circ_mean_p {
  325  my ($x,$opt) = @_;
  326  my ($rad,$sum,$norm);
  327 
  328  if (defined $opt) {
  329    $rad = long PDL::rvals($x,$opt);
  330  }
  331  else {
  332    $rad = long rvals $x;
  333  }
  334  $sum = zeroes($rad->max+1);
  335  PDL::indadd $x->clump(-1), $rad->clump(-1), $sum; # this does the real work
  336  $norm = zeroes($rad->max+1);
  337  PDL::indadd pdl(1), $rad->clump(-1), $norm;       # equivalent to get norm
  338  $sum /= $norm;
  339  return $sum;
  340 }
  341 
  342 =head2 circ_mean
  343 
  344 =for ref
  345 
  346 Smooths an image by applying circular mean.
  347 Optionally takes the center to be used.
  348 
  349 =for usage
  350 
  351    circ_mean($im);
  352    circ_mean($im,{Center => [10,10]});
  353 
  354 =cut
  355 
  356 
  357 sub circ_mean {
  358  my ($x,$opt) = @_;
  359  my ($rad,$sum,$norm,$a1);
  360 
  361  if (defined $opt) {
  362    $rad = long PDL::rvals($x,$opt);
  363  }
  364  else {
  365    $rad = long rvals $x;
  366  }
  367  $sum = zeroes($rad->max+1);
  368  PDL::indadd $x->clump(-1), $rad->clump(-1), $sum; # this does the real work
  369  $norm = zeroes($rad->max+1);
  370  PDL::indadd pdl(1), $rad->clump(-1), $norm;       # equivalent to get norm
  371  $sum /= $norm;
  372  $a1 = $x->clump(-1);
  373  $a1 .= $sum->index($rad->clump(-1));
  374 
  375  return $x;
  376 }
  377 #line 378 "ImageND.pm"
  378 
  379 
  380 
  381 #line 454 "imagend.pd"
  382 
  383 
  384 =head2 kernctr
  385 
  386 =for ref
  387 
  388 `centre' a kernel (auxiliary routine to fftconvolve)
  389 
  390 =for usage
  391 
  392     $kernel = kernctr($image,$smallk);
  393     fftconvolve($image,$kernel);
  394 
  395 kernctr centres a small kernel to emulate the behaviour of the direct
  396 convolution routines.
  397 
  398 =cut
  399 
  400 
  401 *kernctr = \&PDL::kernctr;
  402 
  403 sub PDL::kernctr {
  404     # `centre' the kernel, to match kernel & image sizes and
  405     # emulate convolve/conv2d.  FIX: implement with phase shifts
  406     # in fftconvolve, with option tag
  407     barf "Must have image & kernel for kernctr" if $#_ != 1;
  408     my ($imag, $kern) = @_;
  409     my (@ni) = $imag->dims;
  410     my (@nk) = $kern->dims;
  411     barf "Kernel and image must have same number of dims" if $#ni != $#nk;
  412     my ($newk) = zeroes(double,@ni);
  413     my ($k,$n,$y,$d,$i,@stri,@strk,@b);
  414     for ($i=0; $i <= $#ni; $i++) {
  415     $k = $nk[$i];
  416     $n = $ni[$i];
  417     barf "Kernel must be smaller than image in all dims" if ($n < $k);
  418     $d = int(($k-1)/2);
  419         $stri[$i][0] = "0:$d,";
  420         $strk[$i][0] = (-$d-1).":-1,";
  421         $stri[$i][1] = $d == 0 ? '' : ($d-$k+1).':-1,';
  422         $strk[$i][1] = $d == 0 ? '' : '0:'.($k-$d-2).',';
  423     }
  424     # kernel is split between the 2^n corners of the cube
  425     my ($nchunk) = 2 << $#ni;
  426     CHUNK:
  427       for ($i=0; $i < $nchunk; $i++) {
  428     my ($stri,$strk);
  429     for ($n=0, $y=$i; $n <= $#ni; $n++, $y >>= 1) {
  430         next CHUNK if $stri[$n][$y & 1] eq '';
  431       $stri .= $stri[$n][$y & 1];
  432       $strk .= $strk[$n][$y & 1];
  433     }
  434     chop ($stri); chop ($strk);
  435     (my $t = $newk->slice($stri)) .= $kern->slice($strk);
  436     }
  437     $newk;
  438 }
  439 #line 440 "ImageND.pm"
  440 
  441 
  442 
  443 #line 948 "../../blib/lib/PDL/PP.pm"
  444 
  445 
  446 
  447 =head2 convolveND
  448 
  449 =for sig
  450 
  451   Signature: (k0(); SV *k; SV *aa; SV *a)
  452 
  453 
  454 =for ref
  455 
  456 Speed-optimized convolution with selectable boundary conditions
  457 
  458 =for usage
  459 
  460   $new = convolveND($x, $kernel, [ {options} ]);
  461 
  462 Conolve an array with a kernel, both of which are N-dimensional.
  463 
  464 If the kernel has fewer dimensions than the array, then the extra array
  465 dimensions are broadcasted over.  There are options that control the boundary
  466 conditions and method used.
  467 
  468 The kernel's origin is taken to be at the kernel's center.  If your
  469 kernel has a dimension of even order then the origin's coordinates get
  470 rounded up to the next higher pixel (e.g. (1,2) for a 3x4 kernel).
  471 This mimics the behavior of the earlier L</convolve> and
  472 L<fftconvolve|PDL::FFT/fftconvolve()> routines, so convolveND is a drop-in
  473 replacement for them.
  474 
  475 
  476 The kernel may be any size compared to the image, in any dimension.
  477 
  478 The kernel and the array are not quite interchangeable (as in mathematical
  479 convolution): the code is inplace-aware only for the array itself, and
  480 the only allowed boundary condition on the kernel is truncation.
  481 
  482 convolveND is inplace-aware: say C<convolveND(inplace $x ,$k)> to modify
  483 a variable in-place.  You don't reduce the working memory that way -- only
  484 the final memory.
  485 
  486 OPTIONS
  487 
  488 Options are parsed by PDL::Options, so unique abbreviations are accepted.
  489 
  490 =over 3
  491 
  492 =item boundary (default: 'truncate')
  493 
  494 The boundary condition on the array, which affects any pixel closer
  495 to the edge than the half-width of the kernel.  
  496 
  497 The boundary conditions are the same as those accepted by
  498 L<range|PDL::Slices/range>, because this option is passed directly
  499 into L<range|PDL::Slices/range>.  Useful options are 'truncate' (the
  500 default), 'extend', and 'periodic'.  You can select different boundary 
  501 conditions for different axes -- see L<range|PDL::Slices/range> for more 
  502 detail.
  503 
  504 The (default) truncate option marks all the near-boundary pixels as BAD if
  505 you have bad values compiled into your PDL and the array's badflag is set. 
  506 
  507 =item method (default: 'auto')
  508 
  509 The method to use for the convolution.  Acceptable alternatives are
  510 'direct', 'fft', or 'auto'.  The direct method is an explicit
  511 copy-and-multiply operation; the fft method takes the Fourier
  512 transform of the input and output kernels.  The two methods give the
  513 same answer to within double-precision numerical roundoff.  The fft
  514 method is much faster for large kernels; the direct method is faster
  515 for tiny kernels.  The tradeoff occurs when the array has about 400x
  516 more pixels than the kernel.
  517 
  518 The default method is 'auto', which chooses direct or fft convolution
  519 based on the size of the input arrays.
  520 
  521 =back
  522 
  523 NOTES
  524 
  525 At the moment there's no way to broadcast over kernels.  That could/should
  526 be fixed.
  527 
  528 The broadcasting over input is cheesy and should probably be fixed:
  529 currently the kernel just gets dummy dimensions added to it to match
  530 the input dims.  That does the right thing tersely but probably runs slower
  531 than a dedicated broadcastloop.
  532 
  533 The direct copying code uses PP primarily for the generic typing: it includes
  534 its own broadcastloops.
  535 
  536 
  537 
  538 =for bad
  539 
  540 convolveND does not process bad values.
  541 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  542 
  543 
  544 =cut
  545 #line 546 "ImageND.pm"
  546 
  547 
  548 
  549 #line 949 "../../blib/lib/PDL/PP.pm"
  550 
  551 
  552 use PDL::Options;
  553 
  554 # Perl wrapper conditions the data to make life easier for the PP sub.
  555 
  556 sub PDL::convolveND {
  557   my($a0,$k,$opt0) = @_;
  558   my $inplace = $a0->is_inplace;
  559   my $x = $a0->new_or_inplace;
  560 
  561  
  562   barf("convolveND: kernel (".join("x",$k->dims).") has more dims than source (".join("x",$x->dims).")\n")
  563     if($x->ndims < $k->ndims);
  564   
  565 
  566   # Coerce stuff all into the same type.  Try to make sense.
  567   # The trivial conversion leaves dataflow intact (nontrivial conversions
  568   # don't), so the inplace code is OK.  Non-inplace code: let the existing
  569   # PDL code choose what type is best.
  570   my $type;
  571   if($inplace) {
  572     $type = $a0->get_datatype;
  573   } else {
  574     my $z = $x->flat->index(0) + $k->flat->index(0);
  575     $type = $z->get_datatype;
  576   }
  577   $x = $x->convert($type);
  578   $k = $k->convert($type);
  579     
  580 
  581   ## Handle options -- $def is a static variable so it only gets set up once.
  582   our $def;
  583   unless(defined($def)) {
  584     $def = new PDL::Options( {
  585                               Method=>'a',
  586                               Boundary=>'t'
  587                              }
  588                  );
  589     $def->minmatch(1);
  590     $def->casesens(0);
  591   }
  592 
  593   my $opt = $def->options(PDL::Options::ifhref($opt0));
  594 
  595   ### 
  596   # If the kernel has too few dimensions, we broadcast over the other
  597   # dims -- this is the same as supplying the kernel with dummy dims of
  598   # order 1, so, er, we do that.
  599   $k = $k->dummy($x->dims - 1, 1)
  600     if($x->ndims > $k->ndims);
  601   my $kdims = pdl($k->dims); 
  602 
  603   ###
  604   # Decide whether to FFT or directly convolve: if we're in auto mode,
  605   # choose based on the relative size of the image and kernel arrays.
  606   my $fft = ( ($opt->{Method} =~ m/^a/i) ?
  607            ( $x->nelem > 2500 and ($x->nelem) <= ($k->nelem * 500) ) :
  608            ( $opt->{Method} !~ m/^[ds]/i )
  609           );
  610 
  611   ###
  612   # Pad the array to include boundary conditions
  613   my $adims = pdl($x->dims);
  614   my $koff = ($kdims/2)->ceil - 1;
  615 
  616   my $aa = $x->range( -$koff, $adims + $kdims, $opt->{Boundary} )
  617                ->sever;
  618 
  619   if($fft) {
  620     require PDL::FFT;
  621 
  622     print "convolveND: using FFT method\n" if($PDL::debug);
  623 
  624     # FFT works best on doubles; do our work there then cast back
  625     # at the end.  
  626     $aa = double($aa);
  627     my $aai = $aa->zeroes;
  628 
  629     my $kk = $aa->zeroes;
  630     my $kki = $aa->zeroes;
  631     my $tmp;  # work around new perl -d "feature"
  632     ($tmp = $kk->range( - ($kdims/2)->floor, $kdims, 'p')) .= $k;
  633     PDL::fftnd($kk, $kki);
  634     PDL::fftnd($aa, $aai);
  635 
  636     {
  637       my($ii) = $kk * $aai   +    $aa * $kki;
  638       $aa =     $aa * $kk    -   $kki * $aai;
  639       $aai .= $ii;
  640     }
  641 
  642     PDL::ifftnd($aa,$aai);
  643     $x .= $aa->range( $koff, $adims);
  644 
  645   } else {
  646     print "convolveND: using direct method\n" if($PDL::debug);
  647 
  648     ### The first argument is a dummy to set $GENERIC.  
  649     &PDL::_convolveND_int( $k->flat->index(0), $k, $aa, $x );
  650 
  651   }
  652 
  653 
  654   $x;
  655 }
  656 #line 657 "ImageND.pm"
  657 
  658 
  659 
  660 #line 950 "../../blib/lib/PDL/PP.pm"
  661 
  662 *convolveND = \&PDL::convolveND;
  663 #line 664 "ImageND.pm"
  664 
  665 
  666 
  667 
  668 
  669 #line 34 "imagend.pd"
  670 
  671 
  672 =head1 AUTHORS
  673 
  674 Copyright (C) Karl Glazebrook and Craig DeForest, 1997, 2003
  675 All rights reserved. There is no warranty. You are allowed
  676 to redistribute this software / documentation under certain
  677 conditions. For details, see the file COPYING in the PDL
  678 distribution. If this file is separated from the PDL distribution,
  679 the copyright notice should be included in the file.
  680 
  681 =cut
  682 #line 683 "ImageND.pm"
  683 
  684 
  685 
  686 
  687 # Exit with OK status
  688 
  689 1;