C++ MRIsetVoxVal函数代码示例

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

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

/*!
  \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);
}

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

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

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

int
MRIeraseSmallSegments(MRI_SEGMENTATION *mriseg, MRI *mri_seg, int min_voxels)
{
  int         s, v, nerased = 0, x, y, z ;
  MRI_SEGMENT *mseg ;

  if (DIAG_VERBOSE_ON && 0)
    fprintf(stdout, "compacting segments...\n") ;

  for (s = 0 ; s < mriseg->max_segments ; s++)
  {
    mseg = &mriseg->segments[s] ;
    if (mseg->nvoxels < min_voxels)
    {
      if (mseg->voxels)
	for (v = 0 ; v < mseg->nvoxels ; v++)
	{
	  x = mseg->voxels[v].x ; y = mseg->voxels[v].y ; z = mseg->voxels[v].z ;
	  if (x == Gx && y == Gy && z == Gz)
	    printf("MRIeraseSmallSegments: erasing (%d, %d %d)\n", x, y, z) ;
	  MRIsetVoxVal(mri_seg, x, y, z,0, 0) ;
	  nerased++ ;
	}
    }
  }
  return(nerased) ;
}

/*!
  \fn int CCSegment(MRI *seg, int segid, int segidunknown)
  Constraints a sementation ID to consist of voxels that
  are spatially contiguous (6 face neighbors, not edge
  or corner). The voxels in the largest cluster are not
  changed. The voxels in the other clusters are set to
  segidunknown.
*/
int CCSegment(MRI *seg, int segid, int segidunknown)
{
    MRI_SEGMENTATION *sgmnt;
    int k,kmax,index,c,r,s;

    sgmnt = MRIsegment(seg,segid-.5,segid+.5);
    printf("  Found %d clusters\n",sgmnt->nsegments);

    kmax = 0;
    for (k=0; k < sgmnt->nsegments; k++)
        if (sgmnt->segments[k].nvoxels > sgmnt->segments[kmax].nvoxels)
        {
            kmax = k;
        }

    for (k=0; k < sgmnt->nsegments; k++)
    {
        printf("     %d k %f\n",k,sgmnt->segments[k].area);
        if (k==kmax)
        {
            continue;
        }
        for (index = 0; index < sgmnt->segments[k].nvoxels; index++)
        {
            c = sgmnt->segments[k].voxels[index].x;
            r = sgmnt->segments[k].voxels[index].y;
            s = sgmnt->segments[k].voxels[index].z;
            MRIsetVoxVal(seg,c,r,s,0,segidunknown);
        }
    }
    MRIsegmentFree(&sgmnt);
    return(0);
}

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