A novel multi-modal fundus image fusion method for guiding the laser surgery of central serous chorioretinopathy

1.   Introduction

It is difficult to deal with field theory problems defined in the region $W$ bounded by the non-Jordan curves. However, if a conformal mapping $\Upsilon$ exists that maps region $D$ bounded by smooth Jordan curves onto domain $W$, then we can transform the problems specified in region $W$ into domain $D$ to simply the problems. Meanwhile, we can refer to region $D$ as the preimage domain of the region $W$. Because the preimage domains have better properties, the method for calculating the preimage domain needs to be proposed and extended. Moreover, the preimage domain we discuss is given in the sense of conformal equivalence, which should be emphasized. Conformal mappings are an important theory of analytic function and potential theory, which find in the methods of solving problems in field theory. Conformal mappings can map the complex regions onto the simple domains, which simplifies the issues such as calculating the capacitance of irregular capacitors, harmonic measures ^[1], logarithmic capacity ^[2,3] in potential theory, and so on ^[4,5,6,7]. Many more practical applications of conformal mappings can be found in fluid mechanics ^[8,9] and optics ^[10]. As is known to all, the analytic expressions of conformal mappings are not available in many actual situations, only can be found in some specific regions, and conformal mappings have the superiority in dealing with the problems in field theory. Therefore, the numerical methods to compute conformal mappings have been extensively investigated in the past several decades ^{[11,12,13,14]}. In comparison, the numerical method presented by Nasser, based on the boundary integral equation with the generalized Neumann kernel, has the characteristic of a unified form and can use for thirty-nine canonical slit regions proposed by Koebe ^{[15,16,17,18]}. Using the method presented by Nasser to obtain the numerical solutions of conformal mappings for highly connected regions with piecewise smooth boundaries is quick and accurate ^[19,20].

Several fluid mechanics phenomena are included in external flows past objects. Computing the preimage domains of specified slit regions is of great significance in explaining some fluid phenomena. For example, Crowdy finds several explicit formulas drawn up according to the Schottky-Klein prime function for the problems in fluid dynamics ^{[21,22,23,24]}. According to Ref.^[24] Crowdy (2008), an explicit solution of complex potential for the irrotational incompressible potential fluid driven by an assembly of stirrers moving at specified velocities has been found. Examples presented in ^[24] show streamlines of the irrotational flow generated by circular disk stirrers and paddles, where the paddles can regard as rectilinear slits. The feasible methods to solve these problems are conformally mapping the specified slit region $W$ onto the preimage domain $D$ with the same connectivity and solving these problems in $D$. However, it is difficult to directly calculate conformal mapping that maps a specified slit region onto a preimage domain, but we can find many methods which compute conformal mapping from a preimage domain onto a slit region ^{[11,16,17,25]}. Nasser ^[26] reformulates the problems of stirring fluid as the classical Riemann-Hilbert (R-H) problems and obtains the numerical complex potentials for the flows. When the stirrers are finite rectilinear slits, Nasser ^[26] presents an iterative method, which calculates a preimage domain $D$ bounded by ellipses and conformal mapping $\Upsilon$ that maps $D$ onto the rectilinear slit region. In addition, conformal mappings are also a powerful tool for simulating the flow around bodies, i.e., the phenomenon of boundary layer separation caused by the fluid passing through the surface of objects ^[25,26]. The flow around bodies has been widely studied as a classical problem in fluid mechanics. According to Ref.^[27] Chen et al. (2021), the analytical solution of the potential flow is derived in the region bounded by a circle, an ellipse, and a thin airfoil. They find that the singular equation or the hypersingular equation can be used to solve the symmetric or antisymmetric potential flow problem alone. According to Ref.^[28] Chen et al. (2022), the analytical solution of the potential flow around two identical cylinders is derived by using the boundary integral equation with the bipolar degenerate-kernel. When considering the motion of spiral point vortices in a two-dimensional multiply connected domain with two parallel slits, Aoyama et al. ^[29] first propose an iterative method based on the particle charge simulation method, which can calculate the preimage domain of the slit region. According to Ref.^[29] Aoyama et al. (2013), many equilibria are composed of two spiral point vortices that enhance the downward vertical force and the counterclockwise rotational force.

