自学教程：C++ DiagBreak函数代码示例

static MRI *MRIcomputePriors(MRI *mri_priors, int ndof, MRI *mri_char_priors) {  int     width, height, depth, x, y, z ;  BUFTYPE *pchar_prior, char_prior ;  float   *pprior, prior ;  width = mri_priors->width ;  height = mri_priors->height ;  depth = mri_priors->depth ;  if (!mri_char_priors) {    mri_char_priors = MRIalloc(width, height, depth, MRI_UCHAR) ;  }  for (z = 0 ; z < depth ; z++) {    for (y = 0 ; y < height ; y++) {      pprior = &MRIFvox(mri_priors, 0, y, z) ;      pchar_prior = &MRIvox(mri_char_priors, 0, y, z) ;      for (x = 0 ; x < width ; x++) {        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)          DiagBreak() ;        prior = *pprior++ ;        if (prior > 0)          DiagBreak() ;        if (prior > 10)          DiagBreak() ;        char_prior = (BUFTYPE)nint(100.0f*prior/(float)ndof) ;        if (char_prior > 101)          DiagBreak() ;        *pchar_prior++ = char_prior ;      }    }  }  return(mri_char_priors) ;}

intMRIcomputeMeansAndStds(MRI *mri_mean, MRI *mri_std, int ndof) {  int    x, y, z, width, height, depth ;  float  sum, sum_sq, mean, var ;  width = mri_std->width ;  height = mri_std->height ;  depth = mri_std->depth ;  for (z = 0 ; z < depth ; z++) {    for (y = 0 ; y < height ; y++) {      for (x = 0 ; x < width ; x++) {        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)          DiagBreak() ;        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;        sum = MRIFvox(mri_mean,x,y,z) ;        sum_sq = MRIFvox(mri_std,x,y,z) / ndof ;        mean = MRIFvox(mri_mean,x,y,z) = sum / ndof ;        var = sum_sq - mean*mean;        if (fabs(var) > 1e10 || fabs(mean)  > 1e10)          DiagBreak() ;        MRIFvox(mri_std,x,y,z) = sqrt(var) ;      }    }  }  return(NO_ERROR) ;}

intMRISscaleUp(MRI_SURFACE *mris){  int     vno, n, max_v, max_n ;  VERTEX  *v ;  float   ratio, max_ratio ;  max_ratio = 0.0f ;  max_v = max_n = 0 ;  for (vno = 0 ; vno < mris->nvertices ; vno++)  {    v = &mris->vertices[vno] ;    if (v->ripflag)    {      continue ;    }    if (vno == Gdiag_no)    {      DiagBreak() ;    }    for (n = 0 ; n < v->vnum ; n++)    {      if (FZERO(v->dist[n]))   /* would require infinite scaling */      {        continue ;      }      ratio = v->dist_orig[n] / v->dist[n] ;      if (ratio > max_ratio)      {        max_v = vno ;        max_n = n ;        max_ratio = ratio ;      }    }  }  fprintf(stderr, "max @ (%d, %d), scaling brain by %2.3f/n",          max_v, max_n, max_ratio) ;#if 0  MRISscaleBrain(mris, mris, max_ratio) ;#else  for (vno = 0 ; vno < mris->nvertices ; vno++)  {    v = &mris->vertices[vno] ;    if (v->ripflag)    {      continue ;    }    if (vno == Gdiag_no)    {      DiagBreak() ;    }    for (n = 0 ; n < v->vnum ; n++)    {      v->dist_orig[n] /= max_ratio ;    }  }#endif  return(NO_ERROR) ;}

static intchange_label(MRI *mri_T1,MRI *mri_labeled,int x,int y,int z,int wsize,int left) {  float  wm_mean, hippo_mean, val ;  int    label ;  if (x == 95 && y == 127 && z == 119) /* dark wm (68) */    DiagBreak() ;  if (x == 94 && y == 126 && z == 119)  /* bright hippo (104) */    DiagBreak() ;  val = (float)MRIvox(mri_T1, x, y, z) ;  wm_mean = label_mean(mri_T1, mri_labeled, x, y, z, wsize,                       left ? Left_Cerebral_White_Matter :                       Right_Cerebral_White_Matter) ;  hippo_mean = label_mean(mri_T1, mri_labeled, x, y, z, wsize,                          left ? Left_Hippocampus :                          Right_Hippocampus) ;  if (fabs(wm_mean-val) < fabs(hippo_mean-val))    label = left ?            Left_Cerebral_White_Matter : Right_Cerebral_White_Matter;  else    label = left ?            Left_Hippocampus : Right_Hippocampus;  MRIvox(mri_labeled, x, y, z) = label ;  return(NO_ERROR) ;}

static intnormalize_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) ;}

