**2.2. Geometrical analysis**

In the geometric analysis, there are several stages that can be determined as Figure 5 below:

**Figure 5.** Flowchart of Geometrical Analysis

Users must consider a crater as an ellipse in a real image. This method of consideration will convert an ellipse in a 2D image into a circle on a plane using Conical Projection Analysis. Any ellipse will appear to be a circle from a certain point of views. In other words, an ellipse will be projected into a circle at a certain projection point. At the final stage, this method will be able to calculate the orientation and position of a crater (disc in shape) that is being detected before through the proposed detection algorithm.

#### *2.2.1. Fundamentals of ellipse and rotation matrix*

Mathematically, an ellipse can be defined as the locus of all points on the plane whose distances R1 and R2 (as Figure 6 below) to two fixed points added to the same constant and can be notified as:

$$\mathbf{R1} + \mathbf{R2} = \mathbf{2a} \tag{3}$$

where a = semi major axis and the origin of the coordinate system is at one of the foci (-c,0) and (c,0). These two focis are chosen to be identical with the bounding ellipse algorithm equation. It is sometimes defined as a conical section from cutting a circular conical or cylindrical surface with an oblique plane. There are five vital parameters in ellipse including semi-major axis denoted as **a**, semi-minor axis denoted as **b**, centre, **c** of ellipse in Xcoordinate, **Xc**, centre of ellipse in Y-coordinate, **Yc**, and an angle of rotation denoted as **ω**. The ellipse with all the parameters can be illustrated as below:

**Figure 6.** Ellipse

492 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

technical section and experimental results section.

**2.2. Geometrical analysis** 

**Figure 5.** Flowchart of Geometrical Analysis

detected before through the proposed detection algorithm.

blobs together in a single picture. As a result, the final image will comprise a group of clusters or patches that correspond to craters on the moon's surface. These groups of blobs (light and dark patches) will then be used to measure the minimum distance and angle between each of them. First, the algorithm will calculate all the distances of every patch and will pick only the minimum distance. The light patches with minimum distances that are attached to the dark patches will be considered craters as a first step. Second, every angle between the known input sun vector and pairing patches vector is calculated using a scalar product or dot product. Technically, all these methods will be elaborated further in the

In the geometric analysis, there are several stages that can be determined as Figure 5 below:

Users must consider a crater as an ellipse in a real image. This method of consideration will convert an ellipse in a 2D image into a circle on a plane using Conical Projection Analysis. Any ellipse will appear to be a circle from a certain point of views. In other words, an ellipse will be projected into a circle at a certain projection point. At the final stage, this method will be able to calculate the orientation and position of a crater (disc in shape) that is being An ellipse that lies along the horizontal X-axis with foci points F1 (-c, 0) and F2(c, 0) as can be shown in Figure (6) above, will have an equation of

$$\mathbf{x}^2 \langle \mathbf{a}^2 + \mathbf{y}^2 \rangle \langle \mathbf{a}^2 \cdot \mathbf{c}^2 \rangle = 1 \text{ where } \mathbf{a} \ge \mathbf{c} \text{ for the ellipse} \tag{4}$$

For an ellipse, the distance c between the centre and a focus is less than the distance a between the centre and foci, so a2-c2 is positive and a new constant b>0 is introduced by setting [2]:

$$\mathbf{b}^2 \equiv \mathbf{a}^2 - \mathbf{c}^2 \text{ for ellipses} \tag{5}$$

Hence the equation of an ellipse with F1 (-a, 0) and F2 (a, 0) is simplified to

$$\mathbf{x}^2/\mathbf{a}^2 + \mathbf{y}^2/\mathbf{b}^2 = \mathbf{1} \text{ where } \mathbf{0} \le \mathbf{b} \le \mathbf{a} \tag{6}$$

For both the hyperbola and the ellipse, a number e, called the eccentricity is introduced by setting [2]:

$$\mathbf{e} = \mathbf{c}/\mathbf{a} \text{ or } \mathbf{e} = \sqrt{} \left( \mathbf{a}^2 \mathbf{b}^2 \right) \text{\AA} \tag{7}$$

