float MRIvoxelsInLabelWithPartialVolumeEffects(MRI *mri, MRI *mri_vals, 
                                               int label, 
                                               MRI *mri_mixing_coef, 
                                               MRI *mri_nbr_labels)
{
  float   volume, vox_vol ;
  int     x, y, z, nbr_label_counts[20000];
  int     label_counts[20000], this_label, border;
  int     nbr_label, max_count, vox_label ;
  MRI     *mri_border ;
  float   label_means[20000], pv, mean_label, mean_nbr, val ;
  int     maxlabels = 20000;
  // DNG 6/7/07 : had to use maxlabels instead of MAX_CMA_LABELS here
  // so that segmentations with values > MAX_CMA_LABELS can be
  // accessed. This includes the cortical segmentations as well as
  // white matter segs. Currently, the max seg no is 4181, but this
  // could easily change.
  // NJS 2/17/10 : Indeed, it did change... the Destrieux a2009s atlas has
  // label values up to about 15000.

  if(label >= maxlabels){
    printf("ERROR: MRIvoxelsInLabelWithPartialVolumeEffects()\n");
    printf(" label %d exceeds maximum label number %d\n",label,maxlabels);
    return(-100000);
  }

  vox_vol = mri->xsize*mri->ysize*mri->zsize ;

  /* first find border voxels */
  mri_border = MRImarkLabelBorderVoxels(mri, NULL, label, 1, 1) ;
  if (DIAG_VERBOSE_ON && (Gdiag & DIAG_WRITE))
    MRIwrite(mri_border, "b.mgz") ;
  volume = 0 ;
  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() ;
        vox_label = MRIgetVoxVal(mri, x, y, z, 0) ;
        border = MRIgetVoxVal(mri_border, x, y, z, 0) ;
        if ((vox_label != label) && (border == 0))
          continue ;

        if (border == 0)  volume += vox_vol ;
        else{ /* compute partial volume */
          MRIcomputeLabelNbhd(mri, mri_vals, x, y, z,
           nbr_label_counts, label_means, 1, maxlabels) ;
          MRIcomputeLabelNbhd(mri, mri_vals, x, y, z,
           label_counts, label_means, 7, maxlabels) ;
          val = MRIgetVoxVal(mri_vals, x, y, z, 0) ;  /* compute partial
                                                     volume based on
                                                     intensity */
          mean_label = label_means[vox_label] ;
          nbr_label = -1 ;
          max_count = 0 ;
          /* look for a label that is a nbr and is
             on the other side of val from the label mean */
          for (this_label = 0 ; this_label < maxlabels;  this_label++)  {
            if(this_label == vox_label)  
              continue ;
            if(nbr_label_counts[this_label] == 0) 
              continue ; /* not a nbr */

            if ((label_counts[this_label] > max_count) &&
                ((label_means[this_label] - val) *
                 (mean_label - val) < 0))            {
              max_count = label_means[this_label] ;
              nbr_label = this_label ;
            }
          }
          if (vox_label != label && nbr_label != label)  
            continue ; // this struct not in voxel 

          if (max_count == 0) 
          {
            volume += vox_vol ; // couldn't find an appropriate label
            if (mri_nbr_labels)  // find max nbr label anyway for caller
            {
              for (this_label = 0 ; this_label < maxlabels;  this_label++)  {
                if(this_label == vox_label)  
                  continue ;
                if(nbr_label_counts[this_label] == 0) 
                  continue ; /* not a nbr */
                
                if (label_counts[this_label] > max_count)
                {
                  max_count = label_means[this_label] ;
                  nbr_label = this_label ;
                }
              }
              MRIsetVoxVal(mri_nbr_labels, x, y, z, 0, nbr_label) ;
              if (mri_mixing_coef)
                MRIsetVoxVal(mri_mixing_coef, x, y, z, 0, 1.0) ;
            }
          }
          else{                 // compute partial volume pct 
            mean_nbr = label_means[nbr_label] ;
            pv = (val - mean_nbr) / (mean_label - mean_nbr) ;
            if (pv < 0 || pv > 1)
              DiagBreak() ;
            if (pv > 1)
              pv = 1 ;
            if (pv < 0)
              continue ;  // shouldn't happen 
            if (vox_label != label)
              pv = 1-pv ;
            volume += vox_vol * pv ;
            if (mri_mixing_coef)
              MRIsetVoxVal(mri_mixing_coef, x, y, z, 0, pv);
            if (mri_nbr_labels)  // return nbr label to caller
            {
              if (vox_label != label)
                MRIsetVoxVal(mri_nbr_labels, x, y, z, 0, vox_label);
              else
                MRIsetVoxVal(mri_nbr_labels, x, y, z, 0, nbr_label);
            }
          }
        }
      }
    }
  }

  MRIfree(&mri_border) ;
  return(volume) ;
}