static intcompute_rigid_gradient(MRI_SURFACE *mris, MRI *mri, double *pdx, double *pdy, double *pdz) {  int    vno ;  VERTEX *v ;  Real   val, xw, yw, zw, dx, dy, dz, delV, x, y, z, Ix, Iy, Iz, xv, yv, zv ;  dx = dy = dz = 0.0 ;  for (vno = 0 ; vno < mris->nvertices ; vno++) {    v = &mris->vertices[vno] ;    if (v->ripflag)      continue ;    if (vno == Gdiag_no)      DiagBreak() ;    x = v->x ;    y = v->y ;    z = v->z ;    MRISvertexToVoxel(mris, v, mri, &xw, &yw, &zw);    MRIsampleVolume(mri, xw, yw, zw, &val) ;    MRIsampleVolumeGradient(mri, xw, yw, zw,  &Ix, &Iy, &Iz) ;    // convert back to surface coords    xw += Ix ;    yw += Iy ;    zw += Iz ;    if (mris->useRealRAS)      MRIvoxelToWorld(mri, xw, yw, zw, &xv, &yv, &zv) ;    else      MRIvoxelToSurfaceRAS(mri, xw, yw, zw, &xv, &yv, &zv) ;    Ix = xv-v->x ;    Iy = yv-v->y;    Iz = zv-v->z ;    delV = v->val - val ;    dx += delV * Ix ;    dy += delV * Iy ;    dz += delV * Iz ;    if (!finitep((float)dx))      DiagBreak() ;    if (!finitep((float)dy))      DiagBreak() ;    if (!finitep((float)dz))      DiagBreak() ;  }  dx /= mris->nvertices ;  dy /= mris->nvertices ;  dz /= mris->nvertices ;  *pdx = dx ;  *pdy = dy ;  *pdz = dz ;  return(NO_ERROR) ;}

static intassemble_training_data_and_free_mris(MRI_SURFACE *mris[MAX_SUBJECTS], MRI *mri_overlays[MAX_SUBJECTS][MAX_OVERLAYS], 				     int nsubjects, int noverlays, int **ptraining_classes, double ***ptraining_data, int *pntraining) {  int    sno, i, ntraining, *training_classes, x, y, z, f, vno, tno ;  double **training_data ;  for (ntraining = sno = 0 ; sno < nsubjects ; sno++)    ntraining += MRISvalidVertices(mris[sno]) ;  printf("%2.1fM total training voxels found/n", (float)ntraining/(1024*1024.0)) ;  training_classes = (int *)calloc(ntraining, sizeof(int)) ;  if (training_classes == NULL)    ErrorExit(ERROR_NOFILE, "%s: could not allocate %d-len training class vector", Progname, ntraining) ;  training_data = (double **)calloc(ntraining, sizeof(double *)) ;  if (training_data == NULL)    ErrorExit(ERROR_NOFILE, "%s: could not allocate %d-len training data vector", Progname, ntraining) ;  for (tno = sno = 0 ; sno < nsubjects ; sno++)  {    for (vno = x = 0 ; x < mri_overlays[sno][0]->width ; x++)      for (y = 0 ; y < mri_overlays[sno][0]->height ; y++)	for (z = 0 ; z < mri_overlays[sno][0]->depth ; z++)	  for (f = 0 ; f < mri_overlays[sno][0]->nframes ; f++, vno++)	  {	    if (tno == Gdiag_no)	      DiagBreak() ;	    if (vno == Gdiag_no)	      DiagBreak();	    if (mris[sno]->vertices[vno].ripflag)  	      continue ;   // not in cortex 	    training_classes[tno] = mris[sno]->vertices[vno].marked ;	    training_data[tno] = (double *)calloc(noverlays, sizeof(double)) ;	    if (training_data[tno] == NULL)	      ErrorExit(ERROR_NOFILE, "%s: could not allocate %d-len training data vector [%d]", Progname, noverlays, tno) ;	    for (i = 0 ; i < noverlays ; i++)	      training_data[tno][i] = MRIgetVoxVal(mri_overlays[sno][i], x, y, z, f) ;	    tno++ ;	  }        for (i = 0 ; i < noverlays ; i++)      MRIfree(&mri_overlays[sno][i]) ;     MRISfree(&mris[sno]) ;  }  *pntraining = ntraining ;  *ptraining_classes = training_classes ;  *ptraining_data = training_data ;  return(NO_ERROR) ;}

static intfind_optimal_feature_and_threshold(RANDOM_FOREST *rf, TREE *tree, NODE *parent, NODE *left, NODE *right, int *training_classes, double **training_data, int ntraining){  double  info_gain, best_info_gain, thresh, best_thresh, fmin, fmax, tdata ;  int     f, best_f, fno, nsteps, i, tno;  info_gain = best_info_gain = -1 ; best_f = -1 ; best_thresh = 0 ;  for (f = 0 ; f < tree->nfeatures ; f++)  {    fno = tree->feature_list[f] ;    if (fno == Gdiag_no)      DiagBreak() ;    fmin = rf->feature_max[fno] ; fmax = rf->feature_min[fno] ;    for (tno = 0 ; tno < parent->total_counts ; tno++)    {      i = parent->training_set[tno] ;      tdata = training_data[i][fno] ;      if (tdata < fmin)	fmin = tdata ;      if (tdata > fmax)	fmax = tdata ;    }    nsteps = rf->nsteps ;    do    {      thresh = find_optimal_threshold(rf, tree, parent, left, right, rf->training_data, rf->ntraining,				      fno, &info_gain, nsteps, fmin, fmax);      if (info_gain < 0)	DiagBreak() ;      nsteps *= 5 ;      if (nsteps > rf->max_steps)	break ;  // don't keep trying forever    } while (info_gain <= 0) ;    if (info_gain > best_info_gain)    {      best_info_gain = info_gain ;      best_f = fno ;      best_thresh = thresh ;    }  }  if (best_f < 0)    return(0) ;  info_gain = compute_info_gain(rf, tree, parent, left, right, training_data, best_f,best_thresh) ;  adjust_optimal_threshold(rf, tree, parent, left, right, training_data, best_f,&best_thresh);  parent->thresh = best_thresh ;  parent->feature = best_f ;  return(1) ;}

static MRI *MRIsynthesizeWithFAF(MRI *mri_T1, MRI *mri_PD, MRI *mri_dst, double TR, double alpha, double TE, int nfaf,                     float *faf_coefs[3][2]) {  int   x, y, z, width, height, depth, n ;  Real flash, T1, PD ;  double xb, yb, zb, x0, y0, z0, w0x, w0y, w0z ;  x0 = mri_T1->width/2 ;  y0 = mri_T1->height/2 ;  z0 = mri_PD->depth/2 ;  w0x = 2/x0 ;  w0y = 2/y0 ;  w0z = 2/z0 ;  if (!mri_dst)    mri_dst = MRIclone(mri_T1, NULL) ;  mri_dst->tr = TR ;  mri_dst->flip_angle = alpha ;  mri_dst->te = TE ;  mri_dst->ti = 0 ;  width = mri_T1->width ;  height = mri_T1->height ;  depth = mri_T1->depth ;  for (x = 0 ; x < width ; x++) {    for (xb = 1.0, n=1 ; n <= nfaf ; n++)      xb += faf_coefs[0][0][n-1] * cos(w0x*n*(x-x0)) + faf_coefs[0][1][n-1] * sin(w0x*n*(x-x0)) ;    for (y = 0 ; y < height ; y++) {      for (yb = 1.0, n=1 ; n <= nfaf ; n++)        yb += faf_coefs[1][0][n-1] * cos(w0y*n*(y-y0)) + faf_coefs[1][1][n-1] * sin(w0y*n*(y-y0)) ;      for (z = 0 ; z < depth ; z++) {        for (zb = 1.0, n=1 ; n <= nfaf ; n++)          zb += faf_coefs[2][0][n-1] * cos(w0z*n*(z-z0)) + faf_coefs[2][1][n-1] * sin(w0z*n*(z-z0)) ;        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;        MRIsampleVolume(mri_T1, x, y, z, &T1) ;        if (T1 <= 0)          T1 = 1 ;        if (T1 < 900 && T1 > 600)          DiagBreak() ;        MRIsampleVolume(mri_PD, x, y, z, &PD) ;        flash = FLASHforwardModel(T1, PD, TR, xb*yb*zb*alpha, TE) ;        MRIsetVoxVal(mri_dst, x, y, z, 0, flash) ;      }    }  }  return(mri_dst) ;}

static MRI *MRIupdatePriors(MRI *mri_binary, MRI *mri_priors) {  int     width, height, depth, x, y, z ;  BUFTYPE *pbin ;  float   prob ;  width = mri_binary->width ;  height = mri_binary->height ;  depth = mri_binary->depth ;  if (!mri_priors) {    mri_priors = MRIalloc(width, height, depth, MRI_FLOAT) ;  }  for (z = 0 ; z < depth ; z++) {    for (y = 0 ; y < height ; y++) {      pbin = &MRIvox(mri_binary, 0, y, z) ;      for (x = 0 ; x < width ; x++) {        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)          DiagBreak() ;        prob = *pbin++ / 100.0f ;        MRIFvox(mri_priors, x, y, z) += prob ;      }    }  }  return(mri_priors) ;}

intMRIaccumulateMaskedMeansAndVariances(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) ;}

intMRIcomputeMaskedMeansAndStds(MRI *mri_mean, MRI *mri_std, MRI *mri_dof) {  int    x, y, z, width, height, depth, ndof ;  float  sum, sum_sq, mean, var ;  width = mri_std->width ;  height = mri_std->height ;  depth = mri_std->depth ;  for (z = 0 ; z < depth ; z++) {    for (y = 0 ; y < height ; y++) {      for (x = 0 ; x < width ; x++) {        if (x == DEBUG_X && y == DEBUG_Y && z == DEBUG_Z)          DiagBreak() ;        sum = MRIFvox(mri_mean,x,y,z) ;        ndof = MRIvox(mri_dof, x, y, z) ;        if (!ndof)   /* variance will be 0 in any case */          ndof = 1 ;        sum_sq = MRIFvox(mri_std,x,y,z) / ndof ;        mean = MRIFvox(mri_mean,x,y,z) = sum / ndof ;        var = sum_sq - mean*mean;        MRIFvox(mri_std,x,y,z) = sqrt(var) ;      }    }  }  return(NO_ERROR) ;}

intMRIaccumulateMeansAndVariances(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) ;}

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) ;}

MRI *MRIScomputeDistanceMap(MRI_SURFACE *mris, MRI *mri_distance, int ref_vertex_no){  int    vno ;  VERTEX *v ;  double circumference, angle, distance ;  VECTOR *v1, *v2 ;  if (mri_distance == NULL)    mri_distance = MRIalloc(mris->nvertices, 1, 1, MRI_FLOAT) ;  v1 = VectorAlloc(3, MATRIX_REAL) ; v2 = VectorAlloc(3, MATRIX_REAL) ;  v = &mris->vertices[ref_vertex_no] ;  VECTOR_LOAD(v1, v->x, v->y, v->z) ;  /* radius vector */  circumference = M_PI * 2.0 * V3_LEN(v1) ;  for (vno = 0 ; vno < mris->nvertices ; vno++)  {    v = &mris->vertices[vno] ;    if (vno == Gdiag_no)      DiagBreak() ;    VECTOR_LOAD(v2, v->x, v->y, v->z) ;  /* radius vector */    angle = fabs(Vector3Angle(v1, v2)) ;    distance = circumference * angle / (2.0 * M_PI) ;    MRIsetVoxVal(mri_distance, vno, 0, 0, 0, distance) ;  }  VectorFree(&v1) ; VectorFree(&v2) ;  return(mri_distance) ;}

static intnormalize_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) ;}

/*-----------------------------------------------------  Parameters:  Returns value:  Description  ------------------------------------------------------*/static intmriSegmentReallocateSegments(MRI_SEGMENTATION *mriseg, int max_segments){  MRI_SEGMENT *old_segments ;  int         s ;  old_segments = mriseg->segments ;  mriseg->segments = (MRI_SEGMENT *)calloc(max_segments, sizeof(MRI_SEGMENT)) ;  if (!mriseg->segments)    ErrorExit(ERROR_NOMEMORY, "MRIsegmentAlloc: could not alloc mrisegments");  mriseg->max_segments = max_segments ;  for (s = 0 ; s < mriseg->nsegments ; s++)  {    if (s == 1)      DiagBreak() ;    mriseg->segments[s].area = old_segments[s].area ;    mriseg->segments[s].nvoxels = old_segments[s].nvoxels ;    mriseg->segments[s].max_voxels = old_segments[s].max_voxels ;    mriseg->segments[s].voxels = old_segments[s].voxels ;  }  free(old_segments) ;  for (s = mriseg->nsegments ; s < mriseg->max_segments ; s++)  {    mriseg->segments[s].voxels = (MSV *)calloc(MAX_VOXELS, sizeof(MSV)) ;    if (!mriseg->segments[s].voxels)      ErrorExit(ERROR_NOMEMORY, "mriSegmentRealloc: could not alloc voxels");    mriseg->segments[s].max_voxels = MAX_VOXELS ;  }  return(NO_ERROR) ;}

static intrip_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) ;}

