Changeset 1168 for trunk/src/transform_field

Timestamp:

Dec 9, 2024, 7:40:14 PM (3 months ago)

Author:

sommeria

Message:

bugs repaired in civ TEST mode

Location:

trunk/src/transform_field

Files:

• FFT2_detrend.m (modified) (4 diffs)
• filter2_omitnan.m (added)
• ima_filter.m (modified) (2 diffs)

trunk/src/transform_field/FFT2_detrend.m

@@ -1 @@
-% 'FFT2_detrend': calculate the 2D spectrum of the input scalar after removing the linear trend
+% 'FFT2_detrend': calculate the 2D spectrum of the input scalar or image after removing the linear trend
 %------------------------------------------------------------------------
 %%%%  Use the general syntax for transform fields with a single input %%%%
@@ -34 +34 @@
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%
-[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);
+[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);% detect the input fields
 if ~isempty(errormsg)
     DataOut.Txt=errormsg;
@@ -42 +42 @@
 DataOut.VarDimName={};
 for ilist=1:numel(CellInfo)
-    if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured coordinates
+    if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured 2D coordinates
+        %process scalar
+        ivar=CellInfo{ilist}.VarIndex_scalar(1);
+        VarName=DataIn.ListVarName{ivar};% name o the scalar field
+        z=double(DataIn.(VarName));% transform integer input to double real
+        z(isnan(z))=0;% set NaN values to 0
+        [npy,npx]=size(z);
+        
         %process coordinates
         CoordName=DataIn.ListVarName(CellInfo{ilist}.CoordIndex);
-        x1 = DataIn.(CoordName{2}); y1 = DataIn.(CoordName{1});
-        [x y] = meshgrid(x1,y1);
-        coeff(1,1) = sum(sum(x.^2)); coeff(1,2) = sum(sum(x.*y)); coeff(1,3) = sum(sum(x));
-        coeff(2,1) = sum(sum(x.*y)); coeff(2,2) = sum(sum(y.^2)); coeff(2,3) = sum(sum(y));
-        coeff(3,1) = sum(sum(x)); coeff(3,2) = sum(sum(y)); coeff(3,3) = length(x1)*length(y1);
+        x1 = linspace(DataIn.(CoordName{2})(1),DataIn.(CoordName{2})(end),npx);
+        y1 =linspace(DataIn.(CoordName{1})(1),DataIn.(CoordName{1})(end),npy);
+        [x,y] = meshgrid(x1,y1);% 2D grid of coordinates
+        
+        % prepare the grid of wave vectors
         delta_x = x1(2) - x1(1); delta_y = y1(2) - y1(1);
         Nx = length(x1); Ny = length(y1);
@@ -55 +62 @@
         ssx = find(Nxa<Nx/2); ssy = find(Nya<Ny/2);
         Nxa(Nx-ssx+1) = -Nxa(ssx)+1; Nya(Ny-ssy+1) = -Nya(ssy)+1;
-        [Nxa Ix] = sort(Nxa); [Nya Iy] = sort(Nya);
+        [Nxa, Ix] = sort(Nxa); [Nya, Iy] = sort(Nya);         
         kx1 = (2*pi/delta_x/Nx)*(Nxa-1); ky1 = (2*pi/delta_y/Ny)*(Nya-1);
         ss = find(ky1>=0); ky1 = ky1(ss);
-        [kx ky] = meshgrid(kx1,ky1);
+        [kx,ky] = meshgrid(kx1,ky1);
         DataOut.(CoordName{2}) = kx1; DataOut.(CoordName{1}) = ky1;
         if isfield(DataIn,'CoordUnit')
-            DataOut.CoordUnit=[DataIn.CoordUnit '^{-1}'];
+            DataOut.CoordUnit=[DataIn.CoordUnit '^{-1}'];%unit of wave vectors
             DataOut.ListGlobalAttribute={'CoordUnit'};
         end
-        %process scalar
-        ivar=CellInfo{ilist}.VarIndex_scalar(1);
-        VarName=DataIn.ListVarName{ivar};
-        z=DataIn.(VarName);
-        z(isnan(z))=0;% set to 0 NaN values
+        
+        % get the coeff to calulate the linear trend
+        coeff(1,1) = sum(sum(x.^2)); coeff(1,2) = sum(sum(x.*y)); coeff(1,3) = sum(sum(x));
+        coeff(2,1) = sum(sum(x.*y)); coeff(2,2) = sum(sum(y.^2)); coeff(2,3) = sum(sum(y));
+        coeff(3,1) = sum(sum(x)); coeff(3,2) = sum(sum(y)); coeff(3,3) = length(x1)*length(y1);
         rhs(1) = sum(sum(x.*z)); rhs(2) = sum(sum(y.*z)); rhs(3) = sum(sum(z));
-        lin_coeff = inv(coeff)*rhs';
+       % lin_coeff = inv(coeff)*rhs';
+       lin_coeff = coeff\rhs';
         lin_trend = lin_coeff(1)*x + lin_coeff(2)*y + lin_coeff(3);
-        z2 = z - lin_trend;
-        spec2 = abs(fft2(z2)).^2;
+        z2 = z - lin_trend;% substract the linear trend to the input field
+        
+        spec2 = abs(fft2(z2)).^2;% get spectrum as squared of fft modulus
         spec2 = spec2(Iy,Ix);
         spec2 = spec2(ss,:);
-        DataOut.(VarName) = log10(spec2);
+        DataOut.(VarName) = log10(spec2);% take the log10 of spectrum
         %DataOut.(VarName) = spec2;
-        spec_sum=sum(sum(spec2));
-        kx_mean=sum(sum(spec2.*kx))/spec_sum;
-        ky_mean=sum(sum(spec2.*ky))/spec_sum;
-        theta=atand(ky_mean/kx_mean);
-        lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean));
+%         spec_sum=sum(sum(spec2));
+%         kx_mean=sum(sum(spec2.*kx))/spec_sum;
+%         ky_mean=sum(sum(spec2.*ky))/spec_sum;
+%         theta=atand(ky_mean/kx_mean);
+%         lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean));
         
         DataOut.ListVarName=[CoordName {VarName}];%list of variables

trunk/src/transform_field/ima_filter.m

@@ -39 +39 @@
     dlg_title = 'get the filter size in x and y';
     num_lines= 2;
-    def     = { '20';'20'};
+    def     = { '21';'21'};
     if isfield(Param,'TransformInput')&&isfield(Param.TransformInput,'FilterBoxSize_x')&&...
             isfield(Param.TransformInput,'FilterBoxSize_y')
@@ -74 +74 @@
             else
                 DataOut.(VarName)=filter2(Mfiltre,DataIn.(VarName));
+               %DataOut.(VarName)=filter2_uvmat(Param.TransformInput.FilterBoxSize_x,Param.TransformInput.FilterBoxSize_y,DataIn.(VarName));
                 DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format
             end