Hello,
I have a similar (but hopefully simpler) problem.
I normalized a T1 image to an MNI template (the single subject one). I used
ONLY affine transformations (0 non linear iterations).
How do I calculate the transform W s.t. a point in the 3D image-matrix P
(i,j,k) from my original un-transformed T1 data is sent to Q(i,j,k), the
correspondent point in the SPM-affine-transformed 3D image-matrix?
Basically, I'm trying to extract W, s.t. Q = W*P using the information
stored in the _sn.mat file (Affine, VF.mat, VG.mat). Q and P are point in
the .img of the analyze file, i.e. the simple matrices, no axis
orientations, origins, etc...
Thanks for your help! This issue has been hunting me for almost a week!!
Best regards,
Anna
On Fri, 23 Apr 1999 12:12:45 +0100, John Ashburner <[log in to unmask]>
wrote:
>----------
>X-Sun-Data-Type: text
>X-Sun-Data-Description: text
>X-Sun-Data-Name: text
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 57
>
>
>> We have normalized our T1 image on T1 template from MNI. A sn3d.mat file
>> was created. It contains Affine, Dims, MF, MG, Transform and mgc matrix.
>> We need to convert, coordinates (x y z) of a voxel belonging to a T2
image
>> realigned on our T1 image in the Talairach space (of T1 template). We
want
>> not to modify the intensity of this voxel, just getting the new (x' y'
z')
>> coordinates. We think that it is possible to calculate these coordinates
>> with the parameters of the sn3d. mat file; can you give us the operations
>> to performed and how using these parameters.
>
>I'm afraid that what you want to do may sound very simple, but it is
actually
>a little bit complicated.
>If you are not really interested - then look away now.
>
>
>I have appended a Matlab 5 function for writing out deformation fields
>(spm_write_defs.m). It also requires a few other new matlab functions
>(also attached) as well as some routines from the spm96 distribution.
>You may need to take a look at spm_write_sn.m to see how the deformation
>field could be written for an image that has been co-registered with
>the original image. Note that the values in the deformation fields assume
>that the first voxel in an image is at (0,0,0) rather than (1,1,1) (which
>all the stuff in the official spm releases assume).
>
>This function contains the information about how a voxel in a normalized
>image should map to a voxel in an un-normalized image, but what I think
>you are after is a function that maps from a voxel in an un-normalized
>image to one in a normalized one. This requires that the spatial
>transformations are inverted, which is not a trivial problem. I have
>attached some C routines that can be used to approximately invert a
>deformation field. Compilation is something like:
>
> cc -o invert_def invert_def.c voxels.c penalties.c warpio.c -lm -O2
>
>The functions need the deformations to be stored as contiguous floating
point
>volumes. The arguments are:
>
>invert_def xdim_g,ydim_g,zdim_g/xdim_f,ydim_f,zdim_f y0 y1 y2 iy0 iy1 iy2
>
>where:
> xdim_g, ydim_g and zdim_g are the dimensions of the original
> deformation field.
> xdim_f, ydim_f and zdim_f are the dimensions that the inverted
> deformation field should have.
> y0, y1 and y2 are the names of the files containing the original
> deformation field
> iy0, iy1 and iy2 are the names that the inverted deformation field
> should have.
>
>After running the program, if you sample voxel (x,y,z) of the iy0 file,
then
>it should give you the corresponding x coordinate in the normalized image.
>Sampling the the same point in file iy1 will give you the y coordinate, and
>sampling iy2 will give you the z coordinate.
>
>
>Good luck,
>-John
>----------
>X-Sun-Data-Type: h-file
>X-Sun-Data-Description: h-file
>X-Sun-Data-Name: tets.h
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 9
>
>typedef struct
>{
> REAL ix[4][4];
> REAL dix[4][3];
> REAL volx;
> int o1, o2, o3;
> unsigned long mask;
>} Tettype;
>
>----------
>X-Sun-Data-Type: c-file
>X-Sun-Data-Description: c-file
>X-Sun-Data-Name: invert_def.c
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 169
>
>#ifndef lint
>static char sccsid[]="%W% John Ashburner %E%";
>#endif
>
>#define REAL float
>
>#include<math.h>
>#include<stdio.h>
>#include<malloc.h>
>#define MAXV 2048
>
>extern void invertX(), getM(), invertM(), scan_tetrahedron();
>extern int read_warp(), write_warp();
>
>static REAL x[2][5][4][3] = {
>{{{ 0,0,0},{ 1,0,1},{ 1,0,0},{ 1,1,0}},
> {{ 0,0,0},{ 1,0,1},{ 0,1,1},{ 0,0,1}},
> {{ 0,0,0},{ 0,1,0},{ 0,1,1},{ 1,1,0}},
> {{ 0,0,0},{ 1,0,1},{ 1,1,0},{ 0,1,1}},
> {{ 1,1,1},{ 1,1,0},{ 0,1,1},{ 1,0,1}},},
>
>{{{ 1,0,0},{ 0,0,1},{ 0,0,0},{ 0,1,0}},
> {{ 1,0,0},{ 0,0,1},{ 1,1,1},{ 1,0,1}},
> {{ 1,0,0},{ 1,1,0},{ 1,1,1},{ 0,1,0}},
> {{ 1,0,0},{ 0,0,1},{ 0,1,0},{ 1,1,1}},
> {{ 0,1,1},{ 0,1,0},{ 1,1,1},{ 0,0,1}},},
>};
>
>static REAL ix[2][5][4][4], dix[2][5][4][3];
>static int off[2][4][5];
>
>static void setup_consts(int dim[3])
>{
> int i, k;
> for(k=0; k<2; k++)
> for(i=0; i<5; i++)
> {
> REAL dtx;
> int j;
> invertX(x[k][i], ix[k][i], dix[k][i], &dtx);
> for(j=0; j<4; j++)
> off[k][j][i] = x[k][i][j][0]+dim[0]*(x[k][i]
[j][1]+dim[1]*x[k][i][j][2]);
> }
>}
>
>static void invert_it(int x0, int x1, int x2, float *y0, float *y1, float
*y2, int dim_f[3], float *iy0, float *iy1, float *iy2)
>{
> int i, k;
> REAL d, ix30;
> k = (x0 + x1 + x2 +1)%2;
> for(i=0; i<5; i++)
> {
> REAL Y[4][3], M[4][3], IM[4][3];
> int j, vox[MAXV][3], nvox;
>
> Y[0][0] = y0[off[k][0][i]]; Y[0][1] = y1[off[k][0][i]]; Y[0]
[2] = y2[off[k][0][i]];
> Y[1][0] = y0[off[k][1][i]]; Y[1][1] = y1[off[k][1][i]]; Y[1]
[2] = y2[off[k][1][i]];
> Y[2][0] = y0[off[k][2][i]]; Y[2][1] = y1[off[k][2][i]]; Y[2]
[2] = y2[off[k][2][i]];
> Y[3][0] = y0[off[k][3][i]]; Y[3][1] = y1[off[k][3][i]]; Y[3]
[2] = y2[off[k][3][i]];
>
> getM(Y, ix[k][i], dix[k][i], M, x0, x1, x2, &ix30);
> invertM(M, IM, &d);
>
> scan_tetrahedron(Y, &nvox, vox, MAXV);
> for(j=0; j<nvox; j++)
> {
> if ((vox[j][0]>=0) && (vox[j][0]<dim_f[0]) &&
> (vox[j][1]>=0) && (vox[j][1]<dim_f[1]) &&
> (vox[j][2]>=0) && (vox[j][2]<dim_f[2]))
> {
> int o = vox[j][0]+dim_f[0]*(vox[j][1]+dim_f
[1]*vox[j][2]);
> iy0[o] = IM[0][0]*vox[j][0] + IM[1][0]*vox
[j][1] + IM[2][0]*vox[j][2] + IM[3][0];
> iy1[o] = IM[0][1]*vox[j][0] + IM[1][1]*vox
[j][1] + IM[2][1]*vox[j][2] + IM[3][1];
> iy2[o] = IM[0][2]*vox[j][0] + IM[1][2]*vox
[j][1] + IM[2][2]*vox[j][2] + IM[3][2];
> }
> }
> }
>}
>
>
>
>static void invert_field(int dim_g[3], float y0[], float y1[], float y2
[],
> int dim_f[3], float iy0[], float iy1[], float iy2
[])
>{
> int x2, x1, x0;
> setup_consts(dim_g);
> for(x2=0; x2<dim_g[2]-1; x2++)
> {
> for(x1=0; x1<dim_g[1]-1; x1++)
> for(x0=0; x0<dim_g[0]-1; x0++)
> {
> int o = x0 + dim_g[0]*(x1 + x2*dim_g[1]);
> invert_it(x0, x1, x2, y0+o, y1+o, y2+o,
dim_f, iy0, iy1, iy2);
> }
> (void)printf(".");(void)fflush(stdout);
> }
> (void)printf("\n");(void)fflush(stdout);
>}
>
>static void setnan(float *dat, int n)
>{
> int j;
> float NaN = 0.0/0.0;
> for (j=0; j<n; j++)
> dat[j] = NaN;
>}
>
>#define cleanup {if (y0) {(void)free((void *)y0);y0=0;} if (y1) {(void)free
((void *)y1);y1=0;} if (y2) {(void)free((void *)y2);y2=0;} if (iy0) {(void)
free((void *)iy0);iy0=0;} if (iy1) {(void)free((void *)iy1);iy1=0;} if
(iy2) {(void)free((void *)iy2);iy2=0;}};
>
>main(int argc, char *argv[])
>{
> float *y0=0, *y1=0, *y2=0;
> float *iy0=0, *iy1=0, *iy2=0;
> int dim_g[3], dim_f[3];
>
> if (argc != 8)
> {
> (void)fprintf(stderr, "usage: %s %s %s %s %s %s %s %s\n",
> argv[0],
> "xdim_g,ydim_g,zdim_g/xdim_f,ydim_f,zdim_f",
> "y0", "y1", "y2", "iy0", "iy1", "iy2");
> return(1);
> }
>
> if (sscanf(argv[1], "%d,%d,%d/%d,%d,%d", &dim_g[0],&dim_g[1],&dim_g
[2], &dim_f[0],&dim_f[1],&dim_f[2]) != 6)
> {
> (void)fprintf(stderr,"Can't read dimensions\n");
> return(1);
> }
>
> y0 = (float *)calloc(dim_g[0]*dim_g[1]*dim_g[2],sizeof(float));
> y1 = (float *)calloc(dim_g[0]*dim_g[1]*dim_g[2],sizeof(float));
> y2 = (float *)calloc(dim_g[0]*dim_g[1]*dim_g[2],sizeof(float));
> iy0 = (float *)calloc(dim_f[0]*dim_f[1]*dim_f[2],sizeof(float));
> iy1 = (float *)calloc(dim_f[0]*dim_f[1]*dim_f[2],sizeof(float));
> iy2 = (float *)calloc(dim_f[0]*dim_f[1]*dim_f[2],sizeof(float));
>
> if (!y0 || !y1 || !y2 || !iy0 || !iy1 || !iy2)
> {
> (void)fprintf(stderr,"Out of memory\n");
> cleanup;
> return(1);
> }
>
> if (read_warp(argv[2], dim_g, y0) || read_warp(argv[3], dim_g, y1)
|| read_warp(argv[4], dim_g, y2))
> {
> (void)fprintf(stderr,"Error reading warps\n");
> cleanup;
> return(1);
> }
>
> setnan(iy0, dim_f[0]*dim_f[1]*dim_f[2]);
> setnan(iy1, dim_f[0]*dim_f[1]*dim_f[2]);
> setnan(iy2, dim_f[0]*dim_f[1]*dim_f[2]);
>
> invert_field(dim_g, y0, y1, y2, dim_f, iy0, iy1, iy2);
>
> if (write_warp(argv[5], dim_f, iy0) ||
> write_warp(argv[6], dim_f, iy1) ||
> write_warp(argv[7], dim_f, iy2))
> {
> (void)fprintf(stderr,"Error writing warps\n");
> cleanup;
> return(1);
> }
>
> cleanup;
> return(0);
>}
>----------
>X-Sun-Data-Type: c-file
>X-Sun-Data-Description: c-file
>X-Sun-Data-Name: voxels.c
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 168
>
>#ifndef lint
>static char sccsid[]="%W% John Ashburner %E%";
>#endif
>
>#define EPS 0.0
>#define REAL float
>#include<math.h>
>#include<stdio.h>
>#include<stdlib.h>
>
>static void scan_line(REAL lin[2], int y, int z, int *n, int vox[][3], int
maxv)
>{
> REAL p[2], t;
> int x, xe;
>
> /* sort p into ascending order of x */
> p[0] = lin[0]; p[1] = lin[1];
> if (p[1]<p[0]) {t = p[1]; p[1] = p[0]; p[0] = t;}
>
> /* find voxels where x is integer */
> for(x=ceil(p[0]), xe=floor(p[1]); x<=xe; x++)
> {
> if ((*n)>=maxv-1)
> {
> (void)fprintf(stderr,"Too many voxels inside a
tetrahedran\n");
> (void)exit(1);
> }
> vox[*n][0] = x;
> vox[*n][1] = y;
> vox[*n][2] = z;
> (*n)++;
> }
>}
>
>static void scan_triangle(REAL tri[][2], int z, int *n, int vox[][3], int
maxv)
>{
> REAL *p[3], *t, lin[2];
> REAL x1, x2, y1, y2;
> int y, ye, i;
>
> /* sort p into ascending order of y */
> p[0] = tri[0]; p[1] = tri[1]; p[2] = tri[2];
> if (p[1][1]<p[0][1]) {t = p[1]; p[1] = p[0]; p[0] = t;}
> if (p[2][1]<p[1][1]) {t = p[2]; p[2] = p[1]; p[1] = t;}
> if (p[1][1]<p[0][1]) {t = p[1]; p[1] = p[0]; p[0] = t;}
>
> /* find lower lines cutting triangle where y is integer */
> for(y=ceil(p[0][1]), ye=floor(p[1][1]); y<=ye; y++)
> {
> x1 = p[0][0]; y1 = p[0][1];
> for(i=0; i<2; i++)
> {
> x2 = p[i+1][0]; y2 = p[i+1][1];
> if (y2-y1<=EPS)
> lin[i] = (x1+x2)/2.0;
> else
> lin[i] = (x1*(y2-y)+x2*(y-y1))/(y2-y1);
> }
> scan_line(lin,y,z, n,vox,maxv);
> }
>
> /* find upper lines cutting triangle where y is integer */
> for(y=ceil(p[1][1]), ye=floor(p[2][1]); y<=ye; y++)
> {
> x2 = p[2][0]; y2 = p[2][1];
> for(i=0; i<2; i++)
> {
> x1 = p[i][0]; y1 = p[i][1];
> if (y2-y1<=EPS)
> lin[i] = (x1+x2)/2.0;
> else
> lin[i] = (x1*(y2-y)+x2*(y-y1))/(y2-y1);
> }
> scan_line(lin,y,z, n,vox,maxv);
> }
>}
>
>
>void scan_tetrahedron(REAL Y[4][3], int *n, int vox[][3], int maxv)
>{
> REAL *p[4], *t;
> REAL tri[4][2];
> REAL x1, x2, y1, y2, z1, z2;
> int z, ze, i;
>
> *n = 0;
>
> /* sort p into ascending order of z */
> p[0] = Y[0]; p[1] = Y[1]; p[2] = Y[2]; p[3] = Y[3];
> if (p[1][2]<p[0][2]) {t = p[1]; p[1] = p[0]; p[0] = t;}
> if (p[2][2]<p[1][2]) {t = p[2]; p[2] = p[1]; p[1] = t;}
> if (p[3][2]<p[2][2]) {t = p[3]; p[3] = p[2]; p[2] = t;}
> if (p[1][2]<p[0][2]) {t = p[1]; p[1] = p[0]; p[0] = t;}
> if (p[2][2]<p[1][2]) {t = p[2]; p[2] = p[1]; p[1] = t;}
> if (p[1][2]<p[0][2]) {t = p[1]; p[1] = p[0]; p[0] = t;}
>
> /* find lower triangles that intersect tetrahedron where z is
integer */
> for(z=ceil(p[0][2]), ze=floor(p[1][2]); z<=ze; z++)
> {
> x1 = p[0][0]; y1 = p[0][1]; z1 = p[0][2];
> for(i=0; i<3; i++)
> {
> x2 = p[i+1][0]; y2 = p[i+1][1]; z2 = p[i+1][2];
> if (z2-z1<=EPS)
> {
> tri[i][0] = (x1+x2)/2.0;
> tri[i][1] = (y1+y2)/2.0;
> }
> else
> {
> REAL t2 = z2-z, t1 = z-z1, t = z2-z1;
> tri[i][0] = (x1*t2+x2*t1)/t;
> tri[i][1] = (y1*t2+y2*t1)/t;
> }
> }
> scan_triangle(tri,z, n,vox,maxv);
> }
>
> /* find quadrilaterals that intersect tetrahedron where z is
integer */
> /* each quadrilateral divided into two triangles */
> for(z=ceil(p[1][2]), ze=floor(p[2][2]); z<=ze; z++)
> {
> static int ii[] = {0,1,1,0}, jj[] = {3,3,2,2};
>
> for(i=0; i<4; i++)
> {
> x1 = p[ii[i]][0]; y1 = p[ii[i]][1]; z1 = p[ii[i]]
[2];
> x2 = p[jj[i]][0]; y2 = p[jj[i]][1]; z2 = p[jj[i]]
[2];
> if (z2-z1<=EPS)
> {
> tri[i][0] = (x1+x2)/2.0;
> tri[i][1] = (y1+y2)/2.0;
> }
> else
> {
> REAL t2 = z2-z, t1 = z-z1, t = z2-z1;
> tri[i][0] = (x1*t2+x2*t1)/t;
> tri[i][1] = (y1*t2+y2*t1)/t;
> }
> }
> scan_triangle(tri,z, n,vox,maxv);
> tri[1][0] = tri[3][0];
> tri[1][1] = tri[3][1];
> scan_triangle(tri,z, n,vox,maxv);
> }
>
> /* find upper triangles that intersect tetrahedron where z is
integer */
> for(z=ceil(p[2][2]), ze=floor(p[3][2]); z<=ze; z++)
> {
> x2 = p[3][0]; y2 = p[3][1]; z2 = p[3][2];
> for(i=0; i<3; i++)
> {
> x1 = p[i][0]; y1 = p[i][1]; z1 = p[i][2];
> if (z2-z1<=EPS)
> {
> tri[i][0] = (x1+x2)/2.0;
> tri[i][1] = (y1+y2)/2.0;
> }
> else
> {
> REAL t2 = z2-z, t1 = z-z1, t = z2-z1;
> tri[i][0] = (x1*t2+x2*t1)/t;
> tri[i][1] = (y1*t2+y2*t1)/t;
> }
> }
> scan_triangle(tri,z, n,vox,maxv);
> }
>}
>----------
>X-Sun-Data-Type: c-file
>X-Sun-Data-Description: c-file
>X-Sun-Data-Name: penalties.c
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 493
>
>#ifndef lint
>static char sccsid[] = "%W% John Ashburner %E%";
>#endif
>
>#include<math.h>
>#define REAL float
>
>void invertX(REAL X[4][3], REAL IX[4][4], REAL dIX[4][3], REAL *dt)
>{
> REAL id;
> *dt = X[0][0]*(X[3][1]*(X[1][2]-X[2][2])+X[1][1]*(X[2][2]-X[3][2])+X
[2][1]*(X[3][2]-X[1][2]))+
> X[1][0]*(X[3][2]*(X[0][1]-X[2][1])+X[0][2]*(X[2][1]-X[3][1])+X
[2][2]*(X[3][1]-X[0][1]))+
> X[2][0]*(X[0][1]*(X[1][2]-X[3][2])+X[3][1]*(X[0][2]-X[1][2])+X
[1][1]*(X[3][2]-X[0][2]))+
> X[3][0]*(X[1][2]*(X[2][1]-X[0][1])+X[0][2]*(X[1][1]-X[2][1])+X
[2][2]*(X[0][1]-X[1][1]));
> id = 1/(*dt);
> IX[0][0] = id*(X[1][1]*(X[2][2]-X[3][2])+X[2][1]*(X[3][2]-X[1][2])+X
[3][1]*(X[1][2]-X[2][2]));
> IX[0][1] = id*(X[0][1]*(X[3][2]-X[2][2])+X[2][1]*(X[0][2]-X[3][2])+X
[3][1]*(X[2][2]-X[0][2]));
> IX[0][2] = id*(X[0][1]*(X[1][2]-X[3][2])+X[1][1]*(X[3][2]-X[0][2])+X
[3][1]*(X[0][2]-X[1][2]));
> IX[0][3] = id*(X[0][1]*(X[2][2]-X[1][2])+X[1][1]*(X[0][2]-X[2][2])+X
[2][1]*(X[1][2]-X[0][2]));
> IX[1][0] = id*(X[1][0]*(X[3][2]-X[2][2])+X[2][0]*(X[1][2]-X[3][2])+X
[3][0]*(X[2][2]-X[1][2]));
> IX[1][1] = id*(X[0][0]*(X[2][2]-X[3][2])+X[2][0]*(X[3][2]-X[0][2])+X
[3][0]*(X[0][2]-X[2][2]));
> IX[1][2] = id*(X[0][0]*(X[3][2]-X[1][2])+X[1][0]*(X[0][2]-X[3][2])+X
[3][0]*(X[1][2]-X[0][2]));
> IX[1][3] = id*(X[0][0]*(X[1][2]-X[2][2])+X[1][0]*(X[2][2]-X[0][2])+X
[2][0]*(X[0][2]-X[1][2]));
> IX[2][0] = id*(X[1][0]*(X[2][1]-X[3][1])+X[2][0]*(X[3][1]-X[1][1])+X
[3][0]*(X[1][1]-X[2][1]));
> IX[2][1] = id*(X[0][0]*(X[3][1]-X[2][1])+X[2][0]*(X[0][1]-X[3][1])+X
[3][0]*(X[2][1]-X[0][1]));
> IX[2][2] = id*(X[0][0]*(X[1][1]-X[3][1])+X[1][0]*(X[3][1]-X[0][1])+X
[3][0]*(X[0][1]-X[1][1]));
> IX[2][3] = id*(X[0][0]*(X[2][1]-X[1][1])+X[1][0]*(X[0][1]-X[2][1])+X
[2][0]*(X[1][1]-X[0][1]));
> IX[3][0] = id*(X[1][0]*(X[2][2]*X[3][1]-X[3][2]*X[2][1])+
> X[2][0]*(X[3][2]*X[1][1]-X[1][2]*X[3][1])+
> X[3][0]*(X[1][2]*X[2][1]-X[2][2]*X[1][1]));
> IX[3][1] = id*(X[0][0]*(X[3][2]*X[2][1]-X[2][2]*X[3][1])+
> X[2][0]*(X[0][2]*X[3][1]-X[3][2]*X[0][1])+
> X[3][0]*(X[2][2]*X[0][1]-X[0][2]*X[2][1]));
> IX[3][2] = id*(X[0][0]*(X[1][2]*X[3][1]-X[3][2]*X[1][1])+
> X[1][0]*(X[3][2]*X[0][1]-X[0][2]*X[3][1])+
> X[3][0]*(X[0][2]*X[1][1]-X[1][2]*X[0][1]));
> IX[3][3] = id*(X[0][0]*(X[2][2]*X[1][1]-X[1][2]*X[2][1])+
> X[1][0]*(X[0][2]*X[2][1]-X[2][2]*X[0][1])+
> X[2][0]*(X[1][2]*X[0][1]-X[0][2]*X[1][1]));
>
> dIX[0][0] = -( X[1][1]*X[2][2]-X[1][1]*X[3][2]-X[1][2]*X[2][1]+X[1]
[2]*X[3][1]+X[2][1]*X[3][2]-X[3][1]*X[2][2])*id;
> dIX[0][1] = -(-X[1][0]*X[2][2]+X[1][0]*X[3][2]+X[1][2]*X[2][0]-X[1]
[2]*X[3][0]-X[2][0]*X[3][2]+X[3][0]*X[2][2])*id;
> dIX[0][2] = -( X[1][0]*X[2][1]-X[1][0]*X[3][1]-X[1][1]*X[2][0]+X[1]
[1]*X[3][0]+X[2][0]*X[3][1]-X[3][0]*X[2][1])*id;
>
> dIX[1][0] = ( X[0][1]*X[2][2]-X[0][1]*X[3][2]-X[0][2]*X[2][1]+X[0]
[2]*X[3][1]+X[2][1]*X[3][2]-X[3][1]*X[2][2])*id;
> dIX[1][1] = ( X[3][0]*X[2][2]-X[2][0]*X[3][2]-X[0][0]*X[2][2]+X[0]
[0]*X[3][2]-X[0][2]*X[3][0]+X[0][2]*X[2][0])*id;
> dIX[1][2] = ( X[0][0]*X[2][1]-X[0][0]*X[3][1]-X[0][1]*X[2][0]+X[0]
[1]*X[3][0]+X[2][0]*X[3][1]-X[3][0]*X[2][1])*id;
>
> dIX[2][0] = -( X[0][1]*X[1][2]-X[0][1]*X[3][2]-X[0][2]*X[1][1]+X[0]
[2]*X[3][1]+X[1][1]*X[3][2]-X[1][2]*X[3][1])*id;
> dIX[2][1] = -(-X[0][0]*X[1][2]+X[0][0]*X[3][2]+X[0][2]*X[1][0]-X[0]
[2]*X[3][0]-X[1][0]*X[3][2]+X[1][2]*X[3][0])*id;
> dIX[2][2] = -( X[0][0]*X[1][1]-X[0][0]*X[3][1]-X[0][1]*X[1][0]+X[0]
[1]*X[3][0]+X[1][0]*X[3][1]-X[1][1]*X[3][0])*id;
>
> dIX[3][0] = ( X[0][1]*X[1][2]-X[0][1]*X[2][2]-X[0][2]*X[1][1]+X[0]
[2]*X[2][1]+X[1][1]*X[2][2]-X[1][2]*X[2][1])*id;
> dIX[3][1] = ( X[1][2]*X[2][0]-X[1][0]*X[2][2]+X[0][0]*X[2][2]-X[0]
[0]*X[1][2]-X[0][2]*X[2][0]+X[0][2]*X[1][0])*id;
> dIX[3][2] = ( X[0][0]*X[1][1]-X[0][0]*X[2][1]-X[0][1]*X[1][0]+X[0]
[1]*X[2][0]+X[1][0]*X[2][1]-X[1][1]*X[2][0])*id;
>}
>
>void getM(REAL Y[4][3], REAL IX[4][4], REAL dIX[4][3], REAL M[4][3], int
i, int j, int k, REAL *pix30)
>{
> REAL ix30, ix31, ix32, ix33;
>
> ix30 = IX[3][0] + i*dIX[0][0] + j*dIX[0][1] + k*dIX[0][2];
> ix31 = IX[3][1] + i*dIX[1][0] + j*dIX[1][1] + k*dIX[1][2];
> ix32 = IX[3][2] + i*dIX[2][0] + j*dIX[2][1] + k*dIX[2][2];
> ix33 = IX[3][3] + i*dIX[3][0] + j*dIX[3][1] + k*dIX[3][2];
>
> M[0][0] = Y[0][0]*IX[0][0] + Y[1][0]*IX[0][1] + Y[2][0]*IX[0][2] + Y
[3][0]*IX[0][3];
> M[0][1] = Y[0][1]*IX[0][0] + Y[1][1]*IX[0][1] + Y[2][1]*IX[0][2] + Y
[3][1]*IX[0][3];
> M[0][2] = Y[0][2]*IX[0][0] + Y[1][2]*IX[0][1] + Y[2][2]*IX[0][2] + Y
[3][2]*IX[0][3];
>
> M[1][0] = Y[0][0]*IX[1][0] + Y[1][0]*IX[1][1] + Y[2][0]*IX[1][2] + Y
[3][0]*IX[1][3];
> M[1][1] = Y[0][1]*IX[1][0] + Y[1][1]*IX[1][1] + Y[2][1]*IX[1][2] + Y
[3][1]*IX[1][3];
> M[1][2] = Y[0][2]*IX[1][0] + Y[1][2]*IX[1][1] + Y[2][2]*IX[1][2] + Y
[3][2]*IX[1][3];
>
> M[2][0] = Y[0][0]*IX[2][0] + Y[1][0]*IX[2][1] + Y[2][0]*IX[2][2] + Y
[3][0]*IX[2][3];
> M[2][1] = Y[0][1]*IX[2][0] + Y[1][1]*IX[2][1] + Y[2][1]*IX[2][2] + Y
[3][1]*IX[2][3];
> M[2][2] = Y[0][2]*IX[2][0] + Y[1][2]*IX[2][1] + Y[2][2]*IX[2][2] + Y
[3][2]*IX[2][3];
>
> M[3][0] = Y[0][0]*ix30 + Y[1][0]*ix31 + Y[2][0]*ix32 + Y
[3][0]*ix33;
> M[3][1] = Y[0][1]*ix30 + Y[1][1]*ix31 + Y[2][1]*ix32 + Y
[3][1]*ix33;
> M[3][2] = Y[0][2]*ix30 + Y[1][2]*ix31 + Y[2][2]*ix32 + Y
[3][2]*ix33;
>
> *pix30 = ix30;
>}
>
>static void getMpartial(REAL Y[4][3], REAL IX[4][4], REAL M[4][3])
>{
> M[0][0] = Y[0][0]*IX[0][0] + Y[1][0]*IX[0][1] + Y[2][0]*IX[0][2] + Y
[3][0]*IX[0][3];
> M[0][1] = Y[0][1]*IX[0][0] + Y[1][1]*IX[0][1] + Y[2][1]*IX[0][2] + Y
[3][1]*IX[0][3];
> M[0][2] = Y[0][2]*IX[0][0] + Y[1][2]*IX[0][1] + Y[2][2]*IX[0][2] + Y
[3][2]*IX[0][3];
>
> M[1][0] = Y[0][0]*IX[1][0] + Y[1][0]*IX[1][1] + Y[2][0]*IX[1][2] + Y
[3][0]*IX[1][3];
> M[1][1] = Y[0][1]*IX[1][0] + Y[1][1]*IX[1][1] + Y[2][1]*IX[1][2] + Y
[3][1]*IX[1][3];
> M[1][2] = Y[0][2]*IX[1][0] + Y[1][2]*IX[1][1] + Y[2][2]*IX[1][2] + Y
[3][2]*IX[1][3];
>
> M[2][0] = Y[0][0]*IX[2][0] + Y[1][0]*IX[2][1] + Y[2][0]*IX[2][2] + Y
[3][0]*IX[2][3];
> M[2][1] = Y[0][1]*IX[2][0] + Y[1][1]*IX[2][1] + Y[2][1]*IX[2][2] + Y
[3][1]*IX[2][3];
> M[2][2] = Y[0][2]*IX[2][0] + Y[1][2]*IX[2][1] + Y[2][2]*IX[2][2] + Y
[3][2]*IX[2][3];
>}
>
>void invertM(REAL M[4][3], REAL IM[4][3], REAL *d)
>{
> REAL id;
> *d = M[0][0]*(M[1][1]*M[2][2]-M[1][2]*M[2][1])+
> M[0][1]*(M[1][2]*M[2][0]-M[1][0]*M[2][2])+
> M[0][2]*(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
>
> id = 1.0/(*d);
> IM[0][0] = (M[1][1]*M[2][2]-M[1][2]*M[2][1])*id;
> IM[0][1] = (M[0][2]*M[2][1]-M[0][1]*M[2][2])*id;
> IM[0][2] = (M[0][1]*M[1][2]-M[0][2]*M[1][1])*id;
>
> IM[1][0] = (M[1][2]*M[2][0]-M[1][0]*M[2][2])*id;
> IM[1][1] = (M[0][0]*M[2][2]-M[0][2]*M[2][0])*id;
> IM[1][2] = (M[0][2]*M[1][0]-M[0][0]*M[1][2])*id;
>
> IM[2][0] = (M[1][0]*M[2][1]-M[1][1]*M[2][0])*id;
> IM[2][1] = (M[0][1]*M[2][0]-M[0][0]*M[2][1])*id;
> IM[2][2] = (M[0][0]*M[1][1]-M[0][1]*M[1][0])*id;
>
> IM[3][0] = (M[1][0]*(M[3][1]*M[2][2]-M[2][1]*M[3][2])+
> M[1][1]*(M[2][0]*M[3][2]-M[3][0]*M[2][2])+
> M[1][2]*(M[3][0]*M[2][1]-M[2][0]*M[3][1]))*id;
> IM[3][1] = (M[0][0]*(M[2][1]*M[3][2]-M[3][1]*M[2][2])+
> M[0][1]*(M[3][0]*M[2][2]-M[2][0]*M[3][2])+
> M[0][2]*(M[2][0]*M[3][1]-M[3][0]*M[2][1]))*id;
> IM[3][2] = (M[0][0]*(M[3][1]*M[1][2]-M[1][1]*M[3][2])+
> M[0][1]*(M[1][0]*M[3][2]-M[3][0]*M[1][2])+
> M[0][2]*(M[3][0]*M[1][1]-M[1][0]*M[3][1]))*id;
>}
>
>static void invertJ(REAL J[3][3], REAL IJ[3][3], REAL *dt)
>{
> REAL id;
> *dt = J[0][0]*(J[1][1]*J[2][2]-J[1][2]*J[2][1])+
> J[0][1]*(J[1][2]*J[2][0]-J[1][0]*J[2][2])+
> J[0][2]*(J[1][0]*J[2][1]-J[1][1]*J[2][0]);
> id = 1.0/(*dt);
>
> IJ[0][0] = (J[1][1]*J[2][2]-J[1][2]*J[2][1])*id;
> IJ[0][1] = (J[0][2]*J[2][1]-J[0][1]*J[2][2])*id;
> IJ[0][2] = (J[0][1]*J[1][2]-J[0][2]*J[1][1])*id;
>
> IJ[1][0] = (J[1][2]*J[2][0]-J[1][0]*J[2][2])*id;
> IJ[1][1] = (J[0][0]*J[2][2]-J[0][2]*J[2][0])*id;
> IJ[1][2] = (J[0][2]*J[1][0]-J[0][0]*J[1][2])*id;
>
> IJ[2][0] = (J[1][0]*J[2][1]-J[1][1]*J[2][0])*id;
> IJ[2][1] = (J[0][1]*J[2][0]-J[0][0]*J[2][1])*id;
> IJ[2][2] = (J[0][0]*J[1][1]-J[0][1]*J[1][0])*id;
>}
>
>static void getdIM(REAL M[4][3], REAL ix00, REAL ix10, REAL ix20, REAL
ix30, REAL IM[4][3], REAL dIM[3][4][3])
>{
> REAL dt, id, id2, ddt;
> REAL t00,t10,t20,t30,t01,t11,t21,t31,t02,t12,t22,t32;
>
> dt = M[0][0]*(M[1][1]*M[2][2]-M[1][2]*M[2][1])+
> M[0][1]*(M[1][2]*M[2][0]-M[1][0]*M[2][2])+
> M[0][2]*(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
> id = 1.0/dt;
>
> t00 = M[1][1]*M[2][2]-M[1][2]*M[2][1];
> t01 = M[0][2]*M[2][1]-M[0][1]*M[2][2];
> t02 = M[0][1]*M[1][2]-M[0][2]*M[1][1];
>
> t10 = M[1][2]*M[2][0]-M[1][0]*M[2][2];
> t11 = M[0][0]*M[2][2]-M[0][2]*M[2][0];
> t12 = M[0][2]*M[1][0]-M[0][0]*M[1][2];
>
> t20 = M[1][0]*M[2][1]-M[1][1]*M[2][0];
> t21 = M[0][1]*M[2][0]-M[0][0]*M[2][1];
> t22 = M[0][0]*M[1][1]-M[0][1]*M[1][0];
>
> t30 = M[1][0]*(M[3][1]*M[2][2]-M[2][1]*M[3][2])+
> M[1][1]*(M[2][0]*M[3][2]-M[3][0]*M[2][2])+
> M[1][2]*(M[3][0]*M[2][1]-M[2][0]*M[3][1]);
> t31 = M[0][0]*(M[2][1]*M[3][2]-M[3][1]*M[2][2])+
> M[0][1]*(M[3][0]*M[2][2]-M[2][0]*M[3][2])+
> M[0][2]*(M[2][0]*M[3][1]-M[3][0]*M[2][1]);
> t32 = M[0][0]*(M[3][1]*M[1][2]-M[1][1]*M[3][2])+
> M[0][1]*(M[1][0]*M[3][2]-M[3][0]*M[1][2])+
> M[0][2]*(M[3][0]*M[1][1]-M[1][0]*M[3][1]);
>
> IM[0][0] = t00*id; IM[0][1] = t01*id; IM[0][2] = t02*id;
> IM[1][0] = t10*id; IM[1][1] = t11*id; IM[1][2] = t12*id;
> IM[2][0] = t20*id; IM[2][1] = t21*id; IM[2][2] = t22*id;
> IM[3][0] = t30*id; IM[3][1] = t31*id; IM[3][2] = t32*id;
>
> id2 = id*id;
>
>
> ddt = ix00*(M[1][1]*M[2][2]-M[1][2]*M[2][1])
> +ix10*(M[0][2]*M[2][1]-M[0][1]*M[2][2])
> +ix20*(M[0][1]*M[1][2]-M[0][2]*M[1][1]);
>
> dIM[0][0][0] = -t00*id2*ddt;
> dIM[0][0][1] = -t01*id2*ddt;
> dIM[0][0][2] = -t02*id2*ddt;
> dIM[0][1][0] = (M[1][2]*ix20-ix10*M[2][2])*id - t10*id2*ddt;
> dIM[0][1][1] = (ix00*M[2][2]-M[0][2]*ix20)*id - t11*id2*ddt;
> dIM[0][1][2] = (M[0][2]*ix10-ix00*M[1][2])*id - t12*id2*ddt;
> dIM[0][2][0] = (ix10*M[2][1]-M[1][1]*ix20)*id - t20*id2*ddt;
> dIM[0][2][1] = (M[0][1]*ix20-ix00*M[2][1])*id - t21*id2*ddt;
> dIM[0][2][2] = (ix00*M[1][1]-M[0][1]*ix10)*id - t22*id2*ddt;
> dIM[0][3][0] = (ix10*(M[3][1]*M[2][2]-M[2][1]*M[3][2])+
> M[1][1]*(ix20*M[3][2]-ix30*M[2][2])+
> M[1][2]*(ix30*M[2][1]-ix20*M[3][1]))*id -
t30*id2*ddt;
> dIM[0][3][1] = (ix00*(M[2][1]*M[3][2]-M[3][1]*M[2][2])+
> M[0][1]*(ix30*M[2][2]-ix20*M[3][2])+
> M[0][2]*(ix20*M[3][1]-ix30*M[2][1]))*id -
t31*id2*ddt;
> dIM[0][3][2] = (ix00*(M[3][1]*M[1][2]-M[1][1]*M[3][2])+
> M[0][1]*(ix10*M[3][2]-ix30*M[1][2])+
> M[0][2]*(ix30*M[1][1]-ix10*M[3][1]))*id -
t32*id2*ddt;
>
>
> ddt = ix00*(M[1][2]*M[2][0]-M[1][0]*M[2][2])
> +ix10*(M[0][0]*M[2][2]-M[0][2]*M[2][0])
> +ix20*(M[0][2]*M[1][0]-M[0][0]*M[1][2]);
>
> dIM[1][0][0] = (ix10*M[2][2]-M[1][2]*ix20)*id - t00*id2*ddt;
> dIM[1][0][1] = (M[0][2]*ix20-ix00*M[2][2])*id - t01*id2*ddt;
> dIM[1][0][2] = (ix00*M[1][2]-M[0][2]*ix10)*id - t02*id2*ddt;
> dIM[1][1][0] = -t10*id2*ddt;
> dIM[1][1][1] = -t11*id2*ddt;
> dIM[1][1][2] = -t12*id2*ddt;
> dIM[1][2][0] = (M[1][0]*ix20-ix10*M[2][0])*id - t20*id2*ddt;
> dIM[1][2][1] = (ix00*M[2][0]-M[0][0]*ix20)*id - t21*id2*ddt;
> dIM[1][2][2] = (M[0][0]*ix10-ix00*M[1][0])*id - t22*id2*ddt;
> dIM[1][3][0] = (M[1][0]*(ix30*M[2][2]-ix20*M[3][2])+
> ix10*(M[2][0]*M[3][2]-M[3][0]*M[2][2])+
> M[1][2]*(M[3][0]*ix20-M[2][0]*ix30))*id -
t30*id2*ddt;
> dIM[1][3][1] = (M[0][0]*(ix20*M[3][2]-ix30*M[2][2])+
> ix00*(M[3][0]*M[2][2]-M[2][0]*M[3][2])+
> M[0][2]*(M[2][0]*ix30-M[3][0]*ix20))*id -
t31*id2*ddt;
> dIM[1][3][2] = (M[0][0]*(ix30*M[1][2]-ix10*M[3][2])+
> ix00*(M[1][0]*M[3][2]-M[3][0]*M[1][2])+
> M[0][2]*(M[3][0]*ix10-M[1][0]*ix30))*id -
t32*id2*ddt;
>
>
> ddt = ix00*(M[1][0]*M[2][1]-M[1][1]*M[2][0])
> +ix10*(M[0][1]*M[2][0]-M[0][0]*M[2][1])
> +ix20*(M[0][0]*M[1][1]-M[0][1]*M[1][0]);
>
> dIM[2][0][0] = (M[1][1]*ix20-ix10*M[2][1])*id - t00*id2*ddt;
> dIM[2][0][1] = (ix00*M[2][1]-M[0][1]*ix20)*id - t01*id2*ddt;
> dIM[2][0][2] = (M[0][1]*ix10-ix00*M[1][1])*id - t02*id2*ddt;
> dIM[2][1][0] = (ix10*M[2][0]-M[1][0]*ix20)*id - t10*id2*ddt;
> dIM[2][1][1] = (M[0][0]*ix20-ix00*M[2][0])*id - t11*id2*ddt;
> dIM[2][1][2] = (ix00*M[1][0]-M[0][0]*ix10)*id - t12*id2*ddt;
> dIM[2][2][0] = -t20*id2*ddt;
> dIM[2][2][1] = -t21*id2*ddt;
> dIM[2][2][2] = -t22*id2*ddt;
> dIM[2][3][0] = (M[1][0]*(M[3][1]*ix20-M[2][1]*ix30)+
> M[1][1]*(M[2][0]*ix30-M[3][0]*ix20)+
> ix10*(M[3][0]*M[2][1]-M[2][0]*M[3][1]))*id -
t30*id2*ddt;
> dIM[2][3][1] = (M[0][0]*(M[2][1]*ix30-M[3][1]*ix20)+
> M[0][1]*(M[3][0]*ix20-M[2][0]*ix30)+
> ix00*(M[2][0]*M[3][1]-M[3][0]*M[2][1]))*id -
t31*id2*ddt;
> dIM[2][3][2] = (M[0][0]*(M[3][1]*ix10-M[1][1]*ix30)+
> M[0][1]*(M[1][0]*ix30-M[3][0]*ix10)+
> ix00*(M[3][0]*M[1][1]-M[1][0]*M[3][1]))*id -
t32*id2*ddt;
>}
>
>static void getdIMpartial(REAL M[4][3], REAL ix00, REAL ix10, REAL ix20,
REAL IM[4][3], REAL dIM[3][4][3])
>{
> REAL dt, id, id2, ddt;
> REAL t00,t10,t20,t01,t11,t21,t02,t12,t22;
>
> dt = M[0][0]*(M[1][1]*M[2][2]-M[1][2]*M[2][1])+
> M[0][1]*(M[1][2]*M[2][0]-M[1][0]*M[2][2])+
> M[0][2]*(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
> id = 1.0/dt;
>
> t00 = M[1][1]*M[2][2]-M[1][2]*M[2][1];
> t01 = M[0][2]*M[2][1]-M[0][1]*M[2][2];
> t02 = M[0][1]*M[1][2]-M[0][2]*M[1][1];
>
> t10 = M[1][2]*M[2][0]-M[1][0]*M[2][2];
> t11 = M[0][0]*M[2][2]-M[0][2]*M[2][0];
> t12 = M[0][2]*M[1][0]-M[0][0]*M[1][2];
>
> t20 = M[1][0]*M[2][1]-M[1][1]*M[2][0];
> t21 = M[0][1]*M[2][0]-M[0][0]*M[2][1];
> t22 = M[0][0]*M[1][1]-M[0][1]*M[1][0];
>
> IM[0][0] = t00*id; IM[0][1] = t01*id; IM[0][2] = t02*id;
> IM[1][0] = t10*id; IM[1][1] = t11*id; IM[1][2] = t12*id;
> IM[2][0] = t20*id; IM[2][1] = t21*id; IM[2][2] = t22*id;
>
> id2 = id*id;
>
>
> ddt = ix00*(M[1][1]*M[2][2]-M[1][2]*M[2][1])
> +ix10*(M[0][2]*M[2][1]-M[0][1]*M[2][2])
> +ix20*(M[0][1]*M[1][2]-M[0][2]*M[1][1]);
>
> dIM[0][0][0] = -t00*id2*ddt;
> dIM[0][0][1] = -t01*id2*ddt;
> dIM[0][0][2] = -t02*id2*ddt;
> dIM[0][1][0] = (M[1][2]*ix20-ix10*M[2][2])*id - t10*id2*ddt;
> dIM[0][1][1] = (ix00*M[2][2]-M[0][2]*ix20)*id - t11*id2*ddt;
> dIM[0][1][2] = (M[0][2]*ix10-ix00*M[1][2])*id - t12*id2*ddt;
> dIM[0][2][0] = (ix10*M[2][1]-M[1][1]*ix20)*id - t20*id2*ddt;
> dIM[0][2][1] = (M[0][1]*ix20-ix00*M[2][1])*id - t21*id2*ddt;
> dIM[0][2][2] = (ix00*M[1][1]-M[0][1]*ix10)*id - t22*id2*ddt;
>
> ddt = ix00*(M[1][2]*M[2][0]-M[1][0]*M[2][2])
> +ix10*(M[0][0]*M[2][2]-M[0][2]*M[2][0])
> +ix20*(M[0][2]*M[1][0]-M[0][0]*M[1][2]);
>
> dIM[1][0][0] = (ix10*M[2][2]-M[1][2]*ix20)*id - t00*id2*ddt;
> dIM[1][0][1] = (M[0][2]*ix20-ix00*M[2][2])*id - t01*id2*ddt;
> dIM[1][0][2] = (ix00*M[1][2]-M[0][2]*ix10)*id - t02*id2*ddt;
> dIM[1][1][0] = -t10*id2*ddt;
> dIM[1][1][1] = -t11*id2*ddt;
> dIM[1][1][2] = -t12*id2*ddt;
> dIM[1][2][0] = (M[1][0]*ix20-ix10*M[2][0])*id - t20*id2*ddt;
> dIM[1][2][1] = (ix00*M[2][0]-M[0][0]*ix20)*id - t21*id2*ddt;
> dIM[1][2][2] = (M[0][0]*ix10-ix00*M[1][0])*id - t22*id2*ddt;
>
> ddt = ix00*(M[1][0]*M[2][1]-M[1][1]*M[2][0])
> +ix10*(M[0][1]*M[2][0]-M[0][0]*M[2][1])
> +ix20*(M[0][0]*M[1][1]-M[0][1]*M[1][0]);
>
> dIM[2][0][0] = (M[1][1]*ix20-ix10*M[2][1])*id - t00*id2*ddt;
> dIM[2][0][1] = (ix00*M[2][1]-M[0][1]*ix20)*id - t01*id2*ddt;
> dIM[2][0][2] = (M[0][1]*ix10-ix00*M[1][1])*id - t02*id2*ddt;
> dIM[2][1][0] = (ix10*M[2][0]-M[1][0]*ix20)*id - t10*id2*ddt;
> dIM[2][1][1] = (M[0][0]*ix20-ix00*M[2][0])*id - t11*id2*ddt;
> dIM[2][1][2] = (ix00*M[1][0]-M[0][0]*ix10)*id - t12*id2*ddt;
> dIM[2][2][0] = -t20*id2*ddt;
> dIM[2][2][1] = -t21*id2*ddt;
> dIM[2][2][2] = -t22*id2*ddt;
>}
>
>static void detJ(REAL J[3][3], REAL *d)
>{
> *d = J[0][0]*(J[1][1]*J[2][2]-J[1][2]*J[2][1])+
> J[1][0]*(J[0][2]*J[2][1]-J[0][1]*J[2][2])+
> J[2][0]*(J[0][1]*J[1][2]-J[0][2]*J[1][1]);
>}
>
>static void derivdetJ(REAL J[3][3], REAL scaleJ[3][3], REAL ix00, REAL
ix10, REAL ix20,
> REAL *d, REAL *dd0, REAL *dd1, REAL *dd2)
>{
> *d = J[0][0]*(J[1][1]*J[2][2]-J[1][2]*J[2][1])+
> J[1][0]*(J[0][2]*J[2][1]-J[0][1]*J[2][2])+
> J[2][0]*(J[0][1]*J[1][2]-J[0][2]*J[1][1]);
>
> *dd0 = ix00*scaleJ[0][0]*(J[1][1]*J[2][2]-J[1][2]*J[2][1])+
> ix10*scaleJ[1][0]*(J[0][2]*J[2][1]-J[0][1]*J[2][2])+
> ix20*scaleJ[2][0]*(J[0][1]*J[1][2]-J[0][2]*J[1][1]);
>
> *dd1 = ix00*scaleJ[0][1]*(J[1][2]*J[2][0]-J[1][0]*J[2][2])+
> ix10*scaleJ[1][1]*(J[0][0]*J[2][2]-J[0][2]*J[2][0])+
> ix20*scaleJ[2][1]*(J[0][2]*J[1][0]-J[0][0]*J[1][2]);
>
> *dd2 = ix00*scaleJ[0][2]*(J[1][0]*J[2][1]-J[1][1]*J[2][0])+
> ix10*scaleJ[1][2]*(J[0][1]*J[2][0]-J[0][0]*J[2][1])+
> ix20*scaleJ[2][2]*(J[0][0]*J[1][1]-J[0][1]*J[1][0]);
>}
>
>static REAL sumJ2(REAL J[4][3])
>{
> return(J[0][0]*J[0][0]+J[0][1]*J[0][1]+J[0][2]*J[0][2]+
> J[1][0]*J[1][0]+J[1][1]*J[1][1]+J[1][2]*J[1][2]+
> J[2][0]*J[2][0]+J[2][1]*J[2][1]+J[2][2]*J[2][2]);
>}
>
>static void derivsumJ2(REAL J[3][3], REAL scaleJ[3][3], REAL ix00, REAL
ix10, REAL ix20,
> REAL *d0, REAL *d1, REAL *d2)
>{
> *d0 = 2*(J[0][0]*ix00*scaleJ[0][0] + J[1][0]*ix10*scaleJ[1][0] + J
[2][0]*ix20*scaleJ[2][0]);
> *d1 = 2*(J[0][1]*ix00*scaleJ[0][1] + J[1][1]*ix10*scaleJ[1][1] + J
[2][1]*ix20*scaleJ[2][1]);
> *d2 = 2*(J[0][2]*ix00*scaleJ[0][2] + J[1][2]*ix10*scaleJ[1][2] + J
[2][2]*ix20*scaleJ[2][2]);
>}
>
>static void derivsumIJ2(REAL IJ[3][3], REAL scaleJ[3][3], REAL dIM[3][4]
[3], REAL *d0, REAL *d1, REAL *d2)
>{
> *d0 = 2*(IJ[0][0]*dIM[0][0][0]/scaleJ[0][0] + IJ[0][1]*dIM[0][0]
[1]/scaleJ[1][0] + IJ[0][2]*dIM[0][0][2]/scaleJ[2][0] +
> IJ[1][0]*dIM[0][1][0]/scaleJ[0][1] + IJ[1][1]*dIM[0][1]
[1]/scaleJ[1][1] + IJ[1][2]*dIM[0][1][2]/scaleJ[2][1] +
> IJ[2][0]*dIM[0][2][0]/scaleJ[0][2] + IJ[2][1]*dIM[0][2]
[1]/scaleJ[1][2] + IJ[2][2]*dIM[0][2][2]/scaleJ[2][2]);
> *d1 = 2*(IJ[0][0]*dIM[1][0][0]/scaleJ[0][0] + IJ[0][1]*dIM[1][0]
[1]/scaleJ[1][0] + IJ[0][2]*dIM[1][0][2]/scaleJ[2][0] +
> IJ[1][0]*dIM[1][1][0]/scaleJ[0][1] + IJ[1][1]*dIM[1][1]
[1]/scaleJ[1][1] + IJ[1][2]*dIM[1][1][2]/scaleJ[2][1] +
> IJ[2][0]*dIM[1][2][0]/scaleJ[0][2] + IJ[2][1]*dIM[1][2]
[1]/scaleJ[1][2] + IJ[2][2]*dIM[1][2][2]/scaleJ[2][2]);
> *d2 = 2*(IJ[0][0]*dIM[2][0][0]/scaleJ[0][0] + IJ[0][1]*dIM[2][0]
[1]/scaleJ[1][0] + IJ[0][2]*dIM[2][0][2]/scaleJ[2][0] +
> IJ[1][0]*dIM[2][1][0]/scaleJ[0][1] + IJ[1][1]*dIM[2][1]
[1]/scaleJ[1][1] + IJ[1][2]*dIM[2][1][2]/scaleJ[2][1] +
> IJ[2][0]*dIM[2][2][0]/scaleJ[0][2] + IJ[2][1]*dIM[2][2]
[1]/scaleJ[1][2] + IJ[2][2]*dIM[2][2][2]/scaleJ[2][2]);
>}
>
>
>void IMtimesY(REAL IM[4][3], int y[3], float x[3])
>{
> x[0] = IM[0][0]*y[0] + IM[1][0]*y[1] + IM[2][0]*y[2] + IM[3][0];
> x[1] = IM[0][1]*y[0] + IM[1][1]*y[1] + IM[2][1]*y[2] + IM[3][1];
> x[2] = IM[0][2]*y[0] + IM[1][2]*y[1] + IM[2][2]*y[2] + IM[3][2];
>}
>
>
>static void getJ(REAL M[4][3], REAL scaleJ[3][3], REAL J[3][3])
>{
> J[0][0] = M[0][0]*scaleJ[0][0];
> J[0][1] = M[0][1]*scaleJ[0][1];
> J[0][2] = M[0][2]*scaleJ[0][2];
> J[1][0] = M[1][0]*scaleJ[1][0];
> J[1][1] = M[1][1]*scaleJ[1][1];
> J[1][2] = M[1][2]*scaleJ[1][2];
> J[2][0] = M[2][0]*scaleJ[2][0];
> J[2][1] = M[2][1]*scaleJ[2][1];
> J[2][2] = M[2][2]*scaleJ[2][2];
>}
>
>static void getIJ(REAL IM[4][3], REAL scaleJ[3][3], REAL IJ[3][3])
>{
> IJ[0][0] = IM[0][0]/scaleJ[0][0];
> IJ[0][1] = IM[0][1]/scaleJ[1][0];
> IJ[0][2] = IM[0][2]/scaleJ[2][0];
> IJ[1][0] = IM[1][0]/scaleJ[0][1];
> IJ[1][1] = IM[1][1]/scaleJ[1][1];
> IJ[1][2] = IM[1][2]/scaleJ[2][1];
> IJ[2][0] = IM[2][0]/scaleJ[0][2];
> IJ[2][1] = IM[2][1]/scaleJ[1][2];
> IJ[2][2] = IM[2][2]/scaleJ[2][2];
>}
>
>#include "tets.h"
>
>void fun_ss_im(int x0, int x1, int x2, Tettype *tet, REAL Y[4][3], REAL
scaleJ[3][3], REAL IM[4][3], REAL *hp)
>{
> static REAL M[4][3], J[3][3], IJ[3][3];
> REAL ix30, dt;
> getM(Y, tet->ix, tet->dix, M, x0, x1, x2, &ix30);
> invertM(M, IM, &dt);
> getJ(M, scaleJ, J);
> getIJ(IM, scaleJ, IJ);
> detJ(J, &dt);
> *hp += 0.25*tet->volx*(1+dt)*(sumJ2(J) + sumJ2(IJ) - 6.0);
>}
>
>void dfun_ss_im(int x0, int x1, int x2, Tettype *tet, REAL Y[4][3], REAL
scaleJ[3][3],
> REAL IM[4][3], REAL dIM[3][4][3], REAL *hp, REAL *dhp0, REAL
*dhp1, REAL *dhp2)
>{
> static REAL M[4][3], ix30, J[3][3], IJ[3][3];
> REAL ss, dss0f, dss1f, dss2f, dss0b, dss1b, dss2b, dt, dd0, dd1,
dd2;
>
> getM(Y, tet->ix, tet->dix, M, x0, x1, x2, &ix30);
> getdIM(M, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0], ix30, IM,
dIM);
> getJ ( M, scaleJ, J);
> getIJ(IM, scaleJ, IJ);
>
> derivdetJ(J, scaleJ, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0],
&dt, &dd0, &dd1, &dd2);
> ss = sumJ2(J) + sumJ2(IJ) - 6.0;
> derivsumJ2(J, scaleJ, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0],
&dss0f, &dss1f, &dss2f);
> derivsumIJ2(IJ, scaleJ, dIM, &dss0b, &dss1b, &dss2b);
>
> *hp += 0.25*tet->volx * (1+dt)*ss;
> *dhp0 += 0.25*tet->volx *((1+dt)*(dss0f+dss0b)+dd0*ss);
> *dhp1 += 0.25*tet->volx *((1+dt)*(dss1f+dss1b)+dd1*ss);
> *dhp2 += 0.25*tet->volx *((1+dt)*(dss2f+dss2b)+dd2*ss);
>}
>
>void fun_ss(Tettype *tet, REAL Y[4][3], REAL scaleJ[3][3], REAL *hp)
>{
> static REAL M[4][3], J[3][3], IJ[3][3];
> REAL dt;
> getMpartial(Y, tet->ix, M);
> getJ(M, scaleJ, J);
> invertJ(J, IJ, &dt);
> *hp += 0.25*tet->volx*(1+dt)*(sumJ2(J) + sumJ2(IJ) - 6.0);
>}
>
>void dfun_ss(Tettype *tet, REAL Y[4][3], REAL scaleJ[3][3],
> REAL *hp, REAL *dhp0, REAL *dhp1, REAL *dhp2)
>{
> static REAL M[4][3], IM[4][3], dIM[3][4][3], J[3][3], IJ[3][3];
> REAL ss, dss0f, dss1f, dss2f, dss0b, dss1b, dss2b, dt, dd0, dd1,
dd2;
>
> getMpartial(Y, tet->ix, M);
> getdIMpartial(M, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0], IM,
dIM);
> getJ ( M, scaleJ, J);
> getIJ(IM, scaleJ, IJ);
>
> derivdetJ(J, scaleJ, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0],
&dt, &dd0, &dd1, &dd2);
> ss = sumJ2(J) + sumJ2(IJ) - 6.0;
> derivsumJ2(J, scaleJ, tet->ix[0][0], tet->ix[1][0], tet->ix[2][0],
&dss0f, &dss1f, &dss2f);
> derivsumIJ2(IJ, scaleJ, dIM, &dss0b, &dss1b, &dss2b);
>
> *hp += 0.25*tet->volx * (1+dt)*ss;
> *dhp0 += 0.25*tet->volx *((1+dt)*(dss0f+dss0b)+dd0*ss);
> *dhp1 += 0.25*tet->volx *((1+dt)*(dss1f+dss1b)+dd1*ss);
> *dhp2 += 0.25*tet->volx *((1+dt)*(dss2f+dss2b)+dd2*ss);
>}
>----------
>X-Sun-Data-Type: c-file
>X-Sun-Data-Description: c-file
>X-Sun-Data-Name: warpio.c
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 68
>
>#ifndef lint
>static char sccsid[] = "%W% John Ashburner %E%";
>#endif
>#include<stdio.h>
>
>int read_warp_plane(char *fname, int dim_g[3], float y[])
>{
> FILE *fp;
> fp = fopen(fname,"r");
> if (fp == 0)
> {
> return(1);
> }
> (void)printf("reading warp - %s ", fname);(void)fflush(stdout);
> if (fread(y, sizeof(float), dim_g[0]*dim_g[1]*dim_g[2], fp) != dim_g
[0]*dim_g[1]*dim_g[2])
> {
> (void)fprintf(stderr,"%s: can't read warp data.\n", fname);
> (void)fclose(fp);
> return(1);
> }
> (void)fclose(fp);
> (void)printf("...done\n");(void)fflush(stdout);
> return(0);
>}
>
>int read_warp(char *fname, int dim_g[3], float y[])
>{
> FILE *fp;
> fp = fopen(fname,"r");
> if (fp == 0)
> {
> return(1);
> }
> (void)printf("reading warp - %s ", fname);(void)fflush(stdout);
> if (fread(y, sizeof(float), dim_g[0]*dim_g[1]*dim_g[2], fp) != dim_g
[0]*dim_g[1]*dim_g[2])
> {
> (void)fprintf(stderr,"%s: can't read warp data.\n", fname);
> (void)fclose(fp);
> return(1);
> }
> (void)fclose(fp);
> (void)printf("...done\n");(void)fflush(stdout);
> return(0);
>}
>
>int write_warp(char *fname, int dim_g[3], float y[])
>{
> FILE *fp;
>
> (void)printf("writing warp - %s ", fname);(void)fflush(stdout);
>
> fp = fopen(fname,"w");
> if (fp == 0)
> {
> (void)perror(fname);
> return(1);
> }
> if (fwrite(y, sizeof(float), dim_g[0]*dim_g[1]*dim_g[2], fp) !=
dim_g[0]*dim_g[1]*dim_g[2])
> {
> (void)fprintf(stderr,"%s: can't write warp data.\n", fname);
> (void)fclose(fp);
> return(1);
> }
> (void)fclose(fp);
> (void)printf("...done\n");(void)fflush(stdout);
> return(0);
>}
>
>----------
>X-Sun-Data-Type: default
>X-Sun-Data-Description: default
>X-Sun-Data-Name: spm_write_defs.m
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 108
>
>function img = spm_write_defs(matname, vox,bb)
>% Write deformation field.
>% FORMAT spm_write_defs(matname, vox,bb)
>% matname - the `_sn3d.mat' file containing the spatial normalization
>% parameters.
>% The deformations are stored in y1.img, y2.img and y3.img
>%_______________________________________________________________________
>% %W% John Ashburner %E%
>
>matname = deblank(matname);
>
>load(matname)
>if (exist('mgc') ~= 1)
> error(['Matrix file ' matname ' is the wrong type.']);
>end
>if (mgc ~= 960209)
> error(['Matrix file ' matname ' is the wrong type.']);
>end
>
>if nargin>=3,
> x = (bb(1,1):vox(1):bb(2,1))/Dims(3,1) + Dims(4,1);
> y = (bb(1,2):vox(2):bb(2,2))/Dims(3,2) + Dims(4,2);
> z = (bb(1,3):vox(3):bb(2,3))/Dims(3,3) + Dims(4,3);
>
> dim = [length(x) length(y) length(z)];
>
> origin = round(-bb(1,:)./vox + 1);
> off = -vox.*origin;
> mat = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off
(3) ; 0 0 0 1];
>else,
> dim = Dims(1,:);
> x = 1:dim(1);
> y = 1:dim(2);
> z = 1:dim(3);
>
> mat = MG;
>end;
>
>VX = struct('fname',[matname(1:end-9) '_y1.img'], 'dim',[dim
16], 'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - X');
>VY = struct('fname',[matname(1:end-9) '_y2.img'], 'dim',[dim
16], 'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - Y');
>VZ = struct('fname',[matname(1:end-9) '_y3.img'], 'dim',[dim
16], 'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - Z');
>
>
>X = x'*ones(1,VX.dim(2));
>Y = ones(VX.dim(1),1)*y;
>
>
>if (prod(Dims(2,:)) == 0),
> affine_only = 1;
> basX = 0; tx = 0;
> basY = 0; ty = 0;
> basZ = 0; tz = 0;
> str = 'affine';
>else
> affine_only = 0;
> basX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1);
> basY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1);
> basZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1);
> str = 'nonlinear';
>end
>
>% Start progress plot
>%--------------------------------------------------------------------------
--
>spm_progress_bar('Init',length(z),['Writing ' str ' deformations'],'planes
completed');
>
>spm_create_image(VX);
>spm_create_image(VY);
>spm_create_image(VZ);
>
>% Cycle over planes
>%--------------------------------------------------------------------------
--
>for j=1:length(z)
>
> % Nonlinear deformations
> %-------------------------------------------------------------------
---------
> if (~affine_only)
> % 2D transforms for each plane
> tx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims
(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
> ty = reshape( reshape(Transform(:,2),Dims(2,1)*Dims
(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
> tz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims
(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
>
> X1 = X + basX*tx*basY';
> Y1 = Y + basX*ty*basY';
> Z1 = z(j) + basX*tz*basY';
> end
>
> % Sample each volume
> %-------------------------------------------------------------------
---------
> Mult = Affine;
> if (~affine_only)
> X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
> Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
> Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
> else
> X2= Mult(1,1)*X + Mult(1,2)*Y + (Mult(1,3)*z(j) + Mult
(1,4));
> Y2= Mult(2,1)*X + Mult(2,2)*Y + (Mult(2,3)*z(j) + Mult
(2,4));
> Z2= Mult(3,1)*X + Mult(3,2)*Y + (Mult(3,3)*z(j) + Mult
(3,4));
> end
>
> spm_write_plane(VX,X2-1,j);
> spm_write_plane(VY,Y2-1,j);
> spm_write_plane(VZ,Z2-1,j);
>
> spm_progress_bar('Set',j);
>end
>spm_progress_bar('Clear');
>return;
>%_______________________________________________________________________
>----------
>X-Sun-Data-Type: default
>X-Sun-Data-Description: default
>X-Sun-Data-Name: spm_dctmtx.m
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 59
>
>function C = spm_dctmtx(N,K,n,f)
>% Creates basis functions for Discrete Cosine Transform.
>% FORMAT C = spm_dctmtx(N,K,n)
>% OR C = spm_dctmtx(N,K)
>% OR D = spm_dctmtx(N,K,n,'diff')
>% OR D = spm_dctmtx(N,K,'diff')
>% N - dimension
>% K - order
>% n - optional points to sample
>%
____________________________________________________________________________
>% spm_dctmtx creates a matrix for the first few basis functions of a one
>% dimensional discrete cosine transform.
>% With the 'diff' argument, spm_dctmtx produces the derivatives of the
>% DCT.
>%
>% See: Fundamentals of Digital Image Processing (p 150-154).
>% Anil K. Jain 1989.
>%
____________________________________________________________________________
>% @(#)spm_dctmtx.m 1.3 John Ashburner MRCCU/FIL 96/08/14
>
>d = 0;
>
>if (nargin == 2)
> n = (0:(N-1))';
> if (nargin == 3)
> d = 1;
> end
>elseif (nargin == 3)
> if (strcmp(n,'diff'))
> d = 1;
> n = (0:(N-1))';
> else
> n = n(:);
> end
>elseif (nargin == 4)
> n = n(:);
> if (strcmp(f,'diff'))
> d = 1;
> else
> error('Incorrect Usage');
> end
>else
> error('Incorrect Usage');
>end
>
>C = zeros(size(n,1),K);
>
>if (d == 0)
> C(:,1)=ones(size(n,1),1)/sqrt(N);
> for k=2:K
> C(:,k) = sqrt(2/N)*cos(pi*(2*n+1)*(k-1)/(2*N));
> end
>else
> for k=2:K
> C(:,k) = -2^(1/2)*(1/N)^(1/2)*sin(1/2*pi*(2*n*k-2*n+k-1)/N)
*pi*(k-1)/N;
> end
>end
>
>
>----------
>X-Sun-Data-Type: default
>X-Sun-Data-Description: default
>X-Sun-Data-Name: spm_write_plane.m
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 150
>
>function VO = spm_write_plane(V,A,p)
>% Write a transverse plane of image data.
>% FORMAT spm_write_plane(V,A,p)
>% V - data structure containing image information.
>% - see spm_vol for a description.
>% A - the two dimensional image to write.
>% p - the plane number (beginning from 1).
>%
>% VO - (possibly) modified data structure containing image information.
>% It is possible that future versions of spm_write_plane may
>% modify scalefactors (for example).
>%
>%_______________________________________________________________________
>% @(#)spm_write_plane.m 2.5 John Ashburner 99/03/31
>
>if any(V.dim(1:2) ~= size(A))
> error('Incompatible image dimensions');
>end;
>
>% Write Analyze image by default
>VO = write_analyze_plane(V,A,p);
>if isempty(VO),
> spm_progress_bar('Clear');
> write_error_message(V.fname);
> error(['Error writing ' V.fname '. Check your disk space.']);
>end;
>return;
>%_______________________________________________________________________
>
>%_______________________________________________________________________
>function VO = write_analyze_plane(V,A,p)
>VO = V;
>% Check datatype is OK
>dt = V.dim(4);
>
>% Convert to native datatype
>if dt>256,
> dt = dt/256;
>end;
>s = find(dt == [2 4 8 16 64 128+2 128+4 128+8]);
>if isempty(s)
> VO = [];
> disp(['Unrecognised data type (' num2str(V.dim(4)) ')']);
> return;
>end;
>datasize = [1 2 4 4 8 1 2 4];
>datasize = datasize(s);
>
>% Open the image file
>fname = [spm_str_manip(V.fname,'svd') '.img'];
>fp = fopen(fname,'r+');
>if fp == -1,
> fp = fopen(fname,'w');
> if fp == -1,
> VO = [];
> disp('Can''t open image file');
> return;
> end;
>end;
>
>% Rescale to fit appropriate range
>if any(dt == [2 4 8 128+2 128+4 128+8]),
> msk = find(~finite(A));
> A(msk) = 0.0;
> if dt == 2,
> maxval = max(255*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/255;
> A = round(A/scale);
> A(find(A < 0)) = 0;
> A(find(A > 255)) = 255;
> elseif dt == 4,
> maxval = max(32767*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/32767;
> A = round(A/scale);
> A(find(A > 32767)) = 32767;
> A(find(A < -32768)) = -32768;
> elseif (dt == 8)
> maxval = max((2^31-1)*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/(2^31-1);
> A = round(A/scale);
> A(find(A > 2^31-1)) = 2^31-1;
> A(find(A < -2^31 )) = -2^31;
> elseif dt == 128+2,
> maxval = max(127*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/255;
> A = round(A/scale);
> A(find(A < 0)) = 0;
> A(find(A > 255)) = 255;
> dt = dt - 128;
> elseif dt == 128+4,
> maxval = max(65535*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/32767;
> A = round(A/scale);
> A(find(A > 32767)) = 32767;
> A(find(A < -32768)) = -32768;
> dt = dt - 128;
> elseif (dt == 128+8)
> maxval = max((2^32-1)*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/(2^31-1);
> A = round(A/scale);
> A(find(A > 2^31-1)) = 2^31-1;
> A(find(A < -2^31 )) = -2^31;
> dt = dt - 128;
> end;
>end;
>
>% Seek to the appropriate offset
>off = (p-1)*datasize*prod(V.dim(1:2)); % + V.pinfo(3,1);
>if fseek(fp,off,'bof')==-1,
> % Need this because fseek in Matlab does not
> % seek past the EOF
> fseek(fp,0,'eof');
> curr_off = ftell(fp);
> blanks = zeros(off-curr_off,1);
> if fwrite(fp,blanks,'uchar') ~= prod(size(blanks)),
> fclose(fp);
> VO = [];
> return;
> end;
> if fseek(fp,off,'bof') == -1,
> fclose(fp);
> VO = [];
> return;
> end;
>end;
>
>if fwrite(fp,A,spm_type(dt)) ~= prod(size(A)),
> fclose(fp);
> VO = [];
> return;
>end;
>fclose(fp);
>
>return;
>%_______________________________________________________________________
>
>%_______________________________________________________________________
>function write_error_message(q)
>
>str = {...
> 'Error opening:',...
> [' ',spm_str_manip(q,'k40d')],...
> ' ',...
> 'Check file permissions / disk space / disk quota.'};
>
>msgbox(str,sprintf('%s%s: %s...',spm('ver'),spm('GetUser',' (%
s)'),mfilename),...
> 'error')
>
>return;
>%_______________________________________________________________________
>----------
>X-Sun-Data-Type: default
>X-Sun-Data-Description: default
>X-Sun-Data-Name: spm_create_image.m
>X-Sun-Charset: us-ascii
>X-Sun-Content-Lines: 100
>
>function Vo=spm_create_image(Vi)
>% Create an image file.
>% FORMAT Vo = spm_create_image(Vi)
>% Vi - data structure containing image information.
>% - see spm_vol for a description.
>% Vo - data structure after modification for writing.
>%_______________________________________________________________________
>% @(#)spm_create_image.m 2.3 John Ashburner 98/08/14
>
>[sts,Vo] = create_analyze_image(Vi);
>if sts == -1,
> spm_progress_bar('Clear');
> open_error_message(Vi.fname);
> error(['Error opening ' Vi.fname '. Check that you have write
permission.']);
>end;
>return;
>%_______________________________________________________________________
>
>%_______________________________________________________________________
>function [sts,V] = create_analyze_image(V)
>sts = 0;
>
>% Extract voxel sizes and origin from the 4x4 matrix
>vx = sqrt(sum(V.mat(1:3,1:3).^2));
>orgn = V.mat\[0 0 0 1]';
>orgn = orgn(1:3)';
>
>% Add a description when available
>if ~isfield(V,'descrip'),
> V.descrip = 'SPM compatible';
>end;
>% Check datatype is OK
>dt = V.dim(4);
>
>% Convert to native datatype
>if dt>256,
> dt = dt/256;
>end;
>s = find(dt == [2 4 8 16 64 128+2 128+4 128+8]);
>if isempty(s)
> sts = -1;
> disp(['Unrecognised data type (' num2str(V.dim(4)) ')']);
> return;
>end;
>
>% Compute an appropriate scalefactor
>if dt == 2,
> maxval = max(255*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/255;
>elseif dt == 4,
> maxval = max(32767*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/32767;
>elseif (dt == 8)
> maxval = max((2^31-1)*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/(2^31-1);
>elseif dt == 128+2,
> maxval = max(127*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/255;
> dt = dt - 128;
>elseif dt == 128+4,
> maxval = max(65535*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/32767;
> dt = dt - 128;
>elseif (dt == 128+8)
> maxval = max((2^32-1)*V.pinfo(1,:) + V.pinfo(2,:));
> scale = maxval/(2^31-1);
> dt = dt - 128;
>else
> scale = 1.0;
>end;
>
>V.pinfo = [scale 0 0]';
>V.dim(4)= dt;
>
>% Write the header
>s = spm_hwrite(deblank(V.fname), [V.dim(1:3) 1],...
> vx, scale, dt, 0, orgn, V.descrip);
>if s~= 348,
> sts = -1;
>end;
>
>% Write the matrix
>spm_get_space(deblank(V.fname),V.mat);
>return;
>%_______________________________________________________________________
>
>%_______________________________________________________________________
>function open_error_message(q)
>f=spm_figure('findwin','Graphics');
>if ~isempty(f),
> figure(f);
> spm_figure('Clear','Graphics');
> spm_figure('Clear','Interactive');
> ax=axes('Visible','off','Parent',f);
> text(0,0.60,'Error opening:', 'FontSize',
25, 'Interpreter', 'none');
> text(0,0.55,spm_str_manip(q,'k40d'), 'FontSize',
25, 'Interpreter', 'none');
> text(0,0.40,' Please check that you have write
permission.', 'FontSize', 16, 'Interpreter', 'none');
>end
>return
>%_______________________________________________________________________
>
>
>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>=========================================================================
|