static intcheck_volume(MRI *mri_save, MRI *mri_out, int target_label) {  int   x, y, z, sval, oval ;  for (x = 0 ; x < mri_save->width ; x++)    for (y = 0 ; y < mri_save->height ; y++)      for (z = 0 ; z < mri_save->depth ; z++) {        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;        sval = nint(MRIgetVoxVal(mri_save, x, y, z, 0)) ;        oval = nint(MRIgetVoxVal(mri_out, x, y, z, 0)) ;        if ((oval == target_label && sval == 0) ||            (oval != target_label && sval > 0))          DiagBreak() ;      }  return(NO_ERROR) ;}

static MRI *edit_amygdala(MRI *mri_in_labeled, MRI *mri_T1, MRI *mri_out_labeled) {  int   width, height, depth, x, y, z, nchanged, label, total_changed,  dup, ddown, left ;  MRI   *mri_tmp ;  mri_out_labeled = MRIcopy(mri_in_labeled, mri_out_labeled) ;  mri_tmp = MRIcopy(mri_out_labeled, NULL) ;  width = mri_T1->width ;  height = mri_T1->height ;  depth = mri_T1->depth ;  total_changed = 0 ;  do {    nchanged = 0 ;    /* change gray to wm if near amygdala */    for (z = 0 ; z < depth ; z++) {      for (y = 0 ; y < height ; y++) {        for (x = 0 ; x < width ; x++) {          if (x == 95 && y == 127 && z == 119)            DiagBreak() ;          label = MRIvox(mri_tmp, x, y, z) ;          left = 0 ;          switch (label) {          case Left_Cerebral_Cortex:            left = 1 ;          case Right_Cerebral_Cortex:            dup = distance_to_label(mri_out_labeled,                                    left ?  Left_Amygdala :                                    Right_Amygdala,x,y,z,0,-1,0,2);            ddown = distance_to_label(mri_out_labeled,                                      left ? Left_Cerebral_White_Matter :                                      Right_Cerebral_White_Matter,                                      x,y,z,0,1,0,3);            if (dup <= 1 && ddown <= 1) {              label = left ?                      Left_Cerebral_White_Matter : Right_Cerebral_White_Matter;              MRIvox(mri_tmp, x, y, z) = label ;              nchanged++ ;              continue ;            }          default:            break ;          }        }      }    }    MRIcopy(mri_tmp, mri_out_labeled) ;    total_changed += nchanged ;  } while (nchanged > 0) ;  MRIfree(&mri_tmp) ;  printf("%d amygdala voxels changed./n", total_changed) ;  return(mri_out_labeled) ;}

