Index: /trunk/src/find_field_bounds.m =================================================================== --- /trunk/src/find_field_bounds.m (revision 953) +++ /trunk/src/find_field_bounds.m (revision 954) @@ -54,4 +54,5 @@ CoordMin=zeros(numel(imax),NbDim); Mesh=zeros(1,numel(imax)); +FieldOut.ProjModeRequest='projection';%default for ind=1:numel(imax) if strcmp(CellInfo{imax(ind)}.CoordType,'tps') @@ -90,4 +91,12 @@ Mesh(ind)=min((CoordMax(ind,:)-CoordMin(ind,:))./(NbPoints-1)); end + if isfield(CellInfo{imax(ind)},'ProjModeRequest') + if strcmp(CellInfo{imax(ind)}.ProjModeRequest,'interp_tps') + FieldOut.ProjModeRequest='interp_tps'; + end + if strcmp(CellInfo{imax(ind)}.ProjModeRequest,'interp_lin')&& ~strcmp(FieldOut.ProjModeRequest,'interp_tps') + FieldOut.ProjModeRequest='interp_lin'; + end + end end Mesh=min(Mesh); Index: /trunk/src/mouse_motion.m =================================================================== --- /trunk/src/mouse_motion.m (revision 953) +++ /trunk/src/mouse_motion.m (revision 954) @@ -385,21 +385,22 @@ XYData=AxeData.CurrentOrigin; if isequal(AxeData.Drawing,'create') && isfield(AxeData,'CurrentOrigin') && ~isempty(AxeData.CurrentOrigin) - if strcmp(ObjectData.Type,'line')||strcmp(ObjectData.Type,'polyline')||strcmp(ObjectData.Type,'polygon')||strcmp(ObjectData.Type,'points') - ObjectData.Coord=[ObjectData.Coord ;xy(1,1:2)]; - % ObjectData.Coord(end,:)=xy(1,:); - elseif strcmp(ObjectData.Type,'rectangle')||strcmp(ObjectData.Type,'ellipse')||strcmp(ObjectData.Type,'volume') - ObjectData.Coord=(AxeData.CurrentOrigin+xy(1,1:2))/2;% keep only the first point coordinate + switch ObjectData.Type + case {'line','polyline','polygon','points','plane_z'} + ObjectData.Coord=[ObjectData.Coord ;xy(1,1:2)]; + % ObjectData.Coord(end,:)=xy(1,:); + case {'rectangle','ellipse','volume'} + ObjectData.Coord=(AxeData.CurrentOrigin+xy(1,1:2))/2;% keep only the first point coordinate ObjectData.RangeX=abs(ObjectData.Coord(1,1)-xy(1,1));%rectangle width - ObjectData.RangeY=abs(ObjectData.Coord(1,2)-xy(1,2));%rectangle height - elseif isequal(ObjectData.Type,'plane') %case of 'plane' - DX=(xy(1,1)-ObjectData.Coord(1,1)); - DY=(xy(1,2)-ObjectData.Coord(1,2)); - ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle widt - if isfield(ObjectData,'RangeX') - XMax=sqrt(DX*DX+DY*DY); - if XMax>max(ObjectData.RangeX) - ObjectData.RangeX=[min(ObjectData.RangeX) XMax]; - end - end + ObjectData.RangeY=abs(ObjectData.Coord(1,2)-xy(1,2));%rectangle height + case 'plane' %case of 'plane' + DX=(xy(1,1)-ObjectData.Coord(1,1)); + DY=(xy(1,2)-ObjectData.Coord(1,2)); + ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle widt + if isfield(ObjectData,'RangeX') + XMax=sqrt(DX*DX+DY*DY); + if XMax>max(ObjectData.RangeX) + ObjectData.RangeX=[min(ObjectData.RangeX) XMax]; + end + end end plot_object(ObjectData,ProjObject,AxeData.CurrentObject,'m'); Index: /trunk/src/mouse_up.m =================================================================== --- /trunk/src/mouse_up.m (revision 953) +++ /trunk/src/mouse_up.m (revision 954) @@ -95,5 +95,5 @@ else switch ObjectData.Type - case {'line'} + case {'line','plane_z'} if size(ObjectData.Coord,1)==1 % this is the mouse up for the first point, continue until next click check_multiple=1; @@ -106,5 +106,5 @@ check_multiple=1;% pass to next mous up if width of height=0 end - case 'plane' %case of 'plane' + case 'plane' %case of 'plane', TODO: NOT ACTIVATED DX=(xy(1,1)-ObjectData.Coord(1,1)); DY=(xy(1,2)-ObjectData.Coord(1,2)); @@ -116,4 +116,11 @@ end end + case 'plane_z' + if size(ObjectData.Coord,1)==1 % this is the mouse up for the first point, continue until next click + check_multiple=1; + end + DX=(xy(1,1)-ObjectData.Coord(1,1)); + DY=(xy(1,2)-ObjectData.Coord(1,2)); + ObjectData.Phi=(angle(DX+i*DY))*180/pi;%rectangle width otherwise check_multiple=1; Index: /trunk/src/plot_object.m =================================================================== --- /trunk/src/plot_object.m (revision 953) +++ /trunk/src/plot_object.m (revision 954) @@ -133,131 +133,134 @@ %% determine the coordinates xline, yline,xsup,xinf, yinf,ysup determining the new object plot -test_line= isequal(ObjectData.Type,'points')||isequal(ObjectData.Type,'line')||... - isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon')|| isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume'); -test_patch=isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')||isequal(ObjectData.Type,'volume')... - ||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'); +%test_line= isequal(ObjectData.Type,'points')||isequal(ObjectData.Type,'line')||... +% isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon')|| isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume'); +%test_patch=isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')||isequal(ObjectData.Type,'volume')... +% ||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'); +test_line=ismember(ObjectData.Type,{'points','line','polyline','polygon','plane','plane_z','volume'}); +test_patch=ismember(ObjectData.ProjMode,{'inside','outside','mask_inside','mask_outside'}); if test_line xline=ObjectData.Coord(:,1); yline=ObjectData.Coord(:,2); nbpoints=numel(xline); - if isequal(ObjectData.Type,'line_x') - xline=[xline; ObjectData.RangeX(2)];%creating the line - yline=[yline; ObjectData.RangeY(2)];%creating the line - elseif isequal(ObjectData.Type,'polygon') - xline=[xline; ObjectData.Coord(1,1)];%closing the line - yline=[yline; ObjectData.Coord(1,2)]; - elseif isequal(ObjectData.Type,'plane')|| isequal(ObjectData.Type,'volume') - if ~isfield(ObjectData,'Angle') - ObjectData.Angle=[0 0 0]; - end - phi=ObjectData.Angle(3)*pi/180;%angle in radians - x0=xline(1); y0=yline(1); - xlim=get(haxes,'XLim'); - ylim=get(haxes,'YLim'); - graph_scale=max(abs(xlim(2)-xlim(1)),abs(ylim(2)-ylim(1)))/2;% estimate the length of axes plots - XMax=graph_scale; - YMax=graph_scale; - XMin=-graph_scale; - YMin=-graph_scale; - if ~isempty(XMaxRange) - XMax=XMaxRange; - end - if ~isempty(XMinRange) - XMin=XMinRange; - end - if ~isempty(YMaxRange) - YMax=YMaxRange; - end - if ~isempty(YMinRange) - YMin=YMinRange; - end - % axes lines - xline=NaN(1,13); - xline(1)=x0+min(0,XMin)*cos(phi); % min end of the x axes - yline(1)=y0+min(0,XMin)*sin(phi); - xline(2)=x0+XMax*cos(phi);% max end of the x axes - yline(2)=y0+XMax*sin(phi); - xline(8)=x0-min(0,YMin)*sin(phi);% min end of the y axes - yline(8)=y0+min(0,YMin)*cos(phi); - xline(9)=x0-YMax*sin(phi);% max end of the y axes - yline(9)=y0+YMax*cos(phi); - - %arrows on x axis - arrow_scale=graph_scale/20; - xline(3)=xline(2)-arrow_scale*cos(phi-pi/8); - yline(3)=yline(2)-arrow_scale*sin(phi-pi/8); - xline(5)=xline(2); - yline(5)=yline(2); - xline(6)=xline(2)-arrow_scale*cos(phi+pi/8); - yline(6)=yline(2)-arrow_scale*sin(phi+pi/8); - - %arrows on y axis - xline(10)=xline(9)-arrow_scale*cos(phi+pi/2-pi/8); - yline(10)=yline(9)-arrow_scale*sin(phi+pi/2-pi/8); - xline(12)=xline(9); - yline(12)=yline(9); - xline(13)=xline(9)-arrow_scale*cos(phi+pi/2+pi/8); - yline(13)=yline(9)-arrow_scale*sin(phi+pi/2+pi/8); - %xline=[Xbeg_x Xend_x NaN Ybeg_x Yend_x]; - %yline=[Xbeg_y Xend_y NaN Ybeg_y Yend_y]; - % dashed lines indicating bounds - xsup=NaN(1,5); - ysup=NaN(1,5); - if ~isempty(XMaxRange) - xsup(1)=xline(2)-YMin*sin(phi); - ysup(1)=yline(2)+YMin*cos(phi); - xsup(2)=xline(2)-YMax*sin(phi); - ysup(2)=yline(2)+YMax*cos(phi); - end - if ~isempty(YMaxRange) - xsup(2)=xline(2)-YMax*sin(phi); - ysup(2)=yline(2)+YMax*cos(phi); - xsup(3)=xline(9)+XMin*cos(phi); - ysup(3)=yline(9)+XMin*sin(phi); - end - if ~isempty(XMinRange) - xsup(3)=xline(9)+XMin*cos(phi); - ysup(3)=yline(9)+XMin*sin(phi); - xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); - ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); - end - if ~isempty(YMinRange) - xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); - ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); - xsup(5)=xline(8)-YMin*sin(phi); - ysup(5)=yline(8)+YMin*cos(phi); - end - end - SubLineStyle='none';%default - if isfield(ObjectData,'ProjMode') - if isequal(ObjectData.ProjMode,'projection') - SubLineStyle='--'; %range of projection marked by dash - if isfield (ObjectData,'DX') - ObjectData=rmfield(ObjectData,'DX'); - end - if isfield (ObjectData,'DY') - ObjectData=rmfield(ObjectData,'DY'); - end - elseif isequal(ObjectData.ProjMode,'filter') - SubLineStyle=':';%range of projection not visible - end - end - if isequal(ObjectData.Type,'line')||isequal(ObjectData.Type,'polyline')||isequal(ObjectData.Type,'polygon') - if length(xline)<2 - theta=0; - else - theta=angle(diff(xline)+1i*diff(yline)); - theta(length(xline))=theta(length(xline)-1); - end - xsup(1)=xline(1)+YMax*sin(theta(1)); - xinf(1)=xline(1)-YMax*sin(theta(1)); - ysup(1)=yline(1)-YMax*cos(theta(1)); - yinf(1)=yline(1)+YMax*cos(theta(1)); - for ip=2:length(xline) - xsup(ip)=xline(ip)+YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); - xinf(ip)=xline(ip)-YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); - ysup(ip)=yline(ip)-YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); - yinf(ip)=yline(ip)+YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); - end + switch ObjectData.Type + case 'line_x' + xline=[xline; ObjectData.RangeX(2)];%creating the line + yline=[yline; ObjectData.RangeY(2)];%creating the line + case 'polygon' + xline=[xline; ObjectData.Coord(1,1)];%closing the line + yline=[yline; ObjectData.Coord(1,2)]; + case {'plane','volume'} + if ~isfield(ObjectData,'Angle') + ObjectData.Angle=[0 0 0]; + end + phi=ObjectData.Angle(3)*pi/180;%angle in radians + x0=xline(1); y0=yline(1); + xlim=get(haxes,'XLim'); + ylim=get(haxes,'YLim'); + graph_scale=max(abs(xlim(2)-xlim(1)),abs(ylim(2)-ylim(1)))/2;% estimate the length of axes plots + XMax=graph_scale; + YMax=graph_scale; + XMin=-graph_scale; + YMin=-graph_scale; + if ~isempty(XMaxRange) + XMax=XMaxRange; + end + if ~isempty(XMinRange) + XMin=XMinRange; + end + if ~isempty(YMaxRange) + YMax=YMaxRange; + end + if ~isempty(YMinRange) + YMin=YMinRange; + end + % axes lines + xline=NaN(1,13); + xline(1)=x0+min(0,XMin)*cos(phi); % min end of the x axes + yline(1)=y0+min(0,XMin)*sin(phi); + xline(2)=x0+XMax*cos(phi);% max end of the x axes + yline(2)=y0+XMax*sin(phi); + xline(8)=x0-min(0,YMin)*sin(phi);% min end of the y axes + yline(8)=y0+min(0,YMin)*cos(phi); + xline(9)=x0-YMax*sin(phi);% max end of the y axes + yline(9)=y0+YMax*cos(phi); + + %arrows on x axis + arrow_scale=graph_scale/20; + xline(3)=xline(2)-arrow_scale*cos(phi-pi/8); + yline(3)=yline(2)-arrow_scale*sin(phi-pi/8); + xline(5)=xline(2); + yline(5)=yline(2); + xline(6)=xline(2)-arrow_scale*cos(phi+pi/8); + yline(6)=yline(2)-arrow_scale*sin(phi+pi/8); + + %arrows on y axis + xline(10)=xline(9)-arrow_scale*cos(phi+pi/2-pi/8); + yline(10)=yline(9)-arrow_scale*sin(phi+pi/2-pi/8); + xline(12)=xline(9); + yline(12)=yline(9); + xline(13)=xline(9)-arrow_scale*cos(phi+pi/2+pi/8); + yline(13)=yline(9)-arrow_scale*sin(phi+pi/2+pi/8); + %xline=[Xbeg_x Xend_x NaN Ybeg_x Yend_x]; + %yline=[Xbeg_y Xend_y NaN Ybeg_y Yend_y]; + % dashed lines indicating bounds + xsup=NaN(1,5); + ysup=NaN(1,5); + if ~isempty(XMaxRange) + xsup(1)=xline(2)-YMin*sin(phi); + ysup(1)=yline(2)+YMin*cos(phi); + xsup(2)=xline(2)-YMax*sin(phi); + ysup(2)=yline(2)+YMax*cos(phi); + end + if ~isempty(YMaxRange) + xsup(2)=xline(2)-YMax*sin(phi); + ysup(2)=yline(2)+YMax*cos(phi); + xsup(3)=xline(9)+XMin*cos(phi); + ysup(3)=yline(9)+XMin*sin(phi); + end + if ~isempty(XMinRange) + xsup(3)=xline(9)+XMin*cos(phi); + ysup(3)=yline(9)+XMin*sin(phi); + xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); + ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); + end + if ~isempty(YMinRange) + xsup(4)=x0+XMin*cos(phi)-YMin*sin(phi); + ysup(4)=y0+XMin*sin(phi)+YMin*cos(phi); + xsup(5)=xline(8)-YMin*sin(phi); + ysup(5)=yline(8)+YMin*cos(phi); + end + end +end +SubLineStyle='none';%default +if isfield(ObjectData,'ProjMode') + if isequal(ObjectData.ProjMode,'projection') + SubLineStyle='--'; %range of projection marked by dash + if isfield (ObjectData,'DX') + ObjectData=rmfield(ObjectData,'DX'); + end + if isfield (ObjectData,'DY') + ObjectData=rmfield(ObjectData,'DY'); + end + elseif isequal(ObjectData.ProjMode,'filter') + SubLineStyle=':';%range of projection not visible + end +end +if ismember(ObjectData.Type,{'line','polyline','polygon','plane_z'}) + if length(xline)<2 + theta=0; + else + theta=angle(diff(xline)+1i*diff(yline)); + theta(length(xline))=theta(length(xline)-1); + end + xsup(1)=xline(1)+YMax*sin(theta(1)); + xinf(1)=xline(1)-YMax*sin(theta(1)); + ysup(1)=yline(1)-YMax*cos(theta(1)); + yinf(1)=yline(1)+YMax*cos(theta(1)); + for ip=2:length(xline) + xsup(ip)=xline(ip)+YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); + xinf(ip)=xline(ip)-YMax*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); + ysup(ip)=yline(ip)-YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); + yinf(ip)=yline(ip)+YMax*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2); end end @@ -415,9 +418,9 @@ else% no patch image requested, erase existing ones if isfield(PlotData,'SubObject') - for iobj=1:length(PlotData.SubObject) - if ishandle(PlotData.SubObject(iobj)) && strcmp(get(PlotData.SubObject(iobj),'Type'),'image') - delete(PlotData.SubObject(iobj)) - end - end + for iobj=1:length(PlotData.SubObject) + if ishandle(PlotData.SubObject(iobj)) && strcmp(get(PlotData.SubObject(iobj),'Type'),'image') + delete(PlotData.SubObject(iobj)) + end + end end end @@ -459,5 +462,5 @@ end end - case {'line','polyline','polygon'} + case {'line','polyline','polygon','plane_z'} hh=line(xline,yline,'Color',col); PlotData.SubObject(1)=line(xinf,yinf,'Color',col,'LineStyle',SubLineStyle,'Tag','proj_object');%draw sub-lines Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 953) +++ /trunk/src/proj_field.m (revision 954) @@ -90,6 +90,6 @@ return end -ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes -if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case +% check list of effective projection modes +if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'}) return end @@ -117,5 +117,5 @@ [ProjData,errormsg] = proj_line(FieldData,ObjectData); end - case 'plane' + case {'plane','plane_z'} [ProjData,errormsg] = proj_plane(FieldData,ObjectData); case 'volume' @@ -955,4 +955,12 @@ test90x=0;%=1 for 90 degree rotation alround x axis test90y=0;%=1 for 90 degree rotation alround y axis +if strcmp(ObjectData.Type,'plane_z') + Delta_x=ObjectData.Coord(2,1)-ObjectData.Coord(1,1); + Delta_y=ObjectData.Coord(2,2)-ObjectData.Coord(1,2); + Delta_mod=sqrt(Delta_x*Delta_x+Delta_y*Delta_y); + ObjectData.Angle=[0 0 0]; + ObjectData.Angle(1)=90*Delta_x/Delta_mod; + ObjectData.Angle(2)=90*Delta_y/Delta_mod; +end if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) test90y=isequal(ObjectData.Angle,[0 90 0]); @@ -987,4 +995,8 @@ return end +InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane +if strcmp(ObjectData.Type,'plane_z') + InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation +end %% extrema along each axis @@ -1071,7 +1083,7 @@ end icell_grid=[];% field cell index which defines the grid -if ~strcmp(ObjectData.ProjMode,'projection') +if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize %% define the new coordinates in case of interpolation on a imposed grid - if ~testYMin + if ~testYMin errormsg='min Y value not defined for the projection grid';return end @@ -1119,5 +1131,5 @@ AYName='coord_y'; AXName='coord_x'; - if strcmp(ObjectData.ProjMode,'projection') + if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z') 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 ProjData.coord_x=[FieldData.XMin FieldData.XMax]; @@ -1399,5 +1411,5 @@ DimValue(DimValue==1)=[];%remove singleton dimensions NbDim=numel(DimValue);%update number of space dimensions - nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate + % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate if NbDim>=3 if NbDim>3 @@ -1439,8 +1451,8 @@ Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)}); end - if numel(Coord{NbDim-1})==2 + if numel(Coord{NbDim-1})==2% case of coordinate defined only by the first and last values DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1); end - if numel(Coord{NbDim})==2 + if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1); end @@ -1453,9 +1465,9 @@ end else - YIndexMax=Coord{NbDim-1}(end)/DY; + YIndexMax=numel(Coord{NbDim-1}); YIndexMin=1; end if testXMax - XIndexMax=(XMax-Coord{NbDim}(1))/DY+1;% matrix index corresponding to the max y value for the new field + XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field if testYMin%test_direct(indY) XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field @@ -1464,5 +1476,5 @@ end else - XIndexMax=Coord{NbDim}(end)/DX; + XIndexMax=numel(Coord{NbDim}); XIndexMin=1; end @@ -1627,8 +1639,49 @@ end end - else - errormsg='projection of structured coordinates on oblique plane not yet implemented'; - %TODO: use interp_lin3 - return + else %projection of structured coordinates on oblique plane + % determine the boundaries of the projected field, + % first find the 8 summits of the initial volume in the + % new coordinates + Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates + Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates + Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates + summit=zeros(3,8);% initialize summit coordinates + summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square + summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1) + summit(1:2,5:8)=summit(1:2,1:4); + summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)]; + Mrot=rodrigues(PlaneAngle); + newsummit=zeros(3,8);% initialize the rotated summit coordinates + ObjectData.Coord=[ObjectData.Coord(1,1); ObjectData.Coord(1,2); 0];%TODO: set in set_object + for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data) + newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord)); + end + coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:)); + coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:)); + coord_z_proj=-width:InterpMesh:width; + Mrot_inverse=rodrigues(-PlaneAngle); + Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord; + ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0];% unit vector along the new x coordinates + iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0];% unit vector y coordinates + iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)];% x unit vector z coordinates + [Grid_x,Grid_y,Grid_z]=meshgrid(1:numel(coord_x_proj),1:numel(coord_y_proj),1:numel(coord_z_proj)); + if ismatrix(Grid_x) + Grid_x=shiftdim(Grid_x,-1); + Grid_y=shiftdim(Grid_y,-1); + Grid_z=shiftdim(Grid_z,-1); + end + XI=Origin(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z; + YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z; + ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z; + [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1}); + X=permute(X,[3 1 2]); + Y=permute(Y,[3 1 2]); + Z=permute(Z,[3 1 2]); + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + ListVarName=[ListVarName VarName]; + VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes + ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI); + end end end @@ -2264,4 +2317,6 @@ end else + RotMatrix=rodrigues(om); + errormsg='projection of structured coordinates on oblique plane not yet implemented'; %TODO: use interp3 Index: /trunk/src/series.m =================================================================== --- /trunk/src/series.m (revision 953) +++ /trunk/src/series.m (revision 954) @@ -2292,5 +2292,5 @@ CheckList_1=1;% indicate whether FieldName_1 has been updated handles_coord=[handles.Coord_x handles.Coord_y handles.Coord_z handles.Coord_x_title handles.Coord_y_title handles.Coord_z_title]; -if VelTypeRequest && numel(iview_civ)>=1 +if VelTypeRequest && numel(iview_civ)>=1 menu=set_veltype_display(SeriesData.FileInfo{iview_civ(1)}.CivStage,SeriesData.FileType{iview_civ(1)}); set(handles.VelType,'Value',1)% set first choice by default @@ -2301,5 +2301,5 @@ %CheckList=1; set(handles.FieldName,'Value',1); %velocity vector choice by default - if VelTypeRequest_1 && numel(iview_civ)>=2 + if VelTypeRequest_1 && numel(iview_civ)>=2 menu=set_veltype_display(SeriesData.FileInfo{iview_civ(2)}.CivStage,SeriesData.FileType{iview_civ(2)}); set(handles.VelType_1,'Value',1)% set first choice by default @@ -2317,9 +2317,9 @@ set(handles.VelType,'Visible','off') set(handles.VelType_title,'Visible','off') -end +end %% Detect the types of input files and set menus and default options in 'FieldName' if (FieldNameRequest || VelTypeRequest) && numel(iview_netcdf)>=1 - set(handles.InputFields,'Visible','on') + set(handles.InputFields,'Visible','on')% set the frame InputFields visible if FieldNameRequest && isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName') set(handles.FieldName,'Visible','on') @@ -2346,42 +2346,46 @@ end end - set(handles_coord,'Visible','on') - FieldList=[FieldList;{'get_field...'}]; - if FieldNameRequest_1 && numel(iview_netcdf)>=2 + end + + set(handles_coord,'Visible','on') + FieldList=[FieldList;{'get_field...'}]; + if FieldNameRequest_1 && numel(iview_netcdf)>=2 + set(handles.FieldName_1,'Visible','on') + if CheckList_1==0 % not civ input made + ListVarName=SeriesData.FileInfo{iview_netcdf(2)}.ListVarName; + ind_var=get(handles.FieldName,'Value');%indices of previously selected variables + for ilist=1:numel(ind_var) + if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName))) + FieldList_1={};% previous choice not consistent with new input field + set(handles.FieldName_1,'Value',1) + break + end + end + warn_coord=0; + if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||... + isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName))) + warn_coord=1; + end + if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName))) + FieldList_1={}; + warn_coord=1; + end + if warn_coord + msgbox_uvmat('WARNING','coordiante names do not exist in the second netcdf input file') + end + set(handles.FieldName_1,'Visible','on') - if CheckList_1==0 % not civ input made - ListVarName=SeriesData.FileInfo{iview_netcdf(2)}.ListVarName; - ind_var=get(handles.FieldName,'Value');%indices of previously selected variables - for ilist=1:numel(ind_var) - if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName))) - FieldList_1={};% previous choice not consistent with new input field - set(handles.FieldName_1,'Value',1) - break - end - end - warn_coord=0; - if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||... - isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName))) - warn_coord=1; - end - if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName))) - FieldList_1={}; - warn_coord=1; - end - if warn_coord - msgbox_uvmat('WARNING','coordiante names do not exist in the second netcdf input file') - end - set(handles.FieldName,'String',[FieldList;{'get_field...'}]) - set(handles.FieldName_1,'Visible','on') - set(handles.FieldName_1,'Value',1) - set(handles.FieldName_1,'String',FieldList_1) - end - else - set(handles.FieldName_1,'Visible','off') + set(handles.FieldName_1,'Value',1) + set(handles.FieldName_1,'String',FieldList_1) end else + set(handles.FieldName_1,'Visible','off') + end + if isempty(FieldList) set(handles.FieldName,'Visible','off') - end - set(handles.FieldName,'String',FieldList) + else + set(handles.FieldName,'Visible','on') + set(handles.FieldName,'String',FieldList) + end else set(handles.InputFields,'Visible','off') Index: /trunk/src/series/civ2vel_3C.m =================================================================== --- /trunk/src/series/civ2vel_3C.m (revision 953) +++ /trunk/src/series/civ2vel_3C.m (revision 954) @@ -1,3 +1,15 @@ -%'civ2vel_3C': combine velocity fields from two cameras to get three velocity components +%'civ2vel_3C': combine velocity fields from two cameras to get the three velocity components +% used with the GUI 'series': +% first input line =raw PIV camera 1 (image coordinates) +% second input line=raw PIV camera 2 (image coordinates) +% Three modes: +% 1) no additional input: measurements assumed in the reference plane (laser sheet) +% 2) measurement surface obtained by stereoscopic comparison of the images of the two cameras. +% third input line =surface z(x,y) given by correlation between camera 1 and 2 (expressed in phys apparent coordinates) +% 3) surface z(x,y) given by adding the displacements of each camera with a third intermediate camera (#3) used to reduce the +% to reduce the angle for stereoscopic view. +% third input line =correlation between camera 1 and 3 (expressed in phys apparent coordinates) +% fourth input line =corelation between camera 2 and 3 (expressed in phys apparent coordinates) +% %------------------------------------------------------------------------ % function ParamOut=civ2vel_3C(Param) @@ -7,4 +19,5 @@ % %INPUT: +% % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. % In batch mode, Param is the name of the corresponding xml file containing the same information @@ -15,5 +28,6 @@ % see the current structure Param) % .InputTable: cell of input file names, (several lines for multiple input) -% each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} +% each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}. +% 3 or 4 lines used as described above % .OutputSubDir: name of the subdirectory for data outputs % .OutputDirExt: directory extension for data outputs @@ -22,6 +36,5 @@ % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct % .RUN =0 for GUI input, =1 for function activation -% .RunMode='local','background', 'cluster': type of function use -% +% .RunMode='local','background', 'cluster': type of function use % .IndexRange: set the file or frame indices on which the action must be performed % .InputFields: sub structure describing the input fields withfields @@ -52,5 +65,5 @@ function ParamOut=civ2vel_3C(Param) -disp('test') + %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed if isstruct(Param) && isequal(Param.Action.RUN,0) @@ -124,8 +137,4 @@ %% define the directory for result file (with path=RootPath{1}) OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files -% -% if ~isfield(Param,'InputFields') -% Param.InputFields.FieldName=''; -% end %% calibration data and timing: read the ImaDoc files @@ -160,5 +169,5 @@ return end -ObjectData=Param.ProjObject; +ObjectData=Param.ProjObject;% Object attached to the GUI series xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); @@ -168,5 +177,5 @@ W=zeros(size(XI,1),size(XI,2)); -%% MAIN LOOP ON FIELDS +%%%% MAIN LOOP ON FIELDS %%%%% warning off @@ -175,15 +184,12 @@ CheckOverwrite=Param.CheckOverwrite; end -for index=1:NbField - - update_waitbar(WaitbarHandle,index/NbField) - - - - - %% generating the name of the merged field - i1=i1_series{1}(index); +for field_index=1:NbField + + update_waitbar(WaitbarHandle,field_index/NbField)% waitbar to visualise progress (in RUN mode) + + %% generating the name of the output file for the merged field + i1=i1_series{1}(field_index); if ~isempty(i2_series{end}) - i2=i2_series{end}(index); + i2=i2_series{end}(field_index); else i2=i1; @@ -192,7 +198,7 @@ j2=1; if ~isempty(j1_series{1}) - j1=j1_series{1}(index); + j1=j1_series{1}(field_index); if ~isempty(j2_series{end}) - j2=j2_series{end}(index); + j2=j2_series{end}(field_index); else j2=j1; @@ -201,7 +207,5 @@ OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2); - %% - - + %% program stop if requested on the GUI if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') disp('program stopped by user') @@ -209,127 +213,117 @@ end - if (~CheckOverwrite && exist(OutputFile,'file')) - disp('existing output file already exists, skip to next field') - continue% skip iteration if the mode overwrite is desactivated and the result file already exists - end - + %% check existence of the output file + if (~CheckOverwrite && exist(OutputFile,'file')) + disp('existing output file already exists, skip to next field') + continue% skip iteration if the mode overwrite is desactivated and the result file already exists + end + %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% Data=cell(1,NbView);%initiate the set Data - timeread=zeros(1,NbView); %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken %at the middle between to time step - clear ZItemp - ZItemp=zeros(size(XI,1),size(XI,2),2); - - if index==1 + ZItemp=zeros(size(XI,1),size(XI,2),2); + if field_index==1 first_img=i1_series{1,1}(1,1); %id of the first image of the series - end - - idtemp=0; - for indextemp=index:index+1; - idtemp=idtemp+1; - if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys - - - - [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']); - - if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector - temp=find(Data{3}.Civ3_FF==0); - Zphys=Data{3}.Zphys(temp); - Yphys=Data{3}.Yphys(temp); - Xphys=Data{3}.Xphys(temp); - else - Zphys=Data{3}.Zphys; - Yphys=Data{3}.Yphys; - Xphys=Data{3}.Xphys; + end + + idtemp=0; + % get the surface shape corresponding to the PIV measurements + for indextemp=field_index:field_index+1;%TODO: generalise to field index intervals>1 for PIV + idtemp=idtemp+1; + InputFile_3=fullfile(Param.InputTable{3,1},Param.InputTable{3,2},[Param.InputTable{3,3} '_' int2str(first_img+indextemp-1) '.nc']); + if NbView==3 % if there is only 1 stereo folder (2 cameras only), extract directly Xphys,Yphys and Zphys + % Data{1}: =raw PIV camera 1 only + % Data{2}: =raw PIV camera 2 only + % Data{3}: =correlation between camera 1 and 2 + [Data{3},tild,errormsg] = nc2struct(InputFile_3);%read input file + if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector + temp=find(Data{3}.Civ3_FF==0); + Zphys=Data{3}.Zphys(temp); + Yphys=Data{3}.Yphys(temp); + Xphys=Data{3}.Xphys(temp); + else + Zphys=Data{3}.Zphys; + Yphys=Data{3}.Yphys; + Xphys=Data{3}.Xphys; + end + + elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys + % Data{1}: =raw PIV camera 1 only (left) + % Data{2}: =raw PIV camera 2 only (right) (no PIV done with middle camera) + % Data{3}: =corelation between camera 1 and 3 (left and middle) + % Data{4}: =corelation between camera 2 and 3 (right and middle) + % test if the second camera (3) is the same for both folders + for i=3:4 + indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 + indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1 + camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name + end + if strcmp(camname{3},camname{4})==0 + disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) + return + end + [Data{3},tild,errormsg] = nc2struct(InputFile_3); + if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector + temp=find(Data{3}.Civ3_FF==0); + Xmid3=Data{3}.Xmid(temp); + Ymid3=Data{3}.Ymid(temp); + U3=Data{3}.Uphys(temp); + V3=Data{3}.Vphys(temp); + else + Xmid3=Data{3}.Xmid; + Ymid3=Data{3}.Ymid; + U3=Data{3}.Uphys; + V3=Data{3}.Vphys; + end + %temporary grid of merging the 2 stereos data + [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)); + + %1st folder : interpolate the first camera points on the second (common) camera + %(Dalsa 3) + x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq); + y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq); + + InputFile_4=fullfile(Param.InputTable{4,1},Param.InputTable{4,2},[Param.InputTable{4,3} '_' int2str(first_img+indextemp-1) '.nc']); + [Data{4},tild,errormsg] = nc2struct(InputFile_4); + if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector + temp=find(Data{4}.Civ3_FF==0); + Xmid4=Data{4}.Xmid(temp); + Ymid4=Data{4}.Ymid(temp); + U4=Data{4}.Uphys(temp); + V4=Data{4}.Vphys(temp); + else + Xmid4=Data{4}.Xmid; + Ymid4=Data{4}.Ymid; + U4=Data{4}.Uphys; + V4=Data{4}.Vphys; + end + %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera + %(Dalsa 3) + x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq); + y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq); + + %add the displacements of the two camera pairs + xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); + ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); + u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1); + v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1); + + % get the surface z(x,y) from the combined displacement + [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys + %remove NaN + tempNaN=isnan(Zphys);tempind=find(tempNaN==1); + Zphys(tempind)=[]; + Xphys(tempind)=[]; + Yphys(tempind)=[]; + error(tempind)=[]; + end + ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid - - elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys - - - %test if the seconde camera is the same for both folder - for i=3:4 - indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 - indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1 - camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name - end - - if strcmp(camname{3},camname{4})==0 - disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) - return - end - - - - [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']); - - if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector - temp=find(Data{3}.Civ3_FF==0); - Xmid3=Data{3}.Xmid(temp); - Ymid3=Data{3}.Ymid(temp); - U3=Data{3}.Uphys(temp); - V3=Data{3}.Vphys(temp); - else - Xmid3=Data{3}.Xmid; - Ymid3=Data{3}.Ymid; - U3=Data{3}.Uphys; - V3=Data{3}.Vphys; - end - %temporary gridd of merging the 2 stereos datas - [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)); - - %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera - %(Dalsa 3) - x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq); - y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq); - - - - [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp-1),'.nc']); - if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector - temp=find(Data{4}.Civ3_FF==0); - Xmid4=Data{4}.Xmid(temp); - Ymid4=Data{4}.Ymid(temp); - U4=Data{4}.Uphys(temp); - V4=Data{4}.Vphys(temp); - else - Xmid4=Data{4}.Xmid; - Ymid4=Data{4}.Ymid; - U4=Data{4}.Uphys; - V4=Data{4}.Vphys; - end - - %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera - %(Dalsa 3) - x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq); - y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq); - - xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); - ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); - u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1); - v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1); - - - [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys - %remove NaN - tempNaN=isnan(Zphys);tempind=find(tempNaN==1); - Zphys(tempind)=[]; - Xphys(tempind)=[]; - Yphys(tempind)=[]; - error(tempind)=[]; - - end - - - - - ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd - -end - ZI=mean(ZItemp,3); %mean between two the two time step + end + ZI=mean(ZItemp,3); %mean between two the two times used for surface measurement Vtest=ZItemp(:,:,2)-ZItemp(:,:,1); @@ -337,8 +331,8 @@ [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b - + for iview=1:2 - %% reading input file(s) - [Data{iview},tild,errormsg]=read_civdata(filecell{iview,index},{'vec(U,V)'},'*'); + %% reading PIV input file(s) + [Data{iview},tild,errormsg]=read_civdata(filecell{iview,field_index},{'vec(U,V)'},'*'); if ~isempty(errormsg) disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun) @@ -359,5 +353,5 @@ end end - %remove wrong vector + %remove false vectors temp=find(Data{1}.FF==0); X1=Data{1}.X(temp); @@ -369,8 +363,8 @@ Va=griddata(X1,Y1,V1,Xa,Ya); -% [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X + [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ - %remove wrong vector + %remove false vectors temp=find(Data{2}.FF==0); X2=Data{2}.X(temp); @@ -380,13 +374,13 @@ Ub=griddata(X2,Y2,U2,Xb,Yb); Vb=griddata(X2,Y2,V2,Xb,Yb); - -% [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X + + [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ - + % System to solve S=ones(size(XI,1),size(XI,2),3); D=ones(size(XI,1),size(XI,2),3,3); - + S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; @@ -408,5 +402,5 @@ W(indj,indi)=dxyz(3); end - end + end Error=zeros(size(XI,1),size(XI,2),4); Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua; @@ -415,11 +409,6 @@ Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb; - - - - - %% recording the merged field - if index==1% initiate the structure at first index + if field_index==1% initiate the structure at first index MergeData.ListGlobalAttribute={'Conventions','Time','Dt'}; MergeData.Conventions='uvmat'; @@ -428,5 +417,5 @@ MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'}; MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}... - {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; + {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; MergeData.coord_x=xI; MergeData.coord_y=yI; @@ -437,6 +426,6 @@ MergeData.Z=ZI; -% mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; -% mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; + % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; + % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; MergeData.Error=0.5*sqrt(sum(Error.^2,3)); errormsg=struct2nc(OutputFile,MergeData);%save result file @@ -461,6 +450,7 @@ A(:,:,2,3)=(R(6)-R(9)*Y)./T; -function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U - +function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) +% convert image coordinates to view angles, after removal of quadratic distorsion +% input in pixel, output in radians X1d=Xd-Ud/2; X2d=Xd+Ud/2; @@ -483,5 +473,4 @@ z=0; error=0; - %% first image Index: /trunk/src/set_object.m =================================================================== --- /trunk/src/set_object.m (revision 953) +++ /trunk/src/set_object.m (revision 954) @@ -111,5 +111,5 @@ if isfield(data,'ProjModeMenu') set(handles.ProjMode,'UserData',data.ProjModeMenu)% data.ProjModeMenu as default menu (used in Type_Callback) - end + end errormsg=fill_GUI(data,handles.set_object); if ~isempty(errormsg) @@ -132,5 +132,5 @@ end end - if isfield(data,'RangeX') + if isfield(data,'RangeX')&& ~strcmp(data.Type,'plane_z')%TODO: generalise if ischar(data.RangeX) data.RangeX=str2num(data.RangeX); @@ -139,5 +139,5 @@ set(handles.num_RangeX_1,'String',num2str(min(data.RangeX),3)) end - if isfield(data,'RangeY') + if isfield(data,'RangeY')&& ~strcmp(data.Type,'plane_z')%TODO: generalise if ischar(data.RangeY) data.RangeY=str2num(data.RangeY); @@ -146,5 +146,5 @@ set(handles.num_RangeY_1,'String',num2str(min(data.RangeY),3)) end - if isfield(data,'RangeZ') + if isfield(data,'RangeZ')&& ~strcmp(data.Type,'plane_z')%TODO: generalise if ischar(data.RangeZ) data.RangeZ=str2num(data.RangeZ); @@ -247,6 +247,4 @@ if isempty(get(handles.ProjMode,'UserData')) switch Type - case {'points','line','plane'} - menu_proj={'projection';'interp_lin';'interp_tps';'none'}; case 'polyline' menu_proj={'interp_lin';'interp_tps';'none'}; @@ -328,5 +326,5 @@ set(handles.num_RangeX_2,'TooltipString',['num_RangeX_2: half length of the ' ObjectStyle]) set(handles.num_RangeY_2,'TooltipString',['num_RangeY_2: half width of the ' ObjectStyle]) - case {'plane'} + case {'plane','plane_z'} set(handles.num_Angle_3,'Visible','on') set(handles.num_RangeX_1,'Visible','on') @@ -387,6 +385,6 @@ set(handles.num_RangeY_2,'String',num2str(UvData.Field.YMax)) end - if isempty(get(handles.CoordUnit,'String')) - set(handles.CoordUnit,'String',Field.CoordUnit) + if isempty(get(handles.CoordUnit,'String'))&& isfield(UvData.Field,'CoordUnit') + set(handles.CoordUnit,'String',UvData.Field.CoordUnit) end end Index: /trunk/src/update_imadoc.m =================================================================== --- /trunk/src/update_imadoc.m (revision 953) +++ /trunk/src/update_imadoc.m (revision 954) @@ -34,18 +34,19 @@ testappend=1; backupfile=outputfile; - testexist=2; - while testexist==2 - backupfile=[backupfile '~']; - testexist=exist(backupfile,'file'); - end - [success,message]=copyfile(outputfile,backupfile);%make backup - if success~=1 - errormsg=['errror in xml file backup: ' message]; - return - end t=xmltree(outputfile); %read the file title=get(t,1,'name'); if strcmp(title,'ImaDoc') - testappend=1; + % testappend=1; + %rename the existing file for backup + testexist=2; + while testexist==2 + backupfile=[backupfile '~']; + testexist=exist(backupfile,'file'); + end + [success,message]=movefile(outputfile,backupfile);%make backup + if success~=1 + errormsg=['errror in xml file backup: ' message]; + return + end %if the xml file is ImaDoc uid_calib=find(t,['ImaDoc/' StructName]); Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 953) +++ /trunk/src/uvmat.m (revision 954) @@ -910,4 +910,11 @@ create_object(data,handles) +% -------------------------------------------------------------------- +function Menuplane_z_Callback(hObject, eventdata, handles) +data.Type='plane_z'; +data.ProjMode='projection';%default +data.ProjModeMenu={};% do not restrict ProjMode menus +create_object(data,handles) + %------------------------------------------------------------------------ function Menuvolume_Callback(hObject, eventdata, handles) @@ -949,12 +956,11 @@ check_plot=0; if isfield(UvData,'Field') - Field=UvData.Field; - if isfield(Field,'NbDim')&& isequal(Field.NbDim,3) + if isfield(UvData.Field,'NbDim')&& isequal(UvData.Field.NbDim,3) data.Coord=[0 0 0]; %default end - if isfield(Field,'CoordUnit') - data.CoordUnit=Field.CoordUnit; - end - if isfield(UvData.Field,'CoordMesh')&&~isempty(UvData.Field.CoordMesh) + if isfield(UvData.Field,'CoordUnit') + data.CoordUnit=UvData.Field.CoordUnit; + end + if isfield(UvData.Field,'CoordMesh')&&~isempty(UvData.Field.CoordMesh)&&~strcmp(data.Type,'plane_z') data.RangeX=[UvData.Field.XMin UvData.Field.XMax]; switch data.Type @@ -984,4 +990,7 @@ data.DX=UvData.Field.CoordMesh; data.DY=UvData.Field.CoordMesh; + end + if isfield(UvData.Field,'ProjModeRequest') + data.ProjMode=UvData.Field.ProjModeRequest;%set the request proj mode option by default end end @@ -5882,2 +5891,5 @@ set(handles.TableDisplay,'Visible','off') end + + +