/*!
  \fn MRI *RFstat2P(MRI *rf, RFS *rfs, MRI *binmask, int TwoSided, MRI *p)
  \brief Converts a stat to a p value. If TwoSided, then computes a p value
      based on an unsigned stat, but the sign is still passed to p.
*/
MRI *RFstat2P(MRI *rf, RFS *rfs, MRI *binmask, int TwoSided, MRI *p)
{
  int c,r,s,f=0,m;
  double v,pval;

  if (RFname2Code(rfs) == -1) return(NULL);
  p = MRIclone(rf,p);

  for (c=0; c < rf->width; c++)
  {
    for (r=0; r < rf->height; r++)
    {
      for (s=0; s < rf->depth; s++)
      {
        if (binmask != NULL)
        {
          m = (int)MRIgetVoxVal(binmask,c,r,s,0);
          if (!m) continue;
        }
        for (f=0; f < rf->nframes; f++)
        {
          v = MRIgetVoxVal(rf,c,r,s,f);
	  if(TwoSided) pval = SIGN(v)*2*RFstat2PVal(rfs,fabs(v));
	  else         pval = RFstat2PVal(rfs,v);
          MRIsetVoxVal(p,c,r,s,f,pval);
        }
      }
    }
  }
  return(p);
}

static MRI *
make_atrophy_map(MRI *mri_time1, MRI *mri_time2, MRI *mri_dst, TRANSFORM *transform1, TRANSFORM *transform2,
                 int *gray_labels, int ngray, int *csf_labels, int ncsf) {
  int            x, y, z, label1, label2, n, found, xp, yp, zp, spacing ;
  GCA_MORPH_NODE *gcamn1, *gcamn2 ;
  GCA_MORPH      *gcam1, *gcam2 ;
  float           volume ;

  if (mri_dst == NULL) {
    mri_dst = MRIalloc(mri_time1->width, mri_time1->height, mri_time1->depth, MRI_FLOAT) ;
    MRIcopyHeader(mri_time1, mri_dst) ;
  }

  gcam1 = (GCA_MORPH*)transform1->xform ;
  gcam2 = (GCA_MORPH*)transform2->xform ;
  spacing = gcam1->spacing ;

  for (x = 0 ; x < mri_time1->width ; x++) {
    xp = x / spacing;
    for (y = 0 ; y < mri_time1->height ; y++) {
      yp = y / spacing;
      for (z = 0 ; z < mri_time1->depth ; z++) {
        if (x == Gx && y == Gy && z == Gz)
          DiagBreak() ;
        label1 = MRIgetVoxVal(mri_time1, x, y, z, 0) ;
        label2 = MRIgetVoxVal(mri_time2, x, y, z, 0) ;
        if (label1 == label2)
          continue ;

        /* if label1 was one of the gray types and label2 one of the csf, call it atrophy */
        for (found = n = 0 ; n < ngray ; n++)
          if (label1 == gray_labels[n]) {
            found = 1 ;
            break ;
          }
        if (found == 0)
          continue ;
        for (found = n = 0 ; n < ncsf ; n++)
          if (label2 == csf_labels[n]) {
            found = 1 ;
            break ;
          }
        if (found == 0)
          continue ;
        zp = z / spacing;
        gcamn1 = &gcam1->nodes[xp][yp][zp] ;
        gcamn2 = &gcam2->nodes[xp][yp][zp] ;
        volume = 0 ;
        if (FZERO(gcamn1->area) == 0)
          volume += gcamn1->orig_area / gcamn1->area ;
        if (FZERO(gcamn2->area) == 0)
          volume += gcamn2->orig_area / gcamn2->area ;
        MRIsetVoxVal(mri_dst, x, y, z, 0, volume) ;
      }
    }
  }


  return(mri_dst) ;
}