As a result, computing the preimage domains of the specified slit regions and conformal mappings that map the preimage domains onto the slit regions are problems worthy of attention. According to Ref.^[30] Nasser (2019), the iterative numeral method has successfully computed the preimage domains $D$ and conformal mappings $\Upsilon$ that map $D$ onto $W$, where $W$ is two instances of bounded slit regions called an annulus with radial slit region and the unit disk with radial slit region. However, the thirty-nine canonical slit regions cataloged by Koebe have extensive application prospects. The logarithmic spiral slit regions deserve special attention because when the oblique angle is $\frac{\pi}{2}$ or $0$, the logarithmic spiral will degenerate to a circle with a center at the origin or a ray issuing from the origin, respectively ^{[17,25,31,32]}. Meanwhile, the logarithmic spiral slit regions are closely related to the spiral point vortices.

However, the methods proposed in ^[29,30] can not directly calculate the preimage domains of the logarithmic spiral slit regions. Therefore, we present the improved iterative methods based on the method proposed in ^[30], which calculates the preimage domains $D$ of the specified logarithmic spiral slit regions $W$ and conformal mappings $\Upsilon$ that map $D$ onto $W$. We postulate that the boundaries of the sequence $\left\lbrace D^m\right\rbrace$ of the preimage domains are the ellipses. According to the distribution of logarithmic spiral slits in $W$, the lengths of the long axes and the coordinates of the centers of ellipses, which calculate during iteration, respectively converge to the lengths of the long axes and the coordinates of the centers of the boundaries $\partial D$. In each iteration step, to calculate conformal mappings quickly and accurately, we use the method based on the boundary integral equation with the generalized Neumann kernel to calculate conformal mappings.

The rest of this paper is organized as follows. In Sect.2, we introduce the related knowledge about the boundary integral equation with the generalized Neumann kernel. In Sect.3, we regard two instances of spiral slit regions called the unbounded logarithmic spiral slit region and a unit disk with logarithmic spiral slit region. We briefly review the process of formulating mapping functions as the R-H problems, and then the iterative methods are presented for calculating the preimage domains and conformal mappings that map the preimage domains onto the spiral slit regions. In Sect.4, numerical examples are proposed showing the availability of our iterative methods and simulating the flow around bodies. The results are summed up in the last section.

2.   Preliminaries

In this Section, we introduce the related knowledge about the boundary integral equation with the generalized Neumann kernel. Further details can be found in ^[14,33].

We suppose that $D$ is a $n$-connected region in the extended complex plane $\mathbb{C}\cup\{\infty\}$ containing the point $z = 0$, where the boundary of $D$ is $\varGamma: = \partial D = \mathop{\cup}\limits_{j = 1}^n \varGamma_j.$ $\varGamma_1, \varGamma_2, \dots, \varGamma_n$ are $n$ closed smooth Jordan curves. Denote the orientation of $\varGamma$ to ensure $D$ where the left of $\varGamma$.

Suppose that the curves $\varGamma_1, \varGamma_2, \dots, \varGamma_n$ can be parameterized by the $2\pi$-periodic complex functions $\eta_1(t), \eta_2(t), \dots, \eta_n(t)$ with twice continuously derivative, respectively, in which the first derivative $d\eta_j(t)/dt = \dot{\eta}_j(t)\neq0$ for $t\in J_j = [0, 2\pi], \; j = 1, 2, \dots, n$. The parametrization of the boundary $\varGamma$ can be defined by

Definition 2.1. A real function $f$ is Hölder continuous, when there are nonnegative real constants $C$, $\alpha > 0$, such that

for all $x$ and $y$ in the domain of $f$.

We denote by $H$ the space of real functions on $\varGamma$ whose are Hölder continuous. Thus, $\phi\in H$ indicates that $\phi(t): = \phi(\eta(t))$ is a real Hölder continuous function with $2\pi$-period.

Definition 2.2. A piecewise constant function $h(t)$ is defined on $J_j$ such that

where $h_j\in \mathbb{R}, \; j = 1, 2, \dots, n.$ For convenience, we denote $h(t) = (h_0, h_1, \dots, h_n)$.

Let

