pp_addpm({At=>Top},<<'EOD');

=head1 NAME

PDL::Image2D - Miscellaneous 2D image processing functions


=head1 DESCRIPTION

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

=head1 SYNOPSIS

 use PDL::Image2D;

=cut

EOD

pp_addhdr('

#define IsNaN(x) (x != x)

/* Fast Modulus with proper negative behaviour */

#define REALMOD(a,b) {while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);}

/* Local quick-sort routine for doubles */

void pdl_lqsortD(double* xx, int a, int b) {

   int i,j;

   double t, median;

   i = a; j = b;
   median = xx[(i+j) / 2];
   do {
      while (xx[i] < median)
	 i++;
      while (median < xx[j])
	 j--;
      if (i <= j) {
	 t = xx[i]; xx[i] = xx[j]; xx[j] = t;
	 i++; j--;
      }
   } while (i <= j);

   if (a < j)
      pdl_lqsortD(xx,a,j);
   if (i < b)
      pdl_lqsortD(xx,i,b);
 }

');

pp_def('conv2d', Doc=><<'EOD',
=head2 conv2d

=for ref

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

=for usage

 $new = conv2d $old, $kernel, {OPTIONS}

=for example

  perldl> $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}

=for options

 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

=cut

EOD
        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] : PDL->null;
   &PDL::_conv2d_int($a,$kern,$c,
	(!(defined $opt && exists $$opt{Boundary}))?0:
	(($$opt{Boundary} eq "Reflect") +
	2*($$opt{Boundary} eq "Truncate")));
   return $c;
}

',
        Code => '
           int i,j, i1,j1, i2,j2, poff, qoff;
           double tmp;
           int opt = $COMP(opt);
           int m_size = $COMP(__m_size);
           int n_size = $COMP(__n_size);
           int p_size = $COMP(__p_size);
           int q_size = $COMP(__q_size);
           int *mapi, *mapj;

           mapi = (int *) malloc((p_size+m_size)*sizeof(int));
           mapj = (int *) malloc((q_size+n_size)*sizeof(int));

           if ((mapi==NULL) || (mapj==NULL))
               barf("Out of Memory");

           poff = p_size/2; mapi += p_size-1;
           qoff = q_size/2; mapj += q_size-1;

	   for (i=1-p_size; i<m_size; i++) {
              i2 = i+poff;
              switch (opt) {
	case 1:      /* REFLECT */
                 if (i2<0)
                    i2 = -i2;
                 else if (i2 >= m_size)
                    i2 = 2*m_size-(i2+1);
                 break;
	case 2:      /* TRUNCATE */
                 if (i2<0 || i2 >= m_size)
                    i2 = -1;
                 break;
	default:
                 REALMOD(i2,m_size);
               }
               mapi[i] = i2;
	   }
	   for (j=1-q_size; j<n_size; j++) {
              j2 = j+qoff;
              switch (opt) {
	case 1:      /* REFLECT */
                 if (j2<0)
                    j2 = -j2;
                 else if (j2 >= n_size)
                    j2 = 2*n_size-(j2+1);
                 break;
	case 2:      /* TRUNCATE */
                 if (j2<0 || j2>=n_size)
                    j2 = -1;
                 break;
	default:
                 REALMOD(j2,n_size);
               }
               mapj[j] = j2;
	   }

           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);
                     }
                  }
              }
              $b(m=>i,n=>j) = tmp;
           }} %}
	   free(mapj+1-q_size); free(mapi+1-p_size);
');


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

=head2 med2d

=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)

=for usage

 $new = med2d $old, $kernel, {OPTIONS}

=for example

  perldl> $smoothed = med2d $image, ones(3,3), {Boundary => Reflect}

=for options

 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

=cut

EOD
        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) = @_;
   my $c = $#_ == 2 ? $_[2] : PDL->null;
   &PDL::_med2d_int($a,$kern,$c,
	(!(defined $opt && exists $opt->{Boundary}))?0:
	(($$opt{Boundary} eq "Reflect") +
	2*($$opt{Boundary} eq "Truncate")));
   return $c;
}