static MRI *
compute_bias(MRI *mri_src, MRI *mri_dst, MRI *mri_bias)
{
  int x, y, z ;
  float bias, src, dst ;

  if (!mri_bias)
    mri_bias = MRIalloc
               (mri_src->width, mri_src->height, mri_src->depth, MRI_FLOAT) ;

  MRIcopyHeader(mri_src, mri_bias) ;
  for (x = 0 ; x < mri_src->width ; x++)
  {
    for (y = 0; y < mri_src->height ; y++)
    {
      for (z = 0 ; z < mri_src->depth ; z++)
      {
        src = MRIgetVoxVal(mri_src, x, y, z, 0) ;
        dst = MRIgetVoxVal(mri_dst, x, y, z, 0) ;
        if (FZERO(src))
        {
          bias = 1 ;
        }
        else
        {
          bias = dst/src ;
        }
        MRIsetVoxVal(mri_bias, x, y, z, 0, bias) ;
      }
    }
  }

  return(mri_bias) ;
}

MRI *
MRIcombineDistanceTransforms(MRI *mri_src1, MRI *mri_src2, MRI *mri_dst)
{
  int   x, y, z, f ;
  float val1, val2 ;

  if (mri_dst == NULL)
  {
    mri_dst = MRIclone(mri_src1, NULL) ;
  }

  for (f = 0 ; f < mri_dst->nframes ; f++)
    for (x = 0 ; x < mri_dst->width ; x++)
      for (y = 0 ; y < mri_dst->height ; y++)
        for (z = 0 ; z < mri_dst->depth ; z++)
        {
          val1 = MRIgetVoxVal(mri_src1, x, y, z, f) ;
          val2 = MRIgetVoxVal(mri_src2, x, y, z, f) ;
          if (val2 < 0 && val1 > 0)
          {
            val1 = val2 ;  // in the interior of 1
          }
          else if (val2 > 0 && val1 > val2)  // exterior of both, but closer to border of 2
          {
            val1 = val2 ;
          }
          else if (val2 < 0 && val1 < val2)
          {
            val1 = val2 ;  // interior of both, but closer to border of 2
          }

          MRIsetVoxVal(mri_dst, x, y, z, f, val1) ;
        }
  return(mri_dst) ;
}

static int
normalize_PD(MRI *mri_PD, float target) {
  double mean_PD, scale, val ;
  int    x, y, z ;

  for (mean_PD = 0.0, x = 0 ; x < mri_PD->width ; x++) {
    for (y = 0 ; y < mri_PD->height ; y++) {
      for (z = 0 ; z < mri_PD->depth ; z++) {
        mean_PD += (double)MRIgetVoxVal(mri_PD, x, y, z,0) ;
      }
    }
  }
  mean_PD /= (mri_PD->width * mri_PD->height * mri_PD->depth) ;
  scale = target / mean_PD ;
  printf("mean PD %2.0f, scaling by %2.2f to set mean to %2.0f\n",
         mean_PD, scale, target) ;
  for (mean_PD = 0.0, x = 0 ; x < mri_PD->width ; x++) {
    for (y = 0 ; y < mri_PD->height ; y++) {
      for (z = 0 ; z < mri_PD->depth ; z++) {
        val = (double)MRIgetVoxVal(mri_PD, x, y, z,0) ;
        val *= scale ;
        MRIsetVoxVal(mri_PD, x, y, z,0, val);
      }
    }
  }
  return(NO_ERROR) ;
}

static int
normalize_timepoints_with_parzen_window(MRI *mri, double cross_time_sigma)
{
  int   frame1, frame2, x, y, z ;
  double val0, val, total, g, norm, total_norm ;

  norm = 1 / sqrt(2 * M_PI * SQR(cross_time_sigma)) ;
  for (x = 0 ; x < mri->width ; x++)
    for (y = 0 ; y < mri->height ; y++)
      for (z = 0 ; z < mri->depth ; z++)
      {
        if (x == Gx && y == Gy && z == Gz)
          DiagBreak() ;
        for (frame1 = 0 ; frame1 < mri->nframes ; frame1++)
        {
          val0 = MRIgetVoxVal(mri, x, y, z, frame1) ;
          for (total = total_norm = 0.0, frame2 = 0 ; frame2 < mri->nframes ; frame2++)
          {
            val = MRIgetVoxVal(mri, x, y, z, frame2) ;
            g = norm * exp( - SQR(val-val0) / (2 * SQR(cross_time_sigma))) ;
            total += g*val ; 
            total_norm += g ;
          }
          total /= total_norm ;
          MRIsetVoxVal(mri, x, y, z, frame1, total) ;
        }
      }


  return(NO_ERROR) ;
}

