# On the polar derivative of a polynomial

N. A. Rather, S. H. Ahangar and Suhail Gulzar  Department of Mathematics
University of Kashmir
Srinagar, Hazratbal 190006
India<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Let $P(z)$ be a polynomial of degree $n$ having no zero in $|z|<k$ where
$k\geq 1,$ then for every real or complex number $\alpha$ with $|\alpha|\geq
1$ it is known

$\underset{|z|=1}{\max}|D_{\alpha}P(z)|\leq
n\left(\dfrac{|\alpha|+k}{1+k}\right)\underset{|z|=1}{\max}|P(z)|,$

where $D_{\alpha}P(z)=nP(z)+(\alpha-z)P^{\prime}(z)$ denote the polar
derivative of the polynomial $P(z)$ of degree $n$ with respect to a point
$\alpha\in\mathbb{C}.$ In this paper, by a simple method, a refinement of
above inequality and other related results are obtained.

###### Key words and phrases:

polynomials; Inequalities in the complex domain; Polar derivative; Bernstein’s
inequality.

###### 2010 Mathematics Subject Classification:

30A10, 30C10, 30E10

##  1\.  Introduction and statement of results

If $P(z)$ is a polynomial of degree $n,$ then concerning the estimate of
$|P^{\prime}(z)|$ on the unit disk $|z|=1,$ we have

$\underset{\left|z\right|=1}{\max}\left|P^{\prime}(z)\right|\leq
n\,\underset{\left|z\right|=1}{\max}\left|P(z)\right|.$ (1.1)

Inequality (1.1) is an immediate consequence of Bernstein’s inequality on the
derivative of a trigonometric polynomial (for reference, see [8], [9] or
[10]). Equality in (1.1) holds for $P(z)=az^{n},$ $a\neq 0.$
If we restrict ourselves to the class of polynomials having no zero in
$|z|<1$, then inequality (1.1) can be replaced by

$\underset{\left|z\right|=1}{\max}\left|P^{\prime}(z)\right|\leq\frac{n}{2}\underset{\left|z\right|=1}{\max}\left|P(z)\right|.$
(1.2)

Inequality (1.2) was conjectured by Erdös and later verified by Lax [5].The
result is sharp and equality holds for $P(z)=\alpha z^{n}+\beta,$
$|\alpha|=|\beta|\neq 0.$
As an extension of (1.2), Malik [7] proved that if $P(z)$ is a polynomial of
degree $n$ which does not vanish in $|z|<k$ where $k\geq 1,$ then

$\underset{\left|z\right|=1}{\max}\left|P^{\prime}(z)\right|\leq\frac{n}{1+k}\underset{\left|z\right|=1}{\max}\left|P(z)\right|.$
(1.3)

The result is best possible and equality in (1.3) holds for $P(z)=(z+k)^{n}.$
Let $D_{\alpha}P(z)$ denote the polar derivative of the polynomial $P(z)$ of
degree $n$ with respect to the point $\alpha\in\mathbb{C},$ then

$D_{\alpha}P(z)=nP(z)+(\alpha-z)P^{\prime}(z).$

The polynomial $D_{\alpha}P(z)$ is of degree at most $n-1$ and it generalizes
the ordinary derivative in the sense that

$\underset{\alpha\rightarrow\infty}{\lim}\dfrac{D_{\alpha}P(z)}{\alpha}=P^{\prime}(z).$

A. Aziz [1] extended inequality (1.3) to the polar derivative and proved that
if $P(z)$ is a polynomial of degree $n$ having no zero in $|z|<k$ where $k\geq
1,$ then for $\alpha\in\mathbb{C}$ with $|\alpha|\geq 1$

$\underset{|z|=1}{\max}|D_{\alpha}P(z)|\leq
n\left(\dfrac{|\alpha|+k}{1+k}\right)\underset{|z|=1}{\max}|P(z)|.$ (1.4)

The result is best possible and equality in (1.4) holds for the polynomial
$P(z)=(z+1)^{n}.$
The bound in (1.4) depends only upon the modulus of the zero of smallest
modulus and not on the moduli of other zeros. It is interest to obtain a bound
which depends upon the location of all the zeros rather than just on the
location of the zero of smallest modulus.

In this paper, by a simple method, we first present the following result which
is a refinement of inequality (1.4).

######  Theorem 1.1.

Let $P(z)=a_{n}\prod_{\nu=1}^{n}(z-z_{\nu})$ be a polynomial of degree $n.$ If
$|z_{\nu}|\geq k_{\nu}\geq 1$ where $1\leq\nu\leq n,$ then for
$\alpha\in\mathbb{C}$ with $|\alpha|\geq 1,$

$\underset{|z|=1}{\max}|D_{\alpha}P(z)|\leq
n\left(\dfrac{|\alpha|+t_{0}}{1+t_{0}}\right)\underset{|z|=1}{\max}|P(z)|,$
(1.5)

where

$\displaystyle
t_{0}=\begin{cases}1+\dfrac{n}{\sum_{\nu=1}^{n}\frac{1}{k_{\nu}-1}}&\text{if
}\quad k_{\nu}>1\quad\text{for all}\quad\nu,\,1\leq\nu\leq n\\\
1&\text{if}\quad k_{\nu}=1\quad\text{for some}\quad\nu,\,1\leq\nu\leq
n.\end{cases}$ (1.6)

######  Remark 1.2.

If $k_{\nu}\geq k,$ $k\geq 1$ for $1\leq\nu\leq n,$ then $t_{0}\geq k$ which
implies

$\dfrac{|\alpha|+t_{0}}{1+t_{0}}\leq\dfrac{|\alpha|+k}{1+k}\quad\textnormal{for}\quad|\alpha|\geq
1.$

This shows (1.5) is refinement of inequality (1.4).

######  Remark 1.3.

If we divide the two sides of inequality (1.5) by $|\alpha|$ and letting
$|\alpha|\rightarrow\infty,$ we get a result due to Govil et. al. [4].

Next, as an application of Theorem 1.1, we present the following result.

######  Theorem 1.4.

Let $P(z)=a_{n}\prod_{\nu=1}^{n}(z-z_{\nu})$ be a polynomial of degree $n$
with $P(0)\neq 0.$ If $|z_{\nu}|\leq k_{\nu}\leq 1,$ $1\leq\nu\leq n,$ then
for $\delta\in\mathbb{C}$ with $|\delta|\leq 1,$

$\underset{|z|=1}{\max}|D_{\delta}P(z)|\leq
n\left(\dfrac{1+|\delta|\,s_{0}}{1+s_{0}}\right)\underset{|z|=1}{\max}|P(z)|,$
(1.7)

where

$\displaystyle
s_{0}=\begin{cases}1+\dfrac{n}{\sum_{\nu=1}^{n}\frac{k_{\nu}}{1-k_{\nu}}}&\text{if
}\quad k_{\nu}<1\quad\text{for all}\quad\nu,\,1\leq\nu\leq n\\\
1&\text{if}\quad k_{\nu}=1\quad\text{for some}\quad\nu,\,1\leq\nu\leq
n.\end{cases}$ (1.8)

######  Remark 1.5.

If $k_{\nu}\leq k\leq 1$ for $1\leq\nu\leq n,$ then $1/k\leq s_{0}$ which
implies

$\dfrac{1+|\delta|\,s_{0}}{1+s_{0}}\leq\dfrac{|\delta|+k}{1+k}\quad\textnormal{for}\quad|\delta|\leq
1.$

Therefore, it follows that if $P(z)$ is a polynomial of degree $n$ having all
its zeros in $|z|<k$ where $k\leq 1,$ then for $\delta\in\mathbb{C}$ with
$|\delta|\leq 1,$

$\underset{|z|=1}{\max}|D_{\delta}P(z)|\leq
n\left(\dfrac{|\delta|+k}{1+k}\right)\underset{|z|=1}{\max}|P(z)|.$ (1.9)

The result is sharp.

##  2\. Lemmas

For the proof of these theorems, we need the following lemmas. The first lemma
is due to Gardner and Govil [2].

######  Lemma 2.1.

Let $P(z)=a_{n}\prod_{\nu=1}^{n}(z-z_{\nu})$ be a polynomial of degree $n.$ If
$|z_{\nu}|\geq k_{\nu}\geq 1,$ $1\leq\nu\leq n,$ then for $|z|=1,$

$|Q^{\prime}(z)|\geq t_{0}|P^{\prime}(z)|,$ (2.1)

where $Q(z)=z^{n}\overline{P(1/\overline{z})}$ and $t_{0}$ is given by (1.6).

######  Lemma 2.2.

Let $P(z)$ be the polynomial of degree $n$ and
$Q(z)=z^{n}\overline{P(1/\overline{z})}$ then for $|z|=1,$

$|P^{\prime}(z)|+|Q^{\prime}(z)|\leq n\,\underset{|z|=1}{\max}|P(z)|.$ (2.2)

This is a special case of a result due to Govil and Rahman [6].

##  3\. Proof of Theorems

######  Proof of Theorem 1.1.

Let $Q(z)=z^{n}\overline{P(1/\overline{z})},$ then for $|z|=1$, it can be
easily verified that

$\displaystyle|P^{\prime}(z)|=|nQ(z)-zQ^{\prime}(z)|\quad\textnormal{and}\quad|Q^{\prime}(z)|=|nP(z)-zP^{\prime}(z)|.$
(3.1)

Now, for every real or complex number $\alpha$ and $|z|=1,$ we have by using
(3.1),

$\displaystyle|D_{\alpha}P(z)|$
$\displaystyle=|nP(z)+(\alpha-z)P^{\prime}(z)|$
$\displaystyle\leq|\alpha||P^{\prime}(z)|+|nP(z)-zP^{\prime}(z)|$
$\displaystyle=(|\alpha|-1)|P^{\prime}(z)|+|P^{\prime}(z)|+|Q^{\prime}(z)|.$
(3.2)

Multiplying the two sides of inequality (Proof of Theorem 1.1.) by $t_{0}$ and
using Lemma 2.1, we obtain for $|\alpha|\geq 1,$

$\displaystyle t_{0}|D_{\alpha}P(z)|$
$\displaystyle\leq(|\alpha|-1)t_{0}|P^{\prime}(z)|+t_{0}(|P^{\prime}(z)|+|Q^{\prime}(z)|)$
$\displaystyle\leq(|\alpha|-1)|Q^{\prime}(z)|+t_{0}(|P^{\prime}(z)|+|Q^{\prime}(z)|)\quad\textnormal{for}\quad|z|=1.$
(3.3)

