edges.c

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

	     }
      }

  return(0);
  }


/* --------------------------------------------------------------------
qreflect2() - Modified version of reflect2 that  works properly for 
even-length QMF filters.
*/

int qreflect2(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
  double *filt, *result;
  int x_dim, y_dim, x_pos, y_pos, r_or_e;
  {
  reflect2(filt,x_dim,y_dim,x_pos,y_pos,result,0);
  return(0);
  }

/* --------------------------------------------------------------------
repeat() - repeat edge pixel.  Continuous, with discontinuous first
derivative.
*/

int repeat(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, r_or_e;
  {
  register int y_filt,x_filt, y_res,x_res, y_tmp, x_tmp;
  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 x_overhang = ((x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0));
  int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
  int filt_sz = x_dim*y_dim;
  int mx_pos = (x_dim/2);
  int my_pos = x_dim * (y_dim/2);
  int i;

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

  if (r_or_e IS REDUCE)
    for (y_filt=0, y_res=y_overhang;
	 y_filt<filt_sz;
	 y_filt+=x_dim, y_res+=x_dim)
      for (x_filt=y_filt, x_res=x_overhang;
	   x_filt<y_filt+x_dim;
	   x_filt++, x_res++)
	/* Clip the index: */
	result[((y_res>=0)?((y_res<filt_sz)?y_res:(filt_sz-x_dim)):0) 
	       + ((x_res>=0)?((x_res<x_dim)?x_res:(x_dim-1)):0)]      
		 += filt[x_filt];
  else {
    if ((y_base-y_overhang) ISNT my_pos)
      for (y_res=y_base, y_filt=y_base-ABS(y_overhang);
	   y_filt > y_base-filt_sz;
	   y_filt-=x_dim, y_res-=x_dim)
	if ((x_base-x_overhang) ISNT mx_pos)
	  for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
	       x_filt > x_base-x_dim;
	       x_res--, x_filt--)
	    result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
	else  /* ((x_base-x_overhang) IS mx_pos) */
	  for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
	       x_filt > x_base-x_dim;
	       x_filt--, x_res--)
	    for(x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
	      result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_tmp)];
    else /* ((y_base-y_overhang) IS my_pos) */
      for (y_res=y_base, y_filt=y_base-ABS(y_overhang);
	   y_filt > y_base-filt_sz;
	   y_filt-=x_dim, y_res-=x_dim)
	for (y_tmp=y_filt; y_tmp > y_base-filt_sz; y_tmp-=x_dim)
	  if  ((x_base-x_overhang) ISNT mx_pos)
	    for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
		 x_filt > x_base-x_dim;
		 x_filt--, x_res--)
	      result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_tmp)+ABS(x_filt)];
	  else /* ((x_base-x_overhang) IS mx_pos) */
	    for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
		 x_filt > x_base-x_dim;
		 x_filt--, x_res--)
	      for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
		result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_tmp)+ABS(x_tmp)];
    } /* End, if r_or_e IS EXPAND */

  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,r_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, r_or_e;
  {
  int filt_sz = x_dim*y_dim;
  register int y_filt,x_filt, y_res,x_res, y_tmp, x_tmp;
  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 x_overhang = (x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0);
  int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
  int mx_pos = (x_pos<0)?(x_dim/2):((x_dim-1)/2);
  int my_pos = x_dim * ((y_pos<0)?(y_dim/2):((y_dim-1)/2));
  int i;

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

  /* Modeled on reflect1() */
  if (r_or_e IS REDUCE)
    for (y_filt=0, y_res=y_overhang;
	 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_overhang;
	     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];
	  else
	      {
	      result[y_res+ABS(x_base-ABS(x_res-x_base))] -= filt[x_filt];
	      result[y_res+x_base] += 2*filt[x_filt];
	      }
      else
	for (x_filt=y_filt, x_res=x_overhang;
	     x_filt<y_filt+x_dim;
	     x_filt++, x_res++)
	  if ((x_res>=0) AND (x_res<x_dim))
	      {
	      result[ABS(y_base-ABS(y_res-y_base))+x_res] -= filt[x_filt];
	      result[y_base+x_res] += 2*filt[x_filt];
	      }
	  else
	      {
	      result[ABS(y_base-ABS(y_res-y_base))+ABS(x_base-ABS(x_res-x_base))] 
		-= filt[x_filt];
	      result[y_base+x_base] += 2*filt[x_filt];
	      }
  else {  /* r_or_e ISNT  REDUCE */
      y_overhang = ABS(y_overhang); 
      x_overhang = ABS(x_overhang);
      for (y_res=y_base, y_filt = y_base-y_overhang;
	   y_filt > y_base-filt_sz;
	   y_filt-=x_dim, y_res-=x_dim)
	  {
	  for (x_res=x_base, x_filt=x_base-x_overhang;
	       x_filt > x_base-x_dim;
	       x_res--, x_filt--)
	    result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
	  if (x_pos ISNT 0)
	    if (x_overhang ISNT mx_pos)
	      for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
		   x_filt > x_base-x_dim;
		   x_res--, x_filt--)
		result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
	    else /* x_overhang IS mx_pos */
	      for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
		   x_filt > x_base-x_dim;
		   x_res--, x_filt--)
		for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
		  result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_filt)+ABS(x_tmp)];
	  }
       if (y_pos ISNT 0)
	 if (y_overhang ISNT my_pos)
	   for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang;
		y_filt > y_base-filt_sz;
		y_filt-=x_dim, y_res-=x_dim)
	       {
	       for (x_res=x_base, x_filt=x_base-x_overhang;
		    x_filt > x_base-x_dim;
		    x_res--, x_filt--)
		 result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
	       if ((x_pos ISNT 0) AND (x_overhang ISNT mx_pos))
		 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
		      x_filt > x_base-x_dim;
		      x_res--, x_filt--)
		   result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
	       }
	 else /* y_overhang IS my_pos */
	   for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang-x_dim;
		y_filt > y_base-filt_sz;
		y_res-=x_dim, y_filt-=x_dim)
	     for (y_tmp=y_filt; y_tmp > y_base-filt_sz; y_tmp-=x_dim)
		 {
		 for (x_res=x_base, x_filt=x_base-x_overhang;
		      x_filt > x_base-x_dim;
		      x_res--, x_filt--)
		   result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_tmp)+ABS(x_filt)];
		 if ((x_pos ISNT 0) AND (x_overhang IS mx_pos))
		   for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
			x_filt > x_base-x_dim;
			x_res--, x_filt--)
		     for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
		       result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_tmp)+ABS(x_tmp)];
		 }
      }  /* r_or_e ISNT  REDUCE */

  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,r_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, r_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]);
	  }

  if (r_or_e IS REDUCE)
      {
      /* 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,r_or_e)
  register double *filt, *result;
  register int x_dim;
  int y_dim, x_pos, y_pos, r_or_e;
  {
  register int y_filt,x_filt, y_res,x_res;
  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_overhang = (x_pos>1)?(x_pos-x_dim):((x_pos<-1)?(x_pos+1):0);
  int y_overhang = 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_res=y_overhang;
       y_filt<filt_sz;
       y_filt+=x_dim, y_res+=x_dim)
    for (x_filt=y_filt, x_res=x_overhang;
	 x_filt<y_filt+x_dim;
	 x_filt++, x_res++)
      result[ABS(y_base-ABS(y_res)) + ABS(x_base-ABS(x_res))]
	+= 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: */