/*-------------------------------------------------------------------*/
MRI *RFrescale(MRI *rf, RFS *rfs, MRI *binmask, MRI *rfout)
{
  int c,r,s,f,m;
  double v, gmean, gstddev, gmax;

  if (RFname2Code(rfs) == -1) return(NULL);
  RFexpectedMeanStddev(rfs); // expected
  RFglobalStats(rf, binmask, &gmean, &gstddev, &gmax); //actual

  rfout = MRIclone(rf,rfout);

  for (c=0; c < rf->width; c++)
  {
    for (r=0; r < rf->height; r++)
    {
      for (s=0; s < rf->depth; s++)
      {
        if (binmask != NULL)
        {
          m = (int)MRIgetVoxVal(binmask,c,r,s,0);
          if (!m) continue;
        }
        for (f=0; f < rf->nframes; f++)
        {
          v = MRIgetVoxVal(rf,c,r,s,f);
          v = (v - gmean)*(rfs->stddev/gstddev) + rfs->mean;
          MRIsetVoxVal(rfout,c,r,s,f,v);
        }
      }
    }
  }
  return(rfout);
}

static int
normalize_timepoints_with_samples(MRI *mri, GCA_SAMPLE *gcas, int nsamples, int nsoap)
{
  int   frame, i, x, y, z ;
  double target, val ;
  MRI    *mri_ctrl, *mri_bias, *mri_target, *mri_frame ;

  mri_ctrl = MRIcloneDifferentType(mri, MRI_UCHAR) ;
  mri_bias = MRIcloneDifferentType(mri, MRI_FLOAT) ;
  mri_target = MRIcloneDifferentType(mri, MRI_FLOAT) ;

  for (i = 0 ; i < nsamples ; i++)
  {
    if (i == Gdiag_no)
      DiagBreak() ;
    x = nint(gcas[i].x)  ; y = nint(gcas[i].y) ; z = nint(gcas[i].z) ;
    MRIsetVoxVal(mri_ctrl, x, y, z, 0, CONTROL_MARKED) ;
    for (target = 0.0, frame = 0 ; frame < mri->nframes ; frame++)
      target += MRIgetVoxVal(mri, x, y, z, frame) ;
    target /= mri->nframes ;
    MRIsetVoxVal(mri_target, x, y, z, 0, target) ;
  }

  // build a bias correction for each time point (which each has its own frame)
  for (frame = 0 ; frame < mri->nframes ; frame++)
  {
    MRIclear(mri_bias) ; 
    for (i = 0 ; i < nsamples ; i++)
    {
      if (i == Gdiag_no)
        DiagBreak() ;
      x = nint(gcas[i].x)  ; y = nint(gcas[i].y) ; z = nint(gcas[i].z) ;
      target = MRIgetVoxVal(mri_target, x, y, z, 0) ;
      val = MRIgetVoxVal(mri, x, y, z, frame) ;
      if (FZERO(val))
        val = 1.0 ;
      MRIsetVoxVal(mri_bias, x, y, z, 0, target/val) ;
    }
    MRIbuildVoronoiDiagram(mri_bias, mri_ctrl, mri_bias) ;
    MRIsoapBubble(mri_bias, mri_ctrl, mri_bias, nsoap) ;
    mri_frame = MRIcopyFrame(mri, NULL, frame, 0) ;
    MRImultiply(mri_frame, mri_bias, mri_frame) ;
    if (Gdiag & DIAG_WRITE && DIAG_VERBOSE_ON)
    {
      char fname[STRLEN] ;
      sprintf(fname, "frame%d.mgz", frame) ;
      MRIwrite(mri_frame, fname) ;
      sprintf(fname, "bias%d.mgz", frame) ;
      MRIwrite(mri_bias, fname) ;
      sprintf(fname, "target%d.mgz", frame) ;
      MRIwrite(mri_target, fname) ;
    }
    MRIcopyFrame(mri_frame, mri, 0, frame) ;
  }
  MRIfree(&mri_bias) ; MRIfree(&mri_target) ; MRIfree(&mri_ctrl) ;
  return(NO_ERROR) ;
}

