Rectangular Natural Feature Recognition and Pose Measurement Method for Non-Cooperative Spacecraft

1. Introduction

The advancement of space technology has brought to the forefront a range of critical challenges, such as spacecraft maintenance and rescue in orbit [1], the removal of space debris [2], satellite rendezvous and docking [3], companion flight monitoring, and in-orbit capture [4]. These are urgent problems that demand effective solutions. Key to the successful execution of in-orbit missions is target recognition and pose measurement, which play pivotal roles [5]. Many of the targets in these missions are non-cooperative, lacking known structural, size, or motion information, and often without cooperative markers or communication capabilities with tracking spacecraft. Identifying and measuring non-cooperative targets poses a significant and challenging task [6,7].

Target detection can usually be divided into three stages in space [8]: long-distance (hundreds of meters to tens of kilometers), mid-distance (several kilometers to tens of meters) and close-distance (a few hundred meters to several meters). Generally, the first two stages usually use microwave radar and lidar to detect the approximate location information of the target. They have advantages in long-distance measurement, but the measurement accuracy is low and it is difficult to obtain attitude information [9,10]. For close-distance target detection, visual measurement equipment are usually used to provide the position and attitude information of targets. This article focuses on the pose measurement of non-cooperative spacecraft in short distance.

Some space programs and research departments have adopted a number of methods for close target detection. The Canadian Space Agency and NASA jointly developed triangulation LIDAR for the close-range measurement. But the point cloud matching database needs to be pre-trained [11,12]. Liu used flash laser radar to generate point cloud data, and proposed a pose tracking method based on the known satellite model [13]. Tzschichholz measured the non-cooperative targets by comparing the 3D point cloud features observed by the ToF camera with the actual 3D models [14]. Klionovska used the ToF camera to detect the non-cooperative targets combining target 3D model [15]. Gao et al. proposed an algorithm for measuring the position and attitude of targets using a monocular camera and a laser rangefinder [16]. The ranging accuracies of the pose measurement technique based on lidar and vision are usually respectively centimeter and millimeter. In order to achieve fine work at close range, visions are often chosen as the solution. We typically want to transfer images from on-orbit operations back to the ground for preservation and analysis [17]. Therefore, satellites are usually equipped with cameras (with different optimal imaging distances and observation directions), and no additional equipment is required. It is difficult for monocular vision to measure the pose of the target without prior information. As for the multi-vision, the object information obtained has more redundancy, which can enhance robustness and measurement accuracy of the system, while the cost price is higher and the operation speed is relatively slow. Binocular visions are widely concerned by researchers because of its balance between detection accuracy and operation speed.

There are some examples of binocular stereo vision used to detect non-cooperative targets. There are two main types of target detection methods. The first one is using deep neural network models to achieve target detection and pose estimation [18,19]. Daigo proposed the AI-assisted near-field pose estimation method of spacecraft whose shape and size are completely known. This method detected the keypoints of spacecraft with machine learning techniques and solve the pose by the Efficient Perspective-n-Point (EPnP) algorithm and particle swarm optimizer (PSO) method [20]. Different deep learning approaches have more or less the following problems. Pre-training is expensive and difficult to cover a variety of targets. It requires high computing power of the onboard computer. Many researches are focused on target detection, and there are few researches on target pose solving. The second method relies on image processing techniques to identify the salient features of the target, including but not limited to triangular brackets [21], rectangular solar panels, and circular engine nozzles [22]. Xu introduced a binocular measurement approach for assessing the interface ring and engine nozzle of rockets [23], while Peng devised a virtual stereo vision measurement method reliant on the interplay of triangular and circular features, albeit with the need for manual selection of interior points [24]. Additionally, Qiao proposed a relative pose measurement technique based on point features [25], but it necessitates prior knowledge about the target. Feature-based Target measurement based on feature detection is susceptible to environment and illumination. But it has the advantages of fast computing speed and low cost.

This paper takes the method of identifying the salient features of the target. In theory, larger target features tend to yield higher recognition success rates. Therefore, this paper focuses on researching the aircraft body or rectangular solar panel, which represents the most prominent area within the target. The surface shapes of the aircraft body and the solar panels are usually rectangles. A stereo vision detection method for rectangular features is studied in this paper. The detection of rectangular features in both space and ground applications typically involves processes such as edge line feature extraction, template matching, rectangle fitting, and minimum rectangle envelope extraction. Solar panels consist of numerous small rectangular photo receptors, resulting in a multitude of straight-line features. Extracting these straight-line features is time-consuming, and accurately discerning the target lines can be challenging. Template matching demands prior knowledge of the target, so it is not suitable for the non-cooperative spacecraft detection. Regarding rectangular fitting or minimal rectangular envelope extraction, the accuracy of these methods is limited due to the distortion of rectangular features in Cartesian space when projected onto the image plane [26,27,28]. The error is large in fitting the contour by means of rectangular envelope.

Through the subsequent analysis, it can be found that the projection of the rectangular feature on the image plane is approximately parallelogram. Therefore, we study the parallelogram fitting method of rectangular features and apply it to non-cooperative spacecraft detection. The algorithm has following advantages: (1) The algorithm has strong adaptability and can be used to detect rectangular features that most spacecrafts have; (2) The calculation efficiency of the algorithm is high, because the fine elements and interference elements inside the shape envelope will be gradually eliminated during operation, greatly reducing the amount of calculation; (3) The algorithm has high fitting accuracy. It combines the imaging characteristics of rectangular features. The major contributions of this paper are presented as follows: (1) The imaging characteristics of rectangles, which are the typical natural characters, in the image plane are analyzed. A parallelogram fitting method for rectangular features in the image plane is proposed; (2) The parallelogram fitting method is applied to the measurement of non-cooperative spacecraft. Natural features on non-cooperative targets such as spacecraft bodies and solar panels are detected and measured without adding additional cooperative markers.

