\section{Introduction}
Physical layer security for wireless communications has become a major research topic in recent years because it does not need the computational assumption \cite{InfoTheoreticSec,PhysicalLayerSecurity,Wu2018}.
Different properties of the wireless channel can be exploited using information theoretical tools to prevent leakage of information towards potential eavesdroppers. The classic wiretap model as first proposed by Wyner \cite{Wyner1975} and then generalised by I.~Csisz\'ar and J.~K\"orner \cite{Csiszar1978} was later strengthened to meet cryptographic security standards in~\cite{Csiszar1996} and \cite{hay-wire},
the latter framed within spectrum information-theoretic methods \cite{Han}.
We adopt such approach here: we assume the physical layer security realized by
a stochastic wiretap encoder \cite{Hayashi2011,H13,Bloch2015} based on the privacy amplification method \cite{BBCM,HILL}.
This method decouples reliability and secrecy, enabling the implementation of different security protocols.

In its simplest implementation, the wiretap channel model using privacy amplification can be realized as a one-way security protocol whereby Alice sends a keyless secret message to Bob protected with universal$_2$ hash functions \cite{Wegman1981,Krawczyk,Hayashi-Tsurumaru}.
In the Gaussian wiretap channel,
secrecy capacity is positive as long as
Eve's channel has a worse signal to noise ratio than the channel between Alice and Bob,
i.e., Eve's channel is degraded w.r.t the main channel between Alice and Bob \cite{Hellman1978}.
The same holds to ensure positive secrecy rate in the finite-length \cite{PhySecSat,arXiv}.
However,
when the channel between Alice and Bob has a worse signal to noise ratio than the channel to Eve,
we cannot realize secure communication in this scenario \cite{Liang2006,GLG,EU}.
For example, in the satellite communication,
Eve's satellite usually stays in lower orbit than the orbit of Bob's satellite like in the example scenario in Fig. \ref{ICC},
which implies that Eve has better signal to noise ratio than Bob.
That is, it is hard to realize secure satellite communication
with the above proposals of wiretap codes
if we cannot identify Eve's spatial locations.
This type of attack is often called passive man-in-the-middle attack \cite{Conti2016}.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.3]{TWscenario}
\protect\caption{Illustration for our satellite communication scenario.}
\label{ICC}
\end{figure}

To resolve this problem,
the papers \cite{WGH,FJHWY,ZWPT,TY} introduced two-way protocols, in which,
the channels of both directions are noisy Gaussian channels.
However, when both channels are noisy Gaussian channels,
there still exists a possibility that we cannot realize secure communication
dependently of Eve's and Bob's spatial locations.
To overcome this limitation, Maurer \cite{Maurer1993} proposed
a two-way protocol based on the binary symmetric channel.
In this paper,  inspired by Maurer's idea,
we propose a two-way protocol with Gaussian wiretap channel and
public noiseless feedback, in which,
the feedback channel is given as a public noiseless channel with discrete variable.
This paper assumes that the noise in
the channel of the initial transmission from Bob to Alice
is independent of the noise in
that to Eve while the paper \cite{Correlation} considers the case when these two noises are correlated.
Under this assumption, unless
the channel of the initial transmission to Eve
is noiseless,
this protocol always has positive secure transmission rate regardless Eve's and Bob's spatial locations.
In particular, we focus on the simple but realistic (for fixed user terminals) Gaussian-channel satellite system scenarios.
Our results demonstrate that our two-way protocol greatly outperforms
state-of-the-art one-way and two-way protocols because our two-way protocol always realize secure communication under passive man-in-the-middle attack independently of Eve's spatial location.
Note that extension of our results to channels with fading is straightforward, but we take this purely averaging calculation problem out of the scope of our paper, here we focus on the protocols for the Gaussian channel
without fading. More specifically, on the application to the Gaussian with BPSK satellite channel, which is considered in current satellite communication standards \cite{DVB-S2X}
Further, we can equip authentication in our protocol
by attaching universal$_2$ hash function \cite{Wegman1981,Krawczyk}.
This protocol can prevent active man-in-the-middle attack.
Therefore, the advantages of our protocol are summarized as follows.

Contribution 1) We address wiretap channel security problem (i.e. eavesdropping) and propose for the first time a novel practical Gaussian wiretap protocol implementing theoretical ideas in Maurer's paper \cite{Maurer1993}.

Contribution 2) Our novel two-way protocol greatly outperforms state-of-the-art one-way protocol because our protocol shows always positive secrecy capacity independently of Eve's location (i.e. it does not require channel degradation condition of one-way wiretap channel). This is our main technical result based on novel and rigorous information theoretical proof.

Contribution 3) Our practical two-way protocol outperforms other proposals of two-way protocols \cite{WGH,FJHWY,ZWPT} in the following sense.
First, because while other two-way protocols require several communication rounds, our protocol only requires two rounds. Second, because other two-way protocols while outperforming one-way protocol, they still may have negative secrecy capacity. Finally, because our protocol only requires two rounds, our protocol is highly
suitable to secure communication channels with large delay, e.g. satellite channels.

Contribution 4) We show the performance of our protocol with meaningful numerical results for the realistic Gaussian BPSK modulated satellite channel, which is included in current satellite communication standards \cite{DVB-S2X}. For this, we use our system modeling which allows to evaluate the security capacity as a function of the system parameters. Hence, our method and results are useful for secure communication design.

\if0
Contribution 5) Our graphical numerical results allow to draw consistent (rather than purely theoretical) conclusions about the advantages of our protocol compared to other currently advanced solutions such as quantum key distribution protocols. We also note that, as well as in any other competing solution, our protocol needs authentication, which we show is possible using universal$_2$ hashing \cite{Wegman1981,Krawczyk,Fung,Correlation}.
Moreover, our protocol is clearly compatible with complementary security protocols at upper layers.
\fi

One might consider that the real feedback channel is also a noisy channel.
However, if we choose sufficiently strong intensity and a
suitable error correcting code, the information transmission of
the feedback channel can be regarded as a noiseless channel.
In this case, we can regard the feedback channel as a noiseless
public channel with discrete variable.
In contrast, the noise of
the initial Gaussian channel is essential because its presence
makes the difference between mutual informations from Alice
to Bob and Eve.
Furthermore, the information leakage in the channel during the second phase needs not be considered because information leakage is only relevant on the first phase.
Hence, it is allowed to make the power of the transmission signal very strong
in the second phase
while we cannot use such a strong power in the first phase
to control the secrecy.
Since the information transmission rate with the strong power
is much larger than that with the weak power,
the consuming time of the first phase is dominant
in comparison with that of the second phase is dominant.
That is, the first phase is the bottleneck in this setting.
Therefore, this paper optimizes the amount of
noise in the initial Gaussian channel to maximize the wiretap
capacity in this model.

In fact, the paper \cite{RealisticChannel} considers a similar topic.
It is a follow up of this submission with practical focus.
Hence, it contains only the brief description of the two-way protocol without proper proof.
Also, the analysis of the secure satellite communication in \cite{RealisticChannel} is different from our analysis in Sections IV and V.
That is, it did not consider the optimization while the analysis in this paper is based on the optimization given in Section IV.

Relations to other studies are summarized as follows.
While the paper \cite{TY} discusses two-way wiretap channels,
it considers a new scheme cooperative jamming.
The scheme in \cite{TY} has several users that are cooperative
while this paper has only two cooperative users, the legitimate sender and the legitimate receiver.
Therefore, the method in this paper cannot be compared with \cite{TY}.
In the paper \cite{AFJK}, Bob feeds back some randomness that is used as a secret key,
exactly like the present manuscript.
However, it assumes that the feedback is noiseless and
secure, so Eve does not observe it, which is different from here.
In the paper \cite{BPS},
unlike the present work, Bob does not control the feedback link.
However, it allows all players to observe everything.
Hence, the method in this paper cannot be compared with
the results \cite{AFJK,BPS}.

Our work is structured as follows. In Section II,
we review the results of one-way standard protocol.
In Section III, we propose our two-way protocol.
In Section IV, we make numerical optimization for our obtained secret capacities.
In Section V, we  apply our result to simple realistic models of satellite communication.
Finally in Section VI we discuss the protocol and draw some conclusions.
\section{One-Way Physical Layer Security Protocols}
\subsection{Signal and channel Model}
First, we review the results of one-way standard protocol with Gaussian channels, which model wireless channels in relevant realistic communication scenarios such as satellite radiofrequency communication channels.
The signal variables received by Bob and Eve
can be modeled as
\begin{equation}\label{eq:SysModel-P}
Y'=\sqrt{E_s}V +\sqrt{n_{B}} N_1, \quad
Z'=\gamma_g\sqrt{E_s}V +\sqrt{\gamma_n n_{B}} N_2,
\end{equation}
where $V$ is a variable modeling the transmitted signal and $E_s$ is the energy per symbol expressed in Joules.
$Y'$ is the random variable representing the signal received at the legitimate receiver, $Z'$ the random variable representing the signal received at the eavesdropper's receiver and $N_1$ and $N_2$ are zero-mean circular complex Gaussian random variables with unit variance. $n_B$ is the noise spectral density power of Bob's receiver expressed in Joules per Hertz. The coefficient $\gamma_g$ models the amplitude attenuation of the wiretapper's channel w.r.t. the legitimate channel. The analytical expression of $\gamma_g$ will depend on the system under analysis as well as the channel assumptions and corresponding time scale. The multiplicative coefficient $\gamma_n$ expresses wiretapper's receiver noise with respect to Bob's receiver noise. Denote the signal-to-noise ratio (SNR) for Bob as
$$ \eta_B = \frac{E_s}{n_{B}}=\frac{P}{N_{B}}$$
where $P$ and $N_B$ are the system and noise power at Bob's receiver, respectively both expressed in Watts. We can then rewrite the signal model as
\begin{equation}\label{eq:SysModel}
Y=\sqrt{\eta_B}V+N_1, \quad
Z=\gamma_B\sqrt{\eta_B}V+N_2;
\end{equation}
where $Y = Y'/\sqrt{n_{B}}$, $Z = Z'/\sqrt{\gamma_n n_{B}}$ and $\gamma_B := \gamma_g/\sqrt{\gamma_n}$.
Hence, $N_1$ and $N_2$ are
zero-mean circular complex Gaussian random variables with unit variance.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.3]{OneWayProtocol}
\protect\caption{Graphical illustration of the one way protocol with Gaussian channels.
$N_B$ and $N_E$ are noise powers at Bob's and Eve's receivers.
}
\label{fig:OneWayProtocol}
\end{figure}

\subsection{Secrecy capacity with BPSK modulation and soft decision}
In the one-way model, Alice sends the encoded information to Bob via channel  \eqref{eq:SysModel} as
Fig. \ref{fig:OneWayProtocol}.
Now, we assume the BPSK modulation, in which,
Alice encodes her binary information $A \in {\mathbb{F}}_2$ to $V=(-1)^A$.
In this scenario,
the secrecy capacity $C^{\mathop{\rm OW}}_{\mathop{\rm soft}}(\gamma_B,\eta_B)$
for the one-way protocol is given as \cite[(46)]{arXiv}
\begin{align*}
C^{\mathop{\rm OW}}_{\mathop{\rm soft}}(\gamma_B,\eta_B
=&I (A; Y ) - I (A; Z) \nonumber \\
=&\int_{-\infty}^\infty
u \Bigg[ \frac{1}{\sqrt{8 \pi}}\Big(e^{-\frac{(y-\sqrt{\eta_B} )^2}{2}}
+e^{-\frac{(y+\sqrt{\eta_B} )^2}{2}}\Big)\Bigg]
d y \nonumber \\
&-\int_{-\infty}^\infty
u \Bigg[ \frac{1}{\sqrt{8 \pi}}\Big(
e^{-\frac{(z-\gamma_B\sqrt{\eta_B} )^2}{2}}
\!+\!e^{-\frac{(z+\gamma_B\sqrt{\eta_B} )^2}{2}}
\Big) \Bigg]
d z ,
\end{align*}
where $u(x):= - x \log x$
when
\begin{align}
\gamma_B<1 , \quad {\rm i.e.,} \quad \gamma_g^2 <\gamma_n.
\label{LPK}
\end{align}
Also, when the condition (\ref{LPK}) does not hold,
the capacity is zero. In this case, we cannot realize secure communication in this scheme.
Here,
Bob and Eve are assumed to store the sequence of the received continuous signals
and apply the decoder to them.
This type of information processing is called soft decision decoding \cite{Proakis2001}.

\subsection{Secrecy capacity with BPSK modulation and hard decision}
To save the cost of decoding,
the receiver converts the received continuous signal to a binary signal in the reception
and apply the decoder to the sequence of the binary signals.
This type of information processing is called hard decision decoding \cite{Proakis2001}.
When the receiver applies this method,
it is sufficient to store only binary signals, which saves the memory of the receiver.
Here, as another scenario, we consider the case when Bob and Eve obtain $B$ and $E$ using hard decision detection on their received signals $Y$ and $Z$ as defined in the previous sections, respectively.
The crossover probability between Alice and Bob induced by the Bernoulli random variable $X_1$ is given as $\epsilon_B^*= 0.5 \mbox{erfc}(\sqrt{\eta_B/2})$ and the crossover probability between Alice and Eve induced by the Bernoulli random variable $X_2$ is given as $\epsilon_E^*= 0.5 \mbox{erfc}(\gamma_B \sqrt{\eta_B/2})$ with
\begin{align}
\mbox{erfc}(t):=& \frac{2}{\sqrt{\pi}} \int_{t}^{\infty}e^{-t^2}dt  .\label{LBH}
\end{align}
In this scenario,
by using the binary entropy function $h(x):= -x \log x -(1-x)\log (1-x)$,
the secrecy capacity
for the one-way protocol is given as
\begin{align}
C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_B,\eta_B)
= h(\epsilon_E^*)-h( \epsilon_B^*),
\end{align}
which is positive only when $\epsilon_E^* > \epsilon_B^*$, which is equivalent to \eqref{LPK}.

