Changeset 1137 for trunk/src

Timestamp:

Apr 24, 2024, 7:45:09 PM (10 days ago)

Author:

sommeria

Message:

tps improved by a cosine window in each subdomain. Default size of subdomains reduced for faster computation in the civ interface

Location:

trunk/src

Files:

• calc_field_tps.m (modified) (2 diffs)
• filter_tps.m (modified) (3 diffs)
• read_field.m (modified) (1 diff)
• series/civ_input.m (modified) (1 diff)
• series/civ_series.m (modified) (1 diff)
• series/merge_proj.m (modified) (5 diffs)
• tps_coeff_field.m (modified) (1 diff)
• uvmat.m (modified) (1 diff)

trunk/src/calc_field_tps.m

@@ -108 +108 @@
 
 %% loop on subdomains
+NbSubDomain
 for isub=1:NbSubDomain
     nbvec_sub=NbCentre(isub);
     check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
-    ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range
-    nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
+    ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range    
+    x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;%width of the subdomain/pi
+    y_width=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;%width of the subdomain/pi
+    CentreX=(SubRange(1,2,isub)+SubRange(1,1,isub))/2;%centre of the subdomain
+    CentreY=(SubRange(2,2,isub)+SubRange(2,1,isub))/2;
+    x_dist=(Coord_interp(ind_sel,1)-CentreX)/x_width;% relative x distance to the retangle centre*pi/2
+    y_dist=(Coord_interp(ind_sel,2)-CentreY)/y_width;% relative ydistance to the retangle centre
+    weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
+    nbval(ind_sel)=nbval(ind_sel)+weight;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     if check_grid
         EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
@@ -122 +130 @@
         switch FieldName{ilist}
             case 'vec(U,V)'
-                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
-                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.U(ind_sel)=DataOut.U(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,1);
+                DataOut.V(ind_sel)=DataOut.V(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,2);
             case 'U'
-                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
+                DataOut.U(ind_sel)=DataOut.U(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,1);
             case 'V'
-                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.V(ind_sel)=DataOut.V(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,2);
             case 'norm(U,V)'
-                U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
-                V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
+                U=DataOut.U(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,1);
+                V=DataOut.V(ind_sel)+weight.*EM *FieldVar(1:nbvec_sub+3,isub,2);
                 DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
             case 'curl(U,V)'
-                DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-weight.*EMDY *FieldVar(1:nbvec_sub+3,isub,1)+weight.*EMDX *FieldVar(1:nbvec_sub+3,isub,2);
             case 'div(U,V)'
-                DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.div(ind_sel)=DataOut.div(ind_sel)+weight.*EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
             case 'strain(U,V)'
-                DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+weight.*EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
             case 'DUDX(U,V)'
-                DataOut.DUDX(ind_sel)=DataOut.DUDX(ind_sel)+EMDX *FieldVar(1:nbvec_sub+3,isub,1);
+                DataOut.DUDX(ind_sel)=DataOut.DUDX(ind_sel)+weight.*EMDX *FieldVar(1:nbvec_sub+3,isub,1);
             case 'DUDY(U,V)'
-                DataOut.DUDY(ind_sel)=DataOut.DUDY(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1);
+                DataOut.DUDY(ind_sel)=DataOut.DUDY(ind_sel)+weight.*EMDY*FieldVar(1:nbvec_sub+3,isub,1);
             case 'DVDX(U,V)'