This paper is organized as follows: Section 2 outlines the mission of detecting non-cooperative spacecraft and proposes the framework of this paper’s non-cooperative spacecraft measurement method. Section 3 demonstrates the properties of rectangular features when projected onto the image plane. In the fourth section, we introduce the parallelogram fitting algorithm for rectangular features in Cartesian space. Section five presents the pose solution algorithm for non-cooperative space targets. The sixth section describes the establishment of an experimental system to evaluate our proposed method. Finally, the last section provides a summary of the paper and draws conclusions.

2. Space On-Orbit Measurement Tasks and Version Measurement Algorithm Framework

2.1. Space On-Orbit Measurement Tasks

Figure 1 shows the process of detecting and approaching a non-cooperative spacecraft during on-orbit operation. The servicing satellite detects the target satellite from long distance like hundreds of meters to tens of kilometers or predict the target satellite’s orbit. The servicing satellite monitors the target at long-distance or mid-distance and reconstructs the target point clouds. In the whole space on-orbit service process, this paper focuses on the final tracking phase, that is, the Euler distance between the target and the service star is between 1m to 30m. The stereo camera is mounted on the servicing satellite. This camera is utilized to observe the non-cooperative target satellite and compute its position and orientation, which in turn is used to plan the movements of the control satellite.

In the space environment, lighting conditions are both harsh and subject to rapid changes. To ensure effective detection, it is prudent to choose the most prominent target on the satellite, as depicted in Figure 2. The solar panel, with its expansive surface area and conspicuous visibility when deployed, emerges as the optimal choice for visual detection. Notably, solar panels typically exhibit a rectangular plane surface, rendering them an ideal fit for the parallelogram fitting method introduced in this paper.

2.2. Overview of the On-Orbit Pose Measurement Method

The process of the pose measurement algorithm for non-cooperative targets is shown in Figure 3. It comprises the following steps:

Step 1: After capturing the image of the non-cooperative target using a stereo camera, the image distortion is corrected based on pre-calibrated camera parameters.

Step 2: In the initial detection phase, search for the satellite across the entire image. To enhance algorithm speed, the region of interest (ROI) can be defined based on previous measurement results in subsequent detections.

Step 3: Convert the image to the HLS format and apply median filtering.

Step 4: Set detection criteria based on the HLS values of each pixel, extract the target region, and convert it into a binary image.

Step 5: To minimize interference from silicon wafer gaps on the extracted image’s integrity and improve the extraction of the solar panel, perform morphological closing on the binary map.

Step 6: Search for closed contours in the image and fit these contours with minimum rectangular envelopes.

Step 7: Since the solar panel might have internal interference, it may be divided into multiple contours. Fit rectangles to each contour. Merge contour points based on the center point distance of each rectangle and the rectangle side length. Then, fit the minimum rectangle to the newly generated point set. Repeat this step until no further contour merging is possible.

Step 8: Select target contour points based on contour size, the area relationship between the contour and the rectangle, and the shape of the fitted rectangle.

Step 9: Utilize the parallelogram fitting algorithm proposed in this paper to fit the contour points and identify the corner points of the fitted quadrilateral.

Step 10: Solve the pose of the non-cooperative satellite in Cartesian space relative to the stereo vision system based on the corner points.

3. Planar Projection Properties of Rectangular Features

The majority of cameras operate on a center projection model, where features from Cartesian space are subject to distortion upon projection onto the camera’s image plane. The characteristics of rectangular features when projected onto the camera image plane are illustrated as follows. The projection of a rectangle from Cartesian space onto the image plane of the camera is depicted in Figure 4. The camera coordinate system ${\mathit{O}}_c$-${\mathit{X}}_c$${\mathit{Y}}_c$${\mathit{Z}}_c$ is taken as the world coordinate system, and the coordinate of Cartesian space rectangle ${\mathit{P}}_1$${\mathit{P}}_2$${\mathit{P}}_3$${\mathit{P}}_4$ in the camera coordinate system is (${\mathit{X}}_n$, ${\mathit{Y}}_n$, ${\mathit{Z}}_n$) (n = 1,2,3,4). The coordinates they project into the image plane ${\mathit{o}}_c$-${\mathit{x}}_c$${\mathit{y}}_c$ are ${\mathit{p}}_n$(${\mathit{x}}_n$, ${\mathit{y}}_n$) (n = 1,2,3,4). According to the principle of central projection:

$$\left\{\begin{array}{cc}{x}_n=f\frac{X_n}{Z_n}\hfill & \\ y_n=f\frac{Y_n}{Z_n}\hfill & \end{array}\right.$$

(1)

Where, f is the focal length of the camera.

Rectangular features in 3D space have the following two properties when projected to 2D image plane:

Property 1: Vertical segments in the same plane in Cartesian space are not necessarily vertical in the image plane.

prove: In the camera coordinate system, $\overrightarrow{{\mathit{P}}_1{\mathit{P}}_2}\perp \overrightarrow{{\mathit{P}}_1{\mathit{P}}_3}$, so:

$$\overrightarrow{{\mathit{P}}_1{\mathit{P}}_2}\cdot \overrightarrow{{\mathit{P}}_1{\mathit{P}}_3}=({\mathit{X}}_2-{\mathit{X}}_1)({\mathit{X}}_3-{\mathit{X}}_1)+({\mathit{Y}}_2-{\mathit{Y}}_1)({\mathit{Y}}_3-{\mathit{Y}}_1)+({\mathit{Z}}_2-{\mathit{Z}}_1)({\mathit{Z}}_3-{\mathit{Z}}_1)=0$$

(2)

In the image plane:

$$\begin{array}{cc}\hfill \overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}\cdot \overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}=& ({\mathit{x}}_2-{\mathit{x}}_1)({\mathit{x}}_3-{\mathit{x}}_1)+({\mathit{y}}_2-{\mathit{y}}_1)({\mathit{y}}_3-{\mathit{y}}_1)=\hfill \\ \hfill & {\mathit{f}}^2[(\frac{{\mathit{X}}_2}{{\mathit{Z}}_2}-\frac{{\mathit{X}}_1}{{\mathit{Z}}_1})(\frac{{\mathit{X}}_3}{{\mathit{Z}}_3}-\frac{{\mathit{X}}_1}{{\mathit{Z}}_1})+(\frac{{\mathit{Y}}_2}{{\mathit{Z}}_2}-\frac{{\mathit{Y}}_1}{{\mathit{Z}}_1})(\frac{{\mathit{Y}}_3}{{\mathit{Z}}_3}-\frac{{\mathit{Y}}_1}{{\mathit{Z}}_1})]\hfill \end{array}$$

