mygaussian2(x,y,z,varargin)

function [zc,varargout] = mygaussian2(x,y,z,varargin)
%MYGAUSSIAN2 makes use of the symmetry of the Gaussian weighting function
%   and use 1D convolution myconvper in each direction to filter the 2D 
%   signal z sampled on the domain [x,y] with a Gaussian filter
%   
%   ZC = MYGAUSSIAN2(X,Y,Z) returns the filtered signal ZC corresponding to
%   the original signal Z defined on [X,Y].
%   [ZC,HX] = MYGAUSSIAN2(X,Y,Z) returns also the 1D x-direction Gaussian
%   weighting function
%   [ZC,HX,HY] = MYGAUSSIAN2(X,Y,Z) returns the 1D x-direction and the 1D
%   y-direction Gaussian weighting functions.
%   ZC = MYGAUSSIAN2(X,Y,Z,LCXR) lets the user specify the cut-off length as
%   a fraction of the sample size LX = X(END,1)-Y(1,1) in the x-direction,
%   i.e., LCX = LCXR*LX. According to ISO 11562, the cut-off length LC for
%   Gaussian Profile Filter is L/5.
%   ZC = MYGAUSSIAN2(X,Y,Z,LCXR,LCYR) lets the user specify also the cut-off
%   length in the y-direction.
%
%   X and Y must be produced by;
%   [X,Y] = meshgrid(linspace(ax,bx,Nx),linspace(ay,by,Ny))
%   And size(Z) == [Ny,Nx]

%   Author(s): A. Almqvist
%   Copyright 2009- Andreas Almqvist
%   $Revision: 1$  $Date: 2009/11$
%   Homemade function

lcxr = 1/5;
lcyr = 1/5;
if length(varargin) == 1;
    lcxr = varargin{1};
    lcyr = 1/5;
elseif length(varargin) == 2;
    lcxr = varargin{1};
    lcyr = varargin{2};
end

% Create the domain of def. for the Gaussian weighting function
[Ny,Nx] = size(x);
Lx = x(1,end)-x(1,1); dx = x(1,2)-x(1,1); 
Ly = y(end,1)-y(1,1); dy = y(2,1)-y(1,1);
[xf,yf] = meshgrid(linspace(-Lx/2,Lx/2,Nx),linspace(-Ly/2,Ly/2,Ny)');

lcx   = lcxr*Lx; % Cut-off acc to ISO 11562 Gaussian Profile Filter Lx/5
lcy   = lcyr*Ly; % Cut-off acc to ISO 11562 Gaussian Profile Filter Ly/5
alpha = sqrt(log(2)/pi); % Auxiliary constant

% The Gaussian weighting function
variancex = 1/(pi*sqrt(2/log(2)))*lcx; 
variancey = 1/(pi*sqrt(2/log(2)))*lcy; 
hx = dx/(sqrt(2*pi)*variancex)*...
    exp(-((xf/variancex).^2)/2);
hy = dy/(sqrt(2*pi)*variancey)*...
    exp(-((yf/variancey).^2)/2);

% zcx = myfullrepx(myconvper(myperrepx(z),hx));
zcx = myconvper(z(:,1:end-1),hx);
zcx = [zcx,zcx(:,1)];
% zc  = myfullrepy(myconvper(myperrepy(zcx).',hy.').');
zc  = myconvper(zcx(1:end-1,:).',hy.').';
zc  = [zc;zc(1,:)];
if nargout == 2
    varargout{1} = hx;
elseif nargout == 3
    varargout{2} = hx;
    varargout{3} = hy;
end