In this mathematical and geometrical analysis, the authors started to brief in ellipse equations and rotation matrix first which are going to be analyzed soon in the bounding ellipse algorithm and for the reconstruction of ellipse to a circle on a 2-D plane. These methods are beneficial to determine the orientation and the position of the spacecraft during Entry, Descent and Landing (EDL) applications. In this case, the rotation matrix for an ellipse can be illustrated as Figure 4 below:

**Figure 7.** Rotation Ellipse

In S frame as in figure 7 above, it can be shown/demonstrated that using a standard ellipse equation, the x" and y" can be expressed as

$$(\mathbf{x''})^2/\mathbf{a}^2 + (\mathbf{y''})^2/\mathbf{b}^2 = 1 \text{ where } 0 \le \mathbf{b} \le \mathbf{a} \tag{8}$$

Where

$$
\begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} \cos \phi & \sin \phi \\ -\sin \phi \cos \phi \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \tag{9}
$$

$$= \begin{pmatrix} \ x' \cos \phi & y' \sin \phi \\ -\ x' \sin \phi \ y' \cos \phi \end{pmatrix} \tag{10}$$

and

$$
\begin{pmatrix} \mathbf{x}^\* \\ \mathbf{y}^\* \end{pmatrix} = \begin{pmatrix} \mathbf{x} - \mathbf{x}o \\ \mathbf{y} - \mathbf{y}o \end{pmatrix} \tag{11}
$$

After substitution of these three vital equations, new formula for the rotation ellipse, which is also the bounding ellipse equation:

$$\frac{\left(\left(\mathbf{x} - \mathbf{x}o\right)\cos\alpha + \left(\mathbf{y} - \mathbf{y}o\right)\sin\alpha\right)^2}{a^2} + \frac{\left(\left(-\mathbf{x} + \mathbf{x}o\right)\sin\alpha + \left(\mathbf{y} - \mathbf{y}o\right)\cos\alpha\right)^2}{b^2} = 1\tag{12}$$

Again, the rotation ellipse can also be expressed through this formula:

$$\left(\tilde{\mathbf{x}} - \tilde{\mathbf{c}}\right)^{T} A \{\tilde{\mathbf{x}} - \tilde{\mathbf{c}}\} = \left(\mathbf{x} - \mathbf{x}o\right)^{2} A\_{11} + \left(\left(\mathbf{x} - \mathbf{x}o\right)\left(y - yo\right)\right) \left(A\_{21} + A\_{12}\right) + \left(y - yo\right)^{2} A\_{22} = 1\tag{13}$$

Where 11 21 12 22 *A* ,,, *AAA* is the E matrix that is determined from the bounding ellipse algorithm following the form of

$$A = \begin{pmatrix} A\_{11} \ A\_{12} \\ A\_{21} \ A\_{22} \end{pmatrix} \tag{14}$$

Thus, by comparing the equations of 12 and 13 above, the authors express all the ellipse parameters a,b,c (Xc and Yc) and ω in terms of this 11 21 12 22 *A* ,,, *AAA* entities when drawing the bounding ellipse around the target patch which:

$$\mathfrak{a} = (\frac{2\sin 2\alpha}{(A\_{11} + A\_{22})\sin 2\alpha + 2A\_{12}})^{1/2} \tag{15}$$

$$ab = \frac{2\sin 2\alpha}{(A\_{11} + A\_{22})\sin 2\alpha - 2A\_{12}})^{1/2} \tag{16}$$

$$o = \tan^{-1}(\frac{-2A\_{12}}{A\_{22} - A\_{11}}) / 2 \tag{17}$$

Furthermore, the other two ellipse parameters Xc and Yc values can be determined straight away from the bounding ellipse algorithm.

#### *2.2.2. Bounding ellipse algorithm*

494 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

ellipse can be illustrated as Figure 4 below:

**Figure 7.** Rotation Ellipse

Where

and