(3)

In general, the distance between the target and the camera in the z-axis direction is much larger than the size of the solar panel, so ${\mathit{Z}}_1\approx {\mathit{Z}}_2\approx {\mathit{Z}}_3$, and

$$\begin{array}{cc}\hfill \overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}\cdot \overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}& \approx {\mathit{f}}^2[({\mathit{X}}_2-{\mathit{X}}_1)({\mathit{X}}_3-{\mathit{X}}_1)+({\mathit{Y}}_2-{\mathit{Y}}_1)({\mathit{Y}}_3-{\mathit{Y}}_1)]/{\mathit{Z}}_1\hfill \\ \hfill & \ne {\mathit{f}}^2\overrightarrow{{\mathit{P}}_1{\mathit{P}}_2}\cdot \overrightarrow{{\mathit{P}}_1{\mathit{P}}_3}/{\mathit{Z}}_1=0\hfill \end{array}$$

(4)

So $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}\cdot \overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}\ne 0$, which means $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}$ and $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}$ are not vertical. When ${\mathit{Z}}_1={\mathit{Z}}_2$ or ${\mathit{Z}}_1={\mathit{Z}}_3$, $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}$ or $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}$ is parallel to the plane of the image, and $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_2}$ and $\overrightarrow{{\mathit{p}}_1{\mathit{p}}_3}$ are almost perpendicular.

Property 2: When parallel line segments in Cartesian space are far away from the image plane and the length of line segments is short, their projections in the image plane are approximately parallel.

prove: In the camera coordinate system, $\overrightarrow{{\mathit{P}}_1{\mathit{P}}_2}//\overrightarrow{{\mathit{P}}_3{\mathit{P}}_4}$ and $|\overrightarrow{{\mathit{P}}_1{\mathit{P}}_2}|=|\overrightarrow{{\mathit{P}}_3{\mathit{P}}_4}|$, so:

$${\displaystyle \frac{{\mathit{X}}_2-{\mathit{X}}_1}{{\mathit{X}}_4-{\mathit{X}}_3}}={\displaystyle \frac{{\mathit{Y}}_2-{\mathit{Y}}_1}{{\mathit{Y}}_4-{\mathit{Y}}_3}}={\displaystyle \frac{{\mathit{Z}}_2-{\mathit{Z}}_1}{{\mathit{Z}}_4-{\mathit{Z}}_3}}=1$$

(5)

The following conditions exist for the size of solar panels in Cartesian space:${\mathit{X}}_i-{\mathit{X}}_j,{\mathit{Y}}_i-{\mathit{Y}}_j,{\mathit{Z}}_i-{\mathit{Z}}_j<<{\mathit{Z}}_n,{\mathit{Z}}_1\approx {\mathit{Z}}_2\approx {\mathit{Z}}_3\approx {\mathit{Z}}_4$. In the image plane, the slopes of ${\mathit{p}}_1{\mathit{p}}_2$ and ${\mathit{p}}_3{\mathit{p}}_4$ are ${\mathit{k}}_{p1}{\mathit{k}}_{p2}$ and ${\mathit{k}}_{p3}{\mathit{k}}_{p4}$ respectively:

