mrprior.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
float med9(
  float *inp, 
  int dim
) {
//  xxx
//  xXx
//  xxx
  float w[9], *p;
  int n=0;
 
  /* upper row */
  p=inp-dim-1; 
  w[n]=*p++; if(isfinite(w[n])) n++;
  w[n]=*p++; if(isfinite(w[n])) n++; 
  w[n]=*p; if(isfinite(w[n])) n++;
  /* mid row */
  p=inp-1; 
  w[n]=*p++; if(isfinite(w[n])) n++; 
  w[n]=*p++; if(isfinite(w[n])) n++; 
  w[n]=*p; if(isfinite(w[n])) n++;
  /* lower row */
  p=inp+dim-1; 
  w[n]=*p++; if(isfinite(w[n])) n++; 
  w[n]=*p++; if(isfinite(w[n])) n++; 
  w[n]=*p; if(isfinite(w[n])) n++;
  /* Median */
  return(fmedian(w, n));
}
/*****************************************************************************/
float med21(
  float *inp, 
  int dim
) {
//  oxxxo
//  xxxxx
//  xxXxx
//  xxxxx
//  oxxxo
  float w[21], *p;
 
  /* 1st row */
  p=inp-2*dim-1; w[0]=*p++; w[1]=*p++; w[2]=*p;
  /* 2nd row */
  p=inp-dim-2; w[3]=*p++; w[4]=*p++; w[5]=*p++; w[6]=*p++; w[7]=*p;
  /* 3rd row */
  p=inp-2; w[8]=*p++; w[9]=*p++; w[10]=*p++; w[11]=*p++; w[12]=*p;
  /* 4th row */
  p=inp+dim-2; w[13]=*p++; w[14]=*p++; w[15]=*p++; w[16]=*p++; w[17]=*p;
  /* 5th row */
  p=inp+2*dim-1; w[18]=*p++; w[19]=*p++; w[20]=*p;
  /* Median */
  return(fmedian(w, 21));
}
/*****************************************************************************/
void do_prior(
  float *img, 
  float beta, 
  float *med_coef, 
  int dim, 
  float small,
  int maskdim,
  float *maxm
) {
  int i, j;
  float *iptr, *mptr, med, f, mmax=0.0;
  float one_minus_beta=1.0-beta;
 
  //printf("do_prior(%d, %d)\n", dim, maskdim); fflush(stdout);
  for(i=0; i<dim*dim; i++) med_coef[i]=0.0;
 
  if(maskdim==3) {
    iptr=img+dim+1;
    mptr=med_coef+dim+1;
    for(j=1; j<dim-1; j++) {
      for(i=1; i<dim-1; i++) {
        if(*iptr<=small) {
          *mptr=0.0;
        } else {
          med=med9(iptr, dim);
          if(med==0.0) {
            *mptr=0.0;
          } else {
            f=med/(beta*(*iptr) + one_minus_beta*med);
            if(f>1.0E-08) *mptr=f; else *mptr=0.0;
          }
        }
        if(*mptr>mmax) mmax=*mptr;
        iptr++; mptr++;
      }
      iptr+=2; mptr+=2;
    }
    if(maxm!=NULL) *maxm=mmax;
    //printf("out from do_prior()\n"); fflush(stdout);
    return;
  }
 
  if(maskdim==5) {
    iptr=img+2*dim+2;
    mptr=med_coef+2*dim+2;
    for(j=2; j<dim-2; j++) {
      for(i=2; i<dim-2; i++) {
        if(*iptr<=small) {
          *mptr=0.0;
        } else {
          med=med21(iptr, dim);
          if(med==0.0) {
            *mptr=0.0;
          } else {
            f=med/(beta*(*iptr) + one_minus_beta*med);
            if(f>1.0E-08) *mptr=f; else *mptr=0.0;
          }
        }
        if(*mptr>mmax) mmax=*mptr;
        iptr++; mptr++;
      }
      iptr+=4; mptr+=4;
    }
    if(maxm!=NULL) *maxm=mmax;
    //printf("out from do_prior()\n"); fflush(stdout);
    return;
  }
}
/*****************************************************************************/
 
/*****************************************************************************/