标签：... end plot 绘图 matlab ----- circle venn data
1、离散热力图
function out = scatplot(x,y,method,radius,N,n,po,ms)
% Scatter plot with color indicating data density
%
% USAGE:
% out = scatplot(x,y,method,radius,N,n,po,ms)
% out = scatplot(x,y,dd)
%
% DESCRIPTION:
% Draws a scatter plot with a colorscale
% representing the data density computed
% using three methods
%
% INPUT VARIABLES:
% x,y - are the data points
% method - is the method used to calculate data densities:
% 'circles' - uses circles with a determined area
% centered at each data point
% 'squares' - uses squares with a determined area
% centered at each data point
% 'voronoi' - uses voronoi cells to determin data densities
% default method is 'voronoi'
% radius - is the radius used for the circles or squares
% used to calculate the data densities if
% (Note: only used in methods 'circles' and 'squares'
% default radius is sqrt((range(x)/30)^2 + (range(y)/30)^2)
% N - is the size of the square mesh (N x N) used to
% filter and calculate contours
% default is 100
% n - is the number of coeficients used in the 2-D
% running mean filter
% default is 5
% (Note: if n is length(2), n(2) is tjhe number of
% of times the filter is applied)
% po - plot options:
% 0 - No plot
% 1 - plots only colored data points (filtered)
% 2 - plots colored data points and contours (filtered)
% 3 - plots only colored data points (unfiltered)
% 4 - plots colored data points and contours (unfiltered)
% default is 1
% ms - uses this marker size for filled circles
% default is 4
%
% OUTPUT VARIABLE:
% out - structure array that contains the following fields:
% dd - unfiltered data densities at (x,y)
% ddf - filtered data densities at (x,y)
% radius - area used in 'circles' and 'squares'
% methods to calculate densities
% xi - x coordenates for zi matrix
% yi - y coordenates for zi matrix
% zi - unfiltered data densities at (xi,yi)
% zif - filtered data densities at (xi,yi)
% [c,h] = contour matrix C as described in
% CONTOURC and a handle H to a contourgroup object
% hs = scatter points handles
%
%Copy-Left, Alejandro Sanchez-Barba, 2005
if nargin==0
scatplotdemo
return
end
if nargin<3 | isempty(method)
method = 'vo';
end
if isnumeric(method)
gsp(x,y,method,2)
return
else
method = method(1:2);
end
if nargin<4 | isempty(n)
n = 5; %number of filter coefficients
end
if nargin<5 | isempty(radius)
radius = sqrt((range(x)/30)^2 + (range(y)/30)^2);
end
if nargin<6 | isempty(po)
po = 1; %plot option
end
if nargin<7 | isempty(ms)
ms = 4; %markersize
end
if nargin(x(k)-r) & x(y(k)-r) & y<(y(k)+r) );
end %for
area = (2*r)^2;
dd = dd/area;
case 'ci'
for k=1:Ld
dd(k) = sum( sqrt((x-x(k)).^2 + (y-y(k)).^2) < r );
end
area = pi*r^2;
dd = dd/area;
case 'vo' %----- Using voronoi cells ------
[v,c] = voronoin([x,y]);
for k=1:length(c)
%If at least one of the indices is 1,
%then it is an open region, its area
%is infinity and the data density is 0
if all(c{k}>1)
a = polyarea(v(c{k},1),v(c{k},2));
dd(k) = 1/a;
end %if
end %for
end %switch
return
%~~~~~~~~~~ Graf Scatter Plot ~~~~~~~~~~~
function varargout = gsp(x,y,c,ms)
%Graphs scattered poits
map = colormap;
ind = fix((c-min(c))/(max(c)-min(c))*(size(map,1)-1))+1;
h = [];
%much more efficient than matlab's scatter plot
for k=1:size(map,1)
if any(ind==k)
h(end+1) = line('Xdata',x(ind==k),'Ydata',y(ind==k), ...
'LineStyle','none','Color',map(k,:), ...
'Marker','.','MarkerSize',ms);
end
end
if nargout==1
varargout{1} = h;
end
returnclear,clc; x = normrnd(10,1,1000,1); y = x * 3 + normrnd(10,1,1000,1); scatplot(double(x),double(y))
效果如下：
2、绘制3D球体
function scatter3sph(X,Y,Z,varargin)
%SCATTER3SPH (X,Y,Z) Makes a 3d scatter plot with 3D spheres
% SCATTER3SPH is like scatter3 only drawing spheres instead
% of flat circles, at coordinates specified by vectors X, Y, Z. All three
% vectors have to be of the same length.
% SCATTER3SPH(X,Y,Z) draws the spheres with the default size and color.
% SCATTER3SPH(X,Y,Z,'size',S) draws the spheres with sizes S. If length(S)= 1
% the same size is used for all spheres.
% SCATTER3SPH(X,Y,Z,'color',C) draws the spheres with colors speciffied in a
% N-by-3 matrix C as RGB values.
% SCATTER3SPH(X,Y,Z,'transp',T) applies transparency level 'T' to the spheres
% T= 0 => transparent, T= 1 => opaque.
% Parameter names can be abreviated to 3 letters. For example: 'siz' or
% 'col'. Case is irrelevant.
%
% Example
% %Coordinates
% X= 100*rand(9,1); Y= 100*rand(9,1); Z= 100*rand(9,1);
%
% %Colors: 3 blue, 3 red and 3 green
% C= ones(3,1)*[0 0 1];
% C= [C;ones(3,1)*[1 0 0]];
% C= [C;ones(3,1)*[0 1 0]];
%
% %Sizes
% S= 5+10*rand(9,1);
%
% scatter3sph(X,Y,Z,'size',S,'color',C,'trans',0.3);
% axis vis3d
%-- Some checking...
if nargin < 3 error('Need at least three arguments'); return; end
if mean([length(X),length(Y),length(Z)]) ~= length(X) error ('Imput vectors X, Y, Z are of different lengths'); return; end
%-- Defaults
C= ones(length(X),1)*[0 0 1];
S= 0.1*max([X;Y;Z])*ones(length(X),1);
nfacets= 15;
transp= 0.5;
%-- Extract optional arguments
for j= 1:2:length(varargin)
string= lower(varargin{j});
switch string(1:min(3,length(string)))
case 'siz'
S= varargin{j+1};
if length(S) == 1
S= ones(length(X),1)*S;
elseif length(S) < length(X)
error('The vector of sizes must be of the same length as coordinate vectors (or 1)');
return
end
case 'col'
C= varargin{j+1};
if size(C,2) < 3 error('Colors matrix must have 3 columns'); return; end
if size(C,1) == 1
C= ones(length(X),1)*C(1:3);
elseif size(C,1) < length(X)
error('Colors matrix must have the same number of rows as length of coordinate vectors (or 1)');
return
end
case 'fac'
nfacets= varargin{j+1};
case 'tra'
transp= varargin{j+1};
otherwise
error('Unknown parameter name. Allowed names: ''size'', ''color'', ''facets'', ''transparency'' ');
end
end
%-- Sphere facets
[sx,sy,sz]= sphere(nfacets);
%--- Correct potential distortion
maxax= max([range(X), range(Y), range(Z)]);
ratios= [range(X)/maxax, range(Y)/maxax, range(Z)/maxax];
sx= sx*ratios(1);
sy= sy*ratios(2);
sz= sz*ratios(3);
%-- Plot spheres
hold on
for j= 1:length(X)
surf(sx*S(j)+X(j), sy*S(j)+Y(j), sz*S(j)+Z(j),...
'LineStyle','none',...
'FaceColor',C(j,:),...
'FaceAlpha',transp);
end
daspect([ratios(1), ratios(2), ratios(3)]);
light('Position',[1 1 1],'Style','infinit','Color',[1 1 1]);
lighting gouraud;
view(30,30)clear,clc; X= 100*rand(9,1); Y= 100*rand(9,1); Z= 100*rand(9,1); C= ones(3,1)*[0 0 1]; C= [C;ones(3,1)*[1 0 0]]; C= [C;ones(3,1)*[0 1 0]]; S= 5+10*rand(20,1); scatter3sph(X,Y,Z,'size',S,'color',C,'trans',0.3); axis vis3d grid on
3、matlab绘制复杂矢量图
function hh = quiverwcolorbar(varargin)
% Melissa Day 5/24/2013
% Upgrade of Andrew Diamond's quiverc2wcmap to generate a quiver plot
% with arrows colored according to vector magnitude.
% Functional differences from quiverc2wcmap:
% 1) Allows user to specify colormap boundaries using 'bound': changes
% colorbar axes AND corresponding vector coloring
% (much more useful for intercomparison of datasets)
% 2) Improved fidelity to magnitude colors in large datasets
% (at a cost of increased computation time...)
% Uses some improvements from DS's vfield_color
% 3) Small clarifications added
% Example:
% x = rand(1,50).*100;
% y = rand(1,50).*100;
% u = rand(1,50) .* 10;
% v = rand(1,50) .* 10;
% scale = 0;
% figure; quiverwcolorbar(x',y',u',v',scale); %compare to:
% figure; quiverwcolorbar(x',y',u',v',scale,'bounds',[0 10]);
%----------------
% function hh = quiverc2wcmap(varargin)
% Andrew Diamond 3/17/2005
% This is based off of Tobias H鰂fken which was based off of Bertrand Dano
% keeping with their main intent to show a color w/r to magnitude quiver plot
% while maintaining a large amount of quiver API backwards compatability.
% Functional differences from quiverc2:
% 1) This works under 6.5.1
% 2) It handles NaNs
% 3) It draws a colormap that is w/r to the quiver magnitudes (hard coded to
% 20 segments of the colormap as per quiverc2 - seems fine for a quiver).
% 4) Various bug fixes (I think)
% In order to do this I needed some small hacks on 6.5.1's quiver to take a
% color triplet. I have included as part of this file in a subfunction below.
%----------------
% Comments from quiverc2
% changed Tobias H鰂fken 3-14-05
% totally downstripped version of the former
% split input field into n segments and do a quiver qhich each of them
% Modified version of Quiver to plots velocity vectors as arrows
% with components (u,v) at the points (x,y) using the current colormap
% Bertrand Dano 3-3-03
% Copyright 1984-2002 The MathWorks, Inc.
% changed T. H鰂fken 14.03.05, for high data resolution
% using fixed color "spacing" of 20
%QUIVERC Quiver color plot.
% QUIVERC(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y). The matrices X,Y,U,V must all be the same size
% and contain corresponding position and velocity components (X and Y
% can also be vectors to specify a uniform grid). QUIVER automatically
% scales the arrows to fit within the grid.
%
% QUIVERC(U,V) plots velocity vectors at equally spaced points in
% the x-y plane.
%
% QUIVERC(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
% arrows to fit within the grid and then stretches them by S. Use
% S=0 to plot the arrows without the automatic scaling.
%
% QUIVERC(...,LINESPEC) uses the plot linestyle specified for
% the velocity vectors. Any marker in LINESPEC is drawn at the base
% instead of an arrow on the tip. Use a marker of '.' to specify
% no marker at all. See PLOT for other possibilities.
%
% QUIVERC(...,'filled') fills any markers specified.
%
% H = QUIVERC(...) returns a vector of line handles.
%
% Example:
% [x,y] = meshgrid(-2:.2:2,-1:.15:1);
% z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
% contour(x,y,z), hold on
% quiverc(x,y,px,py), hold off, axis image
%
% See also FEATHER, QUIVER3, PLOT.
% Clay M. Thompson 3-3-94
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.21 $ $Date: 2002/06/05 20:05:16 $
%-------------------------------------------------------------
nin = nargin; %number of inputs
% error(nargchk(2,5,nin));
error(nargchk(2,7,nin)); %added +2 to maxargs to account for 'bounds' add
% Check numeric input arguments
if nin= 5) % quiver(x,y,u,v,s) or quiver(x,y,u,v,s,'bounds',[start end])
if(isscalar(varargin{5}))
scale = varargin{5};
end
end
%-------------Define matrix of vector magnitudes-------------
% Define matrix of vector magnitudes
vr = sqrt(u.^2+v.^2);
% From quiverc2wcmap:
% if data has Nan, don't let it contaminate the computations that segment the
% data; I could just do this with vr and then get clever with the indices but
% this make for easy debugging and as this is a graphics routine the computation
% time is completely irrelevant.
nonNaNind = find(~isnan(vr(:)));
xyuvvrNN = [x(nonNaNind),y(nonNaNind),u(nonNaNind),v(nonNaNind),vr(nonNaNind)];
[xyuvvrNNs, xyuvvrNNsi] = sortrows(xyuvvrNN,5);
% From quiverc2wcmap, no longer necessary
% n = 20; %number of colors
% CC = colormap;
% iCCs=round(linspace(1,size(CC,1),n));
% iData = round(linspace(0,size(xyuvvrNNs,1),n+1));
% figure;
%-------------Generate colorbar-------------
% Includes a clever way to generate (and subsequently hide) an image that
% is required for colorbar without running out of memory for large datasets.
%
% Condensed comments from quiverc2wcmap:
% In 6.5.1 if you ever want a colorbar tick marks to reflect real data ranges
% (versus just the indices of a colormap) then there apparently has to be an
% image somewhere in the figure. Of course, I don't want an image but I figured
% I just make it invisible and then draw the quiver plot over it.
% Unfortunately, it seems that colorbar uses caxis to retrive the data range in
% the image and for invisible images it always seems to be 0 UNLESS you
% explictly reset the caxis.
% This will work but then the axis will be oversized to accomodate the image
% because images have their first and last vitual "pixels" CENTERED around the
% implicit or explict xmin,xmax,ymin,ymax (as per imagesc documentation) but the
% start and end of each of these "pixels" is +/- half a unit where the unit
% corresponds to subdividing the limits by the number of pixels (-1). Given
% that formula and given my invisible 2x2 image for which it is desired to
% diplay in an axis that ISN'T oversized (i.e. min(x), max(x), min(y),max(y)) it
% is possible to solve for the limits (i.e. an artifically reduced limit)
% that need to be specified for imagesc to make its final oversized axis to be the
% non oversized axis that we really want.
% xa,xb,ya,yb compenstates for the axis extention given by imagesc to make
% pixels centered at the limit extents (etc.). Note, this "easy"
% formula would only work for 2x2 pixel images.
xs = min(x); %x(1);
xf = max(x); %x(end)
xa = (3 * xs + xf)/4;
xb = (3 * xf + xs)/4;
ys = min(y); %y(1);
yf = max(y); %y(end)
ya = (3 * ys + yf)/4;
yb = (3 * yf + ys)/4;
% Determine magnitude min/max (which is reflected in colorbar)
colormin = min(xyuvvrNNs(:,5)); %column 5 is NaN-cleared vr
colormax = max(xyuvvrNNs(:,5));
% Allow user to edit bounds using "'bounds',[colormin colormax]" input
for k=1:nin
if (k~=nin) && (length(varargin{k})==6) && strcmp(varargin{k},'bounds')
bounds = varargin{k+1};
if isempty(bounds), error('Specify colormap boundaries'); end
colormin = bounds(1);
colormax = bounds(2);
end
end
% mapbounds = reshape(xyuvvrNNs([1,end;1,end],5),2,2); %from quiverc2wcmap
mapbounds = [colormin colormax; colormin colormax];
h=imagesc([xa,xb],[ya,yb],mapbounds);
set(h,'visible','off');
% Prep colorbar
rang = (colormax-colormin)/colormax;
ticknum = 6; %if you want to toggle number of ticks on colorbar
incr = rang./(ticknum-1);
B = [colormin/colormax:incr:1];
B = B.*colormax;
C = sprintf(['%4.2e',repmat([' \n%4.2e'], 1, ticknum)],B);
C = str2num(C);
caxis([colormin colormax])
colorbar('EastOutside','ytick',B,'yticklabel',C,...
'ticklength',[0.04 0.1],'YLim',[B(1) B(ticknum)],'FontSize',20)
% In quiverc2wcmap this loop plotted for each color level (n=20) and was very
% fast, but I found it did not plot some large data sets with enough color
% accuracy. Switched the loop to examine each data point individually.
% Takes longer but I'm more confident that the colors are correct. Note:
% much of this overhaul was inspired by DS's vfield_color.
hold on;
cmap = jet(64); %toggle type of colormap
CC = colormap(cmap);
cm_stepsize = (colormax-colormin)/length(CC);
for it=1:size(xyuvvrNNs,1) %takes ~13 seconds for a 8730-point dataset
% Some quiverc2wcmap fragments for reference:
% for it=1:n %10x faster, but may not be accurate
% c = CC(iCCs(it),:); %colormap([1:64](it),:) %ie "This RGB color row ="
% si = iData(it)+1; %[1:size(data)](it)+1; %ie "This start row"
% ei = iData(it+1); %[1:size(data)](it+1); %ie "This end row"
% hh=quiver(xyuvvrNNs(si:ei,1),xyuvvrNNs(si:ei,2),...
% xyuvvrNNs(si:ei,3),xyuvvrNNs(si:ei,4),scale*it/n,'Color',c)
cm_index = floor( (xyuvvrNNs(it,5) - colormin) / ( cm_stepsize ) ) + 1;
if cm_index == 1 %in case colormin is zero
c = CC(cm_index,:);
elseif cm_index > length(CC) %in case max(xyuvvrNNs) > colormax
cm_index = length(CC);
c = CC(cm_index,:);
elseif cm_index <= 0 %in case min(xyuvvrNNs) < colormin
cm_index = 1;
c = CC(cm_index,:);
else
cm_index = cm_index-1;
c = CC(cm_index,:);
end
hh=quiver(xyuvvrNNs(it,1),xyuvvrNNs(it,2),...
xyuvvrNNs(it,3),xyuvvrNNs(it,4),scale,'Color',c);
end
%----------Rest of document is from quiverc2wcmap---------------
% This is Matlab's 6.5.1 quiver. I figure that ensures a fair amouint of backward
% compatibility but I needed to hack it to allow a 'Color' property. Obviously
% a person could do more.
function hh = quiver(varargin)
%QUIVER Quiver plot.
% QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y). The matrices X,Y,U,V must all be the same size
% and contain corresponding position and velocity components (X and Y
% can also be vectors to specify a uniform grid). QUIVER automatically
% scales the arrows to fit within the grid.
%
% QUIVER(U,V) plots velocity vectors at equally spaced points in
% the x-y plane.
%
% QUIVER(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
% arrows to fit within the grid and then stretches them by S. Use
% S=0 to plot the arrows without the automatic scaling.
%
% QUIVER(...,LINESPEC) uses the plot linestyle specified for
% the velocity vectors. Any marker in LINESPEC is drawn at the base
% instead of an arrow on the tip. Use a marker of '.' to specify
% no marker at all. See PLOT for other possibilities.
%
% QUIVER(...,'filled') fills any markers specified.
%
% H = QUIVER(...) returns a vector of line handles.
%
% Example:
% [x,y] = meshgrid(-2:.2:2,-1:.15:1);
% z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
% contour(x,y,z), hold on
% quiver(x,y,px,py), hold off, axis image
%
% See also FEATHER, QUIVER3, PLOT.
% Clay M. Thompson 3-3-94
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.21 $ $Date: 2002/06/05 20:05:16 $
% Arrow head parameters
alpha = 0.33; % Size of arrow head relative to the length of the vector
beta = 0.33; % Width of the base of the arrow head relative to the length
autoscale = 1; % Autoscale if ~= 0 then scale by this.
plotarrows = 1; % Plot arrows
sym = '';
filled = 0;
ls = '-';
ms = '';
col = 'b';
nin = nargin;
ColorSpecInd = find(strcmpi(varargin, 'color'));
if(length(ColorSpecInd)==1 & nin > ColorSpecInd)
col = varargin{ColorSpecInd+1};
varargin = varargin([1:ColorSpecInd-1,ColorSpecInd+2:nin]);
nin = nin-2;
end
% Parse the string inputs
while isstr(varargin{nin}),
vv = varargin{nin};
if ~isempty(vv) & strcmp(lower(vv(1)),'f')
filled = 1;
nin = nin-1;
else
[l,c,m,msg] = colstyle(vv);
if ~isempty(msg),
error(sprintf('Unknown option "%s".',vv));
end
if ~isempty(l), ls = l; end
if ~isempty(c), col = c; end
if ~isempty(m), ms = m; plotarrows = 0; end
if isequal(m,'.'), ms = ''; end % Don't plot '.'
nin = nin-1;
end
end
error(nargchk(2,5,nin));
% Check numeric input arguments
if nin0
len = sqrt((u.^2 + v.^2)/del);
maxlen = max(len(:));
else
maxlen = 0;
end
if maxlen>0
autoscale = autoscale*0.9 / maxlen;
else
autoscale = autoscale*0.9;
end
u = u*autoscale; v = v*autoscale;
end
ax = newplot;
next = lower(get(ax,'NextPlot'));
hold_state = ishold;
% Make velocity vectors
x = x(:).'; y = y(:).';
u = u(:).'; v = v(:).';
uu = [x;x+u;repmat(NaN,size(u))];
vv = [y;y+v;repmat(NaN,size(v))];
% h1 = plot(uu(:),vv(:),[col ls]);
h1 = plot(uu(:),vv(:),ls,'Color',col);
if plotarrows,
% Make arrow heads and plot them
hu = [x+u-alpha*(u+beta*(v+eps));x+u; ...
x+u-alpha*(u-beta*(v+eps));repmat(NaN,size(u))];
hv = [y+v-alpha*(v-beta*(u+eps));y+v; ...
y+v-alpha*(v+beta*(u+eps));repmat(NaN,size(v))];
hold on
% h2 = plot(hu(:),hv(:),[col ls]);
h2 = plot(hu(:),hv(:),ls,'Color',col);
else
h2 = [];
end
if ~isempty(ms), % Plot marker on base
hu = x; hv = y;
hold on
% h3 = plot(hu(:),hv(:),[col ms]);
h3 = plot(hu(:),hv(:),ls,'Color',col);
if filled, set(h3,'markerfacecolor',get(h1,'color')); end
else
h3 = [];
end
if ~hold_state, hold off, view(2); set(ax,'NextPlot',next); end
if nargout>0, hh = [h1;h2;h3]; end
function retval = isscalar(m)
retval = prod(size(m)) == 1;
clear，clc; x = rand(1,100).*100; y = rand(1,100).*100; u = rand(1,100) .* 10; v = rand(1,100) .* 10; scale = 0; figure; quiverwcolorbar(x',y',u',v',scale); figure; quiverwcolorbar(x',y',u',v',scale,'bounds',[0 10]);
4、精细化作图工具箱
因为是数据文件只能加群：【matlab&&python资源共享】：获取panel-2.14.rar文件
4、matlab绘制维恩图
function varargout = venn (varargin)
%VENN Plot 2- or 3- circle area-proportional Venn diagram
%
% venn(A, I)
% venn(Z)
% venn(..., F)
% venn(..., 'ErrMinMode', MODE)
% H = venn(...)
% [H, S] = venn(...)
% [H, S] = venn(..., 'Plot', 'off')
% S = venn(..., 'Plot', 'off')
% [...] = venn(..., P1, V1, P2, V2, ...)
%
%venn(A, I) by itself plots circles with total areas A, and intersection
%area(s) I. For two-circle venn diagrams, A is a two element vector of circle
%areas [c1 c2] and I is a scalar specifying the area of intersection between
%them. For three-circle venn diagrams, A is a three element vector [c1 c2 c3],
%and I is a four element vector [i12 i13 i23 i123], specifiying the
%two-circle intersection areas i12, i13, i23, and the three-circle
%intersection i123.
%
%venn(Z) plots a Venn diagram with zone areas specified by the vector Z.
%For a 2-circle venn diagram, Z is a three element vector [z1 z2 z12]
%For a 3-circle venn, Z is a 7 element vector [z1 z2 z3 z12 z13 z23 z123]
%
%venn(..., F) specifies optional optimization options. VENN uses FMINBND to
%locate optimum pair-wise circle distances, and FMINSEARCH to optimize
%overall three-circle alignment. F is a structure with fields specifying
%optimization options for these functions. F may be a two-element array of
%structures, in which case the first structure is used for FMINBND
%function calls, and the second structure is used for FMINSEARCH function
%calls.
%
%venn(..., 'ErrMinMode', MODE)
%Used for 3-circle venn diagrams only. MODE can be 'TotalError' (default),
%'None', or 'ChowRodgers'. When ErrMinMode is 'None', the positions and
%sizes of the three circles are fixed by their pairwise-intersections,
%which means there may be a large amount of error in the area of the three-
%circle intersection. Specifying ErrMinMode as 'TotalError' attempts to
%minimize the total error in all four intersection zones. The area of the
%three circles are kept constant in proportion to their populations. The
%'ChowRodgers' mode uses the the method proposed by Chow and Rodgers
%[Ref. 1] to draw 'nice' three-circle venn diagrams which appear more
%visually representative of the desired areas, although the actual areas of
%the circles are allowed to deviate from requested values.
%
%H = venn(...) returns a two- or three- element vector to the patches
%representing the circles.
%
%[H, S] = venn(...) returns a structure containing descriptive values
%computed for the requested venn diagram. S is a structure with the
%following fields, where C is the number of circles (N = 2 or 3), Z is
%the number of zones (Z = 3 or 7), and I is the number of intersection
%areas (1 or 4)
%
% Radius C-element vector of circle radii
%
% Position C*2 array of circle centers
%
% ZoneCentroid Z*2 array of zone centroids (Can be used for labeling)
%
% CirclePop C-element vector of supplied circle populations.
% (I.e., the 'true' circle areas)
%
% CircleArea C-element of actual circle areas
%
% CircleAreaError = (CircleArea-CirclePop)/CirclePop
%
% IntersectPop I-element vector of supplied intersection populations
% (I.e., the 'true' intersection areas)
%
% IntersectArea I-element vector of actual intersection areas
%
% IntersectError = (IntersectArea-IntersectPop)/IntersectPop
%
% ZonePop Z-element vector of supplied zone populations. (I.e.
% 'true' zone areas
%
% ZoneArea Z-element vector of actual zone areas.
%
% ZoneAreaError = (ZoneArea-ZonePop)/ZonePop
%
%
%[H, S] = venn(..., 'Plot', 'off')
%S = venn(..., 'Plot', 'off')
%Returns a structure of computed values, without plotting the diagram. This
%which can be useful when S is used to draw custom venn diagrams or for
%exporting venn diagram data to another application. When Plot is set to off,
%the handles vector H is returned as an empty array. Alternatively, the command
%S = venn(..., 'Plot', 'off) will return only the output structure.
%
%[...] = venn(..., P1, V1, P2, V2, ...)
%Specifies additional patch settings in standard Matlab parameter/value
%pair syntax. Parameters can be any valid patch parameter. Values for patch
%parameters can either be single values, or a cell array of length LENGTH(A),
%in which case each value in the cell array is applied to the corresponding
%circle in A.
%
%Examples
%
% %Plot a simple 2-circle venn diagram with custom patch properties
% figure, axis equal, axis off
% A = [300 200]; I = 150;
% venn(A,I,'FaceColor',{'r','y'},'FaceAlpha',{1,0.6},'EdgeColor','black')
%
% %Compare ErrMinModes
% A = [350 300 275]; I = [100 80 60 40];
% figure
% subplot(1,3,1), h1 = venn(A,I,'ErrMinMode','None');
% axis image, title ('No 3-Circle Error Minimization')
% subplot(1,3,2), h2 = venn(A,I,'ErrMinMode','TotalError');
% axis image, title ('Total Error Mode')
% subplot(1,3,3), h3 = venn(A,I,'ErrMinMode','ChowRodgers');
% axis image, title ('Chow-Rodgers Mode')
% set([h1 h2], 'FaceAlpha', 0.6)
%
% %Using the same areas as above, display the error optimization at each
% iteration. Get the output structure.
% F = struct('Display', 'iter');
% [H,S] = venn(A,I,F,'ErrMinMode','ChowRodgers','FaceAlpha', 0.6);
%
% %Now label each zone
% for i = 1:7
% text(S.ZoneCentroid(i,1), S.ZoneCentroid(i,2), ['Zone ' num2str(i)])
% end
%
%See also patch, bar, optimset, fminbdn, fminsearch
%
%Copyright (C) 2008 Darik Gamble, University of Waterloo.
%dgamble@engmail.uwaterloo.ca
%
%References
%1. S Chow and P Rodgers. Extended Abstract: Constructing Area-Proportional
% Venn and Euler Diagrams with Three Circles. Presented at Euler Diagrams
% Workshop 2005. Paris. Available online:
% http://www.cs.kent.ac.uk/pubs/2005/2354/content.pdf
%
%2. S Chow and F Ruskey. Drawing Area-Proportional Venn and Euler Diagrams.
% Lecture Notes in Computer Science. 2004. 2912: 466-477. Springer-Verlag.
% Available online: http://www.springerlink.com/content/rxhtlmqav45gc84q/
%
%3. MP Fewell. Area of Common Overlap of Three Circles. Australian Government
% Department of Defence. Defence Technology and Science Organisation. 2006.
% DSTO-TN-0722. Available online:
% http://dspace.dsto.defence.gov.au/dspace/bitstream/1947/4551/4/DSTO-TN-0722.PR.pdf
%Variable overview
% A0, A Desired and actual circle areas
% A = [A1 A2] or [A1 A2 A3]
% I0, I Desired and actual intersection areas
% I = I12 or [I12 I13 I23 I123]
% Z0, Z Desired and actual zone areas
% Z = [Z1 Z2 Z12] or [Z1 Z2 Z3 Z12 Z13 Z23 Z123]
% x, y Circle centers
% x = [x1 x2] or [x1 x2 x3]
% r Circle radii
% r = [r1 r2] or [r1 r2 r3]
% d Pair-wise distances between circles
% d = d12 or [d12 d13 d23]
%Parse input arguments and preallocate settings
[A0, I0, Z0, nCirc, fminOpts, vennOpts, patchOpts] = parseArgsIn (varargin);
[d, x, y, A, I, Z] = preallocVectors (nCirc);
zoneCentroids = []; %Will only be calculated if needed
%Circle Radii
r = sqrt(A0/pi);
%Determine distance between first circle pair
d(1) = circPairDist(r(1), r(2), I0(1), fminOpts(1));
%Position of second circle is now known
x(2) = d(1);
%First intersection area
I(1) = areaIntersect2Circ(r(1), r(2), d(1));
if nCirc==3
%Pairwise distances for remaining pairs 1&3 and 2&3
d(2) = circPairDist(r(1), r(3), I0(2), fminOpts(1)); %d13
d(3) = circPairDist(r(2), r(3), I0(3), fminOpts(1)); %d23
%Check triangle inequality
srtD = sort(d);
if ~(srtD(end)<(srtD(1)+srtD(2)))
error('venn:triangleInequality', 'Triangle inequality not satisfied')
end
%Guess the initial position of the third circle using the law of cosines
alpha = acos( (d(1)^2 + d(2)^2 - d(3)^2) / (2 * d(1) * d(2)) );
x(3) = d(2)*cos(alpha);
y(3) = d(2)*sin(alpha);
%Each pair-wise intersection fixes the distance between each pair
%of circles, so technically there are no degrees of freedom left in
%which to adjust the three-circle intersection. We can either try
%moving the third circle around to minimize the total error, or
%apply Chow-Rodgers
switch vennOpts.ErrMinMode
case 'TotalError'
%Minimize total intersection area error by moving the third circle
pos = fminsearch(@threeCircleAreaError, [x(3) y(3)], fminOpts(2));
x(3) = pos(1);
y(3) = pos(2);
case 'ChowRodgers'
%note that doChowRodgersSearch updates x and y in this
%workspace as a nested fcn
doChowRodgersSearch;
end
%Make sure everything is 'up to date' after optimization
update3CircleData;
end
%Are we supposed to plot?
if vennOpts.Plot
if isempty(vennOpts.Parent)
vennOpts.Parent = gca;
end
hVenn = drawCircles(vennOpts.Parent, x, y, r, patchOpts.Parameters, patchOpts.Values);
else
hVenn = [];
end
%Only determine zone centroids if they're needed
%Needed for output structure
nOut = nargout;
if (nOut==1 && ~vennOpts.Plot) || nOut==2
if nCirc == 2
%Need to calculate new areas
A = A0; %Areas never change for 2-circle venn
Z = calcZoneAreas(2, A, I);
zoneCentroids = zoneCentroids2(d, r, Z);
else
zoneCentroids = zoneCentroids3(x, y, d, r, Z);
end
end
%Figure out output arguments
if nOut==1
if vennOpts.Plot
varargout{1} = hVenn;
else
varargout{1} = getOutputStruct;
end
elseif nOut==2
varargout{1} = hVenn;
varargout{2} = getOutputStruct;
end
function err = threeCircleAreaError (pos)
x3 = pos(1);
y3 = pos(2);
%Calculate distances
d(2) = sqrt(x3^2 + y3^2); %d13
d(3) = sqrt((x3-d(1))^2 + y3^2); %d23
%Calculate intersections
%Note: we're only moving the third circle, so I12 is not changing
I(2:3) = areaIntersect2Circ (r(1:2), r([3 3]), d(2:3)); %I13 and I23
I(4) = areaIntersect3Circ (r, d); %I123
%Replace 0 (no intersection) with infinite error
I(I==0) = Inf;
%Error
err = sum(abs((I-I0)./I0));
end
function doChowRodgersSearch
%Adapted from Ref. [1]
%Initialize an index matrix to select all 7choose2 zone pairs (21 pairs)
idx = nchoosek(1:7, 2);
%Which zone-zone pairs are considered equal?
%Zones within 10% of each other considered equal
zonePairAreas0 = Z0(idx);
%Percent difference in population between the two members of a pair
ar0 = 2*abs(zonePairAreas0(:,1)-zonePairAreas0(:,2))./sum(zonePairAreas0, 2)*100;
eqPairCutoff = 10;
pairIsEq = ar0<=eqPairCutoff;
%Calculate allowable range for pairs of zones considered unequal
if any(~pairIsEq)
%Sort zone areas
[zUneqAreas0, zUneqArea***tIdx] = sort(zonePairAreas0(~pairIsEq,:), 2);
%Make a real index array out of the inconvenient index sort returns
n = sum(~pairIsEq);
zUneqArea***tIdx = sub2ind([n,2], [1:n; 1:n]', zUneqArea***tIdx);
%rp = (largepopulation/smallpopulation)-1
rp = zUneqAreas0(:,2)./zUneqAreas0(:,1)-1;
rpMin = 1 + 0.3*rp;
rpMax = 1 + 2*rp;
end
%Preallocate zone error vector
zoneErr = zeros(1,21);
%Initialize independent parameters to search over
guessParams = [r(1) x(2) r(2) x(3) y(3) r(3)];
%Search!
pp = fminsearch(@chowRodgersErr, guessParams, fminOpts(2));
[r(1) x(2) r(2) x(3) y(3) r(3)] = deal(pp(1), pp(2), pp(3), pp(4), pp(5), pp(6));
function err = chowRodgersErr (p)
%params = [x2 r2 x3 y3 r3]
[r(1), x(2), r(2), x(3), y(3), r(3)] = deal(p(1), p(2), p(3), p(4), p(5), p(6));
%After changing x2, r2, x3, y3, and r3, update circle areas,
%distances, intersection areas, zone areas
update3CircleData;
if any(pairIsEq)
%For zone pairs considered equal, error is equal to square of the
%distance beyond the cutoff; 0 within cutoff
zAreas = Z(idx(pairIsEq,:));
ar = 2*abs(zAreas(:,1)-zAreas(:,2))./sum(zAreas, 2)*100;
isWithinRange = ar<eqPairCutoff;
ar(isWithinRange) = 0;
ar(~isWithinRange) = ar(~isWithinRange) - eqPairCutoff;
%Amplify error for equal zones with unequal areas
eqZoneUneqAreaErrorGain = 10;
ar(~isWithinRange) = ar(~isWithinRange)*eqZoneUneqAreaErrorGain;
zoneErr(pairIsEq) = ar.^2;
end
if any(~pairIsEq)
%For zone pairs considered unequal, error is equal to square of
%the distance from the allowable range of rp
%rp = (largepopulation/smallpopulation)-1
zUneqPairAreas = Z(idx(~pairIsEq,:));
%Sort based on the population sizes (determined by parent
%function doChowRodgersSearch)
zUneqPairAreas = zUneqPairAreas(zUneqArea***tIdx);
rp = zUneqPairAreas(:,2)./zUneqPairAreas(:,1)-1;
lessThanMin = rprpMax;
rp(~lessThanMin & ~moreThanMax) = 0;
%Determine how far out of range errors are
rp(lessThanMin) = rp(lessThanMin) - rpMin(lessThanMin);
rp(moreThanMax) = rp(moreThanMax) - rpMax(moreThanMax);
%Consider the case where rp < rpMin to be more
%erroneous than the case where rp > rpMax
tooSmallErrorGain = 10;
rp(lessThanMin) = rp(lessThanMin)*tooSmallErrorGain;
zoneErr(~pairIsEq) = rp.^2;
end
%Total error
err = sum(zoneErr);
end %chowRodgersErr
end %doChowRodgersSearch
function update3CircleData
%Circle areas
A = pi*r.^2;
%Calculate distances
d(1) = abs(x(2)); %d12
d(2) = sqrt(x(3)^2 + y(3)^2); %d13
d(3) = sqrt((x(3)-d(1))^2 + y(3)^2); %d23
%Calculate actual intersection areas
I(1:3) = areaIntersect2Circ (r([1 1 2]), r([2 3 3]), d); %I12, I13, I23
I(4) = areaIntersect3Circ (r, d); %I123
%Calculate actual zone areas
Z = calcZoneAreas(3, A, I);
end
function S = getOutputStruct
S = struct(...
'Radius' ,r ,...
'Position' ,[x' y'] ,...
'ZoneCentroid' ,zoneCentroids ,...
'CirclePop' ,A0 ,...
'CircleArea' ,A ,...
'CircleAreaError' ,(A-A0)./A0 ,...
'IntersectPop' ,I0 ,...
'IntersectArea' ,I ,...
'IntersectError' ,(I-I0)./I0 ,...
'ZonePop' ,Z0 ,...
'ZoneArea' ,Z ,...
'ZoneAreaError' ,(Z-Z0)./Z0 );
end
end %venn
function D = circPairDist (rA, rB, I, opts)
%Returns an estimate of the distance between two circles with radii rA and
%rB with area of intersection I
%opts is a structure of FMINBND search options
D = fminbnd(@areadiff, 0, rA+rB, opts);
function dA = areadiff (d)
intersectArea = areaIntersect2Circ (rA, rB, d);
dA = abs(I-intersectArea)/I;
end
end
function hCirc = drawCircles(hParent, xc, yc, r, P, V)
hAx = ancestor(hParent, 'axes');
nextplot = get(hAx, 'NextPlot');
%P and V are cell arrays of patch parameter/values
xc = xc(:); yc = yc(:); %Circle centers
r = r(:); %Radii
n = length(r);
%Independent parameter
dt = 0.05;
t = 0:dt:2*pi;
%Origin centered circle coordinates
X = r*cos(t);
Y = r*sin(t);
hCirc = zeros(1,n);
c = {'r', 'g', 'b'}; %default colors
fa = {0.6, 0.6, 0.6}; %default face alpha
tag = {'Circle1', 'Circle2', 'Circle3'}; %default tag
for i = 1:n
xx = X(i,:)+xc(i);
yy = Y(i,:)+yc(i);
hCirc(i) = patch (xx, yy, c{i}, 'FaceAlpha', fa{i}, 'Parent', hParent, 'Tag', tag{i});
if i==1
set(hAx, 'NextPlot', 'add');
end
end
set(hAx, 'NextPlot', nextplot);
%Custom patch parameter values
if ~isempty(P)
c = cellfun(@iscell, V);
%Scalar parameter values -- apply to all circles
if any(~c)
set(hCirc, {P{~c}}, {V{~c}});
end
%Parameters values with one value per circle
if any(c)
%Make sure all vals are column cell arrays
V = cellfun(@(val) (val(:)), V(c), 'UniformOutput', false);
set(hCirc, {P{c}}, [V{:}])
end
end
end %plotCircles
function A = areaIntersect2Circ (r1, r2, d)
%Area of Intersection of 2 Circles
%Taken from [2]
alpha = 2*acos( (d.^2 + r1.^2 - r2.^2)./(2*r1.*d) );
beta = 2*acos( (d.^2 + r2.^2 - r1.^2)./(2*r2.*d) );
A = 0.5 * r1.^2 .* (alpha - sin(alpha)) ...
+ 0.5 * r2.^2 .* (beta - sin(beta));
end
function [A, x, y, c, trngArea] = areaIntersect3Circ (r, d)
%Area of common intersection of three circles
%This algorithm is taken from [3].
% Symbol Meaning
% T theta
% p prime
% pp double prime
%[r1 r2 r3] = deal(r(1), r(2), r(3));
%[d12 d13 d23] = deal(d(1), d(2), d(3));
%Intersection points
[x,y,sinTp,cosTp] = intersect3C (r,d);
if any(isnan(x)) || any(isnan(y))
A = 0;
%No three circle intersection
return
end
%Step 6. Use the coordinates of the intersection points to calculate the chord lengths c1,
%c2, c3:
i1 = [1 1 2];
i2 = [2 3 3];
c = sqrt((x(i1)-x(i2)).^2 + (y(i1)-y(i2)).^2)';
%Step 7: Check whether more than half of circle 3 is included in the circular triangle, so
%as to choose the correct expression for the area
lhs = d(2) * sinTp;
rhs = y(2) + (y(3) - y(2))/(x(3) - x(2))*(d(2)*cosTp - x(2));
if lhs < rhs
sign = [-1 -1 1];
else
sign = [-1 -1 -1];
end
%Calculate the area of the three circular segments.
ca = r.^2.*asin(c/2./r) + sign.*c/4.*sqrt(4*r.^2 - c.^2);
trngArea = 1/4 * sqrt( (c(1)+c(2)+c(3))*(c(2)+c(3)-c(1))*(c(1)+c(3)-c(2))*(c(1)+c(2)-c(3)) );
A = trngArea + sum(ca);
end
function [x, y, sinTp, cosTp] = intersect3C (r, d)
%Calculate the points of intersection of three circles
%Adapted from Ref. [3]
%d = [d12 d13 d23]
%x = [x12; x13; x23]
%y = [y12; y13; y23]
% Symbol Meaning
% T theta
% p prime
% pp double prime
x = zeros(3,1);
y = zeros(3,1);
%Step 1. Check whether circles 1 and 2 intersect by testing d(1)
if ~( ((r(1)-r(2))<d(1)) && (d(1)<(r(1)+r(2))) )
%x = NaN; y = NaN;
%bigfix: no returned values for sinTp, cosTp
[x, y, sinTp, cosTp] = deal(NaN);
return
end
%Step 2. Calculate the coordinates of the relevant intersection point of circles 1 and 2:
x(1) = (r(1)^2 - r(2)^2 + d(1)^2)/(2*d(1));
y(1) = 0.5/d(1) * sqrt( 2*d(1)^2*(r(1)^2 + r(2)^2) - (r(1)^2 - r(2)^2)^2 - d(1)^4 );
%Step 3. Calculate the values of the sines and cosines of the angles tp and tpp:
cosTp = (d(1)^2 + d(2)^2 - d(3)^2) / (2 * d(1) * d(2));
cosTpp = -(d(1)^2 + d(3)^2 - d(2)^2) / (2 * d(1) * d(3));
sinTp = (sqrt(1 - cosTp^2));
sinTpp = (sqrt(1 - cosTpp^2));
%Step 4. Check that circle 3 is placed so as to form a circular triangle.
cond1 = (x(1) - d(2)*cosTp)^2 + (y(1) - d(2)*sinTp)^2 < r(3)^2;
cond2 = (x(1) - d(2)*cosTp)^2 + (y(1) + d(2)*sinTp)^2 >r(3)^2;
if ~(cond1 && cond2)
x = NaN; y = NaN;
return
end
%Step 5: Calculate the values of the coordinates of the relevant intersection points involving
%circle 3
xp13 = (r(1)^2 - r(3)^2 + d(2)^2) / (2 * d(2));
%yp13 = -0.5 / d(2) * sqrt( 2 * d(2)^2 * (r(2)^2 + r(3)^2) - (r(1)^2 - r(3)^2)^2 - d(2)^4 );
yp13 = -0.5 / d(2) * sqrt( 2 * d(2)^2 * (r(1)^2 + r(3)^2) - (r(1)^2 - r(3)^2)^2 - d(2)^4 );
x(2) = xp13*cosTp - yp13*sinTp;
y(2) = xp13*sinTp + yp13*cosTp;
xpp23 = (r(2)^2 - r(3)^2 + d(3)^2) / (2 * d(3));
ypp23 = 0.5 / d(3) * sqrt( 2 * d(3)^2 * (r(2)^2 + r(3)^2) - (r(2)^2 - r(3)^2)^2 - d(3)^4 );
x(3) = xpp23*cosTpp - ypp23*sinTpp + d(1);
y(3) = xpp23*sinTpp + ypp23*cosTpp;
end
function z = calcZoneAreas(nCircles, a, i)
%Uses simple set addition and subtraction to calculate the zone areas
%with circle areas a and intersection areas i
if nCircles==2
%a = [A1 A2]
%i = I12
%z = [A1-I12, A2-I12, I12]
z = [a(1)-i, a(2)-i, i];
elseif nCircles==3
%a = [A1 A2 A3]
%i = [I12 I13 I23 I123]
%z = [A1-I12-I13+I123, A2-I12-I23+I123, A3-I13-I23+I123, ...
% I12-I123, I13-I123, I23-I123, I123];
z = [a(1)-i(1)-i(2)+i(4), a(2)-i(1)-i(3)+i(4), a(3)-i(2)-i(3)+i(4), ...
i(1)-i(4), i(2)-i(4), i(3)-i(4), i(4)];
else
error('')
%This error gets caught earlier in the stack w. better error msgs
end
end
function [Cx, Cy, aiz] = centroid2CI (x, y, r)
%Finds the centroid of the area of intersection of two circles.
%Vectorized to find centroids for multiple circle pairs
%x, y, and r are nCirclePairs*2 arrays
%Cx and Cy are nCirclePairs*1 vectors
%Centroid of the area of intersection of two circles
n = size(x,1);
xic = zeros(n,2);
az = zeros(n,2);
dx = x(:,2)-x(:,1);
dy = y(:,2)-y(:,1);
d = sqrt(dx.^2 + dy.^2);
%Translate the circles so the first is at (0,0) and the second is at (0,d)
%By symmetry, all centroids are located on the x-axis.
%The two circles intersect at (xp, yp) and (xp, -yp)
xp = 0.5*(r(:,1).^2 - r(:,2).^2 + d.^2)./d;
%Split the inner zone in two
%Right side (Area enclosed by circle 1 and the line (xp,yp) (xp,-yp)
%Angle (xp,yp) (X1,Y1) (xp,-yp)
alpha = 2*acos(xp./r(:,1));
%Area and centroid of the right side of the inner zone
[xic(:,1) az(:,1)] = circleChordVals (r(:,1), alpha);
%Angle (xp,yp) (X2,Y2) (xp,-yp)
alpha = 2*acos((d-xp)./r(:,2));
%Area and centroid of the left side of the inner zone
[xic(:,2) az(:,2)] = circleChordVals (r(:,2), alpha);
xic(:,2) = d - xic(:,2);
%Thus the overall centroid & area of the inner zone
aiz = sum(az,2);
Cx = sum(az.*xic,2)./aiz;
%Now translate the centroid back based on the original positions of the
%circles
theta = atan2(dy, dx);
Cy = Cx.*sin(theta) + y(:,1);
Cx = Cx.*cos(theta) + x(:,1);
end
function centroidPos = zoneCentroids2 (d, r, Z)
centroidPos = zeros(3,2);
%Find the centroids of the three zones in a 2-circle venn diagram
%By symmetry, all centroids are located on the x-axis.
%First, find the x-location of the middle (intersection) zone centroid
%Centroid of the inner zone
centroidPos(3,1) = centroid2CI([0 d], [0 0], r);
%Now, the centroid of the left-most zone is equal to the centroid of
%the first circle (0,0) minus the centroid of the inner zone
centroidPos(1,1) = -centroidPos(3,1)*Z(3)/Z(1);
%Similarly for the right-most zone; the second circle has centroid at x=d
centroidPos(2,1) = (d*(Z(2)+Z(3)) - centroidPos(3,1)*Z(3))/Z(2);
end
function centroidPos = zoneCentroids3 (x0, y0, d, r, Z)
Z = Z(:);
%Get area, points of intersection, and chord lengths
[act, xi, yi, c, atr] = areaIntersect3Circ (r, d);
atr = atr(:);
r = r(:);
%Area and centroid of the triangle within the circular triangle is
xtr = sum(xi/3);
ytr = sum(yi/3);
%Now find the centroids of the three segments surrounding the triangle
i = [1 2; 1 3; 2 3];
xi = xi(i); yi = yi(i);
[xcs, ycs, acs] = circSegProps (r(:), x0(:), y0(:), xi, yi, c(:));
%Overall centroid of the circular triangle
xct = (xtr*atr + sum(xcs.*acs))/act;
yct = (ytr*atr + sum(ycs.*acs))/act;
%Now calculate the centroids of the three two-pair intersection zones
%(Zones 12 13 23)
%Entire zone centroid/areas
%x, y, and r are nCirclePairs*2 arrays
%Cx and Cy are nCirclePairs*1 vectors
i = [1 2; 1 3; 2 3];
[x2c, y2c, a2c] = centroid2CI (x0(i), y0(i), r(i));
%Minus the three-circle intersection zone
xZI2C = (x2c.*a2c - xct*act)./(a2c-act);
yZI2C = (y2c.*a2c - yct*act)./(a2c-act);
x0 = x0(:);
y0 = y0(:);
%Finally, the centroids of the three circles minus the intersection
%areas
i1 = [4 4 5]; i2 = [5 6 6];
j1 = [1 1 2]; j2 = [2 3 3];
x1C = (x0*pi.*r.^2 - xZI2C(j1).*Z(i1) - xZI2C(j2).*Z(i2) - xct*act)./Z(1:3);
y1C = (y0*pi.*r.^2 - yZI2C(j1).*Z(i1) - yZI2C(j2).*Z(i2) - yct*act)./Z(1:3);
%Combine and return
centroidPos = [x1C y1C; xZI2C yZI2C; xct yct];
end
function [x, a] = circleChordVals (r, alpha)
%For a circle centered at (0,0), with angle alpha from the x-axis to the
%intersection of the circle to a vertical chord, find the x-centroid and
%area of the region enclosed between the chord and the edge of the circle
%adapted from http://mathworld.wolfram.com/CircularSegment.html
a = r.^2/2.*(alpha-sin(alpha)); %Area
x = 4.*r/3 .* sin(alpha/2).^3 ./ (alpha-sin(alpha)); %Centroid
end
function [xc, yc, area] = circSegProps (r, x0, y0, x, y, c)
%Translate circle to (0,0)
x = x-[x0 x0];
y = y-[y0 y0];
%Angle subtended by chord
alpha = 2*asin(0.5*c./r);
%adapted from http://mathworld.wolfram.com/CircularSegment.html
area = r.^2/2.*(alpha-sin(alpha)); %Area
d = 4.*r/3 .* sin(alpha/2).^3 ./ (alpha-sin(alpha)); %Centroid
%Perpindicular bisector of the chord
m = -(x(:,2)-x(:,1))./(y(:,2)-y(:,1));
%angle of bisector
theta = atan(m);
%centroids
xc = d.*cos(theta);
yc = d.*sin(theta);
%Make sure we're on the correct side
%Point of intersection of the perp. bisector and the circle perimeter
xb = (x(:,1)+x(:,2))/2;
xc(xb<0) = xc(xb<0)*-1;
yc(xb<0) = yc(xb0
if isstruct(args{1})
%FMIN search options
f = args{1};
nIn = nIn - 1;
if nIn>0, args = args(2:end); end
if length(f) == 1
%Just double up
fminOpts = [f f];
elseif length(f) == 2
%ok
fminOpts = f;
else
error('venn:parseArgsIn', 'FMINOPTS must be a 1 or 2 element structure array.')
end
else
%Use defaults
fminOpts = [optimset('fminbnd'), optimset('fminsearch')];
end
else
%Use defaults
fminOpts = [optimset('fminbnd'), optimset('fminsearch')];
end
%If there's an even number of args in remaining
if nIn>0
if mod(nIn, 2)==0
%Parameter/Value pairs
p = args(1:2:end);
v = args(2:2:end);
[vennOpts, patchOpts] = parsePVPairs (p, v, nZones);
else
error('venn:parseArgsIn', 'Parameter/Value options must come in pairs')
end
else
vennOpts = defaultVennOptions;
patchOpts = struct('Parameters', [], 'Values', []);
end
end %parseArgsIn
function [vennOpts, patchOpts] = parsePVPairs (p, v, nZones)
p = lower(p);
%Break up P/V list into Venn parameters and patch parameters
vennParamNames = {'plot', 'errminmode', 'parent'};
[isVennParam, idx] = ismember(p, vennParamNames);
idx = idx(isVennParam);
%vennParams = p(isVennParam);
vennVals = v(isVennParam);
%First do Patch options
patchOpts.Parameters = p(~isVennParam);
patchOpts.Values = v(~isVennParam);
%Now do Venn options
vennOpts = defaultVennOptions;
%PLOT
i = find(idx==1, 1);
if i
plot = lower(vennVals{i});
if islogical(plot)
vennOpts.Plot = plot;
else
if ischar(plot) && any(strcmp(plot, {'on', 'off'}))
vennOpts.Plot = strcmp(plot, 'on');
else
error('venn:parsePVPairs', 'Plot must be ''on'', ''off'', or a logical value.')
end
end
end
%ERRMINMODE
i = find(idx==2, 1);
if i
mode = lower(vennVals{i});
okModes = {'None', 'TotalError', 'ChowRodgers'};
[isOkMode, modeIdx] = ismember(mode, lower(okModes));
if isOkMode
vennOpts.ErrMinMode = okModes{modeIdx};
else
error('venn:parsePVPairs', 'ErrMinMode must be None, TotalError, or ChowRodgers')
end
end
%PARENT
i = find(idx==5, 1);
if i
h = v{i};
if length(h)==1 && ishandle(h)
vennOpts.Parent = h;
else
error('venn:parsePVPairs', 'Parent must be a valid scalar handle')
end
end
end %parsePVPairs
function [A0, I0, Z0] = parseInputAreas (A0, I0, Z0)
%Switch to row vectors
A0 = A0(:)';
I0 = I0(:)';
Z0 = Z0(:)';
if isempty(Z0)
%A0 and I0 supplied
Z0 = calcZoneAreas (length(A0), A0, I0);
else
%Z0 supplied
switch length(Z0)
case 3
A0 = Z0(1:2)+Z0(3);
I0 = Z0(3);
case 7
A0 = Z0(1:3)+Z0([4 4 5])+Z0([5 6 6])+Z0(7);
I0 = [Z0(4:6)+Z0(7) Z0(7)];
otherwise
error('')
end
end
end
function vennOpts = defaultVennOptions
vennOpts = struct(...
'Plot' ,true ,...
'Labels' ,[] ,...
'PopLabels' ,false ,...
'DrawLabels' ,false ,...
'Parent' ,[] ,...
'Offset' ,[0 0] ,...
'ErrMinMode' ,'TotalError' );
end
function [d, x, y, A, I, Z] = preallocVectors (nCirc)
%Initialize position vectors
x = zeros(1, nCirc);
y = zeros(1, nCirc);
if nCirc==2
d = 0;
I = 0;
A = zeros(1,2);
Z = zeros(1,3);
else %nCirc==3
d = zeros(1,3);
I = zeros(1,4);
A = zeros(1,3);
Z = zeros(1,7);
end
endclear,clc
%Plot a simple 2-circle venn diagram with custom patch properties
figure, axis equal, axis off
A = [300 200]; I = 150;
venn(A,I,'FaceColor',{'r','y'},'FaceAlpha',{1,0.6},'EdgeColor','black')
%Compare ErrMinModes
A = [350 300 275]; I = [100 80 60 40];
figure
subplot(1,3,1), h1 = venn(A,I,'ErrMinMode','None');
axis image, title ('No 3-Circle Error Minimization')
subplot(1,3,2), h2 = venn(A,I,'ErrMinMode','TotalError');
axis image, title ('Total Error Mode')
subplot(1,3,3), h3 = venn(A,I,'ErrMinMode','ChowRodgers');
axis image, title ('Chow-Rodgers Mode')
set([h1 h2], 'FaceAlpha', 0.6)
5、颜色条
function varargout = contourfnu(x,y,data,v,cmap,pos_colorbar,overticklabel,method,ninterp,nancolor)
% non-uniform contourf/imagesc/colorbar
%
% Input variables
% x : x-coordinates of grid, vector or 2d matrix
% y : y-coordinates of grid, vector or 2d matrix
% If x and y are vectors, then length(x)==size(z,2) and length(y)==size(Z,1).
% If x and y are 2d matrix, they are generated by meshgrid
% data : 2d matrix to be ploted
%
% Optional:
% v : vector of contour levels (default:linspace(datamin,datamax,10))
% cmap : color map array (default:jet)
% pos_colorbar : 'none', or location with respect to the axes (default:'eastoutside')
% overticklabel : whether or not label the overrange ticks at colorbar (default:true)
% method : imagesc, contourf, contour or pcolor (default:imagesc)
% ninterp : repeatedly interp times in each dimension (default:0)
% nancolor : axis backgroud color
%
% Output variable
% hout : structure with handles
% .h plot handle
% .c contour matrix (method='contourf')
% .hc colorbar handle (pos_colorbar~='none')
datamax = max(data(:));
datamin = min(data(:));
if(~exist('v','var')||isempty(v)), v = linspace(datamin,datamax,10); end
if(~exist('cmap','var')||isempty(cmap)), cmap = jet; end
if(~exist('pos_colorbar','var')||isempty(pos_colorbar)), pos_colorbar = 'eastoutside'; end
if(~exist('overticklabel','var')||isempty(overticklabel)), overticklabel = true; end
if(~exist('method','var')||isempty(method)), method = 'imagesc'; end
if(~exist('ninterp','var')||isempty(ninterp)), ninterp = 0; end
if(ninterp>0)
data=interp2(data,ninterp);
if( strcmp(method,'contourf')||strcmp(method,'contour')||strcmp(method,'pcolor') )
if( isvector(x)&&isvector(y) ) % vector to meshgrid
[x,y] = meshgrid(x,y);
end
x = interp2(x,ninterp);
y = interp2(y,ninterp);
end
end
if( strcmp(method,'imagesc')&&~isvector(x) ) % meshgrid to vector
x = x(1,:);
y = y(:,1);
end
if(isrow(v)), v = v'; end
v = sort(v);
% extendmin = dataminv(end);
if(dataminv(end))
v = [v;inf];
end
nlev = length(v);
% piecewise linear mapping, v -> 1:nlev, data -> 1 to nlev
z = zeros(size(data));
for i=1:nlev-1
if( i==1 && v(i)==-inf )
index = data=v(end-1);
z(index) = i+(data(index)-v(end-1))/(datamax-v(end-1));
else
index = data>=v(i) & data0
varargout{1} = hout;
end实例代码
% Case 1
n = 50;
x = 1:n;
y = 1:n;
data = peaks(n).^4;
contour_levels = [0.5, 10, 100, 500, 1000 ];
contour_levels2 = [0,contour_levels,5000];
% % Case 2
% x=1:0.5:5;
% y=1:0.5:10;
% [x,y] = meshgrid(x,y);
% data = x.^2+exp(y);
% data(end-5:end,end-3:end) = nan;
% data(end-2:end,end-2:end) = inf;
% contour_levels = [5,20,100,1000,10000];
% contour_levels2 = [0,contour_levels,50000];
figure('position',[1,1,1150,800])
subplot(331)
contourf(x,y,data)
colorbar,colormap(jet)
title('original linear cmap')
subplot(332)
contourfnu(x,y,data,contour_levels)
title('imagesc default')
subplot(333)
contourfnu(x,y,data,contour_levels,[],[],[],'contourf')
title('contourf default')
subplot(334)
contourfnu(x,y,data,contour_levels,[],[],false)
title('imagesc overticklabel=false')
subplot(335)
cmap = jet;
cmap(1,:) = [1,1,1]; %white
contourfnu(x,y,data,contour_levels,cmap,[],false)
title('imagesc FirstColorWhite')
subplot(336)
contourfnu(x,y,data,contour_levels2)
title('imagesc AnotherContourLevels')
subplot(337)
contourfnu(x,y,data,contour_levels,[],[],[],[],2)
title('imagesc interp')
subplot(338)
contourfnu(x,y,data,[-inf,contour_levels2,inf],[],[],false)
title('imagesc InfInLevels')
subplot(339)
contourfnu(x,y,data,contour_levels2,[],'southoutside')
title('imagesc PosCbar=''southoutside''')
标签：...,end,plot,绘图,matlab,-----,circle,venn,data 来源： https://blog.51cto.com/u_12355165/2782170
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。