static MRI *edit_cortical_gray_matter(MRI *mri_in_labeled,MRI *mri_T1,MRI *mri_out_labeled) {  int   width, height, depth, x, y, z, nchanged, label, total_changed,  left, niter, change ;  MRI   *mri_tmp ;  mri_out_labeled = MRIcopy(mri_in_labeled, mri_out_labeled) ;  mri_tmp = MRIcopy(mri_out_labeled, NULL) ;  width = mri_T1->width ;  height = mri_T1->height ;  depth = mri_T1->depth ;  niter = total_changed = 0 ;  do {    nchanged = 0 ;    /* change gray to wm if near amygdala */    for (z = 0 ; z < depth ; z++) {      for (y = 0 ; y < height ; y++) {        for (x = 0 ; x < width ; x++) {          if (x == 95 && y == 127 && z == 119)            DiagBreak() ;          label = MRIvox(mri_tmp, x, y, z) ;          change = left = 0 ;          switch (label) {          case Left_Cerebral_Cortex:          case Unknown:          case Right_Cerebral_Cortex:            if (neighborLabel(mri_tmp,x,y,z,1,Left_Cerebral_White_Matter) &&                neighborLabel(mri_tmp,x,y,z,1,Right_Cerebral_White_Matter)) {              left =                sagittalNeighbors(mri_tmp,x,y,z,3,Left_Cerebral_White_Matter)>                sagittalNeighbors(mri_tmp,x,y,z,3,Right_Cerebral_White_Matter);              label = left ?                      Left_Cerebral_White_Matter : Right_Cerebral_White_Matter;              MRIvox(mri_tmp, x, y, z) = label ;              nchanged++ ;              continue ;            }          default:            break ;          }        }      }    }    MRIcopy(mri_tmp, mri_out_labeled) ;    total_changed += nchanged ;  } while ((nchanged > 0) && (niter++ < 3)) ;  MRIfree(&mri_tmp) ;  printf("%d midline voxels changed./n", total_changed) ;  return(mri_out_labeled) ;}