MATRIX *ComputeAdjMatrix(MRI *mri_label, MRI *mri_mask, int minlabel, int maxlabel)
{
  MATRIX *AdjMatrix;
  int i, j, label1, label2, offset, numLabels;
  int depth, width, height;
  int x,y,z,cx,cy,cz;

  numLabels = maxlabel - minlabel + 1;

  depth = mri_label->depth;
  width = mri_label->width;
  height = mri_label->height;

  AdjMatrix = (MATRIX *)MatrixAlloc(numLabels, numLabels, MATRIX_REAL);

  if (!AdjMatrix)
    ErrorExit(ERROR_BADPARM, "%s: unable to allowcate memory.\n", Progname);
  /* The diagnoal entries of AdjMatrix is set to zero and remain zero */
  for (i=1; i <= numLabels;i++)
    for (j=i; j <= numLabels; j++)
    {
      AdjMatrix->rptr[i][j] = 0.0;
      AdjMatrix->rptr[j][i] = 0.0;
    }


  for (z=0; z < depth; z++)
    for (y=0; y< height; y++)
      for (x=0; x < width; x++)
      {
        if (MRIvox(mri_mask, x, y, z) == 0) continue;
        label1 = (int) MRIgetVoxVal(mri_label, x, y, z,0);
        if (label1 < minlabel || label1 > maxlabel) continue;

        /* Find all 6-neighbor with different label */
        for (offset = 0; offset < 6; offset++)
        {
          cx = x + xoff[offset];
          cy = y + yoff[offset];
          cz = z + zoff[offset];
          if (cx < 0 || cx >= width || cy < 0 || cy >= height
              || cz < 0 || cz >= depth) continue;

          label2 = (int) MRIgetVoxVal(mri_label, cx, cy, cz,0);
          if (label2 < minlabel || label2 > maxlabel || label2 == label1)
            continue;

          AdjMatrix->rptr[label1-minlabel+1][label2-minlabel+1] = 1.0;
          AdjMatrix->rptr[label2-minlabel+1][label1-minlabel+1] = 1.0;
        } /* for_offset */
      }

  return (AdjMatrix);
}

/*
 figure out what to do with voxels that were turned 'off' by the
 topology correction. This is a hack, but for now just make them
 the most likely of the nbr voxel labels.
*/
static int
resegment_erased_voxels(MRI *mri_T1, MRI *mri_in, MRI *mri_out, int target_label) {
  int       x, y, z, label_in, label_out, xi, yi, zi, xk, yk, zk, label, changed=0 ;
  HISTOGRAM *histos[MAX_CMA_LABEL+1] ;
  double    p, max_p, val ;

  build_label_histograms(mri_in, mri_T1, histos) ;
  for (x = 0 ; x < mri_in->width ; x++) {
    for (y = 0 ; y < mri_in->height ; y++) {
      for (z = 0 ; z < mri_in->depth ; z++) {
        label_in = nint(MRIgetVoxVal(mri_in, x, y, z, 0)) ;
        label_out = nint(MRIgetVoxVal(mri_out, x, y, z, 0)) ;
        if (label_in == target_label) {
          // find most likely nbr label
          max_p = 0 ;
          label_out = label_in ;
          for (xk = -1 ; xk <= 1 ; xk++) {
            xi = x + xk ;
            if (xi < 0 || xi >= mri_in->width)
              continue ;
            for (yk = -1 ; yk <= 1 ; yk++) {
              yi = y + yk ;
              if (yi < 0 || yi >= mri_in->height)
                continue ;
              for (zk = -1 ; zk <= 1 ; zk++) {
                zi = z + zk ;
                if (zi < 0 || zi >= mri_in->depth)
                  continue ;
                label = nint(MRIgetVoxVal(mri_in, xi, yi, zi, 0)) ;
                if (label == label_in)
                  continue ;  // would be topologically incorrect
                val = MRIgetVoxVal(mri_T1, xi, yi, zi, 0) ;
                p = HISTOvalToCount(histos[label], val) ;
                if (p > max_p) {
                  max_p = p ;
                  label_out = label ;
                }
              }
            }
          }
          changed++ ;
          MRIsetVoxVal(mri_out, x, y, z, 0, label_out) ;
        }
      }
    }
  }
  printf("%d voxels resegmented to be ML\n", changed) ;
  return(NO_ERROR) ;
}

