As can be seen, in matrix $[\Omega remainder]$, there appear eight real unknowns ${fs,fD,\varphi s,\varphi D,\alpha ,\phi ,\chi ,\eta}$, a matrix $[Tremainderi]$, and nine complex observables, since $\Omega remainder(i,j)\u2260\Omega remainder(j,i)$. This formulation leads to a determined nonlinear equation system. Therefore, to determine the rest of the unknown parameters ${fs,fD,\varphi s,\varphi D,\alpha ,\phi ,\chi ,\eta}$ and $[Tremainderi]$ simultaneously, the nonlinear least square optimization method is implemented.^{20} In order to improve the quality of this proposed method, the initial values for orientation angle of single bounce $\kappa $ and double bounce $\eta $ are all equal to $\u2212\pi 4$ (^{10}). With each pair value ($\chi ,\eta $), the parameters ${fs,fD,\varphi s,\varphi D,\alpha ,\phi}$ are determined by employing an singular value decomposition of the $\Omega remainder$. The parameter set ${fs,fD,\varphi s,\varphi D,\alpha ,\phi ,\chi ,\eta}$ in this step is determined from condition minimum of Frobenius norm of matrix $[Tremainder]i$. We show that the parameter set ${fs,fD,fv,\varphi s,\varphi D,\varphi v,\alpha ,\phi ,\chi ,\eta}$ is equivalent to a best fit under condition that the Frobenius norm of remainder matrix becomes zero, where the estimated parameters are perfectly matched to the observations. Finally, we repeat both the above steps for each pixel in the image. The algorithm is summarized in Fig. 2. When the optimal parameters are obtained from GMBD, we can extract the surface topography phase as in Eq. (27) Display Formula
$\phi s=arg{\Omega 11(\varphi )\u2212|A|2B*},$(27)
with Display Formula$A=\Omega 12(\varphi )sin\u20092\kappa +\Omega 13(\varphi )cos\u20092\kappa cos\u20092\eta \u2009sin\u20092\kappa \u2212sin\u20092\eta \u2009cos\u20092\kappa andB=\Omega 22(\varphi )sin\u20094\kappa +2\Omega 23(\varphi )cos2\u20092\kappa cos2\u20092\eta \u2009sin\u20094\kappa \u2212sin\u20094\eta \u2009cos2\u20092\kappa .$(28)