Changeset 954 for trunk


Ignore:
Timestamp:
Jun 22, 2016, 1:36:50 PM (8 years ago)
Author:
sommeria
Message:

update calib modfied + various updates

Location:
trunk/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/find_field_bounds.m

    r924 r954  
    5454CoordMin=zeros(numel(imax),NbDim);
    5555Mesh=zeros(1,numel(imax));
     56FieldOut.ProjModeRequest='projection';%default
    5657for ind=1:numel(imax)
    5758    if strcmp(CellInfo{imax(ind)}.CoordType,'tps')
    … …  
    9091            Mesh(ind)=min((CoordMax(ind,:)-CoordMin(ind,:))./(NbPoints-1));
    9192    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 
    92101end
    93102Mesh=min(Mesh);

  • trunk/src/mouse_motion.m

    r924 r954  
    385385    XYData=AxeData.CurrentOrigin;
    386386    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
    392393                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 widt
    398             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                 end
    403             end
     394                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
    404405        end
    405406        plot_object(ObjectData,ProjObject,AxeData.CurrentObject,'m');

  • trunk/src/mouse_up.m

    r924 r954  
    9595    else
    9696        switch ObjectData.Type
    97             case {'line'}
     97            case {'line','plane_z'}
    9898                if size(ObjectData.Coord,1)==1 % this is the mouse up for the first point, continue until next click
    9999                    check_multiple=1;
    … …  
    106106                    check_multiple=1;% pass to next mous up if width of height=0
    107107                end
    108             case 'plane' %case of 'plane'
     108            case 'plane' %case of 'plane', TODO: NOT ACTIVATED
    109109                DX=(xy(1,1)-ObjectData.Coord(1,1));
    110110                DY=(xy(1,2)-ObjectData.Coord(1,2));
    … …  
    116116                    end
    117117                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
    118125            otherwise
    119126                check_multiple=1;

  • trunk/src/plot_object.m

    r937 r954  
    133133
    134134%% 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');
     139test_line=ismember(ObjectData.Type,{'points','line','polyline','polygon','plane','plane_z','volume'});
     140test_patch=ismember(ObjectData.ProjMode,{'inside','outside','mask_inside','mask_outside'});
    139141if test_line
    140142    xline=ObjectData.Coord(:,1);
    141143    yline=ObjectData.Coord(:,2);
    142144    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
     234end
     235SubLineStyle='none';%default
     236if 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
     248end
     249if 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);
    262265    end
    263266end
    … …  
    415418    else% no patch image requested, erase existing ones
    416419        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             end
    421         end
     420            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
    422425        end
    423426    end
    … …  
    459462                end
    460463            end
    461         case {'line','polyline','polygon'}
     464        case {'line','polyline','polygon','plane_z'}
    462465            hh=line(xline,yline,'Color',col);
    463466                PlotData.SubObject(1)=line(xinf,yinf,'Color',col,'LineStyle',SubLineStyle,'Tag','proj_object');%draw sub-lines

  • trunk/src/proj_field.m

    r937 r954  
    9090    return
    9191end
    92 ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes
    93 if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case
     92% check list of effective projection modes
     93if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'})
    9494    return
    9595end
    … …  
    117117            [ProjData,errormsg] = proj_line(FieldData,ObjectData);
    118118        end
    119     case 'plane'
     119    case {'plane','plane_z'}
    120120        [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
    121121    case 'volume'
    … …  
    955955test90x=0;%=1 for 90 degree rotation alround x axis
    956956test90y=0;%=1 for 90 degree rotation alround y axis
     957if 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;
     964end   
    957965if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
    958966    test90y=isequal(ObjectData.Angle,[0 90 0]);
    … …  
    987995    return
    988996end
     997InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane
     998if strcmp(ObjectData.Type,'plane_z')
     999    InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation
     1000end
    9891001
    9901002%% extrema along each axis
    … …  
    10711083end
    10721084    icell_grid=[];% field cell index which defines the grid
    1073 if ~strcmp(ObjectData.ProjMode,'projection')
     1085if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize
    10741086    %% define the new coordinates in case of interpolation on a imposed grid
    1075     if ~testYMin
     1087    if ~testYMin 
    10761088        errormsg='min Y value not defined for the projection grid';return
    10771089    end
    … …  
    11191131        AYName='coord_y';
    11201132        AXName='coord_x';
    1121         if strcmp(ObjectData.ProjMode,'projection')
     1133        if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z')
    11221134            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
    11231135            ProjData.coord_x=[FieldData.XMin FieldData.XMax];
    … …  
    13991411            DimValue(DimValue==1)=[];%remove singleton dimensions
    14001412            NbDim=numel(DimValue);%update number of space dimensions
    1401             nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
     1413           % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
    14021414            if NbDim>=3
    14031415                if NbDim>3
    … …  
    14391451                        Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});
    14401452                    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
    14421454                        DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1);
    14431455                    end
    1444                     if numel(Coord{NbDim})==2
     1456                    if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values
    14451457                        DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1);
    14461458                    end
    … …  
    14531465                        end         
    14541466                    else
    1455                         YIndexMax=Coord{NbDim-1}(end)/DY;
     1467                        YIndexMax=numel(Coord{NbDim-1});
    14561468                        YIndexMin=1;
    14571469                    end
    14581470                    if testXMax
    1459                          XIndexMax=(XMax-Coord{NbDim}(1))/DY+1;% matrix index corresponding to the max y value for the new field
     1471                         XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field
    14601472                        if testYMin%test_direct(indY)
    14611473                            XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field
    … …  
    14641476                        end         
    14651477                    else
    1466                         XIndexMax=Coord{NbDim}(end)/DX;
     1478                        XIndexMax=numel(Coord{NbDim});
    14671479                        XIndexMin=1;
    14681480                    end
    … …  
    16271639                        end
    16281640                    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
    16331686                end
    16341687            end
    … …  
    22642317                    end
    22652318                else
     2319                    RotMatrix=rodrigues(om);
     2320                   
    22662321                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
    22672322                    %TODO: use interp3

  • trunk/src/series.m

    r951 r954  
    22922292CheckList_1=1;% indicate whether FieldName_1 has been updated
    22932293handles_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 
     2294if VelTypeRequest && numel(iview_civ)>=1
    22952295    menu=set_veltype_display(SeriesData.FileInfo{iview_civ(1)}.CivStage,SeriesData.FileType{iview_civ(1)});
    22962296    set(handles.VelType,'Value',1)% set first choice by default
    … …  
    23012301    %CheckList=1;
    23022302    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
    23042304        menu=set_veltype_display(SeriesData.FileInfo{iview_civ(2)}.CivStage,SeriesData.FileType{iview_civ(2)});
    23052305        set(handles.VelType_1,'Value',1)% set first choice by default
    … …  
    23172317    set(handles.VelType,'Visible','off')
    23182318    set(handles.VelType_title,'Visible','off')
    2319 end   
     2319end
    23202320
    23212321%% Detect the types of input files and set menus and default options in 'FieldName'
    23222322if (FieldNameRequest || VelTypeRequest) && numel(iview_netcdf)>=1
    2323     set(handles.InputFields,'Visible','on')   
     2323    set(handles.InputFields,'Visible','on')% set the frame InputFields visible
    23242324    if FieldNameRequest && isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName')
    23252325        set(handles.FieldName,'Visible','on')
    … …  
    23462346            end
    23472347        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           
    23512377            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)
    23812380        end
    23822381    else
     2382        set(handles.FieldName_1,'Visible','off')
     2383    end
     2384    if isempty(FieldList)
    23832385        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
    23862390else
    23872391    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%               
    214%------------------------------------------------------------------------
    315% function ParamOut=civ2vel_3C(Param)
    … …  
    719%
    820%INPUT:
     21%
    922% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
    1023% In batch mode, Param is the name of the corresponding xml file containing the same information
    … …  
    1528% see the current structure Param)
    1629%    .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
    1832%    .OutputSubDir: name of the subdirectory for data outputs
    1933%    .OutputDirExt: directory extension for data outputs
    … …  
    2236%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
    2337%             .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         
    2639%    .IndexRange: set the file or frame indices on which the action must be performed
    2740%    .InputFields: sub structure describing the input fields withfields
    … …  
    5265
    5366function ParamOut=civ2vel_3C(Param)
    54 disp('test')
     67
    5568%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed
    5669if isstruct(Param) && isequal(Param.Action.RUN,0)
    … …  
    124137%% define the directory for result file (with path=RootPath{1})
    125138OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    126 %
    127 % if ~isfield(Param,'InputFields')
    128 %     Param.InputFields.FieldName='';
    129 % end
    130139
    131140%% calibration data and timing: read the ImaDoc files
    … …  
    160169    return
    161170end
    162 ObjectData=Param.ProjObject;
     171ObjectData=Param.ProjObject;% Object attached to the GUI series
    163172xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
    164173yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
    … …  
    168177W=zeros(size(XI,1),size(XI,2));
    169178
    170 %% MAIN LOOP ON FIELDS
     179%%%% MAIN LOOP ON FIELDS %%%%%
    171180warning off
    172181
    … …  
    175184    CheckOverwrite=Param.CheckOverwrite;
    176185end
    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);
     186for 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);
    186192    if ~isempty(i2_series{end})
    187         i2=i2_series{end}(index);
     193        i2=i2_series{end}(field_index);
    188194    else
    189195        i2=i1;
    … …  
    192198    j2=1;
    193199    if ~isempty(j1_series{1})
    194         j1=j1_series{1}(index);
     200        j1=j1_series{1}(field_index);
    195201        if ~isempty(j2_series{end})
    196             j2=j2_series{end}(index);
     202            j2=j2_series{end}(field_index);
    197203        else
    198204            j2=j1;
    … …  
    201207    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2);
    202208   
    203     %%
    204    
    205    
     209    %% program stop if requested on the GUI
    206210    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    207211        disp('program stopped by user')
    … …  
    209213    end
    210214   
    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   
    216221    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    217222    Data=cell(1,NbView);%initiate the set Data
    218     timeread=zeros(1,NbView);
    219223   
    220224    %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken
    221225    %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
    226228        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           
    247322        end
    248323       
     324        ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid
    249325       
    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
    334328    Vtest=ZItemp(:,:,2)-ZItemp(:,:,1);
    335329   
    … …  
    337331    [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b
    338332   
    339    
     333    
    340334    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)'},'*');
    343337        if ~isempty(errormsg)
    344338            disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)
    … …  
    359353        end
    360354    end
    361     %remove wrong vector
     355    %remove false vectors
    362356    temp=find(Data{1}.FF==0);
    363357    X1=Data{1}.X(temp);
    … …  
    369363    Va=griddata(X1,Y1,V1,Xa,Ya);
    370364   
    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
    372366    [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~
    373367   
    374     %remove wrong vector
     368    %remove false vectors
    375369    temp=find(Data{2}.FF==0);
    376370    X2=Data{2}.X(temp);
    … …  
    380374    Ub=griddata(X2,Y2,U2,Xb,Yb);
    381375    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
    384378    [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~
    385    
     379    
    386380   
    387381    % System to solve
    388382    S=ones(size(XI,1),size(XI,2),3);
    389383    D=ones(size(XI,1),size(XI,2),3,3);
    390 
     384   
    391385    S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
    392386    S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
    … …  
    408402            W(indj,indi)=dxyz(3);
    409403        end
    410     end   
     404    end
    411405    Error=zeros(size(XI,1),size(XI,2),4);
    412406    Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua;
    … …  
    415409    Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb;
    416410   
    417    
    418 
    419    
    420  
    421    
    422411    %% recording the merged field
    423     if index==1% initiate the structure at first index
     412    if field_index==1% initiate the structure at first index
    424413        MergeData.ListGlobalAttribute={'Conventions','Time','Dt'};
    425414        MergeData.Conventions='uvmat';
    … …  
    428417        MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'};
    429418        MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
    430                 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     419            {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
    431420        MergeData.coord_x=xI;
    432421        MergeData.coord_y=yI;
    … …  
    437426    MergeData.Z=ZI;
    438427   
    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;
    441430    MergeData.Error=0.5*sqrt(sum(Error.^2,3));
    442431    errormsg=struct2nc(OutputFile,MergeData);%save result file
    … …  
    461450A(:,:,2,3)=(R(6)-R(9)*Y)./T;
    462451
    463 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X  and Ud to U
    464 
     452function [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
    465455X1d=Xd-Ud/2;
    466456X2d=Xd+Ud/2;
    … …  
    483473z=0;
    484474error=0;
    485 
    486475
    487476%% first image

  • trunk/src/set_object.m

    r924 r954  
    111111    if isfield(data,'ProjModeMenu')
    112112        set(handles.ProjMode,'UserData',data.ProjModeMenu)% data.ProjModeMenu as default menu (used in Type_Callback)
    113     end
     113    end     
    114114    errormsg=fill_GUI(data,handles.set_object);
    115115    if ~isempty(errormsg)
    … …  
    132132        end
    133133    end
    134     if isfield(data,'RangeX')
     134    if isfield(data,'RangeX')&& ~strcmp(data.Type,'plane_z')%TODO: generalise
    135135        if ischar(data.RangeX)
    136136            data.RangeX=str2num(data.RangeX);
    … …  
    139139        set(handles.num_RangeX_1,'String',num2str(min(data.RangeX),3))
    140140    end
    141     if isfield(data,'RangeY')
     141    if isfield(data,'RangeY')&& ~strcmp(data.Type,'plane_z')%TODO: generalise
    142142        if ischar(data.RangeY)
    143143            data.RangeY=str2num(data.RangeY);
    … …  
    146146        set(handles.num_RangeY_1,'String',num2str(min(data.RangeY),3))
    147147    end
    148     if isfield(data,'RangeZ')
     148    if isfield(data,'RangeZ')&& ~strcmp(data.Type,'plane_z')%TODO: generalise
    149149        if ischar(data.RangeZ)
    150150            data.RangeZ=str2num(data.RangeZ);
    … …  
    247247if isempty(get(handles.ProjMode,'UserData'))
    248248    switch Type
    249         case {'points','line','plane'}
    250             menu_proj={'projection';'interp_lin';'interp_tps';'none'};
    251249        case 'polyline'
    252250            menu_proj={'interp_lin';'interp_tps';'none'};
    … …  
    328326        set(handles.num_RangeX_2,'TooltipString',['num_RangeX_2: half length of the ' ObjectStyle])
    329327        set(handles.num_RangeY_2,'TooltipString',['num_RangeY_2: half width of the ' ObjectStyle])
    330     case {'plane'} 
     328    case {'plane','plane_z'} 
    331329        set(handles.num_Angle_3,'Visible','on')
    332330        set(handles.num_RangeX_1,'Visible','on')
    … …  
    387385            set(handles.num_RangeY_2,'String',num2str(UvData.Field.YMax))
    388386        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)
    391389        end       
    392390    end

  • trunk/src/update_imadoc.m

    r924 r954  
    3434    testappend=1;
    3535    backupfile=outputfile;
    36     testexist=2;
    37     while testexist==2
    38         backupfile=[backupfile '~'];
    39         testexist=exist(backupfile,'file');
    40     end
    41     [success,message]=copyfile(outputfile,backupfile);%make backup
    42     if success~=1
    43         errormsg=['errror in xml file backup: ' message];
    44         return
    45     end
    4636    t=xmltree(outputfile); %read the file
    4737    title=get(t,1,'name');
    4838    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
    5051        %if the xml file is  ImaDoc
    5152        uid_calib=find(t,['ImaDoc/' StructName]);

  • trunk/src/uvmat.m

    r951 r954  
    910910create_object(data,handles)
    911911
     912% --------------------------------------------------------------------
     913function Menuplane_z_Callback(hObject, eventdata, handles)
     914data.Type='plane_z';
     915data.ProjMode='projection';%default
     916data.ProjModeMenu={};% do not restrict ProjMode menus
     917create_object(data,handles)
     918
    912919%------------------------------------------------------------------------
    913920function Menuvolume_Callback(hObject, eventdata, handles)
    … …  
    949956check_plot=0;
    950957if 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)
    953959         data.Coord=[0 0 0]; %default
    954960    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')
    959965        data.RangeX=[UvData.Field.XMin UvData.Field.XMax];
    960966        switch data.Type
    … …  
    984990        data.DX=UvData.Field.CoordMesh;
    985991        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
    986995    end
    987996end
    … …  
    58825891    set(handles.TableDisplay,'Visible','off')
    58835892end
     5893
     5894
     5895
Note: See TracChangeset for help on using the changeset viewer.