Adding (Proof of Theorem 1.1.), (Proof of Theorem 1.1.) and using Lemma 2.2,
we get for $|\alpha|\geq 1$ and $|z|=1,$

$\displaystyle(1+t_{0})|D_{\alpha}P(z)|$
$\displaystyle\leq(|\alpha|+t_{0})(|P^{\prime}(z)|+|Q^{\prime}(z)|)$ (3.4)
$\displaystyle\leq n(|\alpha|+t_{0})\underset{|z|=1}{\max}|P(z)|,$

which gives

$\displaystyle|D_{\alpha}P(z)|\leq
n\left(\dfrac{|\alpha|+t_{0}}{1+t_{0}}\right)\underset{|z|=1}{\max}|P(z)|\quad\textnormal{for}\quad|z|=1.$
(3.5)

This completes the proof of theorem 1.1. ∎

######  Proof of Theorem 1.4.

Since $P(z)=a_{n}\prod_{\nu=1}^{n}(z-z_{\nu})$ where $|z_{\nu}|\leq
k_{\nu}\leq 1,$ $\nu=1,2,\cdots,n$ with $P(0)\neq 0,$
$Q(z)=z^{n}\overline{P(1/\overline{z})}=\overline{a_{n}}\prod_{\nu=1}^{n}(1-\bar{z}_{\nu}z)$
is a polynomial of degree $n$ with
$\left|\frac{1}{z_{\nu}}\right|\geq\frac{1}{k_{\nu}}\geq 1.$ Applying Theorem
1.1 to the polynomial $Q(z)$ and noting that $|P(z)|=|Q(z)|$ for $|z|=1,$ we
obtain for $\alpha\in\mathbb{C}$ with $|\alpha|\geq 1$ and $|z|=1,$

$\displaystyle\underset{|z|=1}{\max}|D_{\alpha}Q(z)|\leq
n\left(\dfrac{|\alpha|+s_{0}}{1+s_{0}}\right)\underset{|z|=1}{\max}|P(z)|,$
(3.6)

where

$\displaystyle
s_{0}=\begin{cases}1+\dfrac{n}{\sum_{\nu=1}^{n}\frac{k_{\nu}}{1-k_{\nu}}}&\text{if
}\quad k_{\nu}<1\quad\text{for all}\quad\nu,\,1\leq\nu\leq n\\\
1&\text{if}\quad k_{\nu}=1\quad\text{for some}\quad\nu,\,1\leq\nu\leq
n.\end{cases}$

Again, for $|z|=1$ so that $z\overline{z}=1,$ we have

$\displaystyle|D_{\alpha}Q(z)|$
$\displaystyle=|nQ(z)+(\alpha-z)Q^{\prime}(z)|$
$\displaystyle=\left|nz^{n}\overline{P(1/\overline{z})}+(\alpha-z)\left\\{nz^{n-1}\overline{P(1/\overline{z})}-z^{n-2}\overline{P^{\prime}(1/\overline{z})}\right\\}\right|$
$\displaystyle=\left|\alpha\left\\{nz^{n-1}\overline{P(1/\overline{z})}-z^{n-2}\overline{P^{\prime}(1/\overline{z})}\right\\}+z^{n-1}\overline{P^{\prime}(1/\overline{z})}\right|$
$\displaystyle=\left|\alpha\left(n\overline{P(z)}-\overline{z}\overline{P^{\prime}(z)}\right)+\overline{P^{\prime}(z)}\right|$
$\displaystyle=|\overline{\alpha}nP(z)+(1-\overline{\alpha}z)P^{\prime}(z)|=|\overline{\alpha}||D_{1/\overline{\alpha}}P(z)|.$

Hence, from (3.6), we get

$\displaystyle|\alpha|\underset{|z|=1}{\max}|D_{1/\bar{\alpha}}P(z)|\leq
n\left(\dfrac{|\alpha|+s_{0}}{1+s_{0}}\right)\underset{|z|=1}{\max}|P(z)|.$
(3.7)

Replacing $1/\overline{\alpha}$ by $\delta,$ we obtain for every real or
complex number $\delta$ with $|\delta|\leq 1,$

$\displaystyle\underset{|z|=1}{\max}|D_{\delta}P(z)|\leq
n\left(\dfrac{1+|\delta|\,s_{0}}{1+s_{0}}\right)\underset{|z|=1}{\max}|P(z)|.$
(3.8)

This completes the proof of theorem 1.4. ∎

Acknowledgement. The authors are highly grateful to the referees for careful
reading and their valuable suggestions. The second and third authors are
supported by Council of Scientific and Industrial Research, New Delhi.

## References

  * [1] A. Aziz, Inequalities for the polar derivative of a polynomial, J. Aprox. Theory 55 (1988), 183-193.
  * [2] Robert B. Gardner and N. K. Govil, Inequalities concerning $L^{p}$ norm of a polynomial and its derivative, J. Math. Anal. Appl., 179 (1993), 208-213.
  * [3] Robert B. Gardner and N. K. Govil, An $L^{p}$ inequality for a polynomial and its derivative, J. Math. Anal. Appl., 193 (1995), 490-496.
  * [4] N. K. Govil and G. Labelle, On Bernstein’s inequality, J. Math. Anal. Appl., 126, (1987), 494-500.
  * [5] P. D. Lax, Proof of a conjecture of P. Erdös on the derivative of a polynomial, Bull. Amer. Math. Soc., 50 (1944), 509-513
  * [6] N.K. Govil and Q.I. Rahman, Functions of exponential type not vanishing in a half plane and related polynomials, Trans. Amer. Math. Soc., 137 (1969), 501-517.
  * [7] M. A. Malik, On the derivative of a polynomial, J. Lond. Math. Soc., Second Series 1 (1969), 57-60.
  * [8] G. V. Milovanovic, D. S. Mitrinovic and Th. M. Rassias, Topics in Polynomials: Extremal Properties, Inequalities, Zeros, World scientific Publishing Co., Singapore, (1994).
  * [9] Q. I. Rahman and G. Schmessier, Analytic theory of polynomials, Claredon Press, Oxford, 2002.
  * [10] A. C. Schaffer, Inequalities of A. Markoff and S. Bernstein for polynomials and related functions, Bull. Amer. Math. Soc., 47 (1941), 565-579.


# Mode-Coupling-Driven Frequency Stabilization in Semiconductor Lasers with
Bragg Grating Waveguide

M. R. Mahani  Igor A. Nechepurenko  Yasmin Rahimof  Andreas Wicht  Ferdinand-
Braun-Institut (FBH), Leibniz-Institut für Höchstfrequenztechnik, Gustav-
Kirchhoff-Straße 4, 12489, Berlin, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Precisely stabilizing laser frequency is crucial for advancing laser
technology and unlocking the full potential of various quantum technologies.
Here, we propose a compact device for stabilizing frequency of a semiconductor
laser through mode coupling effects, which provides enhanced sensitivity. Our
proposed architecture features a main ridge waveguide with a Bragg grating,
flanked by two curved ridge waveguides. This configuration exhibits an optical
phenomenon characterized by a transmission crossing at the wavelength of the
Bragg grating. Using particle swarm optimization strategy and employing
efficient figures of merit, we achieve a high transmission crossing. The
observed asymmetric transmission crossing not only holds the promise for an
efficient and compact on-chip laser frequency stabilizer, but also fosters the
development of novel sensing platforms with heightened sensitivity.

##  1 Introduction

Integrated photonics has emerged as a transformative field, revolutionizing
the control of light on a nanoscale platform [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].
One of the most captivating aspects of this domain is the ability to engineer
complex optical structures that offer unprecedented functionalities [11, 12,
13, 14]. In this context, the interplay between waveguide (WG) configurations
and periodic structures, such as Bragg gratings (BGs), has led to advancements
in signal processing [15, 16], modulation [17, 18], and sensing applications
[19]. The effective use of BGs in WG configurations for all-optical signal
processing applications has been demonstrated [20]. The modulation
capabilities of BGs integrated with WGs, exhibiting their potential for high-
speed data modulation has been shown [21]. Moreover, the sensing capabilities
of these structures have been explored [22], where BGs within WGs were
employed for label-free biosensing with enhanced sensitivity [23].

Tilted BGs has recently attracted attention due to their various applications
which include (i) mode coupling and conversion allowing for various signal
processing and modulation techniques [24, 25, 26], (ii) sensitivity to
external perturbations allowing for various sensing applications, including
temperature and strain sensors and chemical or biological sensors [27, 28, 29,
30, 31, 32], (iii) wavelength filtering and tuning allowing for precise
control of the wavelengths of transmitted or reflected light [33].

These gratings, characterized by a tilt angle, offer control over mode
coupling and propagation properties, rendering them indispensable in the realm
of integrated photonics. By inducing controlled variations in the refractive
index along the WG, tilted BGs enable efficient energy transfer between guided
modes of different orders.

Stable laser sources are essential for various quantum optical experiments,
including quantum communication, quantum computation, and quantum sensing [34,
35, 36, 37]. Here, we introduce an optimized on-chip frequency stabilizer
design that leverages the phenomenon of asymmetric transmission crossing
within tilted BG WGs. By utilizing this interplay, the device aims to provide
a robust and compact laser stabilization solution for these quantum optical
applications. The configuration consists of a primary narrow ridge WG
comprising a tilted BG. Positioned on either side of the main WG are curved
narrow ridge WGs, introducing an inherent asymmetry to the system. When light
is launched into the main WG and its transmission is monitored in the flanking
curved WGs, a pronounced transmission crossing is observed exactly at the
Bragg resonance wavelength of the gratings. Through the integration of
particle swarm optimization (PSO) algorithms and finite-difference time-domain
(FDTD) simulation, we have enhanced this effect.

This paper is structured as follows: Section 2 provides an overview of the
theoretical underpinnings of the proposed device, the optimization algorithms
used to maximize the effect, and a discussion of the results, focusing on the
observed asymmetric transmission crossing and its implications. Finally,
Section 3 concludes the paper by highlighting the potential applications of
this phenomenon in integrated photonics devices.

##  2 Simulations and Optimization

###  2.1 Structure