\section{Two-Way Physical Layer Security Protocol with Gaussian channels and BPSK modulation}

\subsection{Signal and channel Model}
One-way model requires the condition that
the mutual information between the sender and the legitimate receiver
is larger than that between the sender and the eavesdropper.
This assumption does not hold when
the eavesdropper performs passive man-in-the-middle attack.
To resolve this problem, we consider two-way protocol for the Gaussian channel and BPSK modulation
as follows.
In an initial step, Bob sends a random variable $V$ to Alice.
In this case, Alice and Eve obtain the variables $Y$ and $Z$, respectively, as follows.
\begin{align}
Y=\sqrt{\eta_A}V+N_1, \quad
Z=\gamma_A\sqrt{\eta_A}V
+N_2,
\label{eq:GaussChannel}
\end{align}
where $N_1$ and $N_2$ are zero-mean circular complex Gaussian random variables with unit variance and the coefficient $\gamma_A$ models the amplitude attenuation of the wiretapper's channel w.r.t. the legitimate channel (which is now between Bob and Alice and not between Alice and Bob).
Notice that the transmitter of the noisy Gaussian channel \eqref{eq:GaussChannel} is not Alice who is
the sender of the secret message of this protocol.
Therefore $\eta_A$ is now
$$ \eta_A = \frac{E_s}{n_{A}}=\frac{P}{N_{A}},$$
where $n_A$ is the noise spectral density power of Alice's receiver expressed in Joules per Hertz.

Bob generates the binary variable $B \in {\mathbb{F}}_2$ subject to the binary uniform distribution, and sends $V = (-1)^{B}$ via the above RF channel. Applying hard decision decoding to $Y$, Alice obtains the binary variable $A$.
In the next step, Alice prepares another binary variable $X$, and sends $X':= X\oplus A$ to Bob via a public channel. When $X$ is regarded as the channel input information, the legitimate receiver's output is $B$ and $X'$ while the eavesdropper's output is $Z$ and $X'$. The overall process along with the generated random variables is shown in Fig. \ref{fig:TwoWay_Gaussian}.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.4]{TwoWay_Gaussian}
\protect\caption{Graphical illustration of the two-way protocol with Gaussian channels.
Phase 1 is shown in red and Phase 2 is shown in black.}
\label{fig:TwoWay_Gaussian}
\end{figure}

\subsection{Protocol}
Based on the above discussion, we fully describe our concrete protocol.
For this aim, we fix an error correction code, i.e., the pair of the encoder $\phi_{e,n}$ and
the decoder $\phi_{d,n}$ with block length $n$.
Then,
combining the error correction code and universal$_2$ hash functions \cite{Wegman1981,Krawczyk,Hayashi-Tsurumaru},
we employ the wiretap code given in
\cite[Appendix A]{arXiv}
using random seed $S$ and we have the wiretap encoder $\phi_{e,n|S}$ and the wiretap decoder $\phi_{d,n|S}.$
This code construction achieves the strong secrecy even in the continuous system \cite{Hayashi2011,PhySecSat,arXiv}.

Then, we propose the following two way protocol.
\begin{description}
\item[(1)]
Bob generates the binary data sequence $B_1, \ldots, B_n \in {\mathbb{F}}_2$.
Then, sends $(-1)^{B_1}, \ldots, (-1)^{B_n}$ via the channel described in (\ref{eq:GaussChannel}).

\item[(2)]
Alice makes the hard decision.
Then, she obtains the binary data $A_1, \ldots, A_n$.
In this case, $A_i$ and $B_i$ are connected via the binary symmetric channel with crossover probability
$\epsilon_A$, whose mathematical expression will be given later.

\item[(3)]
To send the secret message $M$, using an auxiliary variable $L$, Alice applies the code
$X^n=(X_1, \ldots, X_n):= \phi_{e,n|S}(M,L) \in {\mathbb{F}}_2^n$.
Then, Alice sends ${X'}^n=(X_1', \ldots, X_n')$ to Bob via public channel, where
$X_i':= A_i\oplus X_i$.
Here, an auxiliary variable $L$ is a variable independent of the message $M$,
and is used to realize the secrecy of $M$.

\item[(4)]
Bob decodes $M$ by $\phi_{d,n|S}({X^n}'')$, where
${X^n}''=(X_1'', \ldots, X_n'')$ and $X_i'':=X_i'\oplus B_i$.
\end{description}

To realize the public channel from Alice to Bob in the second phase,
they employ the RF channel \eqref{eq:SysModel} with
an error correcting code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$
different from the error correction $\phi_{e,n|S}$
so that the decoding error probability of the code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$ is close to zero (i.e., the bit error rate is e.g. below $10^{-6}$).
Here, if the coding rate of the code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$ is $\hat{R}$,
they use the RF channel \eqref{eq:SysModel} $\hat{n}=n/\hat{R}$ times in Step (3) physically.
That is, in Step (3), Alice sends $\hat{X}^{\hat{n}}=\hat{\phi}_{e,\hat{n}} ({X'}^n) $
to Bob via the RF channel \eqref{eq:SysModel} of coefficient $\eta_B$.
Also, in Step (4), to get ${X'}^n$, Bob applies the decoder $\hat{\phi}_{d,\hat{n}}$
to the received $\hat{n}$ symbols via the RF channel \eqref{eq:SysModel} of coefficient $\eta_B$.
Indeed, when the bit error rate of the public channel (i.e., the bit error rate
of the code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$)
is e.g. below $10^{-6}$,
it can be negligible in comparison with the bit error rate between $A$ and $B$.
Hence, we can consider that
the bit error rate of the channel from $X_i$ to $X_i''$
given in Steps (3) and (4)
almost equals that between $A$ and $B$, in practice.

The above description has no authentication.
However, it is possible by attaching using universal$_2$ hash function \cite{Wegman1981,Krawczyk},
which prevents active man-in-the-middle attack \cite{Conti2016},
while it is difficult to avoid active man-in-the-middle attacks without authentication \cite{Eberz}.
The detail is discussed in \cite{Fung} and the arXiv version of \cite{Correlation}.

\subsection{Secrecy capacity when Eve uses hard decision}
When Eve has limited memory, it is natural that Eve uses hard detection decoding
when receiving $Z$
so that Eve obtains the binary variable $E$.

Indeed, the first phase can be regarded as a preparation step for the secure communication.
In order to prevent Eve to make soft decision, Alice and Bob can consider the following strategy.
Before starting the second phase, Alice and Bob continue the first phase so that the length of their obtained random numbers
$A_1, \ldots, A_n$ and $B_1, \ldots, B_n$ is close to the limitation of their memory.
In this case, the length of their obtained random numbers is across several coding blocks.
Since satellite has a limitation of size of memory due to the limitation of physical space.
it is natural that the size of Eve's memory is similar to that of Alice and Bob.
In this case, it is difficult for Eve to keep the all the outcomes of soft decision, i.e.,
Eve needs to choose hard decision in this case.
For example, the preceding paper \cite{Endo2015}, which is oriented to an application side,
analyzed the security for the Poisson wiretap channel when Eve has limited memory, i.e.,
Eve uses hard detection decoding.

Using two independent Bernoulli random variables $X_1$ and $X_2$ on ${\mathbb{F}}_2$,
we have
\begin{align}
\label{eq:BSC1}
A=B\oplus X_1, \quad
E=B\oplus X_2.
\end{align}
The crossover probability between $A$ and $B$ is $\epsilon_A= 0.5 \mbox{erfc}(\sqrt{\eta_A/2})$
and the crossover probability between $E$ and $B$ is $\epsilon_E= 0.5 \mbox{erfc}(\gamma_A \sqrt{\eta_A/2})$,
where $\mbox{erfc}(t)$ is defined in \eqref{LBH}.
Hence,
the problem is reduced to the case with BSC channels, which was discussed by Maurer \cite{Maurer1993}.
Hence,
the capacity $C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$
when Eve uses hard detection is calculated to
\begin{align}
&C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_A,\eta_A) =
I(A;B|E)=
I(A;B)-I(A;E)\nonumber \\
=&
H(B)-H(B|A)-(H(E)-H(E|A)) \nonumber \\
=&H(E|A)-H(B|A)
= h(\epsilon_E +\epsilon_A -2\epsilon_E \epsilon_A)-h(\epsilon_A)
\end{align}
because $H(B)=H(E)=h(\frac{1}{2})$ and
the probability of $E\neq A$ is
$\epsilon_E(1-\epsilon_A) +(1-\epsilon_E)\epsilon_A
=\epsilon_E +\epsilon_A -2\epsilon_E \epsilon_A $.
When $ \gamma_A=\gamma_B$ and $\eta_A=\eta_B$, i.e.,
$\epsilon_B^*=\epsilon_A$
and
$\epsilon_E^*=\epsilon_E$, we have
\begin{align}
C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)
\ge C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_B,\eta_B)
\end{align}
because $
h(\epsilon_E +\epsilon_A -2\epsilon_E \epsilon_A)
=h(\epsilon_E +\epsilon_A (1-2\epsilon_E))
\ge h(\epsilon_E)=h(\epsilon_E^*)$.

\subsection{Secrecy capacity when Eve uses soft decision}
When Eve has sufficient size of memory
and her computation power is unlimited,
she can employ soft decision decoding.
That is, to consider Eve's best strategy, we need to address the case when Eve
uses the variables $Z$ and $X'$.
To analyze this case, we use the Markov chain $A-B-Z$.
We focus on the wiretap channel composed of the main channel $W_B:= P_{B,X'|X} $ and the eavesdropper channel $W_E:= P_{Z,X'|X} $.
Since $A=X'\oplus X$, we have $P_{B,X'|X} (b,x'|x)= P_{B|A}(b|x'\oplus x)P_A(x'\oplus x)$.
Hence,
the Markov chain $A-B-Z$ condition guarantees that
\begin{align}
& P_{Z,X'|X}(z,x'|x)\nonumber \\
=&
\sum_{b}P_{Z|B}(z|b) P_{B|A}(b|x'\oplus x)P_A(x'\oplus x) \nonumber\\
=&
\sum_{b}P_{Z|B}(z|b) P_{B,X'|X}(b,x'|x).
\end{align}
Thus, the channel $W_E$ is a degraded channel of the channel $W_B$.
Further,
the channels $W_B$ and $W_E$
are symmetric, the wiretap capacity is attained when $P_X$ is the binary uniform distribution,
and the wiretap capacity $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$ is calculated to
\begin{align}
& C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)=I(X;X'B)-I(X;X'Z) \nonumber \\
=& I(X;X')+I(X;B|X')-(I(X;X')+I(X;Z|X'))\nonumber \\
=& I(A\oplus X';B|X')-I(A\oplus X';Z|X')\nonumber \\
=& I(A;B|X')-I(A;Z|X')
\stackrel{(a)}{=} I(A;B)-I(A;Z)\nonumber \\
=& I(A;BZ)-I(A;Z)= I(A;B|Z),
\end{align}
where $(a)$ follows from the independence of $X'$ from $A,B,Z$,
which can be shown by the uniformity of the conditional distribution $P_{X'|A}$.
Therefore, the wiretap capacity $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}$ is always positive regardless of $\gamma_A$,
regardless the condition $\gamma_A <1$ does not hold.
The wiretap capacity expresses the limit of the secure transmission rate
when we use a proper coding under the condition that
the mutual information between the message and Eve's information goes to zero.

Further,
the probability $P_Z$ and
the conditional probability $P_{B|Z}$ are calculated as
\begin{align}
P_{Z}(z)=&
\sum_{b=0}^1 P_{Z,B}(z,b)=
\sum_{b=0}^1 P_{Z|B}(z|b)P_B(b) \nonumber\\
=&
\frac{1}{2\sqrt{2 \pi }}
(e^{-\frac{(z+\gamma_A\sqrt{\eta_A})^2}{2}}
+e^{-\frac{(z-\gamma_A\sqrt{\eta_A})^2}{2}}),\\
P_{B|Z}(b|z)=&
\frac{P_{B,Z}(b,z)}{P_Z(z)}
= \frac{P_{Z|B}(z|b) P_B(b)}{ P_Z(z)}
= \frac{P_{Z|B}(z|b)}{2 P_Z(z)}\nonumber \\
= &
\frac{e^{-\frac{(z-(-1)^b\gamma_A \sqrt{\eta_A})^2}{2}}
}
{
e^{-\frac{(z+\gamma_A\sqrt{\eta_A})^2}{2}}
+e^{-\frac{(z-\gamma_A\sqrt{\eta_A})^2}{2}}
}.
\end{align}
Hence,
due to the Markov chain $Z-B-A$,
we can calculate the conditional mutual information;
\begin{align*}
 I(A;B|Z=z)
 =
 h\Big(P_{B|Z}(0|z) \epsilon_A+ P_{B|Z}(1|z)(1-\epsilon_A)\Big)-h(\epsilon_A).
