Christian Soeller > PDL-2.4.1 > PDL::Image2D

Download:
PDL-2.4.1.tar.gz

Annotate this POD

CPAN RT

Open  0
View/Report Bugs
Source   Latest Release: PDL-2.007_05

NAME ^

PDL::Image2D - Miscellaneous 2D image processing functions

DESCRIPTION ^

Miscellaneous 2D image processing functions - for want of anywhere else to put them.

SYNOPSIS ^

 use PDL::Image2D;

AUTHORS ^

Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams (rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu), and Doug Burke (burke@ifa.hawaii.edu).

All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.

2D convolution of an array with a kernel (smoothing)

For large kernels, using a FFT routine, such as fftconvolve() in PDL::FFT, will be quicker.

 $new = conv2d $old, $kernel, {OPTIONS}
 $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}
 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
            => Default  - periodic boundary conditions 
                          (i.e. wrap around axis)
            => Reflect  - reflect at boundary
            => Truncate - truncate at boundary

EOD BadDoc => 'Unlike the FFT routines, conv2d is able to process bad values.', HandleBad => 1, Pars => 'a(m,n); kern(p,q); [o]b(m,n);', OtherPars => 'int opt;', PMCode => '

sub PDL::conv2d { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' if $#_<1 || $#_>2; my($a,$kern) = @_; my $c = $#_ == 2 ? $_[2] : $a->nullcreate; &PDL::_conv2d_int($a,$kern,$c, (!(defined $opt && exists $$opt{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; }

', Code => init_vars( { vars => 'PDL_Double tmp;' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j<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_map("i") . 
           init_map("j") .
'
           threadloop %{
           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) {
                             if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(kern(p=>i1,q=>j1)) ) {
                                tmp += $a(m=>i2,n=>j2) * $kern(p=>i1,q=>j1);
                                flag = 1;
                             } /* if: good */
                          } /* if: i2 >= 0 */
                       } /* for: i1 */
                    } /* if: j2 >= 0 */
                 } /* for: j1 */
                 if ( flag ) { $b(m=>i,n=>j) = tmp; }
                 else        { $SETBAD(b(m=>i,n=>j)); }
              } /* for: i */
           } /* for: j */
           %}
           free(mapj+1-q_size); free(mapi+1-p_size);',

); # pp_def: conv2d

pp_def('med2d', Doc=> <<'EOD', =for ref

2D median-convolution of an array with a kernel (smoothing)

Note: only points in the kernel >0 are included in the median, other points are weighted by the kernel value (medianing lots of zeroes is rather pointless)

 $new = med2d $old, $kernel, {OPTIONS}
 $smoothed = med2d $image, ones(3,3), {Boundary => Reflect}
 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
            => Default  - periodic boundary conditions (i.e. wrap around axis)
            => Reflect  - reflect at boundary
            => Truncate - truncate at boundary

EOD BadDoc => 'Bad values are ignored in the calculation. If all elements within the kernel are bad, the output is set bad.', HandleBad => 1, Pars => 'a(m,n); kern(p,q); [o]b(m,n);', OtherPars => 'int opt;', PMCode => '

sub PDL::med2d { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\' if $#_<1 || $#_>2; my($a,$kern) = @_; croak "med2d: kernel must contain some positive elements.\n" if all( $kern <= 0 ); my $c = $#_ == 2 ? $_[2] : $a->nullcreate; &PDL::_med2d_int($a,$kern,$c, (!(defined $opt && exists $opt->{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; }

', Code => init_vars( { vars => 'PDL_Double *tmp, kk; int count;', malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));', check => '(tmp==NULL) || ' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j<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 */

                 PDL->qsort_D( tmp, 0, count-1 );
                 $b(m=>i,n=>j) = tmp[(count-1)/2];

              } /* for: i */
           } /* for: j */
           %}
           free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
',
        BadCode => 
           init_vars( { vars => 'PDL_Double *tmp, kk, aa; int count, flag;',
                        malloc => 'tmp = malloc(p_size*q_size*sizeof(PDL_Double));',
                        check => '(tmp==NULL) || ' } ) .
           init_map("i") . 
           init_map("j") .
'
           threadloop %{
           for(j=0; j<n_size; j++) { 
              for(i=0; i<m_size; i++) {
                 count = 0;
                 flag = 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);
                             aa = $a(m=>i2,n=>j2);
                             if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(aa,a) ) {
                                flag = 1;
                                if ( kk > 0 ) {
                                   tmp[count++] = aa * kk;
                                }
                             }
                          } /* if: i2 >= 0 */
                       } /* for: i1 */
                 } /* for: j1 */
                 if ( flag == 0 ) {
                    $SETBAD(b(m=>i,n=>j));
                 } else {

                    PDL->qsort_D( tmp, 0, count-1 );
                    $b(m=>i,n=>j) = tmp[(count-1)/2];

                 }
              } /* for: i */
           } /* for: j */
           %}
           free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
