pattern matching - How can I separate a separable function

I have a separable function $f[x,y]$, and I would like to find two functions $g[x]$ and $h[y]$ with

$f[x,y]=g[x] h[y]$

where $g[x]$ doesn't depend on $y$ and $h[y]$ doesn't depend on $x$. Ideally, $g$ and $h$ should have the same magnitude, to prevent overflows/underflows. I have a hackish approach that works, but involves a lot of manual labor.

Background: $f[x,y]$ is a filter kernel I want to apply to an image, and using two separate 1d-filters is much more efficient.

My first approach was to start with $g[x]=f[x,0]$. But that doesn't work for e.g. $f[x,y]=\frac{e^{-\frac{x^2+y^2}{2 \sigma ^2}} x y}{2 \pi \sigma ^6}$

Currently, I have a function that "removes" $x$ or $y$ from $f[x,y]$ using pattern matching:

removeSymbol[f_, s_] := f //. {s^_ + a_ -> a, s^_.*a_ -> a}

but that means I have to manually adjust this pattern for different f's.

Is there a more elegant way to do this? $f[x,y]$ is usually a derivative of a gaussian, e.g.

gaussian[x_,y_] := 1/(2 π σ^2) Exp[-((x^2 + y^2)/(2 σ^2))]
f[x_,y_] := D[gaussian[x,y], x, y]

Answer

I would take a logarithmic derivative with respect to one variable - it should be then independent of the other one, then integrate it back over the first variable and exponentiate. The second function is found by plain division. Here is the code:

ClearAll[getGX];
getGX[expr_, xvar_, yvar_] :=
  With[{dlogg = D[Log[expr], xvar] // FullSimplify},
     Exp[Integrate[dlogg, xvar]] /; FreeQ[dlogg, yvar]];

Clear[getHY];

getHY[expr_, xvar_, yvar_] := FullSimplify[(#/getGX[#, xvar, yvar]) &[expr]]

A test function:

ftest[x_, y_] := (x^2 + 1)*y^3 *Exp[-x - y] 

Now,

getGX[ftest[x,y],x,y]

(* E^-x (1+x^2)  *)


getHY[ftest[x,y],x,y]

(* E^-y y^3 *)

The integration constant ambiguity translates into an ambiguity of how you split the function, since this operation is only defined up to a multiplicative constant factor by which you can multiply one function, and divide the other one.