-                DataOut.DVDX(ind_sel)=DataOut.DVDX(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.DVDX(ind_sel)=DataOut.DVDX(ind_sel)+weight.*EMDX*FieldVar(1:nbvec_sub+3,isub,2);
             case 'DVDY(U,V)'
-                DataOut.DVDY(ind_sel)=DataOut.DVDY(ind_sel)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
+                DataOut.DVDY(ind_sel)=DataOut.DVDY(ind_sel)+weight.*EMDY *FieldVar(1:nbvec_sub+3,isub,2);
         end
     end

trunk/src/filter_tps.m

@@ -95 +95 @@
             break %  go to next subdomain
         % if too few selected vectors, increase the subrange for next iteration
-        elseif numel(ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous);
+        elseif numel(ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)
             SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
             SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
@@ -112 +112 @@
             % if no value exceeds threshold, the result is recorded
             if isequal(numel(ind_ind_sel),numel(ind_sel))
-                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
-                V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub;
+                x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
+                y_width=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;
+                x_dist=(Coord(ind_sel,1)-CentreX(isub))/x_width;% relative x distance to the retangle centre
+                y_dist=(Coord(ind_sel,2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
+                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
+                U_smooth(ind_sel)=U_smooth(ind_sel)+weight.*U_smooth_sub;
+                V_smooth(ind_sel)=V_smooth(ind_sel)+weight.*V_smooth_sub;
                 NbCentre(isub)=numel(ind_sel);
                 Coord_tps(1:NbCentre(isub),:,isub)=Coord(ind_sel,:);
                 U_tps(1:NbCentre(isub)+3,isub)=U_tps_sub;
                 V_tps(1:NbCentre(isub)+3,isub)=V_tps_sub;
-                nb_select(ind_sel)=nb_select(ind_sel)+1;
+                nb_select(ind_sel)=nb_select(ind_sel)+weight;
                 display(['tps done in subdomain # ' num2str(isub)  ' among ' num2str(NbSubDomain)])
                 break
             % if too few selected vectors, increase the subrange for next iteration
-            elseif numel(ind_ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous);
+            elseif numel(ind_ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)
                 SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
                 SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
@@ -129 +134 @@
                 [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),U(ind_sel(ind_ind_sel)),rho);
                 [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),V(ind_sel(ind_ind_sel)),rho);
-                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub;
-                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;
+                x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
+                y_width=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;
+                x_dist=(Coord(ind_sel(ind_ind_sel),1)-CentreX(isub))/x_width;% relative x distance to the retangle centre
+                y_dist=(Coord(ind_sel(ind_ind_sel),2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
+                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
+                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+weight.*U_smooth_sub;
+                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+weight.*V_smooth_sub;
                 NbCentre(isub)=numel(ind_ind_sel);
                 Coord_tps(1:NbCentre(isub),:,isub)=Coord(ind_sel(ind_ind_sel),:);
                 U_tps(1:NbCentre(isub)+3,isub)=U_tps_sub;
                 V_tps(1:NbCentre(isub)+3,isub)=V_tps_sub;
-                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
+                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+weight;
                 display(['tps redone after elimination of erratic vectors in subdomain # ' num2str(isub) ' among ' num2str(NbSubDomain)])
                 break

trunk/src/read_field.m

@@ -203 +203 @@
                         find(strcmp(ParamIn.Coord_x,VarDimName{ilist}))];
                 end
-                if ~isempty(DimOrder)
+                if numel(DimOrder)>=ndims(Field.(ListVarName{ilist}))
                     Field.(ListVarName{ilist})=permute(Field.(ListVarName{ilist}),DimOrder);
                     VarDimName{ilist}=VarDimName{ilist}(DimOrder);

trunk/src/series/civ_input.m

@@ -81 +81 @@
     set(handles.CheckThreshold,'Visible','on')
     set(handles.CheckDeformation,'Value',0)% desactivate 
+    set(handles.num_SubDomainSize(1),'String','250')
+    set(handles.num_SubDomainSize(2),'String','500')
 end
 switch Param.Action.ActionName

trunk/src/series/civ_series.m

@@ -924 +924 @@
         time_total=toc(tstart);
         disp(['ellapsed time ' num2str(time_total/60,2) ' minutes'])
-        disp(['time image reading ' num2str(time_input/60,2) ' minutes'])
-        disp(['time civ1 ' num2str(time_civ1/60,2) ' minutes'])
-        disp(['time patch1 ' num2str(time_patch1/60,2) ' minutes'])
-        disp(['time civ2 ' num2str(time_civ2/60,2) ' minutes'])
-        disp(['time patch2 ' num2str(time_patch2/60,2) ' minutes'])
-        disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2)/60,2) ' minutes'])
+        disp(['time image reading ' num2str(time_input,2) ' s'])
+        disp(['time civ1 ' num2str(time_civ1,2) ' s'])
+        disp(['time patch1 ' num2str(time_patch1,2) ' s'])
+        disp(['time civ2 ' num2str(time_civ2,2) ' s'])
+        disp(['time patch2 ' num2str(time_patch2,2) ' s'])
+        disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
     end
 end

