Home > bioelectromagnetism > fiducial_coregister.m

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

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

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