'

); # pp_def: med2d

pp_def('med2df', Doc=> <<'EOD', =for ref

2D median-convolution of an array in a pxq window (smoothing)

Note: this routine does the median over all points in a rectangular window and is not quite as flexible as med2d in this regard but slightly faster instead

 $new = med2df $old, $xwidth, $ywidth, {OPTIONS}
 $smoothed = med2df $image, 3, 3, {Boundary => Reflect}
 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
            => Default  - periodic boundary conditions (i.e. wrap around axis)
            => Reflect  - reflect at boundary
            => Truncate - truncate at boundary

EOD Pars => 'a(m,n); [o]b(m,n);', # funny parameter names to avoid special case in 'init_vars' OtherPars => 'int __p_size; int __q_size; int opt;', PMCode => '

sub PDL::med2df { my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\'; die \'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )\' if $#_<2 || $#_>3; my($a,$p,$q) = @_; croak "med2df: kernel must contain some positive elements.\n" if $p == 0 && $q == 0; my $c = $#_ == 3 ? $_[3] : $a->nullcreate; &PDL::_med2df_int($a,$c,$p,$q, (!(defined $opt && exists $opt->{Boundary}))?0: (($$opt{Boundary} eq "Reflect") + 2*($$opt{Boundary} eq "Truncate"))); return $c; }

', Code => init_vars( { vars => '$GENERIC() *tmp, kk; int count;', malloc => 'tmp = malloc(p_size*q_size*sizeof($GENERIC()));', check => '(tmp==NULL) || ' } ) . init_map("i") . init_map("j") . ' threadloop %{ for(j=0; j<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) {
                                tmp[count++] = $a(m=>i2,n=>j2);
                          } /* if: i2 >= 0 */
                       } /* for: i1 */
                 } /* for: j1 */

                 $b(m=>i,n=>j) =
                        quick_select_$TBSULQFD(B,S,U,L,Q,F,D) (tmp, count );

              } /* for: i */
           } /* for: j */
           %}
           free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
',

); # pp_def: med2df

pp_addhdr(<<'EOH'); #define EZ(x) ez ? 0 : (x) EOH pp_def('box2d', Pars => 'a(n,m); [o] b(n,m)', OtherPars => 'int wx; int wy; int edgezero', Code => ' register int nx = 0.5*$COMP(wx); register int ny = 0.5*$COMP(wy); register int xs = $SIZE(n); register int ys = $SIZE(m); register int ez = $COMP(edgezero); double div, sum, lsum; int xx,yy,y,ind1,ind2,first;

                div = 1/((2.0*nx+1)*(2.0*ny+1));
                
                threadloop %{
                  first = 1;
                  for (y=0;y<ys;y++)
                    for (xx=0; xx<nx; xx++) {
                        ind1 = xs-1-xx;  /* rightmost strip */
                        $b(n=>xx,m=>y) = EZ($a(n=>xx,m=>y));
                        $b(n=>ind1,m=>y) = EZ($a(n=>ind1,m=>y));
                  }
                  for (xx=0;xx<xs;xx++)
                    for (y=0; y<ny; y++) {
                        ind1 = ys-1-y;  /* topmost strip */
                        $b(n=>xx,m=>y) = EZ($a(n=>xx,m=>y));
                        $b(n=>xx,m=>ind1) = EZ($a(n=>xx,m=>ind1));
                  }
                  for (y=ny;y<ys-ny;y++) {
                    if (first) {
                       first = 0;
                       lsum = 0;
                       for (yy=y-ny;yy<=y+ny;yy++)   /* initialize sum */
                          for (xx=0;xx<=2*nx;xx++)
                            lsum += $a(n=>xx,m=>yy);
                    } else {
                        ind1 = y-ny-1;
                        ind2 = y+ny;
                        for (xx=0;xx<=2*nx;xx++) {
                            lsum -= $a(n=>xx,m=>ind1); /* remove top pixels */
                            lsum += $a(n=>xx,m=>ind2); /* add bottom pixels */
                        }
                    }
                    sum = lsum;
                    $b(n=>nx,m=>y) = div*sum;     /* and assign */
                    for (xx=nx+1;xx<xs-nx;xx++) {     /* loop along line */
                      ind1 = xx-nx-1;
                      ind2 = xx+nx;
                      for (yy=y-ny;yy<=y+ny;yy++) {
                        sum -= $a(n=>ind1,m=>yy); /* remove leftmost data */
                        sum += $a(n=>ind2,m=>yy); /* and add rightmost */
                      }
                      $b(n=>xx,m=>y) = div*sum;         /* and assign */
                    }
                  }
                %}',
       Doc => << 'EOD',