The main narrow ridge WG, with a width of 250 nm and length of either 600 or
200 $\mu m$, serves as the backbone of the structure (Fig. 1). The main WG
comprises 72 tilted BG grooves with a period of 1 $\mu m$ and a refractive
index of 3.5. These grooves are designed to enable frequency selective
reflection at a Bragg wavelength. In addition, the tilt of the grooves allows
an excitation of a higher order mode. A combination of a symmetrical
fundamental mode and anti-symmetrical first order mode introduces an asymmetry
to the field distribution in the main waveguide. There are also two curved WGs
with different distance to the main WG, listed in the table 1. The different
distance to the main WG imparts an asymmetry to the system. We refer to the
one on the left side of the prorogation direction of light as the left WG and
the other one as the right WG, and light is detected at the end of the side
WGs. The schematic of the device’s concept is shown in Fig.2(a). This concept
is implemented and shown in Fig.2(b), where various parts are introduced.

Figure 1:  (Colour online) The 3D schematic structure of the device. The main
ridge waveguide with a Bragg grating, flanked by two curved ridge waveguides.
Figure 2:  (Colour online) (a) The schematic of the device’s concept, (b) full
device, showing various parts; main WG, curved WGs, source, Bragg grating
section, and the monitors.

Our investigation focuses on the transmission characteristics of this
integrated photonics device. Light is launched from one end of the main WG,
adjacent to the BGs, and the transmitted light is monitored in the left and
right WGs. The intensity of the signal in the right and left WGs is equal when
the frequency of the incident light coincides with the frequency at which the
transmission crossing occurs. As the frequency of incident light detunes, the
coupling of modes in the main and side WGs changes. Depending on the sign of
the change, the feedback signal can be created to tune the device back to the
Bragg reflectance peak. Using optimization described below, we have achieved a
transmission crossing in the spectral domain, occurring precisely at the Bragg
resonance wavelength of the gratings.

###  2.2 Optimization

The simulation of the proposed integrated photonic device, was conducted using
ANSYS Optics (Lumerical) software. The geometry of the device, including the
main WG, tilted BGs and curved WGs was defined in the layout editor using
Lumerical script. A simulation region encompassing the entire device structure
was established. A broadband light source was defined in order to excite a
fundamental mode in the desired spectral range for analysis.

Lumerical’s FDTD solver was chosen for the simulation due to its computational
efficiency [38, 39, 40, 41]. To achieve accurate and reliable results the mesh
for the simulation was 20 nm (wavelength/50). Perfectly matched layer (PML)
was used to absorb scattered light. The simulation time was limited by an
auto-shutoff value which was set to be $10^{-5}$ fraction of the energy left
in the system. The reflectance was calculated for the fundamental mode by
means of the mode projection [42]. We used the built-in tools provided by
ANSYS Optics simulation environment to collect and analyze the data.

The simulation procedure to optimize the proposed device was divided into two
stages. The first set of optimization aimed to achieve 25% reflectivity at the
Bragg resonance (1063 nm) and the second set aimed to maximize transmission
crossing ratio or difference. With such an approach, we split our optimization
problem into two stages, each with smaller dimension. This is possible due to
the low dependency of the Bragg reflectance on its tilt and the geometry of
the side WGs. There is no dependency of the resonance wavelength on the
tilting angle up to 15 degree and the amplitude of the resonance drops by
approximately 10% of the maximum.

For the first optimization step, we choose the figure of merit (FOM) as the
absolute value of the difference between the actual and desired reflectivity
at 1063 nm. The optimization of the BG is rather straightforward, given a good
analytical approximation for the period of the grating. The period is defined
as $P=\frac{{N\lambda}}{2{n_{eff}}}$, where $\lambda$ is the desired
wavelength, $N$ is the order of the grating and $n_{eff}$ is an effective
refractive index. The bandwidth of the Bragg resonance is very narrow, so the
variation for the period should be of the order of 1% from the analytical
approximation. Refractive index variation determine the strength of the
interaction of the WG mode with the grating. The number of the grooves also
determines the strength of the reflection, and in general it can be fixed,
when the refractive index is optimized. We implemented a PSO algorithm in
connection with Lumerical software to optimize the period, the number of
gratings and the refractive index.

For the second set of optimization, to maximize the crossing ratio and
difference at exactly Bragg resonance, we chose seven parameters and FOM that
are described in the section 2.3. The set of seven parameters, which
collectively define the device’s geometry and structural characteristics, are:
BG tilt angle (1), width of the main WG (1), width of the side WGs (2), length
of the side WGs (1), distance between the left/right WGs with the main WG (we
refer to these as left/right gaps) (2). The selection of the parameters is
determined by the need to create an asymmetric coupling to the side WGs. Thus
the geometry of the side WGs should be different.

The PSO algorithm simulates the behavior of a swarm of particles in a multi-
dimensional search space. Each particle, representing a potential solution
with a set of parameter values, collectively explores the solution space to
find the optimal configuration. The algorithm iteratively updates particle
positions over multiple generations, where a generation refers to a specific
iteration in the algorithm’s optimization process. During each generation, the
particle positions are updated based on their historical best and the swarm’s
best solution. The objective function is evaluated for each particle, and the
algorithm converges towards the optimal parameter values. The optimization
process continues until convergence criteria are met.

For the device optimization in the 7-dimensional parameter space, we used 50
particles (simulations) in each generation. We tried various generation sizes
(number of particles in a generation). For larger size it requires many more
simulation to converge, and smaller size is not sufficient to cover the entire
parameter space efficiently. However, the results are not very sensitive to
the exact number of the particles. The optimization procedure converges
relatively quickly, after roughly 22 generations. We ran the simulations up to
27 generations ($27\times 50=1350$ simulations) to ensure that the results are
truly the optimum. We also repeated this procedure multiple times with
different generation size, since the particles are generated randomly at the
start of each optimization.

###  2.3 Equations for the figure of merit

The choice of the FOM is very important to obtain the desired outcome. For the
second set of optimizations (optimizing transmission crossing ratio or
difference), we had to try different FOMs to find the one that yields the best
result.

For crossing ratio, we defined the following FOM (We refer to it as the first
FOM) to evaluate the performance of each particle in PSO, quantifying the
degree of asymmetric transmission crossing achieved by the device. The
objective function takes into account the intensity imbalance between the left
and right WGs at the Bragg resonance wavelength. The parameters in Eq.1 are:
$\lambda_{0}=1.064\times 10^{-6}$, $\lambda_{1}=1.0635\times 10^{-6}$ and
$\lambda_{2}=1.0645\times 10^{-6}$.

$\begin{split}\text{FOM}=(\min(x)-1)+(\min(x)>1.05)\cdot\min(\\{T_{1}(i_{1}),T_{2}(i_{2})\\}),\\\
i_{1}=\min\left(\left|\lambda-\lambda_{2}\right|\right),i_{2}=\min\left(\left|\lambda-\lambda_{1}\right|\right),\\\
x=\left[\frac{T_{1}(i_{1})}{T_{2}(i_{1})+1\times
10^{-6}},\frac{T_{2}(i_{2})}{T_{1}(i_{2})+1\times 10^{-6}}\right].\end{split}$
(1)

The goal was to maximize this imbalance, thereby maximizing the
distinguishability of spectral shift.The range of the parameters and the
optimized values are listed in the table 1, marked with the figure number,
Figs.3(a) and Figs.3(b).

Table 1: Optimized parameters and their respective range. Parameters for 1st optimization | Parameter range | Optimized parameters
---|---|---
Period (P) | [0.99, 1.01] $\mu m$ | 1.0 $\mu m$
Number of gratings (N) | [25, 72] | 72
Refractive index (n) | [3.4, 3.5] | 3.482
Parameters for 2nd optimization (main WG with the length of 600 $\mu$m) |  | Fig.3(a) - Fig.4(a)
Rotation angle ($\theta$) | [0, 15] ∘ | 7.46 ∘ \- 0 ∘
Width of the groovs ($w_{x}$) | [80, 120] $nm$ | 120 $nm$ \- 120 $nm$
Left/right WG length ($l$) | [10, 300] $\mu m$ | 267 $\mu m$\- 271 $\mu m$
Left gap ($d_{l}$) | [50, 200] $nm$ | 200 $nm$ \- 66 $nm$
Right gap ($d_{r}$) | [50, 200] $nm$ | 94 $nm$ \- 165 $nm$
Width of the left WG ($w_{l}$ ) | [100, 400] $nm$ | 368 $nm$ \- 225 $nm$
Width of the right WG ($w_{r}$) | [100, 400] $nm$ | 400 $nm$ \- 300 $nm$
Parameters for 2nd optimization (main WG with the length of 200 $\mu$m) |  | Fig.3(b) - Fig.4(b)
Rotation angle ($\theta$) | [0, 15] ∘ | 15 ∘\- 4.54 ∘
Width of the groovs ($w_{x}$) | [80, 120] $nm$ | 110 $nm$ \- 120 $nm$
Left/right WG length ($l$) | [10, 30] $\mu m$ | 18.5 $\mu m$\- 30 $\mu m$
Left gap ($d_{l}$) | [50, 200] $nm$ | 59 $nm$ \- 142 $nm$
Right gap ($d_{r}$) | [50, 200] $nm$ | 62 $nm$ \- 50 $nm$
Width of the left WG ($w_{l}$ ) | [100, 300] $nm$ | 114 $nm$\- 300 $nm$
Width of the right WG ($w_{r}$) | [100, 300] $nm$ | 100 $nm$ \- 226 $nm$

The optimal structure demonstrates a large crossing ratio (16.8 for Fig.3(b)).
The device has a large change of the signals detected between the left and
right WGs. However, the overall signal is small in comparison with the
transmitted light. In order to optimize for a higher transmission of light to
the side WGs, we developed another FOM. The second FOM was used for Fig.4.

In this case, for crossing difference, we defined the FOM as the following. We
first need to calculate the transmission difference $dT$ and its maximal
variation $x_{1}$.

$\begin{split}dT&=T_{1}-T_{2},\\\
x_{1}&=\min\left(-\min(dT),\max(dT)\right).\\\ \end{split}$ (2)

Then we wrote a code that iterates through the elements of the $dT$ array to
find points where the sign of adjacent elements changes. It keeps track of
these positions in the crossings array. The code then checks if the crossing
is greater than 0.01, it calculates FOM using a specific formula involving
$i_{0}$ and x. Otherwise, it calculates FOM based on $x_{1}$.

