---
abstract: 'Following the presentation and proof of the hypothesis that image features are particularly perceived at points where the Fourier components are maximally in phase, the concept of phase congruency (PC) is introduced. Subsequently, a two-dimensional multi-scale phase congruency (2D-MSPC) is developed, which has been an important tool for detecting and evaluation of image features. However, the 2D-MSPC requires many parameters to be appropriately tuned for optimal image features detection. In this paper, we defined a criterion for parameter optimization of the 2D-MSPC, which is a function of its maximum and minimum moments. We formulated the problem in various optimal and suboptimal frameworks, and discussed the conditions and features of the suboptimal solutions. The effectiveness of the proposed method was verified through several examples, ranging from natural objects to medical images from patients with a neurological disease, multiple sclerosis.'
author:
- 'Seyed Mohammad Mahdi Alavi, and Yunyan Zhang [^1] [^2]'
title: Phase Congruency Parameter Optimization for Enhanced Detection of Image Features for both Natural and Medical Applications
---

Optimization, Phase Congruency, Image Processing, Image Analysis, Image Recognition, Feature Detection, Biomedical image processing, Magnetic Resonance Imaging, Multiple Sclerosis.

Introduction
============

Motivation and Literature Survey
--------------------------------

Image Features Detection (IFD) is an important topic in the image-processing field [@Kovesi1996]. It aims at finding image features including lines, edges, Mach bands, corners, and blobs, by using quantitative methods. Many IFD methodologies have been proposed, which can be classified into two broad categories. First, the IFD methods that are based on dimensional metrics. The majority of the published IFD methods lie into this category including dilation-erosion residue operator [@Noble1989], gradient [@Canny1983], weak membrane [@Blake1987], anisotropic diffusion [@Perona1990], and univalue segment assimilating nucleus [@Smith1997] -based methods. The main issue of the dimensional IFD methods is that they are too sensitive to image contrasts and spatial magnifications [@Kovesi1996]. Second, the IFD methods that are based on non-dimensional metrics. Phase congruency (PC) is such a method [@Morrone1986]. References [@Kovesi1996], [@Ziou1998], [@Patel2014], and [@Noble1989] review many of the proposed IFD techniques in details.

The PC-IFD method is important for several reasons [@Kovesi1996; @Kovesi1999; @Kovesi2003]: 1) It is invariant to the image contrast, because it is not based on intensity gradient; 2) No assumption is made in the PC-IFD formulation and computation; 3) It can detect various features; and 4) It uses the image phase information. The latter might be one of the most significant features of the PC-IFD, because phase is shown to be more informative than magnitude in image processing [@Oppenheim1981]. PC is developed based on the hypothesis that image features are optimally perceived at points where the Fourier series components are maximally in phase, meaning that phases of Fourier series€™ components are similar. This hypothesis was first introduced and verified for Mach bands in [@Morrone1986]. Subsequently, its effectiveness for detecting other features types has been verified in [@Morrone1986; @Venkatesh1990; @Kovesi1996; @Kovesi1999; @Kovesi2003].

It is shown in [@Venkatesh1990] and [@Morrone1987] that the maxima of the local energy occur at the maxima of PC and vice versa. Thus, in practice, the PC measure is often obtained by computing the normalized local energy, [@Morrone1987; @Venkatesh1990; @Kovesi1999; @Kovesi2003]. Moreover, it has been illustrated that the peaks of PC outcome are higher and more distinct when local energy is computed using window-based approaches than the peaks obtained from local energy over the whole signal, [@Kovesi1996]. Consequently, Peter Kovesi proposed a multi-scale PC measurement method based on Gabor wavelets and extended it to solve two dimensions (2D), [@Kovesi1996; @Kovesi1999; @Kovesi2003]. Kovesi further modified the PC formulation to overcome noise and ill-conditioning issues, and introduced a weighting function to penalize PC measures at locations where the spread of frequencies is narrow. The proposed 2D multi-scale PC (2D-MSPC) method has been the basis of many IFD studies in various fields such as history, media, basic sciences, and medicine [@Cao2006; @Struc2009; @Zhang2012; @Obara2012; @Mouats2015; @Rijal2015; @Ziaei2014].

The current computation of 2D-MSPC requires several parameters to be tuned [*a priori*]{}. In [@Kovesi1996; @Kovesi1999; @Kovesi2003], several empirical hints have been provided for parameter tuning. But, the selection of parameters either is not explicitly discussed as in [@Cao2006; @Struc2009; @Zhang2012], or is performed manually based on trail and error [@Obara2012; @Mouats2015; @Rijal2015; @Ziaei2014]. In [@Obara2012] and [@Mouats2015], the parameters were manually optimized to increase the visualization of image features. In [@Rijal2015] the parameters were changed by trial and error, and the study found that an improvement in IFD occurred with parameters that maximized the signal-to-noise ratio of the image. In [@Ziaei2014], a set of experiments was also done based on the trial and error, to determine the best, fixed values for computing the maximum moment of PC covariance. The authors reported that there was no clear rule on how the 2D-MSPC parameters should be tuned. Collectively, there is a critical need for a standard parameter tuning method for 2D-MSPC. Finding a solution for this problem can enhance our understanding of the concept of PC and thereby promoting its applications.

Contribution of This Paper
--------------------------

We have classified the 2D-MSPC optimization problem as follows:

-   What optimization criteria should be used for the tuning of 2D-MSPC parameters?

-   How can we formulate the problem mathematically and compute the parameters of 2D-SMPC optimally and automatically?

This paper aims to address these questions. The 2D-MSPC, its definition, computational method and a list of its fundamental parameters are briefly described in section \[sec:2dmspc\]. The importance of parameter setting in 2D-MSPC IFD, which leads to Q1 and Q2, is discussed in section \[sec:pardisc\]. In section \[sec:opt\], we define a criterion for parameter optimization, based on the maximum and minimum moments of the 2D-MSPC. We formulate the problem in various optimal frameworks, and then describe about a suboptimal solution accordingly. Finally, the effectiveness of the proposed method is verified through several examples from both natural and medical images in section \[sec:ex\]. This includes two examples from magnetic resonance imaging (MRI) of patients with multiple sclerosis (MS). We illustrate that the proposed 2D-MSPC optimization is useful to MS lesions detection.

Nomenclatures
-------------

In this paper the following general notations are used. ${\ensuremath{\mathbb{N}}}$, ${\ensuremath{\mathbb{R}}}$ and ${\ensuremath{\mathbb{C}}}$ denote the integer, real and complex domains, respectively. $\|A\|$ denotes the norm of matrix $A$, where no subscript means that any norm can be used. $\|A\|_F$, $\|A\|_p$, $\|A\|_\infty$, $\|A\|_2$ and $\|A\|_1$ represent Frobenius-norm ($F-$norm), Schatten $p-$norm, $\infty-$norm, 2-norm and 1-norm of matrix $A$. The consistent norm of matrix $A$ is denoted by $\|A\|_c$. The determinant of matrix $A$ is denoted by $\mbox{det}(A)$. The optimal value of the parameter $a$ is denoted by $a^*$ (a superscript asterisk). It should not be confused with the convolution operator $*$, which is used in this paper as $A*B$. Note that there are several other specific notations also defined in the text, which will be introduced when they appear.