trunk/src/series/merge_proj.m

@@ -68 +68 @@
     ParamOut.FieldTransform = 'on';%can use a transform function
     %%%%% list of possible transform functions (needed only for compilation)
-%         ListTransform={'phys','phys_polar','sub_field'};%list of possible transform functions (needed only for compilation)
-        % if 0==1 %never satisfied but trigger compilation with the appropriate transform functions
-        %     phys
-        %     phys_polar
-        %     sub_field
-        try
-            for ilist=1:numel(ListTransform)
-                eval(ListTransform{ilist})
-            end
-        end
+    if 0==1 %never satisfied but trigger compilation with the appropriate transform functions ('eval' inactive for compilation)
+        phys
+        phys_polar
+        sub_field
+    end
 
     ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions 
@@ -192 +187 @@
     currentdir=pwd;
     cd(Param.FieldTransform.TransformPath)
+    if strcmp(Param.FieldTransform.TransformName,'sub_field')
+        checksub=2;% the two into field series will be subtracted
+    end
     transform_fct=str2func(Param.FieldTransform.TransformName);
     cd (currentdir)
@@ -199 +197 @@
         end
     end 
-    checksub=nargin(transform_fct);% number of input arguments for the selected transform fct
     if checksub>2 && NbView>2
         disp_uvmat('WARNING',['only the two first input file series will be combined by ' Param.FieldTransform.TransformName],checkrun)
@@ -321 +318 @@
         
         %% transform the input field iview (e.g; phys) if requested (no transform involving two input fields at this stage)
-        checksub=0;
-        if ~isempty(transform_fct)
+        if ~isempty(transform_fct)&& checksub==0
             checksub=nargin(transform_fct);
-            if checksub>=2
+            if nargin(transform_fct)>=2
                 Data{iview}=transform_fct(Data{iview},XmlData{iview});
-            elseif checksub==1
+            else
                 Data{iview}=transform_fct(Data{iview});
             end
@@ -356 +352 @@
 
     %% merge the NbView fields 
-%     if checksub<=2
+    if checksub==0
         [MergeData,errormsg]=merge_field(Data);%concatene all the input field series by fct merge_data
-%     elseif checksub==3
-%         MergeData=transform_fct(Data{1},XmlData{1},Data{2}); %combine the two input file series
-%     else
-%         MergeData=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});%combine the two input file series with calibration parameters
-%     end
+    else
+        MergeData=transform_fct(Data{1},XmlData{1},Data{2}); %combine the two input file series
+    % else
+    %     MergeData=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});%combine the two input file series with calibration parameters
+    end
     if ~isempty(errormsg)
         disp_uvmat('ERROR',errormsg,checkrun);

trunk/src/tps_coeff_field.m

@@ -43 +43 @@
 Smoothing=0;
 end
-SubDomainNbPoint=1000; %default, estimated nbre of data source points in a subdomain used for tps
+SubDomainNbPoint=300; %default, estimated nbre of data source points in a subdomain used for tps
 if isfield(DataIn,'SubDomain')
     SubDomainNbPoint=DataIn.SubDomain;%old convention

trunk/src/uvmat.m

@@ -3640 +3640 @@
         end
         % case of input vector field, get the scalar used for vector color
-        if ~isempty(regexp(FieldName,'^vec('))
+        if ~isempty(regexp(FieldName,'^vec(','once'))
             list_code=get(handles.ColorCode,'String');% list menu fields
             index_code=get(handles.ColorCode,'Value');% selected string index