$$\left\{\begin{array}{cc}\hfill {\mathit{k}}_{p1}{\mathit{k}}_{p2}=& ({\mathit{y}}_2-{\mathit{y}}_1)/({\mathit{x}}_2-{\mathit{x}}_1)=\mathit{f}({\mathit{Y}}_2/{\mathit{Z}}_2-{\mathit{Y}}_1/{\mathit{Z}}_1)/\left[\mathit{f}({\mathit{X}}_2/{\mathit{Z}}_2-{\mathit{X}}_1/{\mathit{Z}}_1)\right]\hfill \\ \hfill =& ({\mathit{Z}}_1{\mathit{Y}}_2-{\mathit{Z}}_2{\mathit{Y}}_1)/({\mathit{Z}}_1{\mathit{X}}_2-{\mathit{Z}}_2{\mathit{X}}_1)\approx ({\mathit{Y}}_2-{\mathit{Y}}_1)/({\mathit{X}}_2-{\mathit{X}}_1)\hfill \\ \hfill {\mathit{k}}_{p3}{\mathit{k}}_{p4}=& ({\mathit{y}}_4-{\mathit{y}}_3)/({\mathit{x}}_4-{\mathit{x}}_3)\approx ({\mathit{Y}}_4-{\mathit{Y}}_3)/({\mathit{X}}_4-{\mathit{X}}_3)=({\mathit{Y}}_2-{\mathit{Y}}_1)/({\mathit{X}}_2-{\mathit{X}}_1)\hfill \end{array}\right.$$

(6)

So ${\mathit{k}}_{p1p2}\approx {\mathit{k}}_{p3p4}$, which means ${\mathit{p}}_1{\mathit{p}}_2//{\mathit{p}}_3{\mathit{p}}_4$, similarly, ${\mathit{p}}_1{\mathit{p}}_3//{\mathit{p}}_2{\mathit{p}}_4$.

Indeed, parallel lines in Cartesian space will converge at infinity when projected onto the image plane. However, when these three-dimensional parallel lines occupy only a few pixels in the image plane and are distant from the camera plane, they can be approximated as nearly parallel in the image plane. In the on-orbit operation, when the target is close to the camera, the relative attitude error between the target and the camera is generally very small. The proof of property 2 still holds. Furthermore, the algorithm presented in this paper includes a correction step for the fitted parallel lines in the image plane. This correction aligns the fitted lines with the actual target contour rather than maintaining their parallel nature. This adjustment ensures that the fitted lines closely adhere to the contour points, enhancing the accuracy and effectiveness of the fitting process.

4. The Parallelogram Fitting Algorithm for Rectangular Features

4.1. Parallelogram Fitting Algorithm Framework

The comprehensive procedure of the parallelogram fitting algorithm, as proposed in this paper, is delineated in Figure 5. The algorithm can be fed with input data in the form of a closed contour, a point set, or a set of contours. Here is a step-by-step breakdown of the algorithm: Step 1: Linear fitting of contour points, resulting in the centerline $\mathit{L}$ of contour points based on the input point set. Step 2: Division of the point set into ${\mathit{N}}_1$ and ${\mathit{N}}_2$ on both sides of $\mathit{L}$ to facilitate subsequent line fitting. Step 3: The point set for fitting parallel lines is selected and updated as ${\mathit{N}}_1$ and ${\mathit{N}}_2$ based on the distance between $\mathit{L}$ and each point in these sets. Step 4: Use ${\mathit{N}}_1$ and ${\mathit{N}}_2$ to fit parallel lines ${\mathit{l}}_1$ and ${\mathit{l}}_2$, updating ${\mathit{N}}_1$ and ${\mathit{N}}_2$ according to the distance between each point and the parallel lines. Step 5: Solve for the two groups of points at the edges of ${\mathit{N}}_1$ and ${\mathit{N}}_2$. These two sets of points are interconnected to form a starting line for locating another set of points. Step 6: Preliminarily confirm the point set used to fit another group of parallel lines and subsequently fit those parallel lines. Step 7: Iteratively confirm the point sets used to fit the two groups of parallel lines and perform the fitting. Four sets of points and two sets of parallel lines are cyclically updated based on the distance between each point in the point sets and the lines. The distance threshold reduces as the iteration count increases to exclude points that significantly interfere with the parallel line fitting. When the elements of each point set no longer shrink after two iterations, the final point sets for fitting the four edges are derived. Step 8: To enhance the fitting precision, four lines are fitted separately using the four point sets. The intersection of these four lines yields the final fitting quadrilateral, with its four corners precisely determined. This algorithm provides a systematic approach for accurately fitting parallelograms, ensuring robust results for various input scenarios.

4.2. Line Fitting of Contour Points

In order to fit a line based on the input contour point set, it’s necessary to first calculate the coordinates of the contour’s center point. The coordinates of the center point $({\mathit{u}}_0,{\mathit{v}}_0)$ of the contour points $\mathit{M}$ can be computed using the following formula:

$${\mathit{u}}_0=\frac{1}{\mathit{n}}\sum _{i=1}^n{\mathit{u}}_i,{\mathit{v}}_0=\frac{1}{\mathit{n}}\sum _{i=1}^n{\mathit{v}}_i$$

(7)

$\mathit{L}$ fitted by contour points satisfies the following equation:

$$\left\{\begin{array}{cc}au+bv+c=0\hfill & \\ a^2+b^2=1\hfill & \end{array}\right.$$

(8)

The distance between contour points and $\mathit{L}$ are: $a{u}_i+b{v}_i+c$. To solve the center line $\mathit{L}$ of the contour point is to solve a, b and c, so that the minimum value of the following formula can be obtained:

$$\mathit{f}=\sum {(a{u}_i+b{v}_i+c)}^2$$

(9)

Solve the above equation with the Lagrange multiplier method:

$$\left\{\begin{array}{cc}\hfill & \mathit{f}=\sum {(a{u}_i+b{v}_i+c)}^2-\lambda (a^2+b^2-1)\hfill \\ \hfill & \frac{\partial \mathit{f}}{\partial \mathit{a}}=0\frac{\partial \mathit{f}}{\partial \mathit{b}}=0\frac{\partial \mathit{f}}{\partial \mathit{c}}=0\frac{\partial \mathit{f}}{\partial \lambda }=0\hfill \end{array}\right.$$

(10)

Simplified above formula:

$$\left[\begin{array}{cc}\sum {({\mathit{x}}_i-\overline{x})}^2& \sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y})\\ \sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y})& \sum {({\mathit{y}}_i-\overline{y})}^2\end{array}\right]\left[\begin{array}{c}\mathit{a}\\ \mathit{b}\end{array}\right]=\lambda \left[\begin{array}{c}\mathit{a}\\ \mathit{b}\end{array}\right]$$

(11)

It is obtained by solving for the eigenvalue:

$$\lambda ={\displaystyle \frac{\left[\sum {({\mathit{x}}_i-\overline{x})}^2+\sum {({\mathit{y}}_i-\overline{y})}^2-\sqrt{{(\sum {({\mathit{x}}_i-\overline{x})}^2-\sum {({\mathit{y}}_i-\overline{y})}^2)}^2+4{(\sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y}))}^2}\right]}{2}}$$

(12)

