Let's assume I have a set of 3D points that have been sampled on a surface of the form $(x,y,f(x,y))$ where $f$ is a 2D Gaussian function: $$f(x,y)=\frac{1}{2\pi\sigma_x \sigma_y} \exp{\left(-\frac1{2} \left(\frac{(x-\mu_x)^2}{\sigma_x^2}+\frac{(y-\mu_y)^2}{\sigma_y^2}\right)\right)},$$ or, in matricial form, $$f({\bf x})=\tfrac{1}{2\pi\sqrt{|\Sigma|}} e^{\left(-\frac12 \left({\bf x} - {\bf \mu}\right)^T \Sigma^{-1} \left({\bf x} - {\bf \mu}\right)\right)},\quad{\bf x} = \begin{bmatrix} x \\ y \end{bmatrix},\quad{\bf \mu} = \begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix},\quad\Sigma=\begin{bmatrix} \sigma_x^2 && 0 \\ 0 && \sigma_y^2 \end{bmatrix}.$$ Here, $(\mu_x,\mu_y)$ indicates the center of the bell, and $\Sigma$ the way it spreads all over the plane. I'm thinking more from a geometric point of view, but, of course, those would be the mean and the covariance matrix (with uncorrelated variables) of a 2D normal distribution.
I know that, in case I wanted to recover those parameters from the set of points, I can use a nonlinear model fit. For instance, I've tried NonlinearModelFit in Mathematica and it works well.
But now, imagine that my set of points $\{(x_n,y_n,f(x_n,y_n))\}$ is rotated and/or translated to somewhere else in the 3D space. My question is the following one:
Given a set of points $\{(x_k, y_k, z_k)\}_{k=1}^N$, what method can I use to guess the parameters for the rotation ($\theta_x,\theta_y,\theta_z$), the translation ($v_x,v_y,v_z$) and the Gaussian ($\sigma_x,\sigma_y$) such that $$\begin{bmatrix} x_k \\ y_k \\ z_k \end{bmatrix} = R_{\theta_x,\theta_y,\theta_z}\begin{bmatrix} \tilde x_k \\ \tilde y_k \\ f_{\sigma_x,\sigma_y}(\tilde x_k,\tilde y_k) \end{bmatrix} + \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix}$$ for some unknown $(\tilde x_k,\tilde y_k)\in\mathbb{R}^2$?
EDIT:
This is what I came up with so far. I have to solve the following optimization problem: $$\min_{\begin{array}{c}\theta_x,\theta_y,\theta_z\\v_x,v_y,v_z\\\sigma_x,\sigma_y\end{array}} \sum_{k=1}^N \big(\tilde z_k - f_{\sigma_x,\sigma_y}(\tilde x_k,\tilde y_k)\big)^2$$ where $$\begin{bmatrix} \tilde x_k \\ \tilde y_k \\ \tilde z_k \end{bmatrix} = R_{\theta_x,\theta_y,\theta_z}^{-1} \left(\begin{bmatrix} x_k \\ y_k \\ z_k \end{bmatrix} - \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} \right) $$ However, in this way I do not have the classical form of the residuals $r_k=z_k-g(x_k,y_k)$, where $g$ depends on the parameters. Can I go further or is that the final expression I should use?
You new model is still the Gaussian but applied to the $x,y,f$ coordinates after multiplication by the unknown rotation matrix and translation by the unknown vector, which adds $6$ parameters to your equations.
Anyway, this would form a degenerate problem as two translational degrees of freedom have the same effect as $\mu_x,\mu_y$, and one rotational degree of freedom has the effect of nonzero covariance.
So you can recast as the fitting of an $8$-parameter model including