is a continuously differentiable complex function. $\alpha\in D$ is a fixed point, and $\theta(t) = (\theta_1, \theta_2, \dots, \theta_n)$ with constants $\theta_j \in \mathbb{R}$ for $j = 1, 2, \dots, n$ representing the oblique angle of each slit is a specified piecewise function.

The generalized Neumann kernel $N(s, t)$ formed with $B$ and $\eta$ is defined as follows:

The singular kernel $M(s, t)$ is defined as follows:

The Fredholm integral operator ${\bf N}$ and singular integral operator ${\bf M}$ on $H$ can be defined as follows:

and

Theorem 2.1. (^[17])For a specified function $\gamma\in H$, there exists a unique real function $\mu$ and a unique piecewise constant real function $h = (h_1, h_2, \dots, h_n)$ defined on $\varGamma$ such that

are boundary values of an analytic function $g$ in $D$, which $g$ is a unique solution of R-H problem ${\rm Re}\left[Bg\right] = \gamma+h$. Then the function $\mu$ is the solution of the integral equation

and the function $h$ is obtained by

3.   The iterative methods for computing the preimage domains

In this Section, we briefly review the process of formulating mapping functions as the R-H problems, and the iterative methods are presented for calculating the preimage domains $D$ and conformal mappings that map $D$ onto $W$, which $W$ is a specified unbounded or bounded logarithmic spiral slit region (For short: spiral slit region).

In the $w$ plane, the boundaries of the spiral slit regions are composed of the slits that have the origin as the asymptotic point and satisfy

If any ray issuing from the origin intersects a slit, $\theta$ named the oblique angle of the slit is the angle between the ray with the tangent of the slit at the intersection. Especially, for the oblique angle $\theta = 0$, all points of the slit will fall on a ray, which the slit is a radial slit pointing to the origin; For the oblique angle $\theta = \frac{\pi}{2}$, the ray being always perpendicular to the tangent of any point in the slit, and the slit will be a circle slit centered at the origin^{[16,17,25,32]}.

3.1.   The preimage domain of the unbounded logarithmic spiral slit region

In this Subsection, we propose an iterative method based on the boundary integral equation, which calculates a preimage domain of a specified unbounded logarithmic spiral slit region.

Theorem 3.1. (^[25]) Let $\theta(t) = (\theta_1, \theta_2, \dots, \theta_n)$ and a fixed point $\alpha$ in the unbounded multiply connected region $D$ if $\Upsilon$ that conformally maps $D$ onto the unbounded logarithmic spiral slit region satisfies normalization conditions $\Upsilon(\alpha) = 0, \; \Upsilon(\infty) = \infty, \; {\lim\limits_{z \to \infty}}\frac{\Upsilon(z)}{z} = 1$, and then the Laurent series

is uniquely determined.

Suppose $D$ is an unbounded $n$-connected region in the $z$ plane. The unbounded spiral slit region $W$ is the entire $w$ plane with $n$ spiral slits having the oblique angles of $\theta_1, \theta_2, \dots, \theta_n$, respectively. According to the normalization conditions, conformal mapping $\Upsilon(z)$ can be defined by (see ^[17,25,32])

which the analytic function $g(z)$ defines on $D\cup\varGamma$. Conformal mapping $\omega = \Upsilon(z)$ maps $\varGamma_j$ onto the spiral slit $L_j$ for $j = 1, 2, \dots, n$ (see Figure 1). According to the boundary conditions (3.1), the boundary values of conformal mapping $\Upsilon(z)$ should satisfy

Plugging (3.2) into (3.3) leads to

where $R_j\in \mathbb{R}$ for $j = 1, 2, \dots, n$ are undetermined real constants obtained by solving the R-H problem with

According to Theorem 2.1, by solving the integral equation (2.7) and calculating (2.8), we can obtain $\mu(t)$ and $h(t) = (R_1, R_2, \dots, R_n)$ such that

which satisfie the boundary conditions (3.4). For $z\in D$, we use the following Cauchy integral to calculate the values of $g(z)$

Together with (3.2), we can obtain conformal mapping $\Upsilon(z)$.

Note that the boundaries of a specified spiral slit region are fixed. However, $R_j$ calculated by the above method are unable to specify in advance, which affect the location distributions of the boundaries $L_j$ for $j = 1, 2, \dots, n$, and conformal mapping cannot be computed directly if the preimage $D$ is unknown. Thus, the problems will become finding a preimage domain $D$ of a specified spiral slit region $W$ and calculating conformal mapping $\Upsilon(z)$ that maps $D$ onto $W$.

In the $\zeta$ plane,

are $n$ line segments with lengths $l_j^*$ and upper endpoints $c_j$ for $j = 1, 2, \dots, n$ (see Figure 1). Any point on a line segment can be represented as

Therefore, any point on a spiral slit can be written as

In the following, we will propose an iterative method based on the boundary integral equation, which calculates a preimage domain $D$ of a specified unbounded spiral slit region and derives conformal mapping $\Upsilon(z)$ from $D$ onto $W$ at the same time. Each step iteration of our iterative method will generate a preimage domain denoted by $D^m$, and the sequence $\left\lbrace D^m\right\rbrace$ will converge to our needed preimage domain $D.$

Assume that the boundary of the domain $D^m$ is composed of ellipses. In the iterating step $m = 0, 1, 2, \dots,$ the ellipses $\varGamma_j^m$ can are parametrized by

where

● $Q_j^m$ denote the coordinates of the centers of ellipses $\varGamma_j^m$ for $j = 1, 2, \dots, n$,

● $a_j^m$ represent half the lengths of the long axes of elliptical boundaries $\varGamma_j^m$ for $j = 1, 2, \dots, n$,

● $b\in(0, 1)$ is specified uniformly, which denotes the ratio between the lengths of the short axes and the long axes of ellipses $\varGamma_j^m$ for $j = 1, 2, \dots, n$.

$Q_j^m$ and $a_j^m$ can be derived by the following iterative method.

Remark: The displacement of the ellipses centers and the incrementing of the half lengths of the long axes are controlled by giving different values of $\beta$ in iterations. When $\beta = 1$, the Algorithmic 1 is the iteration method proposed in ^[30]. Max and $\varepsilon_1$ are the specified maximum number of iterations allowed and the tolerance, respectively.

3.2.   The preimage domain of the unit disk with logarithmic spiral slit region

In this Subsection, we propose an iterative method based on the boundary integral equation, which calculates a preimage domain of a specified unit disk with logarithmic spiral slit region.

Theorem 3.2. (^[32]) For the specified oblique angles $\theta(t) = (\frac{\pi}{2}, \theta_2, \dots, \theta_n)$ of spiral slits, suppose $D$ is the bounded multiply connected region inside the unit disc. $\varGamma_1 = \left\lbrace |z| = 1\right\rbrace$ is a boundary component of $D$ and $0\in D$. Then there exists a unique analytical function $\Upsilon(z)$ in $D\cup\varGamma$, which satisfies the normalization conditions $\Upsilon(0) = 0, \; \Upsilon^{'}(0) > 0$ and conformally maps $D$ onto the spiral slit region $W$ bounded by $L_1 = \left\lbrace |w| = 1\right\rbrace$. The spiral slits $L_2, \; L_3\; , \dots, \; L_n$ are inside of $L_1$.

Conformal mapping conforming to Theorem 3.2 can be defined by(see ^[17,32])

where the analytic function $g(z)$ defines on $D\cup\varGamma$, and $c$ is an undetermined real constant. $w = \Upsilon(z)$ maps $\varGamma_j$ onto the spiral slit $L_j$ for $j = 1, 2, \dots, n$ (see Figure 2). According to the boundary conditions (3.1), the boundary values of $\Upsilon(z)$ should satisfy

Plugging (3.9) into (3.10) leads to

where $c$ and $h_j = R_j-{\rm Re}\left[{\rm e}^{{\rm i}(\frac{\pi}{2}-\theta_j)}\ln(c)\right]$ for $j = 2, 3, \dots, n$ are undetermined real constants obtained by solving the R-H problem with

According to Theorem 2.1, by solving the integral equation (2.7) and calculating (2.8), we can obtain $\mu(t)$ and $h(t) = (-\ln(c), h_2, \dots, h_n)$ such that

which satisfie the boundary conditions (3.11). For $z\in D$, we also use the Cauchy integral (3.6) to calculate the values of $g(z)$, and conformal mapping $\Upsilon(z)$ is obtained according to (3.9).

The images of the spiral slits $L_j$ under the mapping

are still the line segments $S_j$ for $j = 2, 3, \dots, n$ in the $\zeta$ plane (see Figure 2). Thus we uniformly use $l_j^*$ and $c_j$ for $j = 2, 3, \dots, n$ to denote the lengths and upper endpoints of $S_j$. Any point on a line segment can be represented as

Therefore, any point on a spiral slit still can be written as

As mentioned in Sect.3.1, for $D$ is unknown, we will present an iterative method for computing a preimage domain $D$ of a specified unit disk with spiral slit region $W$ and conformal mapping $\omega = \Upsilon(z)$ that maps $D$ onto $W$. Each step iteration of our iterative method will also generate a preimage domain $D^m$, and the sequence $\left\lbrace D^m\right\rbrace$ of the preimage domains will converge to our needed preimage domain $D$. For $m = 0, 1, 2, \dots,$ we suppose the region $D^m$ is inside the unit circle $\varGamma_1^m$ that is parametrized by

and is outside the ellipses $\varGamma_j^m$ parametrized by

where

● $Q_j^m$ denote the coordinates of the centers of ellipses $\varGamma_j^m$ for $j = 2, 3, \dots, n$,

● $a_j^m$ represent half the lengths of the long axes of elliptical boundaries $\varGamma_j^m$ for $j = 2, 3, \dots, n$,

● $b\in(0, 1)$ is specified uniformly, which denotes the ratio between the lengths of the short axes and the long axes of ellipses $\varGamma_j^m$ for $j = 2, 3, \dots, n$.

$Q_j^m$ and $a_j^m$ can be derived by the following iterative method.

Remark: The displacement of the ellipses centers and the incrementing of the half lengths of the long axes are controlled by giving different values of $\beta$ in iterations. Max and $\varepsilon_2$ are the specified maximum number of iterations allowed and the tolerance, respectively.

4.   Numerical examples

In this Section, we propose three numerical examples to test the availability of our iterative methods.

It is of great significance to reach the flow around bodies in energy and power engineering, aerospace engineering, and ocean engineering. Conformal mappings have unique superiorities in simulating the flow around bodies. In three numerical examples, we will simulate the irrotational plane flow around bodies in region $D$. Let $c(z) = \phi(x, y)+{\rm i}\psi(x, y)\; \text{and }\; C(w)$ are the complex potentials of the flows in the fluid region $D$ and $W$, respectively, where $\phi$ and $\psi$ denote the velocity potential and the stream function. $w = \Upsilon(z)$ conformally maps $D$ onto the slit region $W$. At the same time, the family of equipotential lines and streamlines in $D$ will also be mapped by $w = \Upsilon(z)$ onto the equipotential lines and streamlines in $W$. Thus, when $C(w)$ is known, the complex potential in $D$ is given by

The boundary integral equation (2.7) and the function (2.8) can be calculated with precision by means of the MATLAB function $\texttt{fbie}$. We set the parameters $\texttt{restart} = \left[\; \right]$, $\texttt{iprec} = 5$, $\texttt{maxit} = 100$, and $\texttt{gmrestol} = 0.5\times 10^{-14}$, in $\texttt{fbie}$. $N = 2048$ is the number of equidistant nodes in the discretization of each boundary component. Using the MATLAB function $\texttt{fcau}$, we can rapidly derive the values of $g(z)$ for $z\in D$. For more details, see ^[19].

In numerical examples, the figures illustrate the availability of our iterative methods and discuss the influences of $\beta$ on the number of iterations needed for the iterative methods to converge:

● Figure 3, Figure 6 and Figure 9 show that both the specified slit region and the calculated preimage domain are composed of boundaries and orthogonal grids. Furthermore, using the mapping $\zeta(w)$, the boundaries of the given spiral slit domains are mapped onto the line segments in Figure 3(c), Figure 6(c) and Figure 9(c).

● Figure 4(b) and Figure 8 show the influences generated by the different values of $b$ and $\beta$ on the number of iterations needed for the iterative methods to converge for Example 4.1 and Example 4.2.

● Figure 5 and Figure 10(b) show the influences of several values of the ratio $b$ on the number of iterations needed for the iterative methods to converge and display the error variation in the iterative process, in which the value of $\beta$ is optimal for Example 4.1, Example 4.2, and Example 4.3.

● Figure 4(a), Figure 7 and Figure 10(a) show the results of simulating the irrotational plane flow around bodies.

Example 4.1 Regard the unbounded domain that is the entire $w$-plane with three logarithmic spiral slits for $\theta(t) = (\frac{\pi}{4}, \frac{\pi}{4}, \frac{\pi}{4})$ (see Figure 3(b)). When using the tolerance $\varepsilon_1 = 1\times 10^{-10}$, the iterative method converges for $0.05\leqslant b\leqslant 0.8$ (see Figure 5(a)). Figure 4(a) shows streamlines of the flow around bodies generated by a spiral point vortex $C(w) = \frac{5-5{\rm i}}{2\pi {\rm i}}\ln(w)$.

Example 4.2 Regard the unbounded domain that is the entire $w$-plane with five logarithmic spiral slits for $\theta(t) = (\frac{\pi}{4}, \frac{\pi}{4}, \frac{\pi}{4}, \frac{\pi}{4}, \frac{\pi}{4})$ (see Figure 6). When using tolerance $\varepsilon_1 = 1\times 10^{-10}$ and $0.05\leqslant b\leqslant 0.45$, the iterative method should be convergent (see Figure 8). Figure 7 shows streamlines and equipotential lines of the flows around bodies generated by a uniform stream $C(w) = w$ and a source point $C(w) = \ln(w)$, respectively.

Example 4.3 Regard the unit disk with logarithmic spiral slit region that oblique angles are $\frac{\pi}{4}$ (see Figure 9). When using tolerance $\varepsilon_2 = 1\times 10^{-4}$ and $0.05\leqslant b\leqslant 0.25$, the iterative method should be convergent. At the same time, equipotential lines of the flow around bodies generated by ten source points $C(w) = \mathop{\sum}\limits_{j = 2}^{n} \ln(w-Q_j^*)$ are shown in Figure 10.

Figure 4(b), Figure 5, Figure 8 and Figure 10(b) illustrate the availability and convergence of the iterative methods proposed in Sect.3. Figure 4(b) and Figure 8 show that the values of $\beta$ lead to a significant impact on the number of iterations needed for the iterative methods to converge. In Example 4.1 and Example 4.2, when $\beta = 1$, the iterative methods are non-convergence for some ratio $b$. Therefore, it is necessary to choose an appropriate value of $\beta$, which can improve the convergence of our iterative methods and reduce the number of iterations needed for the iterative methods to converge.

5.   Concluding

In this paper, for a specified multiply connected logarithmic spiral slit region, we propose an iterative method that successfully calculates a conformally equivalent preimage domain and conformal mapping that maps the preimage domain onto the specified spiral slit region. In the proposed iterative method, we introduce parameter $\beta$ to control the changes to the boundaries in the preimage domain, ensuring the convergence of that approach. Through several numerical examples, we also explore how parameters $\beta$ affects the proposed iterative method's convergence and show how many iterations are necessary when parameter $\beta = 2$ and $\beta = 20$ for unbounded and bounded logarithmic spiral slit regions, respectively. Moreover, in physical planes, using the above conformal mappings, we simulate the several flows around bodies, including the spiral point vortex, the uniform stream, and the source point, which develop the applications of conformal mappings in the simulations of the flow around bodies. Since our method can calculate the complex potentials of spiral point vortices and other potential flows at any point in the logarithmic spiral slit regions, it provides a possibility to explore the dynamic behavior of potential flows in that domain in the future.

Acknowledgment

KW is grateful to Prof. Mohamed M. S. Nasser for his valuable suggestions. This work was supported by the National Natural Science Foundation of China (NSFC) (Grant No. 11461037) and Yunnan Fundamental Research Projects (Grant NO. 202101BE070001-050).

Conflict of interest

The authors declare that they have no conflict of interest.