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
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