Changeset 954
- Timestamp:
- Jun 22, 2016, 1:36:50 PM (9 years ago)
- Location:
- trunk/src
- Files:
-
- 11 edited
-
find_field_bounds.m (modified) (2 diffs)
-
mouse_motion.m (modified) (1 diff)
-
mouse_up.m (modified) (3 diffs)
-
plot_object.m (modified) (3 diffs)
-
proj_field.m (modified) (12 diffs)
-
series.m (modified) (4 diffs)
-
series/civ2vel_3C.m (modified) (22 diffs)
-
set_object.m (modified) (7 diffs)
-
update_imadoc.m (modified) (1 diff)
-
uvmat.fig (modified) (previous)
-
uvmat.m (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/find_field_bounds.m
r924 r954 54 54 CoordMin=zeros(numel(imax),NbDim); 55 55 Mesh=zeros(1,numel(imax)); 56 FieldOut.ProjModeRequest='projection';%default 56 57 for ind=1:numel(imax) 57 58 if strcmp(CellInfo{imax(ind)}.CoordType,'tps') … … 90 91 Mesh(ind)=min((CoordMax(ind,:)-CoordMin(ind,:))./(NbPoints-1)); 91 92 end 93 if isfield(CellInfo{imax(ind)},'ProjModeRequest') 94 if strcmp(CellInfo{imax(ind)}.ProjModeRequest,'interp_tps') 95 FieldOut.ProjModeRequest='interp_tps'; 96 end 97 if strcmp(CellInfo{imax(ind)}.ProjModeRequest,'interp_lin')&& ~strcmp(FieldOut.ProjModeRequest,'interp_tps') 98 FieldOut.ProjModeRequest='interp_lin'; 99 end 100 end 92 101 end 93 102 Mesh=min(Mesh); -
trunk/src/mouse_motion.m
r924 r954 385 385 XYData=AxeData.CurrentOrigin; 386 386 if isequal(AxeData.Drawing,'create') && isfield(AxeData,'CurrentOrigin') && ~isempty(AxeData.CurrentOrigin) 387 if strcmp(ObjectData.Type,'line')||strcmp(ObjectData.Type,'polyline')||strcmp(ObjectData.Type,'polygon')||strcmp(ObjectData.Type,'points') 388 ObjectData.Coord=[ObjectData.Coord ;xy(1,1:2)]; 389 % ObjectData.Coord(end,:)=xy(1,:); 390 elseif strcmp(ObjectData.Type,'rectangle')||strcmp(ObjectData.Type,'ellipse')||strcmp(ObjectData.Type,'volume') 391 ObjectData.Coord=(AxeData.CurrentOrigin+xy(1,1:2))/2;% keep only the first point coordinate 387 switch ObjectData.Type 388 case {'line','polyline','polygon','points','plane_z'} 389 ObjectData.Coord=[ObjectData.Coord ;xy(1,1:2)]; 390 % ObjectData.Coord(end,:)=xy(1,:); 391 case {'rectangle','ellipse','volume'} 392 ObjectData.Coord=(AxeData.CurrentOrigin+xy(1,1:2))/2;% keep only the first point coordinate 392 393 ObjectData.RangeX=abs(ObjectData.Coord(1,1)-xy(1,1));%rectangle width 393 ObjectData.RangeY=abs(ObjectData.Coord(1,2)-xy(1,2));%rectangle height 394 elseif isequal(ObjectData.Type,'plane')%case of 'plane'395 DX=(xy(1,1)-ObjectData.Coord(1,1));396 DY=(xy(1,2)-ObjectData.Coord(1,2));397 ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle widt398 if isfield(ObjectData,'RangeX')399 XMax=sqrt(DX*DX+DY*DY);400 if XMax>max(ObjectData.RangeX)401 ObjectData.RangeX=[min(ObjectData.RangeX) XMax];402 end403 end394 ObjectData.RangeY=abs(ObjectData.Coord(1,2)-xy(1,2));%rectangle height 395 case 'plane' %case of 'plane' 396 DX=(xy(1,1)-ObjectData.Coord(1,1)); 397 DY=(xy(1,2)-ObjectData.Coord(1,2)); 398 ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle widt 399 if isfield(ObjectData,'RangeX') 400 XMax=sqrt(DX*DX+DY*DY); 401 if XMax>max(ObjectData.RangeX) 402 ObjectData.RangeX=[min(ObjectData.RangeX) XMax]; 403 end 404 end 404 405 end 405 406 plot_object(ObjectData,ProjObject,AxeData.CurrentObject,'m'); -
trunk/src/mouse_up.m
r924 r954 95 95 else 96 96 switch ObjectData.Type 97 case {'line' }97 case {'line','plane_z'} 98 98 if size(ObjectData.Coord,1)==1 % this is the mouse up for the first point, continue until next click 99 99 check_multiple=1; … … 106 106 check_multiple=1;% pass to next mous up if width of height=0 107 107 end 108 case 'plane' %case of 'plane' 108 case 'plane' %case of 'plane', TODO: NOT ACTIVATED 109 109 DX=(xy(1,1)-ObjectData.Coord(1,1)); 110 110 DY=(xy(1,2)-ObjectData.Coord(1,2)); … … 116 116 end 117 117 end 118 case 'plane_z' 119 if size(ObjectData.Coord,1)==1 % this is the mouse up for the first point, continue until next click 120 check_multiple=1; 121 end 122 DX=(xy(1,1)-ObjectData.Coord(1,1)); 123 DY=(xy(1,2)-ObjectData.Coord(1,2)); 124 ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle width 118 125 otherwise 119 126 check_multiple=1; -
trunk/src/plot_object.m
r937 r954 133 133 134 134 %% determine the coordinates xline, yline,xsup,xinf, yinf,ysup determining the new object plot 135 test_line= isequal(ObjectData.Type,'points')||isequal(ObjectData.Type,'line')||... 136 isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon')|| isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume'); 137 test_patch=isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')||isequal(ObjectData.Type,'volume')... 138 ||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'); 135 %test_line= isequal(ObjectData.Type,'points')||isequal(ObjectData.Type,'line')||... 136 % isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon')|| isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume'); 137 %test_patch=isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')||isequal(ObjectData.Type,'volume')... 138 % ||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'); 139 test_line=ismember(ObjectData.Type,{'points','line','polyline','polygon','plane','plane_z','volume'}); 140 test_patch=ismember(ObjectData.ProjMode,{'inside','outside','mask_inside','mask_outside'}); 139 141 if test_line 140 142 xline=ObjectData.Coord(:,1); 141 143 yline=ObjectData.Coord(:,2); 142 144 nbpoints=numel(xline); 143 if isequal(ObjectData.Type,'line_x') 144 xline=[xline; ObjectData.RangeX(2)];%creating the line 145 yline=[yline; ObjectData.RangeY(2)];%creating the line 146 elseif isequal(ObjectData.Type,'polygon') 147 xline=[xline; ObjectData.Coord(1,1)];%closing the line 148 yline=[yline; ObjectData.Coord(1,2)]; 149 elseif isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume') 150 if ~isfield(ObjectData,'Angle') 151 ObjectData.Angle=[0 0 0]; 152 end 153 phi=ObjectData.Angle(3)*pi/180;%angle in radians 154 x0=xline(1); y0=yline(1); 155 xlim=get(haxes,'XLim'); 156 ylim=get(haxes,'YLim'); 157 graph_scale=max(abs(xlim(2)-xlim(1)),abs(ylim(2)-ylim(1)))/2;% estimate the length of axes plots 158 XMax=graph_scale; 159 YMax=graph_scale; 160 XMin=-graph_scale; 161 YMin=-graph_scale; 162 if ~isempty(XMaxRange) 163 XMax=XMaxRange; 164 end 165 if ~isempty(XMinRange) 166 XMin=XMinRange; 167 end 168 if ~isempty(YMaxRange) 169 YMax=YMaxRange; 170 end 171 if ~isempty(YMinRange) 172 YMin=YMinRange; 173 end 174 % axes lines 175 xline=NaN(1,13); 176 xline(1)=x0+min(0,XMin)*cos(phi); % min end of the x axes 177 yline(1)=y0+min(0,XMin)*sin(phi); 178 xline(2)=x0+XMax*cos(phi);% max end of the x axes 179 yline(2)=y0+XMax*sin(phi); 180 xline(8)=x0-min(0,YMin)*sin(phi);% min end of the y axes 181 yline(8)=y0+min(0,YMin)*cos(phi); 182 xline(9)=x0-YMax*sin(phi);% max end of the y axes 183 yline(9)=y0+YMax*cos(phi); 184 185 %arrows on x axis 186 arrow_scale=graph_scale/20; 187 xline(3)=xline(2)-arrow_scale*cos(phi-pi/8); 188 yline(3)=yline(2)-arrow_scale*sin(phi-pi/8); 189 xline(5)=xline(2); 190 yline(5)=yline(2); 191 xline(6)=xline(2)-arrow_scale*cos(phi+pi/8); 192 yline(6)=yline(2)-arrow_scale*sin(phi+pi/8); 193 194 %arrows on y axis 195 xline(10)=xline(9)-arrow_scale*cos(phi+pi/2-pi/8); 196 yline(10)=yline(9)-arrow_scale*sin(phi+pi/2-pi/8); 197 xline(12)=xline(9); 198 yline(12)=yline(9); 199 xline(13)=xline(9)-arrow_scale*cos(phi+pi/2+pi/8); 200 yline(13)=yline(9)-arrow_scale*sin(phi+pi/2+pi/8); 201 %xline=[Xbeg_x Xend_x NaN Ybeg_x Yend_x]; 202 %yline=[Xbeg_y Xend_y NaN Ybeg_y Yend_y]; 203 % dashed lines indicating bounds 204 xsup=NaN(1,5); 205 ysup=NaN(1,5); 206 if ~isempty(XMaxRange) 207 xsup(1)=xline(2)-YMin*sin(phi); 208 ysup(1)=yline(2)+YMin*cos(phi); 209 xsup(2)=xline(2)-YMax*sin(phi); 210 ysup(2)=yline(2)+YMax*cos(phi); 211 end 212 if ~isempty(YMaxRange) 213 xsup(2)=xline(2)-YMax*sin(phi); 214 ysup(2)=yline(2)+YMax*cos(phi); 215 xsup(3)=xline(9)+XMin*cos(phi); 216 ysup(3)=yline(9)+XMin*sin(phi); 217 end 218 if ~isempty(XMinRange) 219 xsup(3)=xline(9)+XMin*cos(phi); 220 ysup(3)=yline(9)+XMin*sin(phi); 221 xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); 222 ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); 223 end 224 if ~isempty(YMinRange) 225 xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); 226 ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); 227 xsup(5)=xline(8)-YMin*sin(phi); 228 ysup(5)=yline(8)+YMin*cos(phi); 229 end 230 end 231 SubLineStyle='none';%default 232 if isfield(ObjectData,'ProjMode') 233 if isequal(ObjectData.ProjMode,'projection') 234 SubLineStyle='--'; %range of projection marked by dash 235 if isfield (ObjectData,'DX') 236 ObjectData=rmfield(ObjectData,'DX'); 237 end 238 if isfield (ObjectData,'DY') 239 ObjectData=rmfield(ObjectData,'DY'); 240 end 241 elseif isequal(ObjectData.ProjMode,'filter') 242 SubLineStyle=':';%range of projection not visible 243 end 244 end 245 if isequal(ObjectData.Type,'line')||isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon') 246 if length(xline)<2 247 theta=0; 248 else 249 theta=angle(diff(xline)+1i*diff(yline)); 250 theta(length(xline))=theta(length(xline)-1); 251 end 252 xsup(1)=xline(1)+YMax*sin(theta(1)); 253 xinf(1)=xline(1)-YMax*sin(theta(1)); 254 ysup(1)=yline(1)-YMax*cos(theta(1)); 255 yinf(1)=yline(1)+YMax*cos(theta(1)); 256 for ip=2:length(xline) 257 xsup(ip)=xline(ip)+YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 258 xinf(ip)=xline(ip)-YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 259 ysup(ip)=yline(ip)-YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 260 yinf(ip)=yline(ip)+YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 261 end 145 switch ObjectData.Type 146 case 'line_x' 147 xline=[xline; ObjectData.RangeX(2)];%creating the line 148 yline=[yline; ObjectData.RangeY(2)];%creating the line 149 case 'polygon' 150 xline=[xline; ObjectData.Coord(1,1)];%closing the line 151 yline=[yline; ObjectData.Coord(1,2)]; 152 case {'plane','volume'} 153 if ~isfield(ObjectData,'Angle') 154 ObjectData.Angle=[0 0 0]; 155 end 156 phi=ObjectData.Angle(3)*pi/180;%angle in radians 157 x0=xline(1); y0=yline(1); 158 xlim=get(haxes,'XLim'); 159 ylim=get(haxes,'YLim'); 160 graph_scale=max(abs(xlim(2)-xlim(1)),abs(ylim(2)-ylim(1)))/2;% estimate the length of axes plots 161 XMax=graph_scale; 162 YMax=graph_scale; 163 XMin=-graph_scale; 164 YMin=-graph_scale; 165 if ~isempty(XMaxRange) 166 XMax=XMaxRange; 167 end 168 if ~isempty(XMinRange) 169 XMin=XMinRange; 170 end 171 if ~isempty(YMaxRange) 172 YMax=YMaxRange; 173 end 174 if ~isempty(YMinRange) 175 YMin=YMinRange; 176 end 177 % axes lines 178 xline=NaN(1,13); 179 xline(1)=x0+min(0,XMin)*cos(phi); % min end of the x axes 180 yline(1)=y0+min(0,XMin)*sin(phi); 181 xline(2)=x0+XMax*cos(phi);% max end of the x axes 182 yline(2)=y0+XMax*sin(phi); 183 xline(8)=x0-min(0,YMin)*sin(phi);% min end of the y axes 184 yline(8)=y0+min(0,YMin)*cos(phi); 185 xline(9)=x0-YMax*sin(phi);% max end of the y axes 186 yline(9)=y0+YMax*cos(phi); 187 188 %arrows on x axis 189 arrow_scale=graph_scale/20; 190 xline(3)=xline(2)-arrow_scale*cos(phi-pi/8); 191 yline(3)=yline(2)-arrow_scale*sin(phi-pi/8); 192 xline(5)=xline(2); 193 yline(5)=yline(2); 194 xline(6)=xline(2)-arrow_scale*cos(phi+pi/8); 195 yline(6)=yline(2)-arrow_scale*sin(phi+pi/8); 196 197 %arrows on y axis 198 xline(10)=xline(9)-arrow_scale*cos(phi+pi/2-pi/8); 199 yline(10)=yline(9)-arrow_scale*sin(phi+pi/2-pi/8); 200 xline(12)=xline(9); 201 yline(12)=yline(9); 202 xline(13)=xline(9)-arrow_scale*cos(phi+pi/2+pi/8); 203 yline(13)=yline(9)-arrow_scale*sin(phi+pi/2+pi/8); 204 %xline=[Xbeg_x Xend_x NaN Ybeg_x Yend_x]; 205 %yline=[Xbeg_y Xend_y NaN Ybeg_y Yend_y]; 206 % dashed lines indicating bounds 207 xsup=NaN(1,5); 208 ysup=NaN(1,5); 209 if ~isempty(XMaxRange) 210 xsup(1)=xline(2)-YMin*sin(phi); 211 ysup(1)=yline(2)+YMin*cos(phi); 212 xsup(2)=xline(2)-YMax*sin(phi); 213 ysup(2)=yline(2)+YMax*cos(phi); 214 end 215 if ~isempty(YMaxRange) 216 xsup(2)=xline(2)-YMax*sin(phi); 217 ysup(2)=yline(2)+YMax*cos(phi); 218 xsup(3)=xline(9)+XMin*cos(phi); 219 ysup(3)=yline(9)+XMin*sin(phi); 220 end 221 if ~isempty(XMinRange) 222 xsup(3)=xline(9)+XMin*cos(phi); 223 ysup(3)=yline(9)+XMin*sin(phi); 224 xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); 225 ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); 226 end 227 if ~isempty(YMinRange) 228 xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); 229 ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); 230 xsup(5)=xline(8)-YMin*sin(phi); 231 ysup(5)=yline(8)+YMin*cos(phi); 232 end 233 end 234 end 235 SubLineStyle='none';%default 236 if isfield(ObjectData,'ProjMode') 237 if isequal(ObjectData.ProjMode,'projection') 238 SubLineStyle='--'; %range of projection marked by dash 239 if isfield (ObjectData,'DX') 240 ObjectData=rmfield(ObjectData,'DX'); 241 end 242 if isfield (ObjectData,'DY') 243 ObjectData=rmfield(ObjectData,'DY'); 244 end 245 elseif isequal(ObjectData.ProjMode,'filter') 246 SubLineStyle=':';%range of projection not visible 247 end 248 end 249 if ismember(ObjectData.Type,{'line','polyline','polygon','plane_z'}) 250 if length(xline)<2 251 theta=0; 252 else 253 theta=angle(diff(xline)+1i*diff(yline)); 254 theta(length(xline))=theta(length(xline)-1); 255 end 256 xsup(1)=xline(1)+YMax*sin(theta(1)); 257 xinf(1)=xline(1)-YMax*sin(theta(1)); 258 ysup(1)=yline(1)-YMax*cos(theta(1)); 259 yinf(1)=yline(1)+YMax*cos(theta(1)); 260 for ip=2:length(xline) 261 xsup(ip)=xline(ip)+YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 262 xinf(ip)=xline(ip)-YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 263 ysup(ip)=yline(ip)-YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 264 yinf(ip)=yline(ip)+YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); 262 265 end 263 266 end … … 415 418 else% no patch image requested, erase existing ones 416 419 if isfield(PlotData,'SubObject') 417 for iobj=1:length(PlotData.SubObject)418 if ishandle(PlotData.SubObject(iobj)) && strcmp(get(PlotData.SubObject(iobj),'Type'),'image')419 delete(PlotData.SubObject(iobj))420 end421 end420 for iobj=1:length(PlotData.SubObject) 421 if ishandle(PlotData.SubObject(iobj)) && strcmp(get(PlotData.SubObject(iobj),'Type'),'image') 422 delete(PlotData.SubObject(iobj)) 423 end 424 end 422 425 end 423 426 end … … 459 462 end 460 463 end 461 case {'line','polyline','polygon' }464 case {'line','polyline','polygon','plane_z'} 462 465 hh=line(xline,yline,'Color',col); 463 466 PlotData.SubObject(1)=line(xinf,yinf,'Color',col,'LineStyle',SubLineStyle,'Tag','proj_object');%draw sub-lines -
trunk/src/proj_field.m
r937 r954 90 90 return 91 91 end 92 ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes93 if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case92 % check list of effective projection modes 93 if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'}) 94 94 return 95 95 end … … 117 117 [ProjData,errormsg] = proj_line(FieldData,ObjectData); 118 118 end 119 case 'plane'119 case {'plane','plane_z'} 120 120 [ProjData,errormsg] = proj_plane(FieldData,ObjectData); 121 121 case 'volume' … … 955 955 test90x=0;%=1 for 90 degree rotation alround x axis 956 956 test90y=0;%=1 for 90 degree rotation alround y axis 957 if strcmp(ObjectData.Type,'plane_z') 958 Delta_x=ObjectData.Coord(2,1)-ObjectData.Coord(1,1); 959 Delta_y=ObjectData.Coord(2,2)-ObjectData.Coord(1,2); 960 Delta_mod=sqrt(Delta_x*Delta_x+Delta_y*Delta_y); 961 ObjectData.Angle=[0 0 0]; 962 ObjectData.Angle(1)=90*Delta_x/Delta_mod; 963 ObjectData.Angle(2)=90*Delta_y/Delta_mod; 964 end 957 965 if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) 958 966 test90y=isequal(ObjectData.Angle,[0 90 0]); … … 987 995 return 988 996 end 997 InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane 998 if strcmp(ObjectData.Type,'plane_z') 999 InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation 1000 end 989 1001 990 1002 %% extrema along each axis … … 1071 1083 end 1072 1084 icell_grid=[];% field cell index which defines the grid 1073 if ~strcmp(ObjectData.ProjMode,'projection') 1085 if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize 1074 1086 %% define the new coordinates in case of interpolation on a imposed grid 1075 if ~testYMin 1087 if ~testYMin 1076 1088 errormsg='min Y value not defined for the projection grid';return 1077 1089 end … … 1119 1131 AYName='coord_y'; 1120 1132 AXName='coord_x'; 1121 if strcmp(ObjectData.ProjMode,'projection') 1133 if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z') 1122 1134 ProjData.coord_y=[FieldData.YMin FieldData.YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined 1123 1135 ProjData.coord_x=[FieldData.XMin FieldData.XMax]; … … 1399 1411 DimValue(DimValue==1)=[];%remove singleton dimensions 1400 1412 NbDim=numel(DimValue);%update number of space dimensions 1401 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate1413 % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate 1402 1414 if NbDim>=3 1403 1415 if NbDim>3 … … 1439 1451 Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)}); 1440 1452 end 1441 if numel(Coord{NbDim-1})==2 1453 if numel(Coord{NbDim-1})==2% case of coordinate defined only by the first and last values 1442 1454 DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1); 1443 1455 end 1444 if numel(Coord{NbDim})==2 1456 if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values 1445 1457 DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1); 1446 1458 end … … 1453 1465 end 1454 1466 else 1455 YIndexMax= Coord{NbDim-1}(end)/DY;1467 YIndexMax=numel(Coord{NbDim-1}); 1456 1468 YIndexMin=1; 1457 1469 end 1458 1470 if testXMax 1459 XIndexMax=(XMax-Coord{NbDim}(1))/D Y+1;% matrix index corresponding to the max y value for the new field1471 XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field 1460 1472 if testYMin%test_direct(indY) 1461 1473 XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field … … 1464 1476 end 1465 1477 else 1466 XIndexMax= Coord{NbDim}(end)/DX;1478 XIndexMax=numel(Coord{NbDim}); 1467 1479 XIndexMin=1; 1468 1480 end … … 1627 1639 end 1628 1640 end 1629 else 1630 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 1631 %TODO: use interp_lin3 1632 return 1641 else %projection of structured coordinates on oblique plane 1642 % determine the boundaries of the projected field, 1643 % first find the 8 summits of the initial volume in the 1644 % new coordinates 1645 Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates 1646 Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates 1647 Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates 1648 summit=zeros(3,8);% initialize summit coordinates 1649 summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square 1650 summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1) 1651 summit(1:2,5:8)=summit(1:2,1:4); 1652 summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)]; 1653 Mrot=rodrigues(PlaneAngle); 1654 newsummit=zeros(3,8);% initialize the rotated summit coordinates 1655 ObjectData.Coord=[ObjectData.Coord(1,1); ObjectData.Coord(1,2); 0];%TODO: set in set_object 1656 for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data) 1657 newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord)); 1658 end 1659 coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:)); 1660 coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:)); 1661 coord_z_proj=-width:InterpMesh:width; 1662 Mrot_inverse=rodrigues(-PlaneAngle); 1663 Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord; 1664 ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0];% unit vector along the new x coordinates 1665 iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0];% unit vector y coordinates 1666 iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)];% x unit vector z coordinates 1667 [Grid_x,Grid_y,Grid_z]=meshgrid(1:numel(coord_x_proj),1:numel(coord_y_proj),1:numel(coord_z_proj)); 1668 if ismatrix(Grid_x) 1669 Grid_x=shiftdim(Grid_x,-1); 1670 Grid_y=shiftdim(Grid_y,-1); 1671 Grid_z=shiftdim(Grid_z,-1); 1672 end 1673 XI=Origin(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z; 1674 YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z; 1675 ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z; 1676 [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1}); 1677 X=permute(X,[3 1 2]); 1678 Y=permute(Y,[3 1 2]); 1679 Z=permute(Z,[3 1 2]); 1680 for ivar=VarIndex 1681 VarName=FieldData.ListVarName{ivar}; 1682 ListVarName=[ListVarName VarName]; 1683 VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1684 ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI); 1685 end 1633 1686 end 1634 1687 end … … 2264 2317 end 2265 2318 else 2319 RotMatrix=rodrigues(om); 2320 2266 2321 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 2267 2322 %TODO: use interp3 -
trunk/src/series.m
r951 r954 2292 2292 CheckList_1=1;% indicate whether FieldName_1 has been updated 2293 2293 handles_coord=[handles.Coord_x handles.Coord_y handles.Coord_z handles.Coord_x_title handles.Coord_y_title handles.Coord_z_title]; 2294 if VelTypeRequest && numel(iview_civ)>=1 2294 if VelTypeRequest && numel(iview_civ)>=1 2295 2295 menu=set_veltype_display(SeriesData.FileInfo{iview_civ(1)}.CivStage,SeriesData.FileType{iview_civ(1)}); 2296 2296 set(handles.VelType,'Value',1)% set first choice by default … … 2301 2301 %CheckList=1; 2302 2302 set(handles.FieldName,'Value',1); %velocity vector choice by default 2303 if VelTypeRequest_1 && numel(iview_civ)>=2 2303 if VelTypeRequest_1 && numel(iview_civ)>=2 2304 2304 menu=set_veltype_display(SeriesData.FileInfo{iview_civ(2)}.CivStage,SeriesData.FileType{iview_civ(2)}); 2305 2305 set(handles.VelType_1,'Value',1)% set first choice by default … … 2317 2317 set(handles.VelType,'Visible','off') 2318 2318 set(handles.VelType_title,'Visible','off') 2319 end 2319 end 2320 2320 2321 2321 %% Detect the types of input files and set menus and default options in 'FieldName' 2322 2322 if (FieldNameRequest || VelTypeRequest) && numel(iview_netcdf)>=1 2323 set(handles.InputFields,'Visible','on') 2323 set(handles.InputFields,'Visible','on')% set the frame InputFields visible 2324 2324 if FieldNameRequest && isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName') 2325 2325 set(handles.FieldName,'Visible','on') … … 2346 2346 end 2347 2347 end 2348 set(handles_coord,'Visible','on') 2349 FieldList=[FieldList;{'get_field...'}]; 2350 if FieldNameRequest_1 && numel(iview_netcdf)>=2 2348 end 2349 2350 set(handles_coord,'Visible','on') 2351 FieldList=[FieldList;{'get_field...'}]; 2352 if FieldNameRequest_1 && numel(iview_netcdf)>=2 2353 set(handles.FieldName_1,'Visible','on') 2354 if CheckList_1==0 % not civ input made 2355 ListVarName=SeriesData.FileInfo{iview_netcdf(2)}.ListVarName; 2356 ind_var=get(handles.FieldName,'Value');%indices of previously selected variables 2357 for ilist=1:numel(ind_var) 2358 if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName))) 2359 FieldList_1={};% previous choice not consistent with new input field 2360 set(handles.FieldName_1,'Value',1) 2361 break 2362 end 2363 end 2364 warn_coord=0; 2365 if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||... 2366 isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName))) 2367 warn_coord=1; 2368 end 2369 if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName))) 2370 FieldList_1={}; 2371 warn_coord=1; 2372 end 2373 if warn_coord 2374 msgbox_uvmat('WARNING','coordiante names do not exist in the second netcdf input file') 2375 end 2376 2351 2377 set(handles.FieldName_1,'Visible','on') 2352 if CheckList_1==0 % not civ input made 2353 ListVarName=SeriesData.FileInfo{iview_netcdf(2)}.ListVarName; 2354 ind_var=get(handles.FieldName,'Value');%indices of previously selected variables 2355 for ilist=1:numel(ind_var) 2356 if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName))) 2357 FieldList_1={};% previous choice not consistent with new input field 2358 set(handles.FieldName_1,'Value',1) 2359 break 2360 end 2361 end 2362 warn_coord=0; 2363 if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||... 2364 isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName))) 2365 warn_coord=1; 2366 end 2367 if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName))) 2368 FieldList_1={}; 2369 warn_coord=1; 2370 end 2371 if warn_coord 2372 msgbox_uvmat('WARNING','coordiante names do not exist in the second netcdf input file') 2373 end 2374 set(handles.FieldName,'String',[FieldList;{'get_field...'}]) 2375 set(handles.FieldName_1,'Visible','on') 2376 set(handles.FieldName_1,'Value',1) 2377 set(handles.FieldName_1,'String',FieldList_1) 2378 end 2379 else 2380 set(handles.FieldName_1,'Visible','off') 2378 set(handles.FieldName_1,'Value',1) 2379 set(handles.FieldName_1,'String',FieldList_1) 2381 2380 end 2382 2381 else 2382 set(handles.FieldName_1,'Visible','off') 2383 end 2384 if isempty(FieldList) 2383 2385 set(handles.FieldName,'Visible','off') 2384 end 2385 set(handles.FieldName,'String',FieldList) 2386 else 2387 set(handles.FieldName,'Visible','on') 2388 set(handles.FieldName,'String',FieldList) 2389 end 2386 2390 else 2387 2391 set(handles.InputFields,'Visible','off') -
trunk/src/series/civ2vel_3C.m
r927 r954 1 %'civ2vel_3C': combine velocity fields from two cameras to get three velocity components 1 %'civ2vel_3C': combine velocity fields from two cameras to get the three velocity components 2 % used with the GUI 'series': 3 % first input line =raw PIV camera 1 (image coordinates) 4 % second input line=raw PIV camera 2 (image coordinates) 5 % Three modes: 6 % 1) no additional input: measurements assumed in the reference plane (laser sheet) 7 % 2) measurement surface obtained by stereoscopic comparison of the images of the two cameras. 8 % third input line =surface z(x,y) given by correlation between camera 1 and 2 (expressed in phys apparent coordinates) 9 % 3) surface z(x,y) given by adding the displacements of each camera with a third intermediate camera (#3) used to reduce the 10 % to reduce the angle for stereoscopic view. 11 % third input line =correlation between camera 1 and 3 (expressed in phys apparent coordinates) 12 % fourth input line =corelation between camera 2 and 3 (expressed in phys apparent coordinates) 13 % 2 14 %------------------------------------------------------------------------ 3 15 % function ParamOut=civ2vel_3C(Param) … … 7 19 % 8 20 %INPUT: 21 % 9 22 % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. 10 23 % In batch mode, Param is the name of the corresponding xml file containing the same information … … 15 28 % see the current structure Param) 16 29 % .InputTable: cell of input file names, (several lines for multiple input) 17 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} 30 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}. 31 % 3 or 4 lines used as described above 18 32 % .OutputSubDir: name of the subdirectory for data outputs 19 33 % .OutputDirExt: directory extension for data outputs … … 22 36 % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct 23 37 % .RUN =0 for GUI input, =1 for function activation 24 % .RunMode='local','background', 'cluster': type of function use 25 % 38 % .RunMode='local','background', 'cluster': type of function use 26 39 % .IndexRange: set the file or frame indices on which the action must be performed 27 40 % .InputFields: sub structure describing the input fields withfields … … 52 65 53 66 function ParamOut=civ2vel_3C(Param) 54 disp('test') 67 55 68 %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed 56 69 if isstruct(Param) && isequal(Param.Action.RUN,0) … … 124 137 %% define the directory for result file (with path=RootPath{1}) 125 138 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 126 %127 % if ~isfield(Param,'InputFields')128 % Param.InputFields.FieldName='';129 % end130 139 131 140 %% calibration data and timing: read the ImaDoc files … … 160 169 return 161 170 end 162 ObjectData=Param.ProjObject; 171 ObjectData=Param.ProjObject;% Object attached to the GUI series 163 172 xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); 164 173 yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); … … 168 177 W=zeros(size(XI,1),size(XI,2)); 169 178 170 %% MAIN LOOP ON FIELDS179 %%%% MAIN LOOP ON FIELDS %%%%% 171 180 warning off 172 181 … … 175 184 CheckOverwrite=Param.CheckOverwrite; 176 185 end 177 for index=1:NbField 178 179 update_waitbar(WaitbarHandle,index/NbField) 180 181 182 183 184 %% generating the name of the merged field 185 i1=i1_series{1}(index); 186 for field_index=1:NbField 187 188 update_waitbar(WaitbarHandle,field_index/NbField)% waitbar to visualise progress (in RUN mode) 189 190 %% generating the name of the output file for the merged field 191 i1=i1_series{1}(field_index); 186 192 if ~isempty(i2_series{end}) 187 i2=i2_series{end}( index);193 i2=i2_series{end}(field_index); 188 194 else 189 195 i2=i1; … … 192 198 j2=1; 193 199 if ~isempty(j1_series{1}) 194 j1=j1_series{1}( index);200 j1=j1_series{1}(field_index); 195 201 if ~isempty(j2_series{end}) 196 j2=j2_series{end}( index);202 j2=j2_series{end}(field_index); 197 203 else 198 204 j2=j1; … … 201 207 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2); 202 208 203 %% 204 205 209 %% program stop if requested on the GUI 206 210 if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') 207 211 disp('program stopped by user') … … 209 213 end 210 214 211 if (~CheckOverwrite && exist(OutputFile,'file')) 212 disp('existing output file already exists, skip to next field') 213 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 214 end 215 215 %% check existence of the output file 216 if (~CheckOverwrite && exist(OutputFile,'file')) 217 disp('existing output file already exists, skip to next field') 218 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 219 end 220 216 221 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 217 222 Data=cell(1,NbView);%initiate the set Data 218 timeread=zeros(1,NbView);219 223 220 224 %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken 221 225 %at the middle between to time step 222 clear ZItemp 223 ZItemp=zeros(size(XI,1),size(XI,2),2); 224 225 if index==1 226 ZItemp=zeros(size(XI,1),size(XI,2),2); 227 if field_index==1 226 228 first_img=i1_series{1,1}(1,1); %id of the first image of the series 227 end 228 229 idtemp=0; 230 for indextemp=index:index+1; 231 idtemp=idtemp+1; 232 if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys 233 234 235 236 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']); 237 238 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 239 temp=find(Data{3}.Civ3_FF==0); 240 Zphys=Data{3}.Zphys(temp); 241 Yphys=Data{3}.Yphys(temp); 242 Xphys=Data{3}.Xphys(temp); 243 else 244 Zphys=Data{3}.Zphys; 245 Yphys=Data{3}.Yphys; 246 Xphys=Data{3}.Xphys; 229 end 230 231 idtemp=0; 232 % get the surface shape corresponding to the PIV measurements 233 for indextemp=field_index:field_index+1;%TODO: generalise to field index intervals>1 for PIV 234 idtemp=idtemp+1; 235 InputFile_3=fullfile(Param.InputTable{3,1},Param.InputTable{3,2},[Param.InputTable{3,3} '_' int2str(first_img+indextemp-1) '.nc']); 236 if NbView==3 % if there is only 1 stereo folder (2 cameras only), extract directly Xphys,Yphys and Zphys 237 % Data{1}: =raw PIV camera 1 only 238 % Data{2}: =raw PIV camera 2 only 239 % Data{3}: =correlation between camera 1 and 2 240 [Data{3},tild,errormsg] = nc2struct(InputFile_3);%read input file 241 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 242 temp=find(Data{3}.Civ3_FF==0); 243 Zphys=Data{3}.Zphys(temp); 244 Yphys=Data{3}.Yphys(temp); 245 Xphys=Data{3}.Xphys(temp); 246 else 247 Zphys=Data{3}.Zphys; 248 Yphys=Data{3}.Yphys; 249 Xphys=Data{3}.Xphys; 250 end 251 252 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 253 % Data{1}: =raw PIV camera 1 only (left) 254 % Data{2}: =raw PIV camera 2 only (right) (no PIV done with middle camera) 255 % Data{3}: =corelation between camera 1 and 3 (left and middle) 256 % Data{4}: =corelation between camera 2 and 3 (right and middle) 257 % test if the second camera (3) is the same for both folders 258 for i=3:4 259 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 260 indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1 261 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name 262 end 263 if strcmp(camname{3},camname{4})==0 264 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 265 return 266 end 267 [Data{3},tild,errormsg] = nc2struct(InputFile_3); 268 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 269 temp=find(Data{3}.Civ3_FF==0); 270 Xmid3=Data{3}.Xmid(temp); 271 Ymid3=Data{3}.Ymid(temp); 272 U3=Data{3}.Uphys(temp); 273 V3=Data{3}.Vphys(temp); 274 else 275 Xmid3=Data{3}.Xmid; 276 Ymid3=Data{3}.Ymid; 277 U3=Data{3}.Uphys; 278 V3=Data{3}.Vphys; 279 end 280 %temporary grid of merging the 2 stereos data 281 [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)-min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)-min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2)); 282 283 %1st folder : interpolate the first camera points on the second (common) camera 284 %(Dalsa 3) 285 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq); 286 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq); 287 288 InputFile_4=fullfile(Param.InputTable{4,1},Param.InputTable{4,2},[Param.InputTable{4,3} '_' int2str(first_img+indextemp-1) '.nc']); 289 [Data{4},tild,errormsg] = nc2struct(InputFile_4); 290 if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector 291 temp=find(Data{4}.Civ3_FF==0); 292 Xmid4=Data{4}.Xmid(temp); 293 Ymid4=Data{4}.Ymid(temp); 294 U4=Data{4}.Uphys(temp); 295 V4=Data{4}.Vphys(temp); 296 else 297 Xmid4=Data{4}.Xmid; 298 Ymid4=Data{4}.Ymid; 299 U4=Data{4}.Uphys; 300 V4=Data{4}.Vphys; 301 end 302 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 303 %(Dalsa 3) 304 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq); 305 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq); 306 307 %add the displacements of the two camera pairs 308 xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); 309 ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); 310 u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1); 311 v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1); 312 313 % get the surface z(x,y) from the combined displacement 314 [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys 315 %remove NaN 316 tempNaN=isnan(Zphys);tempind=find(tempNaN==1); 317 Zphys(tempind)=[]; 318 Xphys(tempind)=[]; 319 Yphys(tempind)=[]; 320 error(tempind)=[]; 321 247 322 end 248 323 324 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid 249 325 250 251 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 252 253 254 %test if the seconde camera is the same for both folder 255 for i=3:4 256 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 257 indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1 258 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name 259 end 260 261 if strcmp(camname{3},camname{4})==0 262 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 263 return 264 end 265 266 267 268 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']); 269 270 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 271 temp=find(Data{3}.Civ3_FF==0); 272 Xmid3=Data{3}.Xmid(temp); 273 Ymid3=Data{3}.Ymid(temp); 274 U3=Data{3}.Uphys(temp); 275 V3=Data{3}.Vphys(temp); 276 else 277 Xmid3=Data{3}.Xmid; 278 Ymid3=Data{3}.Ymid; 279 U3=Data{3}.Uphys; 280 V3=Data{3}.Vphys; 281 end 282 %temporary gridd of merging the 2 stereos datas 283 [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)-min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)-min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2)); 284 285 %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera 286 %(Dalsa 3) 287 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq); 288 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq); 289 290 291 292 [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp-1),'.nc']); 293 if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector 294 temp=find(Data{4}.Civ3_FF==0); 295 Xmid4=Data{4}.Xmid(temp); 296 Ymid4=Data{4}.Ymid(temp); 297 U4=Data{4}.Uphys(temp); 298 V4=Data{4}.Vphys(temp); 299 else 300 Xmid4=Data{4}.Xmid; 301 Ymid4=Data{4}.Ymid; 302 U4=Data{4}.Uphys; 303 V4=Data{4}.Vphys; 304 end 305 306 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 307 %(Dalsa 3) 308 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq); 309 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq); 310 311 xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); 312 ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); 313 u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1); 314 v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1); 315 316 317 [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys 318 %remove NaN 319 tempNaN=isnan(Zphys);tempind=find(tempNaN==1); 320 Zphys(tempind)=[]; 321 Xphys(tempind)=[]; 322 Yphys(tempind)=[]; 323 error(tempind)=[]; 324 325 end 326 327 328 329 330 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd 331 332 end 333 ZI=mean(ZItemp,3); %mean between two the two time step 326 end 327 ZI=mean(ZItemp,3); %mean between two the two times used for surface measurement 334 328 Vtest=ZItemp(:,:,2)-ZItemp(:,:,1); 335 329 … … 337 331 [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b 338 332 339 333 340 334 for iview=1:2 341 %% reading input file(s)342 [Data{iview},tild,errormsg]=read_civdata(filecell{iview, index},{'vec(U,V)'},'*');335 %% reading PIV input file(s) 336 [Data{iview},tild,errormsg]=read_civdata(filecell{iview,field_index},{'vec(U,V)'},'*'); 343 337 if ~isempty(errormsg) 344 338 disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun) … … 359 353 end 360 354 end 361 %remove wrong vector355 %remove false vectors 362 356 temp=find(Data{1}.FF==0); 363 357 X1=Data{1}.X(temp); … … 369 363 Va=griddata(X1,Y1,V1,Xa,Ya); 370 364 371 % [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 365 [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 372 366 [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ 373 367 374 %remove wrong vector368 %remove false vectors 375 369 temp=find(Data{2}.FF==0); 376 370 X2=Data{2}.X(temp); … … 380 374 Ub=griddata(X2,Y2,U2,Xb,Yb); 381 375 Vb=griddata(X2,Y2,V2,Xb,Yb); 382 383 % [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 376 377 [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 384 378 [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ 385 379 386 380 387 381 % System to solve 388 382 S=ones(size(XI,1),size(XI,2),3); 389 383 D=ones(size(XI,1),size(XI,2),3,3); 390 384 391 385 S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; 392 386 S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; … … 408 402 W(indj,indi)=dxyz(3); 409 403 end 410 end 404 end 411 405 Error=zeros(size(XI,1),size(XI,2),4); 412 406 Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua; … … 415 409 Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb; 416 410 417 418 419 420 421 422 411 %% recording the merged field 423 if index==1% initiate the structure at first index412 if field_index==1% initiate the structure at first index 424 413 MergeData.ListGlobalAttribute={'Conventions','Time','Dt'}; 425 414 MergeData.Conventions='uvmat'; … … 428 417 MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'}; 429 418 MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}... 430 419 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; 431 420 MergeData.coord_x=xI; 432 421 MergeData.coord_y=yI; … … 437 426 MergeData.Z=ZI; 438 427 439 % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;440 % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;428 % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 429 % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 441 430 MergeData.Error=0.5*sqrt(sum(Error.^2,3)); 442 431 errormsg=struct2nc(OutputFile,MergeData);%save result file … … 461 450 A(:,:,2,3)=(R(6)-R(9)*Y)./T; 462 451 463 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U 464 452 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) 453 % convert image coordinates to view angles, after removal of quadratic distorsion 454 % input in pixel, output in radians 465 455 X1d=Xd-Ud/2; 466 456 X2d=Xd+Ud/2; … … 483 473 z=0; 484 474 error=0; 485 486 475 487 476 %% first image -
trunk/src/set_object.m
r924 r954 111 111 if isfield(data,'ProjModeMenu') 112 112 set(handles.ProjMode,'UserData',data.ProjModeMenu)% data.ProjModeMenu as default menu (used in Type_Callback) 113 end 113 end 114 114 errormsg=fill_GUI(data,handles.set_object); 115 115 if ~isempty(errormsg) … … 132 132 end 133 133 end 134 if isfield(data,'RangeX') 134 if isfield(data,'RangeX')&& ~strcmp(data.Type,'plane_z')%TODO: generalise 135 135 if ischar(data.RangeX) 136 136 data.RangeX=str2num(data.RangeX); … … 139 139 set(handles.num_RangeX_1,'String',num2str(min(data.RangeX),3)) 140 140 end 141 if isfield(data,'RangeY') 141 if isfield(data,'RangeY')&& ~strcmp(data.Type,'plane_z')%TODO: generalise 142 142 if ischar(data.RangeY) 143 143 data.RangeY=str2num(data.RangeY); … … 146 146 set(handles.num_RangeY_1,'String',num2str(min(data.RangeY),3)) 147 147 end 148 if isfield(data,'RangeZ') 148 if isfield(data,'RangeZ')&& ~strcmp(data.Type,'plane_z')%TODO: generalise 149 149 if ischar(data.RangeZ) 150 150 data.RangeZ=str2num(data.RangeZ); … … 247 247 if isempty(get(handles.ProjMode,'UserData')) 248 248 switch Type 249 case {'points','line','plane'}250 menu_proj={'projection';'interp_lin';'interp_tps';'none'};251 249 case 'polyline' 252 250 menu_proj={'interp_lin';'interp_tps';'none'}; … … 328 326 set(handles.num_RangeX_2,'TooltipString',['num_RangeX_2: half length of the ' ObjectStyle]) 329 327 set(handles.num_RangeY_2,'TooltipString',['num_RangeY_2: half width of the ' ObjectStyle]) 330 case {'plane' }328 case {'plane','plane_z'} 331 329 set(handles.num_Angle_3,'Visible','on') 332 330 set(handles.num_RangeX_1,'Visible','on') … … 387 385 set(handles.num_RangeY_2,'String',num2str(UvData.Field.YMax)) 388 386 end 389 if isempty(get(handles.CoordUnit,'String')) 390 set(handles.CoordUnit,'String', Field.CoordUnit)387 if isempty(get(handles.CoordUnit,'String'))&& isfield(UvData.Field,'CoordUnit') 388 set(handles.CoordUnit,'String',UvData.Field.CoordUnit) 391 389 end 392 390 end -
trunk/src/update_imadoc.m
r924 r954 34 34 testappend=1; 35 35 backupfile=outputfile; 36 testexist=2;37 while testexist==238 backupfile=[backupfile '~'];39 testexist=exist(backupfile,'file');40 end41 [success,message]=copyfile(outputfile,backupfile);%make backup42 if success~=143 errormsg=['errror in xml file backup: ' message];44 return45 end46 36 t=xmltree(outputfile); %read the file 47 37 title=get(t,1,'name'); 48 38 if strcmp(title,'ImaDoc') 49 testappend=1; 39 % testappend=1; 40 %rename the existing file for backup 41 testexist=2; 42 while testexist==2 43 backupfile=[backupfile '~']; 44 testexist=exist(backupfile,'file'); 45 end 46 [success,message]=movefile(outputfile,backupfile);%make backup 47 if success~=1 48 errormsg=['errror in xml file backup: ' message]; 49 return 50 end 50 51 %if the xml file is ImaDoc 51 52 uid_calib=find(t,['ImaDoc/' StructName]); -
trunk/src/uvmat.m
r951 r954 910 910 create_object(data,handles) 911 911 912 % -------------------------------------------------------------------- 913 function Menuplane_z_Callback(hObject, eventdata, handles) 914 data.Type='plane_z'; 915 data.ProjMode='projection';%default 916 data.ProjModeMenu={};% do not restrict ProjMode menus 917 create_object(data,handles) 918 912 919 %------------------------------------------------------------------------ 913 920 function Menuvolume_Callback(hObject, eventdata, handles) … … 949 956 check_plot=0; 950 957 if isfield(UvData,'Field') 951 Field=UvData.Field; 952 if isfield(Field,'NbDim')&& isequal(Field.NbDim,3) 958 if isfield(UvData.Field,'NbDim')&& isequal(UvData.Field.NbDim,3) 953 959 data.Coord=[0 0 0]; %default 954 960 end 955 if isfield( Field,'CoordUnit')956 data.CoordUnit= Field.CoordUnit;957 end 958 if isfield(UvData.Field,'CoordMesh')&&~isempty(UvData.Field.CoordMesh) 961 if isfield(UvData.Field,'CoordUnit') 962 data.CoordUnit=UvData.Field.CoordUnit; 963 end 964 if isfield(UvData.Field,'CoordMesh')&&~isempty(UvData.Field.CoordMesh)&&~strcmp(data.Type,'plane_z') 959 965 data.RangeX=[UvData.Field.XMin UvData.Field.XMax]; 960 966 switch data.Type … … 984 990 data.DX=UvData.Field.CoordMesh; 985 991 data.DY=UvData.Field.CoordMesh; 992 end 993 if isfield(UvData.Field,'ProjModeRequest') 994 data.ProjMode=UvData.Field.ProjModeRequest;%set the request proj mode option by default 986 995 end 987 996 end … … 5882 5891 set(handles.TableDisplay,'Visible','off') 5883 5892 end 5893 5894 5895
Note: See TracChangeset
for help on using the changeset viewer.
