0001 function plot_contour(md,data,options_structure);
0002
0003
0004
0005
0006
0007
0008
0009
0010 [x y z index is2d]=processmesh(md,options_structure);
0011 [data isongrid]=processdata(md,data,options_structure);
0012
0013
0014 if ~is2d,
0015 error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
0016 end
0017
0018
0019 if ~isongrid
0020 data=averaging(md,data,0);
0021 end
0022
0023
0024 if isnan(options_structure.contouronly)
0025
0026 if isnan(options_structure.contourcolor)
0027 options_structure.contourcolor='y';
0028 end
0029 else
0030
0031 cla
0032 end
0033 color=options_structure.contourcolor;
0034
0035
0036 if isnumeric(options_structure.contourlevels),
0037 levels=round_ice(linspace(max(data),min(data),options_structure.contourlevels),2);
0038 else
0039 levels=[];
0040 for i=1:length(options_structure.contourlevels)
0041 levels(end+1)=options_structure.contourlevels{i};
0042 end
0043 levels=sort(unique(levels),'descend');
0044 end
0045 numlevels=length(levels);
0046
0047
0048 numberofelements=md.numberofelements;
0049 elementslist=1:numberofelements;
0050 c=[];
0051 h=[];
0052
0053
0054
0055 edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
0056
0057 [edges,I,J]=unique(sort(edges,2),'rows');
0058
0059 vec=J;
0060
0061
0062 edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];
0063
0064
0065 Seg1=index(:,[1 2]);
0066 Seg2=index(:,[2 3]);
0067 Seg3=index(:,[3 1]);
0068
0069
0070 Seg1_num=edges_tria(:,1);
0071 Seg2_num=edges_tria(:,2);
0072 Seg3_num=edges_tria(:,3);
0073
0074
0075 Data1=data(Seg1);
0076 Data2=data(Seg2);
0077 Data3=data(Seg3);
0078
0079
0080 Range1=sort(Data1,2);
0081 Range2=sort(Data2,2);
0082 Range3=sort(Data3,2);
0083
0084 for i=1:numlevels
0085
0086 level=levels(i);
0087
0088
0089 pos1=(Range1(:,1)<level & Range1(:,2)>level);
0090 pos2=(Range2(:,1)<level & Range2(:,2)>level);
0091 pos3=(Range3(:,1)<level & Range3(:,2)>level);
0092
0093
0094 poselem12=(pos1 & pos2);
0095 poselem13=(pos1 & pos3);
0096 poselem23=(pos2 & pos3);
0097 poselem=find(poselem12 | poselem13 | poselem23);
0098 numelems=length(poselem);
0099
0100
0101 if numelems==0,
0102 continue,
0103 end
0104
0105
0106 x1=zeros(numelems,1);
0107 x2=zeros(numelems,1);
0108 y1=zeros(numelems,1);
0109 y2=zeros(numelems,1);
0110 edge_l=zeros(numelems,2);
0111
0112 for j=1:numelems,
0113
0114 weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1));
0115 weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1));
0116 weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1));
0117
0118 if poselem12(poselem(j));
0119
0120 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
0121 x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
0122 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
0123 y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
0124 edge_l(j,1)=Seg1_num(poselem(j));
0125 edge_l(j,2)=Seg2_num(poselem(j));
0126
0127 elseif poselem13(poselem(j)),
0128
0129 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
0130 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
0131 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
0132 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
0133 edge_l(j,1)=Seg1_num(poselem(j));
0134 edge_l(j,2)=Seg3_num(poselem(j));
0135
0136 elseif poselem23(poselem(j)),
0137
0138 x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
0139 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
0140 y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
0141 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
0142 edge_l(j,1)=Seg2_num(poselem(j));
0143 edge_l(j,2)=Seg3_num(poselem(j));
0144 else
0145
0146 end
0147 end
0148
0149
0150
0151
0152 while ~isempty(edge_l),
0153
0154
0155 e1=edge_l(1,1); e2=edge_l(1,2);
0156 xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
0157
0158
0159 edge_l(1,:)=[];
0160 x1(1)=[]; x2(1)=[];
0161 y1(1)=[]; y2(1)=[];
0162
0163 [ro1,co1]=find(edge_l==e1);
0164
0165 while ~isempty(ro1)
0166
0167 if co1==1,
0168 xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
0169
0170
0171 e1=edge_l(ro1,2);
0172
0173 else
0174 xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
0175
0176
0177 e1=edge_l(ro1,1);
0178 end
0179
0180
0181 edge_l(ro1,:)=[];
0182 x1(ro1)=[]; x2(ro1)=[];
0183 y1(ro1)=[]; y2(ro1)=[];
0184
0185
0186 [ro1,co1]=find(edge_l==e1);
0187 end
0188
0189
0190 [ro2,co2]=find(edge_l==e2);
0191
0192 while ~isempty(ro2)
0193
0194 if co2==1,
0195 xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
0196
0197
0198 e2=edge_l(ro2,2);
0199 else
0200 xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
0201
0202
0203 e2=edge_l(ro2,1);
0204 end
0205
0206
0207 edge_l(ro2,:)=[];
0208 x1(ro2)=[]; x2(ro2)=[];
0209 y1(ro2)=[]; y2(ro2)=[];
0210
0211
0212 [ro2,co2]=find(edge_l==e2);
0213 end
0214
0215
0216 zc=level*ones(length(xc)+1,1);
0217 if isnan(color),
0218 h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat')];
0219 hold on
0220 else
0221 h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color)];
0222 hold on
0223 end
0224
0225
0226
0227 c = horzcat(c,[level, xc'; length(xc), yc']);
0228
0229 end
0230 end
0231
0232
0233 if ~strcmpi(options_structure.contourticks,'off')
0234 if ~isnan(options_structure.contourcolor)
0235 clabel(c,h,'color',color);
0236 else
0237 clabel(c,h);
0238 end
0239 end