**3.2. Geometrical analysis**

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

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

dot product in the vector analysis which is:

<sup>c</sup>

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

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

 |**Distance| =** √[(**x**2-**x**1)2 + (**y**2-**y**1)2] = |**r**| (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

**rs r s** . cos

(24)

axis-x and axis-y to describe horizontal and vertical coordinates respectively.

the centroids calculated previously are determined using this formula:

<sup>z</sup> A B {z | B A f} (22)

#### *3.2.1. Bounding ellipse algorithm*

The geometrical part is mainly based on the mathematical analysis on vector calculus, the rotation matrix, ellipse equations and also the mathematical applications to determine two major features that is the orientation of a crater that has been modelled as a disc and the position of the disc after projection. These two features make a vital solution to have a safe landing on the moon's surface. Firstly, from the previous codes, the authors choose one of the best eight in a crater's list to run a bounding ellipse algorithm around the targeted crater. The bounding ellipse algorithm input the P matrix which is composed from the targeted crater itself. The tolerance of the bounding image is set to get more precise result. As discussed above in the methodology section, the output of this algorithm will be the A

matrix and c, the centre of the ellipse. A matrix is in the form of 11 12 21 22 <sup>A</sup> *A A A A* . The reason

the authors highlighted this is to show that another important parameter in ellipse is embedded in this A matrix. They are semi-major axis **a**, semi-minor axis **b**, and angle of rotation **ω**. As a result, the final equations (equations 13, 14 and 15) are used to draw the ellipse to bind the targeted patches as derived from the previous methodology section. The next step is to draw the ellipse using another algorithm and all those parameters above and also the centre of the ellipse that can be determined straight away from the bounding ellipse algorithm. Statistically, after the algorithm has calculated all the values of the parameter, the authors have the same value for semi-major axis **a** and semi-minor axis **b** which suggests that a circle is formed. This important knowledge tells the authors that the camera is actually pointing straight around 90 degrees vertically to the targeted crater on the moon's surface. If the authors have an ellipse instead of a circle, it means that the camera is not pointing straight 90 degrees downward away from it (crater). That is the first assumption the authors made before doing the third stage algorithm that is Reconstruction of an Ellipse to a circle in a 2-D plane.

#### *3.2.2. Reconstruction of an ellipse on image plane*

First, the input of this algorithm is the five primary parameters of the ellipse that are **a**, **b**, **Xc**, **Yc**, **ω**, the half-length of the semi major axis is denoted as capital **A**, and the half-length of the semi-minor axis is denoted as capital **B**. This reconstruction vision is to model the craters or projected ellipse as a disc. Mathematically, the concepts are adapted by taking the focal point as an origin of a camera frame and the x-axis is aligned with the focal axis of the camera and x > 0 is what the camera is looking for. Furthermore, the image plane is always assumed to be in front of the focal point rather than those in practice. Taking this objective into account, the half-length of the semi major axis **A** is taken to 0 because it is a disc. In comparison, the half-length of the semi-minor axis **B** is actually the radius of the crater or a disc. This reconstruction or projection ellipse algorithm is based on Stephen's proposed complex method on reconstruction spheroid. The authors will use the equation 17 previously to determine the orientation of a disc or crater that the authors modelled from the reconstruction algorithm.

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

*p* equation describes the orientation of the reconstructed disc (crater) and equation *q* describes the position of a disc (crater) after reconstruction or projection. As the above two equations, there are four possible solutions (regardless of the ambiguity case of <sup>2</sup> as discussed before in the methodology section), two of them are from equation *p* and another two from equation *q* . The equation *p* is independent from the half-length semi minor axis B as can be seen from the equation. Therefore, the authors can determine the orientation of a disc straight away from the algorithm. Unfortunately, the equation *q* is dependent on the half-length of semi-minor axis B which the authors are not aware of the radius of a crater because the authors are not able to determine the distance to the moon's surface. This is the main challenge. But getting the value of *p* will then lead the authors to get the position of *q* by some means because they are related. The position *q* of the disc will then be a greater future work to be explored.

#### **4. Experimental results and discussion**

#### **4.1. Craters detection algorithm result**

In this section, some experimental results are reported which were obtained by applying craters detection algorithm. In particular, the authors choose real images (optical images) of the moon's surface and develop the codes mainly based on the binary image analysis and morphological techniques. As mentioned before, this is the independent algorithm which does not depend on the data elevation map as well as past and future imaging data, which is the advantage compared to other detection algorithms proposed previously. There are seven stages to develop this algorithm as mentioned in the methodology (refer to the flowchart). The authors will discuss the results from the first stage until the last stage by attaching the result image for a better view and understanding of the concept used for analysis.

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

First, the input of this algorithm is the five primary parameters of the ellipse that are **a**, **b**, **Xc**, **Yc**, **ω**, the half-length of the semi major axis is denoted as capital **A**, and the half-length of the semi-minor axis is denoted as capital **B**. This reconstruction vision is to model the craters or projected ellipse as a disc. Mathematically, the concepts are adapted by taking the focal point as an origin of a camera frame and the x-axis is aligned with the focal axis of the camera and x > 0 is what the camera is looking for. Furthermore, the image plane is always assumed to be in front of the focal point rather than those in practice. Taking this objective into account, the half-length of the semi major axis **A** is taken to 0 because it is a disc. In comparison, the half-length of the semi-minor axis **B** is actually the radius of the crater or a disc. This reconstruction or projection ellipse algorithm is based on Stephen's proposed complex method on reconstruction spheroid. The authors will use the equation 17 previously to determine the orientation of a disc or crater that the authors modelled from

> 2 1 1 2 sin cot cos tan tan

*p* equation describes the orientation of the reconstructed disc (crater) and equation *q* describes the position of a disc (crater) after reconstruction or projection. As the above two equations, there are four possible solutions (regardless of the ambiguity case of <sup>2</sup> as discussed before in the methodology section), two of them are from equation *p* and another two from equation *q* . The equation *p* is independent from the half-length semi minor axis B as can be seen from the equation. Therefore, the authors can determine the orientation of a disc straight away from the algorithm. Unfortunately, the equation *q* is dependent on the half-length of semi-minor axis B which the authors are not aware of the radius of a crater because the authors are not able to determine the distance to the moon's surface. This is the main challenge. But getting the value of *p* will then lead the authors to get the position of *q* by some means because they are related. The position *q* of the disc will then be a greater

In this section, some experimental results are reported which were obtained by applying craters detection algorithm. In particular, the authors choose real images (optical images) of the moon's surface and develop the codes mainly based on the binary image analysis and morphological techniques. As mentioned before, this is the independent algorithm which does not depend on the data elevation map as well as past and future imaging data, which is the advantage compared to other detection algorithms proposed previously. There are seven stages to develop this algorithm as mentioned in the methodology (refer to the flowchart). The authors will discuss the results from the first

 

*pR R*

cos

2 2

(25)

 

*3.2.2. Reconstruction of an ellipse on image plane* 

the reconstruction algorithm.

future work to be explored.

**4. Experimental results and discussion** 

**4.1. Craters detection algorithm result** 

Erosion and dilation are two fundamentals in Morphological Image Analysis. In Figures 10 and 11, the erosion and dilation are experimented in a combined process for dark patches and light patches. For an ideal result, the authors just used the eroded image for both light and dark patches by adding them together in a single picture for centroid detection later. If the authors take consideration of both the dilation and erosion process as a combined process, the image itself will produce too much noise as shown in the figure above. For centroid determination, therefore, the authors want the noise to be kept at a minimum level in order to attach valid pairing patches and in order to produce minimum small blobs or patches that do not have their pairs (light and dark patch). This is also to reduce the processing time (running time) and complexity of centroid calculation.

**Figure 9.** Real Image of Craters (Optical Image)

**Figure 10.** Erosion and Dilation for dark patches

**Figure 11.** Erosion and Dilation for light patches

**Figure 12.** Threshold Image

**Figure 10.** Erosion and Dilation for dark patches

**Figure 11.** Erosion and Dilation for light patches

Secondly, after converting the RGB image plane to HSV image plane and taking into consideration only the '*Value*' from the HSV image plane, thresholding is then used for image segmentation proposal. Thresholding is used firstly to differentiate the object (foreground) from the background. The targeted foreground is ideally the craters themselves that are composed by light and dark patches pattern. The light patch occurs when the pixel value is more than the calculated threshold brightness using the same approach that is fully described in the technical section. Whilst the dark patch is detected when the pixel value is less than the same threshold brightness calculated before. The threshold image can be shown in Figure (12) above.

**Figure 13.** Centroid Determination using 'regionprops' function

Centroid determination stage is completed by using the 'regionprops' function in Matlab and this is basically to compute the desired region properties or targeted properties. This step also attached the pairing patches together to form a complete crater for further analysis. Regionprops only takes the label matrix as an input. Therefore, the authors have to preprocess the image by labelling it first using 'bwlabel' function along with regionprops function later. In MATLAB, bwlabel is used to label the connected components in binary image and bwlabel also supports the 2-D inputs only. After labeling the entire target then only the authors can use them in a memory with regionprops. Regionprops takes this labelled component as an input and return them with an image of centroid determination labelled as a cross (\*) symbol as in Figure (13) above. For a smoother image, the authors apply another function in MATLAB called bwareaopen to remove all of the unnecessary objects or all connected components (objects) that are fewer than P pixels set, producing another binary image.

The final stage of this craters detection algorithm is minimum radial distance and angle calculation. The authors have tested the algorithm with two different image conditions (different sun angles and background noises). According to the algorithm, the authors can select the number of craters to be detected by the system. The users can choose how many craters to be detected by the system based on the original image taken. As for this case, it will select the eight best craters that satisfy the minimum distance and angle conditions by assuming that the sun direction is known. An image of the landing site on the moon's surface has to be captured first and the amounts of craters needs to be detected are calculated manually. Based? on the original image, the authors have assumed the sun direction by identifying the shade and sunny pattern locations that formed the craters If the image has less than eight craters, then the system will choose the maximum number of craters on the image. If one deals with the image with many craters, he/she can choose any number that he/she wants to detect based on the original image, which has to be captured first prior to the detection. Although sometimes the system will detect the craters with wrong pairing patches (light patches connected to the wrong dark patch and vice versa) the Lander should understand that one of them might still be one of the hazards that have to be avoided during landing application. The radial distance is calculated for each of the pairs detected in green lines as in Figures (14) and (15) below and using the equation (23) and, the shortest distance between adjacent pairs (light and dark patch) is chosen as a preliminary result for further angle calculations.

After the light patches were connected to the dark patches with a minimum distance between them, then the system will calculate the minimum angle by inputting the direction of the sun (assuming that the authors know the sun's angle) and later comparing it with the pairing patches angle using the dot product (equation 24) as has been thoroughly described in the methodology section. By taking into account both techniques (minimum distance and angle), the authors can determine the best craters that they want on an image. The final craters detected will be denoted in yellow lines as shown in Figures (14) and (15) below.

Ideally, this algorithm will work on any images that have a clear pattern of light and dark patches and the authors do not even have to know the important parameters such as radius, gradient and etc of the craters. Unfortunately, this algorithm will work effectively on the image that has a clear pattern of these light and dark patches only. A crater with a clear pattern in this context will have clear features of light and dark patches that constitute one crater. To determine the accuracy of the algorithm created by the authors, the authors will provide a figure taken from two separate images tested by this algorithm. The final results have to be compared with the real image in order to determine which craters are true. In a real scenario, all craters detected will be considered hazards even though they are connected to the wrong pair of light or dark patches.

Comparing the results taken from Image 1to the original image (real image of the moon's surface), there are eleven valid craters or true craters after pre-processing (after erosion and noise reduction (*bwareaopen* function), labelled by BW label). Valid craters or true craters in this context means that the craters have a clear pattern of light and dark blobs regardless of the size and other features. Over these eleven valid craters, eight of them are successfully detected using this craters detection algorithm; this suggests that the authors have a 73% successful rate. The authors have tested this craters detection algorithm into two different optical images with different sun angles which are 10 degrees, > 10 degrees and with a noisy image as mentioned below. These evaluations below are to measure the accuracy of the algorithm based on the two images proposed and how robust the algorithm is. The determination of accuracy of the algorithm is based on how often craters are detected and the calculations are shown below:

#### **Based on the original Image 1 as can be illustrated in Figure (14) below**

#### **Sun direction: 10 degrees**

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

another binary image.

result for further angle calculations.

Centroid determination stage is completed by using the 'regionprops' function in Matlab and this is basically to compute the desired region properties or targeted properties. This step also attached the pairing patches together to form a complete crater for further analysis. Regionprops only takes the label matrix as an input. Therefore, the authors have to preprocess the image by labelling it first using 'bwlabel' function along with regionprops function later. In MATLAB, bwlabel is used to label the connected components in binary image and bwlabel also supports the 2-D inputs only. After labeling the entire target then only the authors can use them in a memory with regionprops. Regionprops takes this labelled component as an input and return them with an image of centroid determination labelled as a cross (\*) symbol as in Figure (13) above. For a smoother image, the authors apply another function in MATLAB called bwareaopen to remove all of the unnecessary objects or all connected components (objects) that are fewer than P pixels set, producing

The final stage of this craters detection algorithm is minimum radial distance and angle calculation. The authors have tested the algorithm with two different image conditions (different sun angles and background noises). According to the algorithm, the authors can select the number of craters to be detected by the system. The users can choose how many craters to be detected by the system based on the original image taken. As for this case, it will select the eight best craters that satisfy the minimum distance and angle conditions by assuming that the sun direction is known. An image of the landing site on the moon's surface has to be captured first and the amounts of craters needs to be detected are calculated manually. Based? on the original image, the authors have assumed the sun direction by identifying the shade and sunny pattern locations that formed the craters If the image has less than eight craters, then the system will choose the maximum number of craters on the image. If one deals with the image with many craters, he/she can choose any number that he/she wants to detect based on the original image, which has to be captured first prior to the detection. Although sometimes the system will detect the craters with wrong pairing patches (light patches connected to the wrong dark patch and vice versa) the Lander should understand that one of them might still be one of the hazards that have to be avoided during landing application. The radial distance is calculated for each of the pairs detected in green lines as in Figures (14) and (15) below and using the equation (23) and, the shortest distance between adjacent pairs (light and dark patch) is chosen as a preliminary

After the light patches were connected to the dark patches with a minimum distance between them, then the system will calculate the minimum angle by inputting the direction of the sun (assuming that the authors know the sun's angle) and later comparing it with the pairing patches angle using the dot product (equation 24) as has been thoroughly described in the methodology section. By taking into account both techniques (minimum distance and angle), the authors can determine the best craters that they want on an image. The final craters detected will be denoted in yellow lines as shown in Figures (14) and (15) below.

Manual Detection (number of craters after pre-processed): 11 Automatic Detection (number of craters detected): 8 8/11 x 100 = **73% accuracy** 

**Figure 14.** Craters detected from Image 1 in yellow line based on distance and angle measurements

## **Based on the original Image 2 as can be illustrated in Figure (15) below**

**Sun direction: > 10 degrees** 

Manual Detection (number of craters detected after pre-processed): 10 Automatic Detection (number of craters detected): 8 8/10 x 100 = **80% accuracy** 

**Figure 15.** Craters detected from Image 2 in yellow line based on distance and angle measurements

In Figure (15) above, yellow lines, which denote as craters are detected with a minimum distance and angle detection while green lines, which denote as craters are detected with a minimum distance only prior to the minimum angle detection. This angle detection will be a final stage in defining the craters based on the light and dark patch pattern (sometimes denoted as sunny and shady parts) and the final craters are those with yellow lines. By comparison, the accuracy of the algorithm based on these two images with different types of craters, angle (Sun) and lighting condition is said to be 77% and it is quite a satisfactorily accurate.

This accuracy factor can be improved if the authors know exactly the sun elevation angle since in this research; the authors just assumed the angle and the value is not really accurate. In a real application, this sun angle can be measured separately using the satellite, altimeter or radar prior to this detection process and the value will be more accurate. Besides, this algorithm will detect the craters that are above 0.0265 meters in image size (100 pixels). This can be vouched by using the '*bwareaopen'* function where it will remove the entire blobs pixel which is less than 100 pixels.