equation, the x" and y" can be expressed as

'

In this mathematical and geometrical analysis, the authors started to brief in ellipse equations and rotation matrix first which are going to be analyzed soon in the bounding ellipse algorithm and for the reconstruction of ellipse to a circle on a 2-D plane. These methods are beneficial to determine the orientation and the position of the spacecraft during Entry, Descent and Landing (EDL) applications. In this case, the rotation matrix for an

In S frame as in figure 7 above, it can be shown/demonstrated that using a standard ellipse

" cos sin " sin cos

 

'cos 'sin

'

*x y x y* 

*x x xo y y yo* 

'

*x y* 

'

'sin 'cos

 

*x y*

(x")2/a2 + (y")2/b2 = 1 where 0 < b < a (8)

(9)

(10)

(11)

e = c/a or e = √ (a2-b2)/a (7)

This is the first method used in the geometrical analysis and the reason of using it is to get the bounding ellipse around the targeted patch. This is further used in a final stage of this section that is reconstruction or projection of the ellipse to a circle in a 2-D plane. The bounding ellipse takes five inputs of the ellipse parameters described previously, which are semi-major axis a, semi-minor axis b, centre of ellipse in x-coordinate Xc, centre of ellipse in y-coordinate Yc and rotation angle ω. These five parameters are embedded in the A matrix produced by this algorithm as an output along with the Xc and Yc values. In order to draw the bounded ellipse the authors need to express all those ellipse parameters in terms of the entities 11 21 12 22 *A* ,,, *AAA* of the A matrix. The derivations of these terms are comprehensively described in the previous section. The results of this bounding ellipse algorithm will then be used to reconstruct this 2-D ellipse to a 2-D circle in a plane using Reconstruction Ellipse Algorithm. There are two vital equations of reconstructing a disc that will be used in this algorithm to determine the orientation of the spacecraft which is described below [18]; nonetheless in this chapter the authors will only determine the orientation:

$$\tilde{q} = \frac{B}{\sin \alpha} \left| \cos \beta \tilde{\mathcal{R}}\_1 \pm\_1 \cos \alpha \sin \beta \sqrt{\tan^2 \alpha - \tan^2 \beta} \tilde{\mathcal{R}}\_2 \right| \tag{18}$$

$$\tilde{p} = \pm\_2 \cot a \left\{ \frac{\sin \beta}{\cos a} \tilde{R}\_1 \pm\_1 \cos \beta \sqrt{\tan^2 \alpha - \tan^2 \beta} \tilde{R}\_2 \right\} \tag{19}$$

Where:

*q* = is the position of the reconstructed ellipse *p* = is the orientation of the reconstructed ellipse B = is the radius of a disc (craters that the authors model as a disc) α= is the arc length of the semi-major axis β= is the arc length of the semi-minor axis *R*1 = Rotation matrix of the column vector *R*2 =Rotation matrix of the column vector

This is the reconstruction or projection ellipse equations where the authors consider an ellipse as a half-length along the axis symmetry, which is taken to 0 that is **A**=0. In this case, the authors need to model the crater bounded as a disc. That is the reason for the half-length along the axis symmetry **A**, is taken to 0. **A** in this case is not the attributes of the matrix determined previously. The authors have to deal with those two equations above where as can be seen the equation *p* , the orientation of a disc, is independent of parameters B (the radius) which means that the authors are able to determine quickly the reconstruction algorithm in order to identify the orientation of the spacecraft relative to the moon's surface. There is the ambiguity case in equation *p* which is the negative and positive case of 2 sign. In this case, the *p* equations always takes a positive value instead of negative as the crater's orientation is just pointing upward towards a camera or in other words, a disc will only be seen if its outward face points towards the camera rather than away from the camera [18]. This is a discerning case; when one considers how human calculates an object's position, its exact size is needless in finding the direction of neither its centre nor its orientation. In contrast, the equation *q* is dependent on the radius of a disc, B which the authors have no knowledge of the radius of the craters and how far it is from the moon's surface. Therefore, searching for a solution q is actually an interesting future work that needs to be achieved if the authors want to find the position of the spacecraft during the EDL operations.