\end{align*}
Notice that $ I(A;B|Z=z) =0$ if and only if
$P_{B|Z}(0|z)$ is $0$ or $1$.
Hence,
the capacity $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$ with Eve's soft decision
is calculated as a function of $\eta_A,\gamma$ by
\begin{align}
&C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)= I(A;B|Z) \nonumber \\
 =&
\int_{-\infty}^{\infty}
P_Z(z)dz
\Big( h\Big(  P_{B|Z}(0|z) \epsilon_A+ P_{B|Z}(1|z)(1-\epsilon_A)\Big)
-h(\epsilon_A)
\Big).\label{LLF}
\end{align}
Thus, unless $P_{B|Z}(0|z)$ is $0$ or $1$ for all $z$,
\eqref{LLF} is strictly positive.
That is, when $ \gamma_A$ is a finite value,
$P_{B|Z}(0|z)$ is an intermediate value between $0$ and $1$ for all $z$.
Hence, the capacity $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$ is strictly positive.

Even in this scenario,
when $X$ is also subject to the uniform distribution,
$B\oplus X'$ is Bob's sufficient statistics with respect to $X$.
Hence, we have
$I(X;X'\oplus B)-I(X;X'Z)
=I(X;X'B)-I(X;X'Z)
=C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$.
Therefore, even when Bob uses only $B\oplus X' $ for his decoding
while Eve uses the two variables $Z$ and $X'$,
the capacity $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$ can be attained.
Consequently, the channel $W_B$ can be modeled as a computation channel as illustrated in Fig. \ref{fig:ComputationChannel}.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.3]{ComputationChannel}
\protect\caption{Bob's computation channel model for the two-way protocol. Independently of whether Eve uses hard or soft detection, two-way secrecy capacity can be attained with this computation model.}
\label{fig:ComputationChannel}
\end{figure}

Now, we compare
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$
and $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_B,\eta_B)$
when $ \gamma_A=\gamma_B=\gamma$ and $\eta_A=\eta_B=\eta$.
For this comparison, we fix $\eta$ and change $\gamma$.
Due to the above discussion, when $\gamma>1$,
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma,\eta)$ is larger than $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma,\eta)$.
In contrast, under the limit $\gamma\to 0$, we have
\begin{align}
& \lim_{\gamma\to 0}C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma,\eta)
=\! \log 2 - h(\epsilon_A)
=\! I(B_{\mathop{\rm TW}};A_{\mathop{\rm TW}})
<  I (B_{\mathop{\rm TW}};Y_{\mathop{\rm TW}})
\nonumber\\
= &I (A_{\mathop{\rm OW}};Y_{\mathop{\rm OW}}) =I (X_{\mathop{\rm OW}};Y_{\mathop{\rm OW}})
=\lim_{\gamma\to 0}C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma,\eta),
\label{CapPeor}
\end{align}
where the subscripts $\mathop{\rm TW}$ and $\mathop{\rm OW}$ of  the random variables
express the protocol to be considered.
This opposite inequality is caused by the hard decision on Alice's received signal $Y$ in the two way protocol.
Therefore, when $\gamma$ is smaller than a certain threshold, $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma,\eta)$ is smaller than $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma,\eta)$.
That is, the one-way protocol may have greater capacity than the two-way protocol for some threshold of $\gamma$.
We have computed the value of such threshold as a function of the SNR, which is shown in Fig. \ref{fig:umbralesX}.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.17]{gamma_threshold}
\protect\caption{Threshold values of $\gamma$ at which $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma,\eta)$ equals $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma,\eta)$.}
\label{fig:umbralesX}
\end{figure}

\if0

Since we do not care the information leakage of $X_1', \ldots, X_n'$ to Eve,
we do not need to care information leakage in the RF channel \eqref{eq:SysModel}.
Hence,
even when the ratio $\gamma_B$ is fixed,
we can choose so large $\eta_B$ that
$\hat{R}$ is large and the decoding error probability of the code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$ is close to zero (i.e., the bit error rate is $10^{-6}$).
Further, for the code $(\hat{\phi}_{e,\hat{n}},\hat{\phi}_{d,\hat{n}})$, we do not need to
restrict the modulation scheme to the BPSK.
We can choose any modulation to optimize the transmission rate $\hat{R}$.

Hence, it is natural to assume that we can make the transmission rate $\hat{R}$ larger than $1$
by choosing a suitable code and a suitable modulation scheme (e.g., 8PSK etc) with sufficiently large
$\eta_B$.
In this case, we have $\hat{n} \le n$.
If the cost of one signal transmission via the RF channel \eqref{eq:SysModel} from Alice to Bob
equals
that via the RF channel \eqref{eq:GaussChannel} from Bob to Alice,
we do not need to decrease the number $\hat{n}$ and need to increase
the rate of the code $(\phi_{e,n|S},\phi_{d,n|S})$, which can be considered as the bottleneck.
Hence, it is sufficient to optimize
$C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$ and
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_A,\eta_A)$.

\fi

\section{Optimization}
Here, to extract a higher communication speed,
we consider how to optimize the channel parameters
in the RF channel \eqref{eq:GaussChannel}.
When the power of transmitting antenna of Bob increases,
the coefficient ${\eta_A}$ increases and
the ratio between the coefficients of signal in Alice's and Eve's sides is not changed.
For simple analysis, we first assume that
Alice and Bob can know the value of $\eta_A$ by using test transmission,
and control it by changing the power of transmitting antenna of Bob,
where other components (e.g., the receiving antenna gains of Alice and Eve, the directions of antennas etc)
are fixed.
In practice, it is not so easy to know the value of the ratio $\gamma_A$
because it depends on Eve's position.
However, if we know the type of Eve's orbit, we know the range ${\cal G}$ of possible values of $\gamma_A$.
In this case, we consider the worst case for Alice and Bob, i.e.,
$\gamma_{A,\max}:=\max_{\gamma_A \in {\cal G}} \gamma_A$.
In fact, when we have two possible values $\gamma_{A,1}>\gamma_{A,2}$ for the ratio,
the channel to Eve with $\gamma_{A,2}$ is a degraded channel of the channel to Eve with $\gamma_{A,1}$.
Hence \footnote{Note that in a traditional communication scenario, the transmission power is designed so that the link budget can provide a required link quality (e.g. in terms of target bit error rate) is met. However, here we are designing a secure communication scenario and therefore the link budget is constrained to meet the security requirements. Such security requirement means here that the link budget is designed to provide the maximum secrecy capacity.}, a secure code for the channel to Eve with $\gamma_{A,1}$
is also secure for the channel to Eve with $\gamma_{A,2}$.
Therefore, it is sufficient to prepare a code with the largest value
$\gamma_{A,\max}$.
We optimize
$C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A)$ and
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A)$ by changing $\eta_A$.
Here, the parameter $\eta_A=P/N_A$ can be changed by changing the power $P$.
The optimum secret capacities are given as
\begin{align}
C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})&:=\max_{\eta_A}C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A) \\
C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max})&:=\max_{\eta_A}C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A).
\end{align}
Hence, we need to find suitable value for $\eta_A$ dependently of
$\gamma_{A,\max}$.
$\eta_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max}):= \mathbb{\rm argmax}_{\eta_A}C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A)$
and
$\eta_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max}):= \mathbb{\rm argmax}_{\eta_A}C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max},\eta_A)$
are the optimal intensities of $\eta_A$.

In fact, we can apply a similar optimization to the one way case.
In this case, we consider the following optimum secret capacities;
\begin{align}
C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})& :=\max_{\eta_B}C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max},\eta_B)\\
C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max})&:=\max_{\eta_B}C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max},\eta_B),
\end{align}
where $\gamma_{B,\max}$ is the maximum value among possible values of
$\gamma_{B}$.
$\eta_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max}):= \mathbb{\rm argmax}_{\eta_B}C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max},\eta_B)$
and
$\eta_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max}):= \mathbb{\rm argmax}_{\eta_B}C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max},\eta_B)$
are the optimal intensities of $\eta_B$.

Now we show numerical calculations to compare the capacities of the one-way protocol and the two-way protocol. For easy visualization, we calculate numerical values assuming $\eta_B=\eta_A$ and $\gamma_A=\gamma_B$.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.2]{OptimalCaps}
\caption{Comparison among the optimum secret capacities
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})$, and
$C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max})$.
The vertical axis expresses these optimal capacities.
The horizontal axis expresses
$\gamma_{A,\max}$ and $\gamma_{B,\max}$ with log scale, which runs from $10^{-3}$ to $10^2$.}
\label{fig:1_Comparison}
\end{figure}

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.17]{gamma_threshold1}
\protect\caption{Comparison of OW and TW optimal secrecy capacities assuming $\gamma_{B,\max}=\gamma_{A,\max}=\gamma_{\max}$.}
\label{fig:umbrales}
\end{figure}

Fig. \ref{fig:1_Comparison} shows the comparison
among the optimum secret capacities
$C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$C_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})$, and
$C_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max})$.
We can observe that while for the one-way protocol the capacity is zero whenever the eavesdropper has higher SNR than Bob, in the two-way protocol the capacity is always positive and greater than zero. We have computed the difference between optimal secrecy capacities $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})$ and $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})$ for $\gamma_{B,\max}=\gamma_{A,\max}=\gamma_{\max}$
as Fig. \ref{fig:umbrales}.
We observe that as obtained theoretically in \eqref{CapPeor}, the OW secrecy capacity is slightly bigger than the TW secrecy capacity when the channel to Bob is advantageous over that to Eve.
However, also in agreement with the theoretical derivations, the TW secrecy capacity starts to be bigger than
the OW
when the channel to Bob is not so advantageous over that to Eve.
In  Fig. \ref{fig:umbrales}, we observe in the zoom plot that this occurs at $\gamma_{max}=0.3185$.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.22]{OptimalEtas}
\protect\caption{
Optimal intensities
$\eta_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$\eta_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$\eta_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})$, and
$\eta_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max})$.
The vertical axis expresses these optimal intensities with log scale.
The horizontal axis expresses
$\gamma_{A,\max}$ and $\gamma_{B,\max}$ with log scale,
where
$\gamma_{A,\max}$ runs from $10^{-3}$ to $10^2$, but
$\gamma_{B,\max}$ runs from $10^{-3}$ to $1$.
}
\label{fig:2_Comparison}
\end{figure}

Fig. \ref{fig:2_Comparison} shows
the optimal intensities
$\eta_{\mathop{\rm soft}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$\eta_{\mathop{\rm hard}}^{\mathop{\rm TW}}(\gamma_{A,\max})$,
$\eta_{\mathop{\rm soft}}^{\mathop{\rm OW}}(\gamma_{B,\max})$, and
$\eta_{\mathop{\rm hard}}^{\mathop{\rm OW}}(\gamma_{B,\max})$.
These values are the optimal choices for the intensity $\eta_A$ or $\eta_B$ in the respective cases.

\section{Application to real satellite communication}
Now, we apply our analysis to the following two types of real satellite communication scenarios.
\begin{description}
\item[(I)]
The transmitter is the earth station and
the legitimate receiver is the GEO satellite in the noisy Gaussian channels
\eqref{eq:SysModel} and \eqref{eq:GaussChannel}.
That is, Alice is the earth station and Bob is the GEO satellite in the OW,
and
Bob is the earth station and Alice is the GEO satellite in the TW.
Notice that
the noiseless public channel from the GEO satellite to the earth station
is also required by using a proper combination of wireless communication and
outer error correcting code in the TW.
The noisy Gaussian channels \eqref{eq:SysModel} and \eqref{eq:GaussChannel}
of these scenarios
are explained in the two figures on the top in Fig \ref{fig:geometry}, which describe the case when the Earth station is the information data communication source in the noisy Gaussian channel.
In this case, the eavesdropper, Eve is a low Earth orbit (LEO) satellite
or a medium Earth orbit (MEO) satellite.

\item[(II)]
The transmitter is the GEO satellite and
the legitimate receiver is the earth station in the noisy Gaussian channels \eqref{eq:SysModel} and \eqref{eq:GaussChannel}.
That is, Alice is the GEO satellite and Bob is the earth station in the OW,
and
Bob is the GEO satellite and Alice is the earth station in the TW.
The noisy Gaussian channels \eqref{eq:SysModel} and \eqref{eq:GaussChannel} of these scenarios
are explained in the two figures on the down in Fig \ref{fig:geometry},
which describe the case when the GEO satellite is the information data communication source in the noisy Gaussian channel.
In this case, Eve is
an LEO satellite, an MEO satellite, or a GEO satellite.
\end{description}

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.4]{geo1}
\includegraphics[scale=0.4]{geo2}

