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

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

int
MRISscaleUp(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 int
change_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 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) ;
}

static int
compute_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 int
assemble_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 int
find_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) ;
}

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

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

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

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