Index: trunk/src/fill_GUI.m =================================================================== --- trunk/src/fill_GUI.m (revision 1121) +++ trunk/src/fill_GUI.m (revision 1122) @@ -103,5 +103,5 @@ input_string=num2str(input_data(ibox),4); end - else + elseif ischar(input_data) input_string=input_data; end Index: trunk/src/get_file_info.m =================================================================== --- trunk/src/get_file_info.m (revision 1121) +++ trunk/src/get_file_info.m (revision 1122) @@ -58,4 +58,7 @@ FileInfo.FileType='';% input file does not exist FileInfo.FieldType=''; %default output +if ~ischar(fileinput) + return +end % check the existence (not possible for OpenDAP data) if ~isempty(regexp(fileinput,'^http://'))|| exist(fileinput,'file') Index: trunk/src/read_image.m =================================================================== --- trunk/src/read_image.m (revision 1121) +++ trunk/src/read_image.m (revision 1122) @@ -35,4 +35,7 @@ %----------------------------------------------------------------------- if ~exist('FileType','var') + FileType='image'; +end +if isempty(FileType) FileType='image'; end Index: trunk/src/series/civ_input.m =================================================================== --- trunk/src/series/civ_input.m (revision 1121) +++ trunk/src/series/civ_input.m (revision 1122) @@ -114,5 +114,5 @@ end NomTypeNc=NomTypeInput; - ind_opening=FileInfo.CivStage; + ind_opening=SeriesData.FileInfo{1}.CivStage; if isempty(regexp(NomTypeInput,'[ab|AB|-]', 'once')) set(handles.ListCompareMode,'Value',2) %mode displacement advised if the nomencalture does not involve index pairs @@ -120,5 +120,5 @@ set(handles.ListCompareMode,'Value',1) end - [Data,tild,tild,errormsg]=nc2struct(FileInfo.FileName,[]); + [Data,tild,tild,errormsg]=nc2struct(SeriesData.FileInfo{1}.FileName,[]); if ~isempty(errormsg) msgbox_uvmat('ERROR',['error in netcdf input file: ' errormsg]) Index: trunk/src/series/civ_series.m =================================================================== --- trunk/src/series/civ_series.m (revision 1121) +++ trunk/src/series/civ_series.m (revision 1122) @@ -377,4 +377,6 @@ end end + ImageName_A='';ImageName_B='';%default + VideoObject_A=[];VideoObject_B=[]; %% Civ1 % if Civ1 computation is requested @@ -650,7 +652,7 @@ par_civ2.ImageB=[]; if strcmp(Param.ActionInput.ListCompareMode,'displacement') - ImageName_A_Civ2=Param.ActionInput.RefFile; - else - ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_civ2,[],j1_civ2); + ImageName_A_Civ2=Param.ActionInput.RefFile; + else + ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_civ2,[],j1_civ2); end if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield)) @@ -665,4 +667,5 @@ [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield)); end + [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A_Civ2); par_civ2.ImageWidth=FileInfo_A.Width; par_civ2.ImageHeight=FileInfo_A.Height; Index: trunk/src/series/merge_proj_volume.m =================================================================== --- trunk/src/series/merge_proj_volume.m (revision 1121) +++ trunk/src/series/merge_proj_volume.m (revision 1122) @@ -200,5 +200,5 @@ DataVol.Time=0; %TO UPDATE DataVol.ListVarName={'coord_x','coord_y','coord_z','A'}; -DataVol.VarDimName={'coord_x','coord_y','coord_z',{'coord_z','coord_y','coord_x'}}; +DataVol.VarDimName={'coord_x','coord_y','coord_z',{'coord_y','coord_x','coord_z'}}; DataVol.VarAttribute{1}.Role='coord_x'; DataVol.VarAttribute{2}.Role='coord_y'; @@ -257,5 +257,5 @@ for index_i=1:NbField_i - AMerge=zeros(NbField_j,numel(DataVol.coord_y),numel(DataVol.coord_x)); + for index_j=1:NbField_j index_j @@ -268,79 +268,47 @@ %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% - Data=cell(1,NbView);%initiate the set Data + timeread=zeros(1,NbView); - for iview=1:NbView +% for iview=1:NbView %% reading input file(s) - [Data{iview},tild,errormsg] = read_field(filecell{iview,index_j},FileType{iview},ParamIn{iview},frame_index{iview}(index_j)); + [Data,tild,errormsg] = read_field(filecell{1,index_j},FileType{1},ParamIn{1},frame_index{1}(index_j)); if ~isempty(errormsg) disp_uvmat('ERROR',['ERROR in merge_proj/read_field/' errormsg],checkrun) return end - ListVar=Data{iview}.ListVarName; + Data.A=Data.A(1:4:end,1:4:end);% reduce the mage size + ListVar=Data.ListVarName; + %reduce the image size for ilist=1:numel(ListVar) - Data{iview}.(ListVar{ilist})=double(Data{iview}.(ListVar{ilist}));% transform all fields in double before all operations + Data.(ListVar{ilist})=double(Data.(ListVar{ilist}));% transform all fields in double before all operations end % get the time defined in the current file if not already defined from the xml file - if ~isempty(time) && isfield(Data{iview},'Time') - timeread(iview)=Data{iview}.Time; + if ~isempty(time) && isfield(Data,'Time') + timeread(iview)=Data.Time; end - Data{iview}.ZIndex=index_j; - [X,Y]=meshgrid(DataVol.coord_x,DataVol.coord_y);%grid in physical coordinates - Data{iview}=proj_plane(Data{iview},XmlData{iview},X,Y); %project on the common x,y plane + Data.ZIndex=index_j; + [X,Y,Z]=meshgrid(DataVol.coord_x,DataVol.coord_y,DataVol.coord_z);%grid in physical coordinates + %Data{iview}=proj_plane(Data{iview},XmlData{iview},X,Y); %project on the common x,y plane + + if index_j==1 + AMerge=zeros(size(Data.A,1),size(Data.A,2),NbField_j); end - % merge the NbView fields only to merge views from several cameras) - - % [MergeData,errormsg]=merge_field(ProjData); - % if ~isempty(errormsg) - % disp_uvmat('ERROR',errormsg,checkrun); - % return - % end - AMerge(index_j,:,:)=Data{iview}.A; + AMerge(:,:,index_j)=Data.A; %%%%%%%%%%%%%%%% END LOOP FOR VOLUME SCAN %%%%%%%%%%%%%%%% end %interpolate on the vertical grid - Z=zeros(1,NbField_j); - for j_index=1:numel(DataVol.coord_y) - for i_index=1:numel(DataVol.coord_x) - for ZIndex=1:NbField_j - Z(ZIndex)=Zpos(XmlData{iview},ZIndex,X(j_index,i_index),Y(j_index,i_index)); - end - DataVol.A(:,j_index,i_index) = interp1(Z,AMerge(:,j_index,i_index),DataVol.coord_z); - end - end + DataVol.A=proj_volume(AMerge,XmlData{1},X,Y,Z); +% Z=zeros(1,NbField_j); +% for j_index=1:numel(DataVol.coord_y) +% for i_index=1:numel(DataVol.coord_x) +% for ZIndex=1:NbField_j +% Z(ZIndex)=Zpos(XmlData{iview},ZIndex,X(j_index,i_index),Y(j_index,i_index)); +% end +% DataVol.A(:,j_index,i_index) = interp1(Z,AMerge(:,j_index,i_index),DataVol.coord_z); +% end +% end error=struct2nc(OutputFile,DataVol)%save result file - - % if index==1 - % TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos','Time'}; - % TimeData.Conventions='uvmat'; - % TimeData.Project='2016_Circumpolar'; - % TimeData.CoordUnit='cm'; - % TimeData.TimeUnit='s'; - % % TimeData.ZPos=ZPos; - % TimeData.ListVarName={'radius','azimuth','U','V','curl','div'}; - % TimeData.VarDimName={'radius','azimuth',{'radius','azimuth'},{'radius','azimuth'}... - % {'radius','azimuth'},{'radius','azimuth'}}; - % TimeData.VarAttribute{1}.Role=''; - % TimeData.VarAttribute{2}.Role=''; - % TimeData.VarAttribute{3}.Role='vector_x'; - % TimeData.VarAttribute{4}.Role='vector_y'; - % TimeData.VarAttribute{5}.Role='scalar'; - % TimeData.VarAttribute{6}.Role='scalar'; - % - % end - % - % %% append data to the netcdf file for next iterations - % - % - % error=struct2nc(OutputFile,TimeData);%save result file - % if isempty(error) - % disp(['output file ' OutputFile ' written']) - % else - % disp(error) - % end - % ellapsed_time=toc(tstart); - % disp(['ellapsed time since start ' num2str(ellapsed_time/60,2) ' minutes']) end @@ -426,4 +394,46 @@ end + +function ZIndex=XYZtoIndex(XmlData,X,Y,Z)% Z positions corresponding to X,Y positions +Zp=Z-XmlData.Slice.SliceCoord(1,3); +DZ=XmlData.Slice.SliceCoord(end,3)-XmlData.Slice.SliceCoord(1,3)/(size(XmlData.Slice.SliceCoord,1)-1); +if DZ~=0 +ZIndex=(Z-XmlData.Slice.SliceCoord(1,3))/DZ+1;% Z=XmlData.SliceCoord(1,3)+(ZIndex-1)*DZ +end +% effect of angular deviation +SliceAngleMax=XmlData.Slice.SliceAngle(end,:)-XmlData.Slice.SliceAngle(1,:); +normAxe=norm(SliceAngleMax); +if normAxe>0 % case of angle scan +a=-SliceAngleMax(2)/normAxe; +b=-SliceAngleMax(1)/normAxe; +c=-a*XmlData.Slice.SliceCoord(1,1)+b*XmlData.Slice.SliceCoord(1,2);%equation of the axis ax+by+c=0 +DNormal=norm(a*X+b*Y+c); +Ang=atand(Zp./DNormal); +ZIndex=ZIndex+Ang-norm(XmlData.Slice.SliceAngle(1,:))+1; +end + +%------------------------------------------------------------ +% proj_volume: poject each image on a common grid given by coord_x and coord_y +function A=proj_volume(AMerge,XmlData,X,Y,Z) + +A=[]; %default output + +%% initial image coordinates +[npy,npx,npz]=size(AMerge); +xima=0.5:npx-0.5;%image coordinates of corners +%yima=npy-0.5:-1:0.5; +yima=0.5:npy-0.5; +zima=0.5:npz-0.5; +[XIMA_init,YIMA_init,ZIMA_init]=meshgrid(xima,yima,zima);%grid of initial image in px coordinates + +%% projected coordinates +ZIndex=XYZtoIndex(XmlData,X,Y,Z);% Z positions correspoonding to X,Y positions + +%% interpolation on the new grid +[XIMA,YIMA]=px_XYZ(XmlData.GeometryCalib,XmlData.Slice,X,Y,Z);% image coordinates for each point in the real +A=interp3(XIMA_init,YIMA_init,ZIMA_init,AMerge,XIMA,YIMA,ZIndex); + + + % proj_plane: poject each image on a common grid given by coord_x and coord_y function DataOut=proj_plane(DataIn,XmlData,X,Y) Index: trunk/src/series/sliding_average.m =================================================================== --- trunk/src/series/sliding_average.m (revision 1121) +++ trunk/src/series/sliding_average.m (revision 1122) @@ -1,5 +1,5 @@ -%'turb_correlation_time': calculate the time correlation function at each point +%'sliding_average_natalya' %------------------------------------------------------------------------ -% function ParamOut=turb_correlation_time(Param) +% function ParamOut=sliding_average(Param) % %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%% Index: trunk/src/tps_coeff_field.m =================================================================== --- trunk/src/tps_coeff_field.m (revision 1121) +++ trunk/src/tps_coeff_field.m (revision 1122) @@ -1,5 +1,5 @@ %'tps_coeff_field': calculate the thin plate spline (tps) coefficients within subdomains for a field structure %--------------------------------------------------------------------- -% DataOut=tps_coeff_field(DataIn,checkall) +% DataOut=tps_coeff_field(DataIn,checkall,Smoothing) % % OUTPUT: @@ -13,5 +13,6 @@ % checkall:=1 if tps is needed for all fields (a projection mode interp_tps has been chosen), % =0 otherwise (tps only needed to get spatial derivatives of scattered data) -% +% Smoothing parameter: =0 (no smoothing) by default + % called functions: % 'find_field_cells': analyse the input field structure, grouping the variables into 'fields' with common coordinates @@ -37,6 +38,9 @@ %======================================================================= -function [DataOut,errormsg]=tps_coeff_field(DataIn,checkall) +function [DataOut,errormsg]=tps_coeff_field(DataIn,checkall,Smoothing) DataOut=DataIn;%default +if ~exist('Smoothing','var') +Smoothing=0; +end SubDomainNbPoint=1000; %default, estimated nbre of data source points in a subdomain used for tps if isfield(DataIn,'SubDomain') @@ -53,5 +57,5 @@ nbtps=0;% indicate the number of tps coordinate sets in the field structure (in general =1) -for icell=1:numel(CellInfo); +for icell=1:numel(CellInfo) if NbDimArray(icell)>=2 && strcmp(CellInfo{icell}.CoordType,'scattered') %if the coordinates are scattered NbCoord=NbDimArray(icell);% dimension of space @@ -129,5 +133,5 @@ for isub=1:size(SubRange,3) ind_sel=IndSelSubDomain(1:NbCentre(isub),isub);% array indices selected for the subdomain - [tild,Var_tps(1:NbCentre(isub)+NbCoord+1,isub)]=tps_coeff([X(ind_sel) Y(ind_sel)],DataIn.(ListVarInterp{ilist})(ind_sel),0);%calculate the tps coeff in the subdomain + [tild,Var_tps(1:NbCentre(isub)+NbCoord+1,isub)]=tps_coeff([X(ind_sel) Y(ind_sel)],DataIn.(ListVarInterp{ilist})(ind_sel),Smoothing);%calculate the tps coeff in the subdomain end DataOut.(ListNewVar{ilist+3})=Var_tps; Index: trunk/src/transform_field/ima_filter.m =================================================================== --- trunk/src/transform_field/ima_filter.m (revision 1121) +++ trunk/src/transform_field/ima_filter.m (revision 1122) @@ -1,12 +1,15 @@ -% 'ima_filter': example of image transform with input parameters: low-pass filter of an image +% 'ima_filter': low-pass filter of an image or other 2D fields defined on a regular grid +% the size of the filtering window in x and y is interactivement defined %------------------------------------------------------------------------ %%%% Use the general syntax for transform fields with a single input and parameters %%%% -% OUTPUT: -% DataOut: output field structure +% OUTPUT: +% DataOut: output field structure % %INPUT: % DataIn: input field structure +% % Param: matlab structure whose field Param.TransformInput contains the filter parameters +% DataIn_1: variables possibly introduced as a second input field %----------------------------------- @@ -78,20 +81,20 @@ if exist('DataIn_1','var') [CellInfo,NbDim,errormsg]=find_field_cells(DataIn_1); -for icell=1:numel(CellInfo) - if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') - for ivar=1:numel(CellInfo{icell}.VarIndex) - VarName=DataIn_1.ListVarName{CellInfo{icell}.VarIndex(ivar)}; - Atype=class(DataIn_1.(VarName));% detect integer 8 or 16 bits - if numel(size(DataIn_1.(VarName)))==3 - DataOut.(VarName)=filter2(Mfiltre,sum(DataIn_1.(VarName),3));%filter the input image, after summation on the color component (for color images) - DataOut.(VarName)=uint16(DataOut.(VarName)); %transform to 16 bit images - else - DataOut.(VarName)=filter2(Mfiltre,DataIn_1.(VarName)); - DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + for icell=1:numel(CellInfo) + if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') + for ivar=1:numel(CellInfo{icell}.VarIndex) + VarName=DataIn_1.ListVarName{CellInfo{icell}.VarIndex(ivar)}; + Atype=class(DataIn_1.(VarName));% detect integer 8 or 16 bits + if numel(size(DataIn_1.(VarName)))==3 + DataOut.(VarName)=filter2(Mfiltre,sum(DataIn_1.(VarName),3));%filter the input image, after summation on the color component (for color images) + DataOut.(VarName)=uint16(DataOut.(VarName)); %transform to 16 bit images + else + DataOut.(VarName)=filter2(Mfiltre,DataIn_1.(VarName)); + DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + end end end end end -end Index: trunk/src/transform_field/ima_filter_high.m =================================================================== --- trunk/src/transform_field/ima_filter_high.m (revision 1122) +++ trunk/src/transform_field/ima_filter_high.m (revision 1122) @@ -0,0 +1,99 @@ +% 'ima_filter_high': high-pass filter of an image or other 2D fields defined on a regular grid +% the size of the filtering window in x and y is interactivement defined + +%------------------------------------------------------------------------ +%%%% Use the general syntax for transform fields with a single input and parameters %%%% +% OUTPUT: +% DataOut: output field structure +% +%INPUT: +% DataIn: input field structure +% Param: matlab structure whose field Param.TransformInput contains the filter parameters +% DataIn_1: variables possibly introduced as a second input field +%----------------------------------- + +%======================================================================= +% Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France +% http://www.legi.grenoble-inp.fr +% Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr +% +% This file is part of the toolbox UVMAT. +% +% UVMAT is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published +% by the Free Software Foundation; either version 2 of the license, +% or (at your option) any later version. +% +% UVMAT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License (see LICENSE.txt) for more details. +%======================================================================= + +function DataOut=ima_filter_high(DataIn,Param,DataIn_1) + +%% request input parameters +if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) + prompt = {'npx';'npy'}; + dlg_title = 'get the filter size in x and y'; + num_lines= 2; + def = { '20';'20'}; + if isfield(Param,'TransformInput')&&isfield(Param.TransformInput,'FilterBoxSize_x')&&... + isfield(Param.TransformInput,'FilterBoxSize_y') + def={num2str(Param.TransformInput.FilterBoxSize_x);num2str(Param.TransformInput.FilterBoxSize_y)}; + end + answer = inputdlg(prompt,dlg_title,num_lines,def); + DataOut.TransformInput.FilterBoxSize_x=str2num(answer{1}); %size of the filtering window + DataOut.TransformInput.FilterBoxSize_y=str2num(answer{2}); %size of the filtering window + return +end + +DataOut=DataIn; %default + +%definition of the cos shape matrix filter +ix=1/2-Param.TransformInput.FilterBoxSize_x/2:-1/2+Param.TransformInput.FilterBoxSize_x/2;% +iy=1/2-Param.TransformInput.FilterBoxSize_y/2:-1/2+Param.TransformInput.FilterBoxSize_y/2;% +%del=np/3; +%fct=exp(-(ix/del).^2); +fct2_x=cos(ix/((Param.TransformInput.FilterBoxSize_x-1)/2)*pi/2); +fct2_y=cos(iy/((Param.TransformInput.FilterBoxSize_y-1)/2)*pi/2); +%Mfiltre=(ones(5,5)/5^2); +Mfiltre=fct2_y'*fct2_x; +Mfiltre=Mfiltre/(sum(sum(Mfiltre)));%normalize filter + +[CellInfo,NbDim,errormsg]=find_field_cells(DataIn) +for icell=1:numel(CellInfo) + if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') + for ivar=1:numel(CellInfo{icell}.VarIndex) + VarName=DataIn.ListVarName{CellInfo{icell}.VarIndex(ivar)}; + Atype=class(DataIn.(VarName));% detect integer 8 or 16 bits + if numel(size(DataIn.(VarName)))==3 + DataOut.(VarName)=filter2(Mfiltre,sum(DataIn.(VarName),3));%filter the input image, after summation on the color component (for color images) + DataOut.(VarName)=uint16(DataOut.(VarName)); %transform to 16 bit images + else + DataOut.(VarName)=DataIn.(VarName)-filter2(Mfiltre,DataIn.(VarName)); + DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + end + end + end +end +if exist('DataIn_1','var') + [CellInfo,NbDim,errormsg]=find_field_cells(DataIn_1); +for icell=1:numel(CellInfo) + if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') + for ivar=1:numel(CellInfo{icell}.VarIndex) + VarName=DataIn_1.ListVarName{CellInfo{icell}.VarIndex(ivar)}; + Atype=class(DataIn_1.(VarName));% detect integer 8 or 16 bits + if numel(size(DataIn_1.(VarName)))==3 + DataOut.(VarName)=filter2(Mfiltre,sum(DataIn_1.(VarName),3));%filter the input image, after summation on the color component (for color images) + DataOut.(VarName)=uint16(DataOut.(VarName)); %transform to 16 bit images + else + DataOut.(VarName)=filter2(Mfiltre,DataIn_1.(VarName)); + DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + end + end + end +end +end + +