=for ref

fast 2D boxcar average

  $smoothim = $im->box2d($wx,$wy,$edgezero=1);

The edgezero argument controls if edge is set to zero (edgezero=1) or just keeps the original (unfiltered) values.

box2d should be updated to support similar edge options as conv2d and med2d etc.

Boxcar averaging is a pretty crude way of filtering. For serious stuff better filters are around (e.g., use conv2d with the appropriate kernel). On the other hand it is fast and computational cost grows only approximately linearly with window size. EOD ); # pp_def box2d

pp_def('patch2d', Doc=><<'EOD', =for ref

patch bad pixels out of 2D images using a mask

 $patched = patch2d $data, $bad;

$bad is a 2D mask array where 1=bad pixel 0=good pixel. Pixels are replaced by the average of their non-bad neighbours; if all neighbours are bad, the original data value is copied across.

EOD BadDoc => 'This routine does not handle bad values - use patchbad2d instead', HandleBad => 0, Pars => 'a(m,n); int bad(m,n); [o]b(m,n);', Code => 'int m_size, n_size, i,j, i1,j1, i2,j2, norm; double tmp;

         m_size = $COMP(__m_size); n_size = $COMP(__n_size);

      threadloop %{
         
         for(j=0; j<n_size; j++) { 
            for(i=0; i<m_size; i++) {

               $b(m=>i,n=>j) = $a(m=>i,n=>j);

               if ( $bad(m=>i,n=>j)==1 ) {
                  tmp = 0; norm=0;
                  for(j1=-1; j1<=1; j1++) { 
                     j2 = j+j1;
                     if ( j2>=0 && j2<n_size ) {
                        for(i1=-1; i1<=1; i1++) {
                           /* ignore central pixel, which we know is bad */
                           if ( i1!=0 || j1!=0 ) {
                              i2 = i+i1; 
                              if ( i2>=0 && i2<m_size && $bad(m=>i2,n=>j2)!=1 ) {
                                 tmp += $a(m=>i2,n=>j2); 
                                 norm++;
                              }
                           } /* if: i1!=0 || j1!=0 */
                        } /* for: i1 */
                     }
                  } /* for: j1 */

                  if (norm>0) {  /* Patch */
                     $b(m=>i,n=>j) = tmp/norm;
                  }
         
               } /* if: bad() */

            } /* for: i */
         } /* for: j */

      %} /* threadloop */

         ', # Code
);

pp_def('patchbad2d', Doc=><<'EOD', =for ref

patch bad pixels out of 2D images containing bad values

 $patched = patchbad2d $data;

Pixels are replaced by the average of their non-bad neighbours; if all neighbours are bad, the output is set bad. If the input piddle contains no bad values, then a straight copy is performed (see patch2d).

EOD BadDoc => 'patchbad2d handles bad values. The output piddle may contain bad values, depending on the pattern of bad values in the input piddle.', HandleBad => 1, Pars => 'a(m,n); [o]b(m,n);', Code => 'loop(n,m) %{ $b() = $a(); %}', # just copy CopyBadStatusCode => '', # handled by BadCode BadCode => 'int m_size, n_size, i,j, i1,j1, i2,j2, norm, flag; double tmp; $GENERIC(a) a_val;

         flag = 0;
         m_size = $COMP(__m_size); n_size = $COMP(__n_size);

      threadloop %{
         
         for(j=0; j<n_size; j++) { 
            for(i=0; i<m_size; i++) {

               a_val = $a(m=>i,n=>j);   
               if ( $ISGOODVAR(a_val,a) ) {
                  $b(m=>i,n=>j) = a_val;

               } else { 
                  tmp = 0; norm=0;
                  for(j1=-1; j1<=1; j1++) { 
                     j2 = j+j1;
                     if ( j2>=0 && j2<n_size ) {
                        for(i1=-1; i1<=1; i1++) {
                           /* ignore central pixel, which we know is bad */
                           if ( i1!=0 || j1!=0 ) {
                              i2 = i+i1; 
                              if ( i2>=0 && i2<m_size ) {
                                 a_val = $a(m=>i2,n=>j2);
                                 if ( $ISGOODVAR(a_val,a) ) {
                                    tmp += a_val; 
                                    norm++;
                                 }
                              }
                           } /* if: i1!=0 || j1!=0 */
                        } /* for: i1 */
                     }
                  } /* for: j1 */

                  /* Patch */
                  if (norm>0) {
                     $b(m=>i,n=>j) = tmp/norm;
                  } else {
                     $SETBAD(b(m=>i,n=>j));
                     flag = 1;
                  }

               } /* if: ISGOODVAR() */

            } /* for: i */
         } /* for: j */

      %} /* threadloop */

         /* handle bad flag */
         if ( flag ) $PDLSTATESETBAD(b);
         ', # BadCode
);

