Skip to content
Snippets Groups Projects
Select Git revision
  • 143c1ed7dbb8c60ae22be0ac5a9fabd04329a53d
  • master default protected
  • updates-2025
  • updates-2023
  • stereo
5 results

lsqgradient2d.m

Blame
  • lsqgradient2d.m 1.52 KiB
    function [dFdx, dFdy] = lsqgradient2d(F, dx, dy)
    	%% [dFdx, dFdy] = lsqgradient(F, dx, dy)
    	% computes numerical gradient based on second order,
    	% least squares approximation as detailed in PIV: A Practical Guide
    	%
    	% dF/dx_i ~=  2 f_(i+2) + f_(i+1) - f(i-1) - 2f_(i-2)
    	%              ------------------------------------
    	%                          10 dx_i
    	%
    	% F must be a 2D array of size Nx x Ny
    	% dFdx_i an Nx x Ny array
    	% at the edges,
    	% dFdx_i is calculated with either first order central 
    	% differences or at the very edges, forwards/backwards differences 
    	% 
    
    	%% parse input
    	if nargin < 2
    		dx = 1;
    		dy = 1;
    	elseif nargin < 3
    		dy = 1;
    	end
    	
    	%% allocate memory for dFdx
    	Nx		= size(F, 1);
    	Ny    = size(F, 2);
    	dFdx	= nan(size(F), class(F));
    	dFdy  = nan(size(F), class(F));
    	
    	%% forwards and backwards difference at edges
    	dFdx(1, :)		= (F(2,:)	- F(1,:)		) / dx;
    	dFdy(:, 1)		= (F(:,2)	- F(:,1)		) / dy;
    	dFdx(end,:)		= (F(end,:) - F(end-1,:)) / dx; 
    	dFdy(:,end)		= (F(:,end)	- F(:,end-1)) / dy;
    	
    	%% first order central difference near edges
    	if Nx >= 3
    		for i = [2 Nx-1]
    			dFdx(i, :) = (F(i+1,:) - F(i-1,:)) / (2*dx);
    		end
    	end
    	if Ny >= 3
    		for j = [2 Ny-1]
    			dFdy(:, j) = (F(:,j+1) - F(:,j-1)) / (2*dy);
    		end
    	end
    	
    	%% second order least squares near middle
    	% use an expansion of f(x)' which is accurate to O(dx^4)
    	for i = 3:(Nx-2)
    		dFdx(i, :) = (2*F(i+2,:) + F(i+1,:) - F(i-1,:) - 2*F(i-2,:)) / (10*dx);
    	end
    	for j = 3:(Ny-2)
    		dFdy(:, j) = (2*F(:,j+2) + F(:,j+1) - F(:,j-1) - 2*F(:,j-2)) / (10*dy);
    	end
    end