\protect\caption{Geometry considered to illustrate the real satellite communication secrecy analysis.
The two figures on the top describe Case (I), in which the Earth station is the source
in the noisy Gaussian channels \eqref{eq:SysModel} and \eqref{eq:GaussChannel}.
The two figures on the bottom describe Case (II), in which the satellite is the source
in the noisy Gaussian channels \eqref{eq:SysModel} and \eqref{eq:GaussChannel}.
The geometries for OW (left) and TW (right) cases illustrate LEO and MEO orbit heights.
The dashed line expresses the noiseless public channel from Alice to Bob in the TW case.
The eavesdropper can be at any LEO or MEO orbit while the legitimate transmitter and receiver are either the earth station of the GEO satellite. }
\label{fig:geometry}
\end{figure}

Let $\alpha(\theta)$ be the normalised transmitter's antenna radiation's pattern of the earth station
in response to the angle $\theta$ from the boresight axis directed to
the GEO satellite to account for spatial attenuation.
The function $\alpha(\theta)$ can be considered exactly in case the
(normalized) antenna pattern is known,
or otherwise it can be considered in terms of the allowed emission of radiation according to space regulations. A typical analytical expression for $\alpha(\theta)$ is
\begin{equation}\label{H31}
\alpha(\theta) := \frac{J_1(k\sin(\theta))}{2k\sin(\theta)}
+ 36\frac{J_3(k\sin(\theta))}{(k\sin(\theta))^3}
\end{equation}
where $k = 2.0712/\sin(\theta^{3dB})$, with $\theta^{3dB}$
being the one-sided half-power angular beamwidth and $J_1$ and $J_3$ are
the Bessel functions of the first kind, of order one and three respectively.
Our interest is the parameter $\alpha(\theta_E)$ in the specific angle $\theta_E$ between Bob's and Eve's directions.
$g_{L}$ and $g_{E}$ are legitimate and eavesdropper's receiver's antenna gains towards Earth station.
Now we introduce a model for the coefficient that gives Eve's signal strength relative to Bob's signal in \eqref{eq:SysModel} and \eqref{eq:GaussChannel} (see \cite{PhySecSat}\cite{arXiv}\cite{RealisticChannel}). We first introduce the parameter $\mu$ to account for the relative antenna gain, i.e.,
$\mu:= \sqrt{g_{L}/g_{E}}$.
We also define $\beta\left(r,\rho_{E}\right)$ to account for relative propagation losses between Bob and Eve as
\begin{align}\label{H33}
\beta^2\left(r,\rho_{E}\right)=\frac{\rho_{L}^{2}}{\rho_{E}^{r}}.
\end{align}
The exponent $r$ accounts for the power attenuation decay that affects eavesdropper's propagation channel. Different values of the exponent model
correspond to different assumptions about eavesdropper.
Specifically, the eavesdropper can be modeled as a terrestrial, aerial or satellite station.
For example, while for the satellite case $r=2$, in case of aerial eavesdropper,
a good assumption is to consider a large scale two-ray ground multipath model,with $r>2$.
Then, we discuss
the parameter $\gamma_B$ in the OW \eqref{eq:SysModel} and
the parameter $\gamma_A$ in the TW \eqref{eq:GaussChannel}
because
the channel \eqref{eq:SysModel} in the OW case
is the same as
the channel \eqref{eq:GaussChannel} in the TW case
in each scenario (I) or (II).
These parameters are given as
\begin{align}
\gamma(\theta_E, \rho_E, r,\mu,\gamma_n): =
\frac{\alpha\left(\theta_{E}\right)
\mu\beta\left(r,\rho_{E}\right)}{\sqrt{\gamma_n}},
\label{grrt}
\end{align}
where, $\gamma_n$ is the ratio between the powers of the noises in
legitimate receiver's and eavesdropper's detectors.
In doing a secrecy analysis, we assume in which orbit Eve is,
but we don't make any assumption on which angle she is (since she is orbiting).
In this case, to guarantee the security,
we need to consider the worst case.
For this aim, we consider the possible range ${\cal R}$ of the value $(\theta_E, \rho_E)$.
Then, the maximums of $\gamma_B$  and $\gamma_A$ are calculated to
\begin{align}
\gamma_{B,\max}=\gamma_{A,\max}=
\max_{(\theta_E, \rho_E) \in {\cal R}}
\gamma(\theta_E, \rho_E, r,\mu,\gamma_n).
\end{align}

Now, we assume representative values for Eve's possible orbits according to basic orbital mechanics \cite{OrbitMechanics} and usual low or medium orbit terminology. In Case (I),
when Eve is MEO (LEO), we assume that the height of Eve's orbit runs from
$\rho_{\mathop{\rm MEO}\min}=5000$ to $\rho_{\mathop{\rm MEO}\max}=20000$ km
($\rho_{\mathop{\rm LEO}\min}=150$ to $\rho_{\mathop{\rm LEO}\max}=2000$ km).
Also, we assume that the height of our GEO orbit is $\rho_{\mathop{\rm GEO}}= 36000$ km.
Since the maximum of $\alpha\left(\theta_{E}\right)$ is realized by $\theta_{E}=0^{\circ}$,
when Eve is MEO (LEO),
we have
$\gamma(0^{\circ}, \rho_E, r,\mu,\gamma_n)
= \frac{\mu \rho_{\mathop{\rm GEO}}}{ \rho_{\mathop{\rm MEO}\min}^{r/2}\sqrt{\gamma_n}}$
($= \frac{\mu \rho_{\mathop{\rm GEO}}}{ \rho_{\mathop{\rm LEO}\min}^{r/2}\sqrt{\gamma_n}}$).
For illustration, we now assume Eve equally powerful than the legitimate receiver, i.e. $\mu=1$ and $\gamma_n=1$. Also, we have $r=2$ for both LEO and MEO.
Hence, the maximum $\gamma_{B,\max}=\gamma_{A,\max}$
is given as
$\gamma_{{\rm (I)}\mathop{\rm MEO}}:=
\frac{ \rho_{\mathop{\rm GEO}}}{ \rho_{\mathop{\rm MEO}\min} }= 36000 / 5000 = 7.2$
for the case when Eve is MEO,
and it is given as
$\gamma_{{\rm (I)}\mathop{\rm LEO}}:=
\frac{ \rho_{\mathop{\rm GEO}}}{ \rho_{\mathop{\rm LEO}\min} }= 36000/150 = 240$
for the case when Eve is LEO.
Then, we have the capacities of the worst case as Table \ref{T1}.
While these capacities are small in comparison with the conventional communication,
we see that the secure communication is possible in these scenarios only in the TW case.

Applying same reasoning in Case (II), when Eve is MEO,
the minimum of $\rho_E$ is $\rho_{\mathop{\rm GEO}}-\rho_{\mathop{\rm MEO}\max}$ at $\theta_E=0^{\circ}$
and the maximum of $\alpha\left(\theta_{E}\right)$ is realized by $\theta_{E}=0^{\circ}$.
The same observation holds when Eve is LEO or GEO.
In this case, it seems
reasonable to assume the legitimate receiver having a more
powerful antenna gain and less detector noise than the eavesdropper.
Hence, we can again assume $\mu=1$ and $\gamma_n=1$.
Again, we also have $r=2$.
Therefore, when Eve is MEO,
the maximum $\gamma_{B,\max}=\gamma_{A,\max}$
is calculated to be
$\gamma_{{\rm (II)}\mathop{\rm MEO}}:=
\frac{\rho_{\mathop{\rm GEO}}}{ (\rho_{\mathop{\rm GEO}}-\rho_{\mathop{\rm MEO}\max})}=36000/(36000-20000)=9/4$.
When Eve is LEO,
it is calculated to be
$\gamma_{{\rm (II)}\mathop{\rm LEO}}:=
\frac{ \rho_{\mathop{\rm GEO}}}{ (\rho_{\mathop{\rm GEO}}-\rho_{\mathop{\rm LEO}\max}) }=36000/(36000-2000)=18/17$.
When Eve is GEO,
the minimum distance between the transmitter of the initial transmission and Eve
is 1km.
Hence,
it is calculated as
$\gamma_{{\rm (II)}\mathop{\rm GEO}}:=
\frac{\rho_{\mathop{\rm GEO}}}{ 1 }=36000$.

\begin{table}[htpb]
  \caption{Summary of optimum capacities in satellite models}
\label{T1}
\begin{center}
  \begin{tabular}{|l|c|c|c|c|}
\hline
$\gamma_{B,\max}=\gamma_{A,\max}$
 & $C_{\mathop{\rm soft}}^{\mathop{\rm TW}}$ &  $C_{\mathop{\rm hard}}^{\mathop{\rm TW}}$ & $C_{\mathop{\rm soft}}^{\mathop{\rm OW}}$ &  $C_{\mathop{\rm hard}}^{\mathop{\rm OW}}$ \\
\hline
$\gamma_{{\rm (I)}\mathop{\rm MEO}}$ & $0.029$ & $0.036$&  0&0 \\
\hline
$\gamma_{{\rm (I)}\mathop{\rm LEO}}$ & $3.6  \times 10^{-4}$ & $3.9  \times 10^{-4}$ & 0&0 \\
\hline
$\gamma_{{\rm (II)}\mathop{\rm GEO}}$ & $1.4  \times 10^{-12}$ & $1.4  \times 10^{-12}$ &  0&0 \\
\hline
$\gamma_{{\rm (II)}\mathop{\rm MEO}}$ & $0.086$ & $0.108$ &  0&0 \\
\hline
$\gamma_{{\rm (II)}\mathop{\rm LEO}}$ & $0.159$ & $0.198$ &  0&0\\
\hline
  \end{tabular}
\end{center}
\end{table}

\section{Conclusions and further improvements}
We have introduced a two-way protocol for the BPSK modulation
to overcome the limitations of the classic (one-way) wiretap physical layer security protocol.
While the secrecy capacity of the one-way protocol is negative when Eve's channel is better than Bob's channel (i.e $\gamma_B > 1$), we show that the two-way protocol always provides positive capacity with higher gains even for
$\gamma_A \ge 1$ and $\gamma_B \ge 1$ when the noises exist and are independent.
We have shown that
the one-way protocol cannot realize secure communication in realistic scenarios of satellite communication,
while our two-way protocol can realize secure communication in these realistic scenarios.
Notice that this conclusion does not change whenever the maximum of possible values of
$\gamma_A $ and $\gamma_B $ is greater than $1$.

For example, in the scenario (I), we have the transmission rate 3.6$\times 10^{-4}$ with the worst case analysis
in the two-way protocol with Eve's soft decision
while the eavesdropper has an extremely stronger power in the detection process than the legitimate receiver, i.e., the eavesdropper has 240 times power in the receiving signal as the legitimate receiver.
The conventional one-way method cannot realize secure communication in this case.
This numerical analysis shows that even in this case,
we can realize secure communication with the same physical device
if we accept approximately e.g. one thousandth reduction of the speed of the conventional communication without secrecy (e.g., 1 Gbps is reduced to 1 Mbps).
While this cost seems very large,
the cost is still much smaller than quantum key distribution (QKD) \cite{BB84} due to the following reason.
Since QKD requires expensive devices, it is available only for extremely limited users
(big governments and/or big military organizations).
However, since the proposed method is based on the conventional satellite system,
even though the transmission speed is very low,
it is available for ordinary users.
In fact, the user can use the proposed system
when the size of communication is reduced, e.g., the user uses only e-mail instead of video.
On the other hand, in the scenario (II),
we assume that the receiving antenna of the earth station has almost same performance
as that of Eve.
Under such a worst case assumption, we obtain very small transmission speed.
To improve this, the receiving antenna of the earth station needs to be more powerful than that of Eve.
However, to realize this condition,
the earth station needs to prepare an expensive receiving device, which may or may not imply to restrict ordinary users since such higher cost is shared by all service/users sharing the earth station.
In this sense, the scenario (II) may seem less practical for ordinary users.
In any case, to realize secure communication, it is sufficient to share secure keys between two users.
because one-time use of shared secure key realizes secure communication with both directions.
Hence, it is sufficient to establish secure communication only with one direction.
Therefore, the scenario (I) of TW is enough for our purpose.

However, in the scenario (I) of TW, if Eve is an eavesdropping terrestrial node, e.g., a drone near the terrestrial earth station, she has better performance than LEO/MEO satellite.
In this case, the secure transmission rate is worse than 3.6$\times 10^{-4}$
when the angle $\theta_E$ is set to be $0$.
However, the possible minimum angle $\theta_E$ of this case in practice is larger than
that in the case with an eavesdropping LEO/MEO satellite.
Therefore, it is needed to evaluate the secure transmission rate with
an eavesdropping terrestrial node and
the possible minimum angle $\theta_E$.
However, it is not so easy to find the minimum angle $\theta_E$
among practically possible values.
Therefore, this type of analysis is remained as a future study.