static intmark_clusters(MRI_SURFACE *mris, MRI *mri_profiles, CLUSTER *ct, int k) {  int    i, vno, min_i, nsamples, num[MAX_CLUSTERS], nchanged ;  double min_dist, dist ;  VECTOR *v1 ;  VERTEX *v ;  memset(num, 0, sizeof(num)) ;  nsamples = mri_profiles->nframes ;  v1 = VectorAlloc(nsamples, MATRIX_REAL) ;  for (nchanged = vno = 0 ; vno < mris->nvertices ; vno++) {    v = &mris->vertices[vno] ;    if (v->ripflag)      continue ;    if (vno == Gdiag_no)      DiagBreak() ;    for (i = 0 ; i < nsamples ; i++)      VECTOR_ELT(v1, i+1) = MRIgetVoxVal(mri_profiles, vno, 0, 0, i) ;    min_dist = MatrixMahalanobisDistance(ct[0].v_mean, ct[0].m_cov, v1) ;    min_i = 0 ;    for (i = 1 ; i < k ; i++) {      if (i == Gdiag_no)        DiagBreak() ;      dist = MatrixMahalanobisDistance(ct[i].v_mean, ct[i].m_cov, v1) ;      if (dist < min_dist) {        min_dist = dist ;        min_i = i ;      }    }    CTABannotationAtIndex(mris->ct, min_i, &v->annotation) ;    if (v->curv != min_i)      nchanged++ ;    v->curv = min_i ;    num[min_i]++ ;  }  for (i = 0 ; i < k ; i++)    if (num[i] == 0)      DiagBreak() ;  VectorFree(&v1) ;  printf("%d vertices changed clusters.../n", nchanged) ;  return(nchanged) ;}