int
MRIaccumulateMeansAndVariances(MRI *mri, MRI *mri_mean, MRI *mri_std) {
  int    x, y, z, width, height, depth ;
  float  val, *pmean, *pstd ;

  width = mri->width ;
  height = mri->height ;
  depth = mri->depth ;

  for (z = 0 ; z < depth ; z++) {
    for (y = 0 ; y < height ; y++) {
      pmean = &MRIFvox(mri_mean, 0, y, z) ;
      pstd = &MRIFvox(mri_std, 0, y, z) ;
      for (x = 0 ; x < width ; x++) {
        val = MRIgetVoxVal(mri,x,y,z,0) ;
        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)
          DiagBreak() ;
#if 1
        *pmean++ += (float) val ;
        *pstd++ += ((float) val)*((float) val) ;
#else
        MRIFvox(mri_mean,x,y,z) += val ;
        MRIFvox(mri_std,x,y,z) += val*val ;
#endif
      }
    }
  }
  return(NO_ERROR) ;
}

/*-------------------------------------------------------------------*/
int RFsynth(MRI *rf, RFS *rfs, MRI *binmask)
{
  int c,r,s,f;
  double v,m;

  if (RFname2Code(rfs) == -1) return(1);

  for (c=0; c < rf->width; c++)
  {
    for (r=0; r < rf->height; r++)
    {
      for (s=0; s < rf->depth; s++)
      {
        if (binmask != NULL)
        {
          m = MRIgetVoxVal(binmask,c,r,s,0);
          if (m < 0.5) continue;
        }
        for (f=0; f < rf->nframes; f++)
        {
          v = RFdrawVal(rfs);
          MRIsetVoxVal(rf,c,r,s,f,v);
        }
      }
    }
  }
  return(0);
}

static int
apply_bias_field(MRI *mri, int nbias, float *bias_coefs[3][2]) {
  int    x, y, z, n ;
  double xb, yb, zb, x0, y0, z0, w0x, w0y, w0z ;
  float  val ;

  x0 = mri->width/2 ;
  y0 = mri->height/2 ;
  z0 = mri->depth/2 ;
  w0x = 2/x0 ;
  w0y = 2/y0 ;
  w0z = 2/z0 ;
  for (x = 0 ; x < mri->width ; x++) {
    for (xb = 1.0, n=1 ; n <= nbias ; n++)
      xb += bias_coefs[0][0][n-1] * cos(w0x*n*(x-x0)) + bias_coefs[0][1][n-1] * sin(w0x*n*(x-x0)) ;
    for (y = 0 ; y < mri->height ; y++) {
      for (yb = 1.0, n=1 ; n <= nbias ; n++)
        yb += bias_coefs[1][0][n-1] * cos(w0y*n*(y-y0)) + bias_coefs[1][1][n-1] * sin(w0y*n*(y-y0)) ;
      for (z = 0 ; z < mri->depth ; z++) {
        for (zb = 1.0, n=1 ; n <= nbias ; n++)
          zb += bias_coefs[2][0][n-1] * cos(w0z*n*(z-z0)) + bias_coefs[2][1][n-1] * sin(w0z*n*(z-z0)) ;
        val = MRIgetVoxVal(mri, x, y, z, 0) ;
        val = val * xb * yb * zb ;
        MRIsetVoxVal(mri, x, y, z, 0, val) ;
      }
    }
  }

  return(NO_ERROR) ;
}