As the price to pay, the protocol requires higher delay to establish the secure channel when compared to the one-way.
However, this cost is much cheaper than
the previous two-way protocols in the papers \cite{WGH,FJHWY,ZWPT,TY}
because
they require many rounds of communication while our protocol requires only two rounds of communication.
On the other hand, the transmission of information can be over a public channel while for randomness sharing, the channel needs to be previously authenticated like \cite{Correlation,Fung}.
As discussed in \cite{Correlation,Fung}, the required amount of the random numbers shared  between Alice and Bob in advance is the logarithm order of the size of intended secure communication.
Also, this cost is much cheaper than the realization of quantum key distribution. Therefore, considering the cost-benefit performance, we find that our two-way method is useful.
Furthermore, the bias of the variable $B$ may reduce the effectiveness of the protocol and reduce the secrecy capacity gains.
To improve this problem, we often distill uniform random numbers from the thermal noise.
It is known that it is possible to distill uniform random numbers
by applying a hash function to a biased random numbers \cite{BBCM,HILL,Hayashi2011,H13}.
To obtain the ultimate secure uniform random number, we may employ quantum random number generator \cite{HG,MYQZ,HZ},
which requires much cheaper cost than
quantum key distribution because it does not need quantum communication.

Unfortunately, this paper discusses only the asymptotic performance.
Since the implemented communication system has finite-length codes,
we need to evaluate the security of finite-length codes for its practical application \cite{ISIT2013,HTW1,HTW2,YSP,HarXiv}.
Since the finite-length analysis depends on the choice of the security criterion,
we need to be careful of its choice \cite{hay-wire,H13}.
As such a study is beyond the focus on this paper,
it is considered as a future study.
Furthermore, it is well known that a good model for the land mobile satellite (LMS) channel model is a Markov model \cite{Fontan}.
Hence, it is a completely different channel from \eqref{eq:GaussChannel}.
Therefore, follow up studies also include considering different satellite channel models such as fading models accounting for frequency-dependent atmospheric effects and for the case of mobile (legitimate) users.


\section{Introduction}
The world we live in is so rich and varied. But there must be some invariable factors that enable people to identify the same object from the change.
Among these complex transformations in nature, affine transformation is one of the most basic and common issue in computer vision, which is the largest set of linear transformations, including rotation, translation, scaling, et cetera.
Geometric affine transformation in shape space can be used to approximate the projective transformations in images where the camera happened a slight change of viewing perspective, which is the most common situation in the imaging process.
Photometric affine transformation in channel space is the best linear model to describe the color degradation of the outdoor scenes, which has been proved in \cite{geusebroek2001color}.
And affine transformations often occur in both shape and color space simultaneously, which can be termed as \textit{Dual-Affine Transformation} (DAT). DAT is a useful and common transformation model to describe the degradations happened in vision information from our diverse and colorful world, which considering both shape and color degradation.
Invariants under DAT is a kind of important image feature that can be used for object classification and object retrieval.

With the diversification of visual data acquisition means, the formats of visual data are also diversified, which exist in different dimension and with different number of color channels, such as 2D color images, 3D color objects, 4D movies of a 3D color object, even higher-dimensional data and multi-spectrum images with multi-channels. Other data organized in channels (such as vector fields) can also be included in this work. Some examples are shown in Figure \ref{fig:data}.

In common, we can derive dual-affine moment invariants for objects in a particular dimension with a certain number of channels separately. To the best of our knowledge, there is no general framework to derive invariants for all of these data formats. In this paper,
we systematically analyze the properties of invariants under DAT for different data dimensions and present a general framework for deriving invariants under DAT for objects in $M$-dimensional shape space with $N$ color channels.

\subsection{Object definition in $M$-dimensional space with $N$ channels}

An object in $M$-dimensional ($M$-D) shape space, ($X,Y,Z,S4,S5,...$), with $N$ channels, ($R,G,B,C4,C5,...$), can be described as Equation (\ref{equ:fun}):
\begin{align}\label{equ:fun}
O_{M,N} = f_{M,N}(\vec{X}) = (&R(\vec{X}); G(\vec{X});B(\vec{X}); C4(\vec{X});C5(\vec{X});...)
\end{align}
where $\vec{X}$ represents the coordinate $(x,y,z,s4,s5,...)^T$ in $M$-D space.
And $R(\cdot), G(\cdot), B(\cdot), C4(\cdot), C5(\cdot)$ represent the functions in different channels.
Object $O$ can be defined in another form as Equation (\ref{equ:obj}):
\begin{align}\label{equ:obj}
O_{M,N}(id) = &[S_{M}(id),C_N(id)] \notag \\
=&[(x_{id},y_{id},z_{id},s4_{id},s5_{id},...)^T,(r_{id},g_{id},b_{id},c4_{id},c5_{id}...)^T]
\end{align}
where $id$ represents the index of pixels of object $O$. $S_M$ represents the coordinate in shape space and $C_N$ represents the values in channel space.

Lots of data are organized in channel forms, which can represent different meanings in different data, such colors in color objects and vectors in vector fields.
When we specify $M, N$ to different values, we can get different data formats.
For examples, gray image is object in 2D space with 1 channel where $M=2, N=1$.
2D vector field can be regarded as object in 2D space with 2 channels where $M=2, N=2$. The values of 2 channels represent the vector at each point.
Color image is object in 2D space with 3 channels where $M=2, N=3$.
Both 3D vector field and 3D Point Cloud are objects in 3D space with 3 channels. The difference is that values of 3 channels in 3D vector field represent the vector at each point and that in 3D Point Cloud represent the color at each point.
Some objects in different formats are in Figure \ref{fig:data}.

\begin{figure}
	\centering
	\subfigure[]{
		\includegraphics[width=0.18\linewidth]{figures/lena2.jpg}}
	\subfigure[]{
		\includegraphics[width=0.18\linewidth]{figures/vf.png}}
	\subfigure[]{
		\includegraphics[width=0.18\linewidth]{figures/lena.jpg}}
	\subfigure[]{
		\includegraphics[height=0.065\textheight]{figures/vf3.png}}
	\subfigure[]{
		\includegraphics[height=0.065\textheight]{figures/0x.png}}
	\caption{Examples of different data formats: (a) Gray image in 2D space with 1 channel ($M=2, N=1$); (b) Vector field in 2D space ($M=2, N=2$); (c) Color image in 2D space with 3 channels ($M=2, N=3$); (d) Vector field in 3D space ($M=3, N=3$); (e) Color object in Point Cloud in 3D space with 3 channels ($M=3, N=3$).}
	\label{fig:data}
\end{figure}

\subsection{Dual-Affine Transformation (DAT)}

The definition of $M$-D Affine Transformation Matrix (ATM) is listed in Equation (\ref{equ:ts}), which is a $M\times M$ matrix:

\begin{equation} \label{equ:ts}
ATM_M = \left(
\begin{array}{ccccc}
a_{11} & a_{12} & a_{13} & ... & a_{1M}\\
a_{21} & a_{22} & a_{23} & ... & a_{2M}\\
a_{31} & a_{32} & a_{33} & ... & a_{3M}\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
a_{M1} & a_{M2} & a_{M3} & ... & a_{MM}\\
\end{array}
\right)_{M\times M}
\end{equation}

After $M$-D affine transformation, the coordinates of object pixels in space $S_M=(x,y,z,s4,s5,...,sM)^T$ should be transformed to  $S^*_M=(x^*,y^*,z^*, s4^*, s5^*,...sM^*)^T$, which can be expressed by Equation (\ref{equ:tsxyz}):
\begin{equation}\label{equ:tsxyz}
S^*_M = \left(\begin{array}{c}
x^* \\ y^* \\ z^* \\ s4^* \\ s5^* \\ ... \\ sM^*
\end{array}\right)
=ATM_M \cdot
\left(\begin{array}{c}
x \\ y \\ z \\ s4 \\ s5 \\ ... \\ sM
\end{array}\right) = ATM_M \cdot S_M
\end{equation}

After $N$-D affine transformation, the pixel values at different channels $C_N=(r,g,b,c4,c5,...)^T$ will be transformed to $C^*_N=(r^*,g^*,b^*,c4^*,c5^*,...)^T$ in the same way as Equation (\ref{equ:tcrgb}):

\begin{equation}\label{equ:tcrgb}
C^*_N = ATM_N \cdot C_N
\end{equation}

Affine transformations in shape space and channel space often occur simultaneously, which can be called \textit{Dual-Affine Transformation} (DAT). DAT is the transformation we focus on in this article.
The formal definition of DAT for $M$-D objects with $N$ channels can be described as Equation (\ref{equ:dual}):
\begin{equation} \label{equ:dual}
f^*_{M,N}(\vec{X}) = ATM_N \cdot f_{M,N}(ATM_M^{-1} \cdot \vec{X})
\end{equation}
where $ATM_M$ is the \textit{Spatial Transformation} and $ATM_N$ is the \textit{Channel Transformation}. We depict the relationship between them with DAT in Figure \ref{fig:dat}.
DAT is equivalent to applying a \textit{Spatial Transformation} and a \textit{Channel Transformation} on the object independently.

\begin{figure}[t]
	\centering
	\includegraphics[width=\linewidth]{figures/dat3.png}
	\caption{The relationship between \textit{Spatial Transformation}, \textit{Channel Transformation} and DAT with an example of a 3D color Point Cloud ($M=3,N=3$). }
	\label{fig:dat}
\end{figure}

One typical case of DAT is for color images, where $M=2, N=3$. The 2D geometric affine transformation in space represents the deformation of object's shape in the images, and the 3D photometric affine transformation in channel space represents the deformation of object's color. Moment invariants under DAT of color images has been fully discussed in \cite{gong2017naturally}.
Another case is invariants under DAT in 2D vector fields, which have been presented in \cite{schlemmer2007moment,kostkova2018affine,j2019affine}. There are also researches about 3D moment invariants under single affine transformation for 3D objects with one channel in \cite{lo19893,xu2008geometric}.
However, with the continuous advancement of data acquisition methods, more and more forms of data have emerged, like 3D color objects, 4D or higher-dimensional data, multi-spectral images with more than 3 channels, etc. Invariant descriptors under dual-affine transformations in these domains are highly desirable.

The main contribution of this paper is that, we present a general framework for the derivation of moment invariants under DAT for $M$-dimensional objects with $N$ channels.
And with this framework, we get a set of independent and usable DAMI for 3D color objects with 3 channels.
The paper is organized as follows. First, we give a survey of relevant literature about DAMI in Section 2. In Section 3, we illuminate the DAMI's derived framework in detail and get the DAMI generation formula for objects in any data format.
In Section 4,  we instantiate DAMI with data format of 3D color point cloud as an example.
Through a visual approach, we generate a complete set of DAMI equations with $M=3, N=3$.
We perform independent analysis on all DAMI equations and obtain a set of independent DAMI.
In Section 5, we evaluate the ability of DAMI through numerical experiments on several data sets.

\section{Related Works}
Moment invariants were introduced into visual pattern recognition by \cite{hu1962visual} for character recognition, which are invariant under 2D \textit{similarity transformation}.
These seven invariants are considered as a kind of shape descriptor, which plays an important role in the field of pattern recognition. Since then, various kinds of moment invariants have been proposed for different data formats.
Researchers in this field mainly focus on the following three directions: (i) The first is to extend the transform group of invariants to more complex ones. (ii) The second is to derive invariants for data of different dimensions or formats. (iii) The third is to derive invariants including both shape and color degradations into consideration.