#### **3. Technical sections**

There are different mathematical equations and fundamentals applied during this project development. In order to make the project runs smooth as planned; the authors have divided the logical structures of the project into three different sections:

#### **3.1. Craters detection algorithm**

496 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

nonetheless in this chapter the authors will only determine the orientation:

cos

B = is the radius of a disc (craters that the authors model as a disc)

sin *<sup>B</sup> q R*

*q* = is the position of the reconstructed ellipse *p* = is the orientation of the reconstructed ellipse

α= is the arc length of the semi-major axis β= is the arc length of the semi-minor axis

= Rotation matrix of the column vector

=Rotation matrix of the column vector

**3. Technical sections** 

Where:

*R*1

*R*2

used to reconstruct this 2-D ellipse to a 2-D circle in a plane using Reconstruction Ellipse Algorithm. There are two vital equations of reconstructing a disc that will be used in this algorithm to determine the orientation of the spacecraft which is described below [18];

> 

*pR R*

 2 2 1 1 <sup>2</sup> cos cos sin tan tan

2 1 1 2 sin cot cos tan tan

This is the reconstruction or projection ellipse equations where the authors consider an ellipse as a half-length along the axis symmetry, which is taken to 0 that is **A**=0. In this case, the authors need to model the crater bounded as a disc. That is the reason for the half-length along the axis symmetry **A**, is taken to 0. **A** in this case is not the attributes of the matrix determined previously. The authors have to deal with those two equations above where as can be seen the equation *p* , the orientation of a disc, is independent of parameters B (the radius) which means that the authors are able to determine quickly the reconstruction algorithm in order to identify the orientation of the spacecraft relative to the moon's surface. There is the ambiguity case in equation *p* which is the negative and positive case of 2 sign. In this case, the *p* equations always takes a positive value instead of negative as the crater's orientation is just pointing upward towards a camera or in other words, a disc will only be seen if its outward face points towards the camera rather than away from the camera [18]. This is a discerning case; when one considers how human calculates an object's position, its exact size is needless in finding the direction of neither its centre nor its orientation. In contrast, the equation *q* is dependent on the radius of a disc, B which the authors have no knowledge of the radius of the craters and how far it is from the moon's surface. Therefore, searching for a solution q is actually an interesting future work that needs to be achieved if

the authors want to find the position of the spacecraft during the EDL operations.

divided the logical structures of the project into three different sections:

There are different mathematical equations and fundamentals applied during this project development. In order to make the project runs smooth as planned; the authors have

 

 

(19)

2 2

(18)

 *R*

  This project involved the development of the detection algorithm using MATLAB image processing tool and an image of the moon's surface. It is mainly based on the binary image and morphological image analysis. At the first stage, the authors introduce the concept of HSV (Hue, Saturation, and Value) as well as morphology image investigation such as dilation and erosion to exploit the real image. The Hue is expressed as an angle around a colour hexagon mainly using the real axis as 0° axis. The Value is measured along the cone and has two different conditions. If V=0, the end of the axis is represented by black and white if V=1at the end of the axis is represented by white [4]. The Saturation is the purity of the colour and is measured as the distance from the V axis looking at the hue, saturation and value hexagonal cone.

The HSV colour system is based on cylindrical coordinates. Mathematically, converting form RGB (Red, Green and Blue) image to HSV is actually developing the equation from the Cartesian coordinates to cylindrical coordinates. To reduce the complicatedness of analysis on image detection, the authors analysed a 2-D optical image. The threshold for the image is set using intelligent approach from Sawabe, Natsunaga and Rokugawa in classifying the images as can be shown in equations (20) below. By using this approach, images were classified into two components that are light and dark patches or obviously known as 'sunny and shady patches'. Ideally, these two groups of patches are easily recognizable if the image was taken under low sun elevation. These patterns of light and dark patches were detected when all these equations [13] are satisfied:

