edges-orig.c

matlab的steel金字塔小波分解源代码

	  if (x_edge IS 0) x_edge++;
	  result[ABS(y_base-ABS(y_edge)+x_dim) + ABS(x_base-ABS(x_edge)+1)]
	    += filt[x_filt];
	  }
      }
  return(0);
  }

/* --------------------------------------------------------------------
reflect1() - Reflection through the edge pixels.  This is the right thing
to do if you are subsampling by 2, since it maintains parity (even 
pixels positions remain even, odd ones remain odd). (note: procedure differs 
depending on f_or_e parameter).  */	 

int reflect1(filt,x_dim,y_dim,x_pos,y_pos,result,f_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, f_or_e;
  {
  int filt_sz = x_dim*y_dim;
  register int x_start = 0, y_start = 0, x_stop = x_dim, y_stop = filt_sz;
  register int y_filt,x_filt, y_edge,x_edge;
  register int x_base = (x_pos>0)?(x_dim-1):0;
  register int y_base = (y_pos>0)?(x_dim*(y_dim-1)):0; 
  int x_edge_dist = (x_pos>0)?(x_pos-x_dim):((x_pos<0)?(x_pos+1):0);
  int y_edge_dist = x_dim * ((y_pos>0)?(y_pos-y_dim):((y_pos<0)?(y_pos+1):0));
  int i;
  int mx_pos = (x_dim/2)+1;
  int my_pos = (y_dim/2)+1;

  #ifdef DEBUG
    printf("(%d,%d)  ",y_pos,x_pos);
    if (x_pos==0) printf("\n");
  #endif

  for (i=0; i<filt_sz; i++) result[i] = 0.0;

  /* if EXPAND and filter is centered on image edge, do not reflect */
  if (f_or_e IS EXPAND)
      {
      if (x_pos IS mx_pos) x_stop = (x_dim+1)/2;
      else if (x_pos IS -mx_pos) { x_start = x_dim/2; x_edge_dist = 0; }

      if (y_pos IS my_pos) y_stop = x_dim*((y_dim+1)/2);
      else if (y_pos IS -my_pos) { y_start = x_dim*(y_dim/2); y_edge_dist = 0;}
      }

  /* reflect at boundary of image */
  for (y_filt=y_start, y_edge=y_edge_dist;
       y_filt<y_stop;
       y_filt+=x_dim, y_edge+=x_dim)
    for (x_filt=y_filt+x_start, x_edge=x_edge_dist;
	 x_filt<y_filt+x_stop;
	 x_filt++, x_edge++)
	result[ABS(y_base-ABS(y_edge)) + ABS(x_base-ABS(x_edge))]
	  += filt[x_filt];

  /* if EXPAND and filter is not centered on image edge, mult edge by 2 */
    if (f_or_e IS EXPAND)
      {
      if ( (ABS(x_pos) ISNT mx_pos) AND (x_pos ISNT 0) )
	for (y_filt=x_base; y_filt<filt_sz; y_filt+=x_dim)
	  result[y_filt] += result[y_filt];
      if ( (ABS(y_pos) ISNT my_pos) AND (y_pos ISNT 0) )
	for (x_filt=y_base; x_filt<y_base+x_dim; x_filt++)
	  result[x_filt] += result[x_filt];
      }
  return(0);
  }

/* --------------------------------------------------------------------
Extend() - Extend image by reflecting and inverting about edge pixel
value.  Maintains continuity in intensity AND first derivative (but
not higher derivs).
*/

int Extend(filt,x_dim,y_dim,x_pos,y_pos,result,f_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, f_or_e;
  {
  int filt_sz = x_dim*y_dim;
  register int x_start = 0, y_start = 0, x_stop = x_dim, y_stop = filt_sz;
  register int y_filt,x_filt, y_edge,x_edge;
  register int x_base = (x_pos>0)?(x_dim-1):0;
  register int y_base = (y_pos>0)?(x_dim*(y_dim-1)):0; 
  int x_edge_dist = (x_pos>0)?(x_pos-x_dim):((x_pos<-1)?(x_pos+1):0);
  int y_edge_dist = x_dim * ((y_pos>0)?(y_pos-y_dim):((y_pos<-1)?(y_pos+1):0));
  int i;
  int mx_pos = (x_dim/2)+1;
  int my_pos = (y_dim/2)+1;

  for (i=0; i<filt_sz; i++) result[i] = 0.0;

  /* if EXPAND and filter is centered on image edge, do not reflect */
  if (f_or_e IS EXPAND)
      {
      if (x_pos IS mx_pos) x_stop = (x_dim+1)/2;
      else if (x_pos IS -mx_pos) { x_start = x_dim/2; x_edge_dist = 0; }

      if (y_pos IS my_pos) y_stop = x_dim*((y_dim+1)/2);
      else if (y_pos IS -my_pos) { y_start = x_dim*(y_dim/2); y_edge_dist = 0;}
      }

  /* reflect at boundary of image */
  for (y_filt=y_start, y_edge=y_edge_dist;
       y_filt<y_stop;
       y_filt+=x_dim, y_edge+=x_dim)
    for (x_filt=y_filt+x_start, x_edge=x_edge_dist;
	 x_filt<y_filt+x_stop;
	 x_filt++, x_edge++)
      if (((!y_base AND (sgn(y_edge) IS -1)) /* y overhanging */
	   OR
	   (y_base AND (sgn(y_edge) IS 1)))
	  ISNT			             /* XOR */
	  ((!x_base AND (sgn(x_edge) IS -1)) /* x overhanging */
	   OR
	   (x_base AND (sgn(x_edge) IS 1))))
	  {
	  result[ABS(y_base-ABS(y_edge)) + ABS(x_base-ABS(x_edge))]
	    -= filt[x_filt];
	  result[clip(y_base+y_edge,0,y_dim) + clip(x_base+x_edge,0,x_dim)]
	    += filt[x_filt] + filt[x_filt];
	  }
      else result[ABS(y_base-ABS(y_edge)) + ABS(x_base-ABS(x_edge))]
	  += filt[x_filt];
  return(0);
  }

/* --------------------------------------------------------------------
predict() - Simple prediction.  Like zero, but multiplies the result
by the reciprocal of the percentage of filter being used.  (i.e. if
50% of the filter is hanging over the edge of the image, multiply the
taps being used by 2).  */

int predict(filt,x_dim,y_dim,x_pos,y_pos,result,f_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, f_or_e;
  {
  register int y_filt,x_filt, y_res,x_res;
  register double taps_used = 0.0; /* int *** */
  register double fraction = 0.0;
  int filt_sz = x_dim*y_dim;
  int x_start = ((x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0));
  int y_start = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
  int i;

  for (i=0; i<filt_sz; i++) result[i] = 0.0;

  for (y_filt=0, y_res=y_start;
       y_filt<filt_sz;
       y_filt+=x_dim, y_res+=x_dim)
    if ((y_res >= 0) AND (y_res < filt_sz))
      for (x_filt=y_filt, x_res=x_start;
	   x_filt<y_filt+x_dim;
	   x_filt++, x_res++)
	if ((x_res >= 0) AND (x_res < x_dim))
	  {
	    result[y_res+x_res] = filt[x_filt];
	    taps_used += ABS(filt[x_filt]);
	  }
  printf("TU: %f\n",taps_used);
  if (f_or_e IS FILTER)
      {
      /* fraction = ( (double) filt_sz ) / ( (double) taps_used ); */
      for (i=0; i<filt_sz; i++) fraction += ABS(filt[i]);
      fraction = ( fraction / taps_used );
      for (i=0; i<filt_sz; i++) result[i] *= fraction;
      }
  return(0);
  }


/* --------------------------------------------------------------------
Reflect, multiplying tap of filter which is at the edge of the image
by root 2.  This maintains orthogonality of odd-length linear-phase
QMF filters, but it is not useful for most applications, since it
alters the DC level.  */

int ereflect(filt,x_dim,y_dim,x_pos,y_pos,result,f_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, f_or_e;
  {
  register int y_filt,x_filt, y_edge,x_edge;
  register int x_base = (x_pos>0)?(x_dim-1):0;
  register int y_base = x_dim * ( (y_pos>0)?(y_dim-1):0 ); 
  int filt_sz = x_dim*y_dim;
  int x_edge_dist = (x_pos>1)?(x_pos-x_dim):((x_pos<-1)?(x_pos+1):0);
  int y_edge_dist = x_dim * ( (y_pos>1)?(y_pos-y_dim):((y_pos<-1)?(y_pos+1):0) );
  int i;
  double norm,onorm;

  for (i=0; i<filt_sz; i++) result[i] = 0.0;

  /* reflect at boundary */       
  for (y_filt=0, y_edge=y_edge_dist;
       y_filt<filt_sz;
       y_filt+=x_dim, y_edge+=x_dim)
    for (x_filt=y_filt, x_edge=x_edge_dist;
	 x_filt<y_filt+x_dim;
	 x_filt++, x_edge++)
      result[ABS(y_base-ABS(y_edge)) + ABS(x_base-ABS(x_edge))]
	+= filt[x_filt];

  /* now multiply edge by root 2 */
  if (x_pos ISNT 0) 
    for (y_filt=x_base; y_filt<filt_sz; y_filt+=x_dim)
      result[y_filt] *= ROOT2;
  if (y_pos ISNT 0) 
    for (x_filt=y_base; x_filt<y_base+x_dim; x_filt++)
      result[x_filt] *= ROOT2;

  /* now normalize to norm of original filter */
  for (norm=0.0,i=0; i<filt_sz; i++)
    norm += (result[i]*result[i]);
  norm=sqrt(norm);

  for (onorm=0.0,i=0; i<filt_sz; i++)
    onorm += (filt[i]*filt[i]);
  onorm = sqrt(onorm);

  norm = norm/onorm;
  for (i=0; i<filt_sz; i++)
    result[i] /= norm;
  return(0);
  }


/* ------- printout stuff for testing ------------------------------
  printf("Xpos: %d, Ypos: %d", x_pos, y_pos);
    for (y_filt=0; y_filt<y_dim; y_filt++)
      {
      printf("\n");
      for (x_filt=0; x_filt<x_dim; x_filt++)
	printf("%6.1f", result[y_filt*x_dim+x_filt]);
      }
  printf("\n");
*/



/* Local Variables: */
/* buffer-read-only: t */
/* End: */