$\left\\{\begin{array}[]{l}FOM=10\times\frac{1}{{\min\left({\left|{{i_{cross}}-{i_{0}}+{{10}^{-3}}}\right|}\right)}},\,\,\,N>0.\\\
FOM={x_{1}},\,\,\,N=0.\end{array}\right.$ (3)

In the equation above $N$ is the number of crossing in the given wavelength
range, $i_{cross}$ is the index of the wavelength of the crossing closest to
the index of the desired wavelength $i_{0}$. The optimized parameters and
their range are listed in the table 1, marked with the figure number.

When we use this FOM, the optimized transmittance of light to side WGs becomes
higher. The crossing also agrees well with the Bragg resonance (Fig.4). Such
device is efficient to create a differential feedback signal.

(a)  (b)

Figure 3:  (Colour online) Transmission crossing and Bragg resonance for the
first FOM. (a) Main WG with the length of 600 $\mu$m (a crossing ratio of 2.9
and difference of 5.4$\times 10^{-4}$ can be achieved at 1nm away from Bragg
resonance). (b) Main WG with the length of 200 $\mu$m (a crossing ratio of
16.8 and difference of 8.4$\times 10^{-5}$ can be achieved at 1nm away from
Bragg resonance). The reflectance is shown on the right axis (red color). The
transmittance in the left/right WG (solid/dashed blue line) are shown on the
left axis.

###  2.4 Discussions

The device architecture presented in this study, featuring tilted BGs
integrated with curved WGs, has demonstrated an asymmetric transmission
crossing at the Bragg resonance wavelength. By precisely tuning the
parameters, the device can be used to select and filter specific spectral
bands, a fundamental function in spectrometry or stabilizing laser frequency.

The transmission characteristics of the device are closely related to the
spectral properties of the incident light. When light is launched into the
main WG and monitored in the left/right WGs, the occurrence of a transmission
crossing at the Bragg resonance wavelength signifies a specific wavelength’s
presence in the incident light. By measuring the wavelength at the crossing,
the device can be used to analyze the spectral content of the input light.

(a)  (b)

Figure 4:  (Colour online) Transmission crossing and Bragg resonance for the
second FOM. (a) Main WG with the length of 600 $\mu$m (a crossing ratio of 5.7
and difference of 5.1$\times 10^{-2}$ can be achieved at 1nm away from Bragg
resonance). (b) Main WG with the length of 200 $\mu$m (a crossing ratio of 1.1
and difference of 1.4$\times 10^{-2}$ can be achieved at 1nm away from Bragg
resonance). The reflectance is shown on the right axis (red color). The
transmittance in the left/right WG (solid/dashed blue line) are shown on the
left axis.

A high ratio of crossing transmission ($16.8$ in Fig. 3(b)) suggests that the
device can effectively identify specific wavelengths and is effective at
separating wavelengths of interest from others. This can be important in
applications where it is needed to distinguish and identify specific spectral
lines in the presence of background noise.

On the other hand, if we want to increase the efficiency to create a
differential feedback signal, we need to increase the transmittance of light
in the side WGs. A weak signal ($10^{-4}-10^{-3}$ in Fig. 3) could be
difficult to detect. While it might not be highly significant for quantifying
intensity differences, it can still be useful for precise measurements in
applications where the goal is to quantify subtle differences in spectral
features. The small difference is important when we work with very faint
spectral lines. However, this difference can be separately optimized using
second FOM, as shown in Fig.4(a). Here the crossing difference has improved to
reach 0.05 while the ratio is 5.7.

To gain an understanding of the device’s behavior, we examine the electric
field distribution in 2D space across the three WGs. We consider three cases,
where the light pulse with the narrow spectra (few nms wide) is centered on
three distinct wavelengths, Figs.5(a) $\lambda=1063$ nm, 5(b)
$\lambda=1061.5$, and 5(c) $\lambda=1064.5$. If the incident light is centered
on $\lambda=1063$, we observe the equal intensity in the side WGs,
highlighting the transmission crossing at the BG resonance wavelength. At this
wavelength, the transmission is equal in the side WGs, thus the intensity of
the field is roughly the same (Fig.5(a)).

Figure 5:  (Colour online) Electric field distribution in 2D space for a
specific time (2000 fs after the light was launched) in the three WGs for
three different wavelengths of the incident light. (a) $\lambda=1063$ nm, (b)
$\lambda=1061.5$ nm, (c) $\lambda=1064.5$ nm.

Further investigation reveals that as the incident wavelength deviates from
the resonance point, the field distribution changes on the left WG. Figure
4(a) illustrates that the transmission is notably higher in the left WG for
shorter wavelengths ($\lambda=1061.5$). This asymmetry indicates the
wavelength-dependent propagation of the electric field within the device
(Fig.5(b)). Conversely, as the incident wavelength extends to longer values
($\lambda=1064.5$), a reversal in the transmission and respectively the
electric field distribution occurs, with the right WG exhibiting a higher
intensity (Fig.5(c)). We should keep in mind that the distance between the
main WG and the left and right WG is different (66 nm and 165 nm,
respectively). The coupling of the modes of the left WG with the modes of the
central WG is much stronger than that of the right WG. The strong mode-
coupling translates into short coupling length. This makes the system very
sensitive to the optical length of the structure. As a result, the strong
wavelength dependence of the transmission to the top WG is observed on
Fig.4(a). On the contrary, the coupling of the modes of the right WG with the
mode of the central WG is relatively weak, as a result the transmission of
light to the right WG is almost constant.

The wavelength-dependent modulation in the electric field distribution across
the WGs, in agreement with transmission crossing, elucidates the sensitivity
of the proposed architecture, providing potential applications in laser
frequency stabilization and sensing platforms. Since it is capable of
filtering specific spectral bands, it can be used to isolate and stabilize a
particular frequency. The setup’s ability to detect specific wavelengths (high
ratio of crossing transmission) is important for identifying the laser’s
emission wavelength, ensuring that it remains within the desired range.

When the laser is operating at the wavelength close to the Bragg resonance
wavelength (1063 nm), the device will exhibit a transmission crossing at this
wavelength. This crossing is then used as a feedback signal. When the laser
wavelength drifts away from the Bragg resonance, the transmission
characteristics of the device will change. This change in transmission can be
used as an error signal. Using a feedback control system, once the
transmission characteristics change due to the drift, the control system
brings it back to the the Bragg resonance wavelength. We would like to
emphasize that, the setup’s tunable grating parameters provide flexibility,
allowing it to adapt to different laser sources and frequency stabilization
requirements.

It is important to measure the sensitivity of the results with respect to the
changes in the parameters. In order to analyze the sensitivity of the proposed
device, we varied the values of six most relevant parameters. The parameters
are the left/right gap ($d_{l}$ and $d_{r}$), refractive index ($n$), length
of the side WGs ($l$), width of the main WG and the tilting angle ($\theta$).
We varied these parameters and compared the new transmission crossing
wavelength $\lambda_{new}$ with the reference or optimized wavelength
($\lambda_{ref}$). The sensitivity is then calculated by taking the difference
between transmission crossing wavelengths for $1\%$ change in the parameter
values ($s=\lambda_{ref}-\lambda_{new}$). Calculated sensitivities for the
left and right gap are $s=-0.0312$ nm and $s=0.1863$ nm, respectively. The
sensitivities for refractive index, length of the side WG, width of the main
WG and the tilt angle are $s=0.0012$ nm , $s=-0.0046$ nm, $s=0.1149$ nm and
$s=0.0001$ nm, respectively. The most sensitive parameters in our model are
the the right gap and then the width of the main WG. The tilt angle seems to
be the least sensitive parameter.

The device’s performance may also be sensitive to environmental factors such
as temperature fluctuations and vibrations. Ensuring stability in real-world
conditions might be required. Fabrication imperfections and variations in the
device’s physical parameters (WG dimensions, grating properties) during
manufacturing can also impact the device’s performance.

##  3 Conclusion

The device we proposed here shows an asymmetric transmission crossing at the
Bragg resonance wavelength. Through a combination of PSO algorithm and
Lumerical software, we optimized this crossing at the Bragg resonance. We
highlighted the details of our optimization approach and various FOMs and
discussed the ones that gave rise to desired results and made the simulations
converge fast.

Our investigation revealed the device’s capacity to enhance spectral
discrimination and sensitivity, leveraging the distinctive intensity imbalance
at the asymmetric crossing. We obtained a high ratio of crossing transmission
($16.8$) and a large difference in transmitted intensities ($0.05$). We
envisage applications in fields ranging from advanced optical stabilizing
laser frequency to sensing and quantum applications.

Author contributions: MRM, IN, YR, and AW conceived and planned the research.
MRM, IN and YR carried out the simulations and optimizations. MRM, IN wrote
the manuscript with input from all authors. All authors contributed to the
interpretation of the results, provided feedback, and helped shape the
research.

Research funding: This work was supported by the VDI Technologiezentrum GmbH
with funds provided by the Federal Ministry of Education and Research under
grant no. 13N14906 and the DLR Space Administration with funds provided by the
Federal Ministry for Economic Affairs and Climate Action (BMWK) under Grant No
50WK2272.

Conflict of interest statement: The authors declare no conflicts of interest
regarding this article.

Data Availability: The data underlying this study are available from the
corresponding author upon request.

## References

  * [1] D. A. Miller, “Device requirements for optical interconnects to silicon chips,” Proceedings of the IEEE 97, 1166–1185 (2009).
  * [2] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, “Integrated photonic quantum technologies,” Nature Photonics 14, 273–284 (2020).
  * [3] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund _et al._ , “Deep learning with coherent nanophotonic circuits,” Nature photonics 11, 441–446 (2017).
  * [4] D. A. Miller, “Attojoule optoelectronics for low-energy information processing and communications,” Journal of Lightwave Technology 35, 346–396 (2017).
  * [5] D. A. Miller, “Optical interconnects to electronic chips,” Applied optics 49, F59–F70 (2010).
  * [6] R. Jones, P. Doussiere, J. B. Driscoll, W. Lin, H. Yu, Y. Akulova, T. Komljenovic, and J. E. Bowers, “Heterogeneously integrated inp$\backslash$/silicon photonics: fabricating fully functional transceivers,” IEEE Nanotechnology Magazine 13, 17–26 (2019).
  * [7] J. Wang, M. M. Sanchez, Y. Yin, R. Herzer, L. Ma, and O. G. Schmidt, “Silicon-based integrated label-free optofluidic biosensors: latest advances and roadmap,” Advanced Materials Technologies 5, 1901138 (2020).
  * [8] R. Soref, “Silicon photonics: a review of recent literature,” Silicon 2, 1–6 (2010).
  * [9] R. Soref, “Tutorial: Integrated-photonic switching structures,” Apl Photonics 3 (2018).
  * [10] Y.-J. Hung, C.-C. Yin, Z.-R. Wu, C.-Y. Lin, C.-W. Kao, T.-C. Kao, and C.-W. Huang, “Near-infrared spectroscopy using period-chirped si/sio/sio 2-based guided mode resonance filter,” Optics Letters 47, 561–564 (2022).
  * [11] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nature Photonics 12, 659–670 (2018).
  * [12] P. Cheben, R. Halir, J. H. Schmid, H. A. Atwater, and D. R. Smith, “Subwavelength integrated photonics,” Nature 560, 565–572 (2018).
  * [13] R. Halir, A. Ortega-Monux, D. Benedikovic, G. Z. Mashanovich, J. G. Wangüemert-Pérez, J. H. Schmid, I. Molina-Fernandez, and P. Cheben, “Subwavelength-grating metamaterial structures for silicon photonic devices,” Proceedings of the IEEE 106, 2144–2157 (2018).
  * [14] N. I. Zheludev and Y. S. Kivshar, “From metamaterials to metadevices,” Nature materials 11, 917–924 (2012).
  * [15] M. Burla, L. R. Cortés, M. Li, X. Wang, L. Chrostowski, and J. Azaña, “Integrated waveguide bragg gratings for microwave photonics signal processing,” Optics express 21, 25120–25147 (2013).
  * [16] W. Shi, X. Wang, C. Lin, H. Yun, Y. Liu, T. Baehr-Jones, M. Hochberg, N. A. Jaeger, and L. Chrostowski, “Silicon photonic grating-assisted, contra-directional couplers,” Optics express 21, 3633–3650 (2013).
  * [17] B. Stern, H. Chen, K. Kim, and N. K. Fontaine, “Large dispersion silicon bragg grating for full-field 40-gbd qpsk phase retrieval receiver,” Journal of Lightwave Technology 40, 7358–7363 (2022).
  * [18] Q. Ye, R. Huang, Q. Xu, H. Cai, R. Qu, and Z. Fang, “Numerical investigation of ultrashort complex pulse generation based on pulse shaping using a superstructure fiber bragg grating,” Journal of Lightwave Technology 27, 2449–2456 (2009).
  * [19] L. R. Chen, M.-I. Comanici, P. Moslemi, J. Hu, and P. Kung, “A review of recent results on simultaneous interrogation of multiple fiber bragg grating-based sensors using microwave photonics,” Applied Sciences 9, 298 (2019).
  * [20] C. Sima, J. Gates, H. Rogers, B. Snow, C. Holmes, M. Zervas, and P. Smith, “All-optical signal processing using planar bragg gratings,” in _Photonics and Optoelectronics Meetings (POEM) 2011: Optoelectronic Devices and Integration,_ vol. 8333 (SPIE, 2012), pp. 64–68.
  * [21] V. Pilipovich, A. Esman, I. Goncharenko, and V. Kuleshov, “High-speed continuous tuneable fibre and waveguide lasers with controllable bragg grating,” Optics communications 203, 289–294 (2002).
  * [22] E. Luan, H. Yun, S. H. Lin, K. C. Cheung, L. Chrostowski, and N. A. Jaeger, “Phase-shifted bragg grating-based mach-zehnder interferometer sensor using an intensity interrogation scheme,” in _Optical Fiber Communication Conference,_ (Optica Publishing Group, 2020), pp. M1C–5.
  * [23] E. Luan, H. Yun, M. Ma, D. M. Ratner, K. C. Cheung, and L. Chrostowski, “Label-free biosensing with a multi-box sub-wavelength phase-shifted bragg grating waveguide,” Biomedical optics express 10, 4825–4838 (2019).
  * [24] B. Wang, J. Jiang, and G. P. Nordin, “Compact slanted grating couplers,” Optics express 12, 3313–3326 (2004).
  * [25] D. Dai, Y. Tang, and J. E. Bowers, “Mode conversion in tapered submicron silicon ridge optical waveguides,” Optics express 20, 13425–13439 (2012).
  * [26] Y. Yang, K. Chen, W. Jin, and K. S. Chiang, “Widely wavelength-tunable mode converter based on polymer waveguide grating,” IEEE Photonics Technology Letters 27, 1985–1988 (2015).
  * [27] R. Gamal, Y. Ismail, and M. A. Swillam, “Optical biosensor based on a silicon nanowire ridge waveguide for lab on chip applications,” Journal of Optics 17, 045802 (2015).
  * [28] N. Saha and A. Kumar, “Highly sensitive temperature sensor based on a long period grating inscribed metal clad ridge waveguide with pdms surrounding,” OSA Continuum 2, 946–956 (2019).
  * [29] Z. Huang, N. Yang, J. Xie, X. Han, X. Yan, D. You, K. Li, M. I. Vai, K. W. Tam, T. Guo _et al._ , “Improving accuracy and sensitivity of a tilted fiber bragg grating refractometer using cladding mode envelope derivative,” Journal of Lightwave Technology 41, 4123–4129 (2023).
  * [30] O. Butov, K. Tomyshev, I. Nechepurenko, A. Dorofeenko, and S. Nikitov, “Tilted fiber bragg gratings and their sensing applications,” Uspekhi Fiz. Nauk 192, 1385–1398 (2022).
  * [31] E. D. Chubchev, K. A. Tomyshev, I. A. Nechepurenko, A. V. Dorofeenko, and O. V. Butov, “Machine learning approach to data processing of tfbg-assisted spr sensors,” Journal of Lightwave Technology 40, 3046–3054 (2022).
  * [32] H. Jean-Ruel and J. Albert, “Recent advances and current trends in optical fiber biosensors based on tilted fiber bragg gratings,” TrAC Trends in Analytical Chemistry p. 117663 (2024).
  * [33] Y.-J. Hung, K.-H. Lin, C.-J. Wu, C.-Y. Wang, and Y.-J. Chen, “Narrowband reflection from weakly coupled cladding-modulated bragg gratings,” IEEE Journal of Selected Topics in Quantum Electronics 22, 218–224 (2015).
  * [34] L.-M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, “Long-distance quantum communication with atomic ensembles and linear optics,” Nature 414, 413–418 (2001).
  * [35] T. Aoki, R. Sreekantham, B. Sahoo, B. Arora, A. Kastberg, T. Sato, H. Ikeda, N. Okamoto, Y. Torii, T. Hayamizu _et al._ , “Quantum sensing of the electron electric dipole moment using ultracold entangled fr atoms,” Quantum Science and Technology 6, 044008 (2021).
  * [36] M. Lezius, T. Wilken, C. Deutsch, M. Giunta, O. Mandel, A. Thaller, V. Schkolnik, M. Schiemangk, A. Dinkelaker, A. Kohfeldt, A. Wicht, M. Krutzik, A. Peters, O. Hellmig, and H. Duncker, “Space-borne frequency comb metrology,” Optica 3, 1381–1387 (2016).
  * [37] D. Becker, M. D. Lachmann, S. T. Seidel, H. Ahlers, A. N. Dinkelaker, J. Grosse, O. Hellmig, H. Müntinga, V. Schkolnik, T. Wendrich, A. Wenzlawski, B. Weps, R. Corgier, T. Franz, and N. Gaaloul, “Space-borne bose–einstein condensation for precision interferometry,” Nature 562, 391–395 (2018).
  * [38] X. Wang, Y. Wang, J. Flueckiger, R. Bojko, A. Liu, A. Reid, J. Pond, N. A. Jaeger, and L. Chrostowski, “Precise control of the coupling coefficient through destructive interference in silicon waveguide bragg gratings,” Optics letters 39, 5519–5522 (2014).
  * [39] I. Nechepurenko, Y. Rahimof, M. R. Mahani, S. Wenzel, and A. Wicht, “Finite-difference time-domain simulations of surface bragg gratings,” in _2023 International Conference on Numerical Simulation of Optoelectronic Devices (NUSOD),_ (IEEE, 2023), pp. 3–4.
  * [40] M. R. Mahani, Y. Rahimof, S. Wenzel, I. Nechepurenko, and A. Wicht, “Data-efficient machine learning algorithms for the design of surface bragg gratings,” ACS Applied Optical Materials 1, 1474–1484 (2023).
  * [41] M. R. Mahani, I. Nechepurenko, Y. Rahimof, and A. Wicht, “Designing rectangular surface bragg gratings using machine learning models,” in _2023 International Conference on Numerical Simulation of Optoelectronic Devices (NUSOD),_ (IEEE, 2023), pp. 69–70.
  * [42] A. W. Snyder, J. D. Love _et al._ , _Optical waveguide theory_ , vol. 175 (Chapman and hall London, 1983).


# Perspective-consistent multifocus multiview 3D reconstruction of small
objects

Hengjia Li  College of Engineering and Computer Science
Australian National University
Canberra, Australia
<EMAIL_ADDRESS>Chuong Nguyen  Imaging and Computer Vision Group
CSIRO DATA61
Canberra, Australia
<EMAIL_ADDRESS>
###### Abstract

Image-based 3D reconstruction or 3D photogrammetry of small-scale objects
including insects and biological specimens is challenging due to the use of
high magnification lens with inherent limited depth of field, and the object’s
fine structures and complex surface properties. Due to these challenges,
traditional 3D reconstruction techniques cannot be applied without suitable
image pre-processings. One such preprocessing technique is multifocus stacking
that combines a set of partially focused images captured from the same viewing
angle to create a single in-focus image. Traditional multifocus image capture
uses a camera on a macro rail. Furthermore, the scale and shift are not
properly considered by multifocus stacking techniques. As a consequence, the
resulting in-focus images contain artifacts that violate perspective image
formation. A 3D reconstruction using such images will fail to produce an
accurate 3D model of the object. This paper shows how this problem can be
solved effectively by a new multifocus stacking procedure which includes a new
Fixed-Lens Multifocus Capture and camera calibration for image scale and
shift. Initial experimental results are presented to confirm our expectation
and show that the camera poses of fixed-lens images are at least 3-times less
noisy than those of conventional moving lens images.

###### Index Terms:

Fixed-Lens Multifocus Capture, macro imaging, multifocus stacking, perspective
image formation, multiview stereo, 3D reconstruction, image-based 3D
reconstruction, insects, small objects, small specimens

##  I Introduction

3D reconstruction of small objects including insects and biological specimens
is challenging due to the use of high magnification lens with limited depth of
field, fine features and complex surface properties. Recent advancements [7,
9] show that photogrammetry or image-based multiview 3D reconstruction can be
applied with some success to create true-color 3D models of small specimens.
Solutions enabling images of small specimens to be reconstructed in 3D
include: a two-axis turntable combined with macro rail and macro photography
to capture multifocus multiview images [7], calibration target [7, 9],
multifocus image stacking [14], scale-shift calibration and automatic
background masking [9], and modern photogrammetry software [8, 10, 11].
Despite multiple techniques to tackle different issues of image-based 3D
reconstruction of small specimens, obtaining an accurate 3D model is still
difficult. Particularly, camera pose estimation from stacked in-focus images
is still unreliable, and it causes poor quality 3D reconstruction. To better
estimate camera poses, a special calibration object, such as a textured
sphere, needs to be used [9] and camera poses are precomputed. However, this
workaround solution restricts to predetermined camera positions and relies on
the high repeatability of motorized hardware.

In this paper, we identify multifocus image stacking as one major source of
error causing unreliable 3D reconstruction. Multifocus stacking combines
partially in-focus images by automatically detecting in-focus regions from the
images and merging them to generate an in-focus image. The resulting image is
just a good-looking mosaic image which does not accurately represent a
perspective projection. Scale-shift calibration proposed by [9] partially
corrects this problem by scaling and shifting a set of partially in-focus
images before being blended into a stacked in-focus image. However, as the
original image set are captured at different camera lens positions, the
blended in-focus image never accurately represents a perspective projection.
This is illustrated by the bottom of Figure 1. The use of such stacked in-
focus images breaks the major assumption in 3D reconstruction algorithms [8]
that each input image representing a single perspective pinhole projection at
a single camera position and orientation (called camera pose). As a result, a
3D multiview reconstruction using such in-focus images from different viewing
angles fails to generate an accurate 3D model of the object.

This paper proposes an effective solution to capture multifocus images, and
enable multifocus stacking to produce an in-focus image representing a single
perspective projection for accurate 3D reconstruction of small objects. Our
contributions include: a) a newly proposed technique called Fixed-Lens
Multifocus Capture to capture suitable multifocus images, b) a new multifocus
stacking procedure to create an in-focus image that correctly represents a
perspective projection. An overview of the concept of pinhole camera model,
the problem of perspective inconsistency, depth of focus and multifocus
stacking is introduced in section II. Section III describes our proposed
method of multifocus multiview 3D reconstruction: the fixed lens setup, camera
calibration, multifocus stacking based on Guided Filtering, and background
masking. Section IV presents experimental results showing the effectiveness of
the proposed method versus the conventional method. Finally, section V
concludes the paper with our findings.

