Attachment 'mri.c'

Download

   1 /*
   2  *       FILE NAME:   mri.c
   3  *
   4  *       DESCRIPTION: utilities for MRI  data structure
   5  *
   6  *       AUTHOR:      Bruce Fischl
   7  *       DATE:        1/8/97
   8  *
   9  */
  10 // Warning: Do not edit the following four lines.  CVS maintains them.
  11 // Revision Author: $Author: nicks $
  12 // Revision Date  : $Date: 2006/01/17 22:16:34 $
  13 // Revision       : $Revision: 1.325 $
  14 char *MRI_C_VERSION = "$Revision: 1.325 $";
  15 
  16 /*-----------------------------------------------------
  17   INCLUDE FILES
  18   -------------------------------------------------------*/
  19 #define USE_ELECTRIC_FENCE  1
  20 
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <math.h>
  24 #include <string.h>
  25 #include <memory.h>
  26 #include <errno.h>
  27 #include <ctype.h>
  28 #include "error.h"
  29 #include "proto.h"
  30 #include "mri.h"
  31 #include "macros.h"
  32 #include "diag.h"
  33 #include "volume_io.h"
  34 #include "filter.h"
  35 #include "box.h"
  36 #include "region.h"
  37 #include "mri_transform.h"
  38 #include "utils.h"
  39 #include "matrix.h"
  40 #include "pdf.h"
  41 #include "cma.h"
  42 #include "talairachex.h"
  43 #include "voxlist.h"
  44 
  45 extern int errno;
  46 
  47 /*-----------------------------------------------------
  48   MACROS AND CONSTANTS
  49   -------------------------------------------------------*/
  50 
  51 #define DEBUG_POINT(x,y,z)  (((x==8&&y==9) || (x==9&&y==8)) &&((z)==15))
  52 
  53 /*-----------------------------------------------------
  54   STATIC DATA
  55   -------------------------------------------------------*/
  56 
  57 static long mris_alloced = 0 ;
  58 
  59 /*-----------------------------------------------------
  60   STATIC PROTOTYPES
  61   -------------------------------------------------------*/
  62 
  63 /*-----------------------------------------------------
  64   GLOBAL FUNCTIONS
  65   -------------------------------------------------------*/
  66 
  67 /*----------------------------------------------------------
  68   MRIxfmCRS2XYZ() - computes the matrix needed to compute the
  69   XYZ of the center of a voxel at a given Col, Row, and Slice
  70   from the native geometry of the volume
  71 
  72   x         col
  73   y  = T *  row
  74   z        slice
  75   1          1
  76 
  77   T = [Mdc*D Pxyz0]
  78   [0 0 0   1  ]
  79 
  80   Mdc = [Vcol Vrow Vslice]
  81   V<dim> = the direction cosine pointing from the center of one voxel
  82   to the center of an adjacent voxel in the next dim, where
  83   dim is either colum, row, or slice. Vcol = [x_r x_a x_s],
  84   Vrow = [y_r y_a y_s], Vslice = [z_r z_a z_s]. Vcol can also
  85   be described as the vector normal to the plane formed by
  86   the rows and slices of a given column (ie, the column normal).
  87 
  88   D = diag([colres rowres sliceres])
  89   dimres = the distance between adjacent dim, where colres = mri->xsize,
  90   rowres = mri->ysize, and sliceres = mri->zsize.
  91 
  92   Pxyz0 = the XYZ location at CRS=0. This number is not part of the
  93   mri structure, so it is computed here according to the formula:
  94   Pxyz0 = PxyzCenter - Mdc*D*PcrsCenter
  95 
  96   PcrsCenter = the col, row, and slice at the center of the volume,
  97   = [ ncols/2 nrows/2 nslices/2 ]
  98 
  99   PxyzCenter = the X, Y, and Z at the center of the volume and does
 100   exist in the header as mri->c_r, mri->c_a, and mri->c_s,
 101   respectively.
 102 
 103   Note: to compute the matrix with respect to the first voxel being
 104   at CRS 1,1,1 instead of 0,0,0, then set base = 1. This is
 105   necessary with SPM matrices.
 106 
 107   See also: MRIxfmCRS2XYZtkreg, MRItkReg2Native
 108   ------------------------------------------------------*/
 109 MATRIX *MRIxfmCRS2XYZ(MRI *mri, int base)
 110 {
 111   MATRIX *m;
 112   MATRIX *Pcrs, *PxyzOffset;
 113 
 114   m = MatrixAlloc(4, 4, MATRIX_REAL);
 115 
 116   /* direction cosine between columns scaled by
 117      distance between colums */
 118   *MATRIX_RELT(m, 1, 1) = mri->x_r * mri->xsize;
 119   *MATRIX_RELT(m, 2, 1) = mri->x_a * mri->xsize;
 120   *MATRIX_RELT(m, 3, 1) = mri->x_s * mri->xsize;
 121 
 122   /* direction cosine between rows scaled by
 123      distance between rows */
 124   *MATRIX_RELT(m, 1, 2) = mri->y_r * mri->ysize;
 125   *MATRIX_RELT(m, 2, 2) = mri->y_a * mri->ysize;
 126   *MATRIX_RELT(m, 3, 2) = mri->y_s * mri->ysize;
 127 
 128   /* direction cosine between slices scaled by
 129      distance between slices */
 130   *MATRIX_RELT(m, 1, 3) = mri->z_r * mri->zsize;
 131   *MATRIX_RELT(m, 2, 3) = mri->z_a * mri->zsize;
 132   *MATRIX_RELT(m, 3, 3) = mri->z_s * mri->zsize;
 133 
 134   /* Preset the offsets to 0 */
 135   *MATRIX_RELT(m, 1, 4) = 0.0;
 136   *MATRIX_RELT(m, 2, 4) = 0.0;
 137   *MATRIX_RELT(m, 3, 4) = 0.0;
 138 
 139   /* Last row of matrix */
 140   *MATRIX_RELT(m, 4, 1) = 0.0;
 141   *MATRIX_RELT(m, 4, 2) = 0.0;
 142   *MATRIX_RELT(m, 4, 3) = 0.0;
 143   *MATRIX_RELT(m, 4, 4) = 1.0;
 144 
 145   /* At this point, m = Mdc * D */
 146 
 147   /* Col, Row, Slice at the Center of the Volume */
 148   Pcrs = MatrixAlloc(4, 1, MATRIX_REAL);
 149   *MATRIX_RELT(Pcrs, 1, 1) = mri->width/2.0  + base;
 150   *MATRIX_RELT(Pcrs, 2, 1) = mri->height/2.0 + base;
 151   *MATRIX_RELT(Pcrs, 3, 1) = mri->depth/2.0  + base;
 152   *MATRIX_RELT(Pcrs, 4, 1) = 1.0;
 153 
 154   /* XYZ offset the first Col, Row, and Slice from Center */
 155   /* PxyzOffset = Mdc*D*PcrsCenter */
 156   PxyzOffset = MatrixMultiply(m,Pcrs,NULL);
 157 
 158   /* XYZ at the Center of the Volume is mri->c_r, c_a, c_s  */
 159 
 160   /* The location of the center of the voxel at CRS = (0,0,0)*/
 161   *MATRIX_RELT(m, 1, 4) = mri->c_r - PxyzOffset->rptr[1][1];
 162   *MATRIX_RELT(m, 2, 4) = mri->c_a - PxyzOffset->rptr[2][1];
 163   *MATRIX_RELT(m, 3, 4) = mri->c_s - PxyzOffset->rptr[3][1];
 164 
 165   MatrixFree(&Pcrs);
 166   MatrixFree(&PxyzOffset);
 167 
 168   return(m);
 169 }
 170 /*-------------------------------------------------------------
 171   MRIxfmCRS2XYZtkreg() - computes the linear transform between the
 172   column, row, and slice of a voxel and the x, y, z of that voxel as
 173   expected by tkregister (or for when using a tkregister compatible
 174   matrix). For tkregister, the column DC points in the "x" direction,
 175   the row DC points in the "z" direction, and the slice DC points in
 176   the "y" direction. The center of the coordinates is set to the
 177   center of the FOV. These definitions are arbitrary (and more than a
 178   little confusing). Since they are arbitrary, they must be applied
 179   consistently.
 180   -------------------------------------------------------------*/
 181 MATRIX *MRIxfmCRS2XYZtkreg(MRI *mri)
 182 {
 183   MRI *tmp;
 184   MATRIX *K;
 185 
 186   tmp = MRIallocHeader(mri->width, mri->height, mri->depth, mri->type);
 187 
 188   /* Set tkregister defaults */
 189   /* column         row           slice          center      */
 190   tmp->x_r = -1; tmp->y_r =  0; tmp->z_r =  0; tmp->c_r = 0.0;
 191   tmp->x_a =  0; tmp->y_a =  0; tmp->z_a =  1; tmp->c_a = 0.0;
 192   tmp->x_s =  0; tmp->y_s = -1; tmp->z_s =  0; tmp->c_s = 0.0;
 193 
 194   /* Copy the voxel resolutions */
 195   tmp->xsize = mri->xsize;
 196   tmp->ysize = mri->ysize;
 197   tmp->zsize = mri->zsize;
 198 
 199   K = MRIxfmCRS2XYZ(tmp,0);
 200 
 201   MRIfree(&tmp);
 202 
 203   return(K);
 204 }
 205 /*-------------------------------------------------------------------
 206   MRItkReg2Native() - converts a tkregister-compatible registration matrix
 207   R to one that works with the geometry native to the two volumes. The
 208   matrix R maps tkXYZ of the ref volume to tkXYZ of the mov volume. In
 209   a typical application, ref is the anatomical volume and mov is the
 210   functional volume. The purpose of this function is to be able to use
 211   registration matrices computed by (or compatible with) tkregister
 212   without losing the geometries native to the volumes. If R is null,
 213   it is assumed to be the identity. R: MovXYZ = R*RefXYZ. Typically,
 214   Ref is the Anatomical Reference, and Mov is the functional.
 215 
 216   See also: MRItkRegMtx, MRIxfmCRS2XYZtkreg, MRIxfmCRS2XYZ
 217   -------------------------------------------------------------------*/
 218 MATRIX *MRItkReg2Native(MRI *ref, MRI *mov, MATRIX *R)
 219 {
 220   MATRIX *Kref, *Kmov;
 221   MATRIX *Tref, *Tmov, *D;
 222   MATRIX *invKmov, *invTref;
 223 
 224   Tref = MRIxfmCRS2XYZ(ref,0);
 225   Tmov = MRIxfmCRS2XYZ(mov,0);
 226 
 227   Kref = MRIxfmCRS2XYZtkreg(ref);
 228   Kmov = MRIxfmCRS2XYZtkreg(mov);
 229 
 230   /* D = Tmov * inv(Kmov) * R * Kref *inv(Tref) */
 231   invKmov = MatrixInverse(Kmov,NULL);
 232   invTref = MatrixInverse(Tref,NULL);
 233 
 234   D = MatrixMultiply(Tmov,invKmov,NULL);
 235   if(R!=NULL) MatrixMultiply(D,R,D);
 236   MatrixMultiply(D,Kref,D);
 237   MatrixMultiply(D,invTref,D);
 238 
 239   if(0) {
 240     printf("MRITkReg2Native -----------------------------\n");
 241     printf("Tref ----------------\n");
 242     MatrixPrint(stdout,Tref);
 243     printf("Tmov ----------------\n");
 244     MatrixPrint(stdout,Tmov);
 245     printf("Kref ----------------\n");
 246     MatrixPrint(stdout,Kref);
 247     printf("Kmov ----------------\n");
 248     MatrixPrint(stdout,Kmov);
 249     printf("------------------------------------------\n");
 250   }
 251 
 252   MatrixFree(&Kref);
 253   MatrixFree(&Tref);
 254   MatrixFree(&Kmov);
 255   MatrixFree(&Tmov);
 256   MatrixFree(&invKmov);
 257   MatrixFree(&invTref);
 258 
 259   return(D);
 260 }
 261 /*----------------------------------------------------------------
 262   MRItkRegMtx() - creates a tkregsiter-compatible matrix from the
 263   matrix D that aligns the two volumes assuming the native geometry.
 264   This is the counterpart to MRITkReg2Native(). If D is null, it is
 265   assumed to be the identity. R: MovXYZ = R*RefXYZ. Typically,
 266   Ref is the Anatomical Reference, and Mov is the functional.
 267   Use this function with D=NULL when the two volumes have been
 268   aligned with SPM.
 269   ---------------------------------------------------------------*/
 270 MATRIX *MRItkRegMtx(MRI *ref, MRI *mov, MATRIX *D)
 271 {
 272   MATRIX *Kref, *Kmov;
 273   MATRIX *Tref, *Tmov, *R;
 274   MATRIX *invTmov, *invKref;
 275 
 276   /* Native Goemetry */
 277   Tref = MRIxfmCRS2XYZ(ref,0);
 278   Tmov = MRIxfmCRS2XYZ(mov,0);
 279 
 280   /* TkReg Goemetry */
 281   Kref = MRIxfmCRS2XYZtkreg(ref);
 282   Kmov = MRIxfmCRS2XYZtkreg(mov);
 283 
 284   invTmov = MatrixInverse(Tmov,NULL);
 285   invKref = MatrixInverse(Kref,NULL);
 286 
 287   /* R = Kmov * inv(Tmov) * D * Tref *inv(Kref) */
 288   R = MatrixMultiply(Kmov,invTmov,NULL);
 289   if(D != NULL) MatrixMultiply(R,D,R);
 290   MatrixMultiply(R,Tref,R);
 291   MatrixMultiply(R,invKref,R);
 292 
 293   MatrixFree(&Kref);
 294   MatrixFree(&Tref);
 295   MatrixFree(&Kmov);
 296   MatrixFree(&Tmov);
 297   MatrixFree(&invTmov);
 298   MatrixFree(&invKref);
 299 
 300   return(R);
 301 }
 302 /*-------------------------------------------------------------
 303   MRIfixTkReg() - this routine will adjust a matrix created by the
 304   "old" tkregister program. The old program had a problem in the way
 305   it chose the CRS of a voxel in the functional volume based on a
 306   point in the anatomical volume. The functional CRS of a point in
 307   anatomical space rarely (if ever) falls directly on a functional
 308   voxel, so it's necessary to choose a functional voxel given that the
 309   point falls between functional voxels (or it can be interpolated).
 310   The old tkregister program did not interpolate, rather it would
 311   choose the CRS in the following way: iC = floor(fC), iR = ceil(fR),
 312   and iS = floor(fS), where iC is the integer column number and fC is
 313   the floating point column, etc. Unfortunately, this is not nearest
 314   neighbor and it's not invertible. The right way to do it is to do
 315   nearest neighbor (ie, round to the closest integer). Unfortunately,
 316   there are a lot of data sets out there that have been regsitered
 317   with the old program, and we don't want to force poeple to
 318   reregister with the "new" program. This routine attempts to adjust
 319   the matrix created with the old program so that it will work with
 320   code that assumes that pure nearest neighbor was used.
 321 
 322   It does this by randomly sampling the anatomical volume in xyz
 323   and computing the tkreg'ed CRS for each point.
 324 
 325   Pcrs = inv(Tmov)*R*Pxyz
 326   PcrsTkReg = fcf(Pcrs) -- fcf is floor ceiling floor
 327 
 328   We seek a new R (Rfix) define with
 329 
 330   PcrsFix = inv(Tmov)*Rfix*Pxyz
 331 
 332   such that that the difference between PcrsFix and PcrsTkReg are
 333   minimized. To do this, we set
 334 
 335   PcrsFix = PcrsTkReg = inv(Tmov)*Rfix*Pxyz
 336 
 337   and solve for Rfix (this is an LMS solution):
 338 
 339   Rfix = Tmov*(PcrsTkReg*Pxyz')*inv(Pxyz*Pxyz');
 340 
 341   Applications that read in the registration matrix should detect the
 342   truncation method used (see below), fix the matrix if necessary, and
 343   proceed as if nearest neighbor/rounding was used.  The type of
 344   truncation can be determined from the last line of the registration
 345   file (after the matrix itself). If there is nothing there or the
 346   string "tkregister" is there, then the matrix should be
 347   converted. Otherwise, the string "round" should be there. The
 348   function regio_read_register (from registerio.c) will return the
 349   type of matrix in the float2int variable. It will be either
 350   FLT2INT_TKREG or FLT2INT_ROUND (constants defined in resample.h).
 351   ---------------------------------------------------------------*/
 352 MATRIX *MRIfixTkReg(MRI *mov, MATRIX *R)
 353 {
 354   int n, ntest = 1000;
 355   MATRIX *Pxyz, *Pcrs, *PcrsTkReg;
 356   MATRIX *PxyzT, *PxyzPxyzT, *invPxyzPxyzT;
 357   MATRIX *Tmov,*invTmov,*Rfix;
 358   MATRIX *tmp;
 359   float xrange, yrange, zrange;
 360 
 361   /* Assume a COR reference image */
 362   xrange = 256.0;
 363   yrange = 256.0;
 364   zrange = 256.0;
 365 
 366   Tmov = MRIxfmCRS2XYZtkreg(mov);
 367   invTmov = MatrixInverse(Tmov,NULL);
 368 
 369   Pxyz      = MatrixAlloc(4,ntest,MATRIX_REAL);
 370   PcrsTkReg = MatrixAlloc(4,ntest,MATRIX_REAL);
 371 
 372   /* Fill xyz with rand within the reference volume range */
 373   for(n=0;n<ntest;n++) {
 374     Pxyz->rptr[1][n+1] = xrange * (drand48()-0.5);
 375     Pxyz->rptr[2][n+1] = yrange * (drand48()-0.5);
 376     Pxyz->rptr[3][n+1] = zrange * (drand48()-0.5);
 377     Pxyz->rptr[4][n+1] = 1;
 378   }
 379 
 380   /* Compute floating mov CRS from targ XYZ */
 381   /* Pcrs = inv(Tmov)*R*Pxyz */
 382   tmp = MatrixMultiply(R,Pxyz,NULL);
 383   Pcrs = MatrixMultiply(invTmov,tmp,NULL);
 384   MatrixFree(&tmp);
 385 
 386   /* Truncate floating mov CRS using tkregister method*/
 387   for(n=0;n<ntest;n++) {
 388     PcrsTkReg->rptr[1][n+1] = floor(Pcrs->rptr[1][n+1]);
 389     PcrsTkReg->rptr[2][n+1] =  ceil(Pcrs->rptr[2][n+1]);
 390     PcrsTkReg->rptr[3][n+1] = floor(Pcrs->rptr[3][n+1]);
 391     PcrsTkReg->rptr[4][n+1] = 1;
 392   }
 393   MatrixFree(&Pcrs);
 394 
 395   //Rfix = Tmov*(PcrsTkreg*Pxyz')*inv(Pxyz*Pxyz');
 396   PxyzT = MatrixTranspose(Pxyz,NULL);
 397   PxyzPxyzT = MatrixMultiply(Pxyz,PxyzT,NULL);
 398   invPxyzPxyzT = MatrixInverse(PxyzPxyzT,NULL);
 399   tmp = MatrixMultiply(PcrsTkReg,PxyzT,NULL);
 400   MatrixMultiply(Tmov,tmp,tmp);
 401   Rfix = MatrixMultiply(tmp,invPxyzPxyzT,NULL);
 402 
 403 
 404   MatrixFree(&Pxyz);
 405   MatrixFree(&PcrsTkReg);
 406   MatrixFree(&PxyzT);
 407   MatrixFree(&PxyzPxyzT);
 408   MatrixFree(&invPxyzPxyzT);
 409   MatrixFree(&Tmov);
 410   MatrixFree(&invTmov);
 411   MatrixFree(&tmp);
 412 
 413   return(Rfix);
 414 }
 415 
 416 /*-------------------------------------------------------------------
 417   MRIfsl2TkReg() - converts an FSL registration matrix to one
 418   compatible with tkregister.  Note: the FSL matrix is assumed to map
 419   from the mov to the ref whereas the tkreg matrix maps from the ref
 420   to the mov.
 421 
 422   mov voxel --(Tmov)--> tkRegXYZ(mov)
 423   ^                          ^
 424   | inv(Dmov)                |
 425   |                          |
 426   mov' physvox                |
 427   ^                          |
 428   | inv(Mfsl)                | R
 429   |                          |
 430   ref' physvox                |
 431   ^                          |
 432   | Dref                     |
 433   |                          |
 434   ref voxel <--inv(Tref)--tkRegXYZ(ref)
 435   -------------------------------------------------------------------*/
 436 MATRIX *MRIfsl2TkReg(MRI *ref, MRI *mov, MATRIX *FSLRegMat)
 437 {
 438   MATRIX *RegMat=NULL, *invDmov, *Tmov, *Dref;
 439   MATRIX *invFSLRegMat, *invTref, *Tref;
 440 
 441   /* R = Tmov * inv(Dmov) * inv(Mfsl) * Dref * inv(Tref) */
 442   invDmov = MatrixAlloc(4,4,MATRIX_REAL);
 443   invDmov->rptr[1][1] = 1.0/mov->xsize;
 444   invDmov->rptr[2][2] = 1.0/mov->ysize;
 445   invDmov->rptr[3][3] = 1.0/mov->zsize;
 446   invDmov->rptr[4][4] = 1.0;
 447 
 448   Dref = MatrixAlloc(4,4,MATRIX_REAL);
 449   Dref->rptr[1][1] = ref->xsize;
 450   Dref->rptr[2][2] = ref->ysize;
 451   Dref->rptr[3][3] = ref->zsize;
 452   Dref->rptr[4][4] = 1.0;
 453 
 454   Tmov = MRIxfmCRS2XYZtkreg(mov);
 455   Tref = MRIxfmCRS2XYZtkreg(ref);
 456   invTref = MatrixInverse(Tref,NULL);
 457 
 458   invFSLRegMat = MatrixInverse(FSLRegMat,NULL);
 459 
 460   RegMat = MatrixMultiply(Tmov,invDmov,RegMat);
 461   RegMat = MatrixMultiply(RegMat,invFSLRegMat,RegMat);
 462   RegMat = MatrixMultiply(RegMat,Dref,RegMat);
 463   RegMat = MatrixMultiply(RegMat,invTref,RegMat);
 464 
 465   MatrixFree(&invDmov);
 466   MatrixFree(&FSLRegMat);
 467   MatrixFree(&invFSLRegMat);
 468   MatrixFree(&Tmov);
 469   MatrixFree(&Tref);
 470   MatrixFree(&invTref);
 471   MatrixFree(&Dref);
 472 
 473   return(RegMat);
 474 }
 475 /*-------------------------------------------------------------------
 476   MRItkreg2FSL() - converts tkregister registration matrix to one
 477   compatible with FSL. Note: the FSL matrix is assumed to map from the
 478   mov to the ref whereas the tkreg matrix maps from the ref to the
 479   mov.
 480 
 481   mov voxel<--inv(Tmov)-- tkRegXYZ(mov)
 482   |                           ^
 483   |  Dmov                     |
 484   V                           |
 485   mov' physvox                 |
 486   |                           |
 487   |  Mfsl                     | R
 488   V                           |
 489   ref' physvox                 |
 490   |                           |
 491   |  inv(Dref)                |
 492   V                           |
 493   ref voxel -- Tref  -->  tkRegXYZ(ref)
 494 
 495   -------------------------------------------------------------------*/
 496 MATRIX *MRItkreg2FSL(MRI *ref, MRI *mov, MATRIX *tkRegMat)
 497 {
 498   MATRIX *FSLRegMat=NULL, *Dmov, *Tmov, *invTmov, *Tref, *Dref, *invDref;
 499 
 500   /* R = Tmov * inv(Dmov) * inv(Mfsl) * Dref * inv(Tref) */
 501   /* Mfsl =  inv( Dmov * inv(Tmov) * R * Tref * inv(Dref)) */
 502   Dmov = MatrixAlloc(4,4,MATRIX_REAL);
 503   Dmov->rptr[1][1] = mov->xsize;
 504   Dmov->rptr[2][2] = mov->ysize;
 505   Dmov->rptr[3][3] = mov->zsize;
 506   Dmov->rptr[4][4] = 1.0;
 507 
 508   Dref = MatrixAlloc(4,4,MATRIX_REAL);
 509   Dref->rptr[1][1] = ref->xsize;
 510   Dref->rptr[2][2] = ref->ysize;
 511   Dref->rptr[3][3] = ref->zsize;
 512   Dref->rptr[4][4] = 1.0;
 513   invDref = MatrixInverse(Dref,NULL);
 514 
 515   Tmov = MRIxfmCRS2XYZtkreg(mov);
 516   invTmov = MatrixInverse(Tmov,NULL);
 517   Tref = MRIxfmCRS2XYZtkreg(ref);
 518 
 519   FSLRegMat = MatrixMultiply(Dmov,invTmov,FSLRegMat);
 520   FSLRegMat = MatrixMultiply(FSLRegMat,tkRegMat,FSLRegMat);
 521   FSLRegMat = MatrixMultiply(FSLRegMat,Tref,FSLRegMat);
 522   FSLRegMat = MatrixMultiply(FSLRegMat,invDref,FSLRegMat);
 523   FSLRegMat = MatrixInverse(FSLRegMat,FSLRegMat);
 524 
 525   if(0){
 526     printf("--- Dmov ---------------------\n");
 527     MatrixPrint(stdout,Dmov);
 528     printf("--- Tmov ---------------------\n");
 529     MatrixPrint(stdout,Tmov);
 530     printf("--- R ---------------------\n");
 531     MatrixPrint(stdout,tkRegMat);
 532     printf("--- Tref ---------------------\n");
 533     MatrixPrint(stdout,Tref);
 534     printf("--- Dref ---------------------\n");
 535     MatrixPrint(stdout,Dref);
 536     printf("--- Rfsl ---------------------\n");
 537     MatrixPrint(stdout,FSLRegMat);
 538     printf("--- R (from Rfsl) ------------\n");
 539     tkRegMat = MRIfsl2TkReg(ref,mov,FSLRegMat);
 540     MatrixPrint(stdout,tkRegMat);
 541   }
 542 
 543   MatrixFree(&Dmov);
 544   MatrixFree(&Tmov);
 545   MatrixFree(&invTmov);
 546   MatrixFree(&Tref);
 547   MatrixFree(&Dref);
 548   MatrixFree(&invDref);
 549 
 550   return(FSLRegMat);
 551 }
 552 
 553 /*--------------------------------------------------------------------------
 554   MtxCRS1toCRS0() - generates a matrix that will convert 1-based CRS (as
 555   found in SPM matrices) to 0-based CRS, ie, CRS0 = Q*CRS1 (assuming that
 556   CRS1 has been packed with a 1 in the 4th component.
 557   --------------------------------------------------------------------------*/
 558 MATRIX *MtxCRS1toCRS0(MATRIX *Q)
 559 {
 560   int r,c;
 561 
 562   if(Q == NULL) Q = MatrixAlloc(4,4,MATRIX_REAL);
 563   else{
 564     if(Q->rows != 4 || Q->cols != 4){
 565       printf("ERROR: MtxCRS1toCRS0(): input matrix is not 4x4\n");
 566       return(NULL);
 567     }
 568   }
 569 
 570   for(r=1;r<=4;r++){
 571     for(c=1;c<=4;c++){
 572       if(r==c || c==4) Q->rptr[r][c] = 1.0;
 573       else             Q->rptr[r][c] = 0.0;
 574     }
 575   }
 576 
 577   return(Q);
 578 }
 579 
 580 
 581 /*-------------------------------------------------------------------
 582   MRIgetVoxVal() - returns voxel value as a float regardless of
 583   the underlying data type.
 584   -------------------------------------------------------------------*/
 585 float MRIgetVoxVal(MRI *mri, int c, int r, int s, int f)
 586 {
 587   switch(mri->type){
 588   case MRI_UCHAR: return( (float) MRIseq_vox(mri,c,r,s,f)); break;
 589   case MRI_SHORT: return( (float) MRISseq_vox(mri,c,r,s,f)); break;
 590   case MRI_INT:   return( (float) MRIIseq_vox(mri,c,r,s,f)); break;
 591   case MRI_LONG:  return( (float) MRILseq_vox(mri,c,r,s,f)); break;
 592   case MRI_FLOAT: return( (float) MRIFseq_vox(mri,c,r,s,f)); break;
 593   }
 594   return(-10000000000.9);
 595 }
 596 /*-------------------------------------------------------------------
 597   MRIsetVoxVal() - sets voxel value to that passed as the float
 598   voxval, regardless of the underlying data type. If the underlying
 599   type is integer-based, then it is rounded to the nearest integer. No
 600   attempt is made to prevent overflows.
 601   -------------------------------------------------------------------*/
 602 int MRIsetVoxVal(MRI *mri, int c, int r, int s, int f, float voxval)
 603 {
 604   switch(mri->type){
 605   case MRI_UCHAR: MRIseq_vox(mri,c,r,s,f)  = nint(voxval); break;
 606   case MRI_SHORT: MRISseq_vox(mri,c,r,s,f) = nint(voxval); break;
 607   case MRI_INT:   MRIIseq_vox(mri,c,r,s,f) = nint(voxval); break;
 608   case MRI_LONG:  MRILseq_vox(mri,c,r,s,f) = nint(voxval); break;
 609   case MRI_FLOAT: MRIFseq_vox(mri,c,r,s,f) = voxval;       break;
 610   default: return(1); break;
 611   }
 612   return(0);
 613 }
 614 
 615 /*------------------------------------------------------------------
 616   MRIinterpCode() - returns the numeric interpolation code given the
 617   name of the interpolation method.
 618   -----------------------------------------------------------------*/
 619 int MRIinterpCode(char *InterpString)
 620 {
 621   if(!strncasecmp(InterpString,"nearest",3))
 622     return(SAMPLE_NEAREST);
 623   if(!strncasecmp(InterpString,"trilinear",3))
 624     return(SAMPLE_TRILINEAR);
 625   if(!strncasecmp(InterpString,"sinc",3))
 626     return(SAMPLE_SINC);
 627   if(!strncasecmp(InterpString,"cubic",3))
 628     return(SAMPLE_CUBIC);
 629 
 630   return(-1);
 631 }
 632 
 633 /*------------------------------------------------------------------
 634   MRIinterpString() - returns the the name of the interpolation method
 635   given numeric interpolation code
 636   -----------------------------------------------------------------*/
 637 char * MRIinterpString(int InterpCode)
 638 {
 639   switch (InterpCode)
 640     {
 641     case SAMPLE_NEAREST:   return("nearest"); break ;
 642     case SAMPLE_TRILINEAR: return("trilinear"); break ;
 643     case SAMPLE_SINC:      return("sinc"); break ;
 644     case SAMPLE_CUBIC:     return("cubic"); break ;
 645     }
 646   return(NULL);
 647 }
 648 /*------------------------------------------------------------------
 649   MRIprecisionCode() - returns the numeric code given the
 650   name of the precision. This corresponds to the value of the type
 651   field in the MRI structure.
 652   -----------------------------------------------------------------*/
 653 int MRIprecisionCode(char *PrecisionString)
 654 {
 655   if(!strcasecmp(PrecisionString,"uchar"))
 656     return(MRI_UCHAR);
 657   if(!strcasecmp(PrecisionString,"short"))
 658     return(MRI_SHORT);
 659   if(!strcasecmp(PrecisionString,"int"))
 660     return(MRI_INT);
 661   if(!strcasecmp(PrecisionString,"long"))
 662     return(MRI_LONG);
 663   if(!strcasecmp(PrecisionString,"float"))
 664     return(MRI_FLOAT);
 665 
 666   return(-1);
 667 }
 668 
 669 /*------------------------------------------------------------------
 670   MRIprecisionString() - returns the the name of the precision given
 671   numeric precision code. The code corresponds to the value of the
 672   type field in the MRI structure.
 673   -----------------------------------------------------------------*/
 674 char * MRIprecisionString(int PrecisionCode)
 675 {
 676   switch (PrecisionCode)
 677     {
 678     case MRI_UCHAR: return("uchar"); break ;
 679     case MRI_SHORT: return("short"); break ;
 680     case MRI_INT:   return("int"); break ;
 681     case MRI_LONG:  return("long"); break ;
 682     case MRI_FLOAT: return("float"); break ;
 683     }
 684   return(NULL);
 685 }
 686 
 687 /*-----------------------------------------------------
 688   ------------------------------------------------------*/
 689 int
 690 MRImatch(MRI *mri1, MRI *mri2)
 691 {
 692   return(
 693          (mri1->width == mri2->width) &&
 694          (mri1->height == mri2->height) &&
 695          (mri1->depth == mri2->depth) &&
 696          (mri1->type == mri2->type)
 697          ) ;
 698 }
 699 /*-----------------------------------------------------
 700   ------------------------------------------------------*/
 701 MRI *
 702 MRIscalarMul(MRI *mri_src, MRI *mri_dst, float scalar)
 703 {
 704   int     width, height, depth, x, y, z, frame ;
 705   BUFTYPE *psrc, *pdst ;
 706   float   *pfsrc, *pfdst, dval ;
 707   short   *pssrc, *psdst ;
 708 
 709   width = mri_src->width ;
 710   height = mri_src->height ;
 711   depth = mri_src->depth ;
 712   if (!mri_dst)
 713     mri_dst = MRIclone(mri_src, NULL) ;
 714 
 715   for (frame = 0 ; frame < mri_src->nframes ; frame++)
 716     {
 717       for (z = 0 ; z < depth ; z++)
 718         {
 719           for (y = 0 ; y < height ; y++)
 720             {
 721               switch (mri_src->type)
 722                 {
 723                 case MRI_UCHAR:
 724                   psrc = &MRIseq_vox(mri_src, 0, y, z, frame) ;
 725                   pdst = &MRIseq_vox(mri_dst, 0, y, z, frame) ;
 726                   for (x = 0 ; x < width ; x++)
 727                     {
 728                       dval = *psrc++ * scalar ;
 729                       if (dval < 0)
 730                         dval = 0 ;
 731                       if (dval > 255)
 732                         dval = 255 ;
 733                       *pdst++ = dval ;
 734                     }
 735                   break ;
 736                 case MRI_FLOAT:
 737                   pfsrc = &MRIFseq_vox(mri_src, 0, y, z, frame) ;
 738                   pfdst = &MRIFseq_vox(mri_dst, 0, y, z, frame) ;
 739                   for (x = 0 ; x < width ; x++)
 740                     *pfdst++ = *pfsrc++ * scalar ;
 741                   break ;
 742                 case MRI_SHORT:
 743                   pssrc = &MRISseq_vox(mri_src, 0, y, z, frame) ;
 744                   psdst = &MRISseq_vox(mri_dst, 0, y, z, frame) ;
 745                   for (x = 0 ; x < width ; x++)
 746                     *psdst++ = (short)nint((float)*pssrc++ * scalar) ;
 747                   break ;
 748                 default:
 749                   ErrorReturn
 750                     (NULL,
 751                      (ERROR_UNSUPPORTED,
 752                       "MRIscalarMul: unsupported type %d", mri_src->type)) ;
 753                 }
 754             }
 755         }
 756     }
 757   return(mri_dst) ;
 758 }
 759 /*-----------------------------------------------------
 760   ------------------------------------------------------*/
 761 MRI *
 762 MRIscalarMulFrame(MRI *mri_src, MRI *mri_dst, float scalar, int frame)
 763 {
 764   int     width, height, depth, x, y, z ;
 765   BUFTYPE *psrc, *pdst ;
 766   float   *pfsrc, *pfdst, dval ;
 767   short   *pssrc, *psdst ;
 768 
 769   width = mri_src->width ;
 770   height = mri_src->height ;
 771   depth = mri_src->depth ;
 772   if (!mri_dst)
 773     mri_dst = MRIclone(mri_src, NULL) ;
 774 
 775   for (z = 0 ; z < depth ; z++)
 776     {
 777       for (y = 0 ; y < height ; y++)
 778         {
 779           switch (mri_src->type)
 780             {
 781             case MRI_UCHAR:
 782               psrc = &MRIseq_vox(mri_src, 0, y, z, frame) ;
 783               pdst = &MRIseq_vox(mri_dst, 0, y, z, frame) ;
 784               for (x = 0 ; x < width ; x++)
 785                 {
 786                   dval = *psrc++ * scalar ;
 787                   if (dval < 0)
 788                     dval = 0 ;
 789                   if (dval > 255)
 790                     dval = 255 ;
 791                   *pdst++ = dval ;
 792                 }
 793               break ;
 794             case MRI_FLOAT:
 795               pfsrc = &MRIFseq_vox(mri_src, 0, y, z, frame) ;
 796               pfdst = &MRIFseq_vox(mri_dst, 0, y, z, frame) ;
 797               for (x = 0 ; x < width ; x++)
 798                 *pfdst++ = *pfsrc++ * scalar ;
 799               break ;
 800             case MRI_SHORT:
 801               pssrc = &MRISseq_vox(mri_src, 0, y, z, frame) ;
 802               psdst = &MRISseq_vox(mri_dst, 0, y, z, frame) ;
 803               for (x = 0 ; x < width ; x++)
 804                 *psdst++ = (short)nint((float)*pssrc++ * scalar) ;
 805               break ;
 806             default:
 807               ErrorReturn
 808                 (NULL,
 809                  (ERROR_UNSUPPORTED,
 810                   "MRIscalarMulFrame: unsupported type %d", mri_src->type)) ;
 811             }
 812         }
 813     }
 814   return(mri_dst) ;
 815 }
 816 /*-----------------------------------------------------
 817   ------------------------------------------------------*/
 818 int
 819 MRIvalRange(MRI *mri, float *pmin, float *pmax)
 820 {
 821   int      width, height, depth, x, y, z, frame ;
 822   float    fmin, fmax, *pf, val ;
 823   BUFTYPE  *pb ;
 824 
 825   width = mri->width ;
 826   height = mri->height ;
 827   depth = mri->depth ;
 828 
 829   fmin = 10000.0f ;
 830   fmax = -10000.0f ;
 831   switch (mri->type)
 832     {
 833     case MRI_FLOAT:
 834       for (frame = 0 ; frame < mri->nframes ; frame++)
 835         {
 836           for (z = 0 ; z < depth ; z++)
 837             {
 838               for (y = 0 ; y < height ; y++)
 839                 {
 840                   pf = &MRIFseq_vox(mri, 0, y, z, frame) ;
 841                   for (x = 0 ; x < width ; x++)
 842                     {
 843                       val = *pf++ ;
 844                       if (val < fmin)
 845                         fmin = val ;
 846                       if (val > fmax)
 847                         fmax = val ;
 848                     }
 849                 }
 850             }
 851         }
 852       break ;
 853     case MRI_INT:
 854       for (frame = 0 ; frame < mri->nframes ; frame++)
 855         {
 856           for (z = 0 ; z < depth ; z++)
 857             {
 858               for (y = 0 ; y < height ; y++)
 859                 {
 860                   for (x = 0 ; x < width ; x++)
 861                     {
 862                       val = (float)MRIIseq_vox(mri, x, y, z, frame) ;
 863                       if (val < fmin)
 864                         fmin = val ;
 865                       if (val > fmax)
 866                         fmax = val ;
 867                     }
 868                 }
 869             }
 870         }
 871       break ;
 872     case MRI_SHORT:
 873       for (frame = 0 ; frame < mri->nframes ; frame++)
 874         {
 875           for (z = 0 ; z < depth ; z++)
 876             {
 877               for (y = 0 ; y < height ; y++)
 878                 {
 879                   for (x = 0 ; x < width ; x++)
 880                     {
 881                       val = (float)MRISseq_vox(mri, x, y, z, frame) ;
 882                       if (val < fmin)
 883                         fmin = val ;
 884                       if (val > fmax)
 885                         fmax = val ;
 886                     }
 887                 }
 888             }
 889         }
 890       break ;
 891     case MRI_UCHAR:
 892       for (frame = 0 ; frame < mri->nframes ; frame++)
 893         {
 894           for (z = 0 ; z < depth ; z++)
 895             {
 896               for (y = 0 ; y < height ; y++)
 897                 {
 898                   pb = &MRIseq_vox(mri, 0, y, z, frame) ;
 899                   for (x = 0 ; x < width ; x++)
 900                     {
 901                       val = (float)*pb++ ;
 902                       if (val < fmin)
 903                         fmin = val ;
 904                       if (val > fmax)
 905                         fmax = val ;
 906                     }
 907                 }
 908             }
 909         }
 910       break ;
 911     default:
 912       for (frame = 0 ; frame < mri->nframes ; frame++)
 913         {
 914           for (z = 0 ; z < depth ; z++)
 915             {
 916               for (y = 0 ; y < height ; y++)
 917                 {
 918                   for (x = 0 ; x < width ; x++)
 919                     {
 920                       val = (float)MRIgetVoxVal(mri, x, y, z, frame) ;
 921                       if (val < fmin)
 922                         fmin = val ;
 923                       if (val > fmax)
 924                         fmax = val ;
 925                     }
 926                 }
 927             }
 928         }
 929       break ;
 930     }
 931 
 932   *pmin = fmin ;
 933   *pmax = fmax ;
 934   return(NO_ERROR) ;
 935 }
 936 int
 937 MRInonzeroValRange(MRI *mri, float *pmin, float *pmax)
 938 {
 939   int      width, height, depth, x, y, z, frame ;
 940   float    fmin, fmax, val ;
 941 
 942   width = mri->width ;
 943   height = mri->height ;
 944   depth = mri->depth ;
 945 
 946   fmin = 10000.0f ;
 947   fmax = -10000.0f ;
 948   for (frame = 0 ; frame < mri->nframes ; frame++)
 949     {
 950       for (z = 0 ; z < depth ; z++)
 951         {
 952           for (y = 0 ; y < height ; y++)
 953             {
 954               for (x = 0 ; x < width ; x++)
 955                 {
 956                   val = MRIgetVoxVal(mri, x, y, z, 0) ;
 957                   if (FZERO(val))
 958                     continue ;
 959                   if (val < fmin)
 960                     fmin = val ;
 961                   if (val > fmax)
 962                     fmax = val ;
 963                 }
 964             }
 965         }
 966     }
 967 
 968   *pmin = fmin ;
 969   *pmax = fmax ;
 970   return(NO_ERROR) ;
 971 }
 972 /*-----------------------------------------------------
 973   ------------------------------------------------------*/
 974 int
 975 MRIvalRangeFrame(MRI *mri, float *pmin, float *pmax, int frame)
 976 {
 977   int      width, height, depth, x, y, z ;
 978   float    fmin, fmax, *pf, val ;
 979   BUFTYPE  *pb ;
 980 
 981   width = mri->width ;
 982   height = mri->height ;
 983   depth = mri->depth ;
 984 
 985   fmin = 10000.0f ;
 986   fmax = -10000.0f ;
 987   switch (mri->type)
 988     {
 989     case MRI_FLOAT:
 990       for (z = 0 ; z < depth ; z++)
 991         {
 992           for (y = 0 ; y < height ; y++)
 993             {
 994               pf = &MRIFseq_vox(mri, 0, y, z, frame) ;
 995               for (x = 0 ; x < width ; x++)
 996                 {
 997                   val = *pf++ ;
 998                   if (val < fmin)
 999                     fmin = val ;
1000                   if (val > fmax)
1001                     fmax = val ;
1002                 }
1003             }
1004         }
1005       break ;
1006     case MRI_INT:
1007       for (z = 0 ; z < depth ; z++)
1008         {
1009           for (y = 0 ; y < height ; y++)
1010             {
1011               for (x = 0 ; x < width ; x++)
1012                 {
1013                   val = (float)MRIIseq_vox(mri, x, y, z, frame) ;
1014                   if (val < fmin)
1015                     fmin = val ;
1016                   if (val > fmax)
1017                     fmax = val ;
1018                 }
1019             }
1020         }
1021       break ;
1022     case MRI_SHORT:
1023       for (z = 0 ; z < depth ; z++)
1024         {
1025           for (y = 0 ; y < height ; y++)
1026             {
1027               for (x = 0 ; x < width ; x++)
1028                 {
1029                   val = (float)MRISseq_vox(mri, x, y, z, frame) ;
1030                   if (val < fmin)
1031                     fmin = val ;
1032                   if (val > fmax)
1033                     fmax = val ;
1034                 }
1035             }
1036         }
1037       break ;
1038     case MRI_UCHAR:
1039       for (z = 0 ; z < depth ; z++)
1040         {
1041           for (y = 0 ; y < height ; y++)
1042             {
1043               pb = &MRIseq_vox(mri, 0, y, z, frame) ;
1044               for (x = 0 ; x < width ; x++)
1045                 {
1046                   val = (float)*pb++ ;
1047                   if (val < fmin)
1048                     fmin = val ;
1049                   if (val > fmax)
1050                     fmax = val ;
1051                 }
1052             }
1053         }
1054       break ;
1055     default:
1056       ErrorReturn(ERROR_UNSUPPORTED,
1057                   (ERROR_UNSUPPORTED, "MRIvalRange: unsupported type %d",
1058                    mri->type)) ;
1059     }
1060 
1061   *pmin = fmin ;
1062   *pmax = fmax ;
1063   return(NO_ERROR) ;
1064 }
1065 /*-----------------------------------------------------
1066   ------------------------------------------------------*/
1067 int
1068 MRIvalRangeRegion(MRI *mri, float *pmin, float *pmax, MRI_REGION *region)
1069 {
1070   int      width, height, depth, x, y, z, x0, y0, z0 ;
1071   float    fmin, fmax, *pf, val ;
1072   BUFTYPE  *pb ;
1073 
1074   width = region->x + region->dx ;
1075   if (width > mri->width)
1076     width = mri->width ;
1077   height = region->y + region->dy ;
1078   if (height > mri->height)
1079     height = mri->height ;
1080   depth = region->z + region->dz ;
1081   if (depth > mri->depth)
1082     depth = mri->depth ;
1083   x0 = region->x ;
1084   if (x0 < 0)
1085     x0 = 0 ;
1086   y0 = region->y ;
1087   if (y0 < 0)
1088     y0 = 0 ;
1089   z0 = region->z ;
1090   if (z0 < 0)
1091     z0 = 0 ;
1092 
1093   fmin = 10000.0f ;
1094   fmax = -10000.0f ;
1095   switch (mri->type)
1096     {
1097     default:
1098       for (z = z0 ; z < depth ; z++)
1099         {
1100           for (y = y0 ; y < height ; y++)
1101             {
1102               pf = &MRIFvox(mri, x0, y, z) ;
1103               for (x = x0 ; x < width ; x++)
1104                 {
1105                   val = MRIgetVoxVal(mri, x, y, z, 0) ;
1106                   if (val < fmin)
1107                     fmin = val ;
1108                   if (val > fmax)
1109                     fmax = val ;
1110                 }
1111             }
1112         }
1113       break ;
1114 
1115     case MRI_FLOAT:
1116       for (z = z0 ; z < depth ; z++)
1117         {
1118           for (y = y0 ; y < height ; y++)
1119             {
1120               pf = &MRIFvox(mri, x0, y, z) ;
1121               for (x = x0 ; x < width ; x++)
1122                 {
1123                   val = *pf++ ;
1124                   if (val < fmin)
1125                     fmin = val ;
1126                   if (val > fmax)
1127                     fmax = val ;
1128                 }
1129             }
1130         }
1131       break ;
1132     case MRI_UCHAR:
1133       for (z = z0 ; z < depth ; z++)
1134         {
1135           for (y = y0 ; y < height ; y++)
1136             {
1137               pb = &MRIvox(mri, x0, y, z) ;
1138               for (x = x0 ; x < width ; x++)
1139                 {
1140                   val = (float)*pb++ ;
1141                   if (val < fmin)
1142                     fmin = val ;
1143                   if (val > fmax)
1144                     fmax = val ;
1145                 }
1146             }
1147         }
1148       break ;
1149     }
1150 
1151   *pmin = fmin ;
1152   *pmax = fmax ;
1153   return(NO_ERROR) ;
1154 }
1155 /*-----------------------------------------------------
1156   ------------------------------------------------------*/
1157 MRI_REGION *
1158 MRIclipRegion(MRI *mri, MRI_REGION *reg_src, MRI_REGION *reg_clip)
1159 {
1160   int  x2, y2, z2 ;
1161 
1162   x2 = MIN(mri->width-1, reg_src->x + reg_src->dx - 1) ;
1163   y2 = MIN(mri->height-1, reg_src->y + reg_src->dy - 1) ;
1164   z2 = MIN(mri->depth-1, reg_src->z + reg_src->dz - 1) ;
1165   reg_clip->x = MAX(0, reg_src->x) ;
1166   reg_clip->y = MAX(0, reg_src->y) ;
1167   reg_clip->z = MAX(0, reg_src->z) ;
1168   reg_clip->dx = x2 - reg_clip->x + 1 ;
1169   reg_clip->dy = y2 - reg_clip->y + 1 ;
1170   reg_clip->dz = z2 - reg_clip->z + 1 ;
1171   return(reg_clip) ;
1172 }
1173 /*-----------------------------------------------------
1174   ------------------------------------------------------*/
1175 MRI *
1176 MRIvalScale(MRI *mri_src, MRI *mri_dst, float flo, float fhi)
1177 {
1178   int      width, height, depth, x, y, z ;
1179   float    fmin, fmax, *pf_src, *pf_dst, val, scale ;
1180   short    *ps_src, *ps_dst ;
1181   BUFTYPE  *pb_src, *pb_dst ;
1182 
1183   if (!mri_dst)
1184     mri_dst = MRIclone(mri_src, NULL) ;
1185 
1186   MRIvalRange(mri_src, &fmin, &fmax) ;
1187   scale = (fhi - flo) / (fmax - fmin) ;
1188 
1189   width = mri_src->width ;
1190   height = mri_src->height ;
1191   depth = mri_src->depth ;
1192 
1193   switch (mri_src->type)
1194     {
1195     case MRI_FLOAT:
1196       for (z = 0 ; z < depth ; z++)
1197         {
1198           for (y = 0 ; y < height ; y++)
1199             {
1200               pf_src = &MRIFvox(mri_src, 0, y, z) ;
1201               pf_dst = &MRIFvox(mri_dst, 0, y, z) ;
1202               for (x = 0 ; x < width ; x++)
1203                 {
1204                   val = *pf_src++ ;
1205                   *pf_dst++ = (val - fmin) * scale + flo ;
1206                 }
1207             }
1208         }
1209       break ;
1210     case MRI_SHORT:
1211       for (z = 0 ; z < depth ; z++)
1212         {
1213           for (y = 0 ; y < height ; y++)
1214             {
1215               ps_src = &MRISvox(mri_src, 0, y, z) ;
1216               ps_dst = &MRISvox(mri_dst, 0, y, z) ;
1217               for (x = 0 ; x < width ; x++)
1218                 {
1219                   val = (float)(*ps_src++) ;
1220                   *ps_dst++ = (short)nint((val - fmin) * scale + flo) ;
1221                 }
1222             }
1223         }
1224       break ;
1225     case MRI_UCHAR:
1226       for (z = 0 ; z < depth ; z++)
1227         {
1228           for (y = 0 ; y < height ; y++)
1229             {
1230               pb_src = &MRIvox(mri_src, 0, y, z) ;
1231               pb_dst = &MRIvox(mri_dst, 0, y, z) ;
1232               for (x = 0 ; x < width ; x++)
1233                 {
1234                   val = (float)*pb_src++ ;
1235                   *pb_dst++ = (BUFTYPE)nint((val - fmin) * scale + flo) ;
1236                 }
1237             }
1238         }
1239       break ;
1240     default:
1241       ErrorReturn(mri_dst,
1242                   (ERROR_UNSUPPORTED, "MRIvalScale: unsupported type %d",
1243                    mri_src->type)) ;
1244     }
1245 
1246   return(mri_dst) ;
1247 }
1248 /*-----------------------------------------------------
1249   ------------------------------------------------------*/
1250 MRI *
1251 MRIconfThresh(MRI *mri_src, MRI *mri_probs, MRI *mri_classes, MRI *mri_dst,
1252               float thresh, int min_target, int max_target)
1253 {
1254   int      x, y, z, width, height, depth, class ;
1255   float    *pprobs, prob ;
1256   BUFTYPE  *pclasses, *pdst, *psrc, src ;
1257 
1258   if (!mri_dst)
1259     mri_dst = MRIclone(mri_classes, NULL) ;
1260 
1261   width = mri_classes->width ;
1262   height = mri_classes->height ;
1263   depth = mri_classes->depth ;
1264 
1265   for (z = 0 ; z < depth ; z++)
1266     {
1267       for (y = 0 ; y < height ; y++)
1268         {
1269           pprobs = &MRIFvox(mri_probs, 0, y, z) ;
1270           pclasses = &MRIvox(mri_classes, 0, y, z) ;
1271           pdst = &MRIvox(mri_dst, 0, y, z) ;
1272           psrc = &MRIvox(mri_src, 0, y, z) ;
1273           for (x = 0 ; x < width ; x++)
1274             {
1275               src = *psrc++ ;
1276               prob = *pprobs++ ;
1277               class = (int)*pclasses++ ;
1278               if (prob >= thresh &&
1279                   ((class >= min_target) && (class <= max_target)))
1280                 *pdst++ = src ;
1281               else if ((class >= min_target) && (class <= max_target))
1282                 *pdst++ = 25 ;
1283               else
1284                 *pdst++ = 0 ;
1285             }
1286         }
1287     }
1288   return(mri_dst) ;
1289 }
1290 /*-----------------------------------------------------
1291   ------------------------------------------------------*/
1292 int
1293 MRIboundingBoxNbhd(MRI *mri, int thresh, int wsize,MRI_REGION *box)
1294 {
1295   int      width, height, depth, x, y, z, x1, y1, z1, xi, yi, zi, xk, yk, zk,
1296     whalf, in_brain ;
1297   BUFTYPE  *psrc ;
1298   float    *pfsrc ;
1299   short    *pssrc ;
1300 
1301   whalf = (wsize-1)/2 ;
1302   box->dx = width = mri->width ;
1303   box->dy = height = mri->height ;
1304   box->dz = depth = mri->depth ;
1305 
1306   x1 = y1 = z1 = 0 ;
1307   box->x = width-1 ;
1308   box->y = height-1 ;
1309   box->z = depth-1 ;
1310   switch (mri->type)
1311     {
1312     case MRI_UCHAR:
1313       for (z = 0 ; z < depth ; z++)
1314         {
1315           for (y = 0 ; y < height ; y++)
1316             {
1317               psrc = &MRIvox(mri, 0, y, z) ;
1318               for (x = 0 ; x < width ; x++)
1319                 {
1320                   if (*psrc++ > thresh)
1321                     {
1322                       in_brain = 1 ;
1323                       for (zk = -whalf ; in_brain && zk <= whalf ; zk++)
1324                         {
1325                           zi = mri->zi[z+zk] ;
1326                           for (yk = -whalf ; in_brain && yk <= whalf ; yk++)
1327                             {
1328                               yi = mri->yi[y+yk] ;
1329                               for (xk = -whalf ;
1330                                    in_brain && xk <= whalf ;
1331                                    xk++)
1332                                 {
1333                                   xi = mri->xi[x+xk] ;
1334                                   if (MRIvox(mri, xi, yi, zi) < thresh)
1335                                     in_brain = 0 ;
1336                                 }
1337                             }
1338                         }
1339                       if (in_brain)
1340                         {
1341                           if (x < box->x)
1342                             box->x = x ;
1343                           if (y < box->y)
1344                             box->y = y ;
1345                           if (z < box->z)
1346                             box->z = z ;
1347                           if (x > x1)
1348                             x1 = x ;
1349                           if (y > y1)
1350                             y1 = y ;
1351                           if (z > z1)
1352                             z1 = z ;
1353                         }
1354                     }
1355                 }
1356             }
1357         }
1358       break ;
1359     case MRI_SHORT:
1360       for (z = 0 ; z < depth ; z++)
1361         {
1362           for (y = 0 ; y < height ; y++)
1363             {
1364               pssrc = &MRISvox(mri, 0, y, z) ;
1365               for (x = 0 ; x < width ; x++)
1366                 {
1367                   if (*pssrc++ > thresh)
1368                     {
1369                       in_brain = 1 ;
1370                       for (zk = -whalf ; in_brain && zk <= whalf ; zk++)
1371                         {
1372                           zi = mri->zi[z+zk] ;
1373                           for (yk = -whalf ; in_brain && yk <= whalf ; yk++)
1374                             {
1375                               yi = mri->yi[y+yk] ;
1376                               for (xk = -whalf ;
1377                                    in_brain && xk <= whalf ;
1378                                    xk++)
1379                                 {
1380                                   xi = mri->xi[x+xk] ;
1381                                   if (MRISvox(mri, xi, yi, zi) < thresh)
1382                                     in_brain = 0 ;
1383                                 }
1384                             }
1385                         }
1386                       if (in_brain)
1387                         {
1388                           if (x < box->x)
1389                             box->x = x ;
1390                           if (y < box->y)
1391                             box->y = y ;
1392                           if (z < box->z)
1393                             box->z = z ;
1394                           if (x > x1)
1395                             x1 = x ;
1396                           if (y > y1)
1397                             y1 = y ;
1398                           if (z > z1)
1399                             z1 = z ;
1400                         }
1401                     }
1402                 }
1403             }
1404         }
1405       break ;
1406     case MRI_FLOAT:
1407       for (z = 0 ; z < depth ; z++)
1408         {
1409           for (y = 0 ; y < height ; y++)
1410             {
1411               pfsrc = &MRIFvox(mri, 0, y, z) ;
1412               for (x = 0 ; x < width ; x++)
1413                 {
1414                   if (*pfsrc++ > thresh)
1415                     {
1416                       in_brain = 1 ;
1417                       for (zk = -whalf ; in_brain && zk <= whalf ; zk++)
1418                         {
1419                           zi = mri->zi[z+zk] ;
1420                           for (yk = -whalf ; in_brain && yk <= whalf ; yk++)
1421                             {
1422                               yi = mri->yi[y+yk] ;
1423                               for (xk = -whalf ;
1424                                    in_brain && xk <= whalf ;
1425                                    xk++)
1426                                 {
1427                                   xi = mri->xi[x+xk] ;
1428                                   if (MRIFvox(mri, xi, yi, zi) < thresh)
1429                                     in_brain = 0 ;
1430                                 }
1431                             }
1432                         }
1433                       if (in_brain)
1434                         {
1435                           if (x < box->x)
1436                             box->x = x ;
1437                           if (y < box->y)
1438                             box->y = y ;
1439                           if (z < box->z)
1440                             box->z = z ;
1441                           if (x > x1)
1442                             x1 = x ;
1443                           if (y > y1)
1444                             y1 = y ;
1445                           if (z > z1)
1446                             z1 = z ;
1447                         }
1448                     }
1449                 }
1450             }
1451         }
1452       break ;
1453     default:
1454       ErrorReturn(ERROR_UNSUPPORTED,
1455                   (ERROR_UNSUPPORTED, "MRIboundingBoxNbd: unsupported type %d",
1456                    mri->type)) ;
1457       break ;
1458     }
1459   box->x -= whalf+1 ; box->y -= whalf+1 ; box->z -= whalf+1 ;
1460   x1 += whalf+1 ; y1 += whalf+1 ; z1 += whalf+1 ;
1461   box->x = MAX(0,box->x) ; box->y = MAX(0,box->y) ; box->z = MAX(0,box->z) ;
1462   x1 = MIN(width-1,x1) ; y1 = MIN(height-1, y1) ; z1 = MIN(depth-1, z1) ;
1463   box->dx = x1 - box->x + 1 ;
1464   box->dy = y1 - box->y + 1 ;
1465   box->dz = z1 - box->z + 1 ;
1466   return(NO_ERROR) ;
1467 }
1468 /*-----------------------------------------------------
1469   ------------------------------------------------------*/
1470 #define MIN_DARK 10
1471 
1472 int
1473 MRIfindApproximateSkullBoundingBox(MRI *mri, int thresh,MRI_REGION *box)
1474 {
1475   int      width, height, depth, x, y, z, x1, y1, z1;
1476   int      ndark, max_dark, start, nlight, max_light ;
1477   double   means[3] ;
1478   Real     val ;
1479 
1480   width = mri->width ; height = mri->height ; depth = mri->depth ;
1481 
1482   MRIcenterOfMass(mri, means, thresh) ;
1483 
1484 
1485 #define MAX_LIGHT 30   /* don't let there by 3 cm of
1486                           bright stuff 'outside' of brain */
1487 
1488   /* search for left edge */
1489   nlight = ndark = max_dark = 0 ;
1490   y = nint(means[1]) ; z = nint(means[2]) ;
1491   for (start = x1 = x = nint(means[0]) ; x >= 0 ; x--)
1492     {
1493       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1494 
1495       if (val < thresh)
1496         {
1497           if (!ndark)
1498             start = x ;
1499           ndark++ ;
1500           nlight = 0  ;
1501         }
1502       else
1503         {
1504           if (++nlight > MAX_LIGHT)
1505             max_dark = 0 ;
1506           if (ndark > max_dark)
1507             {
1508               max_dark = ndark ; x1 = start ;
1509             }
1510           ndark = 0 ;
1511         }
1512     }
1513   if (ndark > max_dark)
1514     {
1515       max_dark = ndark ;
1516       x1 = start ;
1517     }
1518   if (max_dark < MIN_DARK)
1519     x1 = 0 ;
1520   box->x = x1 ;
1521 
1522   /* search for right edge */
1523   nlight = ndark = max_dark = 0 ;
1524   y = nint(means[1]) ; z = nint(means[2]) ;
1525   for (start = x1 = x = nint(means[0]) ; x < width ; x++)
1526     {
1527       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1528       if (val < thresh)
1529         {
1530           if (!ndark)
1531             start = x ;
1532           ndark++ ;
1533           nlight = 0 ;
1534         }
1535       else
1536         {
1537           if (++nlight > MAX_LIGHT)
1538             max_dark = 0 ;
1539           if (ndark >= max_dark)
1540             {
1541               max_dark = ndark ; x1 = start ;
1542             }
1543           ndark = 0 ;
1544         }
1545     }
1546   if (ndark > max_dark)
1547     {
1548       max_dark = ndark ;
1549       x1 = start ;
1550     }
1551   if (max_dark < MIN_DARK)
1552     x1 = mri->width-1 ;
1553   box->dx = x1 - box->x + 1 ;
1554 
1555   /* search for superior edge */
1556   nlight = ndark = max_dark = max_light = 0 ;
1557   x = MAX(0,nint(means[0])-20) ; // avoid inter-hemispheric fissure
1558   z = nint(means[2]) ;
1559   for (start = y1 = y = nint(means[1]) ; y >= 0 ; y--)
1560     {
1561       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1562       if (val < thresh)
1563         {
1564           if (nlight > max_light)
1565             max_light = nlight ;
1566           if (!ndark)
1567             start = y ;
1568           ndark++ ;
1569           nlight = 0 ;
1570         }
1571       else
1572         {
1573           if (++nlight > MAX_LIGHT)
1574             max_dark = 0 ;
1575           if (ndark >= max_dark)
1576             {
1577               max_dark = ndark ; y1 = start ;
1578               max_light = 0 ;  // max_light is max in a row light above dark
1579             }
1580           ndark = 0 ;
1581         }
1582     }
1583 
1584   /* if we ended on a string of dark voxels, check two things:
1585      1. the string was longer than the previous longest
1586      2. the strong was longer than 1/2 the previous longest, and there
1587      was an intervening string of light voxels indicated it was still in
1588      brain.
1589   */
1590   if ((ndark > max_dark) || (y < 0 &&
1591                              (ndark > max_dark/2) && max_light > MAX_LIGHT/2))
1592     {
1593       max_dark = ndark ;
1594       y1 = start ;
1595     }
1596   if (max_dark < MIN_DARK)
1597     y1 = 0 ;
1598   box->y = y1 ;
1599 
1600   /* search for inferior edge */
1601   nlight = ndark = max_dark = 0 ;
1602   x = nint(means[0]) ; z = nint(means[2]) ;
1603   for (start = y = y1 = nint(means[1]) ; y < height ; y++)
1604     {
1605       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1606       if (val < thresh)
1607         {
1608           if (!ndark)
1609             start = y ;
1610           ndark++ ;
1611           nlight = 0 ;
1612         }
1613       else
1614         {
1615           if (++nlight > MAX_LIGHT)
1616             max_dark = 0 ;
1617           if (ndark >= max_dark)
1618             {
1619               max_dark = ndark ; y1 = start ;
1620             }
1621           ndark = 0 ;
1622         }
1623     }
1624   if (ndark > max_dark)
1625     {
1626       max_dark = ndark ;
1627       y1 = start ;
1628     }
1629   if (max_dark < MIN_DARK)
1630     y1 = mri->height-1 ;
1631   box->dy = y1 - box->y + 1 ;
1632 
1633   /* search for posterior edge */
1634   nlight = ndark = max_dark = 0 ;
1635   x = nint(means[0]) ; y = nint(means[1]) ;
1636   for (z1 = start = z = nint(means[2]) ; z >= 0 ; z--)
1637     {
1638       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1639       if (val < thresh)
1640         {
1641           if (!ndark)
1642             start = z ;
1643           ndark++ ;
1644           nlight = 0 ;
1645         }
1646       else
1647         {
1648           if (++nlight > MAX_LIGHT)
1649             max_dark = 0 ;
1650           if (ndark >= max_dark)
1651             {
1652               max_dark = ndark ; z1 = start ;
1653             }
1654           ndark = 0 ;
1655         }
1656     }
1657   if (ndark > max_dark)
1658     {
1659       max_dark = ndark ;
1660       z1 = start ;
1661     }
1662   if (max_dark < MIN_DARK)
1663     z1 = 0 ;
1664   box->z = z1 ;
1665 
1666   /* search for anterior edge */
1667   nlight = ndark = max_dark = 0 ;
1668   x = nint(means[0]) ; y = nint(means[1]) ;
1669   for (start = z = nint(means[2]) ; z < depth ; z++)
1670     {
1671       MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
1672       if (val < thresh)
1673         {
1674           if (!ndark)
1675             start = z ;
1676           ndark++ ;
1677           nlight = 0 ;
1678         }
1679       else
1680         {
1681           if (++nlight > MAX_LIGHT)
1682             max_dark = 0 ;
1683           if (ndark >= max_dark)
1684             {
1685               max_dark = ndark ; z1 = start ;
1686             }
1687           ndark = 0 ;
1688         }
1689     }
1690   if (ndark > max_dark)
1691     {
1692       max_dark = ndark ;
1693       z1 = start ;
1694     }
1695   if (max_dark < MIN_DARK)
1696     z1 = mri->depth-1 ;
1697   box->dz = z1 - box->z + 1 ;
1698 
1699   return(NO_ERROR) ;
1700 }
1701 /*-----------------------------------------------------
1702   ------------------------------------------------------*/
1703 int
1704 MRIboundingBox(MRI *mri, int thresh, MRI_REGION *box)
1705 {
1706   int      width, height, depth, x, y, z, x1, y1, z1 ;
1707   BUFTYPE  *psrc ;
1708   float    *pfsrc ;
1709   short    *pssrc ;
1710 
1711   box->dx = width = mri->width ;
1712   box->dy = height = mri->height ;
1713   box->dz = depth = mri->depth ;
1714 
1715   x1 = y1 = z1 = 0 ;
1716   box->x = width-1 ;
1717   box->y = height-1 ;
1718   box->z = depth-1 ;
1719   switch (mri->type)
1720     {
1721     case MRI_UCHAR:
1722       for (z = 0 ; z < depth ; z++)
1723         {
1724           for (y = 0 ; y < height ; y++)
1725             {
1726               psrc = &MRIvox(mri, 0, y, z) ;
1727               for (x = 0 ; x < width ; x++)
1728                 {
1729                   if (*psrc++ > thresh)
1730                     {
1731                       if (x < box->x)
1732                         box->x = x ;
1733                       if (y < box->y)
1734                         box->y = y ;
1735                       if (z < box->z)
1736                         box->z = z ;
1737                       if (x > x1)
1738                         x1 = x ;
1739                       if (y > y1)
1740                         y1 = y ;
1741                       if (z > z1)
1742                         z1 = z ;
1743                     }
1744                 }
1745             }
1746         }
1747       break ;
1748     case MRI_FLOAT:
1749       for (z = 0 ; z < depth ; z++)
1750         {
1751           for (y = 0 ; y < height ; y++)
1752             {
1753               pfsrc = &MRIFvox(mri, 0, y, z) ;
1754               for (x = 0 ; x < width ; x++)
1755                 {
1756                   if (*pfsrc++ > thresh)
1757                     {
1758                       if (x < box->x)
1759                         box->x = x ;
1760                       if (y < box->y)
1761                         box->y = y ;
1762                       if (z < box->z)
1763                         box->z = z ;
1764                       if (x > x1)
1765                         x1 = x ;
1766                       if (y > y1)
1767                         y1 = y ;
1768                       if (z > z1)
1769                         z1 = z ;
1770                     }
1771                 }
1772             }
1773         }
1774       break ;
1775     case MRI_SHORT:
1776       for (z = 0 ; z < depth ; z++)
1777         {
1778           for (y = 0 ; y < height ; y++)
1779             {
1780               pssrc = &MRISvox(mri, 0, y, z) ;
1781               for (x = 0 ; x < width ; x++)
1782                 {
1783                   if (*pssrc++ > thresh)
1784                     {
1785                       if (x < box->x)
1786                         box->x = x ;
1787                       if (y < box->y)
1788                         box->y = y ;
1789                       if (z < box->z)
1790                         box->z = z ;
1791                       if (x > x1)
1792                         x1 = x ;
1793                       if (y > y1)
1794                         y1 = y ;
1795                       if (z > z1)
1796                         z1 = z ;
1797                     }
1798                 }
1799             }
1800         }
1801       break ;
1802     default:
1803       ErrorReturn(ERROR_UNSUPPORTED,
1804                   (ERROR_UNSUPPORTED, "MRIboundingBox: unsupported type %d",
1805                    mri->type)) ;
1806       break ;
1807     }
1808   box->dx = x1 - box->x + 1 ;
1809   box->dy = y1 - box->y + 1 ;
1810   box->dz = z1 - box->z + 1 ;
1811   return(NO_ERROR) ;
1812 }
1813 /*-----------------------------------------------------
1814   ------------------------------------------------------*/
1815 int
1816 MRIcheckSize(MRI *mri_src, MRI *mri_check, int width, int height, int depth)
1817 {
1818   if (!mri_check)
1819     return(0) ;
1820 
1821   if (!width)
1822     width = mri_src->width ;
1823   if (!height)
1824     height = mri_src->height ;
1825   if (!depth)
1826     depth = mri_src->depth ;
1827 
1828   if (width != mri_check->width ||
1829       height != mri_check->height ||
1830       depth != mri_check->depth)
1831     return(0) ;
1832 
1833   return(1) ;
1834 }
1835 /*-----------------------------------------------------
1836   ------------------------------------------------------*/
1837 int
1838 MRItransformRegion(MRI *mri_src, MRI *mri_dst, MRI_REGION *src_region,
1839                    MRI_REGION *dst_region)
1840 {
1841   Real  xw, yw, zw, xt, yt, zt, xv, yv, zv ;
1842 
1843   if (getSliceDirection(mri_src) != getSliceDirection(mri_dst))
1844     ErrorReturn(ERROR_UNSUPPORTED,
1845                 (ERROR_UNSUPPORTED,
1846                  "MRItransformRegion(%s): slice directions must match",
1847                  mri_src->fname)) ;
1848 
1849   if (!mri_src->linear_transform)
1850     ErrorReturn(ERROR_UNSUPPORTED,
1851                 (ERROR_UNSUPPORTED,
1852                  "MRItransformRegion(%s): no transform loaded",
1853                  mri_src->fname)) ;
1854   if (!mri_dst->linear_transform)
1855     ErrorReturn(ERROR_UNSUPPORTED,
1856                 (ERROR_UNSUPPORTED,
1857                  "MRItransformRegion(%s): no transform loaded",
1858                  mri_dst->fname)) ;
1859   /*
1860     The convention  is  that  positive xspace coordinates run
1861     from the patient's  left  side  to  right  side,  positive
1862     yspace  coordinates run from patient posterior to anterior
1863     and positive zspace coordinates run from inferior to superior.
1864   */
1865   switch (getSliceDirection(mri_src))
1866     {
1867     case MRI_CORONAL:
1868       break ;
1869     default:
1870       ErrorReturn
1871         (ERROR_UNSUPPORTED,
1872          (ERROR_UNSUPPORTED,
1873           "MRIregionToTalairachRegion: unsupported slice direction %d",
1874           getSliceDirection(mri_src))) ;
1875     }
1876 
1877   xv = (Real)src_region->x ;
1878   yv = (Real)src_region->y ;
1879   zv = (Real)src_region->z ;
1880   MRIvoxelToWorld(mri_src, xv, yv, zv, &xw, &yw, &zw) ;
1881   transform_point(mri_src->linear_transform, xw, yw, zw, &xt, &yt, &zt) ;
1882   transform_point(mri_dst->inverse_linear_transform, xt, yt, zt, &xw,&yw,&zw);
1883   MRIworldToVoxel(mri_dst, xw, yw, zw, &xv, &yv, &zv) ;
1884   dst_region->x = nint(xv) ;
1885   dst_region->y = nint(yv) ;
1886   dst_region->z = nint(zv) ;
1887 
1888   xv = (Real)(src_region->x + src_region->dx - 1) ;
1889   yv = (Real)(src_region->y + src_region->dy - 1) ;
1890   zv = (Real)(src_region->z + src_region->dz - 1) ;
1891   MRIvoxelToWorld(mri_src, xv, yv, zv, &xw, &yw, &zw) ;
1892   transform_point(mri_src->linear_transform, xw, yw, zw, &xt, &yt, &zt) ;
1893   transform_point(mri_dst->inverse_linear_transform, xt, yt, zt, &xw,&yw,&zw);
1894   MRIworldToVoxel(mri_dst, xw, yw, zw, &xv, &yv, &zv) ;
1895   dst_region->dx = nint(xv - (Real)dst_region->x) + 1 ;
1896   dst_region->dy = nint(yv - (Real)dst_region->y) + 1 ;
1897   dst_region->dz = nint(zv - (Real)dst_region->z) + 1 ;
1898 
1899   return(NO_ERROR) ;
1900 }
1901 /*-----------------------------------------------------
1902   ------------------------------------------------------*/
1903 int
1904 MRIvoxelToVoxel(MRI *mri_src, MRI *mri_dst, Real xv, Real yv, Real zv,
1905                 Real *pxt, Real *pyt, Real *pzt)
1906 {
1907   Real  xw, yw, zw, xt, yt, zt ;
1908 
1909 #if 0
1910   if (!mri_src->linear_transform)
1911     ErrorReturn(ERROR_UNSUPPORTED,
1912                 (ERROR_UNSUPPORTED,
1913                  "MRIvoxelToVoxel(%s): no transform loaded", mri_src->fname));
1914 #endif
1915 
1916   /*
1917     The convention  is  that  positive xspace coordinates run
1918     from the patient's  left  side  to  right  side,  positive
1919     yspace  coordinates run from patient posterior to anterior
1920     and positive zspace coordinates run from inferior to superior.
1921   */
1922   switch (getSliceDirection(mri_src))
1923     {
1924     case MRI_CORONAL:
1925       break ;
1926     default:
1927       ErrorReturn(ERROR_UNSUPPORTED,
1928                   (ERROR_UNSUPPORTED,
1929                    "MRIvoxelToVoxel: unsupported slice direction %d",
1930                    getSliceDirection(mri_src))) ;
1931     }
1932 
1933   if (!mri_src->linear_transform || !mri_dst->inverse_linear_transform)
1934     {
1935       /*
1936          if either doesn't have a transform defined, assume they are in
1937          the same coordinate system.
1938       */
1939       *pxt = xv ; *pyt = yv ; *pzt = zv ;
1940     }
1941   else
1942     {
1943       MRIvoxelToWorld(mri_src, xv, yv, zv, &xw, &yw, &zw) ;
1944       if (mri_src->linear_transform)
1945         transform_point(mri_src->linear_transform, xw, yw, zw, &xt, &yt, &zt) ;
1946       else
1947         {
1948           xt = xw ; yt = yw ; zt = zw ;
1949         }
1950       if (mri_dst->inverse_linear_transform)
1951         transform_point
1952           (mri_dst->inverse_linear_transform, xt,yt,zt,&xw,&yw,&zw);
1953       else
1954         {
1955           xw = xt ; yw = yt ; zw = zt ;
1956         }
1957       MRIworldToVoxel(mri_dst, xw, yw, zw, pxt, pyt, pzt) ;
1958     }
1959 
1960   return(NO_ERROR) ;
1961 }
1962 /*-----------------------------------------------------
1963   ------------------------------------------------------*/
1964 int
1965 MRIvoxelToTalairachVoxel(MRI *mri, Real xv, Real yv, Real zv,
1966                          Real *pxt, Real *pyt, Real *pzt)
1967 {
1968   Real  xw, yw, zw, xt, yt, zt ;
1969 
1970 #if 0
1971   if (!mri->linear_transform)
1972     ErrorReturn(ERROR_UNSUPPORTED,
1973                 (ERROR_UNSUPPORTED,
1974                  "MRIvoxelToTalairachVoxel(%s): no transform loaded",
1975                  mri->fname)) ;
1976 #endif
1977   /*
1978     The convention  is  that  positive xspace coordinates run
1979     from the patient's  left  side  to  right  side,  positive
1980     yspace  coordinates run from patient posterior to anterior
1981     and positive zspace coordinates run from inferior to superior.
1982   */
1983   switch (getSliceDirection(mri))
1984     {
1985     case MRI_CORONAL:
1986       break ;
1987     default:
1988       ErrorReturn(ERROR_UNSUPPORTED,
1989                   (ERROR_UNSUPPORTED,
1990                    "MRIvoxelToTalairachVoxel: unsupported slice direction %d",
1991                    getSliceDirection(mri))) ;
1992     }
1993 
1994   MRIvoxelToWorld(mri, xv, yv, zv, &xw, &yw, &zw) ;
1995   if (mri->linear_transform)
1996     transform_point(mri->linear_transform, xw, yw, zw, &xt, &yt, &zt) ;
1997   else
1998     { xt = xw ; yt = yw ; zt = zw ; }
1999   MRIworldToVoxel(mri, xt, yt, zt, pxt, pyt, pzt) ;
2000 
2001   return(NO_ERROR) ;
2002 }
2003 /*-----------------------------------------------------
2004   ------------------------------------------------------*/
2005 int
2006 MRIvoxelToTalairach(MRI *mri, Real xv, Real yv, Real zv,
2007                     Real *pxt, Real *pyt, Real *pzt)
2008 {
2009   Real  xw, yw, zw ;
2010 
2011 #if 0
2012   if (!mri->linear_transform)
2013     ErrorReturn(ERROR_UNSUPPORTED,
2014                 (ERROR_UNSUPPORTED,
2015                  "MRIvoxelToTalairachVoxel(%s): no transform loaded",
2016                  mri->fname)) ;
2017 #endif
2018 
2019   /*
2020     The convention  is  that  positive xspace coordinates run
2021     from the patient's  left  side  to  right  side,  positive
2022     yspace  coordinates run from patient posterior to anterior
2023     and positive zspace coordinates run from inferior to superior.
2024   */
2025   switch (getSliceDirection(mri))
2026     {
2027     case MRI_CORONAL:
2028       break ;
2029     default:
2030       ErrorReturn(ERROR_UNSUPPORTED,
2031                   (ERROR_UNSUPPORTED,
2032                    "MRIvoxelToTalairachVoxel: unsupported slice direction %d",
2033                    getSliceDirection(mri))) ;
2034     }
2035 
2036   MRIvoxelToWorld(mri, xv, yv, zv, &xw, &yw, &zw) ;
2037   if (mri->linear_transform)
2038     transform_point(mri->linear_transform, xw, yw, zw, pxt, pyt, pzt) ;
2039   else
2040     { *pxt = xw ; *pyt = yw ; *pzt = zw ; }
2041 
2042   return(NO_ERROR) ;
2043 }
2044 /*-----------------------------------------------------
2045   ------------------------------------------------------*/
2046 int
2047 MRItalairachToVoxel(MRI *mri, Real xt, Real yt, Real zt,
2048                     Real *pxv, Real *pyv, Real *pzv)
2049 {
2050   Real  xw, yw, zw ;
2051 
2052 #if 0
2053   if (!mri->inverse_linear_transform)
2054     ErrorReturn(ERROR_UNSUPPORTED,
2055                 (ERROR_UNSUPPORTED,
2056                  "MRItalairachToVoxel(%s): no transform loaded",
2057                  mri->fname)) ;
2058 #endif
2059 
2060   /*
2061     The convention  is  that  positive xspace coordinates run
2062     from the patient's  left  side  to  right  side,  positive
2063     yspace  coordinates run from patient posterior to anterior
2064     and positive zspace coordinates run from inferior to superior.
2065   */
2066   switch (getSliceDirection(mri))
2067     {
2068     case MRI_CORONAL:
2069       break ;
2070     default:
2071       ErrorReturn(ERROR_UNSUPPORTED,
2072                   (ERROR_UNSUPPORTED,
2073                    "MRIvoxelToTalairachVoxel: unsupported slice direction %d",
2074                    getSliceDirection(mri))) ;
2075     }
2076 
2077   if (mri->inverse_linear_transform)
2078     transform_point(mri->inverse_linear_transform, xt, yt, zt, &xw, &yw,&zw);
2079   else
2080     { xw = xt ; yw = yt ; zw = zt ; }
2081   MRIworldToVoxel(mri, xw, yw, zw, pxv, pyv, pzv) ;
2082 
2083   return(NO_ERROR) ;
2084 }
2085 /*-----------------------------------------------------
2086   ------------------------------------------------------*/
2087 int
2088 MRItalairachVoxelToVoxel(MRI *mri, Real xtv, Real ytv, Real ztv,
2089                          Real *pxv, Real *pyv, Real *pzv)
2090 {
2091   Real  xw, yw, zw, xt, yt, zt ;
2092 
2093 #if 0
2094   if (!mri->inverse_linear_transform)
2095     ErrorReturn(ERROR_UNSUPPORTED,
2096                 (ERROR_UNSUPPORTED,
2097                  "MRItalairachVoxelToVoxel(%s): no transform loaded",
2098                  mri->fname)) ;
2099 #endif
2100 
2101   /*
2102     The convention  is  that  positive xspace coordinates run
2103     from the patient's  left  side  to  right  side,  positive
2104     yspace  coordinates run from patient posterior to anterior
2105     and positive zspace coordinates run from inferior to superior.
2106   */
2107   switch (getSliceDirection(mri))
2108     {
2109     case MRI_CORONAL:
2110       break ;
2111     default:
2112       ErrorReturn(ERROR_UNSUPPORTED,
2113                   (ERROR_UNSUPPORTED,
2114                    "MRIvoxelToTalairachVoxel: unsupported slice direction %d",
2115                    getSliceDirection(mri))) ;
2116     }
2117 
2118   MRIvoxelToWorld(mri, xtv, ytv, ztv, &xt, &yt, &zt) ;
2119   if (mri->inverse_linear_transform)
2120     transform_point(mri->inverse_linear_transform, xt, yt, zt, &xw, &yw,&zw);
2121   else
2122     { xw = xt ; yw = yt ; zw = zt ; }
2123   MRIworldToVoxel(mri, xw, yw, zw, pxv, pyv, pzv) ;
2124 
2125   return(NO_ERROR) ;
2126 }
2127 /*-----------------------------------------------------
2128   ------------------------------------------------------*/
2129 int
2130 MRItalairachVoxelToWorld(MRI *mri, Real xtv, Real ytv, Real ztv,
2131                          Real *pxw, Real *pyw, Real *pzw)
2132 {
2133   Real  xw, yw, zw, xt, yt, zt ;
2134 
2135 #if 0
2136   if (!mri->inverse_linear_transform)
2137     ErrorReturn(ERROR_UNSUPPORTED,
2138                 (ERROR_UNSUPPORTED,
2139                  "MRItalairachVoxelToVoxel(%s): no transform loaded",
2140                  mri->fname)) ;
2141 #endif
2142 
2143   /*
2144     The convention  is  that  positive xspace coordinates run
2145     from the patient's  left  side  to  right  side,  positive
2146     yspace  coordinates run from patient posterior to anterior
2147     and positive zspace coordinates run from inferior to superior.
2148   */
2149   switch (getSliceDirection(mri))
2150     {
2151     case MRI_CORONAL:
2152       break ;
2153     default:
2154       ErrorReturn(ERROR_UNSUPPORTED,
2155                   (ERROR_UNSUPPORTED,
2156                    "MRIvoxelToTalairachVoxel: unsupported slice direction %d",
2157                    getSliceDirection(mri))) ;
2158     }
2159 
2160   MRIvoxelToWorld(mri, xtv, ytv, ztv, &xt, &yt, &zt) ;
2161   if (mri->inverse_linear_transform)
2162     transform_point(mri->inverse_linear_transform, xt, yt, zt, &xw, &yw,&zw);
2163   else
2164     { xw = xt ; yw = yt ; zw = zt ; }
2165   *pxw = xw ; *pyw = yw ; *pzw = zw ;
2166 
2167   return(NO_ERROR) ;
2168 }
2169 /*-----------------------------------------------------
2170   ------------------------------------------------------*/
2171 #define V4_LOAD(v, x, y, z, r)  (VECTOR_ELT(v,1)=x, VECTOR_ELT(v,2)=y, \
$^
                                  VECTOR_ELT(v,3)=z, VECTOR_ELT(v,4)=r) ;
