Changeset 907 for trunk


Ignore:
Timestamp:
Jun 1, 2015, 10:10:05 PM (10 years ago)
Author:
sommeria
Message:

extract_rdvision updated for the new version and bug corrected in vector display for netcdf files

Location:
trunk/src
Files:
2 deleted
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/read_field.m

    r896 r907  
    8080    case 'civdata'% new format for civ results
    8181        [Field,ParamOut.VelType,errormsg]=read_civdata(FileName,InputField,ParamIn.VelType);
    82         if ~isempty(errormsg),errormsg=['read_civdata / ' errormsg];return,end     
     82        if ~isempty(errormsg),errormsg=['read_civdata / ' errormsg];return,end
    8383        ParamOut.CivStage=Field.CivStage;
    8484    case 'civx'% old (obsolete) format for civ results
    … …  
    105105                end
    106106                if check_colorvar(ilist)
     107                    if isempty(find(strcmp(InputField{ilist},ListVar)))
    107108                    Role{numel(ListVar)}='ancillary';% not projected with interpolation
    108109                    ProjModeRequest{numel(ListVar)}='';
     110                    end
    109111                else
    110112                    Role{numel(ListVar)}='scalar';
    … …  
    146148        elseif isfield(ParamIn,'TimeVarName')% case of reading of a single time  in a multidimensional array
    147149            [Field,var_detect,ichoice,errormsg]=nc2struct(FileName,'TimeVarName',ParamIn.TimeVarName,num,[ParamIn.Coord_x ParamIn.Coord_y ParamIn.Coord_z ListVar]);
    148             NbCoord=NbCoord+1;% adds time coordinate
     150            if numel(num)~=1
     151                NbCoord=NbCoord+1;% adds time coordinate, except if q single time hqs been selected
     152            end
    149153        else
    150154            [Field,var_detect,ichoice,errormsg]=nc2struct(FileName,[ParamIn.Coord_x ParamIn.Coord_y ParamIn.Coord_z ListVar]);
    … …  
    168172        VName='';
    169173        if numel(Field.ListVarName)>NbCoord % if there are variables beyond coord (1 D plots)
    170         for ilist=1:numel(ListVar)
    171             Field.VarAttribute{ilist+NbCoord}.Role=Role{ilist};
    172             Field.VarAttribute{ilist+NbCoord}.ProjModeRequest=ProjModeRequest{ilist};
    173             if isfield(ParamIn,'FieldName')
    174                 Field.VarAttribute{ilist+NbCoord}.FieldName=ListInputField{ilist};
    175             end
    176             r=regexp(ListInputField{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
    177             if ~isempty(r)&& strcmp(r.Operator,'norm')
    178                 NormName='norm';
    179                 if ~isempty(find(strcmp(ListVar,'norm')))
    180                     NormName='norm_1';
    181                 end
    182                 Field.ListVarName=[Field.ListVarName {NormName}];
    183                 ilistmax=numel(Field.ListVarName);
    184                 Field.VarDimName{ilistmax}=Field.VarDimName{ilist+2};
    185                 Field.VarAttribute{ilistmax}.Role='scalar';
    186                 Field.(NormName)=Field.(r.UName).*Field.(r.UName)+Field.(r.VName).*Field.(r.VName);
    187                 Field.(NormName)=sqrt(Field.(NormName));
    188                 UName=r.UName;
    189                 VName=r.VName;
    190             end
    191         end
    192 
    193         if ~isempty(NormName)% remove U and V if norm has been calculated and U and V are not needed as variables
    194             ind_var_U=find(strcmp(UName,ListVar));%check previous listing of variable r.UName
    195             ind_var_V=find(strcmp(VName,ListVar));%check previous listing of variable r.VName
    196             if ~checkU && ~checkV
    197                 Field.ListVarName([ind_var_U+2 ind_var_V+2])=[];
    198                 Field.VarDimName([ind_var_U+2 ind_var_V+2])=[];
    199                 Field.VarAttribute([ind_var_U+2 ind_var_V+2])=[];
    200             elseif ~checkU
    201                 Field.ListVarName(ind_var_U+2)=[];
    202                 Field.VarDimName(ind_var_U+2)=[];
    203                 Field.VarAttribute(ind_var_U+2 )=[];
    204             elseif ~checkV
    205                 Field.ListVarName(ind_var_V+2)=[];
    206                 Field.VarDimName(ind_var_V+2)=[];
    207                 Field.VarAttribute(ind_var_V+2 )=[];
    208             end
    209         end
    210         % insert coordinates as indices in case of plots vs matrix index
    211         if isfield(ParamIn,'CheckCoordIndex') && ParamIn.CheckCoordIndex
    212             Field.ListVarName=[Field.ListDimName Field.ListVarName];
    213             Field.VarDimName=[Field.ListDimName Field.VarDimName];
    214             for idim=1:numel(Field.ListDimName)
    215                 CoordName=Field.ListDimName{idim};
    216                 Field.(CoordName)=1:Field.DimValue(idim);
    217             end
    218             Field.VarAttribute=[cell(1,numel(Field.ListDimName)) Field.VarAttribute]
    219         end
    220                 end
     174            for ilist=1:numel(ListVar)
     175                Field.VarAttribute{ilist+NbCoord}.Role=Role{ilist};
     176                Field.VarAttribute{ilist+NbCoord}.ProjModeRequest=ProjModeRequest{ilist};
     177                if isfield(ParamIn,'FieldName')
     178                    Field.VarAttribute{ilist+NbCoord}.FieldName=ListInputField{ilist};
     179                end
     180                r=regexp(ListInputField{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
     181                if ~isempty(r)&& strcmp(r.Operator,'norm')
     182                    NormName='norm';
     183                    if ~isempty(find(strcmp(ListVar,'norm')))
     184                        NormName='norm_1';
     185                    end
     186                    Field.ListVarName=[Field.ListVarName {NormName}];
     187                    ilistmax=numel(Field.ListVarName);
     188                    Field.VarDimName{ilistmax}=Field.VarDimName{ilist+2};
     189                    Field.VarAttribute{ilistmax}.Role='scalar';
     190                    Field.(NormName)=Field.(r.UName).*Field.(r.UName)+Field.(r.VName).*Field.(r.VName);
     191                    Field.(NormName)=sqrt(Field.(NormName));
     192                    UName=r.UName;
     193                    VName=r.VName;
     194                end
     195            end
     196           
     197            if ~isempty(NormName)% remove U and V if norm has been calculated and U and V are not needed as variables
     198                ind_var_U=find(strcmp(UName,ListVar));%check previous listing of variable r.UName
     199                ind_var_V=find(strcmp(VName,ListVar));%check previous listing of variable r.VName
     200                if ~checkU && ~checkV
     201                    Field.ListVarName([ind_var_U+2 ind_var_V+2])=[];
     202                    Field.VarDimName([ind_var_U+2 ind_var_V+2])=[];
     203                    Field.VarAttribute([ind_var_U+2 ind_var_V+2])=[];
     204                elseif ~checkU
     205                    Field.ListVarName(ind_var_U+2)=[];
     206                    Field.VarDimName(ind_var_U+2)=[];
     207                    Field.VarAttribute(ind_var_U+2 )=[];
     208                elseif ~checkV
     209                    Field.ListVarName(ind_var_V+2)=[];
     210                    Field.VarDimName(ind_var_V+2)=[];
     211                    Field.VarAttribute(ind_var_V+2 )=[];
     212                end
     213            end
     214            % insert coordinates as indices in case of plots vs matrix index
     215            if isfield(ParamIn,'CheckCoordIndex') && ParamIn.CheckCoordIndex
     216                Field.ListVarName=[Field.ListDimName Field.ListVarName];
     217                Field.VarDimName=[Field.ListDimName Field.VarDimName];
     218                for idim=1:numel(Field.ListDimName)
     219                    CoordName=Field.ListDimName{idim};
     220                    Field.(CoordName)=1:Field.DimValue(idim);
     221                end
     222                Field.VarAttribute=[cell(1,numel(Field.ListDimName)) Field.VarAttribute];
     223            end
     224        end
    221225    case 'video'
    222226        if strcmp(class(ParamIn),'VideoReader')
    … …  
    239243        A=permute(A,[3 2 1]);
    240244    case 'multimage'
    241       %  warning 'off'
     245        %  warning 'off'
    242246        A=imread(FileName,num);
    243247    case 'image'

  • trunk/src/series/extract_rdvision.m

    r887 r907  
    153153t=xmltree;
    154154%%% A REMETTREE %%%%%%%%%%%%%%%%%%%%%
    155 %save(t,fullfile(RootPath,'Running.xml'))%create an xml file to indicate that processing takes place
     155save(t,fullfile(RootPath,'Running.xml'))%create an xml file to indicate that processing takes place
    156156
    157157%% calibration data and timing: read the ImaDoc files
    … …  
    166166%%%  loop on the cameras ( #iview)
    167167%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     168% RootPath=Param.InputTable(:,1);
     169% RootFile=Param.InputTable(:,3);
     170% SubDir=Param.InputTable(:,2);
     171% NomType=Param.InputTable(:,4);
     172% FileExt=Param.InputTable(:,5);
     173
     174% [XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
     175% if size(time,1)>1
     176%     diff_time=max(max(diff(time)));
     177%     if diff_time>0
     178%         disp_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)],checkrun)
     179%     end
     180% end
     181%
     182%      nbfield2=size(time,1);
     183checkpreserve=0;% if =1, will npreserve the original images, else it erases them at the end
    168184for iview=1:size(Param.InputTable,1)
    169     filexml=[fullfile(RootPath,Param.InputTable{iview,3}) '.xml'];%new convention: xml at the level of the image folder
     185    filexml=[fullfile(RootPath,Param.InputTable{iview,2},Param.InputTable{iview,3}) '.xml'];%new convention: xml at the level of the image folder
    170186    if ~exist(filexml,'file')
    171187        disp_uvmat('ERROR',[filexml ' missing'],checkrun)
    172188        return
    173189    end
    174     [XmlData,error]=imadoc2struct_special(filexml);
    175     if isfield(XmlData,'Time')
    176         itime=itime+1;
    177         timecell{itime}=XmlData.Time;
    178     end
    179     if isfield(XmlData,'GeometryCalib') && isfield(XmlData.GeometryCalib,'SliceCoord')
    180         NbSlice_calib{1}=size(XmlData.GeometryCalib.SliceCoord,1);%nbre of slices for Zindex in phys transform
    181         if ~isequal(NbSlice_calib{1},NbSlice_calib{1})
    182             msgbox_uvmat('WARNING','inconsistent number of Z indices for the two field series');
    183         end
    184     end
     190    [XmlData,errormsg]=imadoc2struct(filexml);
    185191   
     192    newxml=[fullfile(RootPath,Param.InputTable{iview,3}) '.xml']
    186193   
    187     % correction to RDvision xml file
    188     t=xmltree(filexml);
    189    
    190     % correct Dtj and Dtk
    191     NomTypeNew='_1_1';% new file nomencalture by default
    192     ImageName='img_1_1.png';% first image name
    193     if isfield(XmlData,'NbDtj')
    194         uid_NbDtj=find(t,'ImaDoc/Camera/BurstTiming/NbDtj');
    195         uid_value=children(t,uid_NbDtj);
    196         if ~isempty(uid_value)
    197             t=set(t,uid_value(1),'value',num2str(XmlData.NbDtj));
    198         end
    199     end
    200     if isfield(XmlData,'NbDtk')
    201         uid_NbDtk=find(t,'ImaDoc/Camera/BurstTiming/NbDtk');
    202         uid_value=children(t,uid_NbDtk);
    203         if ~isempty(uid_value)
    204             t=set(t,uid_value(1),'value',num2str(XmlData.NbDtk));
    205         end
    206     end
    207     if isempty(j1_series{1}) && isfield(XmlData,'NbDti')
    208         uid_Dti=find(t,'ImaDoc/Camera/BurstTiming/Dti');
    209         t=add(t,uid_Dti,'chardata',num2str(XmlData.Dti));
    210         uid_NbDti=find(t,'ImaDoc/Camera/BurstTiming/NbDti');
    211         t=add(t,uid_NbDti,'chardata',num2str(XmlData.NbDti));
    212         uid_NbDtj=find(t,'ImaDoc/Camera/BurstTiming/NbDtj');
    213         uid_NbDtk=find(t,'ImaDoc/Camera/BurstTiming/NbDtk');
    214         t=delete(t,uid_NbDtj);
    215         t=delete(t,uid_NbDtk);
    216         uid_Dtj=find(t,'ImaDoc/Camera/BurstTiming/Dtj');
    217         uid_Dtk=find(t,'ImaDoc/Camera/BurstTiming/Dtk');
    218         t=delete(t,uid_Dtj);
    219         t=delete(t,uid_Dtk);
     194    [success,errormsg] = copyfile(filexml,newxml); %copy the xml file in the upper folder
     195       
     196    nbfield2=size(XmlData.Time,2)-1;
     197    if nbfield2>1
     198        NomTypeNew='_1_1';
     199    else
    220200        NomTypeNew='_1';
    221         ImageName='img_1.png';
    222     end
    223    
    224     %update information of 'Heading'
    225     uid_Heading=find(t,'ImaDoc/Heading');
    226     if isempty(uid_Heading)
    227         [t,uid_Heading]=add(t,1,'element','Heading');
    228     end
    229     uid_SubCampaign=find(t,'ImaDoc/Heading/SubCampaign');
    230     if ~isempty(uid_SubCampaign), t=delete(t,uid_SubCampaign); end
    231     uid_Experiment=find(t,'ImaDoc/Heading/Experiment');
    232     if ~isempty(uid_Experiment), t=delete(t,uid_Experiment); end
    233     uid_Device=find(t,'ImaDoc/Heading/Device');
    234     if ~isempty(uid_Device), t=delete(t,uid_Device); end
    235     uid_Record=find(t,'ImaDoc/Heading/Record');
    236     if ~isempty(uid_Record), t=delete(t,uid_Record); end
    237     uid_DateExp=find(t,'ImaDoc/Heading/DateExp');
    238     if ~isempty(uid_DateExp), t=delete(t,uid_DateExp); end
    239    
    240     %indicate the name of the first image (as a check that the xml file is not moved)
    241     uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
    242     if isempty(uid_ImageName)
    243         [t,uid_ImageName]=add(t,uid_Heading,'element','ImageName');
    244     end
    245     uid_value=children(t,uid_ImageName);
    246     if isempty(uid_value)
    247         t=add(t,uid_ImageName,'chardata',ImageName);%indicate  name of the first image, with ;png extension
    248     else
    249         t=set(t,uid_value(1),'value',ImageName);%indicate  name of the first image, with ;png extension
    250     end
    251    
    252     %indicate the date and time of the image acquisition start
    253     % if isfield(FileInfo,'binrepertoire') && isfield(FileInfo,'starttime')
    254     %     sep_pos=regexp(FileInfo.binrepertoire,'T');
    255     %     DateTime=FileInfo.starttime;
    256     %     if ~isempty(sep_pos)
    257     %         DateTime=[FileInfo.binrepertoire(1:sep_pos-1) ' ' DateTime];
    258     %     end
    259     %     uid_DateTime=find(t,'ImaDoc/Heading/DateTime');
    260     %     if isempty(uid_DateTime)
    261     %         [t,uid_DateTime]=add(t,uid_Heading,'element','DateTime');
    262     %     end
    263     %     uid_value=children(t,uid_DateTime);
    264     %     if isempty(uid_value)
    265     %         t=add(t,uid_DateTime,'chardata',DateTime);%indicate  name of the first image, with ;png extension
    266     %     else
    267     %         t=set(t,uid_value(1),'value',DateTime);%indicate  name of the first image, with ;png extension
    268     %     end
    269     % end
    270    
    271     %% backup the previous xml file and save the corrected one
    272 %    [success,message]=copyfile(filexml,[filexml '~']);%make backup
    273 %     if success~=1
    274 %         disp(['errror in xml file backup: ' message]);
    275 %         return
    276 %     end
    277     save(t,filexml)
    278     nbfield2=1;
    279     if isfield(XmlData,'Time')
    280         nbfield2=size(XmlData.Time,2);
    281     end
    282    
     201    end
    283202    %% get the names of .seq and .sqb files
    284203    switch Param.InputTable{iview,5}
    … …  
    286205            filename_seq=fullfile(RootPath,Param.InputTable{iview,2},[Param.InputTable{iview,3} '.seq']);
    287206            filename_sqb=fullfile(RootPath,Param.InputTable{iview,2},[Param.InputTable{iview,3} '.sqb']);
     207            [success,errormsg] = copyfile(filename_seq,[fullfile(RootPath,Param.InputTable{iview,3}) '.seq']); %copy the seq file in the upper folder
     208            [success,errormsg] = copyfile(filename_sqb,[fullfile(RootPath,Param.InputTable{iview,3}) '.sqb']); %copy the sqb file in the upper folder
    288209        otherwise
    289210            errormsg='input file extension must be .seq or .sqb';
    … …  
    291212    if ~exist(filename_seq,'file')
    292213        errormsg=[filename_seq ' does not exist'];
     214    end
     215    if ~isempty(errormsg)
     216        disp_uvmat('ERRROR',errormsg,checkrun);
    293217        return
    294218    end
    295    
     219   
     220
    296221    %% get data from .seq file
    297222    s=ini2struct(filename_seq);
    … …  
    307232        SeqData.binrepertoire=[SeqData.binrepertoire DirExt];
    308233    end
    309 %     PathDir=fileparts(PathDir);
     234   
     235    %% checking consistency with the xml file
     236    [npi,npj]=size(XmlData.Time);
     237    if ~isequal(SeqData.nb_frames,(npi-1)*(npj-1))
     238        disp_uvmat('ERRROR','inconsistent number of images with respect to the xml file',checkrun);
     239        return
     240    end
    310241   
    311242    %% reading the .sqb file
    … …  
    317248   
    318249    %%%%%%%BRICOLAGE in case of unreadable .sqb file: remplace lecture du fichier
    319 %         ind=[111 114:211];%indices of bin files
    320 %         w=1024;%w=width of images in pixels
    321 %         h=1024;%h=height of images in pixels
    322 %         bpp=2;% nbre of bytes per pixel
    323 %         lengthimage=w*h*bpp;% lengthof an image record on the binary file
    324 %         nbimages=32; %nbre of images of each camera in a bin file
    325 %         for ii=1:32*numel(ind)
    326 %             data(ii).offset=mod(ii-1,32)*2*lengthimage+lengthimage;%Dalsa_2
    327 %             %data(ii).offset=mod(ii-1,32)*2*lengthimage;%Dalsa_1
    328 %             data(ii).file_idx=ind(ceil(ii/32));
    329 %             data(ii).timestamp=0.2*(ii-1);
    330 %         end
    331 %         m.Data=data;
     250    %         ind=[111 114:211];%indices of bin files
     251    %         w=1024;%w=width of images in pixels
     252    %         h=1024;%h=height of images in pixels
     253    %         bpp=2;% nbre of bytes per pixel
     254    %         lengthimage=w*h*bpp;% lengthof an image record on the binary file
     255    %         nbimages=32; %nbre of images of each camera in a bin file
     256    %         for ii=1:32*numel(ind)
     257    %             data(ii).offset=mod(ii-1,32)*2*lengthimage+lengthimage;%Dalsa_2
     258    %             %data(ii).offset=mod(ii-1,32)*2*lengthimage;%Dalsa_1
     259    %             data(ii).file_idx=ind(ceil(ii/32));
     260    %             data(ii).timestamp=0.2*(ii-1);
     261    %         end
     262    %         m.Data=data;
    332263    %%%%%%%
    333    
     264    timestamp=zeros(1,numel(m.Data));
    334265    for ii=1: numel(m.Data)
    335266        timestamp(ii)=m.Data(ii).timestamp;
    336     end
    337     %timestamp %todo: check withDt from the xml file
    338     [BinSize,errormsg]=binread_rdv_series(RootPath,SeqData,m.Data,nbfield2,NomTypeNew)
     267        j1=1;
     268        if ~isequal(nbfield2,1)
     269            j1=mod(ii-1,nbfield2)+1;
     270        end
     271        i1=floor((ii-1)/nbfield2)+1;
     272        diff_time(i1,j1)= timestamp(ii)-XmlData.Time(i1+1,j1+1);
     273    end
     274    time_diff_max=max(diff_time');
     275    time_diff_min=min(diff_time');
     276    if max(time_diff_max)>0.005
     277        disp_uvmat('WARNING',['timestamps exceeds xml time by' num2str(max(time_diff_max))],checkrun)
     278        checkpreserve=1;
     279    end
     280    if min(time_diff_min)<-0.005
     281        disp_uvmat('WARNING',['timestamps is lower than xml time by' num2str(min(time_diff_min))],checkrun)
     282        checkpreserve=1;
     283    end
     284    if checkpreserve
     285          disp(  'max and min of timestamp-xml time for each index i:')
     286          disp(time_diff_max)
     287          disp(time_diff_min)           
     288    end
     289    [BinList,errormsg]=binread_rdv_series(RootPath,SeqData,m.Data,nbfield2,NomTypeNew);
    339290    if ~isempty(errormsg)
    340291        disp_uvmat('ERROR',errormsg,checkrun)
    341292        return
    342293    end
    343 end
    344 delete(fullfile(RootPath,'Running.xml'))%delete the  xml file to indicate that processing is finished
     294    % check the existence of the expected output image files (from the xml)
     295    for i1=1:npi-1
     296        for j1=1:npj-1
     297            OutputFile=fullfile_uvmat(RootPath,SeqData.sequencename,'img','.png',NomTypeNew,i1,[],j1);% TODO: set NomTypeNew from SeqData.mode
     298            A=imread(OutputFile);% check image reading (stop if error)
     299        end
     300    end
     301   
     302    % check images
     303   
     304    delete(fullfile(RootPath,'Running.xml'))%delete the  xml file to indicate that processing is finished
     305    if ~checkpreserve
     306        for ibin=1:numel(BinList)
     307            delete(BinList{ibin})
     308        end
     309        rmdir(fullfile(RootPath,Param.InputTable{iview,2}))
     310    end
     311end
    345312
    346313%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    347314%--------- reads a series of bin files
    348315%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    349 function [BinSize,errormsg]=binread_rdv_series(PathDir,SeqData,SqbData,nbfield2,NomTypeNew)
     316function [BinList,errormsg]=binread_rdv_series(PathDir,SeqData,SqbData,nbfield2,NomTypeNew)
    350317% BINREAD_RDV Permet de lire les fichiers bin gￜnￜrￜs par Hiris ￜ partir du
    351318% fichier seq associￜ.
    … …  
    389356    end
    390357end
     358bin_file_counter=0;
    391359for ii=1:SeqData.nb_frames
    392360    j1=[];
    … …  
    406374            [fid,msg]=fopen(fname,'rb');
    407375            if isequal(fid,-1)
    408                 disp(['error in opening ' fname ': ' msg])
     376                errormsg=['error in opening ' fname ': ' msg];
     377                return
    409378            else
    410379                disp([fname ' opened for reading'])
     380                bin_file_counter=bin_file_counter+1;
     381                BinList{bin_file_counter}=fname;
    411382            end
    412383            fseek(fid,SqbData(ii).offset,-1);%look at the right starting place in the bin file
    … …  
    414385            BinSize(NbBinFile)=0;% strat counter for new bin file
    415386        else
    416             %             fclose(fid);%close the previous bin file
    417             %             fid=fopen(fname,'rb');% open the new bin file
    418387            fseek(fid,SqbData(ii).offset,-1);%look at the right starting place in the bin file
    419388        end
Note: See TracChangeset for help on using the changeset viewer.