static int
is_wmsa_border(MRI *mri_seg, int x, int y, int z)
{
  int found_wmsa, found_non_wmsa, label, xi, yi, zi, xk, yk, zk ;

  for (found_wmsa = found_non_wmsa = 0, xk = -1 ; xk <= 1 ; xk++)
  {
    xi = mri_seg->xi[x+xk] ;
    for (yk = -1 ; yk <= 1 ; yk++)
    {
      yi = mri_seg->yi[y+yk] ;
      for (zk = -1 ; zk <= 1 ; zk++)
      {
	if (abs(xk)+abs(yk)+abs(zk) > 1)  // only 6-connected
	  continue ;
	zi = mri_seg->zi[z+zk] ;
	label = (int)MRIgetVoxVal(mri_seg, xi, yi, zi, 0) ;
	if (IS_WMSA(label))
	  found_wmsa++ ;
	else
	  found_non_wmsa++ ;
	if (found_wmsa && found_non_wmsa)
	  return(1) ;
      }
    }
  }
  return(0) ;
}

int
relabel_hypointensities_neighboring_gray(MRI *mri)
{
  int    x, y, z, label, changed, i ;
  MRI    *mri_tmp = NULL ;

  for (changed = i = 0 ; i < 2 ; i++) {
    mri_tmp = MRIcopy(mri, mri_tmp) ;
    for (x = 0 ; x < mri->width ; x++) {
      for (y = 0 ; y < mri->height ; y++) {
        for (z = 0 ; z < mri->depth ; z++) {
          label = MRIgetVoxVal(mri_tmp, x, y, z, 0) ;
          if (label != WM_hypointensities) {
            continue ;
          }
          if (MRIneighbors(mri_tmp, x, y, z, Left_Cerebral_Cortex) > 0) {
            MRIsetVoxVal(mri, x, y, z, 0, Left_Cerebral_Cortex) ;
            changed++ ;
          } else if (MRIneighbors(mri_tmp,x,y,z,Right_Cerebral_Cortex) > 0) {
            MRIsetVoxVal(mri, x, y, z, 0, Right_Cerebral_Cortex) ;
            changed++ ;
          }
        }
      }
    }
  }

  printf("%d hypointense voxels neighboring cortex changed\n", changed) ;
  return(NO_ERROR) ;
}

static int
most_frequent_label(MRI *mri_seg, MRI_SEGMENT *mseg)
{
  int  label_counts[MAX_CMA_LABELS], i, max_count, max_label, label ;

  memset(label_counts, 0, sizeof(label_counts)) ;
  for (i = max_count = max_label = 0 ; i < mseg->nvoxels ; i++)
  {
    label = MRIgetVoxVal(mri_seg, 
                         mseg->voxels[i].x, 
                         mseg->voxels[i].y, 
                         mseg->voxels[i].z, 0) ;
    if (IS_WM(label) == 0 && IS_UNKNOWN(label) == 0)
    {
      label_counts[label]++ ;
    }
  }
  for (i = max_count = max_label = 0 ; i < MAX_CMA_LABELS ; i++)
  {
    if (label_counts[i] > max_count)
    {
      max_count = label_counts[i] ;
      max_label = i ;
    }
  }
  return(max_label) ;
}

static int
rip_vertices_out_of_fov(MRI_SURFACE *mris, MRI *mri_profiles) {
  int  vno, n, i, good, nsamples ;
  VERTEX *v ;

  nsamples = mri_profiles->nframes ;
  for (vno = 0 ; vno < mris->nvertices ; vno++) {
    v = &mris->vertices[vno] ;
    for (good = i = 0 ; good == 0 && i < nsamples ; i++)
      if (!FZERO(MRIgetVoxVal(mri_profiles, vno, 0, 0, i)))
        good = 1 ;
    if (!good) {
      v->ripflag = 1 ;
      v->annotation = 0 ;
    }
  }

  // rip stuff around the identically 0 ones, as they can't be trusted
  for (vno = 0 ; vno < mris->nvertices ; vno++) {
    v = &mris->vertices[vno] ;
    if (v->ripflag != 1)
      continue ;
    for (n = 0 ; n < v->vtotal ; n++) {
      mris->vertices[v->v[n]].ripflag = 2 ;
      mris->vertices[v->v[n]].annotation = 0 ;
    }
  }
  for (vno = 0 ; vno < mris->nvertices ; vno++) {
    v = &mris->vertices[vno] ;
    if (v->ripflag == 2)
      v->ripflag = 1 ;
  }
  return(NO_ERROR) ;
}