',
        Code => '
           int i,j, i1,j1, i2,j2, poff, qoff,
           count;
           PDL_Double *tmp;
           PDL_Double kk;

           int opt = $COMP(opt);
           int m_size = $COMP(__m_size);
           int n_size = $COMP(__n_size);
           int p_size = $COMP(__p_size);
           int q_size = $COMP(__q_size);
           int *mapi, *mapj;

           poff = p_size/2; qoff = q_size/2;

           tmp = malloc(p_size*q_size*sizeof(PDL_Double));
           mapi = (int *) malloc((p_size+m_size)*sizeof(int));
           mapj = (int *) malloc((q_size+n_size)*sizeof(int));
           if ((tmp==NULL) || (mapi==NULL) || (mapj==NULL))
               barf("Out of Memory");

           poff = p_size/2; mapi += p_size-1;
           qoff = q_size/2; mapj += q_size-1;

	   for (i=1-p_size; i<m_size; i++) {
              i2 = i+poff;
              switch (opt) {
	case 1:      /* REFLECT */
                 if (i2<0)
                    i2 = -i2;
                 else if (i2 >= m_size)
                    i2 = 2*m_size-(i2+1);
                 break;
	case 2:      /* TRUNCATE */
                 if (i2<0 || i2 >= m_size)
                    i2 = -1;
                 break;
	default:
                 REALMOD(i2,m_size);
               }
               mapi[i] = i2;
	   }
	   for (j=1-q_size; j<n_size; j++) {
              j2 = j+qoff;
              switch (opt) {
	case 1:      /* REFLECT */
                 if (j2<0)
                    j2 = -j2;
                 else if (j2 >= n_size)
                    j2 = 2*n_size-(j2+1);
                 break;
	case 2:      /* TRUNCATE */
                 if (j2<0 || j2>=n_size)
                    j2 = -1;
                 break;
	default:
                 REALMOD(j2,n_size);
               }
               mapj[j] = j2;
	   }

           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;
                           }
                        }
              }}
              pdl_lqsortD(tmp,0,count-1);

              $b(m=>i,n=>j) = tmp[(count-1)/2];
           }} %}
           free(mapj+1-q_size); free(mapi+1-p_size); free(tmp);
');

pp_def('patch2d', Doc=><<'EOD',
=head2 patch2d

=for ref

patch bad pixels out of 2D images,

=for usage

$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.

=cut

EOD
        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);

           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++) { for(i1=-1; i1<=1; i1++) {
                     i2 = i+i1; j2 = j+j1;
                     if (j2>0 && i2>0 && i2<m_size && j2<n_size &&
                        $bad(m=>i2,n=>j2)!=1)
                        tmp += $a(m=>i2,n=>j2); norm++;
                 }}
                 if (norm>0) {  /* Patch */
                    $b(m=>i,n=>j) = tmp/norm;
                 }

              } /* Next pixel */
           }}
');


pp_def('max2d_ind',Doc=><<'EOD',
=head2 max2d_ind

=for ref

Return value/position of maximum value in 2D image

Contributed by Tim Jeness

=cut

EOD

        Pars => 'a(m,n); [o]b(); int [o]c(); int[o]d();',
        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;
              }
           %}
        %}
        $b() = cur;
        $c() = curind1;
        $d() = curind2;
        ');

pp_def('centroid2d',Doc=><<'EOD',
=head2 centroid2d

=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

=cut

EOD
	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; i1 = i1>=m_size ? m_size-1 : i1;
   j1 = $y() - $box()/2; j1 = j1<0 ? 0 : j1;
   j2 = $y() + $box()/2; j1 = j1>=n_size ? n_size-1 : j1;

   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;
'
);

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

=for example

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

=cut

',
        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_def('bilin2d',
    Pars => 'I(n,m); O(q,p)',
    Doc=><<'EOD',
=head2 bilin2d

=for ref

Bilineary 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.

=cut

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',
    Pars => 'I(n,m); O(q,p)',
    Doc=><<'EOD',
=head2 rescale2d

=for ref

The first piddle is rescaled to the dimensions of the second
(expandind or meaning values as needed) and then added to it.

=cut

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

ix = $SIZE(n);   
iy = $SIZE(m);   
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(n=>j,m=>i);
	  }
	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>"); 
  ');

pp_addpm({At=>Bot},<<'EOD');

=head1 AUTHORS

Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams
(rjrw@ast.leeds.ac.uk) and Tim Jeness (timj@jach.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.


=cut

EOD

pp_done();