For gray-scale images, Flusser and Suk extended Hu's moment invariants into \textit{affine transformation} in \cite{flusser1993pattern}, which are more common in nature than similar transformation.
Relevant works are also demonstrated in \cite{flusser1994affine,suk2004graph,suk2011affine}.
Invariants under \textit{projective transformation}, which is the exact transformation model of the imaging process, has been analyzed in \cite{weiss1988projective,suk2004projective,li2018image}.
There are also some other interesting transformations that not only for image data. Invariants under these transformations are also an important direction in the field, like chiral invariants in \cite{zhang2017fast}, invariants under \textit{M\"obius transformation} in \cite{zhang2018differential}, and even some researches about invariants under \textit{conformal}, \textit{quasi-conformal}, and \textit{diffeomorphism transformations} in \cite{white1973global,willmore2000surfaces,zeng2011registration,rustamov2013map}.
Invariant primitive is the most important part of the derivation process of invariants and is the essential invariant property under a certain transform group. Invariant construction frameworks are always based on primitives.
Under affine transformation in 2D space, the cross-product of the coordinates of two points is the covariant primitive, which is used to derive affine invariant for gray-scale images in \cite{suk2004graph}. In \cite{xu2008geometric}, the distances, areas and volumes of coordinates of points are considered as the covariant primitives for different transformations. The essential primitive under affine transformation in $n$-D space is considered as the determination of the coordinates of $n$ pixels in \cite{mamistvalov1998n}.
It is the same but different format with the area of triangle in 2D space and the volume of tetrahedron in 3D space.
\cite{li2017shape} analysis the generating formula of invariants under different transformations in detail, which can be considered as the Shape-DNA.

There are some researches extend moment invariants into different data formats.
Some works analyze shape features by extracting invariants for curves in 2D space like \cite{mundy1992geometric,zhao1997affine}.
For the 3D objects with one channel, researches about 3D affine moment invariants can be found in \cite{lo19893,novotni20033d,xu2008geometric}.
The invariant for n-D space has been studied in \cite{mamistvalov1998n}.

Color information is important for the recognition of objects.
Some color descriptors are simple and useful for applications, like color histogram in \cite{funt1995color}, and color moment in \cite{stricker1995similarity}. But they do not consider the color degradation model in reality. Some other color descriptors consider the degradation model of illumination, like works in \cite{gevers1999color,li2009illumination,gong2013moment}.
To consider shape and color degradation in the same framework, there are lots of works. \cite{mindru2004moment} proposed a kind of color moment invariant with Lie group method for only color diagonal transformation.
\cite{gong2017naturally} proposed moment invariants under dual-affine transformation for 2D color images.  Some other invariants for color images can be found in \cite{alferez1999geometric,mindru2004moment,suk2009affine,mo2017shape}.
Affine transformation is considered as the best model to describe the outdoor color degradation statistically in \cite{geusebroek2001color}.
In this paper, we will present a general framework for the derivation of moment invariants under \textit{dual-affine transformation} for $M$-D objects with $N$ channels. And with 3D color point cloud as an example, we get a set of independent and usable DAMI with $M=3, N=3$.

\section{Dual-Affine Moment Invariants (DAMI)}
The flowchart of the derivation, instantiation, and usage of DAMI is depicted in Figure \ref{fig:flow}. We will detail the DAMI derivation framework in this section and leave the instantiation and usage part in Section 4 and 5.
The general derivation framework for moment invariants is to first find the \textit{covariant primitive} under a given transformation, and then derive the invariants in the global domain through the integral framework and the normalization factor. Firstly, we give the definition of \textit{covariant primitive} (\textbf{Step 1}) in Section 3.1 and corresponding $SpaceKernel$ and $ChannelKernel$ under $M$-D affine transformation (\textbf{Step 2}) in Section 3.2.
With the proposed $SpaceKernel$ and $ChannelKernel$, we can derive the moment covariants under DAT (\textbf{Step 3}) through integration in Section 3.3. In Section 3.4, we get the generating formula of DAMI after applied normalization (\textbf{Step 4}). Finally, we analysis the null space and a special case of DAMI In Section 3.5 and 3.6.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.85\linewidth]{figures/dami3.png}
	\caption{The flowchart of DAMI derivation, instantiation and usage. }
	\label{fig:flow}
\end{figure}

\subsection{Covariant primitive}
\begin{theorem}
	The matrix determinant of the coordinates of $M$ points in $M$-D space is an affine covariant primitive.
\end{theorem}

The determinant of coordinates of $M$ points in $M$-D space,  $(x_1,y_1,z_1, (s4)_1,..., (sM)_1)^T$,  $(x_2,y_2,z_2,(s4)_2,...,(sM)_2)^T$, ..., $(x_M,y_M,z_M,(s4)_M,...,(sM)_M)^T$,  is the covariant primitive $P_S(M)$ in shape space as Equation (\ref{equ:sc}):
\begin{equation}\label{equ:sc}
P_{S}(M) = \det(MS(M)) = \left|
\begin{array}{cccc}
x_1 & x_2 & \cdots & x_M \\
y_1 & y_2 & \cdots & y_M \\
z_1 & z_2 & \cdots & z_M \\
(s4)_1 & (s4)_2 & \cdots & (s4)_M \\
\vdots & \vdots & \ddots & \vdots \\
(sM)_1 & (sM)_2 & \cdots & (sM)_M \\
\end{array}
\right|
\end{equation}
After $M$-D affine transformation $ATM(M)$, the transformed  $P_{S}^*(M)$ will be Equation (\ref{equ:tsc}):

\begin{align}\label{equ:tsc}
P_S^*(M) &= \left|MS(M)^* \right| = \left|ATM(M) \cdot MS(M) \right| \notag \\
&= \left|ATM(M)\right| \cdot \left|MS(M) \right| = \left|ATM(M)\right| \cdot P_S(M)
\end{align}
where $MS(M)^*$ is the matrix of the transformed coordinates as Equation (\ref{equ:ms}):
\begin{equation}\label{equ:ms}
MS(M)^* = \left(
\begin{array}{cccc}
x_1^* & x_2^* & \cdots & x_M^* \\
y_1^* & y_2^* & \cdots & y_M^* \\
z_1^* & z_2^* & \cdots & z_M^* \\
(s4)_1^* & (s4)_2^* & \cdots & (s4)_M^* \\
\vdots & \vdots & \ddots & \vdots \\
(sM)_1^* & (sM)_2^* & \cdots & (sM)_M^* \\
\end{array}
\right)
\end{equation}
From the result we can see that the determinant before and after the transformation only differs by one scalar coefficient, which is the determinant of the transformation matrix $ATM(M)$. $P_S(M)$ is a covariant primitive in $M$-D space. \textbf{Theorem 1} is proved.

And similarly,  in $N$-D channel space, the matrix determinant of $N$ points with $N$ channels, $(r_1,g_1,b_1,(c4)_1,...,(cN)_1)^T$, $(r_2,g_2,b_2,(c4)_2,...,(cN)_2)^T$, ..., $(r_N,g_N,b_N,(c4)_N,...,(cN)_N)^T$, is covariant primitive in channel space as Equation (\ref{equ:cc}):
\begin{equation}\label{equ:cc}
P_{C}(N)=\det(MC(N))
\end{equation}
where $MC(N)$ is the matrix of the channels' values of $N$ points.
And after $N$-D  affine transformation, the transformed $P^*_{C}(N)$ will be Equation (\ref{equ:tcc}):
\begin{align}\label{equ:tcc}
P_C^*(N) &= \left|MC(N)^* \right| = \left|ATM(N) \cdot MC(N) \right| \notag \\
&= \left|ATM(N)\right| \cdot \left|MC(N) \right| = \left|ATM(N)\right| \cdot P_C(N)
\end{align}

From the results, we can see that the matrix determinant $P_C(N)$ is covariant under affine transformation $ATM(N)$ in $N$-D channel space.

\subsection{Space Kernel and Channel Kernel}

Further, the continued product of several these covariant primitives  is also covariant, as shown in Equation (\ref{equ:kernel}, \ref{equ:kernel2}).

\begin{align} \label{equ:kernel}
SpaceKernel(M,P,O_S) = \prod_{i} P_S(M_i)^{t_i} \\
ChannelKernel(N,Q,O_C) = \prod_{j} P_C(N_j)^{t_j} \label{equ:kernel2}
\end{align}
In which, $M_i$ and $N_j$ represent  $M$ or $N$ different points and different $i$ and $j$ represent different combinations of points. All points involved in building $ SpaceKernel$ constitute the Point Set $ PS_S $. And all  points involved in building $ ChannelKernel$ constitute the Point Set $ PS_C $.
The number of points in the point set,  $P = \left|PS_S\right|$ and $ Q = \left|PS_C\right| $, are called the \textbf{degree}.
And $t_i$ and $t_j$ are the order of corresponding primitives separately.
And the $O_S = \sum t_{i}$ and $O_C = \sum t_{j} $ are called the \textbf{order}.

\subsection{Covariant through integration}

To define the $Covariant$ on object $O_{M,N}$, we make an integration of the $SpaceKernel$ and $ChannelKernel$ as Equation (\ref{equ:covariant}):
\begin{align} \label{equ:covariant}
& Covariant(M,N,P,Q,O_S,O_C)(O) \notag \\
&\quad = \iint...\iint SpaceKernel(M,P,O_S) \cdot \notag \\
&\quad\quad\quad\quad\quad\quad\quad\quad ChannelKernel(N,Q,O_C) \cdot \notag \\
&\quad\quad\quad\quad\quad\quad\quad\quad dx_1dy_1dz_1 .... ... dx_Ldy_Ldz_L...
\end{align}
where $L=|PS_S \cup PS_C|$.

The relationship between $Covarinat$ before and after $M$-D shape affine and $N$-D channel affine transformation is in Equation (\ref{equ:covariantt}):
\begin{align}\label{equ:covariantt}
&Covariant^*(M,N,P,Q,O_S,O_C)(O^*) \notag \\
&\quad = \iint...\iint \prod_i P^*_S(M_i)^{t_i}  \cdot \prod_j P^*_C(N_j)^{t_j} \cdot \notag \\
&\quad\quad\quad\quad\quad\quad\quad\quad dx_1^*dy_1^*dz_1^* .... ... dx_L^*dy_L^*dz_L^*... \notag \\
&\quad = \left| ATM(M) \right|^{O_S+L} \cdot \left| ATM(N) \right|^{O_C} \cdot \notag \\
&\qquad \iint...\iint \prod_i P_S(M_i)^{t_i}  \cdot \prod_j P_C(N_j)^{t_j} \cdot \notag \\
&\quad\quad\quad\quad\quad\quad\quad\quad dx_1dy_1dz_1 .... ... dx_Ldy_Ldz_L...  \notag \\
&\quad = \left| ATM(M) \right|^{O_S+L} \cdot \left| ATM(N) \right|^{O_C} \cdot \notag \\
&\qquad\qquad Covariant(M,N,P,Q,O_S,O_C)(O)
\end{align}
We can see that the difference between $Covariant$ and $Covariant^*$ is only one scalar coefficient, which dependents on the transformation matrix $ATM(M)$ and $ATM(N)$. The $Covariant$ is covariant under DAT.

\subsection{Invariant through normalization}

To achieve $Invariant$ on object $O_{M,N}$, we should eliminate the covariant scalar of the $Covariant$. We use a $NormalizationFactor$ consists of some special $Covariant$ to divide the $Covariant$ in Equation (\ref{equ:covariant}).
	Then we will get $Invariant$ under DAT in $M$-D space with $N$ channels as Equation (\ref{equ:invariant}):
	\begin{align}\label{equ:invariant}
	&Invariant(M,N,P,Q,O_S,O_C)(O) \notag \\
	&\qquad = \frac{
		Covariant(M,N,P,Q,O_S,O_C)
	}{
		NormalizationFactor
	} \notag \\
	&\qquad = \frac{Covariant(M,N,P,Q,O_S,O_C)}{
		\left\{
		\splitdfrac{Covariant(M,N,1,0,0,0)^{O_S+L-N\cdot O_C/2} \cdot}{ Covariant(M,N,1,N,0,2)^{O_C/2}}
		\right\}}
	\end{align}
	The \textbf{degree} is $L=|PS_S \cup PS_C|$, and the \textbf{order} is $K=O_S+O_C$.
	The $NormalizationFactor$ can be any other one that can eliminate the covariant scalar in the numerator.

This definition can be proved in Equation (\ref{equ:invariantt}):

\begin{align}\label{equ:invariantt}
&Invariant^*(M,N,P,Q,O_S,O_C)(O^*) \notag \\
& = \frac{
	Covariant^*(M,N,P,Q,O_S,O_C)
}{
\left\{
	\splitdfrac{
	Covariant^*(M,N,1,0,0,0)^{O_S+L-N \cdot O_C/2} \cdot}{ Covariant^*(M,N,1,N,0,2)^{O_C/2}}
\right\}
} \notag \\
&=  \frac{
	\left\{
	\splitdfrac{
	\left| ATM(M) \right|^{O_S+L} \cdot \left| ATM(N) \right|^{O_C} \cdot }{ Covariant(M,N,P,Q,O_S,O_C)}
	\right\}
}{\left\{
	\splitdfrac{(\left| ATM(M) \right|Covariant(M,N,1,0,0,0))^{O_S+L-N\cdot O_C/2} \cdot}
	{(\left| ATM(M) \right|^{N} \cdot \left| ATM(N) \right|^{2} Covariant(M,N,1,N,0,2))^{O_C/2}}
	\right\} } \notag \\
& = \frac{Covariant(M,N,P,Q,O_S,O_C)}{
\left\{
\splitdfrac{
	Covariant(M,N,1,0,0,0)^{O_S+L-N\cdot O_C/2} \cdot }{ Covariant(M,N,1,N,0,2)^{O_C/2}} \right\}} \notag \\
& = 	Invariant(M,N,P,Q,O_S,O_C)(O)
\end{align}
The invariants before and after the transformation are completely equal.
$Invariant(M,N,P,Q,O_S,O_C)$ is invariant under DAT in the $M$-D space with $N$ channels, which can be called Dual-Affine Moment Invariants ($DAMI(M,N)$).

\subsection{Null space of DAMI}

As we can see from the definition of covariant primitive in Equation (\ref{equ:sc}) and (\ref{equ:cc}), when the determinant of all primitives equal to 0, the DAMI will be no definition, which can be called the null space of DAMI.

When the null space happened in the primitives in shape space, it means matrix of the coordinations is singular. And some axis is linear combination of the other axises. The object will be part of a  low-dimensional hyperplane in the original space, where the DMAI has no definition. For this kind of objects, we can find all the independent axises first, say there are $M'$ , and then derive $DAMI(M', N)$ on these axises. That will make the DAMI work again.

When the null space happened in the primitives in channel space, it will be the same with that in shape space. Some channels are linear dependent on others. We can find all the independent channels, say there are $N'$, and then derive $DAMI(M,N')$ on these independent channels. That will make the DAMI work again.  The literature \cite{kostkova2019null} has discussed some about the null space of color space. This situation is rarely seen in RGB images, but may exist in multi-spectral images.

