## 2 Materials and Methods  
The proposed IMC method comprised five main sections: initialization, movement, reflection, gathering and scattering. Figure 1 shows the flow chart of the IMC method and the specific algorithm process is as follows:  
Initialization: we assume that a beam of Gaussian laser emerges from the XOY plane of Cartesian coordinates, and the limited size beam propagates along the \(z\)-axis from the plane with z = 0. The divergence angle \(\theta_{0}\) and the radius of spot \(r\) follow the Gaussian distribution, while the azimuths of location \(\varphi\) and the azimuths of propagation \(\varphi_{0}\) are uniformly distributed over \([0,\ 2\pi]\). Hence, the initial coordinates and the initial direction of sampled packets can be obtained as:  
\[ \begin{cases}x=r\cos(\varphi)\\ y=r\sin(\varphi)\\ \quad z=0\end{cases}  \tag{1}\]  
\[ \begin{cases}u_{x}=\sin(\theta_{0})\cos(\varphi_{0})\\ u_{y}=\sin(\theta_{0})\sin(\varphi_{0})\\ \quad u_{z}=\cos(\theta_{0})\end{cases}  \tag{2}\]  
Considering the refraction of beam on the surface of the emission window, the direction of refracted light could be determined according to the Law of Snell. Movement: photon packets will travel a random distance \(l_{R}\) along the direction of transmission while the number of scattering is less than the preset maximum, which is called a free path. The sampling function of it is:  
\[l_{R}=-\ln(\varsigma)/c \tag{3}\]  
where \(\varsigma\) is uniformly distributed over \([0,\ 1]\) and \(c\) is the attenuation coefficient of scattering media, which is the sum of absorption coefficient \(a\) and scattering coefficient \(b\). To describe the inherent optical properties of medium more precisely, we solve for the coefficients \(a\), \(b\) and \(c\) according to the Mie scattering theory [31; 32]. The extinction efficiency factor \(Q_{e}\) and the scattering efficiency factor \(Q_{s}\) represent the ratio of extinction and scattering of the energy incident on the geometric cross-section of the particle, respectively. \[\left\{ \begin{array}{c}Q_{e}(x,m)=\frac{2}{x^{2}}\sum\limits_{n=1}^{\infty} {(2n+1)\text{Re}[a_{n}+b_{n}]}\\ Q_{s}(x,m)=\frac{2}{x^{2}}\sum\limits_{n=1}^{\infty}{(2n+1)\Big{(}{{\left|{a_{ n}}\right|}^{2}+{{\left|{b_{n}}\right|}^{2}}\Big{)}}}\end{array} \right. \tag{4}\]  
Factors \(a_{n}\) and \(b_{n}\) are defined by:  
\[\left\{ \begin{array}{c}a_{n}(x,m)=\frac{{{\psi_{n}^{\prime}(mx){\psi_{n}}(x) -m{\psi_{n}}(mx){\psi_{n}^{\prime}}(x)}}}{{{\psi_{n}^{\prime}(mx){\psi_{n}}( x)-m{\psi_{n}}(mx){\psi_{n}^{\prime}}(x)}}}\\ b_{n}(x,m)=\frac{{{\psi_{n}^{\prime}(mx){\psi_{n}}(x)-{\psi_{n}}(mx){\psi_{n}}( x)}}}{{{m{\psi_{n}^{\prime}(mx){\psi_{n}}(x)-{\psi_{n}}(mx){\psi_{n}^{\prime}}(x)}}}} \end{array} \right. \tag{5}\]  
where the size parameter \(x=2\pi r/\lambda\) represents the ratio of the particle perimeter to the incident wavelength and \(m\) is the refractive index. \(\xi_{n}(x)\) and \(\psi_{n}(x)\) are the Riccati-Bessel function described by spherical Bessel Function of the first kind \(J_{n+1/2}\) and half-integer orders Hankel functions of the second kind \(H_{n+1/2}^{(2)}\), respectively. \[\left\{ \begin{array}{c}{{\psi_{n}}(x)=\sqrt{\pi x/2}J_{n+1/2}(x)}\\ {{\xi_{n}}(x)=\sqrt{\pi x/2}H_{n+1/2}^{(2)}(x)}\end{array} \right. \tag{6}\]  
Figure 1: Flow chart for IMC method of underwater imaging. The extinction and scattering efficiency factor of a single scattering particle to photons are obtained, but the key to characterizing the extinction and scattering specialty of the media are the extinction and scattering cross section within the unit volume of the medium, namely attenuation coefficient \(c\) and scattering coefficient \(b\). Suppose the particle concentration is uniformly distributed, and given the number of particles per unit volume \(\rho\) and the mean radius of particles \(r\), the coefficients could be expressed as:  
\[\left\{ \begin{array}{l}c=\rho Q_{e}\pi r^{2}\\ b=\rho Q_{s}\pi r^{2}\end{array} \right.