$$
\mathbf{R\_{min}} \circ \mathbf{R\_n} + \sigma
$$

$$
\mathbf{R\_{max}} \circ \mathbf{R\_n} + \sigma
$$

$$
\mathbf{P\_{min}} \circ \mathbf{P\_{max}}
$$

Where Rmin indicates the minimum pixel value, Rmax indicates the maximum pixel value; Rm indicates the average pixel value and indicates the standard deviation in the small area including the targeted patches. Pmin and Pmax indicate the positions at the minimum and maximum value pixels from the direction of the sun radiation or sun vector [13].

Apart from that, there are two basic fundamentals on morphological operation applied in the algorithm which are dilation and erosion. Dilation is an operation that grows or expands objects in a binary image. This can be represented by a matrix of 1's and 0's but usually it is convenient to show only the 1's group. The dilation between two images A and B, is denoted by A B and is defined as [1]:

$$\mathbf{A}\ddot{A}\mathbf{B} = \{ \mathbf{z} \mid \begin{pmatrix} \mathbf{B} \end{pmatrix}\_{\mathbf{z}} \cap \mathbf{A} \neq \mathbf{f} \}\tag{21}$$

Nevertheless, erosion shrinks or thins objects in a binary image, which is the opposite case of the dilation operation. The mathematical definition of erosion can be expressed as [1] and the result of eroded and dilated image is shown in Figure 8 as below. This can be compared to the original image in Figure (9):

$$\mathbf{A} \oplus \mathbf{B} = \{ \mathbf{z} \mid \begin{pmatrix} \mathbf{B} \end{pmatrix}\_{\mathbf{z}} \cap \mathbf{A}^{\mathbf{c}} \neq \mathbf{f} \} \tag{22}$$

**Figure 8.** Image after dilation and erosion were applied

The previous stages are then followed by another image analysis methodology that is *regionprops* function and this region indicates the patch that the authors want to analyze. The reason for using this region descriptors function is to plot the centroids of each light and dark patch and gather them (light and dark patches) back together in a single picture which will then be used to further calculate the desirable minimum distance and angle between them to the known sun direction. *Regionprops* takes the labelled matrix as an input and uses axis-x and axis-y to describe horizontal and vertical coordinates respectively.

To classify and consider this region of interest as a crater, the authors proposed two ways of detection which are minimum distance measurement and angle detection based on the known sun vector. In the distance measurement, the minimum distances between each of the centroids calculated previously are determined using this formula:

$$\vdash \mathbf{Dístance} \mid = \vee[(\mathbf{x}^2 \cdot \mathbf{x}^1)^2 + (\mathbf{y}^2 \cdot \mathbf{y}^1)^2] = \mid \mathbf{r} \mid \tag{23}$$

where |**Distance|** is the measurement of distance between two pairing patches (light and dark), **x**2 and **x**1 are the x-component of the centroids and **y**2,**y**1 are the centroid's component of the y-axis respectively. In the angle determination, the authors give an input for the sun vector which is known by looking at the sunray effect at those craters (the position of the light and dark shades). This algorithm will then compute each of the angles of every pairing blob with their minimum distances to the sun vector added input using scalar product or dot product in the vector analysis which is:

$$\mathbf{r.s} = \begin{vmatrix} \mathbf{r} & \mathbf{s} \end{vmatrix} \text{ } \text{s} \begin{vmatrix} \mathbf{s} \end{vmatrix} \text{ } \text{ } \tag{24}$$

Where **r** = vector of each pairing blobs **s =** sun vector |**r**| = distance/length of each pairing blobs |**s**| = distance of sun vector = unit vector = 1 θ = angle between sun vector and vector of each pairing blobs

Each of the pairing blobs angle is calculated using the above equation and those who has a minimum angle to with known sun vector and has a minimum distance calculated previously will be considered and stored as a crater. Oppositely, those who are against the direction and have the maximum angles to the sun vector will be scrapped and considered noise. At the end of these two measurements (distance and angle), the authors managed to get the best eight craters using this reliable craters detection algorithm.