Two-Dimensional Multi-Scale Phase Congruency {#sec:2dmspc}
============================================

In this section, fundamental concepts of the 2D-MSPC are briefly described. Interested readers are directed to [@Kovesi1996], [@Kovesi1999] and [@Kovesi2003] for details. The 2D-MSPC is a combination of one-dimensional PC calculated over several orientations. The 2D-MSPC based on the measurement of local energy is given by $$\begin{aligned}
\label{eq:pc}
& \mbox{PC}(x)= \frac{ \sum_{o}  W_o(x) \left \lfloor  E_o(x)-T_o \right \rfloor}{ \sum_{o} \sum_{n}A_{no}(x) + \varepsilon}
\\ \nonumber & o=1,2,\ldots, O
\\ \nonumber & n=1,2,\ldots, N_o\end{aligned}$$ where, $o$ and $n$ denote the orientation and scale indexes, respectively. $O$ is the total number of orientations. $N_o$ is the total number of scales at the $o-$th orientation. $\lfloor \rfloor$ is the floor operator. Note that $\lfloor y \rfloor = y$ if $y>0$, otherwise zero. $E_o$ is the local energy at the $o-$th orientation, which is calculated by using the Hilbert transform as follows $$\begin{aligned}
\label{eq:e}
E_o(x)=\sqrt{F_o^2(x)+H_o^2(x)}\end{aligned}$$ where, $F_o(x)$ and $H_o(x)$ are the AC component and Hilbert transform of the image signal $I(x)$ at the $o-$th orientation. However, because the Hilbert transform operator is an improper integral, $F_o(x)$ and $H_o(x)$ are computed by convolving the image signal with a pair of even and odd wavelets filters in quadrature as follows [@Morrone1987; @Venkatesh1990]: $$\begin{aligned}
\label{eq:fs}
&F_o(x)= \sum_{n} I(x) * M_{no}^{e}\\
\label{eq:hs}
&H_o(s)= \sum_{n} I(x) * M_{no}^{o}$$ where, $*$ denotes the convolution operator. $M_{no}^e$ and $M_{no}^o$ are even and odd filters at the $o-$th orientation and $n-$th scale, generated by the logarithmic Gabor function $$\begin{aligned}
\label{eq:lgabor1}
G_{no}(f)=\exp \left(\frac{-\left(\log\frac{f}{\hat{f}_{no}}\right)^2}{2\left( \log \sigma_{no}\right)^2}\right)$$ where, $\hat{f}_{no}$ and $\sigma_{no}$ denote the centre frequency and bandwidth of the $n-$th Gabor filer at the $o-$th iteration. The filters are scaled as follows $$\begin{aligned}
\label{eq:lgabor2}
\hat{f}_{no}=\frac{1}{\lambda_{\min o} \times \eta_o^{(n-1)}}$$ where, $\lambda_{\min o}$ is the minimum wavelength (maximum centre frequency) of wavelets in banks of $\{M_{1o}^e, \ldots, M_{No}^e\}$ and $\{M_{1o}^o, \ldots, M_{No}^o\}$, and $\eta_o$ represents the distance between successive filters in the bank. It is noted that $M_{no}^e$ and $M_{no}^o$ are identical with a $90^{\circ}$ shift. Figure \[fig:gf\] shows the spectra of 4-scale logarithmic Gabor wavelets at the $o-$th orientation with $\lambda_{\min o}=3$, $\eta_o=3$ and $\sigma_{no}=0.55$, $n=1,\ldots,4$. Note that in the logarithmic frequency scale, the spectra of all Gabor functions are identical.

![Spectra of 4-scale Gabor wavelets.[]{data-label="fig:gf"}](gff.png)

Subsequently, the normalization factor in the 2D-MSPC is computed in the wavelet’s framework as follows $$\begin{aligned}
\label{eq:nf}
\sum_{n}A_{no}=\sum_{n} \sqrt{  (I(x) * M_{no}^{e})^2+(I(x) * M_{no}^{o})^2}.\end{aligned}$$

     Parameter             Tuning hints
  -- --------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     $c_o$                 [Between 0 and 1, typical 0.4 or 0.55]{}
     $g_o$                 [Typical value, 10]{}
     $ \lambda_{\min o}$   [The smallest value is the Nyquist wavelength of 2 pixels. Because of aliasing 3 pixels or above is suggested.]{}
     $\sigma_{no}$         [The smaller $\sigma_{no}$, the larger the bandwidth of the filter. The following sets are suggested: $[\sigma_{no}=0.85,\eta_o=1.3]$, or ]{}
     $\eta_o$              [$[\sigma_{no}=0.75,\eta_o=1.6]$ (Bandwidth $\approx$ 1 octave), or $[\sigma_{no}=0.65,\eta_o=2.1]$, or $[\sigma_{no}=0.55,\eta_o=3]$ (Bandwidth $\approx$ 2 octaves)]{}
     $k_o$                 [Typical value 2, up to 10 or 20 for noisy images]{}
     $\varepsilon$         [Small scalar, typically 0.0001 ]{}
     $N_o$                 [Try values 3 to 6 ]{}
     $O$                   [A filter spacing of $30^{\circ}$ has found good. ]{}

\[table:par\]

As discussed in [@Kovesi1996], the calculation of PC makes sense only in locations where the spread of frequencies is significant. In 2D-MSPC formulation, $W_o(x)$ is a weighting function, which penalizes the information at locations where the spread of frequencies is narrow at the $o-$th orientation. It is defined by a sigmoid function $$\begin{aligned}
W_o(x)=\frac{1}{1+\exp \left(g_o(c_o-s_o(x))\right)}\end{aligned}$$ where, $c_o$ denotes the cutoff (mid-point) point of the sigmoid function; $g_o$ is a gain that controls the rate of weighting; $s_o$ is a measure of frequency spread, which ranges between 0 and 1, and is given by $$\begin{aligned}
s_o(x)=\frac{1}{N_o}\frac{\textstyle \sum_{n}A_{no}(x)}{A_{\max}(x)+\varepsilon}.\end{aligned}$$ where, $\textstyle \sum_{n}A_{no}(x)$ is computed through and $A_{\max}(x)$ is the amplitude of the filter pair having maximum response at $x$.

The parameter $T_o$ is an estimation of noise power that is subtracted from the energy of the signal. It is computed based on the assumptions that the image noise is additive, that the noise power spectrum is constant, and that the image features occur at isolated locations: $$\begin{aligned}
T_o=\mu_{Ro}+k_o \sigma_{Ro}\end{aligned}$$ where, $\mu_{Ro}$ and $\sigma_{Ro}^2$ are the mean and variance of the Rayleigh that describes the noise energy response; $k_o$ is a scaling factor used to estimate the maximum degree of the noise response.

When the spread of frequencies is narrow, $E_o$ and $\sum_n A_{no}$ become very small, making the computations to become ill-conditioned. The parameter $\varepsilon$ is to address this issues.

Table \[table:par\] summarizes the parameters discussed above which are fundamental and need to be fine-tuned in the 2D-MSPC calculations, including some hints from [@Kovesi1996; @Kovesi1999; @Kovesi2003; @Kovesi_url], mainly for manual tuning.

Hereafter, we denote all parameters in a vector format, $$\begin{aligned}
\label{parvec}
\upsilon_o= \left[ \begin{array}{ccccccccc}
    c_o &  g_o & \lambda_{\min o} & \sigma_{no} & \eta_o & k_o & \varepsilon & N_o & O \\ 
  \end{array}\right]\end{aligned}$$ where, $\upsilon_o$ denote the parameter vector at orientation $o$.

As usual in the numerical optimization, we assume that upper and lower limits of the parameter vector are known [*a priori*]{}. The upper and lower limits of $\upsilon$ in the vector format are denoted by $\overline{\upsilon}_o$ and $\underline{\upsilon}_o$ respectively, such that $$\begin{aligned}
\underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\end{aligned}$$ In this paper, $\upsilon_o^*$ denotes the optimal parameter vector and $\mbox{PC}_o^*$ denotes the PC matrix computed with the optimal parameter vector $\upsilon_o^*$, both at orientation $o$.

Description of Q1 and Q2 {#sec:pardisc}
========================

It is illustrated that visualization of the image will significantly vary by changing only one parameter of the 2D-MSPC [@Mouats2015] and [@Rijal2015]. Figure \[fig:lena\] shows the visualization of a lena’s image by using 2D-MSPC, with different values of the cutoff point $c_o$. Based on information from Table \[table:par\], other parameters of the 2D-MSPC are fixed to: $k_o=2$; $\varepsilon=0.0001$; $N_o=4$; $O=6$; $\lambda_{\min o}=3$; $\eta_o=2.1$; and $g_o=10$, for $o=1,\dots, O$ and $n=1,\ldots,N_o$. Figure \[fig:lena-orig\] shows the lena’s grayscale image of lena that was used in this simulation. Figures \[fig:lena-co-p55\] and \[fig:lena-co-p10\] show the 2D-MSPC images for $c_o=0.55$ and $c_o=0.10$, respectively. It is seen that by decreasing the cutoff value, image features become increasingly detectable. There are several features that could not be detected with $c_o=0.55$ (Figure \[fig:lena-co-p55\]). This issue gives rise to questions Q1 and Q2.

[0.25]{} ![Visualization of a lena image by using the 2D-SMPC, with different values of the cutoff point $c_o$, while other parameters are fixed.[]{data-label="fig:lena"}](lena-orig "fig:"){width="\linewidth"}

\

[0.25]{} ![Visualization of a lena image by using the 2D-SMPC, with different values of the cutoff point $c_o$, while other parameters are fixed.[]{data-label="fig:lena"}](lena-pc-co-p55 "fig:"){width="\linewidth"}

[0.25]{} ![Visualization of a lena image by using the 2D-SMPC, with different values of the cutoff point $c_o$, while other parameters are fixed.[]{data-label="fig:lena"}](lena-pc-co-p10 "fig:"){width="\linewidth"}

In the next few sections, these questions are addressed.

Parameter Optimization {#sec:opt}
======================

In this paper, several optimal frameworks are proposed for parameter tuning for the 2D-MSPC, which are based on the maximization of the PC momentums. The maximum and minimum moments of the 2D-MSPC, $M$ and $m$ respectively, are given by [@Kovesi2003]: $$\begin{aligned}
\label{maxmo}
& M=\frac{1}{2}\left(\alpha+\gamma+\sqrt{\beta^2+(\alpha-\gamma)^2}\right)\\
\label{minmo}
& m=\frac{1}{2}\left(\alpha+\gamma-\sqrt{\beta^2+(\alpha-\gamma)^2}\right)\end{aligned}$$ with $$\begin{aligned}
\label{alpha}
& \alpha=\textstyle \sum_{o}\left(\mbox{PC}_o \cos(\theta_o)\right)^2\\
\label{gamma}
& \gamma=\textstyle \sum_{o}\left(\mbox{PC}_o \sin(\theta_o)\right)^2\\
\label{beta}
& \beta=2 \textstyle \sum_{o}\left(\mbox{PC}_o \cos(\theta_o)\right)\left(\mbox{PC}_o \sin(\theta_o)\right).\end{aligned}$$ where, $\mbox{PC}_o$ and $\theta_o$ represent the PC and axis angle at orientation $o$.

The calculation of momentums corresponds to performing a singular value decomposition to the PC covariance matrix, and thus momentums correspond to the singular values, [@Kovesi2003]. On the other hand, singular values indicate the level of increase in energy that can occur between the input and output of a given system. The larger the maximum singular value of a system, the greater the increase in energy, (see Definition 40.2 and descriptions in §40.2.1, p. 652,[@Levine1995]). Therefore, it may be deduced that momentums represent the maximum of local energies. It was also shown that the maxima of local energy are indications of features in an image, [@Venkatesh1990]. The answer to Q1 can now be simplified as such that a potential criterion for parameter tuning of 2D-MSPC should be given in terms of $M$ and $m$.

From the potential relationships between momentums, singular values, local energies, 2D-MSPC and image features, it may be deduced that the image features detection can be enhanced by increasing the maximum and minimum moments of the PC. In other words, optimal values of 2D-MSPC parameters are those that maximize $M$ and $m$. Both $M$ and $m$ are matrices and their size is equal to the size of $\mbox{PC}_o$. Thus, the optimal parameters of 2D-MSPC should be obtained by solving a matrix optimization problem (an answer to Q2).

In order to formulate the problem, we define a cost function $\mathcal{M}$:$$\begin{aligned}
\label{Mm}
\mathcal{M}=\mu_1 \times  M+\mu_2 \times  m\end{aligned}$$ where, $\mu_i \in [0,1]$, $i=1,2$. The cost function enables us to not only maximize both $M$ and $m$, but also to manage their interactions.

As mentioned above, the question of parameter tuning for 2D-MSPC is a matrix optimization problem. Matrix maximization has widely been used in optimal (experimental) designs [@Fredov1972; @Goodwin1977; @Atkinson1982; @Pronzato1985; @Berger1994; @Pukelsheim1993], where the Fisher information matrix [@Fisher1971] is maximized in order to reduce the Cram' er-Rao bound [@Cramer1946; @Rao1945]. The Fisher information matrix is maximized by maximizing some real-valued summary statistics (optimal) criteria), [@Atkinson2007]. The majority of matrix optimization methods are based on the maximization of the determinant of a matrix, referred to as D-optimal problem [@Berger1994; @StJohn1975]. The rest of the methods are mainly based on the maximization of the trace or maximization of the minimum eigenvalue of a matrix.

In the D-optimal framework, maximization of the cost function is given by $$\begin{aligned}
\begin{aligned}\label{D-opt}
& \underset{\upsilon_o}{\text{minimize}}
& & -\mbox{det}\Big(\mathcal{M}(\upsilon_o)\Big) \\
& \text{subject to}
& & \underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\\
\end{aligned}\end{aligned}$$ where $\det(\mathcal{M})$ denotes the determinant of $\mathcal{M}$. To be consistent with the optimization literature, the maximization of $z$ is written as the minimization of $-z$ throughout this paper. The solution to this optimization problem will provide the optimal parameters for 2D-MSPC.

For the PC-based image processing, we propose another optimization approach based on the norm of $\mathcal{M}$ as follows: $$\begin{aligned}
\begin{aligned}\label{norm-opt}
& \underset{\upsilon_o}{\text{minimize}}
& & -\left\| \mathcal{M}(\upsilon_o) \right\| \\
& \text{subject to}
& & \underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\\
\end{aligned}\end{aligned}$$ where $\| \mathcal{M}\|$ denotes the norm of $\mathcal{M}$. The matrix norm definition is given in Appendix \[matrixnorm\].

For $\mu_1=\mu_2=1$, $\mathcal{M}=M+  m= \textstyle \sum_o PC_o^2$, and the D-optimal and norm-optimal problems and are transferred respectively to $$\begin{aligned}
\begin{aligned}\label{D-opt2}
& \underset{\upsilon_o,~o=1,\ldots,O}{\text{minimize}}
& & -\mbox{det}\Big(\mbox{PC}_o^2(\upsilon_o)\Big) \\
& \text{subject to}
& & \underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\\
& & & \mathcal{M}=M+m\\
\end{aligned}\end{aligned}$$ and, $$\begin{aligned}
\begin{aligned}\label{norm-opt2}
& \underset{\upsilon_o,~o=1,\ldots,O}{\text{minimize}}
& & -\left\| \mbox{PC}_o^2(\upsilon_o) \right\| \\
& \text{subject to}
& & \underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\\
& & & \mathcal{M}=M+m\\
\end{aligned}\end{aligned}$$

Differences between the optimization sets { and } and { and } are the followings.

-   In order to solve or , at each iteration of the optimization, the PC matrices of all orientations and $\mathcal{M}$ must be computed. Solving and does not require $\mathcal{M}$ at each iteration of the optimization.

-   The solution to or results in one optimal parameter vector for all orientations. The solution to or results in $O$ numbers of optimal parameter vectors (one optimal vector for each orientation).

A Sub-Optimal Solution to the Norm-Optimal Problem
--------------------------------------------------

In order to avoid the computation of $\mathcal{M}$ at each iteration of the optimization, a sub-optimal criterion of is derived, which optimizes the parameter vector separately at each orientation. Note that the results of this subsection are valid if the matrix norm is consistent and $\mu_i \neq 0$ and $\mu_i \neq 1$ , $i=1,2$. Based on the literature on optimal designs [@Fredov1972; @Goodwin1977; @Atkinson1982; @Pronzato1985; @Berger1994; @Pukelsheim1993], the covariance matrix satisfies the following inequality $$\mbox{cov}(\upsilon) \geq \mathcal{F}^{-1}(\upsilon)$$ where, $\mbox{cov}$ denotes the covariance matrix, $\upsilon$ is the parameter vector and $\mathcal{F}$ denotes the Fisher information matrix. In this case, $\mathcal{F}^{-1}$ is the lower bound of the covariance of estimations (known as the Cram' er-Rao bound [@Cramer1946; @Rao1945]). It is a common approach that in order to reduce the covariance of estimations, the the Cram' er-Rao bound is reduced by maximizing the Fisher information matrix.

In this section, we obtain an upper bound of the norm of $M$ and $m$. We then derive a sub-optimal solution, which is based on the maximization of the upper bound of $M$ and $m$. First, we need the matrix norm to be consistent.

([@Lyche2012], §8.1.1, p. 178) A matrix norm is [*consistent*]{} if it is defined on ${\ensuremath{\mathbb{C}}}^{q\times r}$ for all $q,r\in {\ensuremath{\mathbb{N}}}$ and the sub-multiplicative property $$\begin{aligned}
\label{ableqab}
\| AB \| \leq \|A\| ~\|B\|\end{aligned}$$ holds for all matrices $A$ and $B$ for which the product $AB$ is defined.

In this paper, the consistent norm is denoted by a subscript, $\|.\|_c$.

([@Lyche2012], [@Schatten1960]) \[lemakf\] The Frobenius and all Schatten $p-$norms are sub-multiplicative, i.e., $$\begin{aligned}
\label{ableqab}
\| AB \|_F &\leq \|A\|_F ~\|B\|_F\\
\| AB \|_p & \leq \|A\|_p ~\|B\|_p\end{aligned}$$

\[lempc\] Let $\mbox{PC}_o$ denotes the PC at orientation $o$ and the 2D-MSPC is computed over $O$ orientations. Then, the following relationships hold between the maximum and minimum moments and $\mbox{PC}_o$. $$\begin{aligned}
M&=\frac{1}{2}\Big(\textstyle \sum_{o} \mbox{PC}_o^2 +\sqrt{\textstyle \sum_{o} \mbox{PC}_o^4 }\Big)\\
m&=\frac{1}{2}\Big(\textstyle \sum_{o} \mbox{PC}_o^2-\sqrt{\textstyle \sum_{o} \mbox{PC}_o^4 }\Big)\end{aligned}$$ [*Proof*]{}: Appendix \[p-lempc\].

\[thpcbnd\] Let $\mbox{PC}_o$ denotes the PC at orientation $o$. The 2D-MSPC is computed over $O$ orientations. If a consistent matrix norm is used, then $$\begin{aligned}
\label{Mbnd}
 &\left\| M \right\|_c  \leq  \textstyle \sum_{o} \left\| \mbox{PC}_o^2 \right\|_c\\
 \label{mbnd}
 &\left\| m \right\|_c  \leq  \textstyle \sum_{o} \left\| \mbox{PC}_o^2 \right\|_c\end{aligned}$$ where, $c$ can take $F$ or $p$. If the subscript $c$ is replaced with $F$, the $F-$norm is applied. If the subscript $c$ is replaced with $p$, the Schatten $p-$norm is applied. : Appendix \[p-pcbnd\].

\[thsubopt\] A sub-optimal solution to the optimization problem with the consistent norm is obtained by solving the following optimization problem. $$\begin{aligned}
\begin{aligned}\label{norm-subopt}
& \underset{\upsilon_o,~o=1,\ldots,O}{\text{minimize}}
& & -(|\mu_1|+|\mu_2|)\left\| \mbox{PC}_o^2(\upsilon_o) \right\|_c \\
& \text{subject to}
& & \underline{\upsilon}_o \leq \upsilon_o \leq \overline{\upsilon}_o\\
\end{aligned}\end{aligned}$$ [*Proof*]{}: Appendix \[p-suboptlem\].

The differences between and are as follows:

-   The solution to is obtained by the parameter vector optimization individually at each orientation, while the solution to is a parameter vector optimal for all orientations.

-   The optimization problem does not require the computation of $\mathcal{M}$ at each iteration of optimization. To solve , the cost function $\mathcal{M}$ is computed at each iteration of optimization.

The differences between and are as follows:

-   The norm in can be of any type, but in , it has to be consistent and satisfies the sub-multiplicative property.

-   In , $\mathcal{M}=M+m$, but $\mathcal{M}$ can be any combinations of $M$ and $m$ in .

Norm Selection
--------------

Any norm type can be used in and . However, if $\mu_i \neq 0$ and $\mu_i \neq 1$, $i=1,2$, only can be applied, not . The prerequisite to apply is that the sub-multiplicative property must hold for the norm.

There are also some other features associated with the norm definitions that might be important. If 1- or $\infty$-norms of the matrix is chosen, because they are associated with the maximum absolute column and row of the matrix, certain features, which are not along that row or column, might be overlooked. The 2-norm, Frobenius-norm ($F-$norm), or Schatten $p-$norm might be more efficient, because they are based on the eigenvalues and singular values, which correspond to the energy of a signal. The 2-norm is the maximum singular value of a matrix: $\|A\|_2=\sigma_{\max}$ [@Lyche2012]. The $F-$norm is the square root of the summation of the singular values: $\|A\|_F=(\textstyle \sum_i \sigma_i^2)^{1/2}$ [@Lyche2012], and the Schatten $p-$norm is $\|A\|_p=(\textstyle \sum_i \sigma_i^p)^{1/p}$ [@Schatten1960]. Thus, it can be deduced that the 2-norm may amplify the single most-dominant feature, while the $F-$ and $p-$ norms may strengthen features by taking their summation.

In summary, the 2D-MSPC parameters can be optimally and automatically tuned by maximizing , through or . If $\mathcal{M}=M+m$ is going to be maximized, and or can be used. If $\mu_i \neq 0$ and $\mu_i \neq 1$, $i=1,2$, can provide a sub-optimal solution to . This answers Q2.

     Orientation      $c_o^*$   $g_o^*$   Cost function opt. value
  -- ---------------- --------- --------- ------------------------------------------------------
     $o=1,\dots,6$    0.1212    45.4827   $\mbox{det}(\mathcal{M})^*=7.0088\times 10^{-55}$
     $o=1,\dots,6$    0.0100    50.0000   $\|\mathcal{M}\|_F^*=114.1408$
     $o=1$            0.1001    49.9997   $\mbox{det}(\mbox{PC}_1^2)^*=6.4589\times 10^{-241}$
     $o=2$            0.1021    39.9334   $\mbox{det}(\mbox{PC}_2^2)^*=1.7131\times 10^{-33}$
     $o=3$            0.1216    46.6111   $\mbox{det}(\mbox{PC}_3^2)^*=1.4946\times 10^{-43}$
     $o=4$            0.1946    42.9731   $\mbox{det}(\mbox{PC}_4^2)^*=3.2306\times 10^{-218}$
     $o=5$            0.1566    39.2954   $\mbox{det}(\mbox{PC}_5^2)^*=1.0452\times 10^{-37}$
     $o=6$            0.1492    46.3088   $\mbox{det}(\mbox{PC}_6^2)^*=1.7254\times 10^{-53}$
     $o=1,\ldots,6$   0.1374    44.1869   $\mbox{det}(\mathcal{M})^*=\mbox{inf}$
     $o=1$            0.0100    50.0000   $\|\mbox{PC}_1^2\|_F^*=8.3413\times 10^{3}$
     $o=2$            0.0100    50.0000   $\|\mbox{PC}_2^2\|_F^*=6.3400\times 10^{3}$
     $o=3$            0.0100    50.0000   $\|\mbox{PC}_3^2\|_F^*=3.2968\times 10^{3}$
     $o=4$            0.0100    49.9998   $\|\mbox{PC}_4^2\|_F^*=2.9181\times 10^{3}$
     $o=5$            0.0100    50.0000   $\|\mbox{PC}_5^2\|_F^*=4.7028\times 10^{3}$
     $o=6$            0.0100    50.0000   $\|\mbox{PC}_6^2\|_F^*=7.6949\times 10^{3}$
     $o=1,\ldots,6$   0.0100    50.0000   $\|\mathcal{M}\|_F^*=3.0379\times 10^{4}$

\[table:resultex1\]

Illustrative Examples {#sec:ex}
=====================

Example 1 {#sec:ex1}
---------

Consider the Lena’s image (Figure \[fig:lena-orig\]). The objective in this example is to find the optimal values for the parameters of the weighting function while keeping the other parameters constant, where $$\begin{aligned}
 & \lambda_{\min o} =3,~ \sigma_{no} =0.55,~ \eta_o=2.1,~ k_o =2,\\& \varepsilon=0.0001 ,~ N_o =4,~ O=6. \end{aligned}$$

The parameter vector in this example contains only $c_o$ and $g_o$, thus $$\begin{aligned}
\upsilon_o= [ \begin{array}{cc}
    c_o & g_o \end{array}]
    \end{aligned}$$ We set fixed upper and lower limits $\upsilon_o$ as follows: $$\begin{aligned}
    & \underline{\upsilon}_o=[ \begin{array}{cc}
    0.1 & 1 \end{array}] \\
    & \overline{\upsilon}_o=[ \begin{array}{cc}
    0.9 & 50 \end{array}].\end{aligned}$$ With $c_o < 0.1$, almost all spreads of frequencies (even narrow ones) will be kept, but that is not desirable. The change in the slope of the weighting function is also not significant for $g_o > 50$.

In order to compare , , and , this example focuses on $\mu_1=\mu_2=1$, and therefore, $$\begin{aligned}
\mathcal{M}=M+m.\end{aligned}$$

Table \[table:resultex1\] shows the optimal values of the parameter vector and cost functions. Recall that the solution to and results in one optimal vector for all orientations, and the solution to or results in one optimal vector for each orientation. The results of Frobenius-norm optimization are given in Table \[table:resultex1\]. Recall that the superscript $^*$ represents the optimal value.

Figure \[fig:lena-opt\] shows the image of $\mathcal{M}^*$ for , , and . Features such as lines, edges, Mach bands, corners, and blobs are satisfactorily detected in all images. Visually, it is seen that all four optimization methods result in the same image.

Figure \[fig:diff-2par\] shows the comparison between determinant-based optimization methods and , and norm-based optimization methods and . Differences are demonstrated in the subtraction images between the two approaches. Notably, although the determinant and norm optimal frameworks result in a similar performance visually, their subtraction images are not empty, meaning not identical.

[0.25]{} ![Detection of features from the Lena’s image. This is done by optimizing the parameters of the weighting function for 2D-MSPC in Example 1. The image of $\mathcal{M}^*$ is obtained through methods , , and .[]{data-label="fig:lena-opt"}](pc_lena_det_Mpm_6or "fig:"){width="\linewidth"}

[0.25]{} ![Detection of features from the Lena’s image. This is done by optimizing the parameters of the weighting function for 2D-MSPC in Example 1. The image of $\mathcal{M}^*$ is obtained through methods , , and .[]{data-label="fig:lena-opt"}](pc_lena_Fronorm_Mpm_6or "fig:"){width="\linewidth"}

\

[0.25]{} ![Detection of features from the Lena’s image. This is done by optimizing the parameters of the weighting function for 2D-MSPC in Example 1. The image of $\mathcal{M}^*$ is obtained through methods , , and .[]{data-label="fig:lena-opt"}](pc_lena_det_PC_6or "fig:"){width="\linewidth"}

[0.25]{} ![Detection of features from the Lena’s image. This is done by optimizing the parameters of the weighting function for 2D-MSPC in Example 1. The image of $\mathcal{M}^*$ is obtained through methods , , and .[]{data-label="fig:lena-opt"}](pc_lena_Fronorm_PC_6or "fig:"){width="\linewidth"}

[0.25]{} ![Comparison between determinant-based optimization methods and , also between norm-based optimization methods and .[]{data-label="fig:diff-2par"}](pc_lena_det_Mpm_m_PC_2par "fig:"){width="\linewidth"}

\

[0.25]{} ![Comparison between determinant-based optimization methods and , also between norm-based optimization methods and .[]{data-label="fig:diff-2par"}](pc_lena_Fnorm_Mpm_m_PC_2par "fig:"){width="\linewidth"}

For 1-, 2-, $\infty-$ and $F-$ norms, the optimization method results in the optimal parameter vector $\upsilon_o^*=[0.1 ~ 50]$ with $\|\mathcal{M}\|^*_1=219.5310$, $\|\mathcal{M}\|^*_2=75.4579$, $\|\mathcal{M}\|^*_{\infty}=204.6513$, and $\|\mathcal{M}\|^*_F=114.1408$, respectively. The image of $\mathcal{M}^*$ with optimal settings for the weighting function is shown in Figure \[fig:lena-norm-opt\], obtained through based on 1-, 2-, $\infty-$ and $F-$ norms. Again, it is seen that all four norm types result in the same visual performance. This example shows that the norm-based optimization of 2D-MSPC for the Lena’s image is robust to the norm type.

[0.25]{} ![Impact of using different norm types for the detection of features in Lena’s image. The image of $\mathcal{M}^*$ and the value of the cost function (calculated for the optimal parameter vector) are achieved from based on 1-, 2-, $\infty-$ and $F-$ norms. []{data-label="fig:lena-norm-opt"}](pc_lena_1norm_Mpm_6or "fig:"){width="\linewidth"}

[0.25]{} ![Impact of using different norm types for the detection of features in Lena’s image. The image of $\mathcal{M}^*$ and the value of the cost function (calculated for the optimal parameter vector) are achieved from based on 1-, 2-, $\infty-$ and $F-$ norms. []{data-label="fig:lena-norm-opt"}](pc_lena_2norm_Mpm_6or "fig:"){width="\linewidth"}

\

[0.25]{} ![Impact of using different norm types for the detection of features in Lena’s image. The image of $\mathcal{M}^*$ and the value of the cost function (calculated for the optimal parameter vector) are achieved from based on 1-, 2-, $\infty-$ and $F-$ norms. []{data-label="fig:lena-norm-opt"}](pc_lena_Infnorm_Mpm_6or "fig:"){width="\linewidth"}

[0.25]{} ![Impact of using different norm types for the detection of features in Lena’s image. The image of $\mathcal{M}^*$ and the value of the cost function (calculated for the optimal parameter vector) are achieved from based on 1-, 2-, $\infty-$ and $F-$ norms. []{data-label="fig:lena-norm-opt"}](pc_lena_Fronorm_Mpm_6or "fig:"){width="\linewidth"}

Because the PC is a matrix with elements between 0 and 1, the determinant is expected to be small. Note that, the larger the image, the smaller the determinant is most likely achieved. Thus, from the numerical point of view, solving such an optimization problem should be more sophisticated than solving the norm-based optimization problem.

     $c_o^*$   $g_o^*$   $\lambda_{\min o}^*$   $\sigma_{no}^*$   $\eta_o^*$   $N_o^*$   $O^*$   Cost function opt. value
  -- --------- --------- ---------------------- ----------------- ------------ --------- ------- ---------------------------------------------
     0.1       49.9994   3.055                  0.4               3.8315       4         1       $\|\mathcal{M}\|_F^*=442.0584$
     0.1       50.00     2.6894                 0.4               4.0          4         1       $\|\mbox{PC}_o^2\|_F^*=6.2033\times 10^{4}$
                                                                                                 $\|\mathcal{M}\|_F^*=1.2407\times 10^{5}$

\[table:resultex2\]

[0.25]{} ![The results of the optimization of 7 parameters by solving Frobenius norm-based methods in Example 2. Figures \[fig:pc-Fnorm-mpm-full\] and \[fig:pc-Fnorm-pc-full\] shows the images of $\mathcal{M}^*$ with optimal parameters obtained through and , respectively. Figure \[fig:diff\] shows the difference by subtracting the images. []{data-label="fig:lena-7par"}](pc_lena_Fnorm_Mpm_full "fig:"){width="\linewidth"}

[0.25]{} ![The results of the optimization of 7 parameters by solving Frobenius norm-based methods in Example 2. Figures \[fig:pc-Fnorm-mpm-full\] and \[fig:pc-Fnorm-pc-full\] shows the images of $\mathcal{M}^*$ with optimal parameters obtained through and , respectively. Figure \[fig:diff\] shows the difference by subtracting the images. []{data-label="fig:lena-7par"}](pc_lena_Fnorm_PC_full "fig:"){width="\linewidth"}

\

[0.25]{} ![The results of the optimization of 7 parameters by solving Frobenius norm-based methods in Example 2. Figures \[fig:pc-Fnorm-mpm-full\] and \[fig:pc-Fnorm-pc-full\] shows the images of $\mathcal{M}^*$ with optimal parameters obtained through and , respectively. Figure \[fig:diff\] shows the difference by subtracting the images. []{data-label="fig:lena-7par"}](pc_lena_Fnorm_Mpm_m_PC_full "fig:"){width="\linewidth"}

[0.5]{} ![image](IM-0004-0054_patient6_T2-origPaper){width="1\linewidth"}

[0.5]{} ![image](cost_IM-0004-0054_patient6_T2_Fnorm_M_full_Paper){width="1\linewidth"}

Example 2:
----------

The objective of this example is to find the optimal values for the number of scales, number of orientations, and the parameters of the weighting function parameters and bank of filters, with the noise and ill-conditioning parameters fixed to $$\begin{aligned}
 k_o =2, \varepsilon=0.0001\end{aligned}$$

The parameter vector in this example is $$\begin{aligned}
\upsilon_o= [ \begin{array}{ccccccc}
     c_o &g_o & \lambda_{\min o} & \sigma_{no} & \eta_o & N_o & O \\
     \end{array}]
    \end{aligned}$$ We set fixed upper and lower limits to $\upsilon_o$ as follows $$\begin{aligned}
&\underline{\upsilon}_o= [ \begin{array}{ccccccc}
    0.1 & 1.0 & 2.0 & 0.4 & 1.0 & 1 & 1 \\
        \end{array}]\\
&\overline{\upsilon}_o= [ \begin{array}{ccccccc}
       0.9 & 50.0 & 5.0 & 1.0 & 4.0 & 4 & 6 \\
     \end{array}]
    \end{aligned}$$

Because the Lena’s IFD has shown robustness to the norm type as demonstrated above, we only test the Frobenius norm-based optimization problems and . The resulting optimal values are given in Table \[table:resultex2\]. Notably, the optimal number of orientations where $O^*=1$ is achieved. Another noticeable result is that the optimal values for the weighting function parameters do not change, for which the boundary value at the upper and lower limits are found.

Figures \[fig:pc-Fnorm-mpm-full\] and \[fig:pc-Fnorm-pc-full\] show the images of $\mathcal{M}^*$ with associated optimal parameters $\upsilon_{o}^*$. By optimizing additional parameters of the 2D-MSPC, it is seen that many features are now better detected than Example 1. This is also consistent with the outcome of cost functions: the value of cost functions is larger than that in Example 1. Figure \[fig:diff\] shows the difference between the images obtained through and .

Example 3:
----------

In this example, magnetic resonance imaging (MRI) from the brain of a patient having multiple sclerosis (MS) is studied. MRI plays a key role in diagnosis and management of MS, [@Rodriguez2013]. MS is a disease that causes nerve damage in the brain and spinal cord. Characteristically, multi-focal plaques (lesions) can be seen using MRI, showing areas of brightness compared to the surrounding tissue (Figure \[fig:IM-0004-0054\_patient6\_T2-orig\]) [@Cabezas2014; @Schmidt2012; @Lorenzo2013]. Our experiments show that maximization of the minimum moment can increase the brightness of brain white matter in MRI. Thus, we chose $\mu_2=0$, and the cost function became $\mathcal{M}=M$. We solved the 2D-MSPC optimization for the $F-$norm with $\underline{\upsilon}_o$ and $\overline{\upsilon}_o$ as given in Example 2. The following optimal parameters were achieved with the optimal cost function $\|\mathcal{M}\|^*_F=234.0708$. $$\begin{aligned}
\upsilon_{o}^{*}= [ \begin{array}{ccccccc}
   49.9396 & 0.1001 & 4.7264 & 0.4 & 1.7499 & 3 & 1 \\
     \end{array}]
    \end{aligned}$$ Figure \[fig:IM-0004-0054\_patient6\_T2-pc\] shows the image of $\mathcal{M}$. It is seen that the location and size of lesions are clearly detectable by using the proposed optimal 2D-MSPC.

Example 4
---------

In this example, we show another MR image from a postmortem brain with MS obtained using a high-field MR scanner (Figure \[fig:NR01196s57R1-orig\]). Arrows indicate MS lesions, which have been confirmed by histological analysis.

Because maximization of the minimum moment increases the brightness of the non-lesion brain areas, we chose $\mu_2=0$. The cost function also became $\mathcal{M}=M$. We solve the 2D-MSPC optimization for the $F-$norm, with $\underline{\upsilon}_o$ and $\overline{\upsilon}_o$ as given in Example 2. The following optimal parameters were achieved with the optimal cost function $\|\mathcal{M}\|^*_F=262.1824$. $$\begin{aligned}
\upsilon_{o}^{*}= [ \begin{array}{ccccccc}
   49.9998 & 0.1 & 2 & 0.4 & 3.998 & 4 & 1 \\
     \end{array}]
    \end{aligned}$$ Figure \[fig:NR01196s57R1-pc\] shows the image of $\mathcal{M}^*$. It is seen that the MS lesions are clearly detectable.

[0.5]{} ![image](NR01196s57R1_orig2){width="1\linewidth"}

[0.5]{} ![image](NR01196s57R1_7Par_Fn_M_cost2){width="1\linewidth"}

Conclusions
===========

In this paper, we have focused on the IFD using a 2D-MSPC method. 2D-MSPC is originally proposed by Peter Kovesi and has shown great potential for the detection of various image features, particularly, lines, edges, corners, Mach bands, and blobs. However, the parameter setting of 2D-MSPC is typically performed manually based on trial and error, and studies for tuning of such parameters are limited. To enhance the application of this method, we have proposed several optimization frameworks for optimal and automatic tuning of the 2D-MSPC parameters. Through demonstration of several examples including MR images from patients with MS, we show that the ability of IFD can significantly be enhanced.

Matrix Norm Definition {#matrixnorm}
======================

([@Lyche2012], §8.1, p.177) \[normdef\] A function $\|.\|: {\ensuremath{\mathbb{C}}}^{q\times r}$ is called a matrix norm on $ {\ensuremath{\mathbb{C}}}^{q\times r}$ if for all $A,B \in  {\ensuremath{\mathbb{C}}}^{q\times r}$ and all $c \in {\ensuremath{\mathbb{C}}}$ $$\begin{aligned}
&\mbox{(Positivity)~}\|A\| \geq 0\\
&\mbox{(homogeneity)~}\| c A\| = |c| ~\|A\| \\
&\mbox{(subadditivity)~} \| A+B\| \leq \|A\|+\|B\|\end{aligned}$$

Proof of Lemma \[lempc\] {#p-lempc}
========================

The replacement of , and in results in $$\begin{aligned}
& M=\frac{1}{2}\Big(\textstyle \sum_o \mbox{PC}_o^2 \cos^2(\theta_o)+\sum_o \mbox{PC}_o^2 \sin^2(\theta_o)+\\
     & \sqrt{4\textstyle \sum_o  \mbox{PC}_o^4 \cos^2(\theta_o)\sin^2(\theta_o)+\textstyle \sum_o \mbox{PC}_o^4(\cos^2(\theta_o)-\sin^2(\theta_o))^2}\Big)\\
     &\hspace{5em}=\frac{1}{2}\Big(\textstyle \sum_o \mbox{PC}_o^2 ( \cos^2(\theta_o)+ \sin^2(\theta_o))+\\
   & \hspace{5em}~~ \sqrt{\textstyle \sum_o  \mbox{PC}_o^4 (\cos^2(\theta_o)+\sin^2(\theta_o))^2}\Big)\\&\hspace{5em}=\frac{1}{2}\Big(\textstyle \sum_o \mbox{PC}_o^2+\sqrt{\textstyle \sum_o  \mbox{PC}_o^4}\Big).\end{aligned}$$

The proof for the minimum moment is exactly the same and omitted.

Proof of Proposition \[thpcbnd\] {#p-pcbnd}
================================

We need the following Lemmas from linear algebra.

\[lemapb\] The following property holds for matrix norms. $$\begin{aligned}
\|A-B\| \leq \|A\| +\|B\|\end{aligned}$$ [*Proof:*]{} is straightforward using the matrix norm definition \[normdef\].

([@Lyche2012], §8.1.1, p.179) \[lemak\] For a consistent matrix norm $\|.\|_c$ on ${\ensuremath{\mathbb{C}}}^{q  \times q}$, the following inequality holds $$\begin{aligned}
\| A^k\|_c \leq \|A\|_c^k, \mbox{~for~} k\in {\ensuremath{\mathbb{N}}}.\end{aligned}$$

The proof of Proposition \[thpcbnd\] is as follows. By using Lemma \[lempc\], we have $$\|M\|_c=\frac{1}{2}\left\| \mbox{PC}_1^2+\ldots+\mbox{PC}_O^2+\sqrt{\mbox{PC}_1^4+\ldots+\mbox{PC}_O^4}\right\|_c$$ From Lemma \[lemapb\], thus $$\|M\|_c \leq \frac{1}{2}\Big(\left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\left\|\sqrt{\mbox{PC}_1^4+\ldots+\mbox{PC}_O^4}\right\|_c\Big)$$ By using Lemma \[lemak\], it yields $$\begin{aligned}
\|M\|_c& \leq \frac{1}{2}\Big( \left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\sqrt{\left\|\mbox{PC}_1^4+\ldots+\mbox{PC}_O^4\right\|_c}\Big)\\
& \leq \frac{1}{2}\Big(\left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\sqrt{\left\|\mbox{PC}_1^4\right\|_c+\ldots+\left\|\mbox{PC}_O^4\right\|_c}\Big)\end{aligned}$$ Lemma \[lemak\] is again used, which yields $$\begin{aligned}
&\|M\|_c \leq \\
&\frac{1}{2}\Big(\left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\sqrt{\left\|\mbox{PC}_1^2\right\|_c^2+\ldots+\left\|\mbox{PC}_O^2\right\|_c^2}\big)\end{aligned}$$ From algebra, $\sqrt{a_1^2+\ldots+a_O^2}\leq a_1+\ldots+a_O$ for $a_i\geq0$, thus, $$\begin{aligned}
\|M\|_c &\leq \frac{1}{2}\Big(\left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\left\|\mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c\Big)\\
&=\frac{1}{2}\Big(2\left\| \mbox{PC}_1^2\right\|_c+\ldots+2\left\|\mbox{PC}_O^2\right\|_c\Big)\end{aligned}$$ and the proof of is complete.

By using Lemma \[lempc\], we have $$\|m\|_c=\frac{1}{2}\Big(\left\| \mbox{PC}_1^2+\ldots+\mbox{PC}_O^2-\sqrt{\mbox{PC}_1^4+\ldots+\mbox{PC}_O^4}\right\|_c\Big)$$ From the matrix norm property, thus $$\|m\|_c \leq \frac{1}{2}\Big(\left\| \mbox{PC}_1^2\right\|_c+\ldots+\left\|\mbox{PC}_O^2\right\|_c+\left\|\sqrt{\mbox{PC}_1^4+\ldots+\mbox{PC}_O^4}\right\|_c\Big)$$ and the rest of the proof is exactly the same as above.

Proof of Proposition \[thsubopt\] {#p-suboptlem}
=================================

From the norm definition, , and , it is deduced that $$\begin{aligned}
 \left\| \mathcal{M}\right\|_c  \leq  (|\mu_1|+|\mu_2|) \textstyle \sum_{o} \left\| \mbox{PC}_o^2 \right\|_c\end{aligned}$$ A suboptimal solution to , with the consistent norm, is obtained by increasing the upper found of $\left\| \mathcal{M}\right\|_c$. This completes the proof.

Acknowledgment {#acknowledgment .unnumbered}
==============

First author would like to thank Shrushrita Sharma and Glen Pridham for some discussions. Funding from the MS Society of Canada, Natural Sciences and Engineering Council of Canada, and Alberta Innovates – Health Solutions is also acknewledged.

[1]{}

Peter Kovesi, [*Invariant Measures of Image Features From Phase Information*]{}, PhD Thesis, The University of Western Australia, 1996.

D. Ziou and S. Tabbone, [*Edge detection techniques: An overview*]{}, International Journal of Pattern Recognition and Image Analysis, 8(4), pp. 537–559, 1998.

T.P. Patel, S.R. Panchal, Corner Detection Techniques: An Introductory Survey, [*International Journal of Engineering Development and Research*]{}, Vol. 2, No. 4, pp. 3680–3686, 2014.

J.A. Noble, [*Descriptions of image surfaces*]{}, D.Phil thesis, Department of Engineering Science, University of Oxford, 1989.

J.F. Canny, [*Finding edges and lines in images*]{}, Master’s thesis, Massachusetts Institute of Technology, AI Lab. TR-720, 1983.

A. Blake, A. Zisserman, [*Visual Reconstruction*]{}, MIT Press, Cambridge, MA, 1987.

P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, [*IEEE Transactions on Pattern Analysis and Machine Intelligence*]{}, Volume 12, Issue 7, pp. 629 - 639, 1990.

S.M. Smith, J.M. Brady, SUSAN: A New Approach to Low Level Image Processing, [*International journal of computer vision*]{}, Volume 23, Issue 1, pp 45–78, 1997.

M.C. Morrone, D.C. Burr, J. Ross, R. Owens, Mach bands are phase dependent, [*Nature*]{}, 324, 250–253, 20 November 1986.

P. Kovesi, Image Features From Phase Congruency, [*Videre: Journal of Computer Vision Research, MIT Press*]{}. Volume 1, Number 3, pp. 1-25, Summer 1999.

P. Kovesi, Phase Congruency Detects Corners and Edges, [*The Australian Pattern Recognition Society Conference: DICTA 2003*]{}, pp 309–318, Sydney. December 2003.

A.V. Oppenheim, J.S. Lim, The importance of phase in signals, [*IEEE Proceedings*]{}, vol. 69, pp. 529–541, 1981.

S. Venkatesh, R. Owens, On the classification of image features, [*Pattern Recognition Letters*]{}, Volume 11, Issue 5, Pages 339–349, May 1990.

M.C. Morrone, R. Owens, Feature detection from local energy, [*Pattern Recognition Letters*]{}, 6, pp. 303–313, 1987.

G. Cao, P. Shi, B. Hu, Ultrasonic Liver Discrimination Using 2-D Phase Congruency, [*IEEE Transactions on Biomedical Engineering*]{}, Vol. 53, No. 10, pp. 2116 – 2119, 2006.

V. Štruc N. Paveši' c, Phase congruency features for palm-print verification, [*IET Signal Processing*]{}, Vol. 3, Issue 4, pp. 258 – 268, 2009.

L. Zhang, L. Zhang, D. Zhang, Z. Guo, Phase congruency induced local features for finger-knuckle-print recognition, [*Pattern Recognition*]{}, Volume 45, Issue 7, pp. 2522 – 2531, 2012.

B. Obara, M. Fricker, D. Gavaghan, V. Grau, Contrast-Independent Curvilinear Structure Detection in Biomedical Images, [*IEEE Transactions on Image Processing*]{}, Vol. 21, No. 5, pp. 2572–2581, 2012

T. Mouats, N. Aouf, M.A. Richardson, A Novel Image Representation via Local Frequency Analysis for Illumination Invariant Stereo Matching, [*IEEE Transactions on Image Processing*]{}, Vol. 24, No. 9, pp. 2685-2700, 2015.

O.M. Rijal, H. Ebrahimian, N.M. Noor, A. Hussin, A. Yunus, and A.A. Mahayiddin, Application of Phase Congruency for Discriminating Some Lung Diseases Using Chest Radiograph, [*Computational and Mathematical Methods in Medicine*]{}, Article ID 424970, 15 pages, 2015.

H. Ziaei Nafchi ; R. Farrahi Moghaddam ; M. Cheriet, Phase-Based Binarization of Ancient Document Images: Model and Applications, [*IEEE Transactions on Image Processing*]{}, Vol. 23, No. 7, 2916 – 2930, July 2014.

[ http://www.peterkovesi.com]( http://www.peterkovesi.com). W.S. Levine, [*The Control Handbook*]{}, CRC Press in Cooperation with IEEE Press, 1995.

V. V. Fedorov, [*Theory of optimal experiments*]{}, Elsevier, 1972.

G. C. Goodwin and R. L. Payne, [*Dynamic system identification: experiment design and data analysis*]{}, Academic press, 1977.

A. C. Atkinson, Developments in the design of experiments, [*International Statistical Review*]{}, vol. 50, pp. 161 – 177, 1982.

L. Pronzato and E. Walter, Robust experiment design via stochastic approximation, [*Mathematical Bio-sciences*]{}, vol. 75, no. 1, pp. 103–120, 1985.

M.P.F. Berger, D-optimal sequential sampling designs for item response theory models, [*Journal of Educational Statistics*]{}, vol. 19, no. 1, pp. 43–56, 1994.

F. Pukelsheim, [*Optimal design of experiments, ser. Prob- ability and mathematical statistic*]{}, Wiley, 1993.

R. A. Fisher, [*The Design of Experiments*]{}, 9th ed., Macmillan, 1971 (1935 1st ed.).

H. Cramér, [*Mathematical Methods of Statistics*]{}, Princeton Univ. Press. Princeton, NJ. 1946.

C.R. Rao, Information and the accuracy attainable in the estimation of statistical parameters. [*Bulletin of the Calcutta Mathematical Society*]{}, 37, p. 81–89, 1945.

A. Atkinson, A. Donev, R. Tobias, [*Optimum Experimental Designs, With SAS*]{}, Oxford University Press, 2007.

R.C. St.John and N. R. Draper, D-optimality for regression designs: a review, [*Technometrics*]{}, vol. 17, no. 1, pp. 15 – 23, 1975.

T. Lyche, [*Lecture Notes for Inf-Mat 4350, 2012*]{}, University of Oslo, August 2012.

R. Schatten, [*Norm Ideals of Completely Continuous Operators*]{}, Springer-Verlag, Berlin, Gottingen, Heidelberg, 1960.

M. Rodriguez, O.H. Kantarci, I. Pirko, [*Multiple Sclerosis*]{}, Oxford University Press, 2013.

M. Cabezas, A. Oliver, E. Roura, J. Freixenet, J. C. Vilanova, L. Ramió-Torrentà, À. Rovira, X. Lladó, Automatic multiple sclerosis lesion detection in brain MRI by FLAIR thresholding, [*Computer Methods and Programs in biomedicine*]{}, 115, pp. 147 – 161, 2014.

P. Schmidt, C. Gaser, M. Arsic, D. Buck, A. Förschler, A. Berthele, M. Hoshi, R. Ilga, V. J. Schmid, C. Zimmer, B. Hemmer, M. Mühlau An automated tool for detection of FLAIR-hyperintense white-matter lesions in Multiple Sclerosis, [*NeuroImage*]{}, Volume 59, Issue 4, pp. 3774–3783, 2012.

D. García-Lorenzo, S. Francis, S. Narayanan, D. L. Arnold, D. L. Collins, Review of automatic segmentation methods of multiple sclerosis white matter lesions on conventional magnetic resonance imaging, [*Medical Image Analysis*]{}, Volume 17, Issue 1, Pages 1–18, January 2013.

[^1]: S.M.M. Alavi is with the Hotchkiss Brain Institute and Department of Clinical Neurosciences at the Cumming School of Medicine, University of Calgary, Canada. Email: mahdi.alavi@ucalgary.ca.

[^2]: Y. Zhang is with the Hotchkiss Brain Institute, Department of Clinical Neurosciences, and Department of Radiology at the Cumming School of Medicine, University of Calgary. Email: yunyan.zhang@ucalgary.ca.