\subsection{Special case}
When $M=2,N=2$, $DAMI(2,2)$ can be used on 2D vector fields with 2 channels. The $ATM(M)$ will be the coordinate transformation of the vector field and the $ATM(N)$ will be the vector transformation at each point. The DAMI will be as Equation (\ref{equ:ami22}):

\begin{align}\label{equ:ami22}
&DAMI(2,2)=Invariant(2,2,P,Q,O_S,O_C)(O) = \notag \\
&\qquad = \frac{
	Covariant^*(2,2,P,Q,O_S,O_C)
}{
\left\{
\splitdfrac{
	Covariant^*(2,2,1,0,0,0)^{O_S+L-O_C}\cdot}{ Covariant^*(2,2,1,2,0,2)^{O_C/2}}
\right\}
}
\end{align}

When select the same parameters with that in \cite{kostkova2018affine,j2019affine}, we can get the same invariants as Equation (\ref{equ:case22}).
\begin{align}\label{equ:case22}
&DAMI(2,2)(1) = u_{0101}*u_{1010} - u_{0110}*u_{1001}  \notag \\
&DAMI(2,2)(2) = u_{0301}*u_{3010} - u_{0310}*u_{3001} - 3*u_{1201}*u_{2110} + \notag \\
&\quad 3*u_{1210}*u_{2101} \notag \\
&DAMI(2,2)(3) = u_{0501}*u_{5010} - u_{0510}*u_{5001} - 5*u_{1401}*u_{4110} + \notag \\
&\quad 5*u_{1410}*u_{4101} + 10*u_{2301}*u_{3210} - 10*u_{2310}*u_{3201} \notag \\
&DAMI(2,2)(4) = u_{0201}^2*u_{2010}^2 - 2*u_{0201}*u_{0210}*u_{2001}*u_{2010} - \notag \\
&\quad 4*u_{0201}*u_{1101}*u_{1110}*u_{2010} + 4*u_{0201}*u_{1110}^2*u_{2001} + \notag \\
&\quad u_{0210}^2*u_{2001}^2 + 4*u_{0210}*u_{1101}^2*u_{2010} - \notag \\
&\quad 4*u_{0210}*u_{1101}*u_{1110}*u_{2001}
\end{align}

\section{Instantiation of $DAMI(3,3)$}
In this section, we will use the $DAMI(3,3)$ for 3D color objects as an example to demonstrate how to instantiate DAMI.

\subsection{Instantiation}

To instant the $DAMI(3,3)$, we can select different $SpaceKernel$ and $ChannelKernel$ with different combinations of points. To make the selection clearer and more organized, we use a visualization method, as shown in Figure \ref{fig:selection33}.
\begin{figure}[h]
	\centering
	\subfigure[$L=3$]{
		\includegraphics[width=.25\linewidth]{figures/selection_N3.png}}
	\subfigure[$L=4$]{\label{fig:selectionb33}
		\includegraphics[width=.22\linewidth]{figures/selection_N4.png}}
	\subfigure[$L=5$]{
		\includegraphics[width=.23\linewidth]{figures/selection_N5.png}}
	\subfigure[$L=6$]{
		\includegraphics[width=.2\linewidth]{figures/selection_N6.png}}
	\caption{Visualization of points selection methods with degree $L=3,4,5,6$.}
	\label{fig:selection33}
\end{figure}
In the figure, we use $1,2,3...$ to represent different points' index and each triangle is a \textit{covariant primitive}. $L$ is the degree of the DAMI, which represents how many points has participated in the construction in total.
Then the $primitives$ constructed with these points can be represented with $D_S(i, j, k)$, $D_C(l, m, n)$,..., and  so on.
Each combination of these triangles (with specified order) would be a $ShapeKernel$ or $ChannelKernel$ of $DAMI(3,3)$. There are limited numbers of triangles in the figure.
And all combinations of $ShapeKernels$ and $ChannelKernels$ in each subfigure make up a Complete Set under the specific degree $L$ and specific order $K$. By the way, Complete Set with lower degree is included in that with higher degree.
What we need to do is just specifying the $i,j,k$ and $l,m,n$ as the different indexes of points.

Here, we only consider $SpaceKernels$ and $ChannelKernels$ with order and degree both not greater than 4, which are corresponding to Figure \ref{fig:selectionb33}. There are four triangles in total, $(1,2,3)$, $(1,2,4)$, $(1,3,4)$ and $(2,3,4)$. We use them to construct the $covariants$ and specify different orders which not bigger than 4. Then we can get different $DAMI(3,3)$ with different triangle combinations and different parameters $(P,Q,O_S,O_C)$. All of the combinations are listed in Table \ref{tab:com33}. It should be illuminated that some other combinations with other point indexes  does not appear in the Table, because those combinations have the same topological structure with one in the table. One example is $D_S(2,3,4)D_C(2,3,4)$, which does not appear in table, but has the same topological structure with the first element in table, $D_S(1,2,3)  D_C(1,2,3)$. These two will get the same invariant. So just one of them is listed in table. Finally, we get a Complete Set of all combinations with both degree and order not greater than 4. Each one has a specific $ID$ in Table \ref{tab:com33}. And $G$ in the first column of the table is used just to make our enumeration more convenient and orderliness.
\begin{table}[!t]
	\centering

	\footnotesize
	\caption{All combinations of degree and order not greater than 4, corresponding to Figure \ref{fig:selectionb33}.}
	\label{tab:com33}
	\begin{tabular}{c|c|c|c}
		\toprule
		$G$ & $ID$ &Combination& $(P,Q,O_S,O_C)$ \\
		\midrule
		1&1&$D_S(1,2,3)  D_C(1,2,3)  $  &(3,3,1,1)  \\

		1&2&$D_S(1,2,3)^2D_C(1,2,3)  $  &(3,3,2,1) \\

		1&3&$D_S(1,2,3)  D_C(1,2,3)^2$  &(3,3,1,2)  \\
		1&4&$D_S(1,2,3)^3D_C(1,2,3)  $  &(3,3,3,1) \\
		1&5&$D_S(1,2,3)  D_C(1,2,3)^3$  &(3,3,1,3)  \\
		1&6&$D_S(1,2,3)^2D_C(1,2,3)^2$  &(3,3,2,2)  \\
		\midrule
		\midrule
		2&7&$D_S(1,2,3)  D_C(1,2,4)  $  &(3,3,1,1)  \\

		2&8&$D_S(1,2,3)^2D_C(1,2,4)  $  &(3,3,2,1)  \\

		2&9&$D_S(1,2,3)  D_C(1,2,4)^2$  &(3,3,1,2)  \\
		2&10&$D_S(1,2,3)^3D_C(1,2,4)  $ &(3,3,3,1)  \\
		2&11&$D_S(1,2,3)  D_C(1,2,4)^3$ &(3,3,1,3)   \\
		2&12&$D_S(1,2,3)^2D_C(1,2,4)^2$ &(3,3,2,2)  \\
		\midrule

		\midrule

		3&13&$D_S(1,2,3)  D_S(1,2,4)D_C(1,2,3)  $ &(4,3,2,1) \\
		3&14&$D_S(1,2,3)^2D_S(1,2,4)D_C(1,2,3)  $ &(4,3,3,1) \\
		3&15&$D_S(1,2,3)D_S(1,2,4)^2D_C(1,2,3)  $ &(4,3,3,1) \\
		3&16&$D_S(1,2,3)  D_S(1,2,4)D_C(1,2,3)^2$ &(4,3,2,2) \\

		3&17&$D_S(1,2,3)  D_C(1,2,3)  D_C(1,2,4)$ &(3,4,1,2) \\
		3&18&$D_S(1,2,3)  D_C(1,2,3)^2D_C(1,2,4)$ &(3,4,1,3) \\
		3&19&$D_S(1,2,3)  D_C(1,2,3)D_C(1,2,4)^2$ &(3,4,1,3) \\
		3&20&$D_S(1,2,3)^2D_C(1,2,3)  D_C(1,2,4)$ &(3,4,2,2) \\

		\midrule

		\midrule

		4&21&$D_S(1,2,3)  D_S(1,3,4)D_C(1,2,4)  $ &(4,3,2,1) \\
		4&22&$D_S(1,2,3)^2D_S(1,3,4)D_C(1,2,4)  $ &(4,3,3,1) \\
		4&23&$D_S(1,2,3)  D_S(1,3,4)D_C(1,2,4)^2$ &(4,3,2,2) \\

		4&24&$D_S(1,2,3)  D_C(1,3,4)  D_C(1,2,4)$ &(3,4,1,2) \\
		4&25&$D_S(1,2,3)  D_C(1,3,4)^2D_C(1,2,4)$ &(3,4,1,3) \\
		4&26&$D_S(1,2,3)^2D_C(1,3,4)  D_C(1,2,4)$ &(3,4,2,2) \\

		\midrule
		\midrule
		5&27&$D_S(1,2,3)  D_S(1,2,4)  D_C(1,3,4)  D_C(2,3,4)$ &(4,4,2,2) \\

		5&28&$D_S(1,2,3)  D_S(1,2,4)  D_S(1,3,4)  D_C(2,3,4)$ &(4,4,3,1) \\
		5&29&$D_S(1,2,3)  D_C(1,2,4)  D_C(1,3,4)  D_C(2,3,4)$ &(4,4,1,3) \\
		\bottomrule
	\end{tabular}
\end{table}

\subsection{Calculation}
To calculate the $DAMI(3,3)$ in practice, we should expand the $DAMI(3,3)$ to the form of \textit{Moment}.
The basic geometric moment in 2D shape space is listed in Equation (\ref{equ:moment}).
\begin{equation} \label{equ:moment}
M_{pq}=\iint x^py^q f(x,y) dxdy
\end{equation}
And the central moment is replacing the origin with the centroid, as Equation (\ref{equ:mu}):
\begin{equation}\label{equ:mu}
U_{pq}=\iint (x-\overline{x})^p (y-\overline{y})^q f(x,y) dxdy
\end{equation}
where $\overline{x}$ and $\overline{y}$ represent the center of the image in corresponding axis.

In 3D shape-color space, we define the $M(3,3)$ as Equation (\ref{equ:scm}):
\begin{equation}\label{equ:scm}
M(3,3)_{pqruvw}= \iiint x^py^qz^r R(x,y,z)^u G(x,y,z)^v B(x,y,z)^w dxdydz
\end{equation}
And the central moment is listed in Equation (\ref{equ:scu}):
\begin{align}\label{equ:scu}
U(3,3)_{pqruvw}= \iiint &(x-\overline{x})^p(y-\overline{y})^q(z-\overline{z})^r \cdot \notag \\
& (R - \overline{R})^u (G-\overline{G})^v (B-\overline{B})^w dxdydz
\end{align}
We can expand all the $DAMI(3,3)$ in the form of $U(3,3)$. By the way, some $DAMI(3,3)$ may expand to 0, such as 2,3,8,9,13,17,21,24, in Table \ref{tab:com33}.  Several examples of expanded $DAMI(3,3)$ are listed in Equation (\ref{equ:example1}-\ref{equ:example2}). We use $ID$ listed in Table \ref{tab:com33} to index the 29 $DAMI(3,3)$.
\begin{align}\label{equ:example1}
DA&MI(3,3)(1)\_numerator = \notag \\
&6*u_{001001}* u_{010010}* u_{100100} - 6* u_{001001}* u_{010100}* u_{100010} - \notag \\
&6* u_{001010}*  u_{010001}* u_{100100} + 6* u_{001010}* u_{010100}* u_{100001} + \notag \\
&6* u_{001100}* u_{010001}* u_{100010} - 6* u_{001100}* u_{010010}* u_{100001} \notag \\
DA&MI(3,3)(1)\_denominator = \notag \\
& (- 6* u_{000200}* {u_{000011}}^2 + 12* u_{000011}* u_{000101}* u_{000110} - \notag \\
&6* u_{000020}* {u_{000101}}^2 - 6* u_{000002}* {u_{000110}}^2 + \notag \\
&6* u_{000002}* u_{000020}* u_{000200})^{1/2} * u_{000000}^{5/2} \notag \\
DA&MI(3,3)(1) =\frac{DAMI(3,3)(1)\_numerator}{DAMI(3,3)(1)\_denominator}
\end{align}
\begin{align}\label{equ:example3}
DA&MI(3,3)(7)\_numerator= \notag \\
&2* u_{000001}* u_{001000}* u_{010010}* u_{100100} - 2* u_{000001}* u_{001000}* \notag \\
&u_{010100}* u_{100010} - 2* u_{000001}* u_{001010}* u_{010000}* u_{100100} + \notag \\
&2* u_{000001}* u_{001010}* u_{010100}* u_{100000} + 2* u_{000001}* u_{001100}* \notag \\
&u_{010000}* u_{100010} - 2* u_{000001}* u_{001100}* u_{010010}* u_{100000} - \notag \\
&2* u_{000010}* u_{001000}* u_{010001}* u_{100100} + 2* u_{000010}* u_{001000}* \notag \\
&u_{010100}* u_{100001} + 2* u_{000010}* u_{001001}* u_{010000}* u_{100100}- \notag \\
&2* u_{000010}* u_{001001}* u_{010100}* u_{100000} - 2* u_{000010}* u_{001100}* \notag \\
&u_{010000}* u_{100001} + 2* u_{000010}* u_{001100}* u_{010001}* u_{100000} + \notag \\
& 2* u_{000100}* u_{001000}* u_{010001}* u_{100010} - 2* u_{000100}* u_{001000}* \notag \\
&u_{010010}* u_{100001} - 2* u_{000100}* u_{001001}* u_{010000}* u_{100010} + \notag \\
& 2* u_{000100}* u_{001001}* u_{010010}* u_{100000} + 2* u_{000100}* u_{001010}* \notag \\
& u_{010000}* u_{100001} - 2* u_{000100}* u_{001010}* u_{010001}* u_{100000} \notag \\
DA&MI(3,3)(7)\_denominator = \notag \\
& (- 6* u_{000200}* {u_{000011}}^2 + 12* u_{000011}* u_{000101}* u_{000110} - \notag \\
&6* u_{000020}* {u_{000101}}^2 - 6* u_{000002}* {u_{000110}}^2 + \notag \\
&6* u_{000002}* u_{000020}* u_{000200})^{1/2} * u_{000000}^{5/2} \notag \\
DA&MI(3,3)(7) =\frac{DAMI(3,3)(7)\_numerator}{DAMI(3,3)(7)\_denominator}
\end{align}

