Changeset 77 for trunk/src


Ignore:
Timestamp:
Apr 2, 2010, 8:57:13 AM (15 years ago)
Author:
sommeria
Message:

introduction of volume projection on a regular grid (mode 'interp') using griddata3. It works but very slow.
minor bug corrections
todo: problem with 3D projection when the same field is read twice (e;g. w used for vector color)

Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r72 r77  
    121121%     FieldData.VarAttribute={};
    122122% end
    123 
    124 if isequal(ObjectData.Style,'points')
     123switch ObjectData.Style
     124    case 'points'
    125125    [ProjData,errormsg]=proj_points(FieldData,ObjectData);
    126 elseif isequal(ObjectData.Style,'line')|isequal(ObjectData.Style,'polyline')
     126    case {'line','polyline'}
    127127     [ProjData,errormsg] = proj_line(FieldData,ObjectData);
    128 elseif isequal(ObjectData.Style,'polygon')|isequal(ObjectData.Style,'rectangle')|isequal(ObjectData.Style,'ellipse')
    129     if isequal(ObjectData.ProjMode,'inside')|isequal(ObjectData.ProjMode,'outside')
    130         [ProjData,errormsg] = proj_patch(FieldData,ObjectData);
    131     else
    132         [ProjData,errormsg] = proj_line(FieldData,ObjectData);
    133     end
     128    case {'polygon','rectangle','ellipse'}
     129        if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')
     130            [ProjData,errormsg] = proj_patch(FieldData,ObjectData);
     131        else
     132            [ProjData,errormsg] = proj_line(FieldData,ObjectData);
     133        end
    134134     %A FAIRE : GERER MASK
    135 elseif isequal(ObjectData.Style,'plane')
     135    case 'plane'
    136136%         if isfield(FieldData,'NbDim') & isequal(FieldData.NbDim,3)
    137137%             ProjData= proj_plane3D(FieldData,ObjectData);%
    … …  
    139139            [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
    140140%         end
     141    case 'volume'
     142        [ProjData,errormsg] = proj_volume(FieldData,ObjectData);
    141143end
    142144if exist('IndexObj','var')
    … …  
    154156    width=ObjectData.Range(1,2);
    155157end
    156 if isfield(ObjectData,'RangeX')&~isempty(ObjectData.RangeX)
     158if isfield(ObjectData,'RangeX')&&~isempty(ObjectData.RangeX)
    157159    width=max(ObjectData.RangeX);
    158160end
    159 if isfield(ObjectData,'RangeY')&~isempty(ObjectData.RangeY)
     161if isfield(ObjectData,'RangeY')&&~isempty(ObjectData.RangeY)
    160162    width=max(width,max(ObjectData.RangeY));
    161163end
    162 if isfield(ObjectData,'RangeZ')&~isempty(ObjectData.RangeZ)
     164if isfield(ObjectData,'RangeZ')&&~isempty(ObjectData.RangeZ)
    163165    width=max(width,max(ObjectData.RangeZ));
    164166end
    … …  
    951953test3D=0;
    952954if isfield(FieldData,'nb_dim')
    953     test3D=isequal(FieldData.nb_dim,3)
     955    test3D=isequal(FieldData.nb_dim,3);
    954956end
    955957test3C=test3D; %default 3 vel components
    … …  
    10461048%case of input fields with unstructured coordinates
    10471049    if testX
    1048         XName=FieldData.ListVarName{ivar_X}
    1049         YName=FieldData.ListVarName{ivar_Y}
     1050        XName=FieldData.ListVarName{ivar_X};
     1051        YName=FieldData.ListVarName{ivar_Y};
    10501052        eval(['coord_x=FieldData.' XName ';'])
    10511053        eval(['coord_y=FieldData.' YName ';'])
    10521054        if length(ivar_Z)==1
    1053             ZName=FieldData.ListVarName{ivar_Z}
     1055            ZName=FieldData.ListVarName{ivar_Z};
    10541056            eval(['coord_z=FieldData.' ZName ';'])
    10551057        end
    … …  
    10681070            indcut=find(abs(fieldZ) <= width);
    10691071            for ivar=VarIndex
    1070                 VarName=FieldData.ListVarName{ivar};
     1072                 VarName=FieldData.ListVarName{ivar};
     1073%                 eval(['size(FieldData.' VarName ')'])
    10711074                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
    10721075%                 end     
    … …  
    11301133        % different cases of projection
    11311134        if isequal(ObjectData.ProjMode,'projection')
    1132             %ProjData.ListDimName=[ProjData.ListDimName
    1133             %FieldData.ListDimName(DimIndices(1))];%add the point index to
    11341135            %the list of dimension
    11351136            ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    … …  
    11581159            coord_y_proj=[YMin:DY:YMax];
    11591160            DimCell={'coord_y','coord_x'};
    1160             %ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)];
    11611161            ProjData.ListVarName={'coord_y','coord_x'};
    11621162            ProjData.VarDimName={'coord_y','coord_x'};   
    1163             nbcoord=2;
    1164             %ProjData.VarDimIndex={};   
     1163            nbcoord=2; 
    11651164            ProjData.coord_y=[YMin YMax];
    11661165            ProjData.coord_x=[XMin XMax];
    … …  
    11871186                    ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    11881187                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1189                     %ProjData.VarDimIndex=[ProjData.VarDimIndex {[1 2]}];
    11901188                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    11911189                        ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
    11921190                    end
    1193 %                     ProjData.VarAttribute{ivar_new}.Coord_2=[XMin XMax];
    1194 %                     ProjData.VarAttribute{ivar_new}.Coord_1=[YMin YMax];
    11951191                    if  ~isequal(ivar_FF,0)
    11961192                        eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
    … …  
    12041200                        eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
    12051201                        FF(indnan)=ones(size(indnan));
    1206 %                         eval(['ProjData.' VarName '(indnan)=zeros(size(indnan));'])%put NaN to 0
    1207 %                         FF=FF|FFlag;
    12081202                        testFF=1;
    12091203                    end
    … …  
    12221216                ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    12231217                ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1224               %  ProjData.VarDimIndex=[ProjData.VarDimIndex {[1 2]}];
    12251218               ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    12261219                ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    … …  
    12351228        VarName=FieldData.ListVarName{VarIndex(1)};
    12361229        eval(['DimValue=size(FieldData.' VarName ');'])
    1237        
    1238         %ListDimName=FieldData.ListDimName(DimIndices);
    12391230        ListDimName=FieldData.VarDimName{VarIndex(1)};
    12401231        ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells)
    … …  
    13941385            for ivar=VarIndex
    13951386                VarName=FieldData.ListVarName{ivar};
    1396 %                 if isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')% coordinates common to all fields
    1397 %                     ProjData.ListDimName={'coord_y','coord_x'};
    1398 %                     ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)];
    1399 %                 else
    1400 %                     icoord=icoord+1;
    1401 %                     ProjData.ListDimName=[ProjData.ListDimName {['coord_y_' num2str(icoord)],['coord_x_' num2str(icoord)]}];
    1402 %                     ProjData.DimValue=[ProjData.DimValue length(coord_y_proj) length(coord_x_proj)];
    1403 %                 end
    14041387                ProjData.ListVarName=[ProjData.ListVarName VarName];
    14051388                ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
    … …  
    14071390                    ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
    14081391                end
    1409 %                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=Xbound;
    1410 %                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=Ybound;
    14111392                eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']);
    14121393            end         
    … …  
    14521433                        vec_B(ind_out,icolor)=zeros(size(ind_out));
    14531434                    end
    1454                     % TODO: A REVOIR:
    1455 %                     if isequal(ObjectData.ProjMode,'interp') || isequal(ObjectData.ProjMode,'filter')% coordinates common to all fields
    1456 %                         ProjData.ListDimName={'coord_y','coord_x'};
    1457 %                         ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)];
    1458 %                     else
    1459 %                         icoord=icoord+1;
    1460 %                         ProjData.ListDimName=[ProjData.ListDimName {['coord_y_' num2str(icoord)],['coord_x_' num2str(icoord)]}];
    1461 %                         ProjData.DimValue=[ProjData.DimValue length(coord_y_proj) length(coord_x_proj)];
    1462 %                     end
    14631435                    ProjData.ListVarName=[ProjData.ListVarName VarName];                 
    1464                    % ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
    14651436                    if length(FieldData.VarAttribute)>=ivar
    14661437                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    1467                     end
    1468                     % A REVOIR:
    1469 %                     if test_direct(2)==1
    1470 %                         ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMin XMax];
    1471 %                     else
    1472 %                         ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMax XMin];
    1473 %                     end
    1474 %                     if test_direct(1)==1
    1475 %                         ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMin YMax];
    1476 %                     else
    1477 %                         ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMax YMin];
    1478 %                     end     
     1438                    end     
    14791439                    eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
    14801440                end
    14811441                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
    14821442                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
    1483               %  ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
    14841443                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
    14851444            else %3D case
    … …  
    14941453                            VarName=FieldData.ListVarName{ivar};
    14951454                            ProjData.ListVarName=[ProjData.ListVarName VarName];
    1496                           %  ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-2 nb_dim-1]];
    1497                             ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    1498                             % A REVOIR
    1499 %                             if test_direct_x==1
    1500 %                                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMin XMax];
    1501 %                             else
    1502 %                                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMax XMin];
    1503 %                             end
    1504 %                             if test_direct_y==1
    1505 %                                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMin YMax];
    1506 %                             else
    1507 %                                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMax YMin];
    1508 %                             end     
     1455                            ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
    15091456                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
    15101457                            %TODO : do a vertical average for a thick plane
    … …  
    15451492
    15461493%-----------------------------------------------------------------
     1494%project in a volume
     1495% A FAIRE ....(copie de proj_plane)
     1496 function  [ProjData,errormsg] = proj_volume(FieldData, ObjectData)
     1497%-----------------------------------------------------------------
     1498
     1499%initialisation of the input parameters of the projection plane
     1500%-----------------------------------------------------------------
     1501ProjMode='projection';%direct projection by default
     1502if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
     1503
     1504%axis origin
     1505if isempty(ObjectData.Coord)
     1506    ObjectData.Coord(1,1)=0;%origin of the volume set to [0 0] by default
     1507    ObjectData.Coord(1,2)=0;
     1508    ObjectData.Coord(1,3)=0;
     1509end
     1510
     1511%rotation angles
     1512Phi=0;%default
     1513Theta=0;
     1514Psi=0;
     1515if isfield(ObjectData,'Phi')&& ~isempty(ObjectData.Phi)
     1516    Phi=(pi/180)*ObjectData.Phi;%first Euler angle in radian
     1517end
     1518if isfield(ObjectData,'Theta')&& ~isempty(ObjectData.Theta)
     1519    Theta=(pi/180)*ObjectData.Theta;%second Euler angle in radian
     1520end
     1521if isfield(ObjectData,'Psi')&& ~isempty(ObjectData.Psi)
     1522    Psi=(pi/180)*ObjectData.Psi;%third Euler angle in radian
     1523end
     1524
     1525%components of the unity vector normal to the projection plane
     1526NormVec_X=-sin(Phi)*sin(Theta);
     1527NormVec_Y=cos(Phi)*sin(Theta);
     1528NormVec_Z=cos(Theta);
     1529
     1530% test for 3D fields
     1531test3D=0;
     1532if isfield(FieldData,'nb_dim')
     1533    test3D=isequal(FieldData.nb_dim,3);
     1534end
     1535test3C=test3D; %default 3 vel components
     1536
     1537%mesh sizes DX and DY
     1538DX=0;
     1539DY=0; %default
     1540if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
     1541     DX=abs(ObjectData.DX);%mesh of interpolation points
     1542end
     1543if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
     1544     DY=abs(ObjectData.DY);%mesh of interpolation points
     1545end
     1546if  ~isequal(ProjMode,'projection') & (DX==0|DY==0)
     1547        errormsg='DX or DY missing';
     1548        display(errormsg)
     1549        return
     1550end
     1551if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
     1552     DZ=abs(ObjectData.DZ);%mesh of interpolation points
     1553end
     1554
     1555%extrema along each axis
     1556testXMin=0;
     1557testXMax=0;
     1558testYMin=0;
     1559testYMax=0;
     1560if isfield(ObjectData,'RangeX')
     1561        XMin=min(ObjectData.RangeX);
     1562        XMax=max(ObjectData.RangeX);
     1563        testXMin=XMax>XMin;
     1564        testXMax=1;
     1565end
     1566if isfield(ObjectData,'RangeY')
     1567        YMin=min(ObjectData.RangeY);
     1568        YMax=max(ObjectData.RangeY);
     1569        testYMin=YMax>YMin;
     1570        testYMax=1;
     1571end
     1572width=0;%default width of the projection band
     1573if isfield(ObjectData,'RangeZ')
     1574      ZMin=min(ObjectData.RangeZ);
     1575        ZMax=max(ObjectData.RangeZ);
     1576        testZMin=YMax>YMin;
     1577        testZMax=1;
     1578%         width=max(ObjectData.RangeZ);
     1579end
     1580
     1581% initiate Matlab  structure for physical field
     1582ProjData=proj_heading(FieldData,ObjectData);
     1583ProjData.NbDim=3;
     1584ProjData.ListDimName={};%name of dimension
     1585ProjData.DimValue=[];%values of dimension (nbre of vectors)
     1586ProjData.ListVarName={};
     1587ProjData.VarDimName={};
     1588
     1589error=0;%default
     1590flux=0;
     1591testfalse=0;
     1592ListIndex={};
     1593
     1594%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     1595%group the variables (fields of 'FieldData') in cells of variables with the same dimensions
     1596%-----------------------------------------------------------------
     1597idimvar=0;
     1598[CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_indices(FieldData);
     1599if ~isempty(errormsg)
     1600    errormsg=['error in proj_field/proj_plane:' errormsg];
     1601    return
     1602end
     1603%LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
     1604% CellVarIndex=cells of variable index arrays
     1605ivar_new=0; % index of the current variable in the projected field
     1606icoord=0;
     1607nbcoord=0;%number of added coordinate variables brought by projection
     1608for icell=1:length(CellVarIndex)
     1609    if NbDim(icell)<2
     1610        continue
     1611    end
     1612    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
     1613    VarType=VarTypeCell{icell};
     1614    ivar_X=VarType.coord_x;
     1615    ivar_Y=VarType.coord_y;
     1616    ivar_Z=VarType.coord_z;
     1617    ivar_U=VarType.vector_x;
     1618    ivar_V=VarType.vector_y;
     1619    ivar_W=VarType.vector_z;
     1620    ivar_C=VarType.scalar ;
     1621    ivar_Anc=VarType.ancillary;
     1622    test_anc=zeros(size(VarIndex));
     1623    test_anc(ivar_Anc)=ones(size(ivar_Anc));
     1624    ivar_F=VarType.warnflag;
     1625    ivar_FF=VarType.errorflag;
     1626    testX=~isempty(ivar_X) && ~isempty(ivar_Y);
     1627    DimCell=FieldData.VarDimName{VarIndex(1)};
     1628    if ischar(DimCell)
     1629        DimCell={DimCell};%name of dimensions
     1630    end
     1631   
     1632%case of input fields with unstructured coordinates
     1633    if testX
     1634        XName=FieldData.ListVarName{ivar_X};
     1635        YName=FieldData.ListVarName{ivar_Y};
     1636        eval(['coord_x=FieldData.' XName ';'])
     1637        eval(['coord_y=FieldData.' YName ';'])
     1638        if length(ivar_Z)==1
     1639            ZName=FieldData.ListVarName{ivar_Z};
     1640            eval(['coord_z=FieldData.' ZName ';'])
     1641        end
     1642
     1643        % translate  initial coordinates
     1644        coord_x=coord_x-ObjectData.Coord(1,1);
     1645        coord_y=coord_y-ObjectData.Coord(1,2);
     1646        if ~isempty(ivar_Z)
     1647            coord_z=coord_z-ObjectData.Coord(1,3);
     1648        end
     1649       
     1650        % selection of the vectors in the projection range (3D case)
     1651        if length(ivar_Z)==1 &&  width > 0
     1652            %components of the unitiy vector normal to the projection plane
     1653            fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane           
     1654            indcut=find(abs(fieldZ) <= width);
     1655            for ivar=VarIndex
     1656                VarName=FieldData.ListVarName{ivar};
     1657                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1658%                 end     
     1659                    % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE             
     1660            end
     1661            coord_x=coord_x(indcut);
     1662            coord_y=coord_y(indcut);
     1663            coord_z=coord_z(indcut);
     1664        end
     1665
     1666       %rotate coordinates if needed
     1667        if isequal(Phi,0)
     1668            coord_X=coord_x;
     1669            coord_Y=coord_y;
     1670            if ~isequal(Theta,0)
     1671                coord_Y=coord_Y *cos(Theta);
     1672            end
     1673        else
     1674            coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
     1675            coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
     1676        end
     1677        if ~isempty(ivar_Z)
     1678            coord_Y=coord_Y+coord_z *sin(Theta);
     1679        end
     1680        if ~isequal(Psi,0)
     1681                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1682                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1683        end
     1684        coord_Z=coord_z; %TO UPDATE
     1685        %restriction to the range of x and y if imposed
     1686        testin=ones(size(coord_X)); %default
     1687        testbound=0;
     1688        if testXMin
     1689            testin=testin & (coord_X >= XMin);
     1690            testbound=1;
     1691        end
     1692        if testXMax
     1693            testin=testin & (coord_X <= XMax);
     1694            testbound=1;
     1695        end
     1696        if testYMin
     1697            testin=testin & (coord_Y >= YMin);
     1698            testbound=1;
     1699        end
     1700        if testYMin
     1701            testin=testin & (coord_Y <= YMax);
     1702            testbound=1;
     1703        end
     1704        if testbound
     1705            indcut=find(testin);
     1706            for ivar=VarIndex
     1707                VarName=FieldData.ListVarName{ivar};
     1708                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])           
     1709            end
     1710            coord_X=coord_X(indcut);
     1711            coord_Y=coord_Y(indcut);
     1712            if length(ivar_Z)==1
     1713                coord_Z=coord_Z(indcut);
     1714            end
     1715        end
     1716        % different cases of projection
     1717        if isequal(ObjectData.ProjMode,'projection')
     1718            %the list of dimension
     1719            ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
     1720            ProjData.DimValue=[ProjData.DimValue length(coord_X)];
     1721            nbvar=0;
     1722            for ivar=VarIndex %transfer variables to the projection plane
     1723                VarName=FieldData.ListVarName{ivar};
     1724                if ivar==ivar_X %x coordinate
     1725                    eval(['ProjData.' VarName '=coord_X;'])
     1726                elseif ivar==ivar_Y % y coordinate
     1727                    eval(['ProjData.' VarName '=coord_Y;'])
     1728                elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
     1729                    eval(['ProjData.' VarName '=FieldData.' VarName ';'])
     1730                end
     1731                if isempty(ivar_Z) || ivar~=ivar_Z
     1732                    ProjData.ListVarName=[ProjData.ListVarName VarName];
     1733                    ProjData.VarDimName=[ProjData.VarDimName DimCell];
     1734                    nbvar=nbvar+1;
     1735                    if isfield(FieldData,'VarAttribute') & length(FieldData.VarAttribute) >=ivar
     1736                        ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
     1737                    end
     1738                end
     1739            end 
     1740        elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
     1741            coord_x_proj=[XMin:DX:XMax];
     1742            coord_y_proj=[YMin:DY:YMax];
     1743            coord_z_proj=[ZMin:DZ:ZMax];
     1744            [XI,YI,ZI]=meshgrid(coord_x_proj,coord_y_proj,coord_z_proj);
     1745            DimCell={'coord_y','coord_x'};
     1746            ProjData.ListVarName={'coord_z','coord_y','coord_x'};
     1747            ProjData.VarDimName={'coord_z','coord_y','coord_x'};   
     1748            nbcoord=3; 
     1749            ProjData.coord_z=[ZMin ZMax];
     1750            ProjData.coord_y=[YMin YMax];
     1751            ProjData.coord_x=[XMin XMax];
     1752            if isempty(ivar_X), ivar_X=0; end;
     1753            if isempty(ivar_Y), ivar_Y=0; end;
     1754            if isempty(ivar_Z), ivar_Z=0; end;
     1755            if isempty(ivar_U), ivar_U=0; end;
     1756            if isempty(ivar_V), ivar_V=0; end;
     1757            if isempty(ivar_W), ivar_W=0; end;
     1758            if isempty(ivar_F), ivar_F=0; end;
     1759            if isempty(ivar_FF), ivar_FF=0; end;
     1760            if ~isequal(ivar_FF,0)
     1761                VarName_FF=FieldData.ListVarName{ivar_FF};
     1762                eval(['indsel=find(FieldData.' VarName_FF '==0);'])
     1763                coord_X=coord_X(indsel);
     1764                coord_Y=coord_Y(indsel);
     1765            end
     1766            FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
     1767            testFF=0;
     1768            size(XI)
     1769            size(YI)
     1770            size(ZI)
     1771            for ivar=VarIndex
     1772                VarName=FieldData.ListVarName{ivar};
     1773                if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
     1774                    ivar_new=ivar_new+1;
     1775                    ProjData.ListVarName=[ProjData.ListVarName {VarName}];
     1776                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1777                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
     1778                        ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
     1779                    end
     1780                    if  ~isequal(ivar_FF,0)
     1781                        eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
     1782                    end
     1783                    eval(['ProjData.' VarName '=griddata3(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.' VarName '),XI,YI,ZI);'])
     1784%                     eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));'])
     1785%                     FFlag= isnan(varline); %detect undefined values NaN
     1786%                     indnan=find(FFlag);
     1787%                     if~isempty(indnan)
     1788%                         varline(indnan)=zeros(size(indnan));
     1789%                         eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
     1790%                         FF(indnan)=ones(size(indnan));
     1791%                         testFF=1;
     1792%                     end
     1793%                     if ivar==ivar_U
     1794%                         ivar_U=ivar_new;
     1795%                     end
     1796%                     if ivar==ivar_V
     1797%                         ivar_V=ivar_new;
     1798%                     end
     1799%                     if ivar==ivar_W
     1800%                         ivar_W=ivar_new;
     1801%                     end
     1802                end
     1803            end
     1804            if testFF
     1805                ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
     1806                ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
     1807               ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1808                ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
     1809            end
     1810        end
     1811%case of fields defined on a structured  grid
     1812    else 
     1813        AYName=FieldData.ListVarName{VarType.coord(1)};
     1814        AXName=FieldData.ListVarName{VarType.coord(2)};
     1815        eval(['AX=FieldData.' AXName ';'])
     1816        eval(['AY=FieldData.' AYName ';'])
     1817        VarName=FieldData.ListVarName{VarIndex(1)};
     1818        eval(['DimValue=size(FieldData.' VarName ');'])
     1819        ListDimName=FieldData.VarDimName{VarIndex(1)};
     1820        ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells)
     1821        ProjData.VarDimName=[{AYName} {AXName} ProjData.VarDimName];
     1822        nbcolor=1; %default
     1823        for idim=1:length(ListDimName)
     1824            DimName=ListDimName{idim};
     1825            if isequal(DimName,'rgb')|isequal(DimName,'nb_coord')|isequal(DimName,'nb_coord_i')
     1826               nbcolor=DimValue(idim);
     1827               DimIndices(idim)=[];
     1828               DimValue(idim)=[];
     1829            end
     1830            if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
     1831                DimIndices(idim)=[];
     1832                DimValue(idim)=[];
     1833            end
     1834        end 
     1835        ind_1=find(DimValue==1);
     1836        DimIndices(ind_1)=[]; %suppress singleton dimensions
     1837%         indxy=find(DimVarIndex(DimIndices));%select dimension variables (DimIndices non zero)
     1838        nb_dim=length(DimIndices);%number of space dimensions
     1839        Coord_z=[];
     1840        Coord_y=[];
     1841        Coord_x=[];   
     1842   
     1843        for idim=1:nb_dim %loop on space dimensions
     1844            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
     1845            test_coord(idim)=0;%test for defined coordinates, =0 by default
     1846            ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
     1847            if ~isequal(ivar,0)%  a variable corresponds to the current dimension
     1848                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% position for the first index
     1849                DCoord=diff(Coord{idim});
     1850                DCoord_min(idim)=min(DCoord);
     1851                DCoord_max=max(DCoord);
     1852                test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
     1853                test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
     1854                if ~isequal(test_direct(idim),test_direct_min)
     1855                     msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
     1856                                return
     1857                end               
     1858                test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
     1859                test_coord(idim)=1;
     1860
     1861            else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
     1862                Coord_i_str=['Coord_' num2str(idim)];
     1863                DCoord_min(idim)=1;%default
     1864                Coord{idim}=[0.5 DimValue(idim)-0.5];
     1865                test_direct(idim)=1;
     1866            end
     1867        end
     1868        if nb_dim==2
     1869            if DY==0
     1870                DY=abs(DCoord_min(1));
     1871            end
     1872            npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation
     1873            npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
     1874            if DX==0
     1875                DX=abs(DCoord_min(2));
     1876            end
     1877            npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol 
     1878            npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid
     1879            Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
     1880            test_direct_y=test_direct(1);
     1881            Coord_x=linspace(Coord{2}(1),Coord{2}(end),npX);
     1882            test_direct_x=test_direct(2);
     1883            DAX=DCoord_min(2);
     1884            DAY=DCoord_min(1);
     1885        elseif nb_dim==3
     1886            DZ=abs(DCoord_min(1));
     1887            npz=1+round(abs(Coord{1}(end)-Coord{1}(1))/DZ);%nbre of points after interpolation
     1888            if DY==0
     1889                DY=abs(DCoord_min(2));
     1890            end
     1891            npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol
     1892            npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol
     1893            if DX==0
     1894                DX=abs(DCoord_min(3));
     1895            end
     1896            npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
     1897            npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol
     1898            Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
     1899            test_direct_z=test_direct(1);
     1900            Coord_y=linspace(Coord{2}(1),Coord{2}(end),npY);
     1901            test_direct_y=test_direct(2);
     1902            Coord_x=linspace(Coord{3}(1),Coord{3}(end),npX);
     1903            test_direct_x=test_direct(3);
     1904        end 
     1905        minAX=min(Coord_x);
     1906        maxAX=max(Coord_x);
     1907        minAY=min(Coord_y);
     1908        maxAY=max(Coord_y);
     1909        xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
     1910        ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
     1911        xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame
     1912        ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi);
     1913        if ~testXMax
     1914            XMax=max(xcor_new);
     1915        end
     1916        if ~testXMin
     1917            XMin=min(xcor_new);
     1918        end
     1919        if ~testYMax
     1920            YMax=max(ycor_new);
     1921        end
     1922        if ~testYMin
     1923            YMin=min(ycor_new);
     1924        end
     1925        DXinit=(maxAX-minAX)/(npx-1);
     1926        DYinit=(maxAY-minAY)/(npy-1);
     1927        if DX==0
     1928            DX=DXinit;
     1929        end
     1930        if DY==0
     1931            DY=DYinit;
     1932        end
     1933        npX=floor((XMax-XMin)/DX+1);
     1934        npY=floor((YMax-YMin)/DY+1);   
     1935        if test_direct_y
     1936            coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
     1937        else
     1938            coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
     1939        end
     1940        if test_direct_x
     1941            coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
     1942        else
     1943            coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
     1944        end
     1945       
     1946        % case with no rotation and interpolation
     1947        if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0)
     1948            if test_direct(1)
     1949                min_ind1=ceil((YMin-Coord{1}(1))/DYinit)+1;
     1950                max_ind1=floor((YMax-Coord{1}(1))/DYinit)+1;
     1951                Ybound(1)=Coord{1}(1)+DYinit*(min_ind1-1);
     1952                Ybound(2)=Coord{1}(1)+DYinit*(max_ind1-1);
     1953            else
     1954                min_ind1=ceil((Coord{1}(1)-YMax)/DYinit)+1;
     1955                max_ind1=floor((Coord{1}(1)-YMin)/DYinit)+1;
     1956                Ybound(2)=Coord{1}(1)-DYinit*(max_ind1-1);
     1957                Ybound(1)=Coord{1}(1)-DYinit*(min_ind1-1);
     1958            end             
     1959            if test_direct(2)==1
     1960                min_ind2=ceil((XMin-Coord{2}(1))/DXinit)+1;
     1961                max_ind2=floor((XMax-Coord{2}(1))/DXinit)+1;
     1962                Xbound(1)=Coord{2}(1)+DXinit*(min_ind2-1);
     1963                Xbound(2)=Coord{2}(1)+DXinit*(max_ind2-1);
     1964            else
     1965                min_ind2=ceil((Coord{2}(1)-XMax)/DXinit)+1;
     1966                max_ind2=floor((Coord{2}(1)-XMin)/DXinit)+1;
     1967                Xbound(2)=Coord{2}(1)+DXinit*(max_ind2-1);
     1968                Xbound(1)=Coord{2}(1)+DXinit*(min_ind2-1);
     1969            end
     1970            min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index)
     1971            min_ind2=max(min_ind2,1);
     1972            max_ind1=min(max_ind1,npy);
     1973            max_ind2=min(max_ind2,npx);
     1974            for ivar=VarIndex
     1975                VarName=FieldData.ListVarName{ivar};
     1976                ProjData.ListVarName=[ProjData.ListVarName VarName];
     1977                ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
     1978                if length(FieldData.VarAttribute)>=ivar
     1979                    ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     1980                end
     1981                eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']);
     1982            end         
     1983        else
     1984        % case with rotation and/or interpolation
     1985            if isempty(Coord_z) %2D case
     1986                [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     1987                XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image
     1988                YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi);
     1989                XIMA=(XIMA-minAX)/DXinit+1;% image index along x
     1990                YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
     1991                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
     1992                YIMA=reshape(round(YIMA),1,npX*npY);
     1993                flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     1994                if isequal(ObjectData.ProjMode,'filter')
     1995                    npx_filter=ceil(abs(DX/DAX));
     1996                    npy_filter=ceil(abs(DY/DAY));
     1997                    Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter);
     1998                    test_filter=1;
     1999                else
     2000                    test_filter=0;
     2001                end
     2002                for ivar=VarIndex
     2003                    VarName=FieldData.ListVarName{ivar} ;
     2004                    if test_interp(1) | test_interp(2)%interpolate on a regular grid         
     2005                          eval(['FieldData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
     2006                    end
     2007                    %filter the field (image) if option 'filter' is used
     2008                    if test_filter 
     2009                         Aclass=class(FieldData.A);
     2010                         eval(['FieldData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
     2011                         if ~isequal(Aclass,'double')
     2012                             eval(['FieldData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
     2013                         end
     2014                    end
     2015                    eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line             
     2016                    ind_in=find(flagin);
     2017                    ind_out=find(~flagin);
     2018                    ICOMB=(XIMA-1)*npy+YIMA;
     2019                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     2020                    vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:);
     2021                    for icolor=1:nbcolor
     2022                        vec_B(ind_out,icolor)=zeros(size(ind_out));
     2023                    end
     2024                    ProjData.ListVarName=[ProjData.ListVarName VarName];                 
     2025                    if length(FieldData.VarAttribute)>=ivar
     2026                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
     2027                    end     
     2028                    eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
     2029                end
     2030                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
     2031                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
     2032                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
     2033            else %3D case
     2034                if isequal(Theta,0) & isequal(Phi,0)       
     2035                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
     2036                    iz_sup=find(test_sup);
     2037                    iz=iz_sup(1);
     2038                    if iz>=1 & iz<=npz
     2039                        ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
     2040                        ProjData.DimValue=[ProjData.DimValue npY npX];
     2041                        for ivar=VarIndex
     2042                            VarName=FieldData.ListVarName{ivar};
     2043                            ProjData.ListVarName=[ProjData.ListVarName VarName];
     2044                            ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
     2045                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
     2046                            %TODO : do a vertical average for a thick plane
     2047                            if test_interp(2) | test_interp(3)
     2048                                eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
     2049                            end
     2050                        end
     2051                    end
     2052                else
     2053                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
     2054                    %TODO: use interp3
     2055                    return
     2056                end
     2057            end
     2058        end
     2059    end
     2060    %projection of  velocity components in the rotated coordinates
     2061    if ~isequal(Phi,0) && length(ivar_U)==1
     2062        if isempty(ivar_V)
     2063            msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
     2064            return
     2065        end
     2066        UName=FieldData.ListVarName{ivar_U};
     2067        VName=FieldData.ListVarName{ivar_V};   
     2068        eval(['ProjData.' UName  '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';'])
     2069        eval(['ProjData.' VName  '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');'])
     2070        if ~isempty(ivar_W)
     2071            WName=FieldData.ListVarName{ivar_W};
     2072            eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
     2073            eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
     2074        end
     2075        if ~isequal(Psi,0)
     2076            eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
     2077            eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
     2078        end
     2079    end
     2080end
     2081
     2082%-----------------------------------------------------------------
    15472083%transmit the global attributes
    15482084function [ProjData,errormsg]=proj_heading(FieldData,ObjectData)

  • trunk/src/read_set_object.m

    r19 r77  
    6767            end
    6868        else
    69            data.RangeZ=ZMax;
     69           data.RangeZ(2)=ZMax;
    7070           dimrange=[dimrange(1) 3];
    7171        end

  • trunk/src/set_object.m

    r74 r77  
    238238    end
    239239end
    240 % if desable_open
    241 %     set(handles.OPEN,'Visible','off')
    242 % else
    243 %     set(handles.OPEN,'Visible','on')
    244 % end
    245240if enable_plot
    246241   set(handles.PLOT,'enable','on')
    … …  
    331326        menu_proj={'inside';'outside';'mask_inside';'mask_outside'};
    332327    case 'volume'
    333         menu_proj={'none'};
     328        menu_proj={'interp';'none'};
    334329end   
    335330proj_index=get(handles.ProjMode,'Value');

  • trunk/src/uvmat.m

    r76 r77  
    43244324function MenuExportField_Callback(hObject, eventdata, handles)
    43254325global CurData
    4326 huvmat=findobj(allchild(0),'Name','uvmat');
    4327 UvData=get(huvmat,'UserData');
    4328 CurData=UvData.ProjField_1;
    4329 evalin('base','global CurData')%make CurData global in the workspace
    4330 display(['UserData of view_field :'])
     4326CurData=get(handles.uvmat,'UserData');
     4327% if isfield(UvData,'ProjField')
     4328    CurData=UvData;
     4329    evalin('base','global CurData')%make CurData global in the workspace
     4330    display(['current field :'])
     4331% end
    43314332evalin('base','CurData') %display CurData in the workspace
    43324333commandwindow;
    … …  
    47224723% ------------------------------------------------------------------
    47234724function Menuvolume_Callback(hObject, eventdata, handles)
     4725data.Style='volume';
     4726data.ProjMode='interp';%default
    47244727% set(handles.create,'Visible','on')
    47254728% set(handles.create,'Value',1)
    4726 VOLUME_Callback(hObject,eventdata,handles)
     4729% VOLUME_Callback(hObject,eventdata,handles)
     4730create_object(data,handles)
    47274731
    47284732%------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.