Home > bioelectromagnetism > elec_sphere_refine.m

elec_sphere_refine

PURPOSE ^

elec_sphere_refine - creates spherical vertices between input vertices

SYNOPSIS ^

function [ p ] = elec_sphere_refine(p)

DESCRIPTION ^

 elec_sphere_refine - creates spherical vertices between input vertices

 [ p ] = elec_sphere_refine( p )

 p is the eeg_toolbox struct (see eeg_toolbox_defaults)
 
 This function works with the p.elec.data struct,
 especially the spherical projections created by 
 elec_sphere_project.
 
 The spherical electrode positions are triangulated with
 convhulln.  Then, for each face, 4 new vertices are created 
 at the triangle edge midpoints and the triangle centroid.  Each
 face is divided into 6 faces.
 
 The vertices created in the plane of each face are then projected
 to the spherical surface.

 See also mesh_refine, mesh_laplacian, mesh_laplacian_interp

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [ p ] = elec_sphere_refine(p)
0002 
0003 % elec_sphere_refine - creates spherical vertices between input vertices
0004 %
0005 % [ p ] = elec_sphere_refine( p )
0006 %
0007 % p is the eeg_toolbox struct (see eeg_toolbox_defaults)
0008 %
0009 % This function works with the p.elec.data struct,
0010 % especially the spherical projections created by
0011 % elec_sphere_project.
0012 %
0013 % The spherical electrode positions are triangulated with
0014 % convhulln.  Then, for each face, 4 new vertices are created
0015 % at the triangle edge midpoints and the triangle centroid.  Each
0016 % face is divided into 6 faces.
0017 %
0018 % The vertices created in the plane of each face are then projected
0019 % to the spherical surface.
0020 %
0021 % See also mesh_refine, mesh_laplacian, mesh_laplacian_interp
0022 %
0023 
0024 
0025 % This can be done until some minimal distance (D) of the mean
0026 % distance between vertices of all triangles is achieved.  If
0027 % no D argument is given, the function refines the mesh once.
0028 %
0029 
0030 % $Revision: 1.2 $ $Date: 2005/07/12 21:29:22 $
0031 
0032 % Licence:  GNU GPL, no implied or express warranties
0033 % History:  05/2002, Darren.Weber_at_radiology.ucsf.edu, created
0034 %
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 
0038 % Check the input struct
0039 if ~exist('p','var'),
0040    [p] = elec_open;
0041 end
0042 
0043 
0044 tic;
0045 fprintf('ELEC_SPHERE_REFINE: processing faces...');
0046 
0047 % NOTE
0048 % The centroid is located one third of the way from each vertex to
0049 % the midpoint of the opposite side. Each median divides the triangle
0050 % into two equal areas; all the medians together divide it into six
0051 % equal parts, and the lines from the median point to the vertices
0052 % divide the whole into three equivalent triangles.
0053 
0054 V = [p.elec.data.Xsp, p.elec.data.Ysp, p.elec.data.Zsp];
0055 F = convhulln(V);
0056 
0057 Zmin = min(V(:,3));
0058 
0059 % Initialise a new vertices and faces matrix
0060 Nvert = size(V,1);
0061 Nface = size(F,1);
0062 V2 = zeros(Nface*4,3);
0063 F2 = zeros(Nface*6,3);
0064 
0065 for r=1:size(F,1),
0066     
0067     % Get the triangle vertices
0068     A = V(F(r,1),:);
0069     B = V(F(r,2),:);
0070     C = V(F(r,3),:);
0071     
0072     % Now find the midpoint between all vertices
0073     ABmid = (A + B) ./ 2;
0074     BCmid = (B + C) ./ 2;
0075     CAmid = (C + A) ./ 2;
0076     
0077     % Now find the centroid length of the medians
0078     Ac = BCmid + ( (A - BCmid) ./3 );
0079     %Bc = CAmid + ( (B - CAmid) ./3 );  % Bc = Ac
0080     %Cc = ABmid + ( (C - ABmid) ./3 );  % Cc = Ac
0081     
0082     % Store the midpoints and the centroid vertices
0083     NABmid = Nvert + (r*4-3);
0084     V2(r*4-3,:) = ABmid;
0085     NBCmid = Nvert + (r*4-2);
0086     V2(r*4-2,:) = BCmid;
0087     NCAmid = Nvert + (r*4-1);
0088     V2(r*4-1,:) = CAmid;
0089     Nc     = Nvert + (r*4-0);
0090     V2(r*4-0,:) = Ac;
0091     
0092     % Create new faces for centroid
0093     F2(r*6-5,:) = [F(r,1), NABmid, Nc     ];
0094     F2(r*6-4,:) = [F(r,1), Nc    , NCAmid ];
0095     F2(r*6-3,:) = [F(r,2), Nc    , NABmid ];
0096     F2(r*6-2,:) = [F(r,2), NBCmid, Nc     ];
0097     F2(r*6-1,:) = [F(r,3), NCAmid, Nc     ];
0098     F2(r*6-0,:) = [F(r,3), Nc    , NBCmid ];
0099     
0100     %figure; patch('vertices',[A;B;C],'faces',[1 2 3],'facecolor',[.7 .7 .7]); hold on;
0101     %plot3(A(1),A(2),A(3),'ro');
0102     %plot3(BCmid(1),BCmid(2),BCmid(3),'ro');
0103     %plot3(Ac(1),Ac(2),Ac(3),'bo')
0104     %plot3(Bc(1),Bc(2),Bc(3),'bd')
0105     %plot3(Cc(1),Cc(2),Cc(3),'b.')
0106     %if isequal(r,2), return; end
0107     
0108 end
0109 
0110 % Add the new vertices to the old ones
0111 V = [V;V2];
0112 
0113 % All of the above is plane geometry,
0114 % So now project all points to the sphere
0115 V = proj_sph(V,p.elec.data.Rsp(1),p.elec.data.centroid);
0116 
0117 % Reset all vertices below Zmin to Zmin
0118 % This is a crude modification - needs work
0119 Vlow = find(V(:,3) < Zmin);
0120 V(Vlow,3) = Zmin;
0121 
0122 % Store the results in new elements of p.elec.data
0123 p.elec.data.Vsp = V;
0124 p.elec.data.Fsp = F2;
0125 
0126 
0127 t=toc;
0128 fprintf('...done (%5.2f sec)\n',t);
0129 
0130 return
0131 
0132 
0133 % Find the length of each median
0134 %A2BClen = sqrt ( sum( (A - BCmid).^2, 2 ) );
0135 %B2CAlen = sqrt ( sum( (B - CAmid).^2, 2 ) );
0136 %C2ABlen = sqrt ( sum( (C - ABmid).^2, 2 ) );
0137 
0138 
0139 function V = proj_sph(v,r,c)
0140     
0141     % Find the projection point of X,Y,Z to the fitted sphere radius r
0142     % Cartesian inputs:
0143     % v is the vertex matrix
0144     % r is the sphere radius
0145     % c is the sphere centroid
0146     
0147     X = v(:,1);
0148     Y = v(:,2);
0149     Z = v(:,3);
0150     
0151     xo = c(1);
0152     yo = c(2);
0153     zo = c(3);
0154         
0155     % Convert Cartesian X,Y,Z to spherical (radians)
0156     theta = atan2( (Y-yo), (X-xo) );
0157     phi = atan2( sqrt( (X-xo).^2 + (Y-yo).^2 ), (Z-zo) );
0158     % do not recalc: r = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2);
0159     
0160     %   Recalculate X,Y,Z for constant r, given theta & phi.
0161     R = ones(size(phi)) * r;    
0162     x = R .* sin(phi) .* cos(theta);
0163     y = R .* sin(phi) .* sin(theta);
0164     z = R .* cos(phi);
0165 
0166     V = [x y z];
0167     
0168 return

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