##  II Perspective image formation and depth of focus effect

###  II-A Pinhole camera and perspective image formation

The firstly invented camera is the pinhole camera as shown in the top of
Figure 1 where the pinhole C is the centre of image formation allowing light
rays passing straight through. Modern cameras use lenses instead of the
pinhole to improve image quality, but it adds the depth of focus where a film
or image sensor is located to capture a clear in-focus image. Due to the image
capture process of a camera, a transformation between the 3D world to 2D image
happens where the depth is missing in the captured image. Furthermore, parts
in the scene appear at different scales or magnifications depending on their
distance to the camera lens. This is demonstrated by the bottom of Figure 1
where the images of the same object are captured at two different distances.
As the camera lens moves, the centre C of image formation moves. This leads to
the changes in the relative distances (and magnification) between parts of the
scene to camera centre. As a result, the change of relative scales of
different parts of the scene causes different perspective distortions.

Figure 1: Top: Pinhole camera principle and perspective projection[12] where
all straight light rays go through this hole (image in public domain). Pinhole
represents the camera centre for image formation. Bottom: the pair of oranges
captured by same camera of 5mm focal length at different distances leading to
different perspective image formations (or distortions).

###  II-B Depth of focus, depth of field, and multifocus stacking

Due to the use of lenses in modern cameras, image quality is significantly
better, but the depth of focus is introduced in the image. As shown in the top
of Figure 2, an image of the object at distance $d_{1}$ forms a image at
distance $d_{2}=\frac{d_{1}f}{d_{1}-f}$, where $f$ is the lens focal length.
If an image sensor is placed at this distance (i.e. at back focal plane), the
image will be clear and in focus.

