source: trunk/src/read_field.m

Last change on this file was 1137, checked in by sommeria, 8 days ago

tps improved by a cosine window in each subdomain. Default size of subdomains reduced for faster computation in the civ interface
File size: 14.1 KB
Line 
1%'read_field': read the fields from files in different formats (netcdf files, images, video)
2%--------------------------------------------------------------------------
3%  function [Field,ParamOut,errormsg] = read_field(FileName,FileType,ParamIn,num)
4%
5% OUTPUT:
6% Field: matlab structure representing the field
7% ParamOut: structure representing parameters:
8%        .FieldName; field name
9%        .VelType
10%        .CivStage: stage of civx processing (=0, not Civx, =1 (civ1), =2  (fix1)....
11%        .Npx,.Npy: for images, nbre of pixels in x and y
12% errormsg: error message, ='' by default
13%
14%INPUT
15% FileName: name of the input file
16% FileType: type of file, as determined by the function get_file_info.m
17% ParamIn: movie object or Matlab structure of input parameters
18%     .FieldName: name (char string) of the input field (for Civx data)
19%     .VelType: char string giving the type of velocity data ('civ1', 'filter1', 'civ2'...)
20%     .ColorVar: variable used for vector color
21%     .Npx, .Npy: nbre of pixels along x and y (used for .vol input files)
22%     .TimeDimName: name of the dimension considered as 'time', selected index value then set by input 'num'
23% num: frame number for movies
24%
25% see also read_image.m,read_civxdata.m,read_civdata.m,
26
27%=======================================================================
28% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
29%   http://www.legi.grenoble-inp.fr
30%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
31%
32%     This file is part of the toolbox UVMAT.
33%
34%     UVMAT is free software; you can redistribute it and/or modify
35%     it under the terms of the GNU General Public License as published
36%     by the Free Software Foundation; either version 2 of the license,
37%     or (at your option) any later version.
38%
39%     UVMAT is distributed in the hope that it will be useful,
40%     but WITHOUT ANY WARRANTY; without even the implied warranty of
41%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
42%     GNU General Public License (see LICENSE.txt) for more details.
43%=======================================================================
44
45function [Field,ParamOut,errormsg] = read_field(FileName,FileType,ParamIn,num)
46%% default output and check input
47Field=[];
48if ~exist('num','var')
49    num=1;
50end
51if isempty(num)
52    num=1;
53end
54if ~exist('ParamIn','var')
55    ParamIn=[];
56end
57ParamOut=ParamIn;%default
58errormsg='';
59if isempty(regexp(FileName,'^http://'))&& ~exist(FileName,'file')
60    errormsg=['input file ' FileName ' does not exist'];
61    return
62end
63A=[];
64InputField={};
65check_colorvar=0;
66if isstruct(ParamIn)
67    if isfield(ParamIn,'FieldName')
68        if ischar(ParamIn.FieldName)
69            InputField={ParamIn.FieldName};
70        else
71            InputField= ParamIn.FieldName;
72        end
73    end
74    check_colorvar=zeros(size(InputField));
75    if isfield(ParamIn,'ColorVar')&&~isempty(ParamIn.ColorVar)
76        InputField=[ParamIn.FieldName {ParamIn.ColorVar}];
77        check_colorvar(numel(InputField))=1;
78    end
79end
80
81%% distingush different input file types
82switch FileType
83    case 'civdata'% new format for civ results
84        [Field,ParamOut.VelType,errormsg]=read_civdata(FileName,InputField,ParamIn.VelType);
85        if ~isempty(errormsg),errormsg=['read_civdata / ' errormsg];return,end
86        ParamOut.CivStage=Field.CivStage;
87    case 'pivdata_fluidimage'
88        if ~isfield(ParamIn,'VelType')
89            ParamIn.VelType='';
90        end
91        [Field,ParamOut.VelType,errormsg]=read_pivdata_fluidimage(FileName,InputField,ParamIn.VelType);
92        ParamOut.CivStage=Field.CivStage;
93    case 'civx'% old (obsolete) format for civ results
94        ParamOut.FieldName='velocity';%Civx data found, set .FieldName='velocity' by default
95        [Field,ParamOut.VelType,errormsg]=read_civxdata(FileName,InputField,ParamIn.VelType);
96        if ~isempty(errormsg),errormsg=['read_civxdata / ' errormsg];return,end
97        ParamOut.CivStage=Field.CivStage;
98    case {'netcdf','mat'}% general netcdf file (not recognized as civ)
99        ListVarName={};
100        Role={};
101        ProjModeRequest={};
102        % scan the list InputField
103        Operator=cell(1,numel(InputField));
104        InputVar=cell(1,numel(InputField));
105        for ilist=1:numel(InputField)
106            % look for input variables to read
107            r=regexp(InputField{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
108            if isempty(r)%  no operator used
109                ListVarName=[ListVarName InputField(ilist)];%append the variable name
110                if check_colorvar(ilist)% case of field used for vector color
111                    Role{numel(ListVarName)}='ancillary';% not projected with interpolation
112                    ProjModeRequest{numel(ListVarName)}='';
113                else
114                    Role{numel(ListVarName)}='scalar';
115                    ProjModeRequest{numel(ListVarName)}='interp_lin';%scalar field (requires interpolation for plot)
116                end
117                Operator{numel(ListVarName)}='';
118            else  % an operator 'vec' or 'norm' is used
119                ListVarName=[ListVarName {r.UName}]; % append the variable in the list if not previously listed
120                if  strcmp(r.Operator,'norm')
121                    if check_colorvar(ilist)
122                    Role=[Role {'ancillary'}];
123                    else
124                       Role=[Role {'scalar'}];
125                    end
126                else
127                     Role=[Role {'vector_x'}];
128                end
129                ListVarName=[ListVarName {r.VName}];% append the variable in the list if not previously listed
130                Role=[Role {'vector_y'}];
131                Operator{numel(ListVarName)-1}=r.Operator;
132                Operator{numel(ListVarName)}='';           
133                if ~check_colorvar(ilist) && strcmp(r.Operator,'norm')
134                    ProjModeRequest{numel(ListVarName)}='interp_lin';%scalar field (requires interpolation for plot)
135                    ProjModeRequest{numel(ListVarName)-1}='interp_lin';%scalar field (requires interpolation for plot)
136                else
137                    ProjModeRequest{numel(ListVarName)}='';
138                    ProjModeRequest{numel(ListVarName)-1}='';
139                end
140            end
141        end
142        if ~isfield(ParamIn,'Coord_z')
143            ParamIn.Coord_z=[];
144        end
145        NbCoord=~isempty(ParamIn.Coord_x)+~isempty(ParamIn.Coord_y)+~isempty(ParamIn.Coord_z);
146        if isfield(ParamIn,'TimeDimName')% case of reading of a single time index in a multidimensional array
147            [Field,var_detect,ichoice,errormsg]=nc2struct(FileName,'TimeDimName',ParamIn.TimeDimName,num,[ParamIn.Coord_x ParamIn.Coord_y ParamIn.Coord_z ListVarName]);
148        elseif isfield(ParamIn,'TimeVarName')% case of reading of a single time  in a multidimensional array
149            [Field,var_detect,ichoice,errormsg]=nc2struct(FileName,'TimeVarName',ParamIn.TimeVarName,num,[ParamIn.Coord_x ParamIn.Coord_y ParamIn.Coord_z ListVarName]);
150            if numel(num)~=1
151                NbCoord=NbCoord+1;% adds time coordinate, except if a single time has been selected
152            end
153        else
154            [Field,var_detect,ichoice,errormsg]=nc2struct(FileName,[ParamIn.Coord_x ParamIn.Coord_y ParamIn.Coord_z ListVarName]);
155        end
156        if ~isempty(errormsg)
157            return
158        end
159        CheckStructured=1;
160        %scan all the variables
161        NbCoord=0;
162        if ~isempty(ParamIn.Coord_x)
163            index_Coord_x=find(strcmp(ParamIn.Coord_x,Field.ListVarName));
164            Field.VarAttribute{index_Coord_x}.Role='coord_x';%
165            NbCoord=NbCoord+1;
166        end
167        if ~isempty(ParamIn.Coord_y)
168            if ischar(ParamIn.Coord_y)
169                index_Coord_y=find(strcmp(ParamIn.Coord_y,Field.ListVarName));
170                Field.VarAttribute{index_Coord_y}.Role='coord_y';%
171                NbCoord=NbCoord+1;
172            else
173                for icoord_y=1:numel(ParamIn.Coord_y)
174                    index_Coord_y=find(strcmp(ParamIn.Coord_y{icoord_y},Field.ListVarName));
175                    Field.VarAttribute{index_Coord_y}.Role='coord_y';%
176                    NbCoord=NbCoord+1;
177                end
178            end
179        end
180        NbDim=1;
181        if ~isempty(ParamIn.Coord_z)
182            index_Coord_z=find(strcmp(ParamIn.Coord_z,Field.ListVarName));
183            Field.VarAttribute{index_Coord_z}.Role='coord_z';%
184            NbCoord=NbCoord+1;
185            NbDim=3;
186        elseif ~isempty(ParamIn.FieldName)
187            NbDim=2;
188        end
189        NormName='';
190        UName='';
191        VName='';
192        if numel(Field.ListVarName)>NbCoord % if there are variables beyond coord (exclude 1 D plots)
193            VarAttribute=cell(1,numel(ListVarName));
194            for ilist=1:numel(ListVarName)
195                index_var=find(strcmp(ListVarName{ilist},Field.ListVarName),1);
196                VarDimName{ilist}=Field.VarDimName{index_var};
197                DimOrder=[];
198                if NbDim ==2
199                    DimOrder=[find(strcmp(ParamIn.Coord_y,VarDimName{ilist})) find(strcmp(ParamIn.Coord_x,VarDimName{ilist}))];
200                elseif NbDim ==3
201                    DimOrder=[find(strcmp(ParamIn.Coord_z,VarDimName{ilist}))...
202                        find(strcmp(ParamIn.Coord_y,VarDimName{ilist}))...
203                        find(strcmp(ParamIn.Coord_x,VarDimName{ilist}))];
204                end
205                if numel(DimOrder)>=ndims(Field.(ListVarName{ilist}))
206                    Field.(ListVarName{ilist})=permute(Field.(ListVarName{ilist}),DimOrder);
207                    VarDimName{ilist}=VarDimName{ilist}(DimOrder);
208                end
209                if numel(Field.VarAttribute)>=index_var
210                VarAttribute{ilist}=Field.VarAttribute{index_var};% read var attributes from input if exist
211                end
212            end
213            check_remove=false(1,numel(Field.ListVarName));
214            for ilist=1:numel(ListVarName)
215                VarAttribute{ilist}.Role=Role{ilist};
216                VarAttribute{ilist}.ProjModeRequest=ProjModeRequest{ilist};
217                if isfield(ParamIn,'FieldName')
218                    VarAttribute{ilist}.Operator=Operator{ilist};
219                end
220                if strcmp(Operator{ilist},'norm')
221                    UName=ListVarName{ilist};
222                    VName=ListVarName{ilist+1};
223                    ListVarName{ilist}='norm';
224                    Field.norm=Field.(UName).*Field.(UName)+Field.(VName).*Field.(VName);
225                    Field.norm=sqrt(Field.norm);
226                    check_remove(ilist+1)=true;
227                    VarAttribute{ilist}.Operator='';
228                end
229            end
230            ListVarName(check_remove)=[];
231            VarDimName(check_remove)=[];
232            VarAttribute(check_remove)=[];
233            Field.ListVarName=[Field.ListVarName(1:NbCoord) ListVarName];% complement the list of vqriables, which may be listed twice
234            Field.VarDimName=[Field.VarDimName(1:NbCoord) VarDimName];
235            Field.VarAttribute=[Field.VarAttribute(1:NbCoord) VarAttribute];
236        end 
237    case 'video'
238        if strcmp(class(ParamIn),'VideoReader')
239            A=read(ParamIn,num);
240        else
241            ParamOut=VideoReader(FileName);
242            A=read(ParamOut,num);
243        end
244    case 'mmreader'
245        if strcmp(class(ParamIn),'mmreader')
246            A=read(ParamIn,num);
247        else
248            ParamOut=mmreader(FileName);
249            A=read(ParamOut,num);
250        end
251    case 'vol'
252        A=imread(FileName);
253        Npz=size(A,1)/ParamIn.Npy;
254        A=reshape(A',ParamIn.Npx,ParamIn.Npy,Npz);
255        A=permute(A,[3 2 1]);
256    case 'multimage'
257        %  warning 'off'
258        A=imread(FileName,num);
259    case 'image'
260        A=imread(FileName);
261    case 'rdvision'
262        [A,FileInfo,timestamps,errormsg]=read_rdvision(FileName,num);
263    case 'image_DaVis'
264        Input=readimx(FileName);
265        if numel(Input.Frames)==1
266            num=1;
267        end
268        A=Input.Frames{num}.Components{1}.Planes{1}';
269        for ilist=1:numel(Input.Frames{1}.Attributes)
270            if strcmp(Input.Frames{1}.Attributes{ilist}.Name,'AcqTimeSeries')
271                timestamps=str2num(Input.Frames{1}.Attributes{ilist}.Value(1:end-3))/1000000;
272                break
273            end
274        end
275    case 'cine_phantom'
276        [A,FileInfo] = read_cine_phantom(FileName,num );
277    otherwise
278        errormsg=[ FileType ': invalid input file type for uvmat'];
279       
280end
281
282%% case of image
283if ~isempty(A)
284    if strcmp(FileType,'rdvision')||strcmp(FileType,'image_DaVis')
285        Field.Time=timestamps;
286    end
287    if isstruct(ParamOut)
288        ParamOut.FieldName='image';
289    end
290    Npz=1;%default
291    npxy=size(A);
292    Field.NbDim=2;%default
293    Field.AName='image';
294    Field.ListVarName={'Coord_y','Coord_x','A'}; %
295    Field.VarAttribute{1}.Unit='pixel';
296    Field.VarAttribute{2}.Unit='pixel';
297    Field.VarAttribute{1}.Role='coord_y';
298    Field.VarAttribute{2}.Role='coord_x';
299    Field.VarAttribute{3}.Role='scalar';
300    if ndims(A)==3
301        if Npz==1;%color
302            Field.VarDimName={'Coord_y','Coord_x',{'Coord_y','Coord_x','rgb'}}; %
303            Field.Coord_y=[npxy(1)-0.5 0.5];
304            Field.Coord_x=[0.5 npxy(2)-0.5]; % coordinates of the first and last pixel centers
305            if isstruct(ParamOut)
306                ParamOut.Npx=npxy(2);% display image size on the interface
307                ParamOut.Npy=npxy(1);
308            end
309        else
310            Field.NbDim=3;
311            Field.ListVarName=['AZ' Field.ListVarName];
312            Field.VarDimName={'AZ','Coord_y','Coord_x',{'AZ','Coord_y','Coord_x'}};
313            Field.AZ=[npxy(1)-0.5 0.5];
314            Field.Coord_y=[npxy(2)-0.5 0.5];
315            Field.Coord_x=[0.5 npxy(3)-0.5]; % coordinates of the first and last pixel centers
316            if isstruct(ParamOut)
317                ParamOut.Npx=npxy(3);% display image size on the interface
318                ParamOut.Npy=npxy(2);
319            end
320        end
321    else
322        Field.VarDimName={'Coord_y','Coord_x',{'Coord_y','Coord_x'}}; %
323        Field.Coord_y=[npxy(1)-0.5 0.5];
324        Field.Coord_x=[0.5 npxy(2)-0.5]; % coordinates of the first and last pixel centers
325    end
326    Field.A=A;
327    Field.CoordUnit='pixel'; %used for mouse_motion
328end
329
330
331
Note: See TracBrowser for help on using the repository browser.