So:

$$\left\{\begin{array}{cc}\hfill & \left(\begin{array}{c}\mathit{a}\\ \mathit{b}\end{array}\right)={\displaystyle \frac{\left(\begin{array}{c}\sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y})\\ \lambda -\sum {({\mathit{x}}_i-\overline{x})}^2\end{array}\right)}{\sqrt{{(\sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y}))}^2+{(\lambda -\sum {({\mathit{x}}_i-\overline{x})}^2)}^2}}}\hfill \\ \hfill & c=-a\overline{x}-b\overline{y}={\displaystyle \frac{-(\sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y}))\overline{x}-(\lambda -\sum {({\mathit{x}}_i-\overline{x})}^2)\overline{y}}{\sqrt{{(\sum ({\mathit{x}}_i-\overline{x})({\mathit{y}}_i-\overline{y}))}^2+{(\lambda -\sum {({\mathit{x}}_i-\overline{x})}^2)}^2}}}\hfill \end{array}\right.$$

(13)

4.3. Parallelogram Fitting

In the preceding section, we successfully fitted the centerline L. We then proceeded to segment the set of contour points into two distinct sets: one on each side of L. By considering the distance between each point and the centerline, we initially identified two point sets, ${\mathit{N}}_1$ and ${\mathit{N}}_2$, which were earmarked for fitting a pair of parallel lines, ${\mathit{l}}_1$ and ${\mathit{l}}_2$. Following this, we selected edge points from ${\mathit{N}}_1$ and ${\mathit{N}}_2$ and used them to determine two additional point sets, ${\mathit{N}}_3$ and ${\mathit{N}}_4$, for fitting yet another pair of parallel lines within the parallelogram. These newly determined point sets, ${\mathit{N}}_3$ and ${\mathit{N}}_4$, were then employed to initiate the fitting of another set of parallel lines, ${\mathit{l}}_3$ and ${\mathit{l}}_4$.

4.3.1. Solving of the points sets on both sides of the center line

Firstly, the points in the contour point set are assigned to the point sets on both sides of the line according to their position relation with L. Replace ${\mathit{u}}_i$ in each point m(${\mathit{u}}_i$, ${\mathit{v}}_i$) in the contour point set with line L: $au+bv+c=0$ to solve ${\mathit{v}}_i^{\prime }$. The relationship between each point in the contour point set and L can be summarized into the following four categories, as shown in Figure 6. As shown in Figure 6 (1) and (2), the slope of L is $|-a/b|<1$. If ${\mathit{v}}_i$ < ${\mathit{v}}_i^{\prime }$, $\mathit{m}\in {\mathit{N}}_1$, otherwise $\mathit{m}\in {\mathit{N}}_2$. As shown in Figure 6 (3), the slope of L is $-a/b⩽-1$. If ${\mathit{v}}_i$ < ${\mathit{v}}_i^{\prime }$, $\mathit{m}\in {\mathit{N}}_1$, otherwise $\mathit{m}\in {\mathit{N}}_2$. As shown in Figure 5 (4), the slope of L is $-a/b⩾1$. If ${\mathit{v}}_i$ > ${\mathit{v}}_i^{\prime }$, $\mathit{m}\in {\mathit{N}}_1$, otherwise $\mathit{m}\in {\mathit{N}}_2$.

4.3.2. Initially fitting of a set of parallel lines

A set of parallel lines are initially fitted with ${\mathit{N}}_1$ and ${\mathit{N}}_2$. Solve the farthest distance ${\mathit{d}}_1$ and ${\mathit{d}}_2$ of the points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$ to $\mathit{L}$ with the follow formula:

$$\left\{\begin{array}{cc}\hfill & d_1=\underset{i=1}{\overset{n_1}{max}}|a{u}_i+b{v}_i+c|\hfill \\ \hfill & d_2=\underset{j=1}{\overset{n_2}{max}}|a{u}_j+b{v}_j+c|\hfill \end{array}\right.$$

(14)

Where, ${\mathit{n}}_1$ and ${\mathit{n}}_2$ are the number of elements of ${\mathit{N}}_1$ and ${\mathit{N}}_2$.

We calculate the distances between the points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$ with respect to the best fitting line L. These distances are classified based on their deviation from the values ${\mathit{d}}_1$ and ${\mathit{d}}_2$. In theory, the set of edge points should be the farthest from L. Therefore, the set of points exhibiting the greatest deviation from the line should be used as the basis for fitting the parallel lines. However, if the number of elements in a point set is too small, it suggests that the set mainly comprises noise points. In such cases, we select the point set with the largest number of elements for fitting the parallel lines, which is then updated as ${\mathit{N}}_1$ and ${\mathit{N}}_2$.

With ${\mathit{N}}_1$ and ${\mathit{N}}_2$ at hand, we proceed to fit a set of parallel lines. We utilize Equations (9) through (13) to fit the centerlines ${\mathit{l}}_1$ and ${\mathit{l}}_2$ for ${\mathit{N}}_1$ and ${\mathit{N}}_2$, with respective slopes denoted as ${\mathit{k}}_1$ and ${\mathit{k}}_2$. The inclination angles of these lines are represented by ${\theta }_1$ and ${\theta }_2$, respectively.

$${\theta }_1=arctan{\mathit{k}}_1{\theta }_2=arctan{\mathit{k}}_2$$

(15)

To ensure that the slopes are equal, set the slope angle of parallel lines $\theta =({\theta }_1+{\theta }_2)/2$. Then the ${\mathit{l}}_1$ and ${\mathit{l}}_2$ are:

$$\left\{\begin{array}{cc}sin\theta \mathit{u}-cos\theta \mathit{v}+c_1=a_{1}u+=b_{1}v+c_1=0\hfill & \hfill c_1=-a_1{\overline{u}}_{N_1}-b_1{\overline{v}}_{N_1}\\ sin\theta \mathit{u}-cos\theta \mathit{v}+c_2=a_{1}u+=b_{1}v+c_2=0\hfill & \hfill c_2=-a_1{\overline{u}}_{N_2}-b_1{\overline{v}}_{N_2}\end{array}\right.$$

(16)

Calculate the distance between points in ${\mathit{N}}_1$ and ${\mathit{l}}_1$ and the distance between points in ${\mathit{N}}_2$ and ${\mathit{l}}_2$. The points whose distance is less than the threshold are retained and reconstructed into ${\mathit{N}}_1$ and ${\mathit{N}}_2$.

4.3.3. Initially fitting another set of parallel lines

For fitting parallel lines, it is essential to partition the point set used for this purpose. Based on the ${\mathit{N}}_1$ and ${\mathit{N}}_2$ sets obtained in the previous section, we can determine a line using the edge points of the point set, and then search for the points near this line to form the set required for fitting parallel lines. The parallel lines, ${\mathit{l}}_1$ and ${\mathit{l}}_2$, fitted in the previous section exhibit four possible scenarios, as illustrated in Figure 7.

When $|a_1|<|b_1|$, the target edge points are chosen based on the u coordinates of points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$, as demonstrated in Figure 7(1) and (2). The points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$ with the smallest u coordinates correspond to the leftmost points, denoted as ${\mathit{p}}_1$ and ${\mathit{p}}_3$, while those with the largest u coordinates correspond to the rightmost points, denoted as ${\mathit{p}}_2$ and ${\mathit{p}}_4$. Conversely, when $|a_1|>|b_1|$, the selection of target edge points is based on the v coordinates of points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$, as illustrated in Figure 7(3) and (4). In this case, the points in ${\mathit{N}}_1$ and ${\mathit{N}}_2$ with the smallest v coordinates represent the topmost points, designated as ${\mathit{p}}_1$ and ${\mathit{p}}_3$, while those with the largest v coordinates correspond to the bottommost points, denoted as ${\mathit{p}}_2$ and ${\mathit{p}}_4$. Connect ${\mathit{p}}_1{\mathit{p}}_3$ and ${\mathit{p}}_2{\mathit{p}}_4$ to form two lines. Take ${\mathit{p}}_1{\mathit{p}}_3$ as an example, its equation is:

$$p_1{p}_3:{\displaystyle \frac{{\displaystyle \frac{v_1-v_3}{u_1-u_3}}}{\sqrt{1+{\left({\displaystyle \frac{v_1-v_3}{u_1-u_3}}\right)}^2}}}u-{\displaystyle \frac{1}{\sqrt{1+{\left({\displaystyle \frac{v_1-v_3}{u_1-u_3}}\right)}^2}}}v+{\displaystyle \frac{v_1-{\displaystyle \frac{v_1-v_3}{u_1-u_3}}{u}_1}{\sqrt{1+{\left({\displaystyle \frac{v_1-v_3}{u_1-u_3}}\right)}^2}}}=0$$

(17)

To proceed, we calculate the distance between each point in the contour point set $\mathit{M}$ and the line segments ${\mathit{p}}_1{\mathit{p}}_3$ and ${\mathit{p}}_2{\mathit{p}}_4$. Points whose distance is less than a predetermined threshold are selected to form the point sets ${\mathit{N}}_3$ and ${\mathit{N}}_4$. Referring to equations (16) and (17), we then utilize ${\mathit{N}}_3$ and ${\mathit{N}}_4$ to initiate the fitting of another set of parallel lines, ${\mathit{l}}_3$ represented as $a_{2}u+b_{2}v+c_3=0$, and ${\mathit{l}}_4$ represented as $a_{2}u+b_{2}v+c_4=0$.

4.3.4. Parallelogram fitting

Following the steps described above, we have successfully obtained two sets of parallel lines, ${\mathit{l}}_1$, ${\mathit{l}}_2$, and ${\mathit{l}}_3$, ${\mathit{l}}_4$, along with four point sets ${\mathit{N}}_i$ (i = 1, 2, 3, 4). During the iterative process, points in ${\mathit{N}}_i$ that are beyond the threshold distance from ${\mathit{l}}_1$ are iteratively removed, and ${\mathit{l}}_1$ is re-fitted after each update of point set ${\mathit{N}}_i$. The distance threshold between points and the parallel lines decreases with each iteration. When, in both cycles, the number of elements in ${\mathit{N}}_i$ falls below the threshold, it indicates that each point set has essentially reached its optimal configuration, and the iterative process is terminated.

Rectangles in Cartesian space only approximate parallelograms in the image plane. To further enhance the fitting accuracy, we use equations (9) through (13) to fit their optimal fitting lines, ${\mathit{l}}_i$, respectively. At this point, ${\mathit{l}}_1$ and ${\mathit{l}}_2$, as well as ${\mathit{l}}_3$ and ${\mathit{l}}_4$, need not remain parallel. These four lines intersect to form the final fitted quadrilateral. We employ the following formula to calculate the u and v coordinates of the four corners of the quadrilateral:

$$\left\{\begin{array}{cc}\hfill & u=(b_i{c}_j-b_j{c}_i)/(a_i{b}_j-a_j{b}_i)\hfill \\ \hfill & v=(a_i{c}_j-a_j{c}_i)/(a_j{b}_i-a_i{b}_j)\hfill \end{array}\right.$$

(18)

5. Pose Solution Method of Non-Cooperative Target

As shown in Figure 8, stereo vision is used to detect the corner points of the satellite. The coordinates of the satellite solar panel corners extracted by the left and right cameras in the pixel plane are respectively ${\mathit{p}}_{li}=[u_{li},v_{li]},p_{ri}=[u_{ri},v_{ri}]$ (i = 1, 2...12).

The left camera coordinate system is the world coordinate system. The position $P_i$ of each corner point of the satellite solar panel is [29]:

$${\left[X_i{Y}_i{Z}_i\right]}^T={\left({\mathit{P}}^T\mathit{P}\right)}^{-1}{\mathit{P}}^T\mathit{Q}$$

(19)

where.

$$\left\{\begin{array}{cc}\hfill & \mathit{P}=\left[\begin{array}{ccc}{u}_{li}{m}_{31}^l-m_{11}^l& u_{li}{m}_{32}^l-m_{12}^l& u_{li}{m}_{33}^l-m_{13}^l\\ v_{li}{m}_{31}^l-m_{21}^l& v_{li}{m}_{32}^l-m_{22}^l& v_{li}{m}_{33}^l-m_{23}^l\\ u_{ri}{m}_{31}^r-m_{11}^r& u_{ri}{m}_{32}^r-m_{12}^r& u_{ri}{m}_{33}^r-m_{13}^r\\ v_{ri}{m}_{31}^r-m_{21}^r& v_{ri}{m}_{32}^r-m_{22}^r& v_{ri}{m}_{33}^r-m_{23}^r\end{array}\right]\hfill \\ \hfill & \mathit{Q}=\left[\begin{array}{cccc}{m}_{14}^l-u_{li}& m_{34}^l{m}_{24}^l-v_{li}& m_{34}^l{m}_{14}^r-u_{ri}& m_{34}^r{m}_{24}^r-u_{ri}{m}_{34}^r\end{array}\right]\hfill \end{array}\right.$$

(20)

$m_{ij}^l$ and $m_{ij}^r$ represent the transformation matrix parameters that relate the pixel plane coordinate system to the world coordinate system for the left and right cameras, respectively. To enhance the accuracy of stereo vision measurements, a procedure involving plane fitting to the detected solar panel corners is conducted, where the projected corner points onto the fitted plane serve as the corrected corner points. The specific steps are as follows:

The least square method is used to fit the solar panel plane. In the left camera coordinate system, the depth information $Z_i$ of each corner point is greater than 0, and the fitting plane will not be parallel to the Z-axis of the left camera coordinate system (Because in parallel, the solar panel is a line in the image and cannot be measured). Therefore, the plane equation $\mathit{M}$ can be set as: $ax+by+c-z=0$. For the fitted plane, the sum of distances from each point to the plane is minimum, that is:

$$\left\{\begin{array}{cc}\hfill & mins=\sum _{i=1}^{12}{\left(a{X}_i+b{Y}_i+c-Z_i\right)}^2\hfill \\ \hfill & \frac{\partial s}{\partial a}=\frac{\partial s}{\partial b}=\frac{\partial s}{\partial c}=0\hfill \end{array}\right.$$

(21)

So,

$$\left[\begin{array}{c}a\\ b\\ c\end{array}\right]={\left[\begin{array}{ccc}\sum X_i^2& \sum X_i{Y}_i& \sum X_i\\ \sum X_i{Y}_i& \sum Y_i^2& \sum Y_i\\ \sum X_i& \sum Y_i& 12\end{array}\right]}^{-1}\left[\begin{array}{c}\sum X_i{Z}_i\\ \sum Y_i{Z}_i\\ \sum Z_i\end{array}\right]$$

(22)

The equation of the plane M is:

$$\left\{\begin{array}{cc}\hfill & \frac{a}{m}x+\frac{b}{m}y-\frac{1}{m}z+\frac{c}{m}=Ax+By+Cz+D=0\hfill \\ \hfill & m=\sqrt{a^2+b^2+1}\hfill \end{array}\right.$$

(23)

The projection coordinates $P_{Mi}(X_{Mi},Y_{Mi},Z_{Mi})$ of each corner point of the satellite on the plane M are:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{X}_{Mi}=X_i-At\\ Y_{Mi}=Y_i-Bt\\ Z_{Mi}=Z_i-Ct\end{array}\right.\hfill \\ \hfill & t=\frac{A{X}_i+B{Y}_i+C{Z}_i+D}{A^2+B^2+C^2}\hfill \end{array}$$

(24)

As shown in Figure 9, with the center of the solar panel as the origin of the coordinate system, the translation matrix of the solar panel relative to the camera coordinate system is $\mathit{T}=(P_{M3}+P_{M4}+P_{M9}+P_{M10})/4$. The solar panel coordinate system Z-axis in the camera coordinate system is: ${[r_{M13},r_{M23},r_{M33}]}^T={[A,B,C]}^T$. The X-axis in the camera coordinate system is:

$$\left[\begin{array}{c}{r}_{M11}\\ r_{M12}\\ r_{M13}\end{array}\right]=\frac{\sum _{i=1}^6{P}_{M2i}-\sum _{i=1}^6{P}_{M(2i-1)}}{6\left|\right|\sum _{i=1}^6{P}_{M2i}-\sum _{i=1}^6{P}_{M(2i-1)}\left|\right|}$$

(25)

So the Y-axis in the camera coordinate system is:

$$\left[\begin{array}{c}{r}_{M12}\\ r_{M22}\\ r_{M32}\end{array}\right]=\left[\begin{array}{c}{r}_{M13}\\ r_{M23}\\ r_{M33}\end{array}\right]\times \left[\begin{array}{c}{r}_{M11}\\ r_{M21}\\ r_{M31}\end{array}\right]$$

(26)

6. Experimental Verification

6.1. Satellite Natural Feature Recognition Experiment

To evaluate the success rate and the calculation speed of the satellite detection algorithm proposed in this paper, 100 times of tests which the distance and attitude between stereo vision and satellite varied a lot were conducted. The experimental system operated on a PC equipped with an Intel(R) i5-8500 CPU (3.0GHz) and 8.0GB of RAM. As an example, considering a monocular camera, the processing results of each step are illustrated in Figure 10. These steps successfully eliminate interference contours and points inside and outside the solar panel, while accurately extracting edge and corner points of the solar panel.

Figure 11 displays typical recognition results from these tests. It is evident that across different observation distances and angles, the algorithm presented in this paper consistently and accurately extracts the corners and contours of the satellite’s solar panels. The experiments resulted in a remarkable 98.5% success rate in satellite detection, with an average detection time of only 0.14 seconds. Instances of detection failure were primarily attributed to scenarios where the solar panel plane was nearly parallel to the camera’s optical axis, causing the solar panel to appear as a line in the image.

6.2. Measurement Accuracy Test of Satellite Pose

To assess the measurement accuracy of the algorithm proposed in this paper for determining the satellite’s pose, we established an experimental platform, as depicted in Figure 12. This platform featured air feet beneath both the target satellite and the service satellite, enabling them to move freely on the air platform. A stereo vision system was mounted on the service satellite, capturing images and transmitting them to a central controller that executed the algorithm presented in this paper. To ensure precise measurements, a laser tracker and calibration board were employed to calibrate the relative pose between the stereo vision system and the laser tracker, as well as to measure the pose between the satellite and the camera. The satellite’s poses were measured using stereo vision. These measurements were then compared with the data obtained from the laser tracker to evaluate the measurement accuracy of the algorithm proposed in this paper. The service satellite was equipped with an LIDAR, which can synchronously measure the relative position between the target satellite and the service satellite.

A comprehensive set of 20 experiments was conducted, and the resulting pose measurement errors are depicted in Figure 13. Figure 13(a) provides insight into the position measurement errors and the root mean square (MSE) of the position errors in the x, y, and z directions of the satellite, while Figure 13(b) displays the attitude measurement errors and the root mean square of the attitude errors in the x, y, and z directions of the satellite. The results indicate that the measurements of pose errors remain consistently stable. The position measurement errors in the x and y directions are less than 7mm, while the position measurement error in the z direction is less than 12mm. Furthermore, the three-dimensional attitude errors are all under 1.2°, demonstrating a high level of accuracy in attitude measurement. The average measurement errors for satellite pose are summarized in Table 1, where the average position measurement error is 10.412mm and the average attitude measurement error is 0.661°. These results underscore the precision and reliability of the proposed algorithm for satellite pose measurement.

The performance of various non-cooperative spacecraft detection methods can be evaluated through experimental testing and existing literature data. The detection results are shown in Table 2. In order to measure the detection success rate and detection speed of various methods, 100 experiments were carried out. As for this paper’s method, the successful rate of satellite detection was 98.5%, with an average detection time of only 0.14 seconds. In comparison to existing methods, including stereo vision-based outer rectangle detection, edge line detection, and LiDAR-based detection, the proposed method stands out. The success rate of stereo vision-based detection for solar panel edge was notably lower at only 67%. This reduced success rate can be attributed to the presence of numerous straight-line features in the image due to the solar panel’s composition of small rectangles, making it challenging to extract the target straight lines accurately. Overall, the method presented in this paper demonstrates superior performance in terms of speed and successful rate, surpassing all other methods except the stereo vision-based outer rectangle detection in the image plane, making it a compelling choice for satellite detection applications. As for the method of detecting satellites with deep learning in [20], the success rate was 96%. Because the algorithms are equipped with different hardware, the running time of the algorithms was not compared. Twenty experiments were carried out to test the measurement accuracy of various algorithms from different relative distance and relative attitude measurements of satellites. The test results show that the proposed method is highly accurate. External rectangle detection methods exhibit generally low accuracy, with fitting errors increasing as the deflection angle of the rectangular feature grows. Given the presence of numerous linear features on satellites, these features can significantly interfere with the extraction of the target line, resulting in a poor success rate and low accuracy in satellite edge line detection. LiDAR, for example, installed in this experimental system shown inaa Figure 12, typically achieves a detection accuracy of approximately 30mm, primarily capable of determining the target’s position. However, it often struggles with determining the attitude of the target and often necessitates the cooperation of visual sensors for this purpose. Methods like deep learning require substantial prior information about satellites, making them costly in terms of resource requirements. Compared with [24]’s method of detecting satellite pose with stereo vision, the detection accuracy of this paper is also significantly higher. In contrast, the parallelogram fitting method for rectangular features of non-cooperative spacecraft in the image plane, as proposed in this paper, exhibits exceptional performance. It excels in terms of detection accuracy, speed, and algorithm robustness, making it a promising choice for non-cooperative spacecraft feature detection, especially when other methods face limitations.

7. Conclusion

This paper introduces a novel method for recognizing rectangular natural features and measuring the pose of non-cooperative spacecraft, addressing the challenging task of detecting the poses of such spacecraft in space. The method is applicable to various objects on the spacecraft, including the spacecraft body and solar panels. It employs a parallelogram fitting approach for detecting non-cooperative targets and measuring their poses. Image features are extracted using stereo vision, and the projection of rectangular features onto the image plane is then detected and fitted into parallelograms. An experimental system was established to conduct pose measurement experiments at different distances. The proposed algorithm achieved an average position measurement error of 10.412mm and an average attitude measurement error of 0.661°. These results demonstrate the method’s suitability for applications of in-orbit close-range target measurement. This paper’s method offers valuable insights and reference points for the development of spatial object detection algorithms and the recognition of rectangular features in Cartesian space across various working scenarios. Future research will focus on satellite detection under different illumination conditions and visual servo control.