Home > bioelectromagnetism > elec_coregister.m

elec_coregister

PURPOSE ^

fiducial_coregister - Determine fiducial points coregistration

SYNOPSIS ^

function [T] = fiducial_coregister(P,Q)

DESCRIPTION ^

 fiducial_coregister - Determine fiducial points coregistration
 
 [T] = fiducial_coregister(P,Q)
 
 P is [3x3], three points in Cartesian XYZ
 Q is [3x3], three points in Cartesian XYZ
 
 This function 3D rotates and translates P into
 the space of Q.  The method is described in
 this m file, developed from Mortenson, M. (1985),
 Geometric Modelling, New York: John Wiley & Sons,
 pp. 353-355.
 
 The return 4x4 matrix (T) contains a 3x3 rotation 
 matrix in the upper left and a 1x3 translation
 row vector in the bottom row. The last column contains
 a projection vector in T(1:3,4)(all zeros here) and 
 the scalar at T(4,4) is a homogeneous coordinate scale 
 factor, usually 1.
 
 Using this T matrix:
 to get P values in the space of Q,
 P2Q = [P,ones(size(P,1),1)] * T;
 to get Q values in the space of P,
 Q2P = [Q,ones(size(Q,1),1)] * inv(T);
 
 In this manner, the rotation and translation is
 applied to any Nx3 points in relation to P space.
 
 For use in coregistration of electrodes to a scalp
 mesh, P & Q are the fiducial points:
 P(1,:) is the electrode nasion XYZ
 P(2,:) is the electrode right preauricular XYZ
 P(3,:) is the electrode left preauricular XYZ
 Similarly:
 Q(1,:) is the scalp mesh nasion XYZ
 Q(2,:) is the scalp mesh right preauricular XYZ
 Q(3,:) is the scalp mesh left preauricular XYZ
 The order of these points is essential!
 
 For example, using emse_read_reg to obtain reg, we have:
 T = fiducial_coregister(reg.fiducials.head,reg.fiducials.mri);
 where T is almost equivalent to reg.elec2mri.  This
 function is a hard fiducial transformation, whereas 
 the EMSE coregistration is a bit more sophisticated,
 providing weights to the fiducial points and some MRI
 surface points.

 See also: MESH_FIT_ELEC

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [T] = fiducial_coregister(P,Q)
0002 
0003 % fiducial_coregister - Determine fiducial points coregistration
0004 %
0005 % [T] = fiducial_coregister(P,Q)
0006 %
0007 % P is [3x3], three points in Cartesian XYZ
0008 % Q is [3x3], three points in Cartesian XYZ
0009 %
0010 % This function 3D rotates and translates P into
0011 % the space of Q.  The method is described in
0012 % this m file, developed from Mortenson, M. (1985),
0013 % Geometric Modelling, New York: John Wiley & Sons,
0014 % pp. 353-355.
0015 %
0016 % The return 4x4 matrix (T) contains a 3x3 rotation
0017 % matrix in the upper left and a 1x3 translation
0018 % row vector in the bottom row. The last column contains
0019 % a projection vector in T(1:3,4)(all zeros here) and
0020 % the scalar at T(4,4) is a homogeneous coordinate scale
0021 % factor, usually 1.
0022 %
0023 % Using this T matrix:
0024 % to get P values in the space of Q,
0025 % P2Q = [P,ones(size(P,1),1)] * T;
0026 % to get Q values in the space of P,
0027 % Q2P = [Q,ones(size(Q,1),1)] * inv(T);
0028 %
0029 % In this manner, the rotation and translation is
0030 % applied to any Nx3 points in relation to P space.
0031 %
0032 % For use in coregistration of electrodes to a scalp
0033 % mesh, P & Q are the fiducial points:
0034 % P(1,:) is the electrode nasion XYZ
0035 % P(2,:) is the electrode right preauricular XYZ
0036 % P(3,:) is the electrode left preauricular XYZ
0037 % Similarly:
0038 % Q(1,:) is the scalp mesh nasion XYZ
0039 % Q(2,:) is the scalp mesh right preauricular XYZ
0040 % Q(3,:) is the scalp mesh left preauricular XYZ
0041 % The order of these points is essential!
0042 %
0043 % For example, using emse_read_reg to obtain reg, we have:
0044 % T = fiducial_coregister(reg.fiducials.head,reg.fiducials.mri);
0045 % where T is almost equivalent to reg.elec2mri.  This
0046 % function is a hard fiducial transformation, whereas
0047 % the EMSE coregistration is a bit more sophisticated,
0048 % providing weights to the fiducial points and some MRI
0049 % surface points.
0050 %
0051 % See also: MESH_FIT_ELEC
0052 %
0053 
0054 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:33 $
0055 
0056 % Licence:  GNU GPL, no express or implied warranties
0057 % History:  09/2002, Darren.Weber_at_radiology.ucsf.edu
0058 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0059 
0060 
0061 DB = 0; % debug, 0 off, 1 on
0062 
0063 
0064 if ~exist('P','var') | ~exist('Q','var'),
0065     msg = sprintf('ELEC_COREGISTER: no input P or Q\n');
0066     error(msg);
0067 elseif isempty(P) | isempty(Q),
0068     msg = sprintf('ELEC_COREGISTER: empty input P or Q\n');
0069     error(msg);
0070 elseif ~(isequal(size(P),size(ones(3,3)))),
0071     msg = sprintf('ELEC_COREGISTER: P must be 3x3 matrix\n');
0072     error(msg);
0073 elseif ~(isequal(size(Q),size(ones(3,3)))),
0074     msg = sprintf('ELEC_COREGISTER: Q must be 3x3 matrix\n');
0075     error(msg);
0076 end
0077 
0078 
0079 % Given a geometric model that contains 3D points p1,p2,p3
0080 % and three other points q1,q2,q3, find the total rigid
0081 % body transformation that:
0082 % (1) transforms p1 into q1;
0083 % (2) transforms the vector (p2 - p1) into the
0084 %     vector (q2 - q1), direction only;
0085 % (3) transforms the plane containing the three points
0086 %     p1,p2,p3 into the plane containing q1,q2,q3
0087 
0088 % Note that when p1,p2,p3 are reference points among
0089 % many points of a geometric model, the transformations
0090 % here can be applied to all points to bring the model
0091 % into coincidence with the points in q1,q2,q3
0092 
0093 if DB,
0094     % These 2 point sets can be coregistered with just
0095     % a rotation of 90 degrees around Z
0096     P = [0 0 0; 1 0 0;  0 1 0];
0097     Q = [0 0 0; 0 1 0; -1 0 0];
0098     P = P + 10;
0099 end
0100 
0101 % Step 2 of Mortenson (1985, p. 354)
0102 V1 = P(2,:) - P(1,:);
0103 W1 = Q(2,:) - Q(1,:);
0104 
0105 % Step 3 of Mortenson (1985, p. 354)
0106 % V3 & W3 are orthogonal to the plane containing
0107 % the points p1,p3 and q1,q3 respectively
0108 V3 = cross( V1, [P(3,:) - P(1,:)] );
0109 W3 = cross( W1, [Q(3,:) - Q(1,:)] );
0110 
0111 % Step 4 of Mortenson (1985, p. 354)
0112 % V2 & W2 are orthogonal to the vectors
0113 % V1,V3 and W1,W3 respectively
0114 V2 = cross( V3, V1 );
0115 W2 = cross( W3, W1 );
0116 
0117 % Thus, V1,V2,V3 form a right hand orthogonal system,
0118 % as do W1,W2,W3
0119 
0120 % Step 5 of Mortenson (1985, p. 354)
0121 % Calculate unit vectors
0122 v1 = V1 / sqrt( V1(1)^2 + V1(2)^2 + V1(3)^2 );
0123 v2 = V2 / sqrt( V2(1)^2 + V2(2)^2 + V2(3)^2 );
0124 v3 = V3 / sqrt( V3(1)^2 + V3(2)^2 + V3(3)^2 );
0125 w1 = W1 / sqrt( W1(1)^2 + W1(2)^2 + W1(3)^2 );
0126 w2 = W2 / sqrt( W2(1)^2 + W2(2)^2 + W2(3)^2 );
0127 w3 = W3 / sqrt( W3(1)^2 + W3(2)^2 + W3(3)^2 );
0128 
0129 % Step 6 of Mortenson (1985, p. 354)
0130 % To transform any point p in the v system into the
0131 % w system, use the transformation relationship
0132 % p* = pR + T
0133 % Obviously we need to find R and T first
0134 
0135 % Step 7 of Mortenson (1985, p. 354)
0136 % [w1 w2 w3] = [v1 v2 v3]R, since [w] and [v] are
0137 % the unit vector matrices.
0138 w = [w1;w2;w3];
0139 v = [v1;v2;v3];
0140 % Then the required rotation matrix with respect
0141 % to the w system is simply:
0142 R = v\w;
0143 % Note in matlab A\B is the matrix division
0144 % of A into B, which is roughly the same as INV(A)*B,
0145 % except it is computed by Gaussian elimination.
0146 
0147 % Step 8 of Mortenson (1985, p. 354)
0148 % Obtain the translation matrix by substituting
0149 % R back into step 6 and solving for T, using
0150 % p1 and q1
0151 T = Q(1,:) - [ P(1,:) * R ];
0152 
0153 % Now construct the homogeneous coordinate
0154 % transformation matrix that encorporates both
0155 % of the above R & T transforms
0156 temp1 = [ [ R; zeros(1,3) ], zeros(4,1) ];
0157 temp1(4,4) = 1;
0158 temp1(4,1:3) = T;
0159 
0160 T = temp1;
0161 
0162 if DB,
0163     P2Q = [P,ones(3,1)] * T;
0164     plot3(P(:,1),P(:,2),P(:,3),'bo'); hold on; view(2)
0165     plot3(Q(:,1),Q(:,2),Q(:,3),'ro');
0166     plot3(P2Q(:,1),P2Q(:,2),P2Q(:,3),'b.');
0167     Q2P = [Q,ones(3,1)] * inv(T);
0168     plot3(Q2P(:,1),Q2P(:,2),Q2P(:,3),'r.');
0169 end
0170 
0171 
0172 return

Generated on Mon 15-Aug-2005 15:36:19 by m2html © 2003