int
MRIaccumulateMaskedMeansAndVariances(MRI *mri, MRI *mri_mask, MRI *mri_dof,
                                     float low_val,
                                     float hi_val,MRI *mri_mean,MRI *mri_std) {
  int     x, y, z, width, height, depth ;
  float   val ;
  BUFTYPE *pmask, mask ;

  width = mri->width ;
  height = mri->height ;
  depth = mri->depth ;

  for (z = 0 ; z < depth ; z++) {
    for (y = 0 ; y < height ; y++) {
      pmask = &MRIvox(mri_mask, 0, y, z) ;
      for (x = 0 ; x < width ; x++) {
        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)
          DiagBreak() ;
        mask = *pmask++ ;
        if (mask >= low_val && mask <= hi_val) {
          val = MRIgetVoxVal(mri,x,y,z,0) ;
          MRIFvox(mri_mean,x,y,z) += val ;
          MRIFvox(mri_std,x,y,z) += val*val ;
          MRIvox(mri_dof,x,y,z)++ ;
        }
      }
    }
  }
  return(NO_ERROR) ;
}

/*-------------------------------------------------------------------------------
  MRIbinarize2() - same as MRIbinarize() but passes theshold, low, and hi as
  doubles instead of UCHARs.
  -------------------------------------------------------------------------------*/
MRI * MRIbinarize2(MRI *mri_src, MRI *mri_dst,
                   double threshold, double low_val, double hi_val) {
  int     width, height, depth, x, y, z, f ;
  double  val;

  if (!mri_dst)  mri_dst = MRIclone(mri_src, NULL) ;

  width = mri_src->width ;
  height = mri_src->height ;
  depth = mri_src->depth ;

  for (f = 0 ; f < mri_src->nframes ; f++) {
    for (z = 0 ; z < depth ; z++) {
      for (y = 0 ; y < height ; y++) {
        for (x = 0 ; x < width ; x++) {
          val = MRIgetVoxVal(mri_src, x, y, z, f);
          if (val > threshold) val = hi_val ;
          else                val = low_val ;
          MRIsetVoxVal(mri_dst, x, y, z, f, val) ;
        }
      }
    }
  }

  return(mri_dst) ;
}

static int
rip_bad_vertices(MRI_SURFACE *mris, MRI *mri_profiles) {
  int    unknown_index, bad, i, vno, index ;
  VERTEX *v ;

  CTABfindName(mris->ct, "Unknown", &unknown_index) ;
  printf("unknown index = %d\n", unknown_index) ;
  for (vno = 0 ; vno < mris->nvertices ; vno++) {
    v = &mris->vertices[vno] ;
    if (vno == Gdiag_no)
      DiagBreak() ;
    if (v->ripflag)
      continue ;

    CTABfindAnnotation(mris->ct, v->annotation, &index);
    if (index == unknown_index) {
      v->annotation = 0 ;
      v->ripflag = 1 ;
    } else {
      bad = 1 ;
      for (i = 0 ; i < mri_profiles->nframes ; i++) {
        if (!FZERO(MRIgetVoxVal(mri_profiles, vno, 0, 0, i)))
          bad = 0 ;
      }
      if (bad) {
        v->ripflag = 1 ;
        v->annotation = 0 ;
      }
    }

  }
  MRISripFaces(mris) ;
  MRIScomputeMetricProperties(mris) ;
  return(NO_ERROR) ;
}