static MRI *MRIsynthesize(MRI *mri_T1, MRI *mri_PD, MRI *mri_T2star, MRI *mri_dst, double TR, double alpha, double TE) {  int   x, y, z, width, height, depth ;  Real flash, T1, PD ;  if (!mri_dst)    mri_dst = MRIclone(mri_T1, NULL) ;  mri_dst->tr = TR ;  mri_dst->flip_angle = alpha ;  mri_dst->te = TE ;  mri_dst->ti = 0 ;  width = mri_T1->width ;  height = mri_T1->height ;  depth = mri_T1->depth ;  for (x = 0 ; x < width ; x++) {    for (y = 0 ; y < height ; y++) {      for (z = 0 ; z < depth ; z++) {        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;        MRIsampleVolume(mri_T1, x, y, z, &T1) ;        if (T1 <= 0)          T1 = 1 ;        if (T1 < 900 && T1 > 600)          DiagBreak() ;        MRIsampleVolume(mri_PD, x, y, z, &PD) ;        if (mri_T2star) {          Real T2star ;          MRIsampleVolume(mri_T2star, x, y, z, &T2star) ;          flash = FLASHforwardModelT2star(T1, PD, T2star, TR, alpha, TE) ;        } else          flash = FLASHforwardModel(T1, PD, TR, alpha, TE) ;        MRIsetVoxVal(mri_dst, x, y, z, 0, flash) ;        if (!finite(flash))          DiagBreak() ;      }    }  }  return(mri_dst) ;}