pp_def('max2d_ind', Doc=><<'EOD', =for ref

Return value/position of maximum value in 2D image

Contributed by Tim Jeness

EOD

      BadDoc=><<'EOD',

Bad values are excluded from the search. If all pixels are bad then the output is set bad.

EOD

       HandleBad => 1,
        Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();',
        Code => '
        double cur; int curind1; int curind2;
        curind1=0;
        curind2=0;
        loop(m) %{
           loop(n) %{
           if((!m && !n) || $a() > cur || IsNaN(cur)) {
                cur = $a(); curind1 = m; curind2 = n;
              }
           %}
        %}
        $val() = cur;
        $x()   = curind1;
        $y()   = curind2;
        ',
        BadCode => '
        double cur; int curind1; int curind2;
        curind1 = -1;
        curind2 = -1;
        loop(m) %{
           loop(n) %{
           if( $ISGOOD(a()) && ( (!n && !m) || ($a() > cur) ) ) {
                cur = $a(); curind1 = m; curind2 = n;
              }
           %}
        %}
        if ( curind1 < 0 ) {
          $SETBAD(val());
          $SETBAD(x());
          $SETBAD(y());
        } else {
          $val() = cur;
          $x()   = curind1;
          $y()   = curind2;
        }
        ');

pp_def('centroid2d',

        Doc=><<'EOD',
=for ref

Refine a list of object positions in 2D image by centroiding in a box

$box is the full-width of the box, i.e. the window is +/- $box/2.

EOD

        BadDoc=><<'EOD',
Bad pixels are excluded from the centroid calculation. If all elements are
bad (or the pixel sum is 0 - but why would you be centroiding
something with negatives in...) then the output values are set bad.

EOD

       HandleBad => 1,
        Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();',

        Code => '
   int i,j,i1,i2,j1,j2,m_size,n_size;
   double sum,data,sumx,sumy;

   m_size = $SIZE(m); n_size = $SIZE(n);

   i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;
   i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;
   j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
   j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;

   sum = sumx = sumy = 0;
   for(j=j1; j<=j2; j++) { for(i=i1; i<=i2; i++) {
      data = $im(m=>i,n=>j);
      sum += data;
      sumx += data*i;
      sumy += data*j;
   }}
   $xcen() = sumx/sum;
   $ycen() = sumy/sum;
',

        BadCode => '
   int i,j,i1,i2,j1,j2,m_size,n_size;
   double sum,data,sumx,sumy;

   m_size = $SIZE(m); n_size = $SIZE(n);

   i1 = $x() - $box()/2; i1 = i1<0 ? 0 : i1;
   i2 = $x() + $box()/2; i2 = i2>=m_size ? m_size-1 : i2;
   j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
   j2 = $y() + $box()/2; j2 = j2>=n_size ? n_size-1 : j2;

   sum = sumx = sumy = 0;
   for(j=j1; j<=j2; j++) { 
      for(i=i1; i<=i2; i++) {
         data = $im(m=>i,n=>j);
         if ( $ISGOODVAR(data,im) ) {
            sum += data;
            sumx += data*i;
            sumy += data*j;
         }
      }
   }
   /*
    * if sum == 0 then we will flag as bad -- although it could just mean that
    * there is negative values in the dataset.
    * - should use a better check than != 0.0  ...
    */
   if ( sum != 0.0 ) {
      $xcen() = sumx/sum;
      $ycen() = sumy/sum;
   } else {
      $SETBAD(xcen());
      $SETBAD(ycen());
  }
'
);

pp_addhdr('

/* Add an equivalence to a list - used by pdl_cc8compt */

void AddEquiv ( PDL_Long* equiv, PDL_Long i, PDL_Long j) {

   PDL_Long k, tmp;

   if (i==j)
      return;

    k = j;
    do {
      k = equiv[k];
    } while ( k != j && k != i );

    if ( k == j ) {
       tmp = equiv[i];
       equiv[i] = equiv[j];
       equiv[j] = tmp;
    }
}

');

pp_def('cc8compt',Doc=>' =for ref

Connected 8-component labeling of a binary image.

Connected 8-component labeling of 0,1 image - i.e. find seperate segmented objects and fill object pixels with object number

 $segmented = cc8compt( $image > $threshold );

', HandleBad => 0, # a marker Pars => 'a(m,n); [o]b(m,n);', Code => '

      PDL_Long i,j,k;
      PDL_Long newlabel;
      PDL_Long neighbour[4];
      PDL_Long nfound;
      PDL_Long pass,count,next,this;
      PDL_Long *equiv;
      PDL_Long i1,j1,i2;
      PDL_Long nx = $SIZE(m);
      PDL_Long ny = $SIZE(n);

      loop(n) %{ loop(m) %{ /* Copy */
         $b() = $a();
      %} %}

      /* 1st pass counts max possible compts, 2nd records equivalences */

      for (pass = 0; pass<2; pass++) {

      if (pass==1) {
         equiv = (PDL_Long*) malloc((newlabel+1)*sizeof(PDL_Long));
         if (equiv==(PDL_Long*)0)
            barf("Out of memory");
         for(i=0;i<=newlabel;i++)
             equiv[i]=i;
      }

      newlabel = 1; /* Running label */

      for(j=0; j<ny; j++) { for(i=0; i<nx; i++) { /* Loop over image pixels */

            nfound = 0; /* Number of neighbour >0 */

            i1 = i-1; j1 = j-1; i2 = i+1;

            if ($b(m=>i, n=>j) > 0) { /* Check 4 neighbour already seen */

               if (i>0 && $b(m=>i1, n=>j)>0)
                   neighbour[nfound++] = $b(m=>i1, n=>j); /* Store label of it */
               if (j>0 && $b(m=>i, n=>j1)>0)
                   neighbour[nfound++] = $b(m=>i, n=>j1);
               if (j>0 && i>0  && $b(m=>i1, n=>j1)>0)
                   neighbour[nfound++] = $b(m=>i1, n=>j1);
               if (j>0 && i<(nx-1) && $b(m=>i2, n=>j1)>0)
                   neighbour[nfound++] = $b(m=>i2, n=>j1);

               if (nfound==0)  { /* Assign new label */
                  $b(m=>i, n=>j) = newlabel++;
               }
               else {
                  $b(m=>i, n=>j) =  neighbour[0];
                  if (nfound>1 && pass == 1) {  /* Assign equivalents */
                      for(k=1; k<nfound; k++)
                         AddEquiv( equiv, (PDL_Long)$b(m=>i, n=>j),
                            neighbour[k] );
                  }
               }
            }

            else {  /* No label */

               $b(m=>i, n=>j) = 0;
            }

      }} /* End of image loop */

      } /* Passes */

      /* Replace each cycle by single label */

       count = 0;
       for (i = 1; i <= newlabel; i++)
         if ( i <= equiv[i] ) {
             count++;
             this = i;
             while ( equiv[this] != i ) {
               next = equiv[this];
               equiv[this] = count;
               this = next;
             }
          equiv[this] = count;
         }


      /* Now remove equivalences */

      for(j=0; j<ny; j++) { for(i=0; i<nx; i++) { /* Loop over image pixels */
           $b(m=>i, n=>j)   = equiv[ (PDL_Long) $b(m=>i, n=>j)  ] ;
      }}

      free(equiv); /* Tidy */
');

pp_addhdr(' #define MAXSEC 32 #define line(x1, x2, y) for (k=x1;k<=x2;k++) \ { /* printf("line from %d to %d\n",x1,x2); */ \ image[k+wx*y] = col; } #define PX(n) ps[2*n] #define PY(n) ps[2*n+1]

void polyfill(PDL_Long *image, int wx, int wy, float *ps, int n, PDL_Long col, int *ierr) { int ymin, ymax, xmin, xmax, fwrd = 1, i, j, k, nsect; int x[MAXSEC], temp, l; float s1, s2, t1, t2;

   ymin = PY(0); ymax = PY(0);
   xmin = PX(0); xmax = PX(0);
   *ierr = 0;
   
   for (i=1; i<n; i++) {
     ymin = ymin > PY(i) ? PY(i) : ymin;        
     ymax = ymax < PY(i) ? PY(i) : ymax;
     xmin = xmin > PX(i) ? PX(i) : xmin;        
     xmax = xmax < PX(i) ? PX(i) : xmax;
   }
   if (xmin < 0 || xmax >= wx || ymin < 0 || ymax >= wy) {
        *ierr = 1; /* clipping */
        return;
   }
   s1 = PX(n-1);
   t1 = PY(n-1);
   for (l=ymin; l<= ymax; l++) {
        nsect = 0;
        fwrd = 1;
        for (i=0; i<n; i++) {
          s2 = PX(i);
          t2 = PY(i);
          if ((t1 < l &&  l <= t2) || (t1 >= l && l > t2)) {
                if (nsect > MAXSEC) {
                        *ierr = 2; /* too complex */
                        return;
                }
                x[nsect] = (s1+(s2-s1)*((l-t1)/(t2-t1)));
                nsect += 1;
          }
          s1 = s2;
          t1 = t2;
        }
        /* sort the intersections */
        for (i=1; i<nsect; i++)
                for (j=0; j<i; j++)
                        if (x[j] > x[i]) {
                                temp = x[j];
                                x[j] = x[i];
                                x[i] = temp;
                        }
        if (fwrd) {
                for (i=0; i<nsect-1; i += 2)
                        line(x[i],x[i+1],l);
                fwrd = 0;
        } else {
                for (i=nsect-1; i>0; i -= 2)
                        line(x[i-1],x[i],l);
                fwrd = 1;
        }
   }
}

');

pp_def('polyfill', HandleBad => 0, # a marker Pars => 'int [o,nc] im(m,n); float ps(two=2,np); int col()', Code => 'int ierr = 0, nerr; threadloop %{ polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr); ierr = ierr < nerr ? nerr : ierr; %} if (ierr) warn("errors during polygonfilling"); ', Doc => << 'EOD', =for ref

fill the area inside the given polygon with a given colour

This function works inplace, i.e. modifies im.

EOD );

pp_add_exported('', 'polyfillv'); pp_addpm(<<'EOPM');

polyfillv

return the (dataflown) area of an image within a polygon

  # increment intensity in area bounded by $poly
  $im->polyfillv($pol)++; # legal in perl >= 5.6
  # compute average intensity within area bounded by $poly
  $av = $im->polyfillv($poly)->avg;

rotate an image by given angle

  # rotate by 10.5 degrees with antialiasing, set missing values to 7
  $rot = $im->rot2d(10.5,7,1);

This function rotates an image through an angle between -90 and + 90 degrees. Uses/doesn't use antialiasing depending on the aa flag. Pixels outside the rotated image are set to bg.

Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based on "A Fast Algorithm for General Raster Rotation" by Alan Paeth, Graphics Interface '86, pp. 77-81.

Use the rotnewsz function to find out about the dimension of the newly created image

  ($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;

PDL::Transform offers a more general interface to distortions, including rotation, with various types of sampling; but rot2d is faster.

EOD );

pp_def('bilin2d', HandleBad => 0, Pars => 'I(n,m); O(q,p)', Doc=><<'EOD', =for ref

Bilinearly maps the first piddle in the second. The interpolated values are actually added to the second piddle which is supposed to be larger than the first one.

EOD , Code =>' int i,j,ii,jj,ii1,jj1,num; double x,y,dx,dy,y1,y2,y3,y4,t,u,sum;

  if ($SIZE(q)>=$SIZE(n) && $SIZE(p)>=$SIZE(m)) {
    threadloop %{
      dx = ((double) ($SIZE(n)-1)) / ($SIZE(q)-1);
      dy = ((double) ($SIZE(m)-1)) / ($SIZE(p)-1);
      for(i=0,x=0;i<$SIZE(q);i++,x+=dx) {
        for(j=0,y=0;j<$SIZE(p);j++,y+=dy) {
          ii = (int) floor(x);
          if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2;
          jj = (int) floor(y);
          if (jj>=($SIZE(m)-1)) jj = $SIZE(m)-2;
          ii1 = ii+1;
          jj1 = jj+1;
          y1 = $I(n=>ii,m=>jj);
          y2 = $I(n=>ii1,m=>jj);
          y3 = $I(n=>ii1,m=>jj1);
          y4 = $I(n=>ii,m=>jj1);
          t = x-ii;
          u = y-jj;
          $O(q=>i,p=>j) += (1-t)*(1-u)*y1 + t*(1-u)*y2 + t*u*y3 + (1-t)*u*y4;
        }
      }
      %}
  }
  else { 
    barf("the second matrix must be greater than first! (bilin2d)");
  }
');

pp_def('rescale2d', HandleBad => 0, Pars => 'I(m,n); O(p,q)', Doc=><<'EOD', =for ref

The first piddle is rescaled to the dimensions of the second (expanding or meaning values as needed) and then added to it in place. Nothing useful is returned.

If you want photometric accuracy or automatic FITS header metadata tracking, consider using PDL::Transform::map instead: it does these things, at some speed penalty compared to rescale2d.

EOD , Code =>' int ix,iy,ox,oy,i,j,lx,ly,cx,cy,xx,yy,num; double kx,ky,temp;

ix = $SIZE(m); iy = $SIZE(n); ox = $SIZE(p); oy = $SIZE(q);

if(ox >= ix && oy >= iy) { threadloop %{ kx = ((double) (ox)) / (ix); ky = ((double) (oy)) / (iy); lx = 0; for(i=0;i<ix;i++) { ly = 0; for(j=0;j<iy;j++) { cx = rint((i+1)*kx)-1; cy = rint((j+1)*ky)-1; for(xx=lx;xx<=cx;xx++) for(yy=ly;yy<=cy;yy++) { /* fprintf(stderr,"i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ $O(p=>xx,q=>yy) += $I(m=>i,n=>j); } ly = cy + 1; } lx = cx + 1; } %} } else if(ox < ix && oy < iy) { threadloop %{ kx = ((double) (ix)) / (ox); ky = ((double) (iy)) / (oy); lx = 0; for(i=0;i<ox;i++) { ly = 0; for(j=0;j<oy;j++) { cx = rint((i+1)*kx)-1; cy = rint((j+1)*ky)-1; temp = 0.0; num = 0; for(xx=lx;xx<=cx;xx++) for(yy=ly;yy<=cy;yy++) { /* fprintf(stderr,"i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ temp += $I(n=>yy,m=>xx); num++; } $O(p=>i,q=>j) += temp/num; ly = cy + 1; } lx = cx + 1; } %} } else if(ox >= ix && oy < iy) { threadloop %{ kx = ((double) (ox)) / (ix); ky = ((double) (iy)) / (oy); lx = 0; for(i=0;i<ix;i++) { ly = 0; for(j=0;j<oy;j++) { cx = rint((i+1)*kx)-1; cy = rint((j+1)*ky)-1; temp = 0.0; num = 0; for(yy=ly;yy<=cy;yy++) { /* fprintf(stderr,"1 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ temp += $I(n=>yy,m=>i); num++; } for(xx=lx;xx<=cx;xx++) { /* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ $O(p=>xx,q=>j) += temp/num; } ly = cy + 1; } lx = cx + 1; } %} } else if(ox < ix && oy >= iy) { threadloop %{ kx = ((double) (ix)) / (ox); ky = ((double) (oy)) / (iy); lx = 0; for(i=0;i<ox;i++) { ly = 0; for(j=0;j<iy;j++) { cx = rint((i+1)*kx)-1; cy = rint((j+1)*ky)-1; temp = 0.0; num = 0; for(xx=lx;xx<=cx;xx++) { /* fprintf(stderr,"1 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ temp += $I(n=>j,m=>xx); num++; } for(yy=ly;yy<=cy;yy++) { /* fprintf(stderr,"2 i: %d, j: %d, xx: %d, yy: %d\n",i,j,xx,yy); */ $O(p=>i,q=>yy) += temp/num; } ly = cy + 1; } lx = cx + 1; } %} } else barf("I am not supposed to be here, please report the bug to <chri@infis.univ.ts.it>"); ');

# functions to make handling 2D polynomial mappings a bit easier # pp_add_exported('', 'fitwarp2d applywarp2d');

pp_addpm( '

fitwarp2d

Find the best-fit 2D polynomial to describe a coordinate transformation.

  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf. { options } )

Given a set of points in the output plane ($u,$v), find the best-fit (using singular-value decomposition) 2D polynomial to describe the mapping back to the image plane ($x,$y). The order of the fit is controlled by the $nf parameter (the maximum power of the polynomial is $nf - 1), and you can restrict the terms to fit using the FIT option.

$px and $py are np by np element piddles which describe a polynomial mapping (of order np-1) from the output (u,v) image to the input (x,y) image:

  x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
  y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j

The transformation is returned for the reverse direction (ie output to input image) since that is what is required by the warp2d() routine. The applywarp2d() routine can be used to convert a set of $u,$v points given $px and $py.

Options:

  FIT     - which terms to fit? default ones(byte,$nf,$nf)
  THRESH  - in svd, remove terms smaller than THRESH * max value
            default is 1.0e-5
FIT

FIT allows you to restrict which terms of the polynomial to fit: only those terms for which the FIT piddle evaluates to true will be evaluated. If a 2D piddle is sent in, then it is used for the x and y polynomials; otherwise $fit->slice(":,:,(0)") will be used for $px and $fit->slice(":,:,(1)") will be used for $py.

THRESH

Remove all singular values whose valus is less than THRESH times the largest singular value.

The number of points must be at least equal to the number of terms to fit ($nf*$nf points for the default value of FIT).

  # points in original image
  $x = pdl( 0,   0, 100, 100 );
  $y = pdl( 0, 100, 100,   0 );
  # get warped to these positions
  $u = pdl( 10, 10, 90, 90 );
  $v = pdl( 10, 90, 90, 10 );
  # 
  # shift of origin + scale x/y axis only
  $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
  print "px = ${px}py = $py";
  px = 
  [
   [-12.5  1.25]
   [    0     0]
  ]
  py = 
  [
   [-12.5     0]
   [ 1.25     0]
  ]
  #
  # Compared to allowing all 4 terms
  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
  print "px = ${px}py = $py";
  px = 
  [
   [         -12.5           1.25]
   [  1.110223e-16 -1.1275703e-17]
  ]
  py = 
  [
   [         -12.5  1.6653345e-16]
   [          1.25 -5.8546917e-18]
  ]

applywarp2d

Transform a set of points using a 2-D polynomial mapping

  ( $x, $y ) = applywarp2d( $px, $py, $u, $v )

Convert a set of points (stored in 1D piddles $u,$v) to $x,$y using the 2-D polynomial with coefficients stored in $px and $py. See fitwarp2d() for more information on the format of $px and $py.

warp2d

  Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { options })

Warp a 2D image given a polynomial describing the reverse mapping.

  $out = warp2d( $img, $px, $py, { options } );

Apply the polynomial transformation encoded in the $px and $py piddles to warp the input image $img into the output image $out.

The format for the polynomial transformation is described in the documentation for the fitwarp2d() routine.

At each point x,y, the closest 16 pixel values are combined with an interpolation kernel to calculate the value at u,v. The interpolation is therefore done in the image, rather than Fourier, domain. By default, a tanh kernel is used, but this can be changed using the KERNEL option discussed below (the choice of kernel depends on the frequency content of the input image).

The routine is based on the warping command from the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and for further details on image resampling see Wolberg, G., "Digital Image Warping", 1990, IEEE Computer Society Press ISBN 0-8186-8944-7).

Currently the output image is the same size as the input one, which means data will be lost if the transformation reduces the pixel scale. This will (hopefully) be changed soon.

  $img = rvals(byte,501,501);
  imag $img, { JUSTIFY => 1 };
  # 
  # use a not-particularly-obvious transformation:
  #   x = -10 + 0.5 * $u - 0.1 * $v 
  #   y = -20 + $v - 0.002 * $u * $v
  #
  $px  = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
  $py  = pdl( [ -20, 0 ], [ 1, 0.002 ] );
  $wrp = warp2d( $img, $px, $py );
  #
  # see the warped image
  imag $warp, { JUSTIFY => 1 };

The options are:

  KERNEL - default value is tanh
  NOVAL  - default value is 0

KERNEL is used to specify which interpolation kernel to use (to see what these kernels look like, use the warp2d_kernel() routine). The options are:

tanh

Hyperbolic tangent: the approximation of an ideal box filter by the product of symmetric tanh functions.

sinc

For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle, which produces a sinc interpolation kernel in the spatial domain:

  sinc(x) = sin(pi * x) / (pi * x)

However, it is not ideal for the 4x4 pixel region used here.

sinc2

This is the square of the sinc function.

lanczos

Although defined differently to the tanh kernel, the result is very similar in the spatial domain. The Lanczos function is defined as

  L(x) = sinc(x) * sinc(x/2)  if abs(x) < 2
       = 0                       otherwise
hann

This kernel is derived from the following function:

  H(x) = a + (1-a) * cos(2*pi*x/(N-1))  if abs(x) < 0.5*(N-1)
       = 0                                 otherwise

with a = 0.5 and N currently equal to 2001.

hamming

This kernel uses the same H(x) as the Hann filter, but with a = 0.54.

NOVAL gives the value used to indicate that a pixel in the output image does not map onto one in the input image.

warp2d_kernel

Return the specified kernel, as used by warp2d

  ( $x, $k ) = warp2d_kernel( $name )

The valid values for $name are the same as the KERNEL option of warp2d().

  line warp2d_kernel( "hamming" );
syntax highlighting: