Attachment 'gca.c'

Download

   1 /**
   2  * @file  gca.c
   3  * @brief Core routines implementing the Gaussian Classifier Atlas mechanism
   4  *
   5  * Routines implementing the mechanism of encoding voxel label information
   6  * based on probabilistic information estimated from a training set.
   7  *
   8  * Reference:
   9  * "Whole Brain Segmentation: Automated Labeling of Neuroanatomical
  10  * Structures in the Human Brain", Fischl et al.
  11  * (2002) Neuron, 33:341-355.
  12  */
  13 /*
  14  * Original Author: Bruce Fischl
  15  * CVS Revision Info:
  16  *    $Author: fischl $
  17  *    $Date: 2007/05/11 20:20:05 $
  18  *    $Revision: 1.230 $
  19  *
  20  * Copyright (C) 2002-2007,
  21  * The General Hospital Corporation (Boston, MA). 
  22  * All rights reserved.
  23  *
  24  * Distribution, usage and copying of this software is covered under the
  25  * terms found in the License Agreement file named 'COPYING' found in the
  26  * FreeSurfer source code root directory, and duplicated here:
  27  * https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferOpenSourceLicense
  28  *
  29  * General inquiries: freesurfer@nmr.mgh.harvard.edu
  30  * Bug reports: analysis-bugs@nmr.mgh.harvard.edu
  31  *
  32  */
  33 
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36 #include <math.h>
  37 #include <errno.h>
  38 
  39 #include "mri.h"
  40 #include "error.h"
  41 #include "diag.h"
  42 #include "utils.h"
  43 #include "gca.h"
  44 #include "proto.h"
  45 #include "fio.h"
  46 #include "transform.h"
  47 #include "macros.h"
  48 #include "utils.h"
  49 #include "cma.h"
  50 #include "flash.h"
  51 #include "talairachex.h"
  52 #include "mrimorph.h"
  53 #include "intensity_eig.h"
  54 #include "numerics.h"
  55 #include "mrisegment.h"
  56 
  57 #if WITH_DMALLOC
  58 #include <dmalloc.h>
  59 #endif
  60 
  61 int Ggca_label = -1 ;
  62 int Ggca_nbr_label = -1 ;
  63 int Ggca_x = -1 ;
  64 int Ggca_y = -1 ;
  65 int Ggca_z = -1 ;
  66 int Gxp = -1;
  67 int Gyp = -1;
  68 int Gzp = -1;
  69 int Gxn = -1; // 32;
  70 int Gyn = -1; // 21;
  71 int Gzn = -1; // 32;
  72 char *G_write_probs = NULL ;
  73 
  74 /* this is the hack section */
  75 double PRIOR_FACTOR = 0.1 ;
  76 #define LABEL_UNDETERMINED   255
  77 
  78 static int total_pruned = 0 ;
  79 
  80 #define MIN_DET 1e-7
  81 /* less than this, and the covariance matrix is poorly conditioned */
  82 
  83 #define MAX_LABELS_PER_GCAN         25
  84 #define MAX_DIFFERENT_LABELS        500
  85 #define MIN_VAR                     (5*5)   /* should make this configurable */
  86 #define BIG_AND_NEGATIVE            -10000000.0
  87 #define VERY_UNLIKELY               1e-10
  88 #define UNKNOWN_DIST                4  /* within 4 mm of some known label */
  89 #define GCA_OLD_VERSION             2.0
  90 #define GCA_UCHAR_VERSION           4.0  // labels were uchars in file
  91 #define GCA_INT_VERSION             5.0  // labels are ints in file
  92 #define DEFAULT_MAX_LABELS_PER_GCAN 4
  93 
  94 //static int gcapBrainIsPossible(GCA_PRIOR *gcap) ;
  95 #if 0
  96 static HISTOGRAM *gcaComputeHistogramNormalization(GCA *gca,
  97     HISTOGRAM *h_mri,
  98     int label) ;
  99 static float gcaFindCerebellarScaleFactor(GCA *gca,
 100     HISTOGRAM *h_mri,
 101     int label,
 102     FILE *logfp);
 103 static double GCAcomputeScaledMeanEntropy(GCA *gca, 
 104                                           MRI *mri, 
 105                                           TRANSFORM *transform,
 106     int *labels, float *scales, int nlabels) ;
 107 #endif
 108 static int gcaScale(GCA *gca, int *labels, int *contra_labels, 
 109                     float *scales, int nlabels, int dir);
 110 #if 0
 111 static int gcaMaxPriorLabel(GCA *gca, 
 112                             MRI *mri, 
 113                             TRANSFORM *transform, 
 114                             int x, int y, int z) ;
 115 #endif
 116 static HISTOGRAM *gcaGetLabelHistogram(GCA *gca, int label, int frame) ;
 117 int GCAmaxLabel(GCA *gca) ;
 118 static int gcaRelabelSegment(GCA *gca,
 119                              TRANSFORM *transform,
 120                              MRI *mri_inputs,
 121                              MRI *mri_dst,
 122                              MRI_SEGMENT *mseg) ;
 123 
 124 #define INTERP_PRIOR 0
 125 #if INTERP_PRIOR
 126 static float gcaComputePrior(GCA *gca, MRI *mri, TRANSFORM *transform,
 127                              int x, int y, int z, int label) ;
 128 #endif
 129 
 130 static int gcapGetMaxPriorLabel(GCA_PRIOR *gcap, double *p_prior) ;
 131 double compute_partial_volume_log_posterior(GCA *gca,
 132     GCA_NODE *gcan,
 133     GCA_PRIOR *gcap,
 134     float *vals, int l1, int l2) ;
 135 static double gcaComputeSampleConditionalLogDensity(GCA_SAMPLE *gcas,
 136     float *vals,
 137     int ninputs,
 138     int label) ;
 139 static VECTOR *load_sample_mean_vector(GCA_SAMPLE *gcas,
 140                                        VECTOR *v_means,
 141                                        int ninputs) ;
 142 static MATRIX *load_sample_covariance_matrix(GCA_SAMPLE *gcas,
 143     MATRIX *m_cov,
 144     int ninputs) ;
 145 static double sample_covariance_determinant(GCA_SAMPLE *gcas, int ninputs) ;
 146 static GC1D *findClosestValidGC(GCA *gca,
 147                                 int x, int y, int z,
 148                                 int label, int check_var) ;
 149 static GC1D *findGCInWindow(GCA *gca,
 150                             int x, int y, int z,
 151                             int label, int wsize) ;
 152 
 153 static int    MRIorderIndices(MRI *mri, short *x_indices, short *y_indices,
 154                               short *z_indices) ;
 155 static int gcaCheck(GCA *gca) ;
 156 static double gcaVoxelLogLikelihood(GCA *gca,
 157                                     MRI *mri_labels,
 158                                     MRI *mri_inputs,
 159                                     int x, int y, int z,
 160                                     TRANSFORM *transform);
 161 static double gcaVoxelGibbsLogLikelihood(GCA *gca,
 162     MRI *mri_labels,
 163     MRI *mri_inputs,
 164     int x, int y, int z,
 165     TRANSFORM *transform,
 166     double gibbs_coef) ;
 167 static double gcaNbhdGibbsLogLikelihood(GCA *gca,
 168                                         MRI *mri_labels,
 169                                         MRI *mri_inputs,
 170                                         int x, int y, int z,
 171                                         TRANSFORM *transform,
 172                                         double gibbs_coef) ;
 173 static double gcaGibbsImpossibleConfiguration(GCA *gca,
 174     MRI *mri_labels,
 175     int x, int y, int z,
 176     TRANSFORM *transform) ;
 177 static GCA_SAMPLE *gcaExtractLabelAsSamples(GCA *gca,
 178     MRI *mri_labeled,
 179     TRANSFORM *transform,
 180     int *pnsamples, int label) ;
 181 #if 0
 182 static GCA_SAMPLE *gcaExtractRegionLabelAsSamples(GCA *gca, MRI *mri_labeled,
 183     TRANSFORM *transform,
 184     int *pnsamples, int label,
 185     int xn, int yn, int zn,
 186     int wsize) ;
 187 #endif
 188 static GCA_SAMPLE *gcaExtractThresholdedRegionLabelAsSamples
 189 (GCA *gca,
 190  MRI *mri_labeled,
 191  TRANSFORM *transform,
 192  int *pnsamples,
 193  int label,
 194  int xn,
 195  int yn,
 196  int zn,
 197  int wsize,
 198  float pthresh);
 199 /*static double gcaComputeSampleConditionalDensity(GCA_SAMPLE *gcas,
 200   float *vals, int ninputs, int label) ;*/
 201 static double gcaComputeLogDensity(GC1D *gc, float *vals, \
 202                                    int ninputs, float prior, int label) ;
 203 static double gcaComputeSampleLogDensity(GCA_SAMPLE *gcas,
 204     float *vals, int ninputs) ;
 205 int MRIcomputeVoxelPermutation(MRI *mri, short *x_indices, short *y_indices,
 206                                short *z_indices);
 207 GCA *gcaAllocMax(int ninputs, float prior_spacing, float node_spacing, \
 208                  int width, int height, int depth,
 209                  int max_labels, int flags) ;
 210 static int GCAupdateNode(GCA *gca, MRI *mri, int xn, int yn, int zn,
 211                          float *vals,int label, GCA *gca_prune, int noint);
 212 static int GCAupdateNodeCovariance(GCA *gca, MRI *mri, int xn, int yn, int zn,
 213                                    float *vals,int label);
 214 static int GCAupdatePrior(GCA *gca,
 215                           MRI *mri,
 216                           int xn, int yn, int zn,
 217                           int label);
 218 static int GCAupdateNodeGibbsPriors(GCA *gca,MRI*mri,int xn,int yn,int zn,
 219                                     int x, int y, int z, int label);
 220 #if 0
 221 static int different_nbr_labels(GCA *gca, int x, int y, int z, int wsize,
 222                                 int label, float pthresh) ;
 223 #endif
 224 static int different_nbr_max_labels(GCA *gca, int x, int y, int z, int wsize,
 225                                     int label) ;
 226 static int gcaRegionStats(GCA *gca, int x0, int y0, int z0, int wsize,
 227                           float *priors, float *vars, \
 228                           float means[MAX_DIFFERENT_LABELS][MAX_GCA_INPUTS]) ;
 229 static int gcaFindBestSample(GCA *gca, int x, int y, int z, int best_label,
 230                              int wsize, GCA_SAMPLE *gcas);
 231 #if 0
 232 static int gcaFindClosestMeanSample(GCA *gca, float *means, float min_prior,
 233                                     int x, int y, int z,
 234                                     int label, int wsize, GCA_SAMPLE *gcas);
 235 static GC1D *gcaFindHighestPriorGC(GCA *gca, int x, int y, int z,int label,
 236                                    int wsize) ;
 237 #endif
 238 static double gcaGibbsImageLogLikelihood(GCA *gca,
 239     MRI *mri_labels,
 240     MRI *mri_inputs,
 241     TRANSFORM *transform) ;
 242 
 243 static int mriFillRegion(MRI *mri, int x,int y,int z,int fill_val,int whalf);
 244 static int gcaFindMaxPriors(GCA *gca, float *max_priors) ;
 245 static int gcaFindIntensityBounds(GCA *gca, float *pmin, float *pmax) ;
 246 static int dump_gcan(GCA *gca,
 247                      GCA_NODE *gcan,
 248                      FILE *fp,
 249                      int verbose,
 250                      GCA_PRIOR *gcap) ;
 251 static GCA_NODE *findSourceGCAN(GCA *gca,
 252                                 MRI *mri_src,
 253                                 TRANSFORM *transform,
 254                                 int x,int y,int z);
 255 #if 0
 256 static int getLabelMean(GCA_NODE *gcan, int label, \
 257                         float *pvar, float *means, int ninputs) ;
 258 #endif
 259 static int   borderVoxel(MRI *mri, int x, int y, int z) ;
 260 static int   GCAmaxLikelihoodBorderLabel(GCA *gca, MRI *mri_inputs,
 261     MRI *mri_labels, TRANSFORM *transform,
 262     int x, int y, int z, float min_ratio);
 263 
 264 
 265 float getPrior(GCA_PRIOR *gcap, int label) ;
 266 GCA_PRIOR *getGCAP(GCA *gca,
 267                    MRI *mri,
 268                    TRANSFORM *transform,
 269                    int xv, int yv, int zv) ;
 270 GCA_PRIOR *getGCAPfloat(GCA *gca,
 271                         MRI *mri,
 272                         TRANSFORM *transform,
 273                         float xv, float yv, float zv) ;
 274 GCA_NODE *getGCAN(GCA *gca,
 275                   MRI *mri,
 276                   TRANSFORM *transform,
 277                   int xv, int yv, int zv) ;
 278 static int gcaNodeToPrior(GCA *gca,
 279                           int xn, int yn, int zn,
 280                           int *pxp, int *pyp, int *pzp) ;
 281 static HISTOGRAM *gcaHistogramSamples(GCA *gca, GCA_SAMPLE *gcas, MRI *mri,
 282                                       TRANSFORM *transform, int nsamples,
 283                                       HISTOGRAM *histo, int frame) ;
 284 int GCApriorToNode(GCA *gca,
 285                    int xp, int yp, int zp,
 286                    int *pxn, int *pyn, int *pzn) ;
 287 
 288 /* arrays for indexing 6-connected neighbors */
 289 static int xnbr_offset[] =
 290   {
 291     1, -1, 0, 0,  0,  0
 292   } ;
 293 static int ynbr_offset[] =
 294   {
 295     0, 0,  1, -1, 0,  0
 296   } ;
 297 static int znbr_offset[] =
 298   {
 299     0, 0,  0, 0,  1, -1
 300   } ;
 301 int check_finite(char *where, double what) ;
 302 static int boundsCheck(int *pix, int *piy, int *piz, MRI *mri);
 303 
 304 static void  set_equilavent_classes(int *equivalent_classes);
 305 
 306 static int initialize_ventricle_alignment(MRI *mri_seg, 
 307                                           MRI *mri, 
 308                                           MATRIX *m_L, 
 309                                           char *base_name) ;
 310 
 311 
 312 void GCAsetVolGeom(GCA *gca, VOL_GEOM *vg)
 313 {
 314   vg->width = gca->width;
 315   vg->height = gca->height;
 316   vg->depth = gca->depth;
 317   vg->xsize = gca->xsize;
 318   vg->ysize = gca->ysize;
 319   vg->zsize = gca->zsize;
 320   vg->x_r = gca->x_r;
 321   vg->y_r = gca->y_r;
 322   vg->z_r = gca->z_r;
 323   vg->x_a = gca->x_a;
 324   vg->y_a = gca->y_a;
 325   vg->z_a = gca->z_a;
 326   vg->x_s = gca->x_s;
 327   vg->y_s = gca->y_s;
 328   vg->z_s = gca->z_s;
 329   vg->c_r = gca->c_r;
 330   vg->c_a = gca->c_a;
 331   vg->c_s = gca->c_s;
 332   vg->valid = 1;
 333 }
 334 
 335 void GCAcopyDCToMRI(GCA *gca, MRI *mri)
 336 {
 337   mri->x_r = gca->x_r;
 338   mri->y_r = gca->y_r;
 339   mri->z_r = gca->z_r;
 340   mri->x_a = gca->x_a;
 341   mri->y_a = gca->y_a;
 342   mri->z_a = gca->z_a;
 343   mri->x_s = gca->x_s;
 344   mri->y_s = gca->y_s;
 345   mri->z_s = gca->z_s;
 346   mri->c_r = gca->c_r;
 347   mri->c_a = gca->c_a;
 348   mri->c_s = gca->c_s;
 349   mri->ras_good_flag = 1;
 350     MRIreInitCache(mri) ;
 351   mri->i_to_r__ = extract_i_to_r(mri);
 352   mri->r_to_i__ = extract_r_to_i(mri);
 353 }
 354 void GCAcopyDCToGCA(GCA *gca, GCA *gca_dst)
 355 {
 356   gca_dst->x_r = gca->x_r;
 357   gca_dst->y_r = gca->y_r;
 358   gca_dst->z_r = gca->z_r;
 359   gca_dst->x_a = gca->x_a;
 360   gca_dst->y_a = gca->y_a;
 361   gca_dst->z_a = gca->z_a;
 362   gca_dst->x_s = gca->x_s;
 363   gca_dst->y_s = gca->y_s;
 364   gca_dst->z_s = gca->z_s;
 365   gca_dst->c_r = gca->c_r;
 366   gca_dst->c_a = gca->c_a;
 367   gca_dst->c_s = gca->c_s;
 368   //
 369   gca_dst->node_i_to_r__ =
 370     MatrixCopy(gca->node_i_to_r__, gca_dst->node_i_to_r__);
 371   gca_dst->node_r_to_i__ =
 372     MatrixCopy(gca->node_r_to_i__, gca_dst->node_r_to_i__);
 373   gca_dst->prior_i_to_r__ =
 374     MatrixCopy(gca->prior_i_to_r__, gca_dst->prior_i_to_r__);
 375   gca_dst->prior_r_to_i__ =
 376     MatrixCopy(gca->prior_r_to_i__, gca_dst->prior_r_to_i__);
 377   gca_dst->tal_i_to_r__ = MatrixCopy(gca->tal_i_to_r__, gca_dst->tal_i_to_r__);
 378   gca_dst->tal_r_to_i__ = MatrixCopy(gca->tal_r_to_i__, gca_dst->tal_r_to_i__);
 379   //
 380   gca_dst->mri_prior__ = MRIcopy(gca->mri_prior__, gca_dst->mri_prior__);
 381   gca_dst->mri_node__ = MRIcopy(gca->mri_node__, gca_dst->mri_node__);
 382   gca_dst->mri_tal__ = MRIcopy(gca->mri_tal__, gca_dst->mri_tal__);
 383 }
 384 
 385 void GCAcleanup(GCA *gca)
 386 {
 387   if (gca->mri_node__)
 388   {
 389     MRIfree(&gca->mri_node__);
 390     gca->mri_node__= 0;
 391   }
 392   if (gca->mri_prior__)
 393   {
 394     MRIfree(&gca->mri_prior__);
 395     gca->mri_prior__ = 0;
 396   }
 397   if (gca->mri_tal__)
 398   {
 399     MRIfree(&gca->mri_tal__);
 400     gca->mri_tal__ = 0;
 401   }
 402   if (gca->node_i_to_r__)
 403   {
 404     MatrixFree(&(gca->node_i_to_r__));
 405     gca->node_i_to_r__ = 0;
 406   }
 407   if (gca->node_r_to_i__)
 408   {
 409     MatrixFree(&(gca->node_r_to_i__));
 410     gca->node_r_to_i__ = 0;
 411   }
 412   if (gca->prior_i_to_r__)
 413   {
 414     MatrixFree(&(gca->prior_i_to_r__));
 415     gca->prior_i_to_r__ = 0;
 416   }
 417   if (gca->prior_r_to_i__)
 418   {
 419     MatrixFree(&(gca->prior_r_to_i__));
 420     gca->prior_r_to_i__ = 0;
 421   }
 422   if (gca->tal_i_to_r__)
 423   {
 424     MatrixFree(&(gca->tal_i_to_r__));
 425     gca->tal_i_to_r__ = 0;
 426   }
 427   if (gca->tal_r_to_i__)
 428   {
 429     MatrixFree(&(gca->tal_r_to_i__));
 430     gca->tal_r_to_i__ = 0;
 431   }
 432   if (gca->tmp__)
 433   {
 434     MatrixFree(&gca->tmp__);
 435     gca->tmp__ = 0;
 436   }
 437 }
 438 
 439 // set up mri's used in GCA
 440 void GCAsetup(GCA *gca)
 441 {
 442   // set up node part /////////////////////
 443   if (gca->mri_node__)
 444   {
 445     MRIfree(&gca->mri_node__);
 446     gca->mri_node__ = 0;
 447   }
 448   gca->mri_node__ = \
 449                     MRIallocHeader(gca->node_width,
 450                                    gca->node_height,
 451                                    gca->node_depth,
 452                                    MRI_UCHAR);
 453   /* Copy the voxel resolutions.  Set the defaults */
 454   gca->mri_node__->xsize = gca->xsize * gca->node_spacing;
 455   gca->mri_node__->ysize = gca->ysize * gca->node_spacing;
 456   gca->mri_node__->zsize = gca->zsize * gca->node_spacing;
 457 
 458   // this will recalculate i_to_r__ etc and
 459   // thus xsize etc must be set correctly
 460   GCAcopyDCToMRI(gca, gca->mri_node__);
 461 
 462   // fprintf(stdout, "node voxelToRAS\n");
 463   // MATRIX *mnode = extract_i_to_r(gca->mri_node__);
 464   // MatrixPrint(stdout, mnode);
 465   // MatrixFree(&mnode);
 466 
 467   // setup prior part //////////////////////////////////////
 468   if (gca->mri_prior__)
 469   {
 470     MRIfree(&gca->mri_prior__);
 471     gca->mri_prior__ = 0;
 472   }
 473   gca->mri_prior__ =
 474     MRIallocHeader(gca->prior_width,
 475                    gca->prior_height,
 476                    gca->prior_depth,
 477                    MRI_UCHAR);
 478 
 479   /* Copy the voxel resolutions.  Set the defaults */
 480   gca->mri_prior__->xsize = gca->xsize * gca->prior_spacing;
 481   gca->mri_prior__->ysize = gca->ysize * gca->prior_spacing;
 482   gca->mri_prior__->zsize = gca->zsize * gca->prior_spacing;
 483 
 484   GCAcopyDCToMRI(gca, gca->mri_prior__);
 485 
 486   // fprintf(stdout, "prior voxelToRAS\n");
 487   // MATRIX *mprior = extract_i_to_r(gca->mri_prior__);
 488   // MatrixPrint(stdout, mprior);
 489   // MatrixFree(&mprior);
 490 
 491   // set up the default talairach volume ////////////////
 492   if (gca->mri_tal__)
 493   {
 494     MRIfree(&gca->mri_tal__);
 495     gca->mri_tal__ = 0;
 496   }
 497   gca->mri_tal__ = 
 498     MRIallocHeader(gca->width, gca->height, gca->depth, MRI_UCHAR);
 499 
 500   /* Copy the voxel resolutions.  Set the defaults */
 501   gca->mri_tal__->xsize = gca->xsize;
 502   gca->mri_tal__->ysize = gca->ysize;
 503   gca->mri_tal__->zsize = gca->zsize;
 504 
 505   GCAcopyDCToMRI(gca, gca->mri_tal__);
 506 
 507   //fprintf(stdout, "tal voxelToRAS\n");
 508   //mtal = extract_i_to_r(gca->mri_tal__);
 509   // MatrixPrint(stdout, mtal);
 510   // MatrixFree(&mtal);
 511   if (gca->node_i_to_r__)
 512   {
 513     MatrixFree(&(gca->node_i_to_r__));
 514     gca->node_i_to_r__ = 0;
 515   }
 516   if (gca->node_r_to_i__)
 517   {
 518     MatrixFree(&(gca->node_r_to_i__));
 519     gca->node_r_to_i__ = 0;
 520   }
 521   if (gca->prior_i_to_r__)
 522   {
 523     MatrixFree(&(gca->prior_i_to_r__));
 524     gca->prior_i_to_r__ = 0;
 525   }
 526   if (gca->prior_r_to_i__)
 527   {
 528     MatrixFree(&(gca->prior_r_to_i__));
 529     gca->prior_r_to_i__ = 0;
 530   }
 531   if (gca->tal_i_to_r__)
 532   {
 533     MatrixFree(&(gca->tal_i_to_r__));
 534     gca->tal_i_to_r__ = 0;
 535   }
 536   if (gca->tal_r_to_i__)
 537   {
 538     MatrixFree(&(gca->tal_r_to_i__));
 539     gca->tal_r_to_i__ = 0;
 540   }
 541   if (gca->tmp__)
 542   {
 543     MatrixFree(&(gca->tmp__));
 544     gca->tmp__ = 0;
 545   }
 546   gca->node_i_to_r__ = extract_i_to_r(gca->mri_node__);
 547   gca->node_r_to_i__ = extract_r_to_i(gca->mri_node__);
 548   gca->prior_i_to_r__ = extract_i_to_r(gca->mri_prior__);
 549   gca->prior_r_to_i__ = extract_r_to_i(gca->mri_prior__);
 550   gca->tal_i_to_r__ = extract_i_to_r(gca->mri_tal__);
 551   gca->tal_r_to_i__ = extract_r_to_i(gca->mri_tal__);
 552   gca->tmp__ = MatrixAlloc(4,4, MATRIX_REAL);
 553 }
 554 
 555 // using the values of mri, modify gca
 556 // now we can keep track of voxelToRAS info even for non-standard type
 557 // This has the potential of changing direction cosines
 558 void GCAreinit(MRI *mri, GCA *gca)
 559 {
 560   // Keep the direction cosine the same to avoid used by
 561   // different mri
 562   // modify direction cosines etc.
 563   gca->x_r = mri->x_r; gca->y_r = mri->y_r; gca->z_r = mri->z_r;
 564   gca->x_a = mri->x_a; gca->y_a = mri->y_a; gca->z_a = mri->z_a;
 565   gca->x_s = mri->x_s; gca->y_s = mri->y_s; gca->z_s = mri->z_s;
 566   gca->c_r = mri->c_r; gca->c_a = mri->c_a; gca->c_s = mri->c_s;
 567 
 568   if (Gdiag & DIAG_SHOW && DIAG_VERBOSE_ON)
 569     printf("gca reinit c_(ras) = (%.2f, %.2f, %.2f)\n",
 570            gca->c_r, gca->c_a, gca->c_s);
 571 
 572   // modify width height depth
 573   if (gca->width != mri->width)
 574     fprintf(stdout, "gca width modified from %d to %d\n",
 575             gca->width, mri->width);
 576   if (gca->height != mri->height)
 577     fprintf(stdout, "gca height modified from %d to %d\n",
 578             gca->height, mri->height);
 579   if (gca->depth != mri->depth)
 580     fprintf(stdout, "gca depth modified from %d to %d\n",
 581             gca->depth, mri->depth);
 582   gca->width = mri->width;
 583   gca->height = mri->height;
 584   gca->depth = mri->depth;
 585 
 586   if (gca->xsize != mri->xsize)
 587     fprintf(stdout, "gca xsize modified from %.3f to %.3f\n",
 588             gca->xsize, mri->xsize);
 589   if (gca->ysize != mri->ysize)
 590     fprintf(stdout, "gca ysize modified from %.3f to %.3f\n",
 591             gca->ysize, mri->ysize);
 592   if (gca->zsize != mri->zsize)
 593     fprintf(stdout, "gca zsize modified from %.3f to %.3f\n",
 594             gca->zsize, mri->zsize);
 595   gca->xsize = mri->xsize;
 596   gca->ysize = mri->ysize;
 597   gca->zsize = mri->zsize;
 598 #if 0
 599     // can't do this without reallocating!!
 600   // then must modify node width etc.
 601   gca->node_width = (int)(((float)gca->width/gca->node_spacing)+.99) ;
 602   gca->node_height = (int)((float)gca->height/gca->node_spacing+.99) ;
 603   gca->node_depth = (int)(((float)gca->depth/gca->node_spacing)+.99) ;
 604   gca->prior_width = (int)(((float)gca->width/gca->prior_spacing)+.99) ;
 605   gca->prior_height = (int)((float)gca->height/gca->prior_spacing+.99) ;
 606   gca->prior_depth = (int)(((float)gca->depth/gca->prior_spacing)+.99) ;
 607 #endif
 608   //
 609   GCAsetup(gca);
 610   fflush(stdout) ;
 611 }
 612 
 613 GCA_PRIOR *
 614 getGCAP(GCA *gca, MRI *mri, TRANSFORM *transform, int xv, int yv, int zv)
 615 {
 616   int       xp, yp, zp ;
 617   GCA_PRIOR *gcap=NULL;
 618 
 619   if (!GCAsourceVoxelToPrior(gca, mri, transform, xv, yv, zv, &xp, &yp, &zp))
 620     gcap = &gca->priors[xp][yp][zp] ;
 621 
 622   return(gcap) ;
 623 }
 624 
 625 GCA_PRIOR *
 626 getGCAPfloat(GCA *gca, MRI *mri, TRANSFORM *transform, float xv, float yv, float zv)
 627 {
 628   int       xp, yp, zp ;
 629   GCA_PRIOR *gcap=NULL;
 630 
 631   if (!GCAsourceFloatVoxelToPrior(gca, mri, transform, xv, yv, zv, &xp, &yp, &zp))
 632     gcap = &gca->priors[xp][yp][zp] ;
 633 
 634   return(gcap) ;
 635 }
 636 
 637 GCA_NODE *
 638 getGCAN(GCA *gca, MRI *mri, TRANSFORM *transform, int xv, int yv, int zv)
 639 {
 640   int       xn, yn, zn ;
 641   GCA_NODE *gcan=NULL;
 642 
 643   if (!GCAsourceVoxelToNode(gca, mri, transform, xv, yv, zv, &xn, &yn, &zn))
 644     gcan = &gca->nodes[xn][yn][zn] ;
 645 
 646   return(gcan) ;
 647 }
 648 
 649 float
 650 getPrior(GCA_PRIOR *gcap, int label)
 651 {
 652   int n ;
 653 
 654   if (gcap==NULL)
 655     return (VERY_UNLIKELY);
 656 
 657   // find the label
 658   for (n = 0 ; n < gcap->nlabels ; n++)
 659     if (gcap->labels[n] == label)
 660       break ;
 661   // cannot find it
 662   if (n >= gcap->nlabels)
 663   {
 664     if (gcap->total_training > 0)
 665       return(0.1f/(float)gcap->total_training) ; /* make it unlikely */
 666     else
 667       return(VERY_UNLIKELY) ;
 668   }
 669   // return found one
 670   return(gcap->priors[n]) ;
 671 }
 672 
 673 int
 674 GCAisPossible(GCA *gca, MRI *mri, int label,
 675               TRANSFORM *transform, int x, int y, int z,
 676               int use_mrf)
 677 {
 678   int       n, i, found, nbr_label, xnbr, ynbr, znbr, xn, yn, zn ;
 679   GCA_PRIOR *gcap = NULL;
 680   GC1D      *gc ;
 681 
 682   gcap = getGCAP(gca, mri, transform, x, y, z) ;
 683   if (gcap==NULL)
 684     return(0) ;
 685   // if label found
 686   for (found = n = 0 ; n < gcap->nlabels ; n++)
 687     if (gcap->labels[n] == label)
 688     {
 689       found = 1 ;
 690       break ;
 691     }
 692   if (found == 0)
 693     return(0) ;
 694   if (use_mrf == 0)
 695     return(1) ;
 696 
 697   if (GCAsourceVoxelToNode(gca, mri, transform, x, y, z, &xn, &yn, &zn)
 698       != NO_ERROR)
 699     return(1) ;
 700   gc = GCAfindGC(gca, xn, yn, zn, label) ;
 701   for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
 702   {
 703     found = 0 ;
 704     xnbr = mri->xi[x+xnbr_offset[i]] ;// xnbr_offset 1, -1, 0,  0, 0,  0
 705     ynbr = mri->yi[y+ynbr_offset[i]] ;// ynbr_offset 0,  0, 1, -1, 0,  0
 706     znbr = mri->zi[z+znbr_offset[i]] ;// znbr_offset 0,  0, 0,  0, 1, -1
 707     nbr_label = nint(MRIgetVoxVal(mri, xnbr, ynbr, znbr, 0)) ;
 708     for (n = 0 ; n < gc->nlabels[i] ; n++)
 709       if (gc->labels[i][n] == nbr_label)
 710       {
 711         found = 1 ;
 712         break ;
 713       }
 714     if (found == 0)  // this pair never occurred here
 715       return(0) ;
 716   }
 717   return(1) ;
 718 }
 719 
 720 #if 0
 721 static float get_voxel_prior(GCA *gca, MRI *mri, TRANSFORM *transform,
 722                              int xv, int yv, int zv, int label);
 723 static float
 724 get_voxel_prior(GCA *gca, MRI *mri, TRANSFORM *transform,
 725                 int xv, int yv, int zv, int label)
 726 {
 727   int       xn, yn, zn ;
 728   GCA_PRIOR *gcap ;
 729 
 730   gcap = getGCAP(gca, mri, transform, xv, yv, zv) ;
 731   // this could return NULL
 732   if (gcap==NULL)
 733     return (VERY_UNLIKELY)
 734 
 735            return(getPrior(gcap, label)) ;
 736 }
 737 #endif
 738 
 739 static float
 740 get_node_prior(GCA *gca, int label, int xn, int yn, int zn)
 741 {
 742   int       xp, yp, zp ;
 743   GCA_PRIOR *gcap ;
 744 
 745   if (gcaNodeToPrior(gca, xn, yn, zn, &xp, &yp, &zp)==NO_ERROR)
 746   {
 747     gcap = &gca->priors[xp][yp][zp] ;
 748     if (gcap == NULL)
 749       return (VERY_UNLIKELY);
 750     return(getPrior(gcap, label)) ;
 751   }
 752   else
 753     return (VERY_UNLIKELY);
 754 }
 755 
 756 // use the same function for bounds checking
 757 // if (ix, iy, iz) is within mri volume, returns NO_ERROR
 758 //                                       otherwise ERROR_BADPARAM
 759 static int boundsCheck(int *pix, int *piy, int *piz, MRI *mri)
 760 {
 761   int ix = *pix;
 762   int iy = *piy;
 763   int iz = *piz;
 764   int errCode = NO_ERROR; // initialize
 765 
 766   if (ix < 0)
 767     errCode = ERROR_BADPARM;
 768   else if (iy < 0)
 769     errCode = ERROR_BADPARM;
 770   else if (iz < 0)
 771     errCode = ERROR_BADPARM;
 772   else if (ix > mri->width-1)
 773     errCode = ERROR_BADPARM;
 774   else if (iy > mri->height-1)
 775     errCode = ERROR_BADPARM;
 776   else if (iz > mri->depth-1)
 777     errCode = ERROR_BADPARM;
 778 
 779 #if 0
 780   // give up returning error
 781   if (ix < 0)
 782     *pix = 0;
 783   if (iy < 0)
 784     *piy = 0;
 785   if (iz < 0)
 786     *piz = 0;
 787   if (ix > mri->width -1)
 788     *pix = mri->width -1;
 789   if (iy > mri->height-1)
 790     *piy = mri->height -1;
 791   if (iz > mri->depth-1)
 792     *piz = mri->depth -1;
 793 #endif
 794 
 795   return errCode;
 796 }
 797 
 798 static int
 799 gcaNodeToPrior(GCA *gca, int xn, int yn, int zn, int *pxp, int *pyp, int *pzp)
 800 {
 801   // initialize errCode
 802   int errCode = NO_ERROR;
 803 
 804   // see GCApriorToNode comment
 805   // node_spacing > prior_spacing
 806   // integer operation is floor
 807   *pxp = xn*gca->node_spacing/gca->prior_spacing;
 808   *pyp = yn*gca->node_spacing/gca->prior_spacing;
 809   *pzp = zn*gca->node_spacing/gca->prior_spacing;
 810   // no error should occur
 811   return errCode ;
 812 }
 813 
 814 int
 815 GCApriorToNode(GCA *gca, int xp, int yp, int zp, int *pxn, int *pyn, int *pzn)
 816 {
 817   int errCode = NO_ERROR;
 818   /////////////////////////////////////////////////////////////////////
 819   // Note that prior and node share the common direction cosines/translation
 820   //
 821   //      prior has voxelToRAS = [ R | T ][prior_size | 0] = X [prior_size | 0]
 822   //                             [ 0 | 1 ][    0      | 1]     [    0      | 1]
 823   //      node  has voxelToRAS = [ R | T ][ node_size | 0] = X [ node_size | 0]
 824   //                             [ 0 | 1 ][    0      | 1]     [    0      | 1]
 825   //
 826   //        prior   --->  RAS
 827   //          |            |
 828   //          |            | identity
 829   //          V            V
 830   //         node    ---> RAS
 831   //                                  -1   -1
 832   //    priorToNode = [ node_size | 0]  * X   * X * [ prior_size | 0]
 833   //                  [     0     | 1]              [     0      | 1]
 834   //
 835   //                = [         -1             |     ]
 836   //                  [ node_size * prior_size |  0  ]
 837   //                  [           0            |  1  ]
 838   //
 839   // Note that node_spacing > prior_spacing and thus the following will do
 840   // .5 becomes 0.
 841   // but this is OK, since the correspondence is
 842   //               |  0  |   1  |   prior
 843   //               |     0      |   node
 844   // i.e. 1 becomes 0
 845   *pxn = xp*gca->prior_spacing/gca->node_spacing;
 846   *pyn = yp*gca->prior_spacing/gca->node_spacing;
 847   *pzn = zp*gca->prior_spacing/gca->node_spacing;
 848   // no error should occur
 849   return errCode;
 850 }
 851 
 852 static int
 853 dump_gcan(GCA *gca, GCA_NODE *gcan, FILE *fp, int verbose, GCA_PRIOR *gcap)
 854 {
 855   int       n, i, j, n1 ;
 856   GC1D      *gc ;
 857   float     prior ;
 858   VECTOR    *v_means = NULL ;
 859   MATRIX    *m_cov = NULL ;
 860 
 861   if (gcap==NULL)
 862     return -1;
 863   // print out gcan labels
 864   fprintf(fp,
 865           "node labels at this point -----"
 866           "---------------------------------\n");
 867   for (n = 0 ; n < gcan->nlabels ; n++)
 868   {
 869     // find the same label in gcap
 870     for (n1 = 0 ; n1 < gcap->nlabels ; n1++)
 871       if (gcap->labels[n1] == gcan->labels[n])
 872         break ;
 873     // the following won't happen
 874     if (n1 >= gcap->nlabels)
 875       continue ;
 876     prior = getPrior(gcap, gcan->labels[n]) ;
 877     if (FZERO(prior))
 878       continue ;
 879     // get gc for this label
 880     gc = &gcan->gcs[n] ;
 881     v_means = load_mean_vector(gc, v_means, gca->ninputs) ;
 882     m_cov = load_covariance_matrix(gc, m_cov, gca->ninputs) ;
 883     fprintf(fp, "%d: label %s (%d), prior %2.3f, ntr %d, means " ,
 884             n, cma_label_to_name(gcan->labels[n]),
 885             gcan->labels[n], prior, gc->ntraining) ;
 886     for (i = 0 ; i < gca->ninputs ; i++)
 887       fprintf(fp, "%2.1f ", VECTOR_ELT(v_means,i+1)) ;
 888     fprintf(fp, ", cov:\n") ;
 889     MatrixPrint(fp, m_cov) ;
 890 
 891     if (verbose)
 892     {
 893       for (i = 0 ; i < gca->ninputs ; i++)
 894       {
 895         for (j = 0 ; j < gca->ninputs ; j++)
 896         {
 897           fprintf(fp, "%2.1f\t", *MATRIX_RELT(m_cov,i+1,j+1)) ;
 898         }
 899         fprintf(fp, "\n") ;
 900       }
 901       // 6 neighbors
 902       for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
 903       {
 904         fprintf(fp, "\tnbr %d (%d,%d,%d): %d labels\n",
 905                 i, xnbr_offset[i], ynbr_offset[i], znbr_offset[i],
 906                 gc->nlabels[i]) ;
 907         fflush(fp);
 908         for (j = 0 ; j < gc->nlabels[i] ; j++)
 909         {
 910           fprintf(fp, "\t\tlabel %s, prior %2.3f\n",
 911                   cma_label_to_name(gc->labels[i][j]),
 912                   gc->label_priors[i][j]) ;
 913           fflush(fp);
 914         }
 915       }
 916     }
 917   }
 918   fprintf(fp, "--------------------------------"
 919           "--------------------------------\n");
 920   MatrixFree(&m_cov) ;
 921   VectorFree(&v_means);
 922   return(NO_ERROR) ;
 923 }
 924 
 925 ////////////////////////////////////////////////////////////////
 926 // transform from template -> node
 927 ////////////////////////////////////////////////////////////////
 928 int GCAvoxelToNodeReal(GCA *gca, MRI *mri, Real xv, Real yv, Real zv,
 929                        Real *pxn, Real *pyn, Real *pzn)
 930 {
 931   //               i_to_r
 932   //        mri    ----->    RAS
 933   //         |                |
 934   //         | we need        | identity
 935   //         V                V
 936   //        node   <-----    RAS
 937   //               r_to_i
 938   MATRIX *rasFromVoxel = mri->i_to_r__; // extract_i_to_r(mri);
 939   MATRIX *nodeFromRAS = gca->node_r_to_i__; // extract_r_to_i(gca->mri_node__);
 940   MATRIX *voxelToNode = gca->tmp__;
 941   MatrixMultiply(nodeFromRAS, rasFromVoxel, gca->tmp__);
 942 
 943   TransformWithMatrix(voxelToNode, xv, yv, zv, pxn, pyn, pzn);
 944 
 945   // MatrixFree(&rasFromVoxel);
 946   // MatrixFree(&nodeFromRAS);
 947   // MatrixFree(&voxelToNode);
 948 
 949   return NO_ERROR;
 950 }
 951 
 952 int
 953 GCAvoxelToNode(GCA *gca, MRI *mri, int xv, int yv, int zv, int *pxn,
 954                int *pyn, int *pzn)
 955 {
 956   Real xn, yn, zn;
 957   int   ixn, iyn, izn;
 958   int errCode = NO_ERROR;
 959 
 960   GCAvoxelToNodeReal(gca, mri, xv, yv, zv, &xn, &yn, &zn);
 961 
 962   // xn, yn, zn are double.  we use voxel center as integer
 963   ixn = (int) floor(xn);
 964   iyn = (int) floor(yn);
 965   izn = (int) floor(zn);
 966   // if outofbounds, tell it
 967   errCode = boundsCheck(&ixn, &iyn, &izn, gca->mri_node__);
 968   //
 969   *pxn = ixn;
 970   *pyn = iyn;
 971   *pzn = izn;
 972 
 973   return errCode;
 974 }
 975 
 976 ////////////////////////////////////////////////////////////////////
 977 // transform from template -> prior
 978 ////////////////////////////////////////////////////////////////////
 979 int GCAvoxelToPriorReal(GCA *gca, MRI *mri, Real xv, Real yv, Real zv,
 980                         Real *pxp, Real *pyp, Real *pzp)
 981 {
 982   MATRIX *rasFromVoxel = mri->i_to_r__; //extract_i_to_r(mri);
 983   MATRIX *priorFromRAS = gca->prior_r_to_i__;
 984   //extract_r_to_i(gca->mri_prior__);
 985   MATRIX *voxelToPrior = gca->tmp__;
 986   MatrixMultiply(priorFromRAS, rasFromVoxel, gca->tmp__);
 987 
 988   TransformWithMatrix(voxelToPrior, xv, yv, zv, pxp, pyp, pzp);
 989 
 990   // MatrixFree(&rasFromVoxel);
 991   // MatrixFree(&priorFromRAS);
 992   // MatrixFree(&voxelToPrior);
 993 
 994   return NO_ERROR;
 995 }
 996 
 997 int
 998 GCAvoxelToPrior(GCA *gca, MRI *mri, int xv, int yv, int zv,
 999                 int *pxp, int *pyp, int *pzp)