In reality, an image is considered in focus if the image of a point remains
smaller than a circular dot called the circle of confusion (COC) with an
empirical diameter $\phi_{coc}$. From [13], $\phi_{coc}$ is chosen
approximately 0.1% of the mean of the width and height of the image sensor.
For example, for a 35mm sensor format, $\phi_{coc}$ is chosen to be 0.025mm.
Strictly speaking, $\phi_{coc}$ should be selected as the larger value of the
size of an image pixel and optical resolution of the lens.

The distance range where the image remains in focus is called the depth of
focus $DOFocus$. Using the proportional relationship, one can obtain the depth
of focus as:

$\displaystyle\frac{DOFocus}{2d_{2}}$
$\displaystyle=\frac{\phi_{coc}}{\phi_{a}}$ (1) $\displaystyle\Rightarrow
DOFocus$ $\displaystyle=2\frac{\phi_{coc}}{\phi_{a}}d_{2}$ (2)
$\displaystyle\text{or }DOFocus$
$\displaystyle=2\frac{f_{number}\phi_{coc}}{f}d_{2}$ (3)

where $\phi_{a}$ is the diameter of the lens aperture, and
$f_{number}=\frac{f}{\phi_{a}}$.

This depth of focus translates to the depth of field $DOField$ where objects
stay within so that their image is in-focus. Assuming the circle of confusion
$\phi_{coc}$ is back projected into the scene to the size
$\frac{\phi_{coc}}{M}$ where $M=\frac{d_{2}}{d_{1}}$ is lens magnification,
the proportional relationship gives:

$\displaystyle\frac{DOField}{2d_{1}}$
$\displaystyle=\frac{\phi_{coc}}{M\phi_{a}}$ (4) $\displaystyle\Rightarrow
DOField$ $\displaystyle=2\frac{\phi_{coc}}{M\phi_{a}}d_{1}$ (5)
$\displaystyle\text{or }DOField$
$\displaystyle=2\frac{d_{1}^{2}\phi_{coc}}{d_{2}\phi_{a}}=2\frac{d_{1}f\phi_{coc}}{(d_{1}-f)\phi_{a}}$
(6)

To capture scenes at different distances, the position of the lens and/or the
image sensor needs to be adjusted to put the image into focus.

Figure 2: Top: image formation through a lens of focal length $f$. After
passing the lens, rays parallel to optical axis go through focal point and
rays going through focal point becomes parallel to optical axis. Bottom: image
sensor is used to capture an object image at a position within the depth of
focus where a point source grows into a circle of confusion [13] of diameter
$\phi_{coc}$. A captured image is considered as out of focus if image sensor
is placed outside the depth of focus or the object is located outside of the
depth of field.

For an object with parts at different distances, only some parts are in focus
and the other parts are out-of-focus. This is especially prevalent in macro
imaging where the size of the scene is tens of centimeter or smaller and high
magnification lens is used. Figure 3 shows images captured as a camera (and
lens) moves towards the specimen at constant incremental depth at the same
viewing angle. Multifocus image stacking produced an in-focus image
representing that single viewing angle.

Figure 3: A set of 20 partially out-of-focus images (only 5 of them shown
here) captured at different camera positions are stacked to produce an in-
focus image. The stacked image however does not represent a single perspective
image formation due to moving lens.

Conventional multifocus capture is to move camera and lens together. The step
size for macro rail is equivalent to the depth of field. However, the
recommended macro rail step size is 50% of the depth of field in this case to
allow 50% overlapping between in-focus regions of successive images to help
guide image registration for optimisation-based multifocus stacking
algorithms.

By using such multifocus stacked images captured at different angles, 3D
reconstruction of small objects becomes possible as reported by [7, 9].
However, such in-focus stacked images do not correctly represent perspective
image formation, leading to reconstruction errors or failures. This becomes
more severe when the size of the objects becomes smaller, i.e. a few
millimeters long or smaller.

##  III Multifocus image stacking with fixed perspective image formation

###  III-A Fixed-Lens Multifocus Capture

Conventional camera setup captures multifocus images by moving the camera and
lens together as shown in the top of Figure 4. As the whole camera moves, the
camera centre moves, and perspective view moves with it. To avoid moving the
perspective view of the camera, we proposed that the lens is fixed and only
the image sensor moves during scanning as shown in the bottom of Figure 4. In
this case, the centre of perspective image formation does not change. Although
the image size of the whole scene becomes larger when the film is moving away
from the lens, the relative scale or magnification of different parts of the
scene stays the same.

Figure 4: Conventional multifocus capture (top) and proposed Fixed-Lens
Multifocus Capture (bottom) for multifocus image capture. For the conventional
scanning, the movement of the camera and lens is the same as that of the front
focal plane. For the proposed scanning, the movement of the camera is smaller
than that of front focal plane if the lens magnification is smaller than 1,
and vice versa.

An example of camera lens setup for the proposed Fixed-Lens Multifocus Capture
is shown on the right of Figure 5 as compared to the conventional one (left).
The lens is fixed to the upper frame while the camera body is moved by a macro
rail to capture multiple partially focused images. A rubber duct connects the
lens and the camera body to prevent ambient light and dusts entering the
camera sensor chamber. A small flower is mounted on a 2-axis turntable to
capture images at different pan-tilt angles.

Figure 5: Lens and pan-tilt setup for multifocus multiview image capture of
small objects. The camera is mounted on a macro rail fixed to the lower frame
for multifocus capture. A macro lens is attached to a camera for conventional
moving lens setup (left) or to a stationary frame for our proposed fixed lens
setup (right). For the fixed lens setup, an expandable rubber duct connects
the macro lens and the camera body. The object/specimen is pinned to a fridge
magnet on a pan-tilt turntable for multiview capture.

###  III-B Camera calibration to measure image scale and shift

Images captured at different camera distances from the lens have different
magnifications or scales. Furthermore, the direction of the camera’s movement
is not generally aligned with the lens optical axis, causing a relative shift
in the image. To account for such changes, a camera calibration is performed
to measure the relative scale change and shift from one image to another. This
can be carried out by capturing multifocus images of a known calibration
target placed perpendicular to the optical axis of the lens. The relative
scale and shift can be measured from the relative positions of the control
points (circles) between the images. Once the relative scales of the scene in
partially focused images are determined, image mapping is applied to remove
the scale change and shift between images captured from the same viewing
angle.

Figure 6 shows calibration images captured as a camera moves perpendicularly
towards the calibration target at constant incremental depth at the same
viewing angle, and this camera motion will be repeated during the scan of a
specimen. The image of calibration pattern represents the scale and location
of the object image as the camera moves to capture multifocus images. A small
in-plane shift between images can be observed in additional to the change of
scale.

