inter_area.m 2.71 KB
 Julius Welzel committed Mar 12, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 ``````%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % inter_area (fun) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This function calucate the area of intersection of a 3D topographie and % and plane in z-axis % Input Variables: - chanlocs: Standard EEGLab chanlocs format with x,y,z % coordinates of channel locations % - map: Vector [1 *length(chanlocs)] with value of % topographie % - lvl_inter: Vector of position of z-level of % intersection % % % Output Variables: - area_inter: Vector with area for every intersection % level [1 *length(lvl_inter)] % Author: Julius Welzel, Mareike Daeglau, University of Oldenburg % Contact: julius.welzel@gmail.com % Version: 1.0 // 11.03.19 // inital setup function area_inter = inter_area(chanlocs, map, n_inter, ID, PATHOUT,c_map) % Calculate and interpolate X & Y coordinates from chanlocs p_interp = 200; % points to interpolate (should be higher for less channels % get cartesian coordinates by MIKE X COHEN :D and grid data [elocsX,elocsY] = pol2cart(pi/180*[chanlocs.theta],[chanlocs.radius]); elocsX = elocsX*-1; %swap left/ right xlin = linspace(min(elocsX),max(elocsX),p_interp); ylin = linspace(min(elocsY),max(elocsY),p_interp); [Xgrid, Ygrid] = meshgrid(xlin,ylin); %% Get 3D representation of indv map min_z = min(min(map)); max_z = max(max(map)); lvl_inter = linspace(max_z,min_z,n_inter); Z = griddata(elocsX,elocsY,map,Xgrid,Ygrid,'cubic'); % estimate Z values for map_st area_inter = NaN(1,length(lvl_inter)); p = numSubplots(length(lvl_inter)); % for subplot config figure for i = 1:length(lvl_inter) try ax = subplot(p(1),p(2),i) mesh(Xgrid,Ygrid,Z) % plot map 3D colormap(ax,c_map) hold on surf(Xgrid,Ygrid,lvl_inter(i)*ones(size(Z)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none') hold on % Take the difference between the two surface heights and find the contour % where that surface is zero. zdiff = Z - (lvl_inter(i)*ones(size(Z))); C = contours(Xgrid, Ygrid, zdiff,[0,0]); %claculate the intersection area [a, ca] = cont_area(C,Xgrid,Ygrid,Z); % Visualize the line. cl = line(ca.xL, ca.yL, ca.zL, 'Color', 'k', 'LineWidth', 3); title(['Area : ' num2str(round(a,2)) ' at z = ' num2str(round(lvl_inter(i),2))]) area_inter(i) = a; catch break area_inter(i) = 0; end end sgtitle (['Subj ' num2str(ID)]) save_fig(gcf,PATHOUT,['SurfArea_subj_' num2str(ID)],'figtype','.png','figsize',[0 0 45 25]) end ``````