1000 {
1001   Real xp, yp, zp;
1002   int ixp, iyp, izp;
1003   int errCode = NO_ERROR;
1004 
1005   GCAvoxelToPriorReal(gca, mri, xv, yv, zv, &xp, &yp, &zp);
1006 
1007   ixp = (int) floor(xp);
1008   iyp = (int) floor(yp);
1009   izp = (int) floor(zp);
1010   // bound check
1011   // if outofbounds, tell it
1012   errCode = boundsCheck(&ixp, &iyp, &izp, gca->mri_prior__);
1013   //
1014   *pxp = ixp;
1015   *pyp = iyp;
1016   *pzp = izp;
1017 
1018   return errCode;
1019 }
1020 
1021 /////////////////////////////////////////////////////////////////////
1022 // transform node->template
1023 /////////////////////////////////////////////////////////////////////
1024 int GCAnodeToVoxelReal(GCA *gca, MRI *mri, Real xn, Real yn, Real zn,
1025                        Real *pxv, Real *pyv, Real *pzv)
1026 {
1027   //               r_to_i
1028   //        mri    <-----    RAS
1029   //         ^                ^
1030   //         | we need        | identity
1031   //         |                |
1032   //        node   ----->    RAS
1033   //               i_to_r
1034   MATRIX *rasFromNode = gca->node_i_to_r__;
1035   //  extract_i_to_r(gca->mri_node__);
1036   MATRIX *voxelFromRAS = mri->r_to_i__; //  extract_r_to_i(mri);
1037   MATRIX *nodeToVoxel = gca->tmp__;
1038   MatrixMultiply(voxelFromRAS, rasFromNode, gca->tmp__);
1039 
1040   TransformWithMatrix(nodeToVoxel, xn, yn, zn, pxv, pyv, pzv);
1041 
1042   // MatrixFree(&rasFromNode);
1043   // MatrixFree(&voxelFromRAS);
1044   // MatrixFree(&nodeToVoxel);
1045 
1046   return NO_ERROR;
1047 }
1048 
1049 int
1050 GCAnodeToVoxel(GCA *gca, MRI *mri, int xn, int yn, int zn,
1051                int *pxv, int *pyv, int *pzv)
1052 {
1053   Real xv, yv, zv;
1054   int ixv, iyv, izv;
1055   int errCode = NO_ERROR;
1056 
1057   GCAnodeToVoxelReal(gca, mri, xn, yn, zn, &xv, &yv, &zv);
1058 
1059   // addition makes overall error smaller
1060   //             0 picks  1 out of 0, 1, 2, 3 possible choices
1061   // without it, 0 pickes 0 out of 0, 1, 2, 3 possible choices
1062   // if you used float, then 0 picks 1.5.
1063   ixv = (int) floor(xv + gca->node_spacing/2.);
1064   iyv = (int) floor(yv + gca->node_spacing/2.);
1065   izv = (int) floor(zv + gca->node_spacing/2.);
1066 
1067   // bound check
1068   // if outofbounds, tell it
1069   errCode = boundsCheck(&ixv, &iyv, &izv, mri);
1070   //
1071 
1072   *pxv = ixv;
1073   *pyv = iyv;
1074   *pzv = izv;
1075   return errCode;
1076 }
1077 
1078 //////////////////////////////////////////////////////////////////////
1079 // transform from prior-> template
1080 //////////////////////////////////////////////////////////////////////
1081 int GCApriorToVoxelReal(GCA *gca, MRI *mri, Real xp, Real yp, Real zp,
1082                         Real *pxv, Real *pyv, Real *pzv)
1083 {
1084   MATRIX *rasFromPrior = gca->prior_i_to_r__;
1085   // extract_i_to_r(gca->mri_prior__);
1086   MATRIX *voxelFromRAS = mri->r_to_i__; // extract_r_to_i(mri);
1087   MATRIX *priorToVoxel = gca->tmp__;
1088   MatrixMultiply(voxelFromRAS, rasFromPrior, gca->tmp__);
1089 
1090   // TransformWithMatrix(priorToVoxel, xp, yp, zp, pxv, pyv, pzv);
1091   TransformWithMatrix(priorToVoxel, xp, yp, zp, pxv, pyv, pzv);
1092 
1093   // MatrixFree(&rasFromPrior);
1094   // MatrixFree(&voxelFromRAS);
1095   // MatrixFree(&priorToVoxel);
1096 
1097   return NO_ERROR;
1098 }
1099 
1100 int
1101 GCApriorToVoxel(GCA *gca, MRI *mri, int xp, int yp, int zp,
1102                 int *pxv, int *pyv, int *pzv)
1103 {
1104   Real xv, yv, zv;
1105   int ixv, iyv, izv;
1106   int errCode = NO_ERROR;
1107 
1108   GCApriorToVoxelReal(gca, mri, xp, yp, zp, &xv, &yv, &zv);
1109   // addition makes overall error smaller
1110   // without it, 0 pickes 0 out of 0, 1 possible choices
1111   ixv = (int) floor(xv + gca->prior_spacing/2.);
1112   iyv = (int) floor(yv + gca->prior_spacing/2.);
1113   izv = (int) floor(zv + gca->prior_spacing/2.);
1114   // bound check
1115   // if outofbounds, tell it
1116   errCode = boundsCheck(&ixv, &iyv, &izv, mri);
1117   //
1118   *pxv = ixv;
1119   *pyv = iyv;
1120   *pzv = izv;
1121 
1122   errCode = NO_ERROR;
1123 
1124   return errCode;
1125 }
1126 
1127 ///////////////////////////////////////////////////////////////////////
1128 // transform from source -> template space -> prior
1129 //////////////////////////////////////////////////////////////////////
1130 int
1131 GCAsourceVoxelToPrior(GCA *gca, MRI *mri, TRANSFORM *transform,
1132                       int xv, int yv, int zv, int *pxp, int *pyp, int *pzp)
1133 {
1134   float   xt, yt, zt ;
1135   Real    xrt, yrt, zrt, xrp, yrp, zrp;
1136 
1137   LTA *lta;
1138   if (transform->type != MORPH_3D_TYPE)
1139   {
1140     if (transform->type == LINEAR_VOX_TO_VOX)
1141     {
1142       lta = (LTA *) transform->xform;
1143       // transform point to talairach volume point
1144       TransformWithMatrix(lta->xforms[0].m_L,
1145                           xv, yv, zv, &xrt, &yrt, &zrt);
1146       xt = xrt;
1147       yt = yrt;
1148       zt = zrt;
1149       // TransformSample(transform, xv, yv, zv, &xt, &yt, &zt) ;
1150     }
1151     else
1152       ErrorExit(ERROR_BADPARM, \
1153                 "GCAsourceVoxelToPrior: needs vox-to-vox transform") ;
1154   }
1155   else // morph 3d type can go directly from source to template
1156   {
1157     TransformSample(transform, xv, yv, zv, &xt, &yt, &zt);
1158   }
1159   // get the position in gca from talairach volume
1160   GCAvoxelToPriorReal(gca, gca->mri_tal__, xt, yt, zt, &xrp, &yrp, &zrp) ;
1161   *pxp = nint(xrp) ;
1162   *pyp = nint(yrp) ;
1163   *pzp = nint(zrp) ;
1164   if (*pxp < 0 || *pyp < 0 || *pzp < 0 ||
1165       *pxp >= gca->prior_width ||
1166       *pyp >= gca->prior_height ||
1167       *pzp >= gca->prior_depth)
1168     return(ERROR_BADPARM) ;
1169   return (NO_ERROR) ;
1170 }
1171 
1172 int
1173 GCAsourceFloatVoxelToPrior(GCA *gca, MRI *mri, TRANSFORM *transform,
1174                            float xv, float yv, float zv, int *pxp, int *pyp, int *pzp)
1175 {
1176   float   xt, yt, zt ;
1177   Real    xrt, yrt, zrt, xrp, yrp, zrp;
1178 
1179   LTA *lta;
1180   if (transform->type != MORPH_3D_TYPE)
1181   {
1182     if (transform->type == LINEAR_VOX_TO_VOX)
1183     {
1184       lta = (LTA *) transform->xform;
1185       // transform point to talairach volume point
1186       TransformWithMatrix(lta->xforms[0].m_L,
1187                           xv, yv, zv, &xrt, &yrt, &zrt);
1188       xt = xrt;
1189       yt = yrt;
1190       zt = zrt;
1191       // TransformSample(transform, xv, yv, zv, &xt, &yt, &zt) ;
1192     }
1193     else
1194       ErrorExit(ERROR_BADPARM, \
1195                 "GCAsourceVoxelToPrior: needs vox-to-vox transform") ;
1196   }
1197   else // morph 3d type can go directly from source to template
1198   {
1199     TransformSample(transform, xv, yv, zv, &xt, &yt, &zt);
1200   }
1201   // get the position in gca from talairach volume
1202   GCAvoxelToPriorReal(gca, gca->mri_tal__, xt, yt, zt, &xrp, &yrp, &zrp) ;
1203   *pxp = nint(xrp) ;
1204   *pyp = nint(yrp) ;
1205   *pzp = nint(zrp) ;
1206   if (*pxp < 0 || *pyp < 0 || *pzp < 0 ||
1207       *pxp >= gca->prior_width ||
1208       *pyp >= gca->prior_height ||
1209       *pzp >= gca->prior_depth)
1210     return(ERROR_BADPARM) ;
1211   return (NO_ERROR) ;
1212 }
1213 
1214 /////////////////////////////////////////////////////////////////////
1215 // transform from source -> template space -> prior
1216 /////////////////////////////////////////////////////////////////////
1217 int
1218 GCAsourceVoxelToNode(GCA *gca, MRI *mri, TRANSFORM *transform,
1219                      int xv, int yv, int zv,
1220                      int *pxn, int *pyn, int *pzn)
1221 {
1222   float xt, yt, zt;
1223   Real  xrt, yrt, zrt ;
1224   LTA *lta;
1225 
1226   if (transform->type != MORPH_3D_TYPE)
1227   {
1228     if (transform->type == LINEAR_VOX_TO_VOX) // from src to talairach volume
1229     {
1230       lta = (LTA *) transform->xform;
1231       // get the talairach position
1232       TransformWithMatrix(lta->xforms[0].m_L,
1233                           xv, yv, zv, &xrt, &yrt, &zrt);
1234       // TransformSample(transform, xv, yv, zv, &xt, &yt, &zt) ;
1235       xt = xrt;
1236       yt = yrt;
1237       zt = zrt;
1238     }
1239     else
1240       ErrorExit(ERROR_BADPARM,
1241                 "GCAsourceVoxelToNode: needs vox-to-vox transform") ;
1242   }
1243   else
1244   {
1245     TransformSample(transform, xv, yv, zv, &xt, &yt, &zt);
1246   }
1247   if (Ggca_x == xv && Ggca_y == yv && Ggca_z == zv && DIAG_VERBOSE_ON)
1248     fprintf(stdout, "source (%d, %d, %d) to talposition (%.2f, %.2f, %.2f)\n",
1249             xv, yv, zv, xt, yt, zt);
1250   fflush(stdout) ;
1251 
1252   // get the position in node from the talairach position
1253   return GCAvoxelToNode(gca, gca->mri_tal__, xt, yt, zt, pxn, pyn, pzn) ;
1254 }
1255 
1256 //////////////////////////////////////////////////////////////////////
1257 // transform from node->template->source
1258 ////////////////////////////////////////////////////////////////////
1259 int
1260 GCApriorToSourceVoxelFloat(GCA *gca, MRI *mri, TRANSFORM *transform,
1261                            int xp, int yp, int zp,
1262                            float *pxv, float *pyv, float *pzv)
1263 {
1264   int   width, height, depth;
1265   Real  xt, yt, zt ;
1266   float  xv, yv, zv ;
1267   Real  xc, yc, zc;
1268   int errCode = NO_ERROR;
1269   LTA *lta;
1270   width = mri->width ;
1271   height = mri->height ;
1272   depth = mri->depth ;
1273   // go to the template voxel position
1274   GCApriorToVoxelReal(gca, gca->mri_tal__, xp, yp, zp, &xt, &yt, &zt);
1275   // got the point in gca->mri_tal__ position
1276   if (transform->type != MORPH_3D_TYPE)
1277   {
1278     if (transform->type == LINEAR_VOX_TO_VOX) // from src to talairach volume
1279     {
1280       lta = (LTA *) transform->xform;
1281       // get the talairach to orig
1282       TransformWithMatrix(lta->inv_xforms[0].m_L,
1283                           xt, yt, zt, &xc, &yc, &zc);
1284       // TransformSampleInverse(transform, xt, yt, zt, &xc, &yc, &zc);
1285       if (xc < 0)
1286         errCode = ERROR_BADPARM;
1287       else if (yc < 0)
1288         errCode = ERROR_BADPARM;
1289       else if (zc < 0)
1290         errCode = ERROR_BADPARM;
1291       else if (xc > (width-1))
1292         errCode = ERROR_BADPARM;
1293       else if (yc > (height-1))
1294         errCode = ERROR_BADPARM;
1295       else if (zc > (depth-1))
1296         errCode = ERROR_BADPARM;
1297       xv = xc;
1298       yv = yc;
1299       zv = zc;
1300     }
1301     else
1302       ErrorExit(ERROR_BADPARM,
1303                 "GCApriorToSourceVoxelFloat: needs vox-to-vox transform") ;
1304   }
1305   else // go directly from template to source
1306   {
1307     TransformSampleInverse(transform, xt, yt, zt, &xv, &yv, &zv);
1308   }
1309   *pxv = xv ;
1310   *pyv = yv ;
1311   *pzv = zv ;
1312   return errCode ;
1313 }
1314 
1315 int
1316 GCAnodeToSourceVoxelFloat(GCA *gca, MRI *mri, TRANSFORM *transform,
1317                           int xn, int yn, int zn,
1318                           float *pxv, float *pyv, float *pzv)
1319 {
1320   int   width, height, depth;
1321   Real  xt, yt, zt ;
1322   float  xv, yv, zv ;
1323   Real  xc, yc, zc ;
1324   int errCode = NO_ERROR;
1325   LTA *lta;
1326   width = mri->width ;
1327   height = mri->height ;
1328   depth = mri->depth ;
1329   // get template voxel position
1330   GCAnodeToVoxelReal(gca, gca->mri_tal__, xn, yn, zn, &xt, &yt, &zt) ;
1331   if (transform->type != MORPH_3D_TYPE)
1332   {
1333     lta = (LTA *) transform->xform;
1334     // get the talairach to orig
1335     TransformWithMatrix(lta->inv_xforms[0].m_L, xt, yt, zt, &xc, &yc, &zc);
1336     // TransformSampleInverse(transform, xt, yt, zt, &xc, &yc, &zc);
1337     if (xc < 0)
1338       errCode = ERROR_BADPARM;
1339     else if (yc < 0)
1340       errCode = ERROR_BADPARM;
1341     else if (zc < 0)
1342       errCode = ERROR_BADPARM;
1343     else if (xc > (width-1))
1344       errCode = ERROR_BADPARM;
1345     else if (yc > (height-1))
1346       errCode = ERROR_BADPARM;
1347     else if (zc > (depth-1))
1348       errCode = ERROR_BADPARM;
1349     xv = xc;
1350     yv = yc;
1351     zv = zc;
1352   }
1353   else  // template to source
1354   {
1355     TransformSampleInverse(transform, xt, yt, zt, &xv, &yv, &zv);
1356   }
1357   *pxv = xv ;
1358   *pyv = yv ;
1359   *pzv = zv ;
1360   return errCode;
1361 }
1362 
1363 int
1364 GCAnodeToSourceVoxel(GCA *gca, MRI *mri, TRANSFORM *transform,
1365                      int xn, int yn, int zn,
1366                      int *pxv, int *pyv, int *pzv)
1367 {
1368   float  xf, yf, zf ;
1369   int errCode = NO_ERROR;
1370   errCode = GCAnodeToSourceVoxelFloat(gca, mri, transform,
1371                                       xn, yn, zn, &xf, &yf, &zf) ;
1372   if (xf < 0)
1373     errCode = ERROR_BADPARM;
1374   else if (yf <0)
1375     errCode = ERROR_BADPARM;
1376   else if (zf < 0)
1377     errCode = ERROR_BADPARM;
1378   else if (xf > (mri->width-1))
1379     errCode = ERROR_BADPARM;
1380   else if (yf > (mri->height-1))
1381     errCode = ERROR_BADPARM;
1382   else if (zf > (mri->depth-1))
1383     errCode = ERROR_BADPARM;
1384 
1385   *pxv = nint(xf) ;
1386   *pyv = nint(yf) ;
1387   *pzv = nint(zf) ;
1388   return errCode;
1389 }
1390 
1391 int
1392 GCApriorToSourceVoxel(GCA *gca, MRI *mri, TRANSFORM *transform,
1393                       int xp, int yp, int zp,
1394                       int *pxv, int *pyv, int *pzv)
1395 {
1396   float  xf, yf, zf ;
1397   int errCode = NO_ERROR;
1398 
1399   errCode = GCApriorToSourceVoxelFloat(gca, mri, transform,
1400                                        xp, yp, zp, &xf, &yf, &zf) ;
1401   if (xf < 0)
1402     errCode = ERROR_BADPARM;
1403   else if (yf < 0)
1404     errCode = ERROR_BADPARM;
1405   else if (zf < 0)
1406     errCode = ERROR_BADPARM;
1407   else if (xf > (mri->width-1))
1408     errCode = ERROR_BADPARM;
1409   else if (yf > (mri->height-1))
1410     errCode = ERROR_BADPARM;
1411   else if (zf > (mri->depth-1))
1412     errCode = ERROR_BADPARM;
1413 
1414   *pxv = nint(xf) ;
1415   *pyv = nint(yf) ;
1416   *pzv = nint(zf) ;
1417   return errCode;
1418 }
1419 
1420 GCA  *
1421 GCAalloc(int ninputs,
1422          float prior_spacing, float node_spacing,
1423          int width, int height, int depth, int flags)
1424 {
1425   int max_labels ;
1426 
1427   if (flags & GCA_NO_GCS)
1428     max_labels = 0 ;
1429   else
1430     max_labels = DEFAULT_MAX_LABELS_PER_GCAN ;
1431   return(gcaAllocMax(ninputs, prior_spacing, node_spacing, width,height,depth,
1432                      max_labels, flags));
1433 }
1434 
1435 GCA *
1436 gcaAllocMax(int ninputs,
1437             float prior_spacing, float node_spacing,
1438             int width, int height, int depth,
1439             int max_labels, int flags)
1440 {
1441   GCA       *gca ;
1442   GCA_NODE  *gcan ;
1443   GCA_PRIOR *gcap ;
1444   int       x, y, z ;
1445 
1446   gca = calloc(1, sizeof(GCA)) ;
1447   if (!gca)
1448     ErrorExit(ERROR_NOMEMORY, "GCAalloc: could not allocate struct") ;
1449 
1450   gca->ninputs = ninputs ;
1451   gca->prior_spacing = prior_spacing ;
1452   gca->node_spacing = node_spacing ;
1453   gca->type = GCA_UNKNOWN; // mark it as unknown
1454 
1455   // setup default direction cosines
1456   gca->x_r = -1;
1457   gca->y_r =  0;
1458   gca->z_r = 0;
1459   gca->c_r = 0;
1460   gca->x_a =  0;
1461   gca->y_a =  0;
1462   gca->z_a = 1;
1463   gca->c_a = 0;
1464   gca->x_s =  0;
1465   gca->y_s = -1;
1466   gca->z_s = 0;
1467   gca->c_s = 0;
1468   gca->xsize = 1;
1469   gca->ysize = 1;
1470   gca->zsize = 1;
1471   //
1472   gca->width = width;
1473   gca->height = height;
1474   gca->depth = depth;
1475 
1476   /* ceil gives crazy results, I don't know why */
1477   gca->node_width = (int)(((float)width/node_spacing)+.99) ;
1478   gca->node_height = (int)((float)height/node_spacing+.99) ;
1479   gca->node_depth = (int)(((float)depth/node_spacing)+.99) ;
1480   gca->prior_width = (int)(((float)width/prior_spacing)+.99) ;
1481   gca->prior_height = (int)((float)height/prior_spacing+.99) ;
1482   gca->prior_depth = (int)(((float)depth/prior_spacing)+.99) ;
1483   gca->flags = flags ;
1484 
1485   gca->nodes = (GCA_NODE ***)calloc(gca->node_width, sizeof(GCA_NODE **)) ;
1486   if (!gca->nodes)
1487     ErrorExit(ERROR_NOMEMORY, "GCAalloc: could not allocate nodes") ;
1488 
1489   // setting vlaues gca->nodes volume
1490   for (x = 0 ; x < gca->node_width ; x++)
1491   {
1492     gca->nodes[x] =
1493       (GCA_NODE **)calloc(gca->node_height, sizeof(GCA_NODE *)) ;
1494     if (!gca->nodes[x])
1495       ErrorExit(ERROR_NOMEMORY, "GCAalloc: could not allocate %dth **",x) ;
1496 
1497     for (y = 0 ; y < gca->node_height ; y++)
1498     {
1499       gca->nodes[x][y] =
1500         (GCA_NODE *)calloc(gca->node_depth, sizeof(GCA_NODE)) ;
1501       if (!gca->nodes[x][y])
1502         ErrorExit(ERROR_NOMEMORY,
1503                   "GCAalloc: could not allocate %d,%dth *",x,y);
1504       for (z = 0 ; z < gca->node_depth ; z++)
1505       {
1506         gcan = &gca->nodes[x][y][z] ;
1507         // set max_labels
1508         gcan->max_labels = max_labels ;
1509 
1510         if (max_labels > 0)
1511         {
1512           /* allocate new ones */
1513           gcan->gcs =
1514             alloc_gcs(gcan->max_labels, flags, gca->ninputs) ;
1515           if (!gcan->gcs)
1516             ErrorExit(ERROR_NOMEMORY,
1517                       "GCANalloc: couldn't allocate gcs to %d",
1518                       gcan->max_labels) ;
1519           // allocate label storage up to max_labels
1520           gcan->labels =
1521             (unsigned short *)calloc(gcan->max_labels,
1522                                     sizeof(unsigned short)) ;
1523           if (!gcan->labels)
1524             ErrorExit(ERROR_NOMEMORY,
1525                       "GCANalloc: couldn't allocate labels to %d",
1526                       gcan->max_labels) ;
1527         }
1528       }
1529     }
1530   }
1531 
1532   gca->priors = (GCA_PRIOR ***)calloc(gca->prior_width, sizeof(GCA_PRIOR **)) ;
1533   if (!gca->priors)
1534     ErrorExit(ERROR_NOMEMORY, "GCAalloc: could not allocate priors") ;
1535   // setting values to gca->prior volume
1536   for (x = 0 ; x < gca->prior_width ; x++)
1537   {
1538     gca->priors[x] =
1539       (GCA_PRIOR **)calloc(gca->prior_height, sizeof(GCA_PRIOR *)) ;
1540     if (!gca->priors[x])
1541       ErrorExit(ERROR_NOMEMORY, "GCAalloc: could not allocate %dth **",x) ;
1542 
1543     for (y = 0 ; y < gca->prior_height ; y++)
1544     {
1545       gca->priors[x][y] =
1546         (GCA_PRIOR *)calloc(gca->prior_depth, sizeof(GCA_PRIOR)) ;
1547       if (!gca->priors[x][y])
1548         ErrorExit(ERROR_NOMEMORY,
1549                   "GCAalloc: could not allocate %d,%dth *",x,y);
1550       for (z = 0 ; z < gca->prior_depth ; z++)
1551       {
1552         gcap = &gca->priors[x][y][z] ;
1553         if (gcap==NULL)
1554           continue;
1555         gcap->max_labels = max_labels ;
1556 
1557         if (max_labels > 0)
1558         {
1559           /* allocate new ones */
1560           // allocate label space
1561           gcap->labels =
1562             (unsigned short *)calloc(max_labels,
1563                                     sizeof(unsigned short)) ;
1564           if (!gcap->labels)
1565             ErrorExit(ERROR_NOMEMORY, \
1566                       "GCANalloc: couldn't allocate labels to %d",
1567                       gcap->max_labels) ;
1568           // create prior space
1569           gcap->priors = (float *)calloc(max_labels, sizeof(float)) ;
1570           if (!gcap->priors)
1571             ErrorExit(ERROR_NOMEMORY,
1572                       "GCANalloc: couldn't allocate priors to %d",
1573                       max_labels) ;
1574         }
1575       }
1576     }
1577   }
1578   // setup
1579   gca->mri_node__ = 0;
1580   gca->mri_prior__ = 0;
1581   gca->mri_tal__ = 0;
1582   // initialize
1583   gca->node_i_to_r__ = gca->node_r_to_i__ = 0;
1584   gca->prior_i_to_r__ = gca->prior_r_to_i__ = 0;
1585   gca->tal_i_to_r__ = gca->tal_r_to_i__ = 0;
1586 
1587   GCAsetup(gca);
1588 
1589   return(gca) ;
1590 }
1591 
1592 int
1593 GCAfree(GCA **pgca)
1594 {
1595   GCA  *gca ;
1596   int  x, y, z ;
1597 
1598   gca = *pgca ;
1599   *pgca = NULL ;
1600 
1601   for (x = 0 ; x < gca->node_width ; x++)
1602   {
1603     for (y = 0 ; y < gca->node_height ; y++)
1604     {
1605       for (z = 0 ; z < gca->node_depth ; z++)
1606       {
1607         GCANfree(&gca->nodes[x][y][z], gca->ninputs) ;
1608       }
1609       free(gca->nodes[x][y]) ;
1610     }
1611     free(gca->nodes[x]) ;
1612   }
1613   free(gca->nodes) ;
1614 
1615   for (x = 0 ; x < gca->prior_width ; x++)
1616   {
1617     for (y = 0 ; y < gca->prior_height ; y++)
1618     {
1619       for (z = 0 ; z < gca->prior_depth ; z++)
1620       {
1621         free(gca->priors[x][y][z].labels) ;
1622         free(gca->priors[x][y][z].priors) ;
1623       }
1624       free(gca->priors[x][y]) ;
1625     }
1626     free(gca->priors[x]) ;
1627   }
1628 
1629   free(gca->priors) ;
1630   GCAcleanup(gca);
1631 
1632   free(gca) ;
1633 
1634   return(NO_ERROR) ;
1635 }
1636 
1637 int
1638 GCANfree(GCA_NODE *gcan, int ninputs)
1639 {
1640   if (gcan->nlabels)
1641   {
1642     free(gcan->labels) ;
1643     free_gcs(gcan->gcs, gcan->nlabels, ninputs) ;
1644   }
1645   return(NO_ERROR) ;
1646 }
1647 
1648 void PrintInfoOnLabels(GCA *gca, int label, int xn, int yn, int zn,
1649                        int xp, int yp, int zp,
1650                        int x, int y, int z)
1651 {
1652   GCA_NODE  *gcan ;
1653   GCA_PRIOR *gcap ;
1654   GC1D *gc ;
1655   int  i ;
1656   // using node to find the label
1657   gc = GCAfindGC(gca, xn, yn, zn, label) ;
1658   if (gc)
1659   {
1660     gcan = &gca->nodes[xn][yn][zn];
1661     fprintf(stdout,
1662             "\n Node (%3d, %3d, %3d) pos (%3d, %3d, %3d) label=%d, labels:",
1663             xn, yn, zn, x, y, z, label);
1664     for (i=0; i < gcan->nlabels; ++i)
1665       fprintf(stdout, "%4d ", gcan->labels[i]);
1666     fprintf(stdout, "\n");
1667     gcap = &gca->priors[xp][yp][zp];
1668     if (gcap==NULL)
1669       return;
1670     fprintf(stdout,
1671             "Prior (%3d, %3d, %3d) pos (%3d, %3d, %3d) label=%d\n",
1672             xp, yp, zp, x, y, z, label);
1673     fprintf(stdout, "prior label histogram  (label):");
1674     for (i=0; i < gcap->nlabels; ++i)
1675       fprintf(stdout, "%4d ", gcap->labels[i]);
1676     fprintf(stdout, "\n");
1677     fprintf(stdout, "                     :(counts):");
1678     for (i=0; i < gcap->nlabels; ++i)
1679       fprintf(stdout, "%4.f ", gcap->priors[i]);
1680     fprintf(stdout, "\n");
1681     fprintf(stdout, "mean: ");
1682     for (i = 0 ; i < gca->ninputs ; i++)
1683       fprintf(stdout, "%2.1f ", gc->means[i] / gc->ntraining) ;
1684     fprintf(stdout, "\n");
1685   }
1686   fflush(stdout) ;
1687 }
1688 
1689 int
1690 GCAtrainCovariances(GCA *gca,
1691                     MRI *mri_inputs,
1692                     MRI *mri_labels,
1693                     TRANSFORM *transform)
1694 {
1695   int    x, y, z, width, height, depth, label, xn, yn, zn ;
1696   float  vals[MAX_GCA_INPUTS] ;
1697   int xp, yp, zp;
1698 
1699   /* convert transform to voxel coordinates */
1700 
1701   /* go through each voxel in the input volume and find the canonical
1702      voxel (and hence the classifier) to which it maps. Then update the
1703      classifiers statistics based on this voxel's intensity and label.
1704   */
1705   width = mri_labels->width ;
1706   height = mri_labels->height;
1707   depth = mri_labels->depth ;
1708   for (x = 0 ; x < width ; x++)
1709   {
1710     for (y = 0 ; y < height ; y++)
1711     {
1712       for (z = 0 ; z < depth ; z++)
1713       {
1714         if (x == Gx && y == Gy && z == Gz)
1715           DiagBreak() ;
1716         // get the segmented value
1717         label = nint(MRIgetVoxVal(mri_labels, x, y, z,0)) ;
1718 
1719 #if 0
1720         if (!label)
1721           continue ;
1722 #endif
1723         // get all input volume values at this point
1724         load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
1725 
1726         // src -> talairach -> node
1727         if (!GCAsourceVoxelToNode(gca, mri_inputs, transform,
1728                                   x, y, z, &xn, &yn, &zn))
1729         {
1730           if (!GCAsourceVoxelToPrior(gca, mri_inputs,
1731                                      transform, x, y, z,
1732                                      &xp, &yp, &zp))
1733           {
1734             ///////////////// debug code ////////////////////////////
1735             if ((xp == Gxp && yp == Gyp && zp == Gzp) &&
1736                 (Ggca_label < 0 || Ggca_label == label))
1737             {
1738               printf("src (%d, %d, %d), "
1739                      "prior (%d, %d, %d), "
1740                      "node (%d, %d, %d), label = %d\n",
1741                      x,y,z, xp, yp, zp, xn, yn, zn, label);
1742             }
1743             ////////////////////////////////////////////////////////
1744             // update the value
1745             GCAupdateNodeCovariance(gca, mri_inputs,
1746                                     xn, yn, zn, vals, label) ;
1747 
1748             //////////////debug code ////////////////////////////
1749             if (xn == Gxn && yn == Gyn && zn == Gzn)
1750             {
1751               fprintf(stdout, "Train Covariance\n");
1752               PrintInfoOnLabels(gca, label,
1753                                 xn,yn,zn,
1754                                 xp,yp,zp,
1755                                 x,y,z);
1756             }
1757             if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z &&
1758                 (label == Ggca_label || Ggca_label < 0))
1759             {
1760               GC1D *gc ;
1761               int  i, nsamples ;
1762               MATRIX *m ;
1763               gc = GCAfindGC(gca, xn, yn, zn, label) ;
1764               if (gc)
1765               {
1766                 nsamples = gc->ntraining - gc->n_just_priors ;
1767                 /* for no-intensity training */
1768                 if (nsamples < 1)
1769                   nsamples = 1  ;
1770                 printf("voxel(%d,%d,%d) = ", x, y, z) ;
1771                 for (i = 0 ; i < gca->ninputs ; i++)
1772                   printf("%d ", nint(vals[i])) ;
1773 
1774                 printf(" --> node(%d,%d,%d), "
1775                        "label %s (%d), mean ",
1776                        xn, yn, zn,
1777                        cma_label_to_name(label),label) ;
1778                 for (i = 0 ; i < gca->ninputs ; i++)
1779                   printf("%2.1f ", gc->means[i]) ;
1780                 printf("\ncovariances (det=%f):\n",
1781                        covariance_determinant
1782                        (gc,gca->ninputs)/pow
1783                        ((double)nsamples,
1784                         (double)gca->ninputs)) ;
1785                 m =
1786                   load_covariance_matrix(gc,
1787                                          NULL,
1788                                          gca->ninputs) ;
1789                 MatrixScalarMul(m, 1.0/(nsamples), m) ;
1790                 MatrixPrint(stdout, m) ;
1791                 MatrixFree(&m) ;
1792               }
1793             }
1794             /////////////////////////////////////////////////
1795           }
1796         } // if (!GCA...)
1797       }
1798     }
1799   }
1800 
1801   fflush(stdout) ;
1802   return(NO_ERROR) ;
1803 }
1804 #include <unistd.h>
1805 
1806 int
1807 GCAtrain(GCA *gca, MRI *mri_inputs, MRI *mri_labels,
1808          TRANSFORM *transform, GCA *gca_prune,
1809          int noint)
1810 {
1811   int    x, y, z, width, height, depth, label, xn, yn, zn,holes_filled,
1812   /*i, xnbr, ynbr, znbr, xn_nbr, yn_nbr, zn_nbr,*/ xp, yp, zp ;
1813   float  vals[MAX_GCA_INPUTS] ;
1814   static int first_time = 1 ;
1815   FILE *logfp = NULL ;
1816   static int logging = 0 ;
1817   GCA_PRIOR *gcap ;
1818   GCA_NODE  *gcan ;
1819   MRI       *mri_mapped ;
1820 
1821     gca->total_training++ ;
1822   mri_mapped = MRIalloc(gca->prior_width, gca->prior_height, gca->prior_depth,
1823                         MRI_UCHAR) ;
1824   if (first_time)
1825   {
1826     first_time = 0 ;
1827     logging = getenv("GCA_LOG") != NULL ;
1828     if (logging)
1829     {
1830       printf("logging image intensities to GCA log file 'gca*.log'\n") ;
1831       for (label = 0 ; label <= MAX_CMA_LABEL ; label++)
1832       {
1833         char fname[STRLEN] ;
1834         sprintf(fname, "gca%d.log", label) ;
1835         if (FileExists(fname))
1836           unlink(fname) ;
1837       }
1838     }
1839   }
1840 
1841 
1842   /* convert transform to voxel coordinates */
1843 
1844   /* go through each voxel in the input volume and find the canonical
1845      voxel (and hence the classifier) to which it maps. Then update the
1846      classifiers statistics based on this voxel's intensity and label.
1847   */
1848   // segmented volume
1849   width = mri_labels->width ;
1850   height = mri_labels->height;
1851   depth = mri_labels->depth ;
1852   for (x = 0 ; x < width ; x++)
1853   {
1854     for (y = 0 ; y < height ; y++)
1855     {
1856       for (z = 0 ; z < depth ; z++)
1857       {
1858         /// debugging /////////////////////////////////////
1859         if (x == Gx && y == Gy && z == Gz)
1860           DiagBreak() ;
1861         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
1862           DiagBreak() ;
1863 
1864         ///////////////////////////////////////////////////
1865 
1866         // get the segmented voxel label
1867         label = nint(MRIgetVoxVal(mri_labels, x, y, z,0)) ;
1868 #if 0
1869         if (!label)
1870           continue ;
1871 #endif
1872         // get all the values to vals[] in inputs at this point
1873         // mri_inputs are T1, PD etc.
1874         load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
1875         // segmented volume->talairach volume->node
1876         if (!GCAsourceVoxelToNode(gca, mri_inputs, transform,
1877                                   x, y, z, &xn, &yn, &zn))
1878           if (!GCAsourceVoxelToPrior(gca, mri_inputs, transform,
1879                                      x, y, z, &xp, &yp, &zp))
1880           {
1881             // got node point (xn, yn. zn) and
1882             // prior point (xp, yp, zp) for
1883             // this label volume point (x, y, z)
1884             // update the value at this prior point
1885             if (xp == Gxp && yp == Gyp && zp == Gzp)
1886             {
1887               DiagBreak() ;
1888               if (Ggca_label < 0 || Ggca_label == label)
1889               {
1890                 printf("src (%d, %d, %d), "
1891                        "prior (%d, %d, %d), "
1892                        "node (%d, %d, %d), label = %d\n",
1893                        x,y,z, xp, yp, zp, xn, yn, zn, label);
1894               }
1895             }
1896             MRIvox(mri_mapped, xp, yp, zp) = 1 ;
1897             GCAupdatePrior(gca, mri_inputs, xp, yp, zp, label) ;
1898             if ((GCAupdateNode(gca, mri_inputs, xn, yn, zn,
1899                                vals,label,gca_prune, noint) ==
1900                  NO_ERROR) &&
1901                 !(gca->flags & GCA_NO_MRF))
1902               //                 node        label point
1903               GCAupdateNodeGibbsPriors(gca, mri_labels,
1904                                        xn, yn, zn, x, y,z, label);
1905 
1906             /// debugging code //////////////////////////////////////
1907             if (xn == Gxn && yn == Gyn && zn == Gzn)
1908             {
1909               PrintInfoOnLabels(gca, label,
1910                                 xn,yn,zn,
1911                                 xp,yp,zp,
1912                                 x,y,z);
1913             }
1914             if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z &&
1915                 (label == Ggca_label || (Ggca_label < 0))
1916                 && (Ggca_nbr_label < 0))
1917             {
1918               GC1D *gc ;
1919               int  i ;
1920               gc = GCAfindGC(gca, xn, yn, zn, label) ;
1921               if (gc)
1922               {
1923                 if (logging)
1924                 {
1925                   char fname[STRLEN] ;
1926                   sprintf(fname, "gca%d.log", label) ;
1927                   logfp = fopen(fname, "a") ;
1928                 }
1929                 printf("voxel(%d,%d,%d) = ", x, y, z) ;
1930                 for (i = 0 ; i < gca->ninputs ; i++)
1931                 {
1932                   printf("%2.1f ", (vals[i])) ;
1933                   if (logging)
1934                     fprintf(logfp, "%2.1f ", (vals[i])) ;
1935                 }
1936 
1937                 printf(" --> node(%d,%d,%d), "
1938                        "label %s (%d), mean ",
1939                        xn, yn, zn,
1940                        cma_label_to_name(label),label) ;
1941                 for (i = 0 ; i < gca->ninputs ; i++)
1942                   printf("%2.1f ", gc->means[i] / gc->ntraining) ;
1943                 printf("\n") ;
1944                 gcan = &gca->nodes[xn][yn][zn];
1945                 printf("   node labels:");
1946                 for (i=0; i < gcan->nlabels; ++i)
1947                   printf("%d ", gcan->labels[i]);
1948                 printf("\n");
1949                 printf(" --> prior (%d,%d,%d)\n", xp, yp, zp);
1950                 gcap = &gca->priors[xp][yp][zp];
1951                 if (gcap==NULL)
1952                   continue;
1953                 printf("   prior labels:");
1954                 for (i=0; i < gcap->nlabels; ++i)
1955                   printf("%d ", gcap->labels[i]);
1956                 printf("\n");
1957                 if (logging)
1958                 {
1959                   fprintf(logfp, "\n") ;
1960                   fclose(logfp) ;
1961                 }
1962               }
1963               ///////////////////////////////////////////
1964             }
1965           }
1966       }
1967       if (gca->flags & GCA_NO_MRF)
1968         continue ;
1969 #if 0
1970       for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
1971       {
1972         xnbr = x+xnbr_offset[i] ;
1973         ynbr = y+ynbr_offset[i] ;
1974         znbr = z+znbr_offset[i] ;
1975         if (!GCAsourceVoxelToNode(gca, mri_inputs, transform,
1976                                   xnbr, ynbr, znbr,
1977                                   &xn_nbr,&yn_nbr,&zn_nbr))
1978         {
1979           if (xn_nbr == xn && yn_nbr == yn && zn_nbr == zn)
1980             continue ;  /* only update if it is a different node */
1981 
1982           if (GCAupdateNode(gca, mri_inputs, xn_nbr, yn_nbr, zn_nbr,
1983                             vals,label,gca_prune,noint) == NO_ERROR)
1984             GCAupdateNodeGibbsPriors(gca, mri_labels,
1985                                      xn_nbr, yn_nbr, zn_nbr,
1986                                      x, y,z, label);
1987         }
1988       }
1989 #endif
1990     }
1991   }
1992 
1993   for (holes_filled = xp = 0 ; xp < gca->prior_width; xp++)
1994   {
1995     for (yp = 0 ; yp < gca->prior_height; yp++)
1996     {
1997       for (zp = 0 ; zp < gca->prior_depth; zp++)
1998       {
1999         if (MRIvox(mri_mapped, xp, yp, zp) > 0)
2000           continue ;
2001         if (!GCApriorToSourceVoxel(gca, mri_inputs, transform,
2002                                    xp, yp, zp, &x, &y, &z))
2003         {
2004           GC1D *gc ;
2005           label = nint(MRIgetVoxVal(mri_labels, x, y, z, 0)) ;
2006           GCAupdatePrior(gca, mri_inputs, xp, yp, zp, label) ;
2007           GCApriorToNode(gca, xp, yp, zp, &xn, &yn, &zn) ;
2008           gc = GCAfindGC(gca, xn, yn, zn, label) ;
2009           if (gc == NULL)  // label doesn't exist at this position
2010           {
2011             load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
2012             if ((GCAupdateNode(gca, mri_inputs, xn, yn, zn,
2013                                vals,label,gca_prune, noint) ==
2014                  NO_ERROR) && !(gca->flags & GCA_NO_MRF))
2015               GCAupdateNodeGibbsPriors(gca, mri_labels,
2016                                        xn, yn, zn, x, y,z, label);
2017 
2018           }
2019           holes_filled++ ;
2020         }
2021       }
2022     }
2023   }
2024   if (Gdiag & DIAG_WRITE && DIAG_VERBOSE_ON)
2025     MRIwrite(mri_mapped, "m.mgz") ;
2026   if (holes_filled > 0)
2027     printf("%d prior holes filled\n", holes_filled) ;
2028   MRIfree(&mri_mapped) ;
2029   return(NO_ERROR) ;
2030 }
2031 
2032 // declare function pointer
2033 static int (*myclose)(FILE *stream);
2034 
2035 int
2036 GCAwrite(GCA *gca, char *fname)
2037 {
2038   FILE      *fp ;
2039   int       x, y, z, n, i, j ;
2040   GCA_NODE  *gcan ;
2041   GCA_PRIOR *gcap ;
2042   GC1D      *gc ;
2043 
2044   if (strstr(fname, ".gcz"))
2045   {
2046     char command[STRLEN];
2047     myclose = pclose;
2048     strcpy(command, "gzip -f -c > ");
2049     strcat(command, fname);
2050     fp = popen(command, "w");
2051     if (errno)
2052     {
2053       pclose(fp);
2054       errno = 0;
2055       ErrorReturn(ERROR_BADPARM,
2056                   (ERROR_BADPARM,
2057                    "GCAwrite(%s): gzip encountered error",
2058                    fname)) ;
2059     }
2060   }
2061   else
2062   {
2063     myclose = fclose;
2064     fp  = fopen(fname, "wb") ;
2065   }
2066   if (!fp)
2067     ErrorReturn(ERROR_NOFILE,
2068                 (ERROR_NOFILE,
2069                  "GCAwrite: could not open GCA %s for writing",fname)) ;
2070 
2071   fwriteFloat(GCA_INT_VERSION, fp) ;
2072   fwriteFloat(gca->prior_spacing, fp) ;
2073   fwriteFloat(gca->node_spacing, fp) ;
2074   fwriteInt(gca->prior_width,fp);
2075   fwriteInt(gca->prior_height,fp);
2076   fwriteInt(gca->prior_depth,fp);
2077   fwriteInt(gca->node_width,fp);
2078   fwriteInt(gca->node_height,fp);
2079   fwriteInt(gca->node_depth,fp);
2080   fwriteInt(gca->ninputs,fp) ;
2081   fwriteInt(gca->flags, fp) ;
2082 
2083   for (x = 0 ; x < gca->node_width ; x++)
2084   {
2085     for (y = 0 ; y < gca->node_height ; y++)
2086     {
2087       for (z = 0 ; z < gca->node_depth ; z++)
2088       {
2089         if (x == 139 && y == 103 && z == 139)
2090           /* wm should be pallidum */
2091           DiagBreak() ;
2092         gcan = &gca->nodes[x][y][z] ;
2093         fwriteInt(gcan->nlabels, fp) ;
2094         fwriteInt(gcan->total_training, fp) ;
2095         for (n = 0 ; n < gcan->nlabels ; n++)
2096         {
2097           int  r, c ;
2098           gc = &gcan->gcs[n] ;
2099                     fwriteInt(gcan->labels[n], fp) ;
2100           for (r = 0 ; r < gca->ninputs ; r++)
2101             fwriteFloat(gc->means[r], fp) ;
2102           for (r = i = 0 ; r < gca->ninputs ; r++)
2103             for (c = r ;  c < gca->ninputs ; c++, i++)
2104               fwriteFloat(gc->covars[i], fp) ;
2105 
2106           if (gca->flags & GCA_NO_MRF)
2107             continue ;
2108           for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
2109           {
2110             fwriteInt(gc->nlabels[i], fp) ;
2111             for (j = 0 ; j < gc->nlabels[i] ; j++)
2112             {
2113               fwriteInt((int)gc->labels[i][j], fp) ;
2114               fwriteFloat(gc->label_priors[i][j], fp) ;
2115             }
2116           }
2117         }
2118       }
2119     }
2120   }
2121 
2122   for (x = 0 ; x < gca->prior_width ; x++)
2123   {
2124     for (y = 0 ; y < gca->prior_height ; y++)
2125     {
2126       for (z = 0 ; z < gca->prior_depth ; z++)
2127       {
2128         if (x == 139 && y == 103 && z == 139)
2129           /* wm should be pallidum */
2130           DiagBreak() ;
2131         gcap = &gca->priors[x][y][z] ;
2132         if (gcap==NULL)
2133           continue;
2134         fwriteInt(gcap->nlabels, fp) ;
2135         fwriteInt(gcap->total_training, fp) ;
2136         for (n = 0 ; n < gcap->nlabels ; n++)
2137         {
2138                     fwriteInt((int)gcap->labels[n], fp) ;
2139           fwriteFloat(gcap->priors[n], fp) ;
2140         }
2141       }
2142     }
2143   }
2144 
2145   // if (gca->type == GCA_FLASH || gca->type == GCA_PARAM)
2146   // always write gca->type
2147   {
2148     int  n ;
2149 
2150     fwriteInt(FILE_TAG, fp) ;  /* beginning of tagged section */
2151 
2152     /* all tags are format: <int: tag> <int: num> <parm> <parm> .... */
2153     fwriteInt(TAG_GCA_TYPE, fp) ;
2154     fwriteInt(1, fp) ;
2155     fwriteInt(gca->type, fp) ;
2156 
2157     if (gca->type == GCA_FLASH)
2158     {
2159       fwriteInt(TAG_PARAMETERS, fp) ;
2160       fwriteInt(3, fp) ;   /* currently only storing 3 parameters */
2161       for (n = 0 ; n < gca->ninputs ; n++)
2162       {
2163         fwriteFloat(gca->TRs[n], fp) ;
2164         fwriteFloat(gca->FAs[n], fp) ;
2165         fwriteFloat(gca->TEs[n], fp) ;
2166       }
2167     }
2168   }
2169 
2170   // write direction cosine information
2171   fwriteInt(TAG_GCA_DIRCOS, fp);
2172   fwriteFloat(gca->x_r, fp);
2173   fwriteFloat(gca->x_a, fp);
2174   fwriteFloat(gca->x_s, fp);
2175   fwriteFloat(gca->y_r, fp);
2176   fwriteFloat(gca->y_a, fp);
2177   fwriteFloat(gca->y_s, fp);
2178   fwriteFloat(gca->z_r, fp);
2179   fwriteFloat(gca->z_a, fp);
2180   fwriteFloat(gca->z_s, fp);
2181   fwriteFloat(gca->c_r, fp);
2182   fwriteFloat(gca->c_a, fp);
2183   fwriteFloat(gca->c_s, fp);
2184   fwriteInt(gca->width, fp);
2185   fwriteInt(gca->height, fp);
2186   fwriteInt(gca->depth, fp);
2187   fwriteFloat(gca->xsize, fp);
2188   fwriteFloat(gca->ysize, fp);
2189   fwriteFloat(gca->zsize, fp);
2190 
2191   // fclose(fp) ;
2192   myclose(fp);
2193 
2194   return(NO_ERROR) ;
2195 }
2196 
2197 GCA *
2198 GCAread(char *fname)
2199 {
2200   FILE      *fp ;
2201   int       x, y, z, n, i, j ;
2202   GCA       *gca ;
2203   GCA_NODE  *gcan ;
2204   GCA_PRIOR *gcap ;
2205   GC1D      *gc ;
2206   float     version, node_spacing, prior_spacing ;
2207   int       node_width, node_height, node_depth,
2208   prior_width, prior_height, prior_depth,
2209   ninputs, flags ;
2210   int       tag;
2211 
2212   if (strstr(fname, ".gcz"))
2213   {
2214     char command[STRLEN];
2215     myclose = pclose;
2216 #if defined(Darwin) || defined(SunOS)
2217     // zcat on Max OS always appends and assumes a .Z extention,
2218     // whereas we want .gcz
2219     strcpy(command, "gunzip -c ");
2220 #else
2221     strcpy(command, "zcat ");
2222 #endif
2223     strcat(command, fname);
2224     errno = 0;
2225     fp = popen(command, "r");
2226     if (errno)
2227     {
2228       pclose(fp);
2229       errno = 0;
2230       ErrorReturn(NULL, (ERROR_BADPARM,
2231                          "GCAread: encountered error executing: '%s'",
2232                          command)) ;
2233     }
2234   }
2235   else
2236   {
2237     myclose = fclose;
2238     fp  = fopen(fname, "rb") ;
2239   }
2240   if (!fp)
2241     ErrorReturn(NULL,
2242                 (ERROR_NOFILE,
2243                  "GCAread: could not open GCA %s for reading",fname)) ;
2244 
2245   if (!freadFloatEx(&version, fp))
2246     ErrorReturn(NULL, (ERROR_BADPARM,"GCAread(%s): could not read file",
2247                        fname)) ;
2248 
2249   if (version < GCA_UCHAR_VERSION)
2250   {
2251     node_spacing = freadFloat(fp) ;
2252     node_width = freadInt(fp);
2253     node_height = freadInt(fp);
2254     node_depth = freadInt(fp);
2255     ninputs = freadInt(fp) ;
2256     if (version == 3.0)
2257       flags = freadInt(fp) ;
2258     else
2259       flags = 0 ;
2260 
2261     gca = gcaAllocMax(ninputs, node_spacing, \
2262                       node_spacing, node_spacing*node_width,
2263                       node_spacing*node_height, \
2264                       node_spacing*node_depth, 0, flags) ;
2265     if (!gca)
2266       ErrorReturn(NULL, (Gerror, NULL)) ;
2267 
2268     for (x = 0 ; x < gca->node_width ; x++)
2269     {
2270       for (y = 0 ; y < gca->node_height ; y++)
2271       {
2272         for (z = 0 ; z < gca->node_depth ; z++)
2273         {
2274           if (x == 28 && y == 39 && z == 39)
2275             DiagBreak() ;
2276           gcan = &gca->nodes[x][y][z] ;
2277           gcan->nlabels = freadInt(fp) ;
2278           gcan->total_training = freadInt(fp) ;
2279           if (gcan->nlabels)
2280           {
2281             gcan->labels = \
2282                            (unsigned short *)calloc(gcan->nlabels,
2283                                                    sizeof(unsigned short)) ;
2284             if (!gcan->labels)
2285               ErrorExit(ERROR_NOMEMORY,
2286                         "GCAread(%s): could not allocate %d "
2287                         "labels @ (%d,%d,%d)",
2288                         fname, gcan->nlabels, x, y, z) ;
2289             gcan->gcs = alloc_gcs(gcan->nlabels,
2290                                   flags, gca->ninputs) ;
2291             if (!gcan->gcs)
2292               ErrorExit(ERROR_NOMEMORY,
2293                         "GCAread(%s); could not allocated %d gcs "
2294                         "@ (%d,%d,%d)",
2295                         fname, gcan->nlabels, x, y, z) ;
2296           }
2297           else // no labels assigned to this node
2298           {
2299             gcan->labels = 0;
2300             gcan->gcs = 0;
2301           }
2302           for (n = 0 ; n < gcan->nlabels ; n++)
2303           {
2304             int r, c;
2305             gc = &gcan->gcs[n] ;
2306             gcan->labels[n] = (unsigned short)fgetc(fp) ;
2307             for (r = 0 ; r < gca->ninputs ; r++)
2308               gc->means[r] = freadFloat(fp) ;
2309             for (i = r = 0 ; r < gca->ninputs ; r++)
2310               for (c = r ; c < gca->ninputs ; c++, i++)
2311                 gc->covars[i] = freadFloat(fp) ;
2312             if (gca->flags & GCA_NO_MRF)
2313               continue ;
2314             for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
2315             {
2316               gc->nlabels[i] = freadInt(fp) ;
2317 
2318 #if 1
2319               /* allocate new ones */
2320               gc->label_priors[i] =
2321                 (float *)calloc(gc->nlabels[i],sizeof(float));
2322               if (!gc->label_priors[i])
2323                 ErrorExit(ERROR_NOMEMORY, "GCAread(%s): "
2324                           "couldn't expand gcs to %d",
2325                           fname,gc->nlabels) ;
2326               gc->labels[i] =
2327                 (unsigned short *)calloc(gc->nlabels[i], \
2328                                         sizeof(unsigned short)) ;
2329               if (!gc->labels)
2330                 ErrorExit(ERROR_NOMEMORY,
2331                           "GCAread(%s): couldn't expand "
2332                           "labels to %d",
2333                           fname, gc->nlabels[i]) ;
2334 #endif
2335               for (j = 0 ; j < gc->nlabels[i] ; j++)
2336               {
2337                 gc->labels[i][j] = (unsigned short)freadInt(fp) ;
2338                 gc->label_priors[i][j] = freadFloat(fp) ;
2339               }
2340             }
2341           }
2342         }
2343       }
2344     }
2345   }
2346   else   /* current version - stores priors at different
2347                             resolution than densities */
2348   {
2349     if (!FEQUAL(version, GCA_UCHAR_VERSION) &&
2350                 !FEQUAL(version, GCA_INT_VERSION))
2351     {
2352       // fclose(fp) ;
2353       myclose(fp);
2354       ErrorReturn(NULL, (ERROR_BADFILE,
2355                          "GCAread(%s), version #%2.1f found, "
2356                          "%2.1f expected",
2357                          fname, version, GCA_INT_VERSION)) ;
2358     }
2359     prior_spacing = freadFloat(fp) ;
2360     node_spacing = freadFloat(fp) ;
2361     prior_width = freadInt(fp);
2362     prior_height = freadInt(fp);
2363     prior_depth = freadInt(fp);
2364     node_width = freadInt(fp);
2365     node_height = freadInt(fp);
2366     node_depth = freadInt(fp);
2367     ninputs = freadInt(fp) ;
2368     flags = freadInt(fp) ;
2369 
2370     gca = gcaAllocMax(ninputs, prior_spacing, node_spacing,
2371                       node_spacing*node_width,
2372                       node_spacing*node_height,
2373                       node_spacing*node_depth, 0, flags) ;
2374     if (!gca)
2375       ErrorReturn(NULL, (Gdiag, NULL)) ;
2376 
2377     for (x = 0 ; x < gca->node_width ; x++)
2378     {
2379       for (y = 0 ; y < gca->node_height ; y++)
2380       {
2381         for (z = 0 ; z < gca->node_depth ; z++)
2382         {
2383           if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
2384             DiagBreak() ;
2385           gcan = &gca->nodes[x][y][z] ;
2386           gcan->nlabels = freadInt(fp) ;
2387           gcan->total_training = freadInt(fp) ;
2388           if (gcan->nlabels)
2389           {
2390             gcan->labels =
2391               (unsigned short *)calloc(gcan->nlabels,
2392                                       sizeof(unsigned short)) ;
2393             if (!gcan->labels)
2394               ErrorExit(ERROR_NOMEMORY, "GCAread(%s): could not "
2395                         "allocate %d "
2396                         "labels @ (%d,%d,%d)",
2397                         fname, gcan->nlabels, x, y, z) ;
2398             gcan->gcs = alloc_gcs(gcan->nlabels,
2399                                   flags,
2400                                   gca->ninputs) ;
2401             if (!gcan->gcs)
2402               ErrorExit(ERROR_NOMEMORY,
2403                         "GCAread(%s); could not allocated %d gcs "
2404                         "@ (%d,%d,%d)",
2405                         fname, gcan->nlabels, x, y, z) ;
2406           }
2407           else // no labels at this node
2408           {
2409             gcan->labels = 0;
2410             gcan->gcs = 0;
2411           }
2412           for (n = 0 ; n < gcan->nlabels ; n++)
2413           {
2414             int  r, c ;
2415 
2416             gc = &gcan->gcs[n] ;
2417 
2418                         if (version == GCA_UCHAR_VERSION)
2419                             gcan->labels[n] = (unsigned short)fgetc(fp) ;
2420                         else
2421                             gcan->labels[n] = (unsigned short)freadInt(fp) ;
2422 
2423             for (r = 0 ; r < gca->ninputs ; r++)
2424               gc->means[r] = freadFloat(fp) ;
2425             for (i = r = 0 ; r < gca->ninputs ; r++)
2426               for (c = r ; c < gca->ninputs ; c++, i++)
2427                 gc->covars[i] = freadFloat(fp) ;
2428             if (gca->flags & GCA_NO_MRF)
2429               continue ;
2430             for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
2431             {
2432               gc->nlabels[i] = freadInt(fp) ;
2433 
2434               /* allocate new ones */
2435               gc->label_priors[i] =
2436                 (float *)calloc(gc->nlabels[i],sizeof(float));
2437               if (!gc->label_priors[i])
2438                 ErrorExit(ERROR_NOMEMORY, "GCAread(%s): "
2439                           "couldn't expand gcs to %d",
2440                           fname,gc->nlabels) ;
2441               gc->labels[i] =
2442                 (unsigned short *)calloc(gc->nlabels[i], \
2443                                         sizeof(unsigned short)) ;
2444               if (!gc->labels)
2445                 ErrorExit(ERROR_NOMEMORY,
2446                           "GCAread(%s): couldn't expand "
2447                           "labels to %d",
2448                           fname, gc->nlabels[i]) ;
2449               for (j = 0 ; j < gc->nlabels[i] ; j++)
2450               {
2451                 gc->labels[i][j] = (unsigned short)freadInt(fp) ;
2452                 gc->label_priors[i][j] = freadFloat(fp) ;
2453               }
2454             }
2455           }
2456         }
2457       }
2458     }
2459 
2460     for (x = 0 ; x < gca->prior_width ; x++)
2461     {
2462       for (y = 0 ; y < gca->prior_height ; y++)
2463       {
2464         for (z = 0 ; z < gca->prior_depth ; z++)
2465         {
2466           if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
2467             DiagBreak() ;
2468           gcap = &gca->priors[x][y][z] ;
2469           if (gcap==NULL)
2470             continue;
2471           gcap->nlabels = freadInt(fp) ;
2472           gcap->total_training = freadInt(fp) ;
2473           if (gcap->nlabels)
2474           {
2475             gcap->labels =
2476               (unsigned short *)calloc(gcap->nlabels, \
2477                                       sizeof(unsigned short)) ;
2478             if (!gcap->labels)
2479               ErrorExit(ERROR_NOMEMORY, "GCAread(%s): could not "
2480                         "allocate %d "
2481                         "labels @ (%d,%d,%d)",
2482                         fname, gcap->nlabels, x, y, z) ;
2483             gcap->priors = (float *)calloc(gcap->nlabels,
2484                                            sizeof(float)) ;
2485             if (!gcap->priors)
2486               ErrorExit(ERROR_NOMEMORY, "GCAread(%s): could "
2487                         "not allocate %d "
2488                         "priors @ (%d,%d,%d)",
2489                         fname, gcap->nlabels, x, y, z) ;
2490           }
2491           else // no labels assigned to this priors
2492           {
2493             gcap->labels = 0;
2494             gcap->priors = 0;
2495           }
2496           for (n = 0 ; n < gcap->nlabels ; n++)
2497           {
2498                         if (version == GCA_UCHAR_VERSION)
2499                             gcap->labels[n] = (unsigned short)fgetc(fp) ;
2500                         else
2501                             gcap->labels[n] = (unsigned short)freadInt(fp) ;
2502             gcap->priors[n] = freadFloat(fp) ;
2503           }
2504         }
2505       }
2506     }
2507   }
2508 
2509   for (x = 0 ; x < gca->node_width ; x++)
2510   {
2511     for (y = 0 ; y < gca->node_height ; y++)
2512     {
2513       for (z = 0 ; z < gca->node_depth ; z++)
2514       {
2515         int xp, yp, zp ;
2516 
2517         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
2518           DiagBreak() ;
2519         gcan = &gca->nodes[x][y][z] ;
2520         if (gcaNodeToPrior(gca, x, y, z, &xp, &yp, &zp)==NO_ERROR)
2521         {
2522           gcap = &gca->priors[xp][yp][zp] ;
2523           if (gcap==NULL)
2524             continue;
2525           for (n = 0 ; n < gcan->nlabels ; n++)
2526           {
2527             gc = &gcan->gcs[n] ;
2528             gc->ntraining =
2529               gcan->total_training * getPrior(gcap,gcan->labels[n]) ;
2530           }
2531         }
2532       }
2533     }
2534   }
2535 
2536   // if (!feof(fp))  // this does not work ;-)
2537   // feof(fp) check does not work, since feof is not signaled until you read
2538   while (freadIntEx(&tag, fp))
2539   {
2540     int  n, nparms ;
2541 
2542     // tag = freadInt(fp) ;
2543     if (tag == FILE_TAG) /* beginning of tagged section */
2544     {
2545       // while (!feof(fp))
2546       while (freadIntEx(&tag, fp))
2547       {
2548         /* all tags are format:
2549            <int: tag> <int: num> <parm> <parm> .... */
2550         // tag = freadInt(fp) ;
2551         switch (tag)
2552         {
2553         case TAG_GCA_TYPE:
2554           freadInt(fp) ;   /* skip num=1 */
2555           gca->type = freadInt(fp) ;
2556           if (DIAG_VERBOSE_ON) switch (gca->type)
2557             {
2558             case GCA_NORMAL:
2559               printf("setting gca type = Normal gca type\n");
2560               break;
2561             case GCA_PARAM:
2562               printf("setting gca type = T1/PD gca type\n");
2563               break;
2564             case GCA_FLASH:
2565               printf("setting gca type = FLASH gca type\n");
2566               break;
2567             default:
2568               printf("setting gca type = Unknown\n");
2569               gca->type=GCA_UNKNOWN;
2570               break;
2571             }
2572           break ;
2573         case TAG_PARAMETERS:
2574           nparms = freadInt(fp) ;
2575           /* how many MR parameters are stored */
2576           printf("reading %d MR parameters out of GCA header...\n",
2577                  nparms) ;
2578           for (n = 0 ; n < gca->ninputs ; n++)
2579           {
2580             gca->TRs[n] = freadFloat(fp) ;
2581             gca->FAs[n] = freadFloat(fp) ;
2582             gca->TEs[n] = freadFloat(fp) ;
2583             printf("input %d: TR=%2.1f msec, FA=%2.1f deg, "
2584                    "TE=%2.1f msec\n",
2585                    n,
2586                    gca->TRs[n],
2587                    DEGREES(gca->FAs[n]), gca->TEs[n]) ;
2588           }
2589           break ;
2590         case TAG_GCA_DIRCOS:
2591           gca->x_r = freadFloat(fp);
2592           gca->x_a = freadFloat(fp);
2593           gca->x_s = freadFloat(fp);
2594           gca->y_r = freadFloat(fp);
2595           gca->y_a = freadFloat(fp);
2596           gca->y_s = freadFloat(fp);
2597           gca->z_r = freadFloat(fp);
2598           gca->z_a = freadFloat(fp);
2599           gca->z_s = freadFloat(fp);
2600           gca->c_r = freadFloat(fp);
2601           gca->c_a = freadFloat(fp);
2602           gca->c_s = freadFloat(fp);
2603           gca->width = freadInt(fp);
2604           gca->height = freadInt(fp);
2605           gca->depth = freadInt(fp);
2606           gca->xsize = freadFloat(fp);
2607           gca->ysize = freadFloat(fp);
2608           gca->zsize = freadFloat(fp);
2609 
2610           if (Gdiag & DIAG_SHOW && DIAG_VERBOSE_ON)
2611           {
2612             printf("Direction cosines read:\n");
2613             printf(" x_r = % .4f, y_r = % .4f, z_r = % .4f\n",
2614                    gca->x_r, gca->y_r, gca->z_r);
2615             printf(" x_a = % .4f, y_a = % .4f, z_a = % .4f\n",
2616                    gca->x_a, gca->y_a, gca->z_a);
2617             printf(" x_s = % .4f, y_s = % .4f, z_s = % .4f\n",
2618                    gca->x_s, gca->y_s, gca->z_s);
2619             printf(" c_r = % .4f, c_a = % .4f, c_s = % .4f\n",
2620                    gca->c_r, gca->c_a, gca->c_s);
2621           }
2622           break;
2623         default:
2624           ErrorPrintf(ERROR_BADFILE,
2625                       "GCAread(%s): unknown tag %x\n", fname, tag) ;
2626           break ;
2627         }
2628       }
2629     }
2630   }
2631 
2632   GCAsetup(gca);
2633 
2634   // fclose(fp) ;
2635   myclose(fp);
2636 
2637   return(gca) ;
2638 }
2639 
2640 static int
2641 GCAupdatePrior(GCA *gca, MRI *mri, int xn, int yn, int zn, int label)
2642 {
2643   int       n ;
2644   GCA_PRIOR *gcap ;
2645 
2646   if (label >= MAX_CMA_LABEL)
2647     ErrorReturn(ERROR_BADPARM,
2648                 (ERROR_BADPARM,
2649                  "GCAupdatePrior(%d, %d, %d, %d): label out of range",
2650                  xn, yn, zn, label)) ;
2651 
2652   if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z)
2653     DiagBreak() ;
2654 
2655   gcap = &gca->priors[xn][yn][zn] ;
2656   if (gcap==NULL)
2657     return -1;
2658   // find the label histogram index n
2659   for (n = 0 ; n < gcap->nlabels ; n++)
2660   {
2661     if (gcap->labels[n] == label)
2662       break ;
2663   }
2664   // if index is beyond what we have, then
2665   if (n >= gcap->nlabels)  /* have to allocate a new classifier */
2666   {
2667 
2668     if (n >= gcap->max_labels)
2669     {
2670       int  old_max_labels ;
2671       unsigned short *old_labels ;
2672       float *old_priors ;
2673 
2674       old_max_labels = gcap->max_labels ;
2675       gcap->max_labels += 2 ;
2676       old_labels = gcap->labels ;
2677       old_priors = gcap->priors ;
2678 
2679       /* allocate new ones */
2680       gcap->priors = (float *)calloc(gcap->max_labels, sizeof(float)) ;
2681 
2682       if (!gcap->priors)
2683         ErrorExit(ERROR_NOMEMORY,
2684                   "GCANupdatePriors: couldn't expand priors to %d",
2685                   gcap->max_labels) ;
2686       gcap->labels =
2687         (unsigned short *)calloc(gcap->max_labels, sizeof(unsigned short)) ;
2688       if (!gcap->labels)
2689         ErrorExit(ERROR_NOMEMORY,
2690                   "GCANupdatePriors: couldn't expand labels to %d",
2691                   gcap->max_labels) ;
2692 
2693       /* copy the old ones over */
2694       memmove(gcap->priors, old_priors, old_max_labels*sizeof(float)) ;
2695       memmove(gcap->labels, old_labels,
2696               old_max_labels*sizeof(unsigned short)) ;
2697 
2698       /* free the old ones */
2699       free(old_priors) ;
2700       free(old_labels) ;
2701     }
2702     // add one
2703     gcap->nlabels++ ;
2704   }
2705 
2706   /* these will be updated when training is complete */
2707   // give the value at n
2708   gcap->priors[n] += 1.0f ; // increment counter
2709   gcap->total_training++ ;  // increment training counter
2710   gcap->labels[n] = label ; // save the label value
2711 
2712   return(NO_ERROR) ;
2713 }
2714 
2715 static int
2716 GCAupdateNodeCovariance(GCA *gca, MRI *mri,
2717                         int xn, int yn, int zn, float *vals,int label)
2718 {
2719   int      n, r, c, v ;
2720   GCA_NODE *gcan ;
2721   GC1D     *gc ;
2722 
2723 
2724   if (label >= MAX_CMA_LABEL)
2725     ErrorReturn(ERROR_BADPARM,
2726                 (ERROR_BADPARM,
2727                  "GCAupdateNodeCovariance(%d, %d, %d, %d): label out of range",
2728                  xn, yn, zn, label)) ;
2729 
2730   ///////////// debug ////////////////////////////////
2731   if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z && label == Ggca_label)
2732     DiagBreak() ;
2733   if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z)
2734     DiagBreak() ;
2735   ////////////////////////////////////////////////////
2736 
2737   gcan = &gca->nodes[xn][yn][zn] ;
2738 
2739   for (n = 0 ; n < gcan->nlabels ; n++)
2740   {
2741     if (gcan->labels[n] == label)
2742       break ;
2743   }
2744   if (n >= gcan->nlabels)
2745     ErrorExit(ERROR_BADPARM,
2746               "GCAupdateNodeCovariance(%d, %d, %d, %d): could not find label",
2747               xn, yn, zn, label) ;
2748 
2749   gc = &gcan->gcs[n] ;
2750 
2751   /* remove means from vals */
2752   for (r = 0 ; r < gca->ninputs ; r++)
2753     vals[r] -= gc->means[r] ;
2754 
2755   /* these will be normalized when training is complete */
2756   for (v = r = 0 ; r < gca->ninputs ; r++)
2757   {
2758     for (c = r ; c < gca->ninputs ; c++, v++)
2759       gc->covars[v] += vals[r]*vals[c] ;
2760   }
2761 
2762   /* put means back in vals so caller isn't confused */
2763   for (r = 0 ; r < gca->ninputs ; r++)
2764     vals[r] += gc->means[r] ;
2765 
2766   return(NO_ERROR) ;
2767 }
2768 
2769 static int
2770 GCAupdateNode(GCA *gca, MRI *mri,
2771               int xn, int yn, int zn, float *vals, int label,
2772               GCA *gca_prune, int noint)
2773 {
2774   int      n, i ;
2775   GCA_NODE *gcan ;
2776   GC1D     *gc ;
2777 
2778 
2779   if (label >= MAX_CMA_LABEL)
2780     ErrorReturn(ERROR_BADPARM,
2781                 (ERROR_BADPARM,
2782                  "GCAupdateNode(%d, %d, %d, %d): label out of range",
2783                  xn, yn, zn, label)) ;
2784 
2785   if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z && label == Ggca_label)
2786     DiagBreak() ;
2787   if (xn == Ggca_x && yn == Ggca_y && zn == Ggca_z)
2788     DiagBreak() ;
2789 
2790   // if non-zero label and gca_prune is there
2791   if (label > 0 && gca_prune != NULL && !noint)
2792   {
2793     GCA_NODE *gcan_prune ;
2794     GC1D     *gc_prune ;
2795     int      nprune ;
2796 
2797     gcan_prune = &gca_prune->nodes[xn][yn][zn] ;
2798 
2799     for (nprune = 0 ; nprune < gcan_prune->nlabels ; nprune++)
2800     {
2801       if (gcan_prune->labels[nprune] == label)
2802         break ;
2803     }
2804     if (nprune >= gcan_prune->nlabels)
2805       ErrorPrintf(ERROR_BADPARM,
2806                   "WARNING: pruning GCA at (%d,%d,%d) doesn't "
2807                   "contain label %d", xn, yn, zn, label) ;
2808     gc_prune = &gcan_prune->gcs[nprune] ;
2809 
2810     if (sqrt(GCAmahDist(gc_prune, vals, gca->ninputs)) > 2)
2811       /* more than 2 stds from mean */
2812     {
2813 #if 0
2814       if (xn == 23 && yn == 30 && zn == 32)
2815       {
2816         printf(
2817           "pruning val %2.0f, label %d @ "
2818           "(%d,%d,%d),u=%2.1f, std=%2.1f\n"
2819           ,val, label, xn,yn,zn,mean, std) ;
2820         DiagBreak() ;
2821       }
2822 #endif
2823       total_pruned++ ;
2824       return(ERROR_BAD_PARM) ;
2825     }
2826   }
2827 
2828   // get one at this point
2829   gcan = &gca->nodes[xn][yn][zn] ;
2830 
2831   // look for this label in the array
2832   for (n = 0 ; n < gcan->nlabels ; n++)
2833   {
2834     if (gcan->labels[n] == label)
2835       break ;
2836   }
2837   // if not found
2838   if (n >= gcan->nlabels)  /* have to allocate a new classifier */
2839   {
2840 
2841     if (n >= gcan->max_labels)
2842     {
2843       int  old_max_labels ;
2844       unsigned short *old_labels ;
2845       GC1D *old_gcs ;
2846 
2847       old_max_labels = gcan->max_labels ;
2848       gcan->max_labels += 2 ;
2849       old_labels = gcan->labels ;
2850       old_gcs = gcan->gcs ;
2851 
2852       /* allocate new ones */
2853 #if 0
2854       gcan->gcs = (GC1D *)calloc(gcan->max_labels, sizeof(GC1D)) ;
2855 #else
2856       gcan->gcs = alloc_gcs(gcan->max_labels, gca->flags, gca->ninputs) ;
2857 #endif
2858 
2859       if (!gcan->gcs)
2860         ErrorExit(ERROR_NOMEMORY,
2861                   "GCANupdateNode: couldn't expand gcs to %d",
2862                   gcan->max_labels) ;
2863       gcan->labels =
2864         (unsigned short *)calloc(gcan->max_labels, sizeof(unsigned short)) ;
2865       if (!gcan->labels)
2866         ErrorExit(ERROR_NOMEMORY,
2867                   "GCANupdateNode: couldn't expand labels to %d",
2868                   gcan->max_labels) ;
2869 
2870       /* copy the old ones over */
2871 #if 0
2872       memmove(gcan->gcs, old_gcs, old_max_labels*sizeof(GC1D)) ;
2873 #else
2874       copy_gcs(old_max_labels, old_gcs, gcan->gcs, gca->ninputs) ;
2875 #endif
2876       memmove(gcan->labels, old_labels,
2877               old_max_labels*sizeof(unsigned short)) ;
2878 
2879       /* free the old ones */
2880       free(old_gcs) ;
2881       free(old_labels) ;
2882     }
2883     gcan->nlabels++ ;
2884   }
2885 
2886   gc = &gcan->gcs[n] ;
2887 
2888   /* these will be updated when training is complete */
2889   if (noint)
2890     gc->n_just_priors++ ;
2891   else
2892   {
2893     // vals[] array is the values of inputs at this point
2894     for (i = 0 ; i < gca->ninputs ; i++)
2895       gc->means[i] += vals[i] ;
2896     // get the mean valu (note it is not divided by ninputs!)
2897     /*gc->var += val*val ; */
2898   }
2899   if (gc->n_just_priors >= gc->ntraining)
2900     DiagBreak() ;
2901   gcan->total_training++ ;
2902 
2903   gcan->labels[n] = label ;
2904   gc->ntraining++ ;
2905 
2906   return(NO_ERROR) ;
2907 }
2908 int
2909 GCAcompleteMeanTraining(GCA *gca)
2910 {
2911   int       x, y, z, n, total_nodes, \
2912   total_gcs, i, j, holes_filled, total_brain_gcs,\
2913   total_brain_nodes, r ;
2914   float     nsamples ;
2915   GCA_NODE  *gcan ;
2916   GCA_PRIOR *gcap ;
2917   GC1D      *gc ;
2918 
2919   total_nodes = gca->node_width*gca->node_height*gca->node_depth ;
2920   total_brain_nodes = total_gcs = total_brain_gcs = 0 ;
2921   for (x = 0 ; x < gca->node_width ; x++)
2922   {
2923     for (y = 0 ; y < gca->node_height ; y++)
2924     {
2925       for (z = 0 ; z < gca->node_depth ; z++)
2926       {
2927         gcan = &gca->nodes[x][y][z] ;
2928         total_gcs += gcan->nlabels ;
2929         if (gcan->nlabels > 1 || !IS_UNKNOWN(gcan->labels[0]))
2930         {
2931           total_brain_gcs += gcan->nlabels ;
2932           total_brain_nodes++ ;
2933         }
2934 
2935         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
2936           DiagBreak() ;
2937         for (n = 0 ; n < gcan->nlabels ; n++)
2938         {
2939           gc = &gcan->gcs[n] ;
2940           nsamples = gc->ntraining ;
2941           if ((gca->flags & GCA_NO_MRF) == 0)
2942           {
2943             for (i = 0 ; i < GIBBS_NEIGHBORS ; i++)
2944             {
2945               for (j = 0 ; j < gc->nlabels[i] ; j++)
2946               {
2947                 gc->label_priors[i][j] /= (float)nsamples ;
2948                 check_finite(
2949                   "GCAcompleteMeanTraining: "
2950                   "label_priors",
2951                   gc->label_priors[i][j]) ;
2952               }
2953             }
2954           }
2955 
2956           nsamples -= gc->n_just_priors ;
2957           /* for no-intensity training */
2958           if (nsamples > 0)
2959           {
2960             for (r = 0 ; r < gca->ninputs ; r++)
2961             {
2962               gc->means[r] /= nsamples ;
2963               check_finite("GCAcompleteMeanTraining: mean",
2964                            gc->means[r]) ;
2965             }
2966           }
2967           else
2968           {
2969             int r ;
2970             for (r = 0 ; r < gca->ninputs ; r++)
2971               gc->means[r] = -1 ;  /* mark it for later processing */
2972           }
2973         }
2974       }
2975     }
2976   }
2977 
2978   for (x = 0 ; x < gca->prior_width ; x++)
2979   {
2980     for (y = 0 ; y < gca->prior_height ; y++)
2981     {
2982       for (z = 0 ; z < gca->prior_depth ; z++)
2983       {
2984         gcap = &gca->priors[x][y][z] ;
2985         if (gcap==NULL)
2986           continue;
2987         for (n = 0 ; n < gcap->nlabels ; n++)
2988         {
2989           gcap->priors[n] /= (float)gcap->total_training ;
2990           check_finite("GCAcompleteMeanTraining: priors",
2991                        gcap->priors[n]) ;
2992         }
2993 
2994       }
2995     }
2996   }
2997   printf("filling holes in the GCA...\n") ;
2998   for (holes_filled = x = 0 ; x < gca->node_width ; x++)
2999   {
3000     for (y = 0 ; y < gca->node_height ; y++)
3001     {
3002       for (z = 0 ; z < gca->node_depth ; z++)
3003       {
3004         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3005           DiagBreak() ;
3006         gcan = &gca->nodes[x][y][z] ;
3007         for (n = 0 ; n < gcan->nlabels ; n++)
3008         {
3009           gc = &gcan->gcs[n] ;
3010           nsamples = gc->ntraining - gc->n_just_priors ;
3011           if (nsamples <= 0)
3012           {
3013             GC1D *gc_nbr ;
3014             int  r, i ;
3015 
3016             gc_nbr = findClosestValidGC(gca, x, y, z,
3017                                         gcan->labels[n], 0) ;
3018             if (!gc_nbr)
3019             {
3020               ErrorPrintf(ERROR_BADPARM,
3021                           "gca(%d,%d,%d,%d) - could not "
3022                           "find valid nbr label "
3023                           "%s (%d)", x,  y, z, n,
3024                           cma_label_to_name(gcan->labels[n]),
3025                           gcan->labels[n]) ;
3026               continue ;
3027             }
3028             holes_filled++ ;
3029             for (i = r = 0 ; r < gca->ninputs ; r++)
3030             {
3031               gc->means[r] = gc_nbr->means[r] ;
3032             }
3033             if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3034               printf("filling hole @ (%d, %d, %d)\n", x, y, z) ;
3035           }
3036         }
3037       }
3038     }
3039   }
3040 
3041   printf("%d classifiers: %2.1f per node, %2.2f in brain (%d holes filled)\n",
3042          total_gcs, (float)total_gcs/(float)total_nodes,
3043          (float)total_brain_gcs/(float)total_brain_nodes,holes_filled) ;
3044   if (total_pruned > 0)
3045   {
3046     printf("%d samples pruned during training\n", total_pruned) ;
3047     total_pruned = 0 ;
3048   }
3049   return(NO_ERROR) ;
3050 }
3051 
3052 int
3053 GCAcompleteCovarianceTraining(GCA *gca)
3054 {
3055   int       x, y, z, n, holes_filled, r, c,v, nregs = 0, nsamples ;
3056   GCA_NODE  *gcan ;
3057   GC1D      *gc ;
3058 
3059   for (x = 0 ; x < gca->node_width ; x++)
3060   {
3061     for (y = 0 ; y < gca->node_height ; y++)
3062     {
3063       for (z = 0 ; z < gca->node_depth ; z++)
3064       {
3065         gcan = &gca->nodes[x][y][z] ;
3066 
3067         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3068           DiagBreak() ;
3069         for (n = 0 ; n < gcan->nlabels ; n++)
3070         {
3071           if (x == Ggca_x && y == Ggca_y && z == Ggca_z &&
3072               (Ggca_label == gcan->labels[n] || Ggca_label < 0))
3073             DiagBreak() ;
3074           gc = &gcan->gcs[n] ;
3075           nsamples = gc->ntraining - gc->n_just_priors ;
3076           /* for no-intensity training */
3077           if (nsamples > 0)
3078           {
3079             for (r = v = 0 ; r < gca->ninputs ; r++)
3080             {
3081               check_finite("GCAcompleteCovarianceTraining: mean",
3082                            gc->means[r]) ;
3083               if (nsamples > 1)
3084               {
3085                 for (c = r ; c < gca->ninputs ; c++, v++)
3086                 {
3087                   gc->covars[v] /= (float)(nsamples-1) ;
3088                   check_finite("GCAcompleteCovarianceTraining:"
3089                                " covar",
3090                                gc->covars[v]) ;
3091                   if (r == c)
3092                     /* diagonal should be positive definite */
3093                   {
3094                     if (gc->covars[v] < -0.1)
3095                       DiagBreak() ;
3096                   }
3097                 }
3098               }
3099             }
3100           }
3101           if (x == Ggca_x && y == Ggca_y && z == Ggca_z &&
3102               (gcan->labels[n] == Ggca_label || Ggca_label < 0))
3103           {
3104             MATRIX *m ;
3105             double det ;
3106 
3107             det = covariance_determinant(gc, gca->ninputs) ;
3108             printf("final covariance matrix for %s "
3109                    "(nsamples=%d, det=%f):\n",
3110                    cma_label_to_name(gcan->labels[n]),
3111                    nsamples, det) ;
3112             m = load_covariance_matrix(gc, NULL, gca->ninputs) ;
3113             MatrixPrint(stdout, m) ;
3114             MatrixFree(&m) ;
3115             fflush(stdout) ;
3116           }
3117         }
3118       }
3119     }
3120   }
3121 
3122   holes_filled = 0 ;
3123 #if 0
3124   printf("filling holes in the GCA...\n") ;
3125   for (holes_filled = x = 0 ; x < gca->node_width ; x++)
3126   {
3127     for (y = 0 ; y < gca->node_height ; y++)
3128     {
3129       for (z = 0 ; z < gca->node_depth ; z++)
3130       {
3131         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3132           DiagBreak() ;
3133         gcan = &gca->nodes[x][y][z] ;
3134         for (n = 0 ; n < gcan->nlabels ; n++)
3135         {
3136           if (x == Ggca_x && y == Ggca_y && z == Ggca_z &&
3137               (Ggca_label == gcan->labels[n] || Ggca_label < 0))
3138             DiagBreak() ;
3139           gc = &gcan->gcs[n] ;
3140           nsamples = gc->ntraining - gc->n_just_priors ;
3141           /* for no-intensity training */
3142 #define MIN_GC_SAMPLES(gca) ((gca->ninputs*(gca->ninputs+1))/2)
3143           if (gc->means[0] < 0 || nsamples < MIN_GC_SAMPLES(gca))
3144           {
3145             GC1D *gc_nbr ;
3146             int  r, c, i ;
3147 
3148             gc_nbr = findClosestValidGC(gca, x, y, z,
3149                                         gcan->labels[n], 1) ;
3150             if (!gc_nbr)
3151             {
3152               ErrorPrintf(ERROR_BADPARM,
3153                           "gca(%d,%d,%d,%d) - could not "
3154                           "find valid nbr label "
3155                           "%s (%d)", x,  y, z, n,
3156                           cma_label_to_name(gcan->labels[n]),
3157                           gcan->labels[n]) ;
3158               continue ;
3159             }
3160             holes_filled++ ;
3161             for (i = r = 0 ; r < gca->ninputs ; r++)
3162             {
3163               gc->means[r] = gc_nbr->means[r] ;
3164               for (c = r ; c < gca->ninputs ; c++, i++)
3165                 gc->covars[i] = gc_nbr->covars[i] ;
3166             }
3167             if (x == Ggca_x && y == Ggca_y && z == Ggca_z &&
3168                 (Ggca_label == gcan->labels[n] || Ggca_label < 0))
3169               printf("filling hole @ (%d, %d, %d)\n", x, y, z) ;
3170           }
3171         }
3172       }
3173     }
3174   }
3175 #endif
3176 
3177   GCAfixSingularCovarianceMatrices(gca) ;
3178   /* shouldn't need the code that follows, but haven't made sure yet */
3179 
3180   /* find and fix singular covariance matrices */
3181   for (x = 0 ; x < gca->node_width ; x++)
3182   {
3183     for (y = 0 ; y < gca->node_height ; y++)
3184     {
3185       for (z = 0 ; z < gca->node_depth ; z++)
3186       {
3187         gcan = &gca->nodes[x][y][z] ;
3188         for (n = 0 ; n < gcan->nlabels ; n++)
3189         {
3190           MATRIX *m_cov, *m_cov_inv ;
3191           double det ;
3192 
3193           if (x == Ggca_x && y == Ggca_y && z == Ggca_z &&
3194               (Ggca_label == gcan->labels[n] || Ggca_label < 0))
3195             DiagBreak() ;
3196           gc = &gcan->gcs[n] ;
3197           m_cov = load_covariance_matrix(gc, NULL, gca->ninputs) ;
3198           m_cov_inv = MatrixInverse(m_cov, NULL) ;
3199           det = covariance_determinant(gc, gca->ninputs) ;
3200           if (det <= 0 )
3201           {
3202             MatrixFree(&m_cov_inv) ;
3203             m_cov_inv = NULL ;
3204           }
3205           if (m_cov_inv == NULL)  /* singular */
3206           {
3207             MATRIX *m_I ;
3208             m_I = MatrixIdentity(gca->ninputs, NULL) ;
3209             MatrixScalarMul(m_I, MIN_VAR, m_I) ;
3210             MatrixAdd(m_cov, m_I, m_cov) ;
3211             m_cov_inv = MatrixInverse(m_cov, NULL) ;
3212             if (m_cov_inv == NULL)
3213               ErrorExit(ERROR_BADPARM,
3214                         "GCAcompleteCovarianceTraining: cannot "
3215                         "regularize singular covariance matrix at "
3216                         "(%d,%d,%d):%d",
3217                         x, y, z, n) ;
3218             det = covariance_determinant(gc, gca->ninputs) ;
3219             if (det < 0)
3220               DiagBreak() ;
3221             for (v = r = 0 ; r < gca->ninputs ; r++)
3222               for (c = r ; c < gca->ninputs ; c++, v++)
3223                 gc->covars[v] = *MATRIX_RELT(m_cov, r+1, c+1) ;
3224             MatrixFree(&m_I) ;
3225             nregs++ ;
3226           }
3227 
3228           MatrixFree(&m_cov_inv) ;
3229           MatrixFree(&m_cov) ;
3230         }
3231       }
3232     }
3233   }
3234   gcaCheck(gca) ;
3235 
3236   if (total_pruned > 0)
3237   {
3238     printf("%d samples pruned during training\n", total_pruned) ;
3239     total_pruned = 0 ;
3240   }
3241   return(NO_ERROR) ;
3242 }
3243 
3244 MRI  *
3245 GCAlabel(MRI *mri_inputs, GCA *gca, MRI *mri_dst, TRANSFORM *transform)
3246 {
3247   int       x, y, z, width, height, depth, label, xn, yn, zn, n, num_pv,
3248   use_partial_volume_stuff ;
3249   GCA_NODE  *gcan ;
3250   GCA_PRIOR *gcap ;
3251   GC1D      *gc ;
3252   float    /*dist,*/ max_p, p, vals[MAX_GCA_INPUTS] ;
3253 #if INTERP_PRIOR
3254   float     prior ;
3255 #endif
3256 
3257   use_partial_volume_stuff = (getenv("USE_PARTIAL_VOLUME_STUFF") != NULL);
3258   if (use_partial_volume_stuff)
3259     printf("using partial volume calculations in labeling...\n") ;
3260 
3261   // labeled volume has the same property of the inputs
3262   if (!mri_dst)
3263   {
3264     mri_dst = MRIalloc(mri_inputs->width,
3265                        mri_inputs->height,
3266                        mri_inputs->depth,
3267                        MRI_UCHAR) ;
3268     if (!mri_dst)
3269       ErrorExit(ERROR_NOMEMORY, "GCAlabel: could not allocate dst") ;
3270     MRIcopyHeader(mri_inputs, mri_dst) ;
3271   }
3272 
3273 
3274   /* go through each voxel in the input volume and find the canonical
3275      voxel (and hence the classifier) to which it maps. Then update the
3276      classifiers statistics based on this voxel's intensity and label.
3277   */
3278   width = mri_inputs->width ;
3279   height = mri_inputs->height;
3280   depth = mri_inputs->depth ;
3281   for (x = 0 ; x < width ; x++)
3282   {
3283     for (y = 0 ; y < height ; y++)
3284     {
3285       for (z = 0 ; z < depth ; z++)
3286       {
3287         if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3288           DiagBreak() ;
3289 
3290         if (x == width/2 && y == height/2 && z == depth/2)
3291           DiagBreak() ;
3292 
3293         if (!GCAsourceVoxelToNode(gca, mri_inputs,
3294                                   transform, x, y, z, &xn, &yn, &zn))
3295         {
3296           load_vals(mri_inputs, x, y, z, vals, gca->ninputs);
3297 
3298 #if 0
3299           if (x == 153 && y == 119 && z == 117)
3300             /* wm should be hippo (1484) */
3301           {
3302             Gx = xn ;
3303             Gy = yn ;
3304             Gz = zn ;
3305             DiagBreak() ;
3306           }
3307 #endif
3308 
3309           gcan = &gca->nodes[xn][yn][zn] ;
3310           gcap = getGCAP(gca, mri_inputs, transform, x, y, z) ;
3311           if (gcap==NULL)
3312             continue;
3313           label = 0 ;
3314           max_p = 2*GIBBS_NEIGHBORS*BIG_AND_NEGATIVE ;
3315           // going through gcap labels
3316           for (n = 0 ; n < gcap->nlabels ; n++)
3317           {
3318 #if 0
3319             p = GCAcomputePosteriorDensity(gcap, gcan, n, vals,
3320                                            gca->ninputs) ;
3321 #else
3322             gc = GCAfindGC(gca, xn, yn, zn, gcap->labels[n]) ;
3323             if (gc == NULL)
3324             {
3325               MRIsetVoxVal(mri_dst, x, y, z,0,0); // unknown
3326               continue ;
3327             }
3328 #if INTERP_PRIOR
3329             prior = gcaComputePrior(gca,
3330                                     mri_inputs,
3331                                     transform,
3332                                     x, y, z,
3333                                     gcap->labels[n]) ;
3334             p = gcaComputeLogDensity(gc, vals,
3335                                      gca->ninputs,
3336                                      prior,
3337                                      gcap->labels[n]) ;
3338 #else
3339             p = gcaComputeLogDensity(gc, vals,
3340                                      gca->ninputs,
3341                                      gcap->priors[n],
3342                                      gcap->labels[n]) ;
3343 #endif
3344 #endif
3345             // look for largest p
3346             if (p > max_p)
3347             {
3348               max_p = p ;
3349               label = gcap->labels[n] ;
3350             }
3351           }
3352           if (use_partial_volume_stuff)
3353             //////////// start of partial volume stuff
3354           {
3355             int n1, l1, l2, max_l1, max_l2, max_n1, max_n2 ;
3356             double max_p_pv ;
3357 
3358             max_p_pv = -10000  ;
3359             if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3360               DiagBreak() ;
3361             max_l1 = label ;
3362             max_l2 = max_n1 = max_n2 = 0 ;
3363             for (n = 0 ; n < gcap->nlabels ; n++)
3364               for (n1 = n+1 ; n1 < gcap->nlabels ; n1++)
3365               {
3366                 l1 = gcap->labels[n] ;
3367                 l2 = gcap->labels[n1] ;
3368                 p = compute_partial_volume_log_posterior
3369                     (gca, gcan, gcap, vals, l1, l2) ;
3370                 if (p > max_p_pv)
3371                 {
3372                   max_l1 = l1 ;
3373                   max_l2 = l2 ;
3374                   max_p_pv = p ;
3375                   max_n1 = n ;
3376                   max_n2 = n1 ;
3377                 }
3378                 if (p > max_p && l1 != label && l2 != label)
3379                   DiagBreak() ;
3380               }
3381 
3382             /* not the label picked before - change it */
3383             if (max_p_pv > max_p &&
3384                 max_l1 != label &&
3385                 max_l2 != label)
3386             {
3387               double p1, p2 ;
3388 
3389               gc = GCAfindGC(gca, xn, yn, zn, max_l1) ;
3390               p1 =
3391                 gcaComputeLogDensity(gc, vals, gca->ninputs,
3392                                      gcap->priors[max_n1],max_l1) ;
3393               gc = GCAfindGC(gca, xn, yn, zn, max_l2) ;
3394               p2 =
3395                 gcaComputeLogDensity(gc, vals, gca->ninputs,
3396                                      gcap->priors[max_n2],max_l2) ;
3397               num_pv++ ;
3398               if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3399                 printf("label @ %d, %d, %d: partial volume "
3400                        "from %s to %s\n",
3401                        x, y, z, cma_label_to_name(label),
3402                        cma_label_to_name
3403                        (p1 > p2 ? max_l1 : max_l2)) ;
3404               label = p1 > p2 ? max_l1 : max_l2 ;
3405               DiagBreak() ;
3406             }
3407           }
3408           //////////// end of partial volume stuff
3409 
3410           // found the label
3411           ///////////////////////// debug code /////////////////////
3412           if (x == Ggca_x && y == Ggca_y && z == Ggca_z)
3413           {
3414             int i ;
3415             printf("(%d, %d, %d): inputs=", x, y, z) ;
3416             for (i = 0 ; i < gca->ninputs ; i++)
3417               printf("%2.1f ", vals[i]) ;
3418 
3419             printf("\nprior label %s (%d), log(p)=%2.2e, "
3420                    "node (%d, %d, %d)\n",
3421                    cma_label_to_name(label), label, max_p,
3422                    xn, yn, zn) ;
3423             dump_gcan(gca, gcan, stdout, 1, gcap) ;
3424           }
3425           /////////////////////////////////////////////
3426           // set the value
3427           MRIsetVoxVal(mri_dst, x, y, z, 0, label) ;
3428         }
3429         else
3430           MRIsetVoxVal(mri_dst, x, y, z, 0, 0) ; // unknown
3431       } // z loop
3432     } // y loop
3433   } // x loop
3434 
3435   return(mri_dst) ;
3436 }
3437 
3438 MRI  *
3439 GCAlabelProbabilities(MRI *mri_inputs,
3440                       GCA *gca,
3441                       MRI *mri_dst,
3442                       TRANSFORM *transform)
3443 {
3444   int       x, y, z, width, height, depth, label,
3445   xn, yn, zn, n ;
3446   GCA_NODE  *gcan ;
3447   GCA_PRIOR *gcap ;
3448   GC1D      *gc ;
3449   double    max_p, p, total_p ;
3450   float     vals[MAX_GCA_INPUTS] ;
3451 
3452   width = mri_inputs->width ;
3453   height = mri_inputs->height;
3454   depth = mri_inputs->depth ;
3455   if (!mri_dst)
3456   {
3457     mri_dst = MRIalloc(width, height, depth, MRI_UCHAR) ;
3458     if (!mri_dst)
3459       ErrorExit(ERROR_NOMEMORY, "GCAlabel: could not allocate dst") ;
3460     MRIcopyHeader(mri_inputs, mri_dst) ;
3461   }
3462 
3463   /* go through each voxel in the input volume and find the canonical
3464      voxel (and hence the classifier) to which it maps. Then update the
3465      classifiers statistics based on this voxel's intensity and label.
3466   */
3467   for (x = 0 ; x < width ; x++)
3468   {
3469     for (y = 0 ; y < height ; y++)
3470     {
3471       for (z = 0 ; z < depth ; z++)
3472       {
3473         /// debug code /////////////////////////
3474         if (x == 152 && y == 126 && z == 127)
3475           DiagBreak() ;
3476         if (x == 63 && y == 107 && z == 120)
3477           DiagBreak() ;
3478         ///////////////////////////////////////
3479 
3480         load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
3481         if (!GCAsourceVoxelToNode(gca, mri_inputs,
3482                                   transform, x, y, z, &xn, &yn, &zn))
3483         {
3484           gcan = &gca->nodes[xn][yn][zn] ;
3485           gcap = getGCAP(gca, mri_inputs, transform, x, y, z) ;
3486           if (gcap == NULL || gcap->nlabels <= 0)
3487             continue;
3488           label = 0 ;
3489           max_p = 2*GIBBS_NEIGHBORS*BIG_AND_NEGATIVE ;
3490           // go through labels and find the one with max probability
3491           for (total_p = 0.0, n = 0 ; n < gcan->nlabels ; n++)
3492           {
3493             gc = &gcan->gcs[n] ;
3494 
3495             /* compute 1-d Mahalanobis distance */
3496             p = GCAcomputePosteriorDensity(gcap, gcan, n, vals,
3497                                            gca->ninputs) ;
3498             if (p > max_p)
3499             {
3500               max_p = p ;
3501               label = gcan->labels[n] ;
3502             }
3503             total_p += p ;
3504           }
3505           max_p = 255.0* max_p / total_p ;
3506           if (max_p > 255) max_p = 255 ;
3507           MRIsetVoxVal(mri_dst, x, y, z, 0, (BUFTYPE)max_p) ;
3508         }
3509         else
3510           MRIsetVoxVal(mri_dst, x, y, z, 0,255); // 0;
3511       }
3512     }
3513   }
3514 
3515   return(mri_dst) ;
3516 }
3517 
3518 MRI  *
3519 GCAcomputeProbabilities(MRI *mri_inputs, GCA *gca, MRI *mri_labels,
3520                         MRI *mri_dst, TRANSFORM *transform)
3521 {
3522   int       x, y, z, width, height, depth, label,
3523   xn, yn, zn, n ;
3524   GCA_NODE  *gcan ;
3525   GCA_PRIOR *gcap ;
3526   GC1D      *gc ;
3527   double    label_p, p, total_p ;
3528   float     vals[MAX_GCA_INPUTS] ;
3529 
3530   width = mri_inputs->width ;
3531   height = mri_inputs->height;
3532   depth = mri_inputs->depth ;
3533   if (!mri_dst)
3534   {
3535     mri_dst = MRIalloc(width, height, depth, MRI_UCHAR) ;
3536     if (!mri_dst)
3537       ErrorExit(ERROR_NOMEMORY, "GCAlabel: could not allocate dst") ;
3538     MRIcopyHeader(mri_inputs, mri_dst) ;
3539   }
3540 
3541 
3542   /* go through each voxel in the input volume and find the canonical
3543      voxel (and hence the classifier) to which it maps. Then update the
3544      classifiers statistics based on this voxel's intensity and label.
3545   */
3546   for (x = 0 ; x < width ; x++)
3547   {
3548     for (y = 0 ; y < height ; y++)
3549     {
3550       for (z = 0 ; z < depth ; z++)
3551       {
3552 #if 0
3553         if (x == 152 && y == 126 && z == 127)
3554           DiagBreak() ;
3555         if (x == 63 && y == 107 && z == 120)
3556           DiagBreak() ;
3557 #endif
3558         load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
3559         if (!GCAsourceVoxelToNode(gca, mri_inputs,
3560                                   transform, x, y, z, &xn, &yn, &zn))
3561         {
3562           label = nint(MRIgetVoxVal(mri_labels, x, y, z, 0)) ;
3563 
3564           gcan = &gca->nodes[xn][yn][zn] ;
3565           gcap = getGCAP(gca, mri_inputs, transform, x, y, z) ;
3566           if (gcap==NULL)
3567             continue;
3568           for (label_p = total_p = 0.0, n = 0 ;
3569                n < gcan->nlabels ;
3570                n++)
3571           {
3572             gc = &gcan->gcs[n] ;
3573 
3574             p = GCAcomputePosteriorDensity(gcap, gcan, n, vals,
3575                                            gca->ninputs) ;
3576             if (label == gcan->labels[n])
3577             {
3578               label_p = p ;
3579             }
3580             total_p += p ;
3581           }
3582           label_p = 255.0* label_p / total_p ;
3583           if (label_p > 255) label_p = 255 ;
3584           MRIsetVoxVal(mri_dst, x, y, z,0,(BUFTYPE)label_p) ;
3585         }
3586       }
3587     }
3588   }
3589 
3590   return(mri_dst) ;
3591 }
3592 
3593 #define STARTING_T   500
3594 
3595 MRI  *
3596 GCAannealUnlikelyVoxels(MRI *mri_inputs,
3597                         GCA *gca,
3598                         MRI *mri_dst,
3599                         TRANSFORM *transform,
3600                         int max_iter,
3601                         MRI  *mri_fixed)
3602 {
3603   int       x, y, z, width, depth, height, *x_indices, *y_indices, *z_indices,
3604   nindices, index, iter, nchanged, xn, yn, zn, n, nbad,
3605   old_label  ;
3606   double    log_likelihood, T, delta_E, p, rn, new_ll, old_ll,
3607   total_likelihood ;
3608   GCA_NODE  *gcan ;
3609   MRI       *mri_bad ;
3610 
3611   width = mri_inputs->width ;
3612   height = mri_inputs->height ;
3613   depth = mri_inputs->depth ;
3614   mri_bad = MRIalloc(width, height, depth, MRI_UCHAR) ;
3615   MRIcopyHeader(mri_inputs, mri_bad);
3616 
3617   for (nindices = x = 0 ; x < width ; x++)
3618   {
3619     for (y = 0 ; y < height ; y++)
3620     {
3621       for (z = 0 ; z < depth ; z++)
3622       {
3623         if (mri_fixed && MRIvox(mri_fixed, x, y, z) > 0)
3624           continue ;
3625 
3626         if (!GCAsourceVoxelToNode(gca, mri_inputs,
3627                                   transform, x, y, z, &xn, &yn, &zn))
3628         {
3629           log_likelihood =
3630             gcaVoxelGibbsLogLikelihood(gca, mri_dst,
3631                                        mri_inputs, x, y, z,transform,
3632                                        PRIOR_FACTOR);
3633           gcan = &gca->nodes[xn][yn][zn] ;
3634           if (log_likelihood < log(1.0f/(float)gcan->total_training))
3635           {
3636             MRIvox(mri_bad, x, y, z) = 1 ;
3637             nindices++ ;
3638           }
3639         }
3640       }
3641     }
3642   }
3643 
3644   printf("annealing %d voxels...\n", nindices) ;
3645   x_indices = (int *)calloc(nindices, sizeof(int)) ;
3646   y_indices = (int *)calloc(nindices, sizeof(int)) ;
3647   z_indices = (int *)calloc(nindices, sizeof(int)) ;
3648   for (index = x = 0 ; x < width ; x++)
3649   {
3650     for (y = 0 ; y < height ; y++)
3651     {
3652       for (z = 0 ; z < depth ; z++)
3653       {
3654         if (MRIvox(mri_bad, x, y, z) > 0)
3655         {
3656           x_indices[index] = x ;
3657           y_indices[index] = y ;
3658           z_indices[index] = z ;
3659           index++ ;
3660         }
3661       }
3662     }
3663   }
3664 
3665   MRIfree(&mri_bad) ;
3666   T = STARTING_T ;
3667   iter = 0 ;
3668   do
3669   {
3670     total_likelihood = 0.0 ;
3671     for (nbad = nchanged = index = 0 ; index < nindices ; index++)
3672     {
3673       x = x_indices[index] ;
3674       y = y_indices[index] ;
3675       z = z_indices[index] ;
3676       if (x == 155 && y == 126 && z == 128)
3677         DiagBreak() ;
3678 
3679       /* find the node associated with this coordinate and classify */
3680       if (!GCAsourceVoxelToNode(gca, mri_inputs, transform, x, y, z,
3681                                 &xn, &yn, &zn))
3682       {
3683         gcan = &gca->nodes[xn][yn][zn] ;
3684 
3685         if (gcan->nlabels == 1)
3686           continue ;
3687         n = (int)randomNumber(0.0, (double)gcan->nlabels-0.0001) ;
3688         if (gcan->labels[n] == nint(MRIgetVoxVal(mri_dst, x, y, z,0)))
3689           continue ;
3690         old_ll =
3691           gcaNbhdGibbsLogLikelihood(gca, mri_dst, mri_inputs,
3692                                     x, y, z, transform,
3693                                     PRIOR_FACTOR) ;
3694         old_label = nint(MRIgetVoxVal(mri_dst, x, y, z,0)) ;
3695         MRIsetVoxVal(mri_dst, x, y, z, 0, gcan->labels[n]) ;
3696         new_ll =
3697           gcaNbhdGibbsLogLikelihood(gca, mri_dst, mri_inputs,
3698                                     x, y, z, transform,
3699                                     PRIOR_FACTOR) ;
3700         delta_E = new_ll - old_ll ;
3701         p = exp(delta_E / T) ;
3702         rn = randomNumber(0.0, 1.0) ;
3703 
3704         if (p > rn)
3705         {
3706           if (new_ll < log(1.0f/(float)gcan->total_training))
3707             nbad++ ;
3708           nchanged++ ;
3709           total_likelihood += new_ll ;
3710         }
3711         else
3712         {
3713           total_likelihood += old_ll ;
3714           if (old_ll < log(1.0f/(float)gcan->total_training))
3715             nbad++ ;
3716           MRIsetVoxVal(mri_dst, x, y, z,0,old_label) ;
3717         }
3718       }
3719     }
3720     T = T * 0.99 ;
3721     fprintf(stdout, "%03d: T = %2.2f, nchanged %d, nbad = %d, ll=%2.2f\n",
3722             iter, T, nchanged, nbad, total_likelihood/(double)nindices) ;
3723     if (!nchanged)
3724       break ;
3725   }
3726   while (iter++ < max_iter) ;
3727 
3728   free(x_indices) ;
3729   free(y_indices) ;
3730   free(z_indices) ;
3731   fflush(stdout) ;
3732   return(mri_dst) ;
3733 }
3734 
3735 GCA  *
3736 GCAreduce(GCA *gca_src)
3737 {
3738 #if 0
3739   GCA       *gca_dst ;
3740   int       xs, ys, zs, xd, yd, zd, swidth, sheight, sdepth,
3741   ns, nd, dwidth, dheight, ddepth, dof ;
3742   float     spacing = gca_src->spacing ;
3743   GCA_NODE  *gcan_src, *gcan_dst ;
3744   GC1D      *gc_src, *gc_dst ;
3745 
3746   swidth = gca_src->width ;
3747   sheight = gca_src->height;
3748   sdepth = gca_src->depth;
3749 
3750   gca_dst =
3751     GCAalloc(gca_src->ninputs, spacing*2, spacing*swidth,
3752              spacing*gca_src->height,spacing*sdepth, gca_src->flags) ;
3753 
3754 
3755   dwidth = gca_dst->width ;
3756   dheight = gca_dst->height;
3757   ddepth = gca_dst->depth;
3758 
3759   for (zs = 0 ; zs < sdepth ; zs++)
3760   {
3761     for (ys = 0 ; ys < sheight ; ys++)
3762     {
3763       for (xs = 0 ; xs < swidth ; xs++)
3764       {
3765         xd = xs/2 ;
3766         yd = ys/2 ;
3767         zd = zs/2 ;
3768         if (xd == 15 && yd == 22 && zd == 35)
3769           DiagBreak() ;
3770         gcan_src = &gca_src->nodes[xs][ys][zs] ;
3771         gcan_dst = &gca_dst->nodes[xd][yd][zd] ;
3772 
3773         for (ns = 0 ; ns < gcan_src->nlabels ; ns++)
3774         {
3775           gc_src = &gcan_src->gcs[ns] ;
3776           dof = gc_src->prior * gcan_src->total_training ;
3777 
3778           /* find label in destination */
3779           for (nd = 0 ; nd < gcan_dst->nlabels ; nd++)
3780           {
3781             if (gcan_dst->labels[nd] == gcan_src->labels[ns])
3782               break ;
3783           }
3784           if (nd >= gcan_dst->nlabels)  /* couldn't find it */
3785           {
3786             if (nd >= gcan_dst->max_labels)
3787             {
3788               int  old_max_labels ;
3789               char *old_labels ;
3790               GC1D *old_gcs ;
3791 
3792               old_max_labels = gcan_dst->max_labels ;
3793               gcan_dst->max_labels += 2 ;
3794               old_labels = gcan_dst->labels ;
3795               old_gcs = gcan_dst->gcs ;
3796 
3797               /* allocate new ones */
3798 #if 0
3799               gcan_dst->gcs =
3800                 (GC1D *)calloc(gcan_dst->max_labels,
3801                                sizeof(GC1D)) ;
3802 #else
3803               gcan_dst->gcs =
3804                 alloc_gcs(gcan_dst->max_labels,
3805                           gca_dst->flags, gca_dst->ninputs) ;
3806 #endif
3807               if (!gcan_dst->gcs)
3808                 ErrorExit(ERROR_NOMEMORY,
3809                           "GCANreduce: couldn't expand gcs to %d",
3810                           gcan_dst->max_labels) ;
3811               gcan_dst->labels =
3812                 (unsigned short *)calloc(gcan_dst->max_labels,
3813                                         sizeof(unsigned short)) ;
3814               if (!gcan_dst->labels)
3815                 ErrorExit(ERROR_NOMEMORY,
3816                           "GCANupdateNode: couldn't expand "
3817                           "labels to %d",
3818                           gcan_dst->max_labels) ;
3819 
3820               /* copy the old ones over */
3821 #if 0
3822               memmove(gcan_dst->gcs, old_gcs,
3823                       old_max_labels*sizeof(GC1D));
3824 #else
3825               copy_gcs(old_max_labels, old_gcs,
3826                        gcan_dst->gcs, gca->ninputs) ;
3827 #endif
3828 
3829               memmove(gcan_dst->labels, old_labels,
3830                       old_max_labels*sizeof(unsigned short)) ;
3831 
3832               /* free the old ones */
3833               free(old_gcs) ;
3834               free(old_labels) ;
3835             }
3836             gcan_dst->nlabels++ ;
3837           }
3838           gc_dst = &gcan_dst->gcs[nd] ;
3839           gc_dst->mean += dof*gc_src->mean ;
3840           gc_dst->var += dof*(gc_src->var) ;
3841           gc_dst->prior += dof ;
3842           gcan_dst->total_training += dof ;
3843           gcan_dst->labels[nd] = gcan_src->labels[ns] ;
3844         }
3845       }
3846     }
3847   }
3848 
3849   for (xd = 0 ; xd < gca_dst->width ; xd++)
3850   {
3851     for (yd = 0 ; yd < gca_dst->height ; yd++)
3852     {
3853       for (zd = 0 ; zd < gca_dst->depth ; zd++)
3854       {
3855         gcan_dst = &gca_dst->nodes[xd][yd][zd] ;
3856         for (nd = 0 ; nd < gcan_dst->nlabels ; nd++)
3857         {
3858           float var ;
3859 
3860           gc_dst = &gcan_dst->gcs[nd] ;
3861           dof = gc_dst->prior ;
3862           /* prior is count of # of occurences now */
3863           gc_dst->mean /= dof ;
3864           var = gc_dst->var / dof ;
3865           if (var < -0.1)
3866             DiagBreak() ;
3867           if (var < MIN_VAR)
3868             var = MIN_VAR ;
3869           gc_dst->var = var ;
3870           gc_dst->prior /= (float)gcan_dst->total_training ;
3871         }
3872       }
3873     }
3874   }
3875   return(gca_dst) ;
3876 #else
3877   /* have to update to include priors at different resolution than nodes */
3878   ErrorReturn(NULL, (ERROR_UNSUPPORTED,
3879                      "GCAreduce: not currently supported") ;) ;
3880 #endif
3881 }
3882 
3883 
3884 typedef struct
3885 {
3886   int   label ;
3887   float prob ;
3888   int   index ;
3889 }
3890 LABEL_PROB ;
3891 
3892 static int compare_sort_probabilities(const void *plp1, const void *plp2);
3893 MRI *
3894 GCAclassify(MRI *mri_inputs,GCA *gca,MRI *mri_dst,
3895             TRANSFORM *transform,int max_labels)
3896 {
3897   int        x, y, z, width, height, depth,
3898   xn, yn, zn, n ;
3899   GCA_NODE   *gcan ;
3900   GCA_PRIOR  *gcap ;
3901   GC1D       *gc ;
3902   float      max_p, p, total_p ;
3903   LABEL_PROB label_probs[1000] ;
3904   float      vals[MAX_GCA_INPUTS] ;
3905 
3906   if (max_labels > MAX_LABELS_PER_GCAN || max_labels <= 0)
3907     max_labels = MAX_LABELS_PER_GCAN ;
3908 
3909   if (!mri_dst)
3910   {
3911     int width = mri_inputs->width, height = mri_inputs->height,
3912                                             depth = mri_inputs->depth ;
3913 
3914     mri_dst = MRIallocSequence(width, height, depth,
3915                                MRI_UCHAR, 2*max_labels) ;
3916     if (!mri_dst)
3917       ErrorExit(ERROR_NOMEMORY, "GCAlabel: could not allocate dst") ;
3918     MRIcopyHeader(mri_inputs, mri_dst) ;
3919   }
3920 
3921   /* go through each voxel in the input volume and find the canonical
3922      voxel (and hence the classifier) to which it maps. Then update the
3923      classifiers statistics based on this voxel's intensity and label.
3924   */
3925   width = mri_inputs->width ;
3926   height = mri_inputs->height;
3927   depth = mri_inputs->depth ;
3928   for (x = 0 ; x < width ; x++)
3929   {
3930     for (y = 0 ; y < height ; y++)
3931     {
3932       for (z = 0 ; z < depth ; z++)
3933       {
3934         if (x == 67 && y == 87 && z == 114)
3935           DiagBreak() ;
3936 
3937         load_vals(mri_inputs, x, y, z, vals, gca->ninputs) ;
3938 
3939         /* find the node associated with this coordinate and classify */
3940         if (!GCAsourceVoxelToNode(gca, mri_inputs,
3941                                   transform, x, y, z, &xn, &yn, &zn))
3942         {
3943           gcan = &gca->nodes[xn][yn][zn] ;
3944           gcap = getGCAP(gca, mri_inputs, transform, x, y, z) ;
3945           if (gcap==NULL)
3946             continue;
3947           max_p = 2*GIBBS_NEIGHBORS*BIG_AND_NEGATIVE ;
3948           for (total_p = 0.0, n = 0 ; n < gcan->nlabels ; n++)
3949           {
3950             gc = &gcan->gcs[n] ;
3951 
3952             p = GCAcomputePosteriorDensity(gcap, gcan, n, vals,
3953                                            gca->ninputs) ;
3954             total_p += p ;
3955             label_probs[n].prob = p ;
3956             label_probs[n].label = gcan->labels[n] ;
3957           }
3958           /* now sort the labels and probabilities */
3959           qsort(label_probs, gcan->nlabels, sizeof(LABEL_PROB),
3960                 compare_sort_probabilities) ;
3961 
3962           for (n = 0 ; n < max_labels ; n++)
3963           {
3964             if (n < gcan->nlabels)
3965             {
3966               MRIseq_vox(mri_dst, x, y, z, n*2) =
3967                 label_probs[n].label ;
3968               MRIseq_vox(mri_dst, x, y, z, n*2+1) =
3969                 (BUFTYPE)nint(255.0*label_probs[n].prob/total_p) ;
3970             }
3971             else
3972             {
3973               MRIseq_vox(mri_dst, x, y, z, n*2) = 255 ;
3974               MRIseq_vox(mri_dst, x, y, z, n*2+1) = 0 ;
3975             }
3976           }
3977         }
3978       }
3979     }
3980   }
3981 
3982   return(mri_dst) ;
3983 }
3984 
3985 
3986 static int