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