intcheck_finite(char *where, double what){  if (!finite(what))  {    ErrorPrintf(ERROR_BADPARM, "%s not finite!/n",where) ;    DiagBreak() ;    return(0) ;  }  return(1) ;}

static intmle_label(MRI *mri_T1, MRI *mri_labeled, int x, int y, int z, int wsize, int l1, int l2) {  float  l1_mean, l2_mean, val ;  int    label ;  if (x == 95 && y == 127 && z == 119) /* dark wm (68) */    DiagBreak() ;  if (x == 94 && y == 126 && z == 119)  /* bright hippo (104) */    DiagBreak() ;  val = (float)MRIvox(mri_T1, x, y, z) ;  l1_mean = label_mean(mri_T1, mri_labeled, x, y, z, wsize, l1) ;  l2_mean = label_mean(mri_T1, mri_labeled, x, y, z, wsize, l2) ;  if (fabs(l1_mean-val) < fabs(l2_mean-val))    label = l1 ;  else    label = l2 ;  MRIvox(mri_labeled, x, y, z) = label ;  return(NO_ERROR) ;}

MRI *MRIcropVolumeToLabel(MRI *mri_src,                      MRI *mri_dst,                      LABEL *area,                      MRI_SURFACE *mris_white,                      MRI_SURFACE *mris_pial){  MHT    *mht_white, *mht_pial ;  int    x, y, z ;  VERTEX *v_white, *v_pial ;  Real   xs, ys, zs ;  mht_white = MHTfillVertexTableRes(mris_white, NULL, CURRENT_VERTICES, 5.0) ;  mht_pial = MHTfillVertexTableRes(mris_pial, NULL, CURRENT_VERTICES, 5.0) ;  if (mri_dst == NULL)    mri_dst = MRIclone(mri_src, NULL) ;  MRISclearMarks(mris_white) ; MRISclearMarks(mris_pial) ;  LabelMarkSurface(area, mris_white) ; LabelMarkSurface(area, mris_pial) ;  for (x = 0 ; x < mri_src->width ; x++)  {    for (y = 0 ; y < mri_src->height ; y++)    {      for (z = 0 ; z < mri_src->depth ; z++)      {        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;                      if (MRIgetVoxVal(mri_src, x, y, z, 0) > 0)        {          MRISsurfaceRASFromVoxel(mris_white, mri_dst,                                   x, y, z,                                   &xs, &ys, &zs) ;          v_white = MHTfindClosestVertexInTable(mht_white, mris_white,                                                 xs, ys, zs, 1) ;          v_pial = MHTfindClosestVertexInTable(mht_pial, mris_pial,                                                xs, ys, zs, 1) ;          if ((v_white && v_white->marked == 1) ||              (v_pial && v_pial->marked == 1))          {            MRIsetVoxVal(mri_dst, x, y, z, 0, 255) ;          }          else            MRIsetVoxVal(mri_dst, x, y, z, 0, 0) ;        }      }    }  }  MHTfree(&mht_white) ; MHTfree(&mht_pial) ;  return(mri_dst) ;}

static doublecvector_average_in_label(float *v, LABEL *area, int num) {  int    i ;  double avg ;  for (avg = 0.0, i = 0 ; i < area->n_points ; i++) {    if (!finite(v[area->lv[i].vno]))      DiagBreak() ;    avg += v[area->lv[i].vno] ;  }  avg /= (double)area->n_points ;  return(avg) ;}