Figure 6: A set of 9 target images (from 0 to 8) captured at different camera
positions for camera calibration. The camera moves toward the calibration
target and its focal plane passes through the calibration target. Thus, the
captured images are from out-of-focus (blurry) to in-focus (sharp) and become
out-of-focus (blurry) again. The scale of the grids gets smaller as the
cameras comes toward the target. There is also a small shift from top left to
bottom right of the image.

Instead of measuring directly the scale and the in-plane shift between an
image and a reference image as in [9], we compute a homography matrix from the
points correspondingly and use this to transform an image to align with the
reference image before multifocus stacking. Figure 7 shows the calibration
dots (A) are detected as in (B) and compared with a set from its reference
image (C). A homography is computed from point corresponding between that
image and the reference image (usually the one at the middle of the set). The
homography matrix $H$ is defined as a linear transformation between the
positions $X_{i}$ of the detected calibration dots on the i-th image and those
$X_{r}$ on the reference image:

$\displaystyle X_{r}$ $\displaystyle=HX_{i}$ (7)
$\displaystyle\begin{bmatrix}x_{r}\\\ y_{r}\\\ w_{r}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}h_{1}^{T}X_{i}\\\ h_{2}^{T}X_{i}\\\
h_{3}^{T}X_{i}\end{bmatrix}$ (8)

Rewrite this as:

$\displaystyle X_{r}\times
HX_{i}=\begin{bmatrix}y_{r}h_{3}^{T}X_{i}-w_{r}h_{2}^{T}X_{i}\\\
w_{r}h_{1}^{T}X_{i}-x_{r}h_{3}^{T}X_{i}\\\
x_{r}h_{2}^{T}X_{i}-y_{r}h_{1}^{T}X_{i}\end{bmatrix}=0$ (9)

Factorise this to:

$\displaystyle\begin{bmatrix}0&-w_{r}X_{i}^{T}&y_{r}X_{i}^{T}\\\
w_{r}X_{i}^{T}&0&-x_{r}X_{i}^{T}\\\
-y_{r}X_{i}^{T}&x_{r}X_{i}^{T}&0\end{bmatrix}\begin{bmatrix}h_{1}\\\ h_{2}\\\
h_{3}\end{bmatrix}=0$ (10)

The homography matrix can be solved by singular value decomposition of
Equation 10.

Figure 7: Calibration process for image stacking. (A) image of calibration
target with black circular dots. (B) detected dots shown as blue (C) Overlay
of circular dots detected from an image (red) and those from a reference image
(green). A homography matrix is computed from the corresponding points and
used to transform the images captured at a similar camera depth position of
red dots to align with the reference image.

After the captured images from the same viewing angle are transformed to align
to the reference image using the homography obtained from corresponding
calibration images, these transformed images can be stacked or blended to
create an in-focus image.

###  III-C Multifocus image stacking of transformed images

To stack the multifocus images, we apply a guided-filter based blending
algorithm focusing on local sharpness/saliency of images [2, 4]. The local
sharpness can be computed by convolving the source image with a Laplacian of
Gaussian filter. As shown in Figure 8, the saliency values in images S1 and S2
reveal the in-focus regions in source images I1 and I2. Therefore, the in-
focus regions in each image can be selected according to the saliency values
and a linear combination of the detected in-focus regions leads to the all-in-
focus images as blending results.

The linear combination of in-focus regions requires weight maps to blend
images in particular proportions. For each pixel in the stack images, its
corresponding weight is 1 if the pixel has the highest sharpness value among
all stack images, otherwise its weight is zero. Thus, the weight map can
indicate the index of image that has an in-focus point among all pixels, and
it can be written as:

$\displaystyle
P_{n}^{k}=\left\\{\begin{matrix}1,ifS_{n}^{k}=max(S_{1}^{k},S_{2}^{k},......,S_{N}^{k})\\\
0,otherwise\end{matrix}\right.$ (11)

where ${N}$ is the number of source images and $S_{n}^{k}$ represents the
saliency value of the pixel $K$ in the $n$-th image. Figure 8 shows example
weight maps (P1 and P2) computed from a pair of images having the different
region in focus. The generated weight maps are noisy at this stage and they
can be cleaned by guided filtering.

Figure 8:  Example outputs of weight-map construction from a pair of input
images that have different regions in focus. From top to bottom: a) The source
images I1 and I2 with different in-focus regions. b) The saliency maps S1 and
S2 highlight image regions in focus. c) The weight maps P1 and P2 from the
comparison between the saliency maps. d) Guided weight maps Wb1 and Wb2 for
blending the base layer. e) Guided weight maps Wd1 and Wd2 for blending the
detail layer.

Guided filtering [2] has a strong edge-preserving property. It enhances the
edges in an input image according to another guidance image. In theory, the
guided filter is defined as linear translation-variant. The filtering output
$F$ is a linear transformation of the guidance image $I$ in a window $w_{k}$
of size $(2r+1)(2r+1)$ and centred at the pixel $k$:

$\displaystyle\ F_{i}=a_{k}I_{i}+b_{k},\forall i\in w_{k}$ (12)

where ($a_{k},b_{k}$) are linear coefficients and $i$ denotes a pixel in
window $w_{k}$. The linear coefficients are determined by minimizing the
difference between the output image $F$ and the input image $P$. Therefore, a
cost function can be written as:

$\displaystyle\ E(a_{k},b_{k})=\Sigma_{\forall i\in
w_{k}}((a_{k}I_{i}+b_{k}-P_{i})^{2}+\varepsilon a_{k}^{2}$ (13)

where $\varepsilon$ denotes a regularization parameter preventing $a_{k}$ from
being too large. By linear regression, the solution of the above equation can
be obtained as:

$\displaystyle\ a_{k}=\frac{\frac{1}{|w|}\Sigma_{i\in
w_{k}}I_{i}P_{i}-\mu_{k}\overline{P_{k}}}{\sigma_{k}^{2}+\varepsilon}$ (14)
$\displaystyle\ b_{k}=\overline{P_{k}}-a_{k}\mu_{k}$ (15)

Here $\mu_{k}$ and $\sigma_{k}$ are the mean and variance of I in $w_{k}$
correspondingly. $|w|$ denotes the number of pixels in the window $w_{k}$ and
$\overline{P_{k}}$ is the mean value of $F_{i}$ in the window. Therefore, the
filtering output can be solved. However, there are multiple windows $w_{k}$
that contains pixel i, so the value of filtering output is not constant while
computing for different windows. The strategy applied here is to average all
possible values of the output and the output of guided filtering becomes:

$\displaystyle\ F_{i}=\frac{1}{|w|}\Sigma_{k:i\in
w_{k}}(a_{k}I_{i}+b_{k})=\overline{a_{i}}I_{i}+\overline{b_{i}}$ (16)
$\displaystyle\overline{a_{i}}=\frac{1}{|w|}\Sigma_{k\in w_{k}}a_{k}$ (17)
$\displaystyle\overline{b_{i}}=\frac{1}{|w|}\Sigma_{k\in w_{k}}b_{k}$ (18)

To use guided filtering to refine the weight maps, we apply a two-scale image
decomposition as defined in [4]. Each source image is decomposed into a base
and a detail layer. The base layer is obtained by average filtering on the
source image, and it extracts large-scale variations in image intensity. The
base layer can be computed as a convolution:

$\displaystyle B_{n}=I_{n}*Z$ (19)

where $I_{n}$ denotes the $n$-th source image and $Z$ is the average filter.
The detail layer is obtained by subtracting the base layer from the source
image, and thus it keeps the small-scale texture information. The detail layer
can be computed as:

$\displaystyle D_{n}=I_{n}-B_{n}$ (20)

Then, the obtained weight maps $P_{n}$ (equation 11) are fed into the guided
filter as input, while the source images $I_{n}$ are applied as guidance.
Therefore, the guided weight maps ${W_{n}}^{B}$ and ${W_{n}}^{D}$ can be
generated for the base and dense layers:

$\displaystyle{W_{n}}^{B}=G_{r_{1},\varepsilon_{1}}(P_{n},I_{n})$ (21)
$\displaystyle{W_{n}}^{D}=G_{r_{2},\varepsilon_{2}}(P_{n},I_{n})$ (22)

where $r_{1},r_{2},\varepsilon_{1}$ and $\varepsilon_{2}$ are parameters of
the guided filter. Meanwhile, the values of the $N$ guided weight maps are
normalized such that they sum to one at each pixel $K$. These normalized
weight maps preserve clearer edges of sharp regions from the source images
comparing to the original weight maps, see Wb1, Wb2 Wd1 and Wd2 in Figure 8
for examples. Then, the base and detail layers of different source images can
be blended together according to the weight maps:

$\displaystyle\bar{B}=\sum_{n=1}^{N}W_{n}^{B}B_{n}$ (23)
$\displaystyle\bar{D}=\sum_{n=1}^{N}W_{n}^{D}D_{n}$ (24)

In the end, the fusing result $O$ can be generated by combining the blended
base layer $\bar{B}$ and the blended detail layer $\bar{D}$:

$\displaystyle O=\bar{B}+\bar{D}.$ (25)

Figure 8 shows an example weight map construction process for two source
images. However, it usually needs multiple multifocus images for recovering a
successful blending result. In that case, the guided-filter based blending
algorithm still works in the same manner, and the only difference is that it
takes more iterations to complete fusion.

###  III-D Background masking for multiview 3D reconstruction

Multiview 3D reconstruction optionally accepts masks or object silhouettes to
ignore image region belong to the background to speed up the process and to
improve the reconstruction accuracy. Automatic background segmentation can
lead to errors and often requires manual input to correct them. To improve
automatic background segmentation, [9] proposed to use two images, one in
normal front lighting and one with a strong back light. The back light is a
uniform light source such as a light box to produce a clear contrast between
background area and object area. As a result, automatic segmentation of image
with strong back light is very accurate and efficient without the need for
manual correction. The only drawback of capturing separate images with back
light is the doubling of the data storage and the time to capture and
preprocess images.

##  IV Experiments and discussion

###  IV-A Experiment summary

An overview of the multifocus multiview 3D reconstruction experiment can be
summarized as follows:

  * •

Capture a set of multifocus single-view images of a calibration dot target
facing perpendicular to the optical axis of the lens.

  * •

Calibrate camera by detecting dot and computing homography matrices for each
image relative to the reference image predefined in the set.

  * •

Capture multifocus multiview images of an object of interest without and with
back light.

  * •

Stack/blend the multifocus images using the guided filtering algorithm. Apply
separately for images without and with back light source.

  * •

Threshold the stacked image with back light to create a background mask and
add this as a transparent channel of the stacked image without back light.

  * •

Obtain the 3D model by feeding the cropped images into the open-source 3D
reconstruction software Meshroom [3, 6, 1]. The 3D reconstruction process
includes following steps:

    * –

Camera intrinsic estimation

    * –

Feature extraction using SIFT [5]

    * –

Image matching and feature matching from SIFT descriptors

    * –

Structure from motion to generate dense scene and camera poses

    * –

Depth map computation

    * –

Meshing and Texturing to produce the 3D model

###  IV-B Camera setup

Figure 5 demonstrates the lens and pan-tilt setup for the multifocus multiview
imaging experiment. The macro lens (Canon MP-E65mm f/2.8 1-5x) in the setup is
connected to a fixed upper-frame, and the camera body (Canon 5DS) is mounted
on a macro rail, which is fixed to the lower frame (right of Figure 5).
Therefore, the lens remains stationary during capturing multifocus images,
while the camera body is moving. An expandable rubber duct connects the lens
and the camera body to protect the camera sensors from ambient illumination or
dust. The target specimen is attached to a pan-tilt turntable (Stackshot 3X)
for multiview capturing. For comparison, another experiment setup with the
conventional moving lens is also tested (left of Figure 5 where the macro lens
is installed on the camera body, and thus it moves with the camera during
capturing).

In the experiment, a dry flower was chosen as a scanning target. The two
experimental setups were performed on the specimen to capture multifocus
multiview images and to reconstruct 3D models. Some parts of the flower
unfortunately broke off during scanning, so it looks slightly different
between the two image acquisitions with different lens setups.

###  IV-C Results and discussion

To reconstruct the 3D model of the target specimen, image regions belong to
backgrounds need to be removed. Therefore, we capture and extract the
silhouettes of the specimen by turning on a back-light source. As illustrated
in the previous section, the silhouette can be converted into a binary mask by
simple thresholding and by adding it to the blended image as a transparency
channel. Figure 9 shows an example of in-focus foreground images (left), in-
focus back light images (middle) and the final images with background masks
(right), for the fixed lens (top) and moving lens (bottom) setups. From the
blending results, it can be seen that the fixed lens setup produces in-focus
image with the better surface texture of the flower.

Figure 9: Examples of in-focus images from fixed lens (top) and moving lens
(bottom). Left to right: blended in-focus image, in-focus mask and in-focus
image with background mask as transparent channel. The in-focus image by fixed
lens is sharper than that by moving lens. All images are cropped to show
mostly the flower.

The blended images (with transparent backgrounds) are then fed into 3D
reconstruction software, and the dense clouds and camera poses are generated
as shown in Figure 10. Figure 10 provides a qualitative comparison of the
dense clouds and camera poses obtained from the blended images. Although there
is a difference in light intensity in input images between two cases, the
visual accuracy of the flow model and camera poses using fixed lens is
significantly better than that of moving lens. From the dense clouds, the
model generated from the moving lens is sparser, and it implies that the
conventional moving lens camera setup recovers fewer details of the surface
geometry than that of our fixed lens setup. From the estimated camera poses,
we can see that the poses reconstructed by the moving lens setup are not well-
aligned and missing (only 82 out of 100 camera positions detected). On the
contrary, the estimated camera poses generated by the fixed lens can reproduce
completed camera trajectory (all 92 camera positions detected), and the
recovered camera positions are also well-aligned as expected on a spherical
surface.

Figure 10: A comparison between structure-from-motion reconstruction from
images captured using the proposed Fixed-Lens Multifocus Capture (left) and
moving-lens multifocus capture (right). Top row and middle row: camera poses
from the reconstruction. Bottom row: sparse point cloud from reconstruction.
More uniform camera poses and denser point cloud from fixed lens indicate
better image quality than that of moving lens.

For a quantitative comparison, we use the estimated camera positions to
compute the error in radial position and pan-tilt rotation. Due to pan tilt
scanning motion, the camera positions are expected to lie on a spherical
surface. A sphere-fit is performed on the estimated camera positions, and the
radial distances to the sphere center and pan-tilt angles are obtained. Note
that only detected poses are included and carefully processed to avoid
introducing artificial outliers. As shown in Table IV-C, the mean and standard
deviation of the radial distances and angles of the poses are provided for
both the moving lens and fixed lens. The results show that the estimated
camera positions of fixed lens are much closer to the expected values, with
estimated panning motion 3.7-times less noisy, than the those of the moving
lens.

o — X[c] — X[c] — X[c] — X[c] — | Mean $\pm$ STD of radial distance | Mean $\pm$ STD of pan step | Mean $\pm$ STD of tilt step
---|---|---|---
Moving lens: pan step = 18°, tilt step = 6° | 144.270mm $\pm$ 0.266mm | 18.066°$\pm$1.060° | 5.708°$\pm$0.398°
Fixed lens: pan step = 20°, tilt step = 6° | 141.042mm $\pm$ 0.112mm |  20.017°$\pm$0.290° |  6.273°$\pm$0.313°

TABLE I: A comparison between expected and estimated camera centres

Figure 11 shows a pair of in-focus image and snapshots of reconstructed models
from a similar angle. From image 2 and 4 in Figure 11, it can be seen that the
model obtained by the fixed lens has a significantly better visual accuracy on
the reconstructed surface geometry of the target specimen. Finer details and
textures of the specimen’s surface are preserved in the model from the fixed
lens than those from moving lens. This matches with the 3D reconstruction
result shown in Figure 10, where the fixed lens setup produces a more
completed and more accurate estimation of camera positions and dense clouds.
However, deviations between the source image and reconstructed model still
exist. For example, either of the models cannot reconstruct the gaps between
leaves in the source images. The deviations might be caused by limitations of
the 3D reconstruction software, as micro-structures (such as transparent wings
of insects or the tiny gaps in the images) are usually very difficult to be
recovered. This defective reconstruction indicates the direction of future
research, which could focus on recovering the micro-structures of specimens by
testing different 3D reconstruction technique or software.

Figure 11: A qualitative comparison between reconstructed 3D models, obtained
by fixed lens and moving lens. Top (1 and 2): blended image and reconstructed
3D model obtained by fixed lens. Bottom (3 and 4): blended image and
reconstructed 3D model obtained by moving lens.

##  V Conclusion

In this paper, we proposed an image-based 3D reconstruction setup using the
proposed Fixed-Lens Multifocus Capture and the homography-based scale-shift
calibration for accurate 3D reconstruction of small-scale objects. Currently,
image-based 3D reconstruction devices using lens moving with the camera to
capture multifocus images suffer from perspective inconsistency that reduces
the accuracy of image stacking and 3D reconstruction. With the proposed Fixed-
Lens Multifocus Capture setup, the lens remains stationary while the camera
and image sensor move during multifocus image capturing. The scale-shift image
calibration via homography estimation is performed to account for the change
of relative scale and in-plane shift between images captured at different
depths. Multifocus images corrected for scale-shift change are stacked/blended
to create an in-focus image that is consistent with the perspective image
formation.

Initial multifocus multiview 3D reconstruction experiments were performed to
demonstrate the effectiveness of our proposed fixed lens setup as compared
with conventional moving lens setup. The reconstruction results show that the
proposed Fixed-Lens Multifocus Capture generated a more accurate 3D model and
more accurate camera poses (at least three times less noisy) than those by the
conventional moving lens. This can be explained by the fact that the
perspective inconsistency in in-focus images captured using moving lens
violates pinhole camera model assumption of 3D reconstruction and this leads
to reconstruction errors and quite a few camera positions are missing. By
using proposed Fixed-Lens Multifocus Capture to generate in-focus images that
preserve perspective image formation, the multiview 3D reconstruction works as
expected, detected all camera positions and created a more accurate 3D model.

## References

  * [1] AliceVision.  Meshroom - 3D reconstruction software.  https://alicevision.github.io/#meshroom, 2019.  [Online; accessed 13-July-2019].
  * [2] Kaiming He, Jian Sun, and Xiaoou Tang.  Guided image filtering.  IEEE transactions on pattern analysis and machine intelligence, 35:1397–1409, 06 2013.
  * [3] Michal Jancosek and Tomas Pajdla.  Multi-view reconstruction preserving weakly-supported surfaces.  In CVPR 2011. IEEE, jun 2011.
  * [4] Shutao Li, Xudong Kang, and Jianwen Hu.  Image fusion with guided filtering.  IEEE transactions on image processing : a publication of the IEEE Signal Processing Society, 22, 01 2013.
  * [5] David G. Lowe.  Distinctive image features from scale-invariant keypoints.  International Journal of Computer Vision, 60:91–110, 2004.
  * [6] Pierre Moulon, Pascal Monasse, and Renaud Marlet.  Adaptive structure from motion with a contrario model estimation.  In Proceedings of the Asian Computer Vision Conference (ACCV 2012), pages 257–270. Springer Berlin Heidelberg, 2012.
  * [7] Chuong V Nguyen, David R Lovell, Matt Adcock, and John La Salle.  Capturing natural-colour 3d models of insects for species discovery and diagnostics.  PloS one, 9(4):e94346, 2014.
  * [8] Steven M Seitz, Brian Curless, James Diebel, Daniel Scharstein, and Richard Szeliski.  A comparison and evaluation of multi-view stereo reconstruction algorithms.  In 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), volume 1, pages 519–528. IEEE, 2006.
  * [9] Bernhard Ströbel, Sebastian Schmelzle, Nico Blüthgen, and Michael Heethoff.  An automated device for the digitization and 3d modelling of insects, combining extended-depth-of-field and all-side multi-view imaging.  ZooKeys, (759):1, 2018.
  * [10] Wikipedia.  3dsom - 3d models from photos.  http://www.3dsom.com/, 2019.  [Online; accessed 14-July-2019].
  * [11] Wikipedia.  Agisoft - metashape (photoscan).  https://www.agisoft.com/, 2019.  [Online; accessed 14-July-2019].
  * [12] Wikipedia.  Camera obscura.  https://en.wikipedia.org/wiki/Camera_obscura, 2019.  [Online; accessed 04-July-2019].
  * [13] Wikipedia.  Circle of confusion.  https://en.wikipedia.org/wiki/Circle_of_confusion, 2019.  [Online; accessed 04-July-2019].
  * [14] Wikipedia.  Focus stacking.  https://en.wikipedia.org/wiki/Focus_stacking, 2019.  [Online; accessed 14-July-2019].