A new polarization direction measurement via local Radon transform and error correction

Vectorial optical field-based spatially polarization modulation has been widely studied for polarization measurement due to its simple system structure. In this system, the polarization information is encoded in the irradiance image, and polarization measurement can be realized by image processing. The classical image processing methods could not meet the increasing demand of practical applications due to their poor computational efficiency. To address this issue, a new image processing method, combining the rapidity of local radon transform (LRT) and the precision of error correction (EC), was proposed in this paper. Firstly, the polarization direction of the light was coarsely estimated from pixels on several circles. Then, the LRT of the input image was completed while the coarsely estimated direction was the center angle for LRT. Finally, the EC was conducted to get the accurate direction depending on the quantitative link between the error of the coarse estimation and the correlation between the LRTs. Experiments on synthetic and real data demonstrate that, compared to the other state-of-the-art methods, our proposed algorithm is more robust and less time-consuming.

spatial modulated polarimetry methods are highly dependent on spatial modulation devices. They were difficult to be deployed in the realistic scenarios.
To address the issue of the classical spatial modulated polarimetry methods, vectorial optical field-based spatially polarization modulated polarimetry methods have been proposed recently [6,18,[21][22][23]. In the strategy, the polarization information is recorded by the modulated intensity pattern of the input light. In this way, when the pattern of the light is captured by a camera, the polarization measurement can be transformed into the problem of analyzing the irradiance image. In [22,23], when a zero-order vortex quarter-wave retarder was used as a space variant birefringence device to achieve spatial modulation for all polarization components, a normalized least-squares method and a hybrid gradient descent algorithm were proposed, respectively, to calculate the polarization state from the irradiance images. Researchers have found that, when the input light was analyzed by an azimuthally (a radially) spatial modulator, the irradiance image has hourglass-shaped gray distribution. In other words, the darkest line of the irradiance image is parallel (or perpendicular) to the polarization direction [21]. Consequently, the polarization direction of the input light can be captured by extracting the darkest line from the image. In [6,10], the global Radon transform (GRT) was adopted. Gao and Lei [18] also chosen GRT to get the intensity modulation curve from which the four Stokes parameters of the input light can be measured. Lei and Liu [21] compared the accuracy and cost time of different image processing algorithms such as interesting area detection (IAD), local correlation (LC), GRT and so on. They found that the precision of the IAD was low, and the Radon transform was quite sensitive to image noise. The LC had more stable and higher accuracy, but it was time-consuming like IAD and Radon transform.
To measure the polarization direction more robust and faster, a novel method is presented in this paper. Motivated by GRT and LC, the method contains three stages: coarse estimation, local Radon transform (LRT), and error compensation (EC). At the first, the coarse direction is estimated based on threshold segmentation. Then, LRT is performed in a local angle range while the coarse estimated direction is taken as the center angle. In the end, the accurate darkest direction (parallel or perpendicular to the polarization direction) is gotten by EC, which establishes a quantitative link between the error of coarse estimation and the correlation between LRTs. The advantages of our algorithm are fourfold.
Firstly, the proposed method is robust to noise owing to the gray integral operation in LRT.
Secondly, the utilization of EC makes the proposed method highly precise. Thirdly, different to the GRT with small angle interval [6,10,18], LRT is only need to be computed in a local angle domain with large angle interval. It is therefore more computationally efficient.
Finally, since most of the processing can be done by looking up tables generated offline, our algorithm is suitable for real-time task for its high speed.
The outline of this article is as follows. In Sect. 2, the coarse estimation, LRT, and EC are introduced, followed by a flowchart to summarize our method. Section 3 shows several experimental results. Finally, some conclusions are drawn in Sect. 4.

Method
It has been verified that, when the input light was analyzed by an azimuthally (or a radially spatial modulator, the hourglass-shaped intensity pattern of the modulated light satisfies Malus's law [21]. In other words, the gray distribution of the irradiance image, as shown in Fig. 1, is directly proportional to the square of the cosine of the angle between the azimuthal angle and the darkest direction. The darkest direction, which is parallel (or perpendicular) to polarization direction, has the minimum radial integral value of the image. To capture the darkest direction accurately and quickly, our method include three stages: coarse estimation, LRT and EC. They are introduced as follow.

Coarse estimation
In our algorithm, the darkest direction is first coarsely estimated based on threshold segmentation. To reduce the computational complexity, threshold segmentation is processed on the pixels on the circles with certain radiuses rather than all the pixels in the image. Given a set of radiuses (e. g, r 1 , r 2 , r 3 , . . . , r N ), the pixels on the circles with different radiuses are collected. Then, the pixels are divided into two parts (i.e., bright area and dark area) based on a predefined threshold T . The average azimuthal angle of pixels in the dark area, denoted by θ c , is treated as the coarse darkest direction, i.e., where I(r, θ) is the gray value of the pixel with the coordinate (r, θ).

Local Radon transform
In this stage, Radon transform [25] is adopted to compute the integral of an image along specified directions. Suppose that f is a 2-D function, the integral of f along the radial line l(θ i ) = x, y : x sin θ i − y cos θ i = 0 is given by For digital images, Eq. (2) can be transferred as In Eq. (3), I(x, y) is the gray of the pixel with the rectangular coordinate (x, y) . W (ξ ) , the weight of the pixel (x, y) for integration along l(θ i ) , can be obtained by d is the distance threshold to determine whether the pixel (x, y) is on the line l(θ i ). Obviously, GRT needs to compute the integral of the image along radial lines orientated from 0 • to 180 • . Moreover, to have the accurate result, the angle interval that the GRT adopts should be as small as possible. Different from GRT, LRT only needs to capture the integral of the image in a local angle range, in which, the coarse darkest direction (i.e., θ c ) is taken as the center angle. For example, assuming the angle range and angle interval for LRT are ±θ T and θ s , LRT is gotten while the radial integral values are arranged in azimuth order.
As illustrated in Fig. 1, the actual darkest direction of the irradiance image is 25 • . As the image is disturbed by Gaussian white noise ( µ = σ 2 = 0.01 ), the darkest direction calculated by coarse estimation is 25.06 • (the white solid line in Fig. 1), LRT is composed of the normalized integral values of the image along the radial lines (the white dotted line in Fig. 1) counterclockwise oriented from 145.06 • to 85.06 • . Here, θ T is set to be 60 • .

Error correction
Theoretically, the darkest direction has the minimum value in LRT. It is regrettable that, the radial integral value of the image is always disturbed by the noise. For instance, the LRT of the image (shown in Fig. 1) is displayed in Fig. 2. The actual darkest direction of the image is 25 • , yet the direction that has the minimum value in LRT is 25.6 • . Apparently, the direction with the minimum value is not the actual darkest direction under the noise. To address this issue, EC is developed to explore the error of coarse estimation.
Assuming we have two modulate irradiance images ( Im 1 and Im 2 ) with hourglass-shaped gray distribution, and the darkest directions of two images are θ d1 and θ d2 , respectively, G 1 (θ d1 − θ a ) and G 2 (θ d2 − θ a ) has the best correlation. That is, and G 2 (θ d2 − θ a ) denote the LRTs of Im 1 and Im 2 while θ d1 − θ a and θ d 2 − θ a are the centers of the local angle ranges for integration, i.e., Similarly, G 1 (θ) is the LRT of the image Im 1 while the center of the local integral angle range is θ . θ a is an arbitrary angle.
Let the coarsely estimated darkest direction for Im 2 is θ c , the error of the coarse estimation is θ e . From Eq. (5), we can infer that, the LRT of Im 1 that has the best correlation with G 2 (θ c ) is G 1 (θ d1 − θ e ) . This inference can be represented as [corr(G 1 (θ), G 2 (θ c ))]. In Eq. (5), the range for θ is [0, π) . In fact, the optimal θ fluctuates around θ d1 as a result of the small error of coarse estimation. To reduce calculation, the range for θ can be decreased to [θ d1 − θ M , θ d1 + θ M ] . Substituting θ e = θ d2 − θ c into Eq. (5), we can have Equation (6) explores the link between the error of coarse estimation and the correction between LRTs. Based on Eq. (6), the actual darkest direction of Im 2 can be captured by In Eq. (7), the range for θ is [0, π) . In fact, the optimal θ fluctuates around θ d1 as a result of the small error of coarse estimation. To reduce calculation, the range for θ can be decreased to In practice, according to Malus's law, Im 1 can be generated and treated as the model image. Apparently, LRTs of Im 1 also satisfies Malus's law. That is, the integral of the image at the direction θ i is A is a coefficient decided by the image brightness. θ d1 is the darkest direction of the image. Depending on Eq. (8), a set of LRTs of Im 1 (i.e., G 1 (θ i ),θ i = θ d1 − θ M + (i − 1)θ r ) can be gotten. For the input image Im 2 , substituting G 2 (θ c ) into Eq. (7), the corrected darkest direction can be captured. Taking the image in Fig. 1 as an example, the working mechanism is illustrated in Fig. 3. In this experiment, the darkest direction of the model image is 0 • . According to Eq. (8), a set of LRTs of the model image are generated while the center angle changes from −5 • (175 • ) to 5 • , in 0.01 • increment. For the input image shown in Fig. 1, the predicted darkest direction estimated by coarse estimation is θ c = 25.6 • . In EC stage, we found that, G 1 (0.6 • ) has the best correlation value with G 2 (θ c ) . G 1 (□) and G 2 (□) denote the LRTs of the model image and the input image, respectively. Finally, according to Eq. (7), the estimated darkest direction of the input image is corrected to be 25 • .

Implementation details
In practice, once the parameters of the algorithm are given, some intermediate data including the coordinates of pixels used for the coarse estimation, the coordinates and weights of pixels for LRT computation keep unchanged while different input images are treated. Hence, these data can be computed ahead and saved in tables which are named as circle pixel coordinate table (CPCT), integral pixel coordinate table (IPCT), and integral pixel weight table (IPWT), respectively. It should be noted that, due to the different darkest direction of the input images, the coordinates and weights of pixels for gray integration should be saved while the azimuth angle changes from 0 • to 180 • . In addition, the LRTs of the model images with different center angles, which are independent to the input image, also can be captured offline using Eq. (8) and saved.
The flow chart and pseudo code of our method are shown in Fig. 4.

Results and discussion
In this section, the performance of algorithms is tested on synthetic data and real data. We first evaluate the sensitivity of our algorithm (LRT + EC) to the angle range (i.e., θ T ), the angle interval (i.e., θ s ) and the image noise, and then the performance of LRT + EC is compared with two state-of-the-art methods, including GRT [18] and IAD + LC [21]. The CUP of our computer is Intel (R) Core (TM) i7-10710U@1.10G, the RAM is 1.61 GHz.
In practice, the marginal and central parts are not well modulated owing to the vignetting effect of optical system and the imperfection of optical components [18], so the algorithms are performed on the ring area with the inner and outer radiuses are 300 pixels and 500 pixels, respectively. In the following experiments, the radiuses of circles (denoted by r in Eq. (1)) used for coarse estimation changes from 300 to 450 pixels in steps of 5 pixels, and the gray threshold for segmentation is 0.4. The image, in which the darkest direction of hourglass-shaped gray distribution pattern is 0 • , is

Experiments on synthetic data
In this section, experiments are performed on synthetic data. The images are all generated according to Malus's law [26]. We first evaluate the sensitivity of our algorithm (LRT + EC) to some parameters and noise. Moreover, the computational complexity is tested. All experimental results in this section are counted based on 500 Monte Carlo simulations.

Performance to angle range and angle interval of LRT
Generally, the larger the angle range ( θ T ) and the smaller the angle interval ( θ s ), the algorithm will get the higher accuracy. However, the complexity of the algorithm will increase significantly. To get appropriate parameters to guarantee the algorithm acts well in both accuracy and runtime, the performance of the algorithm to θ T and θ s are tested, respectively. In this set of experiments, the Gaussian white noise ( µ = 0.01 and σ 2 =0.005 ) is added on the synthetic images. Figure 5 compared the mean absolute error (MAE) of the results when θ s is 0.6 and θ T changes from 30 • to 85 • . It reviews that the MAE of the algorithm increases rapidly with the increase of θ T , and then becomes flat when θ T is larger than 60 • . Hence, to have low computation complexity, θ T is suggested to be 60 • . Moreover, the MAE of our algorithm while θ s changes from 0.1 • to 2 • is summarized in Fig. 6, and experimental results imply that the smaller θ s is, the better our algorithm performs. As described above, in the following experiments, θ T and θ s are setted to be 60 • and 0.6 • , respectively.

Performance to image noise
Apparently, it is easy to extract the darkest direction precisely from the images without noise, but it is more difficult to have the same precision on the real images which are usually disturbed by the noise. As the noise including shot noise, thermal noise, and dark current noise generally obey the Gaussian distribution [10], to effect the performance of our algorithm to noise, the Gaussian white noise is added while µ = 0.01 and σ 2 changes from 0 to 0.01. MAE of different algorithms are compared in Fig. 7. It can be observed that the performance of all algorithms decrease when the σ 2 increases, and LRT + EC and IAD + LC outperform GRT. For fair comparison, in this set of experiments, the gray threshold for ROI extraction in IAD + LC is also setted to be 0.4. Furthermore, comparing the performance of coarse estimation (the green line shown in Fig. 7) and LRT + EC, we can note that, the MAE of LRT + EC is much smaller. To test the effect of integration on the performance of our algorithm, in EC, we choose the gray of pixels on the circle ( r = 400 ) to replace LRT. The result of this method (the rose line in Fig. 7) explores that the accuracy can be highly improved by gray integration.

Processing time
To make an overall comparison, it is necessary to analyze the processing time of each algorithm. Table 1 indicates that, LRT + EC has the lowest computation complexity.

Experiments on real data
To verify the accuracy of algorithms on real data, 20 modulated intensity images are captured continuously in the same state. One of these images is shown in Fig. 8. The 20 images are analyzed by three different algorithms, and the calculated results are illustrated in Fig. 9. Apparently, LRT + EC is most robust. Furthermore, this conclusion also can be verified by the root mean square error (RMSE) of the results (given in Table 2).

Conclusions
In some spatial polarization modulated polarimetry schemes, the polarization direction of the input light can be achieved by image processing algorithms. In this paper, an efficient image processing algorithm, named LRT + EC, is proposed to extract the polarization direction from the irradiance image of the modulated input light. Different to GRT which performs the Radon transform in the global angle range, LRT is obtained by integrating the image along the radial lines orientated in a local angle range. In addition, LRT and EC can be completed by looking up tables generated offline. Therefore, time consumption of our algorithm is reduced to less than 0.01 s, which meet the realtime requirement well. Moreover, owing to the EC which establish the link between the error of coarse estimation and the correlation between LRTs, the accuracy and robustness of the algorithm are highly improved. The experimental results on synthesized and real data verify that, our proposed algorithm outperforms the state-of-the-art methods including GRT and IAD + LC.