2172 
2173 int
2174 MRIvoxelToWorld(MRI *mri, Real xv, Real yv, Real zv,
2175                 Real *pxw, Real *pyw, Real *pzw)
2176 {
2177   VECTOR *vw, *vv;
2178   MATRIX *RfromI;
2179 
2180   // if the transform is not cached yet, then
2181   if (!mri->i_to_r__)
2182     mri->i_to_r__ = extract_i_to_r(mri);
2183   if (!mri->r_to_i__)
2184     mri->r_to_i__ = extract_r_to_i(mri);
2185 
2186   RfromI = mri->i_to_r__; // extract_i_to_r(mri);
2187 
2188   vv = VectorAlloc(4, MATRIX_REAL) ;
2189   V4_LOAD(vv, xv, yv, zv, 1.) ;
2190   vw = MatrixMultiply(RfromI, vv, NULL) ;
2191   *pxw = V3_X(vw);
2192   *pyw = V3_Y(vw);
2193   *pzw = V3_Z(vw);
2194 
2195   // MatrixFree(&RfromI);
2196   VectorFree(&vv);
2197   VectorFree(&vw);
2198   return(NO_ERROR) ;
2199 }
2200 /*-----------------------------------------------------
2201   ------------------------------------------------------*/
2202 int
2203 MRIworldToTalairachVoxel(MRI *mri, Real xw, Real yw, Real zw,
2204                          Real *pxv, Real *pyv, Real *pzv)
2205 {
2206   Real  xt, yt, zt ;
2207 
2208   transform_point(mri->linear_transform, xw, yw, zw, &xt, &yt, &zt) ;
2209   MRIworldToVoxel(mri, xt, yt, zt, pxv, pyv, pzv) ;
2210   return(NO_ERROR) ;
2211 }
2212 /*-----------------------------------------------------
2213   ------------------------------------------------------*/
2214 int MRIworldToVoxelIndex(MRI *mri, Real xw, Real yw, Real zw,
2215                          int *pxv, int *pyv, int *pzv)
2216 {
2217   Real xv, yv, zv;
2218   MRIworldToVoxel(mri, xw, yw, zw, &xv, &yv, &zv);
2219   *pxv = (int) xv;
2220   *pyv = (int) yv;
2221   *pzv = (int) zv;
2222 
2223   /*
2224     switch (getSliceDirection(mri))
2225     {
2226     case MRI_CORONAL:
2227     #if 0
2228     *pxv = ((Real)mri->xend - xw) / mri->xsize ;
2229     *pzv = (yw - (Real)mri->zstart) / mri->zsize ;
2230     *pyv = (-zw - (Real)mri->ystart) / mri->ysize ;
2231     #else
2232     trans_SetBounds ( mri->xstart, mri->xend, mri->ystart, mri->yend,
2233     mri->zstart, mri->zend );
2234     trans_SetResolution ( mri->xsize, mri->ysize, mri->zsize );
2235     trans_RASToVoxelIndex(xw, yw, zw, pxv, pyv, pzv) ;
2236     #endif
2237     break ;
2238     default:
2239     ErrorReturn(ERROR_UNSUPPORTED,
2240     (ERROR_UNSUPPORTED,
2241     "MRIworldToVoxel: unsupported slice direction %d",
2242     getSliceDirection(mri))) ;
2243     break ;
2244     }
2245   */
2246   return(NO_ERROR) ;
2247 }
2248 
2249 /*
2250    int MRIvoxelToSurfaceRAS and int MRIsurfaceRASToVoxel
2251 
2252    get the surfaceRAS values from original voxel
2253    Note that this is different from MRIvoxelToWorld().
2254 
2255    Note that currently MATRIX uses float** to store data.
2256    Going around the circle of transform causes error accumulation quickly.
2257    I noticed that 7 x 10^(-6) is very common.
2258 
2259 */
2260 
2261 MATRIX *surfaceRASFromVoxel_(MRI *mri)
2262 {
2263   // Derivation (avoid expensive matrix calculation)
2264   //
2265   //     orig -----(rasFromVoxel)------> RAS (scanner RAS or physical RAS)
2266   //       |                              |
2267   //(conformedVoxelFromVoxel)            (1)
2268   //       |                              |
2269   //       V                              V
2270   //conformed--(RASfromConformed)----->RAS (scanner RAS or physical RAS)
2271   //       |                              |
2272   //      (1)                      (surfaceRASFromRAS)
2273   //       |                              |
2274   //       V                              V
2275   //conformed-(surfaceRASFromConformed)->surfaceRAS
2276   //
2277   //
2278   // where  RASFromConformed = [ -1  0  0   s1 ] where s1 = c_r + 128
2279   //                           [  0  0  1   s2 ]       s2 = c_a - 128
2280   //                           [  0 -1  0   s3 ]       s3 = c_s + 128
2281   //                           [  0  0  0    1 ]
2282   //
2283   //  surfaceRASFromConformed= [ -1  0  0  128 ]
2284   //                           [  0  0  1 -128 ]
2285   //                           [  0 -1  0  128 ]
2286   //                           [  0  0  0   1  ]
2287   // Therefore
2288   //        surfaceRASFromRAS= [  1  0  0  -c_r]  just a translation matrix
2289   //                           [  0  1  0  -c_a]
2290   //                           [  0  0  1  -c_s]
2291   //                           [  0  0  0    1 ]
2292   //
2293   //  surfaceRASFromVoxel = surfaceRASFromRAS (x) rasFromVoxel
2294   //
2295   //  i.e.       applying translation on RASFromVoxel
2296   //
2297   // This means tha
2298   //
2299   //  if RASFromVoxel = (  X  | T ), then
2300   //                    (  0  | 1 )
2301   //
2302   //    surfaceRASFromVoxel =  ( 1 | -C) * ( X | T ) = ( X | T - C )
2303   //                           ( 0 | 1 )   ( 0 | 1 )   ( 0 |   1   )
2304   //
2305   MATRIX *rasFromVoxel;
2306   MATRIX *sRASFromVoxel;
2307   double m14, m24, m34;
2308 
2309   // if the transform is not cached yet, then
2310   if (!mri->i_to_r__)
2311     mri->i_to_r__ = extract_i_to_r(mri);
2312   if (!mri->r_to_i__)
2313     mri->r_to_i__ = extract_r_to_i(mri);
2314 
2315   rasFromVoxel = mri->i_to_r__; // extract_i_to_r(mri);
2316 
2317   sRASFromVoxel = MatrixCopy(rasFromVoxel, NULL);
2318   // MatrixFree(&rasFromVoxel);
2319   // modify
2320   m14 = *MATRIX_RELT(sRASFromVoxel, 1,4);
2321   *MATRIX_RELT(sRASFromVoxel, 1,4) = m14 - mri->c_r;
2322   m24 = *MATRIX_RELT(sRASFromVoxel, 2,4);
2323   *MATRIX_RELT(sRASFromVoxel, 2,4) = m24 - mri->c_a;
2324   m34 = *MATRIX_RELT(sRASFromVoxel, 3,4);
2325   *MATRIX_RELT(sRASFromVoxel, 3,4) = m34 - mri->c_s;
2326 
2327   return sRASFromVoxel;
2328 }
2329 
2330 MATRIX *surfaceRASFromRAS_(MRI *mri)
2331 {
2332   MATRIX *sRASFromRAS;
2333   sRASFromRAS = MatrixAlloc(4, 4, MATRIX_REAL);
2334   MatrixIdentity(4, sRASFromRAS);
2335   *MATRIX_RELT(sRASFromRAS, 1,4) = - mri->c_r;
2336   *MATRIX_RELT(sRASFromRAS, 2,4) = - mri->c_a;
2337   *MATRIX_RELT(sRASFromRAS, 3,4) = - mri->c_s;
2338   return sRASFromRAS;
2339 }
2340 
2341 MATRIX *RASFromSurfaceRAS_(MRI *mri)
2342 {
2343   MATRIX *RASFromSRAS;
2344   RASFromSRAS = MatrixAlloc(4, 4, MATRIX_REAL);
2345   MatrixIdentity(4, RASFromSRAS);
2346   *MATRIX_RELT(RASFromSRAS, 1,4) = mri->c_r;
2347   *MATRIX_RELT(RASFromSRAS, 2,4) = mri->c_a;
2348   *MATRIX_RELT(RASFromSRAS, 3,4) = mri->c_s;
2349   return RASFromSRAS;
2350 }
2351 
2352 int MRIRASToSurfaceRAS(MRI *mri, Real xr, Real yr, Real zr,
2353                        Real *xsr, Real *ysr, Real *zsr)
2354 {
2355   MATRIX *surfaceRASFromRAS=0;
2356   VECTOR *v, *sr;
2357   v = VectorAlloc(4, MATRIX_REAL);
2358   V4_LOAD(v, xr, yr, zr, 1.);
2359   surfaceRASFromRAS = surfaceRASFromRAS_(mri);
2360   sr = MatrixMultiply(surfaceRASFromRAS, v, NULL);
2361   *xsr = V3_X(sr);
2362   *ysr = V3_Y(sr);
2363   *zsr = V3_Z(sr);
2364   MatrixFree(&surfaceRASFromRAS);
2365   VectorFree(&v);
2366   VectorFree(&sr);
2367   return (NO_ERROR);
2368 }
2369 
2370 int MRIsurfaceRASToRAS(MRI *mri, Real xsr, Real ysr, Real zsr,
2371                        Real *xr, Real *yr, Real *zr)
2372 {
2373   MATRIX *RASFromSurfaceRAS=0;
2374   VECTOR *v, *r;
2375   v = VectorAlloc(4, MATRIX_REAL);
2376   V4_LOAD(v, xsr, ysr, zsr, 1.);
2377   RASFromSurfaceRAS = RASFromSurfaceRAS_(mri);
2378   r = MatrixMultiply(RASFromSurfaceRAS, v, NULL);
2379   *xr = V3_X(r);
2380   *yr = V3_Y(r);
2381   *zr = V3_Z(r);
2382   MatrixFree(&RASFromSurfaceRAS);
2383   VectorFree(&v);
2384   VectorFree(&r);
2385   return (NO_ERROR);
2386 }
2387 
2388 
2389 int MRIvoxelToSurfaceRAS(MRI *mri, Real xv, Real yv, Real zv,
2390                          Real *xs, Real *ys, Real *zs)
2391 {
2392   MATRIX *sRASFromVoxel;
2393   VECTOR *vv, *sr;
2394 
2395   sRASFromVoxel = surfaceRASFromVoxel_(mri);
2396   // calculate the surface ras value
2397   vv = VectorAlloc(4, MATRIX_REAL);
2398   V4_LOAD(vv, xv, yv, zv, 1.);
2399   sr = MatrixMultiply(sRASFromVoxel, vv, NULL);
2400   *xs = V3_X(sr);
2401   *ys = V3_Y(sr);
2402   *zs = V3_Z(sr);
2403 
2404   MatrixFree(&sRASFromVoxel);
2405   VectorFree(&vv);
2406   VectorFree(&sr);
2407 
2408   return (NO_ERROR);
2409 }
2410 
2411 MATRIX *voxelFromSurfaceRAS_(MRI *mri)
2412 {
2413   MATRIX *voxelFromSRAS =0;
2414   //////////////////////////////////////////////////////////////////
2415   // it turned out that this can be done easily without taking inverse
2416   // Note the surfaceRASFromVoxel_ is given by
2417   //
2418   //       ( X | T - C)
2419   //       ( 0 |  1   )
2420   // Note, however, that we define C by
2421   //
2422   //      (  X | T ) (S/2) = ( C )  where S = (w/2, h/2, d/2)^t
2423   //      (  0 | 1 ) ( 1 )   ( 1 )
2424   // Thus
2425   //        X*S/2 + T = C
2426   // or
2427   //        T - C = - X*S/2
2428   // or
2429   //     surfaceRASFromVoxel = ( X | - X*S/2 )
2430   //                           ( 0 |    1    )
2431   // whose inverse is given by
2432   //
2433   //     voxelFromSurfaceRAS = ( X ^(-1)| S/2 )
2434   //                           (  0     |  1  )
2435   //
2436   // since
2437   //           ( X^(-1)  S/2) ( X  -X*S/2) = ( 1   S/2 - S/2) = 1
2438   //           (   0      1 ) ( 0     1  )   ( 0      1     )
2439   //
2440   // thus get r_to_i__ and set translation part to S/2 is the quickest way
2441   /////////////////////////////////////////////////////////////////////
2442   // if the transform is not cached yet, then
2443   if (!mri->i_to_r__)
2444     mri->i_to_r__ = extract_i_to_r(mri);
2445   if (!mri->r_to_i__)
2446     mri->r_to_i__ = extract_r_to_i(mri);
2447 
2448   voxelFromSRAS = MatrixCopy(mri->r_to_i__, NULL);
2449   // modify translation part
2450   *MATRIX_RELT(voxelFromSRAS, 1,4) = mri->width/2;
2451   *MATRIX_RELT(voxelFromSRAS, 2,4) = mri->height/2;
2452   *MATRIX_RELT(voxelFromSRAS, 3,4) = mri->depth/2;
2453 
2454 #if 0
2455   // no more expensive inverse
2456   MATRIX *sRASFromVoxel = surfaceRASFromVoxel_(mri);
2457   MATRIX *voxelFromSRAS = MatrixInverse(sRASFromVoxel, NULL);
2458   MatrixFree(&sRASFromVoxel) ;
2459 #endif
2460 
2461   return voxelFromSRAS;
2462 }
2463 
2464 /* extract the RASToVoxel Matrix */
2465 MATRIX *GetSurfaceRASToVoxelMatrix(MRI *mri){
2466   return voxelFromSurfaceRAS_(mri);
2467 }
2468 
2469 int MRIsurfaceRASToVoxel(MRI *mri, Real xr, Real yr, Real zr,
2470                          Real *xv, Real *yv, Real *zv)
2471 {
2472   MATRIX *voxelFromSRAS;
2473   static VECTOR *sr = NULL, *vv = NULL;
2474 
2475   voxelFromSRAS=voxelFromSurfaceRAS_(mri);
2476   if (sr == NULL)
2477     sr = VectorAlloc(4, MATRIX_REAL);
2478   V4_LOAD(sr, xr, yr, zr, 1.);
2479   vv = MatrixMultiply(voxelFromSRAS, sr, vv);
2480   *xv = V3_X(vv);
2481   *yv = V3_Y(vv);
2482   *zv = V3_Z(vv);
2483 
2484   MatrixFree(&voxelFromSRAS);
2485   //  VectorFree(&sr);
2486   //  VectorFree(&vv);
2487 
2488   return (NO_ERROR);
2489 }
2490 
2491 /*------------------------------------------------------*/
2492 int
2493 MRIworldToVoxel(MRI *mri, Real xw, Real yw, Real zw,
2494                 Real *pxv, Real *pyv, Real *pzv)
2495 {
2496   VECTOR *vv, *vw;
2497   MATRIX *IfromR;
2498 
2499   // if transform is not cached yet, then
2500   if (!mri->r_to_i__)
2501     mri->r_to_i__ = extract_r_to_i(mri);
2502   if (!mri->i_to_r__)
2503     mri->i_to_r__ = extract_i_to_r(mri);
2504 
2505   IfromR = mri->r_to_i__;
2506   vw = VectorAlloc(4, MATRIX_REAL) ;
2507   V4_LOAD(vw, xw, yw, zw, 1.) ;
2508   vv = MatrixMultiply(IfromR, vw, NULL) ;
2509   *pxv = V3_X(vv);
2510   *pyv = V3_Y(vv);
2511   *pzv = V3_Z(vv);
2512 
2513   VectorFree(&vv);
2514   VectorFree(&vw);
2515 
2516   return(NO_ERROR) ;
2517 }
2518 /*-----------------------------------------------------
2519   ------------------------------------------------------*/
2520 int
2521 MRIinitHeader(MRI *mri)
2522 {
2523   mri->ptype = 2 ;
2524 
2525   /* most of these are in mm */
2526   mri->imnr0 = 1 ;
2527   mri->imnr1 = mri->depth ;
2528   mri->fov = mri->width ;
2529   mri->thick = 1.0 ;
2530   mri->xsize = 1.0 ;
2531   mri->ysize = 1.0 ;
2532   mri->zsize = 1.0 ;
2533   mri->ps = 1 ;
2534   mri->xstart = -mri->width/2.0 ;
2535   mri->xend = mri->width/2.0 ;
2536   mri->ystart = -mri->height/2.0 ;
2537   mri->yend = mri->height/2.0 ;
2538   mri->zstart = -mri->depth/2.0 ;
2539   mri->zend = mri->depth/2 ;
2540   // coronal
2541   mri->x_r = -1;
2542   mri->x_a = 0.;
2543   mri->x_s = 0.;
2544   //
2545   mri->y_r = 0.;
2546   mri->y_a = 0.;
2547   mri->y_s = -1;
2548   //
2549   mri->z_r = 0.;
2550   mri->z_a = 1;
2551   mri->z_s = 0.;
2552   //
2553   mri->c_r = mri->c_a = mri->c_s = 0.0;
2554 
2555   mri->ras_good_flag = 1;
2556   mri->gdf_image_stem[0] = '\0';
2557   mri->tag_data = NULL;
2558   mri->tag_data_size = 0;
2559 
2560   if (!mri->i_to_r__)
2561     mri->i_to_r__ = extract_i_to_r(mri);
2562   if (!mri->r_to_i__)
2563     mri->r_to_i__ = extract_r_to_i(mri);
2564 
2565   return(NO_ERROR) ;
2566 }
2567 
2568 /**
2569  * MRIreInitCache
2570  *
2571  * @param mri MRI* whose header information was modified
2572  *
2573  * @return NO_ERROR
2574  */
2575 int MRIreInitCache(MRI *mri)
2576 {
2577   if (mri->i_to_r__)
2578     {
2579       MatrixFree(&mri->i_to_r__);
2580       mri->i_to_r__ = 0;
2581     }
2582   if (mri->r_to_i__)
2583     {
2584       MatrixFree(&mri->r_to_i__);
2585       mri->r_to_i__ = 0;
2586     }
2587   mri->i_to_r__ = extract_i_to_r(mri);
2588   mri->r_to_i__ = extract_r_to_i(mri);
2589 
2590   return (NO_ERROR);
2591 }
2592 
2593 /*-----------------------------------------------------
2594   Parameters:
2595 
2596   Returns value:
2597 
2598   Description
2599   change the direction of slices
2600   ------------------------------------------------------*/
2601 MRI *
2602 MRIextract(MRI *mri_src, MRI *mri_dst, int x0, int y0, int z0,
2603            int dx, int dy, int dz)
2604 {
2605   return(MRIextractInto(mri_src, mri_dst, x0, y0, z0, dx, dy, dz, 0, 0, 0)) ;
2606 }
2607 /*-----------------------------------------------------
2608   Parameters:
2609 
2610   Returns value:
2611 
2612   Description
2613   Extract a cubic region of an MR image and return it to the caller
2614   ------------------------------------------------------*/
2615 MRI *
2616 MRIextractRegion(MRI *mri_src, MRI *mri_dst, MRI_REGION *region)
2617 {
2618   return(MRIextractInto(mri_src, mri_dst, region->x, region->y, region->z,
2619                         region->dx, region->dy, region->dz, 0, 0, 0)) ;
2620 }
2621 /*-----------------------------------------------------
2622   Parameters:
2623 
2624   Returns value:
2625 
2626   Description
2627   Extract a cubic region of an MR image and return it to the caller
2628   ------------------------------------------------------*/
2629 MRI *
2630 MRIextractIntoRegion(MRI *mri_src, MRI *mri_dst, int x0, int y0, int z0,
2631                      MRI_REGION *region)
2632 {
2633   return(MRIextractInto(mri_src, mri_dst, x0, y0, z0, region->dx, region->dy,
2634                         region->dz, region->x, region->y, region->z)) ;
2635 }
2636 /*-----------------------------------------------------
2637   Parameters:
2638 
2639   Returns value:
2640 
2641   Description
2642   Extract a cubic region of an MR image and return it to the caller
2643   ------------------------------------------------------*/
2644 MRI *
2645 MRIextractInto(MRI *mri_src, MRI *mri_dst, int x0, int y0, int z0,
2646                int dx, int dy, int dz, int x1, int y1, int z1)
2647 {
2648   int  width, height, depth, ys, zs, yd, zd, bytes, frame, xsize,ysize,zsize,
2649     dst_alloced = 0 ;
2650   Real c_r, c_a, c_s;
2651 
2652   width = mri_src->width ;
2653   depth = mri_src->depth ;
2654   height = mri_src->height ;
2655 
2656   if (z0 >= depth || y0 >= height || x0 >= width)
2657     ErrorReturn(NULL,
2658                 (ERROR_BADPARM,
2659                  "MRIextractInto: bad src location (%d, %d, %d)", x0,y0,z0));
2660   // validation
2661   if (x0 < 0)
2662     x0 = 0 ;
2663   if (y0 < 0)
2664     y0 = 0 ;
2665   if (z0 < 0)
2666     z0 = 0 ;
2667   if (x0+dx > width)
2668     dx = (width - x0) ;
2669   if (y0+dy > height)
2670     dy = (height - y0) ;
2671   if (z0+dz > depth)
2672     dz = (depth - z0) ;
2673   if (x1 < 0)
2674     x1 = 0 ;
2675   if (y1 < 0)
2676     y1 = 0 ;
2677   if (z1 < 0)
2678     z1 = 0 ;
2679 
2680   if (!mri_dst)
2681     {
2682       mri_dst = MRIallocSequence(dx, dy, dz, mri_src->type, mri_src->nframes) ;
2683       MRIcopyHeader(mri_src, mri_dst) ;
2684       mri_dst->imnr0 = z0 + mri_src->imnr0 - z1 ;
2685       mri_dst->imnr1 = mri_dst->imnr0 + dz - 1 ;
2686       dst_alloced = 1 ;
2687     }
2688 
2689   if (mri_src->type != mri_dst->type)
2690     {
2691       MRIfree(&mri_dst) ;
2692       ErrorReturn(NULL,
2693                   (ERROR_BADPARM,
2694                    "MRIextractInto: src and dst types must match"));
2695     }
2696 
2697   if (x1+dx > mri_dst->width)
2698     dx = (mri_dst->width - x1) ;
2699   if (y1+dy > mri_dst->height)
2700     dy = (mri_dst->height - y1) ;
2701   if (z1+dz > mri_dst->depth)
2702     dz = (mri_dst->depth - z1) ;
2703 
2704   xsize = mri_src->xsize ;
2705   ysize = mri_src->ysize ;
2706   zsize = mri_src->zsize ;
2707 
2708   if (dst_alloced)
2709     {
2710       mri_dst->xstart += x0*xsize ;
2711       mri_dst->xend = mri_dst->xstart + dx*xsize ;
2712       mri_dst->ystart += y0*ysize ;
2713       mri_dst->yend = mri_dst->ystart + dy*ysize ;
2714       mri_dst->zstart += z0*zsize ;
2715       mri_dst->zend = mri_dst->zstart + dz*zsize  ;
2716     }
2717 
2718   bytes = dx ;
2719   switch (mri_src->type)
2720     {
2721     case MRI_FLOAT:
2722       bytes *= sizeof(float) ;
2723       break ;
2724     case MRI_LONG:
2725       bytes *= sizeof(long) ;
2726       break ;
2727     case MRI_INT:
2728       bytes *= sizeof(int) ;
2729       break ;
2730     case MRI_SHORT:
2731       bytes *= sizeof(short) ;
2732       break ;
2733     default:
2734       break ;
2735     }
2736 
2737   for (frame = 0 ; frame < mri_src->nframes ; frame++)
2738     {
2739       for (zd = z1, zs = z0 ; zs < z0+dz ; zs++, zd++)
2740         {
2741           for (yd = y1, ys = y0 ; ys < y0+dy ; ys++, yd++)
2742             {
2743               switch (mri_src->type)
2744                 {
2745                 case MRI_UCHAR:
2746                   memcpy(&MRIseq_vox(mri_dst, x1, yd, zd,frame),
2747                          &MRIseq_vox(mri_src,x0,ys,zs,frame), bytes);
2748                   break ;
2749                 case MRI_FLOAT:
2750                   memcpy(&MRIFseq_vox(mri_dst, x1, yd, zd,frame),
2751                          &MRIFseq_vox(mri_src,x0,ys,zs,frame), bytes);
2752                   break ;
2753                 case MRI_SHORT:
2754                   memcpy(&MRISseq_vox(mri_dst, x1, yd, zd,frame),
2755                          &MRISseq_vox(mri_src,x0,ys,zs,frame), bytes);
2756                   break ;
2757                 case MRI_LONG:
2758                   memcpy(&MRILseq_vox(mri_dst, x1, yd, zd,frame),
2759                          &MRILseq_vox(mri_src,x0,ys,zs,frame), bytes);
2760                   break ;
2761                 }
2762             }
2763         }
2764     }
2765   // calculate c_ras
2766   MRIcalcCRASforExtractedVolume
2767     (mri_src, mri_dst, x0, y0, z0, x1, y1, z1, &c_r, &c_a, &c_s);
2768   mri_dst->c_r = c_r;
2769   mri_dst->c_a = c_a;
2770   mri_dst->c_s = c_s;
2771   // initialize cached transform
2772   MRIreInitCache(mri_dst);
2773 
2774   return(mri_dst) ;
2775 }
2776 /*-----------------------------------------------------
2777   Parameters:
2778 
2779   Returns value:
2780 
2781   Description
2782   change the direction of slices
2783   ------------------------------------------------------*/
2784 MRI *
2785 MRIreslice(MRI *mri_src, MRI *mri_dst, int slice_direction)
2786 {
2787   int     width, height, depth, x1, x2, x3 ;
2788   BUFTYPE *psrc, val, *pdst ;
2789 
2790   int src_slice_direction = getSliceDirection(mri_src);
2791   if (slice_direction == src_slice_direction)
2792     {
2793       mri_dst = MRIcopy(mri_src, NULL) ;
2794       return(mri_dst) ;
2795     }
2796 
2797   width = mri_src->width ;
2798   height = mri_src->height ;
2799   depth = mri_src->depth ;
2800 
2801 
2802   if ((src_slice_direction == MRI_SAGITTAL &&
2803        slice_direction == MRI_CORONAL) ||
2804       (src_slice_direction == MRI_CORONAL &&
2805        slice_direction == MRI_SAGITTAL))
2806     {
2807       /*
2808         coronal images are back to front of the head, thus the depth axis
2809         points from the nose to the back of the head, with x from neck to
2810         crown, and y from ear to ear.
2811       */
2812       /* x1 --> x3
2813          x2 --> x2
2814          x3 --> x1
2815       */
2816       if (!mri_dst)
2817         {
2818           mri_dst = MRIalloc(depth, height, width, mri_src->type) ;
2819           MRIcopyHeader(mri_src, mri_dst) ;
2820         }
2821       else if (depth != mri_dst->width || height != mri_dst->height ||
2822                width != mri_dst->depth)
2823         {
2824           ErrorReturn(NULL,
2825                       (ERROR_BADPARM,
2826                        "MRIreslice: invalid destination size (%d, %d, %d)",
2827                        mri_dst->width, mri_dst->height, mri_dst->depth)) ;
2828         }
2829 
2830       for (x3 = 0 ; x3 < depth ; x3++)
2831         {
2832           for (x2 = 0 ; x2 < height ; x2++)
2833             {
2834               psrc = &MRIvox(mri_src, 0, x2, x3) ;
2835               for (x1 = 0 ; x1 < width ; x1++)
2836                 {
2837                   /* swap so in place transformations are possible */
2838                   mri_dst->slices[x1][x2][x3] = *psrc++ ;
2839                 }
2840             }
2841         }
2842     }
2843   else
2844     if ((src_slice_direction == MRI_HORIZONTAL &&
2845          slice_direction == MRI_CORONAL) ||
2846         (src_slice_direction == MRI_CORONAL &&
2847          slice_direction == MRI_HORIZONTAL))
2848       {
2849         /*
2850           horizontal images are top to bottom of the head, thus the depth axis
2851           points from the top of the head to the neck, with x from ear to ear
2852           and y from nose to back of head.
2853         */
2854         /* x3 --> x2
2855            x2 --> x3
2856            x1 --> x1
2857         */
2858         if (!mri_dst)
2859           {
2860             mri_dst = MRIalloc(width, depth, height, mri_src->type) ;
2861             MRIcopyHeader(mri_src, mri_dst) ;
2862           }
2863         else if (depth != mri_dst->height || height != mri_dst->depth ||
2864                  width != mri_dst->width)
2865           ErrorReturn(NULL,
2866                       (ERROR_BADPARM,
2867                        "MRIreslice: invalid destination size (%d, %d, %d)",
2868                        mri_dst->width, mri_dst->height, mri_dst->depth)) ;
2869 
2870         for (x3 = 0 ; x3 < depth ; x3++)
2871           {
2872             for (x2 = 0 ; x2 < height ; x2++)
2873               {
2874                 psrc = &MRIvox(mri_src, 0, x2, x3) ;
2875                 pdst = &MRIvox(mri_dst, 0, x3, x2) ;
2876                 for (x1 = 0 ; x1 < width ; x1++)
2877                   {
2878                     /* swap so in place transformations are possible */
2879                     *pdst++ = *psrc++ ;
2880                   }
2881               }
2882           }
2883       }
2884     else
2885       if ((src_slice_direction == MRI_SAGITTAL &&
2886            slice_direction == MRI_HORIZONTAL))
2887         {
2888           /*
2889             horizontal images are top to bottom of the head,
2890             thus the depth axis
2891             points from the top of the head to the neck, with x from
2892             ear to ear
2893             and y from nose to back of head.
2894           */
2895           /* x3 --> x2
2896              x1 --> x3
2897              x2 --> x1
2898           */
2899           if (!mri_dst)
2900             {
2901               mri_dst = MRIalloc(width, depth, height, mri_src->type) ;
2902               MRIcopyHeader(mri_src, mri_dst) ;
2903             }
2904           else if (depth != mri_dst->height || height != mri_dst->depth ||
2905                    width != mri_dst->width)
2906             ErrorReturn(NULL,
2907                         (ERROR_BADPARM,
2908                          "MRIreslice: invalid destination size (%d, %d, %d)",
2909                          mri_dst->width, mri_dst->height, mri_dst->depth)) ;
2910 
2911           for (x3 = 0 ; x3 < depth ; x3++)
2912             {
2913               for (x2 = 0 ; x2 < height ; x2++)
2914                 {
2915                   psrc = &MRIvox(mri_src, 0, x2, x3) ;
2916                   for (x1 = 0 ; x1 < width ; x1++)
2917                     {
2918                       /* swap so in place transformations are possible */
2919                       mri_dst->slices[x2][x1][x3] = *psrc++ ;
2920                     }
2921                 }
2922             }
2923         }
2924       else
2925         if (src_slice_direction == MRI_HORIZONTAL &&
2926             slice_direction == MRI_SAGITTAL)
2927           {
2928             /*
2929               horizontal images are top to bottom of the head,
2930               thus the depth axis
2931               points from the top of the head to the neck,
2932               with x from ear to ear
2933               and y from nose to back of head.
2934             */
2935             /* x2 --> x3
2936                x3 --> x1
2937                x1 --> x2
2938             */
2939             if (!mri_dst)
2940               {
2941                 mri_dst = MRIalloc(width, depth, height, mri_src->type) ;
2942                 MRIcopyHeader(mri_src, mri_dst) ;
2943               }
2944             else if (depth != mri_dst->height || height != mri_dst->depth ||
2945                      width != mri_dst->width)
2946               ErrorReturn(NULL,
2947                           (ERROR_BADPARM,
2948                            "MRIreslice: invalid destination size (%d, %d, %d)",
2949                            mri_dst->width, mri_dst->height, mri_dst->depth)) ;
2950 
2951             for (x3 = 0 ; x3 < depth ; x3++)
2952               {
2953                 for (x2 = 0 ; x2 < height ; x2++)
2954                   {
2955                     psrc = &MRIvox(mri_src, 0, x2, x3) ;
2956                     for (x1 = 0 ; x1 < width ; x1++)
2957                       {
2958                         /* swap so in place transformations are possible */
2959                         mri_dst->slices[x1][x3][x2] = *psrc++ ;
2960                       }
2961                   }
2962               }
2963           }
2964         else
2965           switch (src_slice_direction)
2966             {
2967             default:
2968               MRIfree(&mri_dst) ;
2969               ErrorReturn(NULL,
2970                           (ERROR_BADPARM,
2971                            "MRIreslice: mri_src unknown slice direction %d",
2972                            src_slice_direction)) ;
2973               break ;
2974             case MRI_CORONAL:
2975               /*
2976                 coronal images are back to front of the head,
2977                 thus the depth axis
2978                 points from the nose to the back of the head,
2979                 with x from neck to
2980                 crown, and y from ear to ear.
2981               */
2982               switch (slice_direction)
2983                 {
2984                 case MRI_SAGITTAL:
2985                   /* x1 --> x3
2986                      x2 --> x2
2987                      x3 --> x1
2988                   */
2989                   if (!mri_dst)
2990                     {
2991                       mri_dst = MRIalloc(depth, height, width, mri_src->type) ;
2992                       MRIcopyHeader(mri_src, mri_dst) ;
2993                     }
2994                   else if (depth != mri_dst->width ||
2995                            height != mri_dst->height ||
2996                            width != mri_dst->depth)
2997                     ErrorReturn
2998                       (NULL,
2999                        (ERROR_BADPARM,
3000                         "MRIreslice: invalid destination size (%d, %d, %d)",
3001                         mri_dst->width, mri_dst->height, mri_dst->depth)) ;
3002 
3003                   for (x3 = 0 ; x3 < depth ; x3++)
3004                     {
3005                       for (x2 = 0 ; x2 < height ; x2++)
3006                         {
3007                           psrc = &MRIvox(mri_src, 0, x2, x3) ;
3008                           for (x1 = 0 ; x1 < width ; x1++)
3009                             {
3010                               /* swap so in place transformations
3011                                  are possible */
3012                               val = *psrc++ ;
3013 #if 0
3014                               mri_dst->slices[x3][x2][x1] =
3015                                 mri_src->slices[x1][x2][x3] ;
3016 #endif
3017                               mri_dst->slices[x1][x2][x3] = val ;
3018                             }
3019                         }
3020                     }
3021                   break ;
3022                 case MRI_HORIZONTAL:
3023                   break ;
3024                 }
3025               break ;
3026             case MRI_SAGITTAL:
3027               /*
3028                 sagittal images are slices in
3029                 the plane of the nose, with depth going
3030                 from ear to ear.
3031               */
3032               break ;
3033             }
3034   setDirectionCosine(mri_dst, slice_direction);
3035   mri_dst->ras_good_flag = 0;
3036   return(mri_dst) ;
3037 }
3038 /*-----------------------------------------------------
3039   Parameters:
3040 
3041   Returns value:
3042 
3043   Description
3044   Set an MRI intensity values to 0
3045   ------------------------------------------------------*/
3046 int
3047 MRIclear(MRI *mri)
3048 {
3049   int   width, depth, height, bytes, y, z, frame, nframes ;
3050 
3051   width = mri->width ;
3052   height = mri->height ;
3053   depth = mri->depth ;
3054   nframes = mri->nframes ;
3055   bytes = width ;
3056 
3057   switch (mri->type)
3058     {
3059     case MRI_UCHAR:
3060       bytes *= sizeof(unsigned char) ;
3061       break ;
3062     case MRI_BITMAP:
3063       bytes /= 8 ;
3064       break ;
3065     case MRI_FLOAT:
3066       bytes *= sizeof(float) ;
3067       break ;
3068     case MRI_LONG:
3069       bytes *= sizeof(long) ;
3070       break ;
3071     case MRI_INT:
3072       bytes *= sizeof(int) ;
3073       break ;
3074     case MRI_SHORT:
3075       bytes *= sizeof(short) ;
3076       break ;
3077     default:
3078       ErrorReturn(ERROR_UNSUPPORTED,
3079                   (ERROR_UNSUPPORTED,
3080                    "MRIclear: unsupported input type %d", mri->type)) ;
3081       break ;
3082     }
3083 
3084   for (frame = 0 ; frame < nframes ; frame++)
3085     {
3086       for (z = 0 ; z < depth ; z++)
3087         {
3088           for (y = 0 ; y < height ; y++)
3089             memset(mri->slices[z+frame*depth][y], 0, bytes) ;
3090         }
3091     }
3092 
3093   return(NO_ERROR) ;
3094 }
3095 /*-----------------------------------------------------
3096   Parameters:
3097 
3098   Returns value:
3099 
3100   Description
3101   find the principle components of a (binary) MRI. The
3102   eigenvectors are the columns of the matrix mEvectors, the
3103   eigenvalues are returned in the array evalues and the means
3104   in means (these last two must be three elements long.)
3105   ------------------------------------------------------*/
3106 int
3107 MRIcenterOfMass(MRI *mri,double *means, BUFTYPE threshold)
3108 {
3109   int     width, height, depth, x, y, z ;
3110   long    npoints ;
3111   double  mx, my, mz, weight ;
3112   Real    val ;
3113 
3114   width = mri->width ;
3115   height = mri->height ;
3116   depth = mri->depth ;
3117 
3118   mx = my = mz = weight = 0.0f ; npoints = 0L ;
3119 
3120   for (z = 0 ; z < depth ; z++)
3121     {
3122       for (y = 0 ; y < height ; y++)
3123         {
3124           for (x = 0 ; x < width ; x++)
3125             {
3126               MRIsampleVolumeType(mri, x,  y, z, &val, SAMPLE_NEAREST) ;
3127               if (val > threshold)
3128                 {
3129                   weight += val ;
3130                   mx += (float)x*val ;
3131                   my += (float)y*val ;
3132                   mz += (float)z*val ;
3133                   npoints++ ;
3134                 }
3135             }
3136         }
3137     }
3138 
3139   if (weight > 0.0)
3140     {
3141       mx /= weight ;
3142       my /= weight  ;
3143       mz /= weight ;
3144       means[0] = mx ;
3145       means[1] = my ;
3146       means[2] = mz ;
3147     }
3148   else
3149     means[0] = means[1] = means[2] = 0.0f ;
3150 
3151 
3152   return(NO_ERROR) ;
3153 }
3154 /*-----------------------------------------------------
3155   Parameters:
3156 
3157   Returns value:
3158 
3159   Description
3160   find the principle components of a (binary) MRI. The
3161   eigenvectors are the columns of the matrix mEvectors, the
3162   eigenvalues are returned in the array evalues and the means
3163   in means (these last two must be three elements long.)
3164   ------------------------------------------------------*/
3165 int
3166 MRIprincipleComponents(MRI *mri, MATRIX *mEvectors, float *evalues,
3167                        double *means, BUFTYPE threshold)
3168 {
3169   int     width, height, depth, x, y, z ;
3170   BUFTYPE *psrc, val ;
3171   long    npoints ;
3172   MATRIX  *mCov, *mX, *mXT, *mTmp ;
3173   double  mx, my, mz, weight ;
3174 
3175   if (mri->type != MRI_UCHAR)
3176     ErrorReturn(ERROR_UNSUPPORTED,
3177                 (ERROR_UNSUPPORTED,
3178                  "MRIprincipleComponents: unsupported input type %d",
3179                  mri->type)) ;
3180 
3181   width = mri->width ;
3182   height = mri->height ;
3183   depth = mri->depth ;
3184 
3185   mx = my = mz = weight = 0.0f ; npoints = 0L ;
3186 
3187   for (z = 0 ; z < depth ; z++)
3188     {
3189       for (y = 0 ; y < height ; y++)
3190         {
3191           psrc = &MRIvox(mri, 0, y, z) ;
3192           for (x = 0 ; x < width ; x++)
3193             {
3194               val = *psrc++ ;
3195               if (val > threshold)
3196                 {
3197                   weight += val ;
3198                   mx += (float)x*val ;
3199                   my += (float)y*val ;
3200                   mz += (float)z*val ;
3201                   npoints++ ;
3202                 }
3203             }
3204         }
3205     }
3206 
3207   if (weight > 0.0)
3208     {
3209       mx /= weight ;
3210       my /= weight  ;
3211       mz /= weight ;
3212       means[0] = mx ;
3213       means[1] = my ;
3214       means[2] = mz ;
3215     }
3216   else
3217     means[0] = means[1] = means[2] = 0.0f ;
3218 
3219   mX = MatrixAlloc(3, 1, MATRIX_REAL) ;     /* zero-mean coordinate vector */
3220   mXT = NULL ;                              /* transpose of above */
3221   mTmp = MatrixAlloc(3, 3, MATRIX_REAL) ;   /* tmp matrix for covariance */
3222   mCov = MatrixAlloc(3, 3, MATRIX_REAL) ;   /* covariance matrix */
3223 
3224   for (z = 0 ; z < depth ; z++)
3225     {
3226       for (y = 0 ; y < height ; y++)
3227         {
3228           psrc = &MRIvox(mri, 0, y, z) ;
3229           for (x = 0 ; x < width ; x++)
3230             {
3231               val = *psrc++ ;
3232               if (val > threshold)
3233                 {
3234                   mX->rptr[1][1] = (x - (int)mx)*val ;
3235                   mX->rptr[2][1] = (y - (int)my)*val ;
3236                   mX->rptr[3][1] = (z - (int)mz)*val ;
3237                   mXT = MatrixTranspose(mX, mXT) ;
3238                   mTmp = MatrixMultiply(mX, mXT, mTmp) ;
3239                   mCov = MatrixAdd(mTmp, mCov, mCov) ;
3240                 }
3241             }
3242         }
3243     }
3244 
3245   if (weight > 0)
3246     MatrixScalarMul(mCov, 1.0f/weight, mCov) ;
3247 
3248   MatrixEigenSystem(mCov, evalues, mEvectors) ;
3249 
3250   return(NO_ERROR) ;
3251 }
3252 /*-----------------------------------------------------
3253   Parameters:
3254 
3255   Returns value:
3256 
3257   Description
3258   find the principle components of a (binary) MRI. The
3259   eigenvectors are the columns of the matrix mEvectors, the
3260   eigenvalues are returned in the array evalues and the means
3261   in means (these last two must be three elements long.)
3262   ------------------------------------------------------*/
3263 int
3264 MRIbinaryPrincipleComponents(MRI *mri, MATRIX *mEvectors, float *evalues,
3265                              double *means, BUFTYPE threshold)
3266 {
3267   int     width, height, depth, x, y, z ;
3268   BUFTYPE *psrc, val ;
3269   long    npoints ;
3270   MATRIX  *mCov, *mX, *mXT, *mTmp ;
3271   double  mx, my, mz, weight ;
3272 
3273   if (mri->type != MRI_UCHAR)
3274     ErrorReturn(ERROR_UNSUPPORTED,
3275                 (ERROR_UNSUPPORTED,
3276                  "MRIprincipleComponents: unsupported input type %d",
3277                  mri->type)) ;
3278 
3279   width = mri->width ;
3280   height = mri->height ;
3281   depth = mri->depth ;
3282 
3283   mx = my = mz = weight = 0.0f ; npoints = 0L ;
3284 
3285   for (z = 0 ; z < depth ; z++)
3286     {
3287       for (y = 0 ; y < height ; y++)
3288         {
3289           psrc = &MRIvox(mri, 0, y, z) ;
3290           for (x = 0 ; x < width ; x++)
3291             {
3292               val = *psrc++ ;
3293               if (val > threshold)
3294                 {
3295                   weight++ ;
3296                   mx += (float)x ;
3297                   my += (float)y ;
3298                   mz += (float)z ;
3299                   npoints++ ;
3300                 }
3301             }
3302         }
3303     }
3304 
3305   if (weight > 0.0)
3306     {
3307       mx /= weight ;
3308       my /= weight  ;
3309       mz /= weight ;
3310       means[0] = mx ;
3311       means[1] = my ;
3312       means[2] = mz ;
3313     }
3314   else
3315     means[0] = means[1] = means[2] = 0.0f ;
3316 
3317   mX = MatrixAlloc(3, 1, MATRIX_REAL) ;     /* zero-mean coordinate vector */
3318   mXT = NULL ;                              /* transpose of above */
3319   mTmp = MatrixAlloc(3, 3, MATRIX_REAL) ;   /* tmp matrix for covariance */
3320   mCov = MatrixAlloc(3, 3, MATRIX_REAL) ;   /* covariance matrix */
3321 
3322   for (z = 0 ; z < depth ; z++)
3323     {
3324       for (y = 0 ; y < height ; y++)
3325         {
3326           psrc = &MRIvox(mri, 0, y, z) ;
3327           for (x = 0 ; x < width ; x++)
3328             {
3329               val = *psrc++ ;
3330               if (val > threshold)
3331                 {
3332                   mX->rptr[1][1] = (x - (int)mx) ;
3333                   mX->rptr[2][1] = (y - (int)my) ;
3334                   mX->rptr[3][1] = (z - (int)mz) ;
3335                   mXT = MatrixTranspose(mX, mXT) ;
3336                   mTmp = MatrixMultiply(mX, mXT, mTmp) ;
3337                   mCov = MatrixAdd(mTmp, mCov, mCov) ;
3338                 }
3339             }
3340         }
3341     }
3342 
3343   if (weight > 0)
3344     MatrixScalarMul(mCov, 1.0f/weight, mCov) ;
3345 
3346   MatrixEigenSystem(mCov, evalues, mEvectors) ;
3347 
3348   return(NO_ERROR) ;
3349 }
3350 /*-----------------------------------------------------
3351   Parameters:
3352 
3353   Returns value:
3354 
3355   Description
3356   threshold an MRI.
3357   ------------------------------------------------------*/
3358 MRI *
3359 MRIthresholdRangeInto(MRI *mri_src,MRI *mri_dst,BUFTYPE low_val,BUFTYPE hi_val)
3360 {
3361   int     width, height, depth, x, y, z ;
3362   BUFTYPE *psrc, *pdst, val ;
3363   float   *pfsrc, *pfdst, fval ;
3364 
3365   if (!mri_dst)
3366     mri_dst = MRIclone(mri_src, NULL) ;
3367 
3368   width = mri_src->width ;
3369   height = mri_src->height ;
3370   depth = mri_src->depth ;
3371 
3372   switch (mri_src->type)
3373     {
3374     case MRI_UCHAR:
3375       for (z = 0 ; z < depth ; z++)
3376         {
3377           for (y = 0 ; y < height ; y++)
3378             {
3379               psrc = &MRIvox(mri_src, 0, y, z) ;
3380               pdst = &MRIvox(mri_dst, 0, y, z) ;
3381               for (x = 0 ; x < width ; x++, pdst++)
3382                 {
3383                   val = *psrc++ ;
3384                   if (val >=  low_val && val <= hi_val)
3385                     *pdst = val ;
3386                   else
3387                     *pdst = 0 ;
3388                 }
3389             }
3390         }
3391       break ;
3392     case MRI_FLOAT:
3393       for (z = 0 ; z < depth ; z++)
3394         {
3395           for (y = 0 ; y < height ; y++)
3396             {
3397               pfsrc = &MRIFvox(mri_src, 0, y, z) ;
3398               pfdst = &MRIFvox(mri_dst, 0, y, z) ;
3399               for (x = 0 ; x < width ; x++, pdst++)
3400                 {
3401                   fval = *pfsrc++ ;
3402                   if (fval >=  low_val && fval <= hi_val)
3403                     *pfdst = fval ;
3404                   else
3405                     *pfdst = 0 ;
3406                 }
3407             }
3408         }
3409       break ;
3410     default:
3411       ErrorReturn(mri_dst,
3412                   (ERROR_UNSUPPORTED,
3413                    "MRIthresholdRangeInto: unsupported type %d",
3414                    mri_src->type)) ;
3415       break ;
3416     }
3417   return(mri_dst) ;
3418 }
3419 /*-----------------------------------------------------
3420   Parameters:
3421 
3422   Returns value:
3423 
3424   Description
3425   threshold an MRI.
3426   ------------------------------------------------------*/
3427 MRI *
3428 MRIthreshold(MRI *mri_src, MRI *mri_dst, float threshold)
3429 {
3430   int     width, height, depth, x, y, z ;
3431   float   val ;
3432 
3433   if (!mri_dst)
3434     mri_dst = MRIclone(mri_src, NULL) ;
3435 
3436   width = mri_src->width ;
3437   height = mri_src->height ;
3438   depth = mri_src->depth ;
3439 
3440   for (z = 0 ; z < depth ; z++)
3441     {
3442       for (y = 0 ; y < height ; y++)
3443         {
3444           for (x = 0 ; x < width ; x++)
3445             {
3446               val = MRIgetVoxVal(mri_src, x, y, z, 0) ;
3447               if (val < threshold)
3448                 val = 0 ;
3449               MRIsetVoxVal(mri_dst, x, y, z, 0, val) ;
3450             }
3451         }
3452     }
3453 
3454   return(mri_dst) ;
3455 }
3456 /*-----------------------------------------------------
3457   Parameters:
3458 
3459   Returns value:
3460 
3461   Description
3462   threshold an MRI.
3463   ------------------------------------------------------*/
3464 MRI  *
3465 MRIinvertContrast(MRI *mri_src, MRI *mri_dst, float threshold)
3466 {
3467   int     width, height, depth, x, y, z ;
3468   BUFTYPE *psrc, *pdst, val ;
3469 
3470   if (!mri_dst)
3471     mri_dst = MRIclone(mri_src, NULL) ;
3472 
3473   width = mri_src->width ;
3474   height = mri_src->height ;
3475   depth = mri_src->depth ;
3476 
3477   for (z = 0 ; z < depth ; z++)
3478     {
3479       for (y = 0 ; y < height ; y++)
3480         {
3481           psrc = &MRIvox(mri_src, 0, y, z) ;
3482           pdst = &MRIvox(mri_dst, 0, y, z) ;
3483           for (x = 0 ; x < width ; x++)
3484             {
3485               val = *psrc++ ;
3486               if (val > threshold)
3487                 val = 255 - val ;
3488               *pdst++ = val ;
3489             }
3490         }
3491     }
3492 
3493   return(mri_dst) ;
3494 }
3495 /*-----------------------------------------------------
3496   Parameters:
3497 
3498   Returns value:
3499 
3500   Description
3501   threshold an MRI.
3502   ------------------------------------------------------*/
3503 MRI *
3504 MRIbinarize(MRI *mri_src, MRI *mri_dst, BUFTYPE threshold, BUFTYPE low_val,
3505             BUFTYPE hi_val)
3506 {
3507   int     width, height, depth, x, y, z, f ;
3508   Real    val ;
3509 
3510   if (!mri_dst)
3511     mri_dst = MRIclone(mri_src, NULL) ;
3512 
3513   width = mri_src->width ;
3514   height = mri_src->height ;
3515   depth = mri_src->depth ;
3516 
3517   for (f = 0 ; f < mri_src->nframes ; f++)
3518     {
3519       for (z = 0 ; z < depth ; z++)
3520         {
3521           for (y = 0 ; y < height ; y++)
3522             {
3523               for (x = 0 ; x < width ; x++)
3524                 {
3525                   MRIsampleVolumeFrameType
3526                     (mri_src, x, y, z, f, SAMPLE_NEAREST, &val) ;
3527                   if (val < threshold)
3528                     val = low_val ;
3529                   else
3530                     val = hi_val ;
3531                   MRIsetVoxVal(mri_dst, x, y, z, f, val) ;
3532                 }
3533             }
3534         }
3535     }
3536 
3537   return(mri_dst) ;
3538 }
3539 
3540 /*-----------------------------------------------------*/
3541 MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst)
3542 {
3543   int     nframes, width, height, depth, x, y, z, f, s ;
3544   float v1, v2, v=0.0;
3545   BUFTYPE *p1=NULL, *p2=NULL, *pdst=NULL ;
3546   short *ps1=NULL, *ps2=NULL, *psdst=NULL ;
3547   int   *pi1=NULL, *pi2=NULL, *pidst=NULL ;
3548   long  *pl1=NULL, *pl2=NULL, *pldst=NULL ;
3549   float *pf1=NULL, *pf2=NULL, *pfdst=NULL ;
3550 
3551   width = mri1->width ;
3552   height = mri1->height ;
3553   depth = mri1->depth ;
3554   nframes = mri1->nframes ;
3555   if(nframes == 0) nframes = 1;
3556 
3557   if (!mri_dst){
3558     mri_dst = MRIallocSequence(width, height, depth, mri1->type,nframes) ;
3559     MRIcopyHeader(mri1, mri_dst) ;
3560   }
3561 
3562   if(mri1->type != mri2->type){
3563     /* Generic but slow */
3564     for (f = 0 ; f < nframes ; f++){
3565       for (z = 0 ; z < depth ; z++){
3566         for (y = 0 ; y < height ; y++){
3567           for (x = 0 ; x < width ; x++){
3568             v1 = MRIgetVoxVal(mri1,x,y,z,f);
3569             v2 = MRIgetVoxVal(mri2,x,y,z,f);
3570             v = v1-v2;
3571             MRIsetVoxVal(mri_dst,x,y,z,f,v);
3572           }
3573         }
3574       }
3575     }
3576     return(mri_dst) ;
3577   }
3578 
3579 
3580   s = 0;
3581   for (f = 0 ; f < nframes ; f++){
3582     for (z = 0 ; z < depth ; z++){
3583       for (y = 0 ; y < height ; y++){
3584 
3585         switch(mri_dst->type){
3586         case MRI_UCHAR: pdst  =           mri_dst->slices[s][y] ; break;
3587         case MRI_SHORT: psdst = (short *) mri_dst->slices[s][y] ; break;
3588         case MRI_INT:   pidst = (int *)   mri_dst->slices[s][y] ; break;
3589         case MRI_LONG:  pldst = (long *)  mri_dst->slices[s][y] ; break;
3590         case MRI_FLOAT: pfdst = (float *) mri_dst->slices[s][y] ; break;
3591         }
3592         switch(mri1->type){
3593         case MRI_UCHAR:
3594           p1 = mri1->slices[s][y] ;
3595           p2 = mri2->slices[s][y] ;
3596           break;
3597         case MRI_SHORT:
3598           ps1 = (short *) mri1->slices[s][y] ;
3599           ps2 = (short *) mri2->slices[s][y] ;
3600           break;
3601         case MRI_INT:
3602           pi1 = (int *) mri1->slices[s][y] ;
3603           pi2 = (int *) mri2->slices[s][y] ;
3604           break;
3605         case MRI_LONG:
3606           pl1 = (long *) mri1->slices[s][y] ;
3607           pl2 = (long *) mri2->slices[s][y] ;
3608           break;
3609         case MRI_FLOAT:
3610           pf1 = (float *) mri1->slices[s][y] ;
3611           pf2 = (float *) mri2->slices[s][y] ;
3612           break;
3613         }
3614 
3615         for (x = 0 ; x < width ; x++){
3616 
3617           switch(mri1->type){
3618           case MRI_UCHAR: v = (float)(*p1++)  - (float)(*p2++);  break;
3619           case MRI_SHORT: v = (float)(*ps1++) - (float)(*ps2++); break;
3620           case MRI_INT:   v = (float)(*pi1++) - (float)(*pi2++); break;
3621           case MRI_LONG:  v = (float)(*pl1++) - (float)(*pl2++); break;
3622           case MRI_FLOAT: v = (float)(*pf1++) - (float)(*pf2++); break;
3623           }
3624 
3625           switch(mri_dst->type){
3626           case MRI_UCHAR: (*pdst++)  = (BUFTYPE) v; break;
3627           case MRI_SHORT: (*psdst++) = (short)   v; break;
3628           case MRI_INT:   (*pidst++) = (int)     v; break;
3629           case MRI_LONG:  (*pldst++) = (long)    v; break;
3630           case MRI_FLOAT: (*pfdst++) = (float)   v; break;
3631           }
3632 
3633         }
3634       }
3635       s++;
3636     }
3637   }
3638 
3639   return(mri_dst) ;
3640 }
3641 #if 0
3642 /*------------------------------------------------------
3643   MRIsubtract(mri1,mri2,mridiff) - computes mri1-mri2.
3644   ------------------------------------------------------*/
3645 MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst)
3646 {
3647   int     nframes, width, height, depth, x, y, z, f ;
3648   float v1, v2, vdiff;
3649 
3650   width = mri1->width ;
3651   height = mri1->height ;
3652   depth = mri1->depth ;
3653   nframes = mri1->nframes ;
3654   if(nframes == 0) nframes = 1;
3655 
3656   if (!mri_dst){
3657     mri_dst = MRIallocSequence(width, height, depth, mri1->type,nframes) ;
3658     MRIcopyHeader(mri1, mri_dst) ;
3659   }
3660 
3661   for (z = 0 ; z < depth ; z++){
3662     for (y = 0 ; y < height ; y++){
3663       for (x = 0 ; x < width ; x++){
3664         for (f = 0 ; f < nframes ; f++){
3665           v1 = MRIgetVoxVal(mri1,x,y,z,f);
3666           v2 = MRIgetVoxVal(mri2,x,y,z,f);
3667           vdiff = v1-v2;
3668           MRIsetVoxVal(mri_dst,x,y,z,f,vdiff);
3669         }
3670       }
3671     }
3672   }
3673   return(mri_dst) ;
3674 }
3675 #endif
3676 /*-----------------------------------------------------
3677   ------------------------------------------------------*/
3678 MRI *
3679 MRIabsdiff(MRI *mri1, MRI *mri2, MRI *mri_dst)
3680 {
3681   int     width, height, depth, x, y, z ;
3682   BUFTYPE *p1, *p2, *pdst, v1, v2 ;
3683   float   f1, f2 ;
3684 
3685   width = mri1->width ;
3686   height = mri1->height ;
3687   depth = mri1->depth ;
3688 
3689   if (!mri_dst)
3690     {
3691       mri_dst = MRIalloc(width, height, depth, mri1->type) ;
3692       MRIcopyHeader(mri1, mri_dst) ;
3693     }
3694 
3695   if (mri1->type == MRI_UCHAR && mri2->type == MRI_UCHAR)
3696     {
3697       for (z = 0 ; z < depth ; z++)
3698         {
3699           for (y = 0 ; y < height ; y++)
3700             {
3701               p1 = mri1->slices[z][y] ;
3702               p2 = mri2->slices[z][y] ;
3703               pdst = mri_dst->slices[z][y] ;
3704               for (x = 0 ; x < width ; x++)
3705                 {
3706                   v1 = *p1++ ;
3707                   v2 = *p2++ ;
3708                   if (v1 > v2)
3709                     *pdst++ = v1 - v2 ;
3710                   else
3711                     *pdst++ = v2 - v1 ;
3712                 }
3713             }
3714         }
3715     }
3716   else
3717     {
3718       for (z = 0 ; z < depth ; z++)
3719         {
3720           for (y = 0 ; y < height ; y++)
3721             {
3722               for (x = 0 ; x < width ; x++)
3723                 {
3724                   f1 = MRIgetVoxVal(mri1, x, y, z, 0)  ;
3725                   f2 = MRIgetVoxVal(mri2, x, y, z, 0)  ;
3726                   MRIsetVoxVal(mri_dst, x, y, z, 0, fabs(f1 - f2)) ;
3727                 }
3728             }
3729         }
3730     }
3731   return(mri_dst) ;
3732 }
3733 /*-----------------------------------------------------
3734   ------------------------------------------------------*/
3735 MRI *MRIabs(MRI *mri_src, MRI *mri_dst)
3736 {
3737   int   width, height, depth, nframes, x, y, z,f ;
3738   float val;
3739 
3740   width = mri_src->width ;
3741   height = mri_src->height ;
3742   depth = mri_src->depth ;
3743   nframes = mri_src->nframes;
3744 
3745   if (!mri_dst){
3746     mri_dst = MRIallocSequence(width, height, depth, mri_src->type, nframes) ;
3747     MRIcopyHeader(mri_src, mri_dst) ;
3748   }
3749 
3750   for (z = 0 ; z < depth ; z++){
3751     for (y = 0 ; y < height ; y++){
3752       for (x = 0 ; x < width ; x++){
3753         for (f = 0 ; f < nframes ; f++){
3754           val = fabs(MRIgetVoxVal(mri_src,x,y,z,f));
3755           MRIsetVoxVal(mri_src,x,y,z,f,val);
3756         }
3757       }
3758     }
3759   }
3760   return(mri_dst) ;
3761 }
3762 /*-----------------------------------------------------*/
3763 MRI * MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst)
3764 {
3765   int     nframes, width, height, depth, x, y, z, f, s ;
3766   float v1, v2, v=0.0;
3767   BUFTYPE *p1=NULL, *p2=NULL, *pdst=NULL ;
3768   short *ps1=NULL, *ps2=NULL, *psdst=NULL ;
3769   int   *pi1=NULL, *pi2=NULL, *pidst=NULL ;
3770   long  *pl1=NULL, *pl2=NULL, *pldst=NULL ;
3771   float *pf1=NULL, *pf2=NULL, *pfdst=NULL ;
3772 
3773   width = mri1->width ;
3774   height = mri1->height ;
3775   depth = mri1->depth ;
3776   nframes = mri1->nframes ;
3777   if(nframes == 0) nframes = 1;
3778 
3779   if (!mri_dst){
3780     mri_dst = MRIallocSequence(width, height, depth, mri1->type,nframes) ;
3781     MRIcopyHeader(mri1, mri_dst) ;
3782   }
3783 
3784   if(mri1->type == MRI_UCHAR || (mri1->type != mri2->type)){
3785     /* Generic but slow */
3786     for (f = 0 ; f < nframes ; f++){
3787       for (z = 0 ; z < depth ; z++){
3788         for (y = 0 ; y < height ; y++){
3789           for (x = 0 ; x < width ; x++){
3790             v1 = MRIgetVoxVal(mri1,x,y,z,f);
3791             v2 = MRIgetVoxVal(mri2,x,y,z,f);
3792             v = v1+v2;
3793             if (mri_dst->type == MRI_UCHAR && v > 255)
3794               v = 255 ;
3795             MRIsetVoxVal(mri_dst,x,y,z,f,v);
3796           }
3797         }
3798       }
3799     }
3800     return(mri_dst) ;
3801   }
3802 
3803 
3804   s = 0;
3805   for (f = 0 ; f < nframes ; f++){
3806     for (z = 0 ; z < depth ; z++){
3807       for (y = 0 ; y < height ; y++){
3808 
3809         switch(mri_dst->type){
3810         case MRI_UCHAR: pdst  = mri_dst->slices[s][y] ; break;
3811         case MRI_SHORT: psdst = (short *) mri_dst->slices[s][y] ; break;
3812         case MRI_INT:   pidst = (int *)   mri_dst->slices[s][y] ; break;
3813         case MRI_LONG:  pldst = (long *)  mri_dst->slices[s][y] ; break;
3814         case MRI_FLOAT: pfdst = (float *) mri_dst->slices[s][y] ; break;
3815         }
3816 
3817         switch(mri1->type){
3818         case MRI_UCHAR:
3819           p1 = mri1->slices[s][y] ;
3820           p2 = mri2->slices[s][y] ;
3821           break;
3822         case MRI_SHORT:
3823           ps1 = (short *) mri1->slices[s][y] ;
3824           ps2 = (short *) mri2->slices[s][y] ;
3825           break;
3826         case MRI_INT:
3827           pi1 = (int *) mri1->slices[s][y] ;
3828           pi2 = (int *) mri2->slices[s][y] ;
3829           break;
3830         case MRI_LONG:
3831           pl1 = (long *) mri1->slices[s][y] ;
3832           pl2 = (long *) mri2->slices[s][y] ;
3833           break;
3834         case MRI_FLOAT:
3835           pf1 = (float *) mri1->slices[s][y] ;
3836           pf2 = (float *) mri2->slices[s][y] ;
3837           break;
3838         }
3839 
3840         for (x = 0 ; x < width ; x++){
3841 
3842           switch(mri_dst->type){
3843           case MRI_UCHAR:
3844             switch(mri1->type){
3845             case MRI_UCHAR: (*pdst++) = (BUFTYPE) ((*p1++)+(*p2++));   break;
3846             case MRI_SHORT: (*pdst++) = (BUFTYPE) ((*ps1++)+(*ps2++)); break;
3847             case MRI_INT:   (*pdst++) = (BUFTYPE) ((*pi1++)+(*pi2++)); break;
3848             case MRI_LONG:  (*pdst++) = (BUFTYPE) ((*pl1++)+(*pl2++)); break;
3849             case MRI_FLOAT: (*pdst++) = (BUFTYPE) ((*pf1++)+(*pf2++)); break;
3850             }
3851             break;
3852           case MRI_SHORT:
3853             switch(mri1->type){
3854             case MRI_UCHAR: (*psdst++) = ((short)(*p1++)+(*p2++));   break;
3855             case MRI_SHORT: (*psdst++) = (short) ((*ps1++)+(*ps2++)); break;
3856             case MRI_INT:   (*psdst++) = (short) ((*pi1++)+(*pi2++)); break;
3857             case MRI_LONG:  (*psdst++) = (short) ((*pl1++)+(*pl2++)); break;
3858             case MRI_FLOAT: (*psdst++) = (short) ((*pf1++)+(*pf2++)); break;
3859             }
3860             break;
3861           case MRI_INT:
3862             switch(mri1->type){
3863             case MRI_UCHAR: (*pidst++) = ((int)(*p1++)+(*p2++));   break;
3864             case MRI_SHORT: (*pidst++) = ((int)(*ps1++)+(*ps2++)); break;
3865             case MRI_INT:   (*pidst++) = (int) ((*pi1++)+(*pi2++)); break;
3866             case MRI_LONG:  (*pidst++) = (int) ((*pl1++)+(*pl2++)); break;
3867             case MRI_FLOAT: (*pidst++) = (int) ((*pf1++)+(*pf2++)); break;
3868             }
3869             break;
3870           case MRI_LONG:
3871             switch(mri1->type){
3872             case MRI_UCHAR: (*pldst++) = ((long)(*p1++)+(*p2++));   break;
3873             case MRI_SHORT: (*pldst++) = ((long)(*ps1++)+(*ps2++)); break;
3874             case MRI_INT:   (*pldst++) = ((long)(*pi1++)+(*pi2++)); break;
3875             case MRI_LONG:  (*pldst++) = (long) ((*pl1++)+(*pl2++)); break;
3876             case MRI_FLOAT: (*pldst++) = (long) ((*pf1++)+(*pf2++)); break;
3877             }
3878             break;
3879           case MRI_FLOAT:
3880             switch(mri1->type){
3881             case MRI_UCHAR: (*pfdst++) = ((float)(*p1++)+(*p2++)); break;
3882             case MRI_SHORT: (*pfdst++) = ((float)(*ps1++)+(*ps2++)); break;
3883             case MRI_INT:   (*pfdst++) = ((float)(*pi1++)+(*pi2++)); break;
3884             case MRI_LONG:  (*pfdst++) = ((float)(*pl1++)+(*pl2++)); break;
3885             case MRI_FLOAT: (*pfdst++) =         (*pf1++)+(*pf2++);  break;
3886             }
3887             break;
3888           }
3889         }
3890       }
3891       s++;
3892     }
3893   }
3894 
3895   return(mri_dst) ;
3896 }
3897 /*-----------------------------------------------------------
3898   MRIaverage() - computes average of source and destination.
3899   ------------------------------------------------------*/
3900 MRI *
3901 MRIaverage(MRI *mri_src, int dof, MRI *mri_dst)
3902 {
3903   int     width, height, depth, x, y, z, f ;
3904   Real    src, dst ;
3905 
3906   width = mri_src->width ;
3907   height = mri_src->height ;
3908   depth = mri_src->depth ;
3909 
3910   if (!mri_dst)
3911     {
3912       mri_dst = MRIalloc(width, height, depth, mri_src->type) ;
3913       MRIcopyHeader(mri_src, mri_dst) ;
3914     }
3915 
3916   if (!MRIcheckSize(mri_src, mri_dst,0,0,0))
3917     ErrorReturn(NULL,
3918                 (ERROR_BADPARM,"MRIaverage: incompatible volume dimensions"));
3919 #if 0
3920   if ((mri_src->type != MRI_UCHAR) || (mri_dst->type != MRI_UCHAR))
3921     ErrorReturn(NULL,
3922                 (ERROR_UNSUPPORTED,
3923                  "MRISaverage: unsupported voxel format %d",mri_src->type));
3924 #endif
3925 
3926   for (f = 0 ; f < mri_src->nframes ; f++)
3927     {
3928       for (z = 0 ; z < depth ; z++)
3929         {
3930           for (y = 0 ; y < height ; y++)
3931             {
3932               for (x = 0 ; x < width ; x++)
3933                 {
3934                   MRIsampleVolumeFrameType
3935                     (mri_src, x,  y, z, f, SAMPLE_NEAREST, &src) ;
3936                   MRIsampleVolumeFrameType
3937                     (mri_dst, x,  y, z, f, SAMPLE_NEAREST, &dst) ;
3938                   MRIsetVoxVal
3939                     (mri_dst, x, y, z, f, (dst*dof+src)/(Real)(dof+1))  ;
3940                 }
3941             }
3942         }
3943     }
3944   return(mri_dst) ;
3945 }
3946 /*-----------------------------------------------------
3947   Parameters:
3948 
3949   Returns value:
3950 
3951   Description
3952 
3953   ------------------------------------------------------*/
3954 MRI *
3955 MRImultiply(MRI *mri1, MRI *mri2, MRI *mri_dst)
3956 {
3957   int     width, height, depth, x, y, z ;
3958   float   f1, f2 ;
3959 
3960   width = mri1->width ;
3961   height = mri1->height ;
3962   depth = mri1->depth ;
3963 
3964   if (!mri_dst)
3965     {
3966       mri_dst = MRIalloc(width, height, depth, mri1->type) ;
3967       MRIcopyHeader(mri1, mri_dst) ;
3968     }
3969 
3970   for (z = 0 ; z < depth ; z++)
3971     {
3972       for (y = 0 ; y < height ; y++)
3973         {
3974           for (x = 0 ; x < width ; x++)
3975             {
3976               f1 = MRIgetVoxVal(mri1, x, y, z, 0) ;
3977               f2 = MRIgetVoxVal(mri2, x, y, z, 0) ;
3978               MRIsetVoxVal(mri_dst, x, y, z, 0, f1*f2) ;
3979             }
3980         }
3981     }
3982   return(mri_dst) ;
3983 }
3984 /*-----------------------------------------------------
3985   Parameters:
3986 
3987   Returns value:
3988 
3989   Description
3990 
3991   ------------------------------------------------------*/
3992 MRI *
3993 MRIscaleAndMultiply(MRI *mri1, float scale, MRI *mri2, MRI *mri_dst)
3994 {
3995   int     width, height, depth, x, y, z ;
3996   BUFTYPE *p1, *p2, *pdst ;
3997   float   out_val ;
3998 
3999   width = mri1->width ;
4000   height = mri1->height ;
4001   depth = mri1->depth ;
4002 
4003   if (!mri_dst)
4004     {
4005       mri_dst = MRIalloc(width, height, depth, mri1->type) ;
4006       MRIcopyHeader(mri1, mri_dst) ;
4007     }
4008 
4009   for (z = 0 ; z < depth ; z++)
4010     {
4011       for (y = 0 ; y < height ; y++)
4012         {
4013           p1 = mri1->slices[z][y] ;
4014           p2 = mri2->slices[z][y] ;
4015           pdst = mri_dst->slices[z][y] ;
4016           for (x = 0 ; x < width ; x++)
4017             {
4018               out_val = *p1++ * (*p2++/scale) ;
4019               if (out_val > 255)
4020                 out_val = 255 ;
4021               else if (out_val < 0)
4022                 out_val = 0 ;
4023               *pdst++ = (BUFTYPE)nint(out_val) ;
4024             }
4025         }
4026     }
4027   return(mri_dst) ;
4028 }
4029 /*-----------------------------------------------------
4030   Parameters:
4031 
4032   Returns value:
4033 
4034   Description
4035 
4036   ------------------------------------------------------*/
4037 MRI *
4038 MRIdivide(MRI *mri1, MRI *mri2, MRI *mri_dst)
4039 {
4040   int     width, height, depth, x, y, z ;
4041   BUFTYPE *p1, *p2, *pdst ;
4042 
4043   width = mri1->width ;
4044   height = mri1->height ;
4045   depth = mri1->depth ;
4046 
4047   if (!mri_dst)
4048     {
4049       mri_dst = MRIalloc(width, height, depth, mri1->type) ;
4050       MRIcopyHeader(mri1, mri_dst) ;
4051     }
4052 
4053   if (mri1->type != MRI_UCHAR || mri2->type != MRI_UCHAR)
4054     {
4055       Real val1, val2, dst ;
4056 
4057       for (z = 0 ; z < depth ; z++)
4058         {
4059           for (y = 0 ; y < height ; y++)
4060             {
4061               pdst = mri_dst->slices[z][y] ;
4062               for (x = 0 ; x < width ; x++)
4063                 {
4064                   if (x == Gx && y == Gy && z==Gz)
4065                     DiagBreak() ;
4066                   val1 = MRIgetVoxVal(mri1, x, y, z, 0) ;
4067                   val2 = MRIgetVoxVal(mri2, x, y, z, 0) ;
4068                   if  (FZERO(val2))
4069                     {
4070                       dst = FZERO(val1) ? 0 : 255 ;
4071                     }
4072                   else
4073                     dst = val1 / val2 ;
4074                   if (abs(dst) > 1000)
4075                     DiagBreak() ;
4076                   MRIsetVoxVal(mri_dst, x, y, z, 0, dst) ;
4077                 }
4078             }
4079         }
4080     }
4081   else   /* both UCHAR volumes */
4082     {
4083       for (z = 0 ; z < depth ; z++)
4084         {
4085           for (y = 0 ; y < height ; y++)
4086             {
4087               p1 = mri1->slices[z][y] ;
4088               p2 = mri2->slices[z][y] ;
4089               pdst = mri_dst->slices[z][y] ;
4090               for (x = 0 ; x < width ; x++)
4091                 {
4092                   if  (!*p2)
4093                     {
4094                       *pdst = FZERO(*p1) ? 0 : 255 ;
4095                       p2++ ;
4096                     }
4097                   else
4098                     *pdst++ = *p1++ / *p2++ ;
4099                 }
4100             }
4101         }
4102     }
4103   return(mri_dst) ;
4104 }
4105 /*-----------------------------------------------------
4106   MRIclone() - create a copy of an mri struct. Copies
4107   header info and allocs the pixel space (but does not
4108   copy pixel data).
4109   ------------------------------------------------------*/
4110 MRI *MRIclone(MRI *mri_src, MRI *mri_dst)
4111 {
4112   if (!mri_dst)
4113     mri_dst =
4114       MRIallocSequence(mri_src->width, mri_src->height,mri_src->depth,
4115                        mri_src->type, mri_src->nframes);
4116 
4117   MRIcopyHeader(mri_src, mri_dst) ;
4118   return(mri_dst) ;
4119 }
4120 /*---------------------------------------------------------------------
4121   MRIcloneSpace() - create a copy of an mri struct but allows user to
4122   set nframes (ie, all the spatial stuff is copied). Copies header
4123   info and allocs the pixel space (but does not copy pixel data).
4124   -------------------------------------------------------------------*/
4125 MRI *MRIcloneBySpace(MRI *mri_src, int nframes)
4126 {
4127   MRI *mri_dst;
4128   mri_dst = MRIallocSequence(mri_src->width, mri_src->height,mri_src->depth,
4129                              mri_src->type, nframes);
4130   MRIcopyHeader(mri_src, mri_dst) ;
4131   mri_dst->nframes = nframes;
4132   return(mri_dst) ;
4133 }
4134 /*-----------------------------------------------------
4135   Description
4136   Copy one MRI into another (including header info)
4137   ------------------------------------------------------*/
4138 MRI *
4139 MRIcloneRoi(MRI *mri_src, MRI *mri_dst)
4140 {
4141   int  w, h, d ;
4142 
4143   w = mri_src->width - mri_src->roi.x ;
4144   h = mri_src->height - mri_src->roi.y ;
4145   d = mri_src->depth - mri_src->roi.z ;
4146   mri_dst = MRIallocSequence(w, h, d, MRI_FLOAT, mri_src->nframes) ;
4147   MRIcopyHeader(mri_src, mri_dst) ;
4148   mri_dst->xstart = mri_src->xstart + mri_src->roi.x * mri_src->xsize ;
4149   mri_dst->ystart = mri_src->ystart + mri_src->roi.y * mri_src->ysize ;
4150   mri_dst->zstart = mri_src->zstart + mri_src->roi.z * mri_src->zsize ;
4151   mri_dst->xend = mri_src->xstart + w * mri_src->xsize ;
4152   mri_dst->yend = mri_src->ystart + h * mri_src->ysize ;
4153   mri_dst->zend = mri_src->zstart + d * mri_src->zsize ;
4154   return(mri_dst) ;
4155 }
4156 /*-----------------------------------------------------
4157   Parameters:
4158 
4159   Returns value:
4160 
4161   Description
4162   Copy one MRI into another (including header info and data)
4163   ------------------------------------------------------*/
4164 MRI *
4165 MRIcopy(MRI *mri_src, MRI *mri_dst)
4166 {
4167   int     width, height, depth, bytes, x, y, z, frame, val ;
4168   float   *fdst, *fsrc ;
4169   BUFTYPE *csrc, *cdst ;
4170   int     dest_ptype, *isrc;
4171   short   *ssrc, *sdst ;
4172 
4173   if (mri_src == mri_dst)
4174     return(mri_dst) ;
4175   width = mri_src->width ;
4176   height = mri_src->height ;
4177   depth = mri_src->depth ;
4178 
4179   if (!mri_dst)
4180     {
4181       if(mri_src->slices)
4182         mri_dst = MRIallocSequence(width, height, depth, mri_src->type,
4183                                    mri_src->nframes) ;
4184       else
4185         mri_dst = MRIallocHeader(width, height, depth, mri_src->type);
4186     }
4187   dest_ptype = mri_dst->ptype;
4188   MRIcopyHeader(mri_src, mri_dst) ;
4189   mri_dst->ptype = dest_ptype;
4190 
4191   if(!mri_src->slices)
4192     return(mri_dst);
4193 
4194   if (mri_src->type == mri_dst->type)
4195     {
4196       bytes = width ;
4197       switch (mri_src->type)
4198         {
4199         case MRI_UCHAR