function [dudx, dudy, dudz, hessian] = lsqgradient2(u, dx, dy, dz)
	% lsqgradient2(ux, dx, dy, dz)
	% calculates gradients in the 3D scalar field u(x,y,z)
	% according to the second-order polynomial regression technique used by
	% [1] and [2]
	% 
	% in a 5x5x5 region surrounding each point, a polynomial
	% u(x,y,z) = u0 + ux*x + uy*y + uz*z 
	%          + uxy*x*y + uxz*x*z + uyz*y*z
	%			  + uxx*x^2 + uyy*y^2 + uzz*z^2
	% is fitted to the neighbourhood
	%
	% [1] Elsinga GE, Adrian RJ, Oudheusden BW, Scarano F,
	% Threedimensional vortex organization in a high-Reynolds-number 
	% supersonic turbulent boundary layer
	% 
	% [2] Tokgoz S, Elsinga GE, Delfos R, Westerweel J, 2012
	% Spatial resolution and dissipation rate estimation in Taylor-Couette
	% flow for tomographic PIV 

	%% if dx, dy, dz not specified, set dx=dy=dz=1
	if nargin < 2
		dx			= 1;
		dy			= 1;
		dz			= 1;
	end
	
	%% get dimensions
	Nx				= size(u, 1);
	Ny				= size(u, 2);
	Nz				= size(u, 3);	
	
	%% construct kernel
	% model variables U are:
	% [u0 ux uy uz uxx uyy uzz uxy uxz uyz]
	% if u is a 5x5x5 array in the neighbourhood of a point, then the model
	% for U should be the best U that fits
	% u(:) = M * U
	xvec					= (-2 : 2) * dx;
	yvec					= (-2 : 2) * dy;
	zvec					= (-2 : 2) * dz;
	[xmat,ymat,zmat]	= ndgrid(xvec,yvec,zvec);
	xmat					= xmat(:);
	ymat					= ymat(:);
	zmat					= zmat(:);
	M						= [ones(125,1) xmat ymat zmat xmat.^2/2 ymat.^2/2 zmat.^2/2 xmat.*ymat xmat.*zmat ymat.*zmat];
	% least squares solution of u(:) = M * U for u is
	% U = M+ * u(:)
	% where M+ is the moore-penrose pseudoinverse of M
	Mplus					= pinv(M);
	K_dx					= reshape(Mplus(2 , :), [5 5 5]);
	K_dy					= reshape(Mplus(3 , :), [5 5 5]);
	K_dz					= reshape(Mplus(4 , :), [5 5 5]);
	K_dxdx				= reshape(Mplus(5 , :), [5 5 5]);
	K_dydy				= reshape(Mplus(6 , :), [5 5 5]);
	K_dzdz				= reshape(Mplus(7 , :), [5 5 5]);
	K_dxdy				= reshape(Mplus(8 , :), [5 5 5]);
	K_dxdz				= reshape(Mplus(9 , :), [5 5 5]);
	K_dydz				= reshape(Mplus(10, :), [5 5 5]);
	
	K_dxdx				= flip(flip(flip(K_dxdx,3),2),1);
	K_dydy				= flip(flip(flip(K_dydy,3),2),1);
	K_dzdz				= flip(flip(flip(K_dzdz,3),2),1);	
	K_dxdy				= flip(flip(flip(K_dxdy,3),2),1);
	K_dxdz				= flip(flip(flip(K_dxdz,3),2),1);
	K_dydz				= flip(flip(flip(K_dydz,3),2),1);
	
	%% separate kernels
	K_dx					= K_dx(end:-1:1, 1, 1);
	K_dy					= K_dy(1, end:-1:1, 1);
	K_dz					= K_dz(1, 1, end:-1:1);
	K_dxdx				= K_dxdx(:, 1, 1);
	K_dydy				= K_dydy(1, :, 1);
	K_dzdz				= K_dzdz(1, 1, :);
	K_dxdy				= K_dxdy(:, :, 1);
	K_dxdz				= K_dxdz(:, 1, :);
	K_dydz				= K_dydz(1, :, :);
	
	%% apply kernel
	u_box_x				= convn(u,			ones(5,1,1), 'same');
	u_box_y				= convn(u,			ones(1,5,1), 'same');
	u_box_z				= convn(u,			ones(1,1,5), 'same');
	u_box_xy				= convn(u_box_x,	ones(1,5,1), 'same');
	u_box_xz				= convn(u_box_x,  ones(1,1,5), 'same');
	u_box_yz				= convn(u_box_y,  ones(1,1,5), 'same');
	
	dudx					= convn(u_box_yz, K_dx, 'same');
	dudy					= convn(u_box_xz, K_dy, 'same');
	dudz					= convn(u_box_xy, K_dz, 'same');
	
	%% remove edges, as these are invalid
	xr						= [1 2 Nx-1 Nx];
	yr						= [1 2 Ny-1 Ny];
	zr						= [1 2 Nz-1 Nz];
	dudx(xr,:,:)		= nan;
	dudx(:,yr,:)		= nan;
	dudx(:,:,zr)		= nan;
	dudy(xr,:,:)		= nan;
	dudy(:,yr,:)		= nan;
	dudy(:,:,zr)		= nan;
	dudz(xr,:,:)		= nan;
	dudz(:,yr,:)		= nan;
	dudz(:,:,zr)		= nan;
	
	%% if hessian requested, also calculate this
	if nargout > 3
		hessian				= zeros([Nx Ny Nz 3 3]);
		hessian(:,:,:,1,1)= convn(u_box_yz, K_dxdx, 'same');
		hessian(:,:,:,2,2)= convn(u_box_xz, K_dydy, 'same');
		hessian(:,:,:,3,3)= convn(u_box_xy, K_dzdz, 'same');
		hessian(:,:,:,1,2)= convn(u_box_z , K_dxdy, 'same');
		hessian(:,:,:,1,3)= convn(u_box_y , K_dxdz, 'same');
		hessian(:,:,:,2,3)= convn(u_box_x , K_dydz, 'same');
		
		%% remove bad values
		hessian(xr,:,:,:)	=nan;
		hessian(:,yr,:,:)	=nan;
		hessian(:,:,zr,:)	=nan;
		hessian				= permute(hessian, [4 5 1 2 3]);
	end
end