static intadjust_optimal_threshold(RF *rf, TREE *tree, NODE *parent, NODE *left, NODE *right, double **training_data, int fno, double *pbest_thresh) {  double previous_thresh, next_thresh, info_gain, best_info_gain, step, thresh, best_thresh ;  best_thresh = *pbest_thresh ;  best_info_gain = compute_info_gain(rf, tree, parent, left, right, training_data, fno,best_thresh) ;  if (rf->min_step_size > 0)    step = rf->min_step_size ;  else    step = (rf->feature_max[fno] - rf->feature_min[fno]) / (10*rf->nsteps-1);  for (thresh = best_thresh ; thresh <= rf->feature_max[fno] ; thresh += step)  {    info_gain = compute_info_gain(rf, tree, parent, left, right, training_data, fno,thresh) ;    if (info_gain < 0)      DiagBreak() ;    if (info_gain > best_info_gain)    {      best_thresh = thresh ;      best_info_gain = info_gain ;    }    if (info_gain < best_info_gain)      break ;  }  next_thresh = thresh - step ;  for (thresh = best_thresh ; thresh >= rf->feature_min[fno] ; thresh -= step)  {    info_gain = compute_info_gain(rf, tree, parent, left, right, training_data, fno,thresh) ;    if (info_gain > best_info_gain)    {      best_thresh = thresh ;      best_info_gain = info_gain ;    }    if (info_gain < best_info_gain)      break ;  }  previous_thresh = thresh+step ;  thresh = (next_thresh + previous_thresh) / 2 ;  // maximize margin  info_gain = compute_info_gain(rf, tree, parent, left, right, training_data, fno,thresh) ;  if (info_gain < best_info_gain) // don't use it    compute_info_gain(rf, tree, parent, left, right, training_data, fno,best_thresh) ;  else    best_thresh = thresh ;  *pbest_thresh = best_thresh ;    return(NO_ERROR) ;}

static MRI *apply_bias(MRI *mri_orig, MRI *mri_norm, MRI *mri_bias) {  MATRIX   *m_vox2vox;  VECTOR   *v1, *v2;  int      x, y, z ;  double   xd, yd, zd, bias, val_orig, val_norm ;  if (mri_norm == NULL)    mri_norm = MRIclone(mri_orig, NULL) ;  m_vox2vox = MRIgetVoxelToVoxelXform(mri_orig, mri_bias) ;  v1 = VectorAlloc(4, MATRIX_REAL);  v2 = VectorAlloc(4, MATRIX_REAL);  VECTOR_ELT(v1, 4) = 1.0 ;  VECTOR_ELT(v2, 4) = 1.0 ;  for (x = 0 ; x < mri_orig->width ; x++) {    V3_X(v1) = x ;    for (y = 0 ; y < mri_orig->height ; y++) {      V3_Y(v1) = y ;      for (z = 0 ; z < mri_orig->depth ; z++) {        V3_Z(v1) = z ;        if (x == Gx && y == Gy && z == Gz)          DiagBreak() ;        val_orig = MRIgetVoxVal(mri_orig, x, y, z, 0) ;        MatrixMultiply(m_vox2vox, v1, v2) ;        xd = V3_X(v2) ;        yd = V3_Y(v2) ;        zd = V3_Z(v2);        MRIsampleVolume(mri_bias, xd, yd, zd, &bias) ;        val_norm = val_orig * bias ;        if (mri_norm->type == MRI_UCHAR) {          if (val_norm > 255)            val_norm = 255 ;          else if (val_norm < 0)            val_norm = 0 ;        }        MRIsetVoxVal(mri_norm, x, y, z, 0, val_norm) ;      }    }  }  MatrixFree(&m_vox2vox) ;  VectorFree(&v1) ;  VectorFree(&v2) ;  return(mri_norm) ;}

static intcopy_ctrl_points_to_volume(GCA_SAMPLE *gcas, int nsamples, MRI *mri_ctrl, int frame){  int   i, xv, yv, zv, nregions ;  nregions = mri_ctrl->nframes/2 ;  for (i = 0 ; i < nsamples ; i++)  {    xv = gcas[i].x ; yv = gcas[i].y ; zv = gcas[i].z ;    if (xv == Ggca_x && yv == Ggca_y && zv == Ggca_z)      DiagBreak() ;    MRIsetVoxVal(mri_ctrl, xv, yv, zv, frame, gcas[i].label) ;    MRIsetVoxVal(mri_ctrl, xv, yv, zv, frame+nregions, gcas[i].means[0]) ;  }  return(NO_ERROR) ;}