\begin{equation}\label{equ:example2}
DAMI(3,3)(2,3,8,9,13,17,21,24)=0
\end{equation}

Considering the length, the others are not listed here.

\section{Experiments for $DAMI(3,3)$}
As listed in Table \ref{tab:com33}, except the 8 DAMIs that are expanded to 0, we get 21 usable DAMIs. In order to give a comprehensive analysis of these DAMIs, we conduct several experiments in this section.
First, we select a typical object and design an experiment to evaluate the characteristics of DAMI for different types of transformations in Section 5.1.
Then, we conduct experiments on some special objects with 1 or more symmetric planes in Section 5.2, which may affect the stability of some DAMIs. Through these cases, we can get the characteristics of each DAMI.
Based on this analysis, we give some guidelines for selecting a high-order invariant kernels.
Finally, we conduct a classification experiments for the dataset with dual-affine transformations in Section 5.3.

\subsection{Invariance to different transformations}
To evaluate the invariance of the $DAMI(3,3)$, we use one object from 3D-MNIST dataset\footnote{\url{https://www.kaggle.com/daavoo/3d-mnist}}, which is in 3D points cloud format. Then we assign color for every pixels based the $y$ coordination and generate one 3D color object, as shown in Figure \ref{fig:3dmnist}(a). We transform it with some shape affine transformations, including translation, rotation, scaling, and affine, channel affine transformations and dual-affine transformations. Examples of the 3D color object and its transformed versions are shown in Figure \ref{fig:3dmnist}.

\begin{figure}[h]
	\centering
	\subfigure[original]{
		\includegraphics[height=0.11\textheight]{figures/pc0000_1.png}}
	\subfigure[translate]{
		\includegraphics[height=0.11\textheight]{figures/pc0000_3.png}}
	\subfigure[rotation]{
		\includegraphics[height=0.11\textheight]{figures/pc0000_5.png}}
	\subfigure[scaling]{
		\includegraphics[height=0.1\textheight]{figures/pc0000_7.png}}

	\subfigure[affine]{
		\includegraphics[height=0.12\textheight]{figures/pc0000_10.png}}
	\subfigure[color affine]{
		\includegraphics[height=0.12\textheight]{figures/pc0000_c2.png}}
	\subfigure[dual-Affine]{
		\includegraphics[height=0.12\textheight]{figures/pc0000_sc3.png}}
	\caption{A 3D color object of point cloud and its  transformed versions.}
	\label{fig:3dmnist}
\end{figure}

We calculate the 21 $DAMI(3,3)$ for these objects.
To measure the stability of DAMIs, we calculate the coefficient of variation (CV), as Equation (\ref{equ:cv}), for every DAMI to each type of the transformation. The results are listed in Table \ref{tab:cv}. The smaller of the CV is, the more robust of the invariant is. And 0 is the best.

\begin{equation}\label{equ:cv}
c_v = \frac{\sigma}{\mu}
\end{equation}

\begin{table}[h]
	\tiny
	\centering
	\caption{The CVs of invariants for each type of transformation.}
	\label{tab:cv}
	\begin{tabular}{c|rrrrrrr}
		\toprule
		DAMI 		&   translation    &   rotation   &   scaling   &   affine    &   color-affine    &   DAT & ALL      \\
		\midrule
		1  &  30.4598 &  1.9720 &  68.2011 &  1.9537 &  1.0516 &  2.4390 &   2.9668  \\
		4  &  49.0937 &  1.7006 &   4.1107 &  2.3175 &  0.5987 &  0.9042 &   5.4031  \\
		5  &   0.3933 &  1.5842 &   1.7944 &  2.0706 &  1.8723 &  1.0928 &   3.2553  \\
		6  &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.0130} &  \textbf{0.0130} &   \textbf{0.0084}  \\
		7  & 647.1033 &  1.2743 &   1.6939 &  1.9432 &  4.7739 &  2.0432 &   3.5687  \\
		10 &   1.0987 &  1.8283 &   2.0222 &  2.5839 &  2.8689 &  0.9317 &   2.9907  \\
		11 &  97.7044 &  1.3899 &   1.5922 &  1.9368 &  0.0407 &  0.9426 &   3.3357  \\
		12 &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.0082} &  \textbf{0.0082} &   \textbf{0.0052}  \\
		14 & 525.2069 &  1.2780 &   1.7163 &  0.7668 &  0.0169 &  0.9621 &   1.5313  \\
		15 &   8.1589 &  1.3549 &   2.0697 &  2.5449 &  0.6408 &  0.8617 &  16.5947  \\
		16 &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.0765} &  \textbf{0.0765} &   \textbf{0.0550}  \\
		18 &   1.0154 &  2.4479 &   1.3882 &  2.0443 &  1.1687 &  1.1411 &   3.8883  \\
		19 &   0.3031 &  1.6126 &   1.8729 &  2.0763 &  1.8352 &  1.0757 &   3.0969  \\
		20 &   \textbf{0.0000} & \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  1.5923 &  1.5923 &   4.8507  \\
		22 &   2.3639 &  1.8489 &   0.8868 &  2.5441 &  2.4992 &  1.3321 &   2.8626  \\
		23 &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.1781} &  \textbf{0.1781} &   \textbf{0.1076}  \\
		25 &   0.2547 &  1.7778 &  12.4229 &  2.0604 &  1.5818 &  2.1393 &   5.5918  \\
		26 &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.0079} &  \textbf{0.0079} &   \textbf{0.0064}  \\
		27 &   \textbf{0.0000} &  \textbf{0.0000} &   \textbf{0.0000} &  \textbf{0.0000} &  \textbf{0.0516} &  \textbf{0.0516} &   \textbf{0.0316} \\
		28 &   3.3809 &  1.8488 &   0.5721 &  1.1201 &  1.5067 &  0.8754 &   5.1531  \\
		29 &   1.5408 &  1.2632 &   1.2428 &  2.1340 &  5.3576 &  1.2356 &   3.0848  \\
		\bottomrule
	\end{tabular}
\end{table}

From the results, we can see that not all $DAMI(3,3)$ are robust to all transformations. In general, invariants with $cv$ less than 1 are considered unstable. We can see that for all transformations of the object, there are 6 DAMIs robust, 6, 12, 16, 23, 26, 27. And No. 20 DAMI is only robust to shape affine transformations.
We will conduct an in-depth analysis of the causes of the errors through experiments in the next subsection.

\subsection{Guidelines for kernel selection}
We analyzed the relationship between the $ShapeKernel$ and symmetric planes of the objects in shape space. In color space, the relationship is the same.
However, the symmetry in color space rarely occurs on objects.
Therefore, some invariants, with odd-order for some color kernels, are still robust to most objects.
Another phenomenon we can see from Table \ref{tab:cv} is that the 20th DAMI is not robust to color affine invariants. Its $ColorKernel$ is characterized by the fact that points (1, 2, 3) participate in both the shape $primitive$ and the color $primitive$ construction, and the sum of their orders is odd.

From the above experiments, we can conclude that the calculation of the central moment produces a tiny error for the central moment of all orders due to the tiny error in the calculation of the center of the object.
When the object's symmetric plane is exactly the coordinate plane, the odd moments will be symmetry too and the values of the moments will be close to 0. The absolute values of the error and the moment are on the same order of magnitude and this error will be manifested.
In other cases, the magnitude of the error and the absolute value of the central moment differ greatly in magnitude, the error will be recessive, and the value of the central moment will not be affected.
For some DAMIs, although the expression is not 0, the sum of the addend of the invariant is also close to 0, which will also causes this error to be manifested.
In fact, the necessary and sufficient condition for the invariant to be stable for all objects is that when the invariant is finally expanded, there is one addend that all the central moments in it are even-order.

Considering the complexity of the actual application, we summarize the DAMI that can be used in different objects.
First, when the objects in the dataset has no more than one symmetric plane, the available invariants are (6, 12, 16, 23, 26, 27). Second, when the objects in the dataset has three or more symmetrical planes, the available DAMIs are (6, 12, 26). Third, for very few cases where colors of the object are also symmetrical in color space, 26 may be unstable, and the available invariants are (6, 12). Fourth, for the case of objects with a single color, the numerator and denominator of DAMI are all 0. The DAMIs are not available. We can use the AMI only in shape space instead.

In order to be able to guide the derivation of higher-order stable invariants, we present some relationship between  invariant kernels and invariant stability: First, A sufficient condition for the stability of DAMI is that the order of each shape kernel and color kernel is even. Second, a necessary condition for DAMI stability is that each point participating in the shape kernel and color kernel has an even number of participations. With these two guideline, we can derive DAMI in high-oder more easily.

\subsection{Classification Experiment}

In order to evaluate the discriminability of invariants, we conduct classification experiments. We select 10 objects from 3D-MNIST dataset, which are shown in Figure \ref{fig:class}.
Then we conduct 13 DAT on each of them. Finally, we get a dataset with 140 objects of 10 classes.
With the 6 $DAMI(3,3)$, (6, 12, 16, 23, 26, 27), we conduct 10-folds cross validation on the dataset with \textit{K-NearestNeighbor} (KNN). The average accuracy is 100\%. We can conclude that these six invariants have good discriminability.

\begin{figure}[h]
	\centering
	\subfigure[0]{
		\includegraphics[height=.1\textheight]{figures/pc10_1_2.png}}
	\subfigure[1]{
		\includegraphics[height=.1\textheight]{figures/pc10_2_2.png}}
	\subfigure[2]{
		\includegraphics[height=.1\textheight]{figures/pc10_3_2.png}}
	\subfigure[3]{
		\includegraphics[height=.1\textheight]{figures/pc10_4_2.png}}
	\subfigure[4]{
		\includegraphics[height=.1\textheight]{figures/pc10_5_2.png}}

	\subfigure[5]{
		\includegraphics[height=0.09\textheight]{figures/pc10_6_2.png}}
	\subfigure[6]{
		\includegraphics[height=0.09\textheight]{figures/pc10_7_2.png}}
	\subfigure[7]{
		\includegraphics[height=0.09\textheight]{figures/pc10_8_2.png}}
	\subfigure[8]{
		\includegraphics[height=0.09\textheight]{figures/pc10_9_2.png}}
	\subfigure[9]{
		\includegraphics[height=0.09\textheight]{figures/pc10_10_2.png}}
	\caption{10 objects from 3D-MNIST dataset.}
	\label{fig:class}
\end{figure}

\section{Conclusion}

In this paper, we propose a general framework to derive Moment Invariants under Dual-Affine Transformation for objects in $M$-dimensional space with $N$ channels. With this framework, we can not only  unify the invariants of objects with different dimensions into one framework, but also get the DAMI of objects with more complex formats.
This framework integrate the DAMI of different data format into one framework, such as color image in 2D space with 3 channels, 2D vector field in 2D space with 2 channel, 3D color object in 3D shape with 3 channels and objects in higher dimensional space with multi-channels.
This framework is consist with SCAMI in color images in \cite{gong2017naturally}, and 2D vector field \cite{kostkova2018affine}.
We use the 3D color objects as an example and give out a complete set of $DAMI(3,3)$ with both order and degree not greater than 4.
Then, we analyzed the root cause of the error in the moment invariant. At the same time, the invariant stability and the choice of the core in the construction process are given. This can play a very good role in the construction of higher moments.
Through invariance analysis, we select six robust DAMIs, which are robust for the objects with no more than 1 symmetric plane. And we evaluated the discriminable of them with classification experiments. The results shown that the DAMI are robust for objects under DAT. We can use this framework to derive DAMI in higher-dimensional space with multiple channels.

\section*{Acknowledgments}
This work is supported by National Key R\&D Program of China (No. 2017YFB1002703),
National Key Basic Research Program of China (No. 2015CB554507), and National Natural Science Foundation of China (No. 61379082). This work is partially supported by China Scholarship Council (CSC).

\bibliographystyle{model2-names}