Circular Slit Maps of Multiply Connected Regions with Application to Brain Image Processing

In this paper, we present a fast boundary integral equation method for the numerical conformal mapping and its inverse of bounded multiply connected regions onto a disk and annulus with circular slits regions. The method is based on two uniquely solvable boundary integral equations with Neumann-type and generalized Neumann kernels. The integral equations related to the mappings are solved numerically using combination of Nyström method, GMRES method, and fast multipole method. The complexity of this new algorithm is O((M+1)n)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O((M + 1)n)$$\end{document}, where M+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M+1$$\end{document} stands for the multiplicity of the multiply connected region and n refers to the number of nodes on each boundary component. Previous algorithms require O((M+1)3n3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O((M+1)^3 n^3)$$\end{document} operations. The numerical results of some test calculations demonstrate that our method is capable of handling regions with complex geometry and very high connectivity. An application of the method on medical human brain image processing is also presented.


Introduction
Conformal mapping has been applied to the study of the visual cortex.Frederick and Schwartz [1] presented the conformal mapping of retina as a case of identification of a single point (the representation of the blind spot, or optic disk) in the eye.They mapped one half of the retina conformally to the half unit disk, which conformally is mapped to a unit disk preserving the point in the half disk that mapped to the origin of the target disk as a parameter of this mapping.Finally, the unit disk is mapped conformally to an arbitrary 2D region, the flattened data.More applications can be found in [2,3].
Explicit formulae for conformal mappings of multiply circular regions onto canonical slit regions in terms of Schottky-Klein prime function and their applications are described in [4][5][6].Numerical conformal mappings of general multiply connected regions onto canonical slit regions have been described in [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].The obtained boundary integral equation methods in [16][17][18][19][20][23][24][25] were discretized by the Nyström method with trapezoidal rule to obtain a dense and nonsymmetric linear system.Gauss elimination method of order O((M + 1) 3 n 3 ) operations was used for solving the obtained linear system, where M + 1 refers to the multiplicity of the multiply connected region and n stands for the number of nodes in the discretization of each boundary component.Hence, it is very expensive to solve the linear system for very large values of M and n.
This paper illustrates a new integral equation method using the adjoint generalized Neumann and Neumann-type kernels for conformal mapping and its inverse of a bounded multiply connected region onto a disk and annulus with circular slits, extending the work presented in Sangawi et al. [16,17].Unlike the methods presented in [16,17] which require solving three integral equations, the method shown in this paper requires only two integral equations.Discretization of the integral equation yields a system of linear equations which are solved by an iterative method based on the restarted version of the generalized minimal residual (GMRES) method of Saad and Schultz [26].
The GMRES method will converge significantly faster without the need to use preconditioning for solving our linear system.See [27][28][29][30][31] for details on FMM operation and memory requirement.This paper is an improvement of the methods presented in [16,17].The matrix-vector product function can be computed using the function gmres.The function zfmm2dpart in the MATLAB toolbox FMMLIB2D developed by Greengard and Gimbutas [32] is then used for the matrix-vector product function for the coefficient matrix of our linear system.Thus, the obtained linear system can be solved in O((M + 1)n) operations.We then apply the proposed method to brain medical image processing which gives high-quality information that helps in diagnosing diseases and making the treatment easier.
The disposition of this paper is as follows: Sect. 2 presents some auxiliary materials.Section 3 discusses an integral equation method to calculate the modulus of conformal mapping function f .Derivations of integral equations related to f for canonical regions of type disk with circular slits Ω 1 and annulus with circular slits Ω 2 are given in Sects.4 and 5 .Section 6 presents a method to compute the inverse mapping functions from the canonical slits region onto the original region.Section 7 describes the fast

Notations and Auxiliary Materials
Let Ω be a bounded multiply connected region of connectivity M + 1.The boundary Γ consists of M + 1 smooth Jordan curves Γ j , j = 0, 1, . . ., M as shown in Fig. 1.
The curves Γ j are parametrized by 2π -periodic twice continuously differentiable complex functions z j (t) with non-vanishing first derivatives The total parameter domain J is the disjoint union of M + 1 intervals J 0 , . . ., J M .We define a parametrization z(t) of the whole boundary Γ on J by z(t) = z j (t), t ∈ J j . ( Let H * be the space of all real Hölder continuous 2π -periodic functions ω(t) such that Suppose that b(z) and H (z) are complex-valued functions defined on Γ such that b(z) = 0, H (z) = 0 and H (z)/T (z) 2 satisfies the Hölder condition on Γ .Then, the interior relationship is defined as follows: A complex-valued function P(z) is said to satisfy the interior relationship if P(z) is analytic in Ω and satisfies the non-homogeneous boundary relationship where G(z) is analytic in Ω, Hölder continuous on Γ , and G(z) = 0 on Γ .The boundary relationship (2) also has the following equivalent form: The following theorem from [18] gives an integral equation for an analytic function satisfying the interior non-homogeneous boundary relationship (2) or (3). Res where The symbol "conj" in the superscript denotes complex conjugate, while the minus sign in the superscript denotes limit from the exterior.The sum in (4) is over all those zeros a 1 , a 2 , . . ., a M of G that lie inside Ω .If G has no zeros in Ω , then the term containing the residue in (4) will not appear.

Compute the Piecewise Real Function h j
Let A(t) be a complex continuously differentiable 2π -periodic function for all t ∈ J .The generalized Neumann kernel formed with A is defined by [13,18] The adjoint function to the function Ã is given by The generalized Neumann kernel Ñ (s, t) formed with Ã is given by where N * (s, t) = N (t, s) is the adjoint kernel of the generalized Neumann kernel N (s, t) (see [13] for more details).Define the Fredholm integral operator N * by It is known that λ = 1 is an eigenvalue of the kernel N with multiplicity 1 and λ = −1 is an eigenvalue of the kernel N with multiplicity M [13].The eigenfunctions of N corresponding to the eigenvalue λ = −1 are χ [0] , χ [1] , . . ., χ [M] , where χ [ j] (ξ ) = 1, ξ ∈ Γ j , 0, otherwise, j = 0, 1, . . ., M.
Theorem 2 [16] Suppose that the function γ ∈ H * and functions h, μ ∈ S such that are boundary values of an analytic function g(z) in Ω.Then, the function h = (h 0 , h 1 , . . ., h M ) is given by where φ [ j] are solutions of the following integral equations

Disk with Circular Slit Region
Assume that an analytic function w = f (z) maps Ω conformally onto a disk with circular slits Ω 1 (see Fig. 1b).This canonical region Ω 1 consists of the disk |w| ≤ μ 0 slit along M arcs of circles |w| = μ j , j = 1, . . ., M. We assume that w = f (z) maps Γ 0 onto the circle |w| = μ 0 , and the curves Γ 1 , . . ., Γ M onto the arcs of circles |w| = μ j , j = 1, . . ., M, where μ 0 is known, and μ j , j = 1, . . ., M, are unknowns and have to be determined.The function is uniquely determined by prescribing that f (a) = 0 and f (a) > 0. The boundary values of f can be represented in the form where θ j are the boundary correspondence functions of Γ j and μ j are the radii of the circular slits.The unit tangent to Γ at z(t) is denoted by . Thus, it can be shown from (12) that Note that for j = 0, θ j > 0 and for j = 1, 2, . . ., M, θ j may be positive or negative since each circular slit f (Γ j ) is traversed twice (see Fig. 1c).Thus θ j (t) θ j (t) = ±1.
Squaring both sides of (13), we get The mapping function can also be expressed as [11,16] where ĥ(z) is an analytic function and c = f (a) is an undetermined real constant.From ( 14) and (15), we obtain the boundary relationship By taking logarithm on both sides of (15), we obtain which implies where By obtaining h j , j = 0, 1, . . ., m, from (10), we obtain and Comparison of ( 16) and ( 3) leads to a choice of where Note that a is a pole of order one for f (w) . Hence, Combining the integral equation ( 21) and ( 23), we get 123 Write the integral equation (24) as By using single-valuedness of the mapping function f leads to the following conditions Thus, the integral equation (25) with condition ( 26) should give a unique solution provided the parameters c and | f (z q )| = μ q , q = 1, 2, . . ., M that appear in (25) are known.By solving the integral equation ( 11), we get φ [ j] , j = 0, 1, . . ., M, which gives h j through (10) which in turn gives c and μ j through (19) and (20).By solving (25) with (26) and the known values of c and μ j , we get the values of f .Then, the boundary values of ĥ(z) are given by ĥ The logarithm in (27) indicates complex logarithm which define as log(z) = ln |z| + iArg(z).
In This paper, we define the range of Arg(z) as [0, 2π) to avoid the multiple value of Arg(z).Finally, the approximate boundary values of f (z) are obtained from (15), i.e., The approach presented here is an improvement over [16].Unlike [16], no integral equation for θ (t) is required here for the computation of f (z).
Note that if Ω is a simply connected region, then (25) reduced to where For simply connected region, we can just solve integral equation (22) for f (z(t)) and compute f (z(t)) using ( 14).

Annulus with Circular Slit Region
Let w = f (z) be the analytic function which maps Ω conformally onto an annulus with circular slits Ω 2 (see Fig. 1c).This canonical region Ω 2 consists of a circular annulus where μ 0 is known, and μ j , j = 1, . . ., M, are unknowns and have to be determined.
The boundary values of f can be represented in the form where θ j are the boundary correspondence functions of Γ j and μ j are the radii of the circular slits.The mapping function can be expressed as [11,17] where ĥ(z) is an analytic function and let 0 ∈ Ω and b be an interior point in Γ M .The mapping is made uniquely determined by prescribing that f (b) = 0 and c = f (0) > 0 where c is an undetermined real constant.From (32), we get which implies The above equation means either both f (0) and ĥ(0)− 1 b are reals or both pure imaginary complex numbers.The representation (31) also satisfy the boundary relationship (14).Multiplying both sides of ( 14) by z gives the following boundary relationship Comparison of ( 35) and (3) leads to a choice of where Note that zero is a pole of order one for Thus, integral equation, (36) after dividing both sides by f (0), becomes By using single-valuedness of the mapping function, f leads to the following conditions and Thus, the integral equation (39) with conditions ( 40) and ( 41) should give a unique solution provided the parameters c and | f (z q )| = μ q , q = 1, 2, . . ., M that appear in (39) are known.The values of μ can be computed by using the integral equation ( 11), ( 10), ( 19) and (20) with γ (t) = ln |b| − ln |b − z(t)| and A(t) = z(t).By solving (39) with the known values of c and μ, we get the boundary values of f .Then, the following equation gives the values of ĥ(z The logarithm in (42) indicates complex logarithm.Section 4 (after 27) explained the definition of complex algorithm and range of argument to avoid multiple computation.
If f (0) and ( ĥ(0) − 1 b ) are both pure imaginary, then the boundary values of ĥ(z) are given by ĥ Finally, the approximate boundary values of f (z) are given by

Computing Values of Mapping Functions Interior and Inverse Mapping Functions
The approximate interior values of the function f (z) for both canonical slits regions are calculated by the Cauchy integral formula in the form of numerically where Γ 1 w − z dw = 1, the integral in the numerator has the advantage that the denominator in this formula compensates for the error in the numerator (see [33]).The integrals in (45) are approximated by the trapezoidal rule.
For computing the inverse maps, note that the mapping function f −1 (w) = z is analytic in the respective canonical slits region.Then by means of Cauchy's integral formula, we obtain z ∈ Ω by

Numerical Implementation
Since the functions z j (t) are 2π -periodic, a reliable procedure for solving the integral equations ( 11), ( 25) and ( 39) numerically is by using the Nyström's method with the trapezoidal rule [34].The trapezoidal rule is the most accurate method for integrating periodic functions numerically [35, pp. 134-142].The details for solving integral equation ( 11) is given in [36,37].

Solving Integral Equation with Neumann-Type Kernel for Conformal Mapping onto Unit Disk with Circular Slits
For solving integral equation ( 25), multiplying both sides by |z (t)| gives We choose n equidistant collocation points for Γ j t j = (( j − 1)modn)2π Using the parametric representation z(t) on Γ and discretize integral equation ( 47) and (26) gives We define vector hpx as Then, (49) can be written as Adding ( 51) and (48) gives Rearrange (52) yields Let E and G be n × n matrices with the elements and Take t to be the vector with element t k , where k = 1, 2, . . ., (M +1)n, a be (M +1)n×1 vector with element a and Then, (53) can be written in matrix form as From (56), multiplication of matrix G by a vector can be done by FMM of order O((M + 1)n) operations.Let the (M + 1)n × 1 vector p be the vector f (z(t))z (t) and the (M + 1)n × 1 vector q be the vector Let also b and Hence the products Ep and Gq can be computed by using the MATLAB function zfmm2dpart as follows [32]: and In this paper, we take i prec = 4, which implies the tolerance of the FMM is 0.5 × 10 −12 .The solution of the system (56) can be solve by using GMRES method.We use the MATLAB function gmres to solve the system (56) with relative residual < 10 −12 as the stopping criteria.Once the solution f (z(t))z (t) f (0) has been computed, the mapping function f (z(t)) is computed, by the formulas ( 27) and (28).The computation of interior values and inverse value via Cauchy integral formula (45) and ( 46) is referred to Nasser [37].For computing Cauchy integral formula, we also discretized integrals using Nyström method with trapezoidal rule.The summation can be done by using FMM.

Solving Integral Equation with Neumann-Type Kernel for Conformal Mapping onto the Annulus with Circular Slits
For solving integral equation (39), multiplying both sides of equation ( 39) by |z (t)| gives We choose n equidistant collocation points by choosing the same equidistant collocation points as in previous subsection.Using the parametric representation z(t) on Γ and discretizing integral equations ( 61), ( 40) and (41) gives and We define vector hp as Then, (49) can be written as Adding (66), ( 64) and (62) gives Rearrange (67) yields Let E and G be (M + 1)n × (M + 1)n matrices as in ( 54) and ( 55), take t be vector with element t k , where k = 1, 2, . . ., (M + 1)n and Then, (68) can be written in matrix form as Let the (M + 1)n × 1 vector p be the vector | f (z(t))z (t) f (0) and the (M + 1)n × 1 vector q be the vector Let also b and c be 2 × (M + 1)n real vectors as given in ( 57) and ( 58).The products Ep and Gq can be computed by using the MATLAB function zfmm2dpart according to ( 59) and (60), respectively.We take i prec = 4, which implies the tolerance of the FMM is 0.5×10 −12 .The solution of the system (69) can be solved by using GMRES method.We use the MATLAB function gmres to solve the system (69) with relative residual < 10 −12 as the stopping criteria.
For computing Cauchy integral formula, we also discretized integrals using Nyström method with trapezoidal rule.The summation arising from discretization can be done by using FMM.

Numerical Examples
For numerical experiments, we have used some test regions of connectivity one, two and fifteen.All the computations were done using MATLAB 7.12.0.Then, the exact mapping function is given by [8, p. 340] See Figs. 2, 3 and Table 1 for numerical results.

Example 3 Ellipse with Two Circles:
Let Ω be the region bounded by [41,42] Γ 0 : {z(t) = 2 cos t + i sin t}, This example is also discussed in Reichel [41] and Kokkinos et al. [42] which allows for comparison of the radii μ 1 , μ 2 .Since our condition problem is different from them, we should multiply our radii by 2.5 to compare with the result of Reichels and Kokkinos et al.We choose the symbols μ 1R and μ 2R for the radii in Reichel [41] and μ 1K and μ 2K for the radii in Kokkinos et al. [42].Figures 6 and 7 show the region and its images based on our method.See Tables 4, 5 and 6 for comparisons between   our computed values of μ 1 and μ 2 with those computed values of Reichel [41] and Kokkinos et al. [42] and computed values of c D , c A , μ A1 and μ A2 .
Example 4 Consider the region Ω of connectivity fifteen [43], z j (t) = ξ j + e iσ j (a j cos t + ib j sin t), j = 0, 1, . . ., 14, where a = −0.5 − i0.5, b = −0.833− i2.165 and the values of the complex constants ξ j and the real constants a j , b j and σ j are as in Table 7.We chose the symbols μ Dj and μ A j for radii in previous two canonical slits regions, respectively.The numerical results are presented in Figs. 8 and 9. See  where a = 3.9 + i0.2 and b = −4.6871− i1.3355 and the values of the complex constants ξ j are as in Table 9. Mapping function from the original region onto the disk and annulus with slits region and the inverse mapping functions from the disk and

Medical Image Processing Applications
Medical image processing has experienced a dramatic expansion.Experts from applied mathematics, engineering, computer sciences, biology and medicine have been attracted by biomedical research.One of the important parts of clinical routine is computer analytic processing.Image processing is a very important process for analyzing images and gives high-quality information that helps in diagnosing diseases and making treatment easier.A conformal mapping is a one-one and onto transformation w = f (z) that preserves both local angles and shape of infinitesimal small figures but not necessarily their size which is necessary factor in medical image processing.
Conformal mapping can be effectively applied in the field of brain surface mapping.Parameterization of anatomical surfaces in brain imaging is valuable to help analyze anatomical shape.Conformal mapping involving slits has been applied by Wang et al. [3] to study the difference in cortical surface of the brain.It helps in disease diagnosis related to cortical surface such as Alzheimer's disease.Conformal mapping can be used to map an irregular surface onto a disk and annulus while preserving the angle, which is useful for visualization of magnetic resonance imaging (MRI).MRI is a tomographic imaging procedure which uses strong magnetic fields and radio-waves to create images of internal physical and chemical characteristics of an object.An MRI scanner is shown in Fig. 12.
The conformal mapping method of a simply and multiply connected region is applied to a brain surface image onto a disk and annulus with circular slit regions.The results based on our method are presented in Figs. 13 and 14 .Different values of a has the same effect as "warping" and "zooming" to get better and clearer pictures of different parts of the brain.Since conformal mappings are one-to-one, no information is lost.One is able to see how the zoomed part of the brain relates with the other parts.Thus conformal maps help preserve some essential features of visual information of the brain.

Conclusions
In this paper, we have constructed new boundary integral equations for conformal mapping of bounded multiply connected regions onto the disk and annulus with circular slits regions.The advantage of the proposed method is that the boundary integral equations are all linear.Also in contrast to the boundary integral equation used in Nasser [10,11], the right-hand sides of newly constructed boundary integral equations are much simpler and do not contain any singular operator.Moreover, several mappings of the test regions of connectivity one, two, three, fifteen and ninety-two were computed numerically using the proposed method.After computing the boundary values of the mapping function, the interior values of the mapping function and its inverse were calculated by means of Cauchy integral formula.The numerical examples presented have illustrated that the innovative boundary integral equation method has high accuracy.Luckily, this rate of accuracy is highly vital for medical application of image analysis.This paper used MATLAB functions graythresh, imfill, bwareaopen and bwboundaries to recognize and extract object boundary from medical images reproducibly and accurately.Consequently, this led to the achievement of better and clearer images.In the end, one should bear this in mind that such results were achieved after applying the proposed method to the medical image objective that helps in diagnosing diseases and making treatment easier.

Fig. 1
Fig. 1 Mappings of the bounded multiply connected region onto a disk Ω 1 and annulus with circular slits Ω 2

Fig. 2
Fig. 2 Mapping a region bounded by unit circle onto a disk for Example 1

Fig. 13 Fig. 14
Fig. 13 Application of the proposed method on the human brain.Computation times with n = 2 12 are b 2.239811 s, c 2.260020 s, d 2.270281 s, e 2.258278 s, f 2.262168 s

Fig. 15 Fig. 16
Fig. 15 Original MR image and the transformed images based on conformal mapping by choosing different point in the original image.Computation times with n = 2 12 are b 2.470122 s, c 2.515067 s, d 2.313659 s, e 2.414645 s, f 2.563450 s

Table 1
Error norm (unit circle)

Table 5
Radii comparison forExample 3 for disk with slits Ω 1

Table 7
Values of constants a j , b j , ξ j and σ j for Example 4

Table 9
annulus with slits region onto the original region.The numerical results are presented in Figs.10 and 11.