Three-Stage Interpolation Method for Demosaicking Monochrome Polarization DoFP Images

The emergence of polarization image sensors presents both opportunities and challenges for real-time full-polarization reconstruction in scene imaging. This paper presents an innovative three-stage interpolation method specifically tailored for monochrome polarization image demosaicking, emphasizing both precision and processing speed. The method introduces a novel linear interpolation model based on polarization channel difference priors in the initial two stages. To enhance results through bidirectional interpolation, a continuous adaptive edge detection method based on variance differences is employed for weighted averaging. In the third stage, a total intensity map, derived from the previous two stages, is integrated into a residual interpolation process, thereby further elevating estimation precision. The proposed method undergoes validation using publicly available advanced datasets, showcasing superior performance in both global parameter evaluations and local visual details when compared with existing state-of-the-art techniques.


Introduction
As optical imaging and its applications have been rapidly evolving, the study of light's polarization properties emerges as an indispensable facet, providing higher-dimensional insights into the orientation of light wave oscillations.Th polarization of light contains rich information to facilitate a more profound comprehension and interpretation of the physical world.For example, some of the species in nature exhibit acute responsiveness to polarization cues and utilize polarized light for positioning and detection [1].
Notwithstanding the considerable potential of polarized imaging, it remains beyond the perceptual scope of the human eye, necessitating the deployment of specialized polarized imaging apparatus.Recent years have witnessed the rapid proliferation of polarized imaging and measurement technologies across diverse academic disciplines and practical applications, including image recovery [2,3], surface inspection [4], and 3D shape reconstruction [5].
The acquisition of polarization information from target images typically entails the capture of a set of in-site images employing polarizers placed at different angular orientations in front of the camera [6].This enables the computational assessment and analysis of polarization attributes at the level of individual pixels within the acquired images.However, such a multi-shot approach is unsuitable for dynamic scenes and video acquisition.Recently, a novel category of polarization image sensors using the division of focal plane (DoFP) technique emerged as a pioneering sensing solution [7].Analogous to conventional color cameras employing color filter arrays (CFAs) [8], these sensors integrate micro-polarizer arrays into the pixels of focal plane array sensors.Such polarization filter arrays (PFAs) allow for the formation of polarizers at different angles on the chip during the Sensors 2024, 24, 3018 2 of 13 semiconductor-manufacturing process, enabling highly precise alignment with the pixels.The typical monochrome polarizing filter array (MPFA) used in DoFP polarization cameras consists of a 2 × 2 period of four polarizers with angles set at 0, 45, 90, and 135 degrees, as shown in Figure 1.Additionally, the color polarization filter array (CPFA) combines the principles of color imaging with polarization imaging, and features a fundamental 4 × 4 periodic structure with unique color and polarization orientation at each position.
Sensors 2024, 24, x FOR PEER REVIEW 2 of 13 arrays (PFAs) allow for the formation of polarizers at different angles on the chip during the semiconductor-manufacturing process, enabling highly precise alignment with the pixels.The typical monochrome polarizing filter array (MPFA) used in DoFP polarization cameras consists of a 2 × 2 period of four polarizers with angles set at 0, 45, 90, and 135 degrees, as shown in Figure 1.Additionally, the color polarization filter array (CPFA) combines the principles of color imaging with polarization imaging, and features a fundamental 4 × 4 periodic structure with unique color and polarization orientation at each position.The emergence of DoFP polarization cameras provides better opportunities for capturing polarization information, but it also presents new challenges that need to be addressed for effective full-polarization reconstruction.Image demosaicking, also referred to as mosaic interpolation, is a crucial process in digital imaging with CFAs and also PFAs.In the case of MPFAs, where each pixel location holds only one polarization orientation, retrieving data of the other polarization channels requires interpolation or reconstruction for accurate estimation.MPFA demosaicking based on interpolation-based methods is also an integral component of CPFA interpolation-based demosaicking.This process is oriented towards the completion of absent polarization information at each pixel location by extrapolating values of the other three polarization orientations from the pixel's neighborhood.
Compared to CFA demosaicking [10,11], that with MPFAs inherently presents greater complexity due to the increased number of channel information within an array and the challenges in addressing intricate polarization characteristics.The effective polarization estimation and management of inter-channel correlations become imperative in this context.In tackling the demosaicking challenge in monochrome polarization images, some conventional interpolation methods employed in CFA demosaicking, including bilinear interpolation, bicubic interpolation, and cubic spline interpolation, have been adapted for MPFA demosaicking [12].However, these methods primarily rely on spatial correlations and overlook inter-channel correlations.Given the disparities in array types, many well-structured CFA demosaicking algorithms may also not achieve optimal results, owing to inherent distinctions in color and polarization mechanisms [13].Zhang et al. [14] proposed an interpolation method that enhances the demosaicking quality of DoFP polarization images using intensity correlations.They performed dual bicubic interpolation with edge detection for both oblique and vertical pixels, leading to improved results.Recognizing that high-frequency information has lower polarization direction differences, Li et al. [15] proposed a polarization differential domain interpolation method using Newton polynomial interpolation.Morimatsu et al. [16] presented a method using intensityguided edge-aware residual interpolation based on channel differences, demonstrating notable improvements despite the requirement for complex computations.In a recent development, Wu et al. [17] introduced a simple yet effective polarization image demosaicking technique, leveraging prior information about polarization channel differences to The emergence of DoFP polarization cameras provides better opportunities for capturing polarization information, but it also presents new challenges that need to be addressed for effective full-polarization reconstruction.Image demosaicking, also referred to as mosaic interpolation, is a crucial process in digital imaging with CFAs and also PFAs.In the case of MPFAs, where each pixel location holds only one polarization orientation, retrieving data of the other polarization channels requires interpolation or reconstruction for accurate estimation.MPFA demosaicking based on interpolation-based methods is also an integral component of CPFA interpolation-based demosaicking.This process is oriented towards the completion of absent polarization information at each pixel location by extrapolating values of the other three polarization orientations from the pixel's neighborhood.
Compared to CFA demosaicking [10,11], that with MPFAs inherently presents greater complexity due to the increased number of channel information within an array and the challenges in addressing intricate polarization characteristics.The effective polarization estimation and management of inter-channel correlations become imperative in this context.In tackling the demosaicking challenge in monochrome polarization images, some conventional interpolation methods employed in CFA demosaicking, including bilinear interpolation, bicubic interpolation, and cubic spline interpolation, have been adapted for MPFA demosaicking [12].However, these methods primarily rely on spatial correlations and overlook inter-channel correlations.Given the disparities in array types, many wellstructured CFA demosaicking algorithms may also not achieve optimal results, owing to inherent distinctions in color and polarization mechanisms [13].Zhang et al. [14] proposed an interpolation method that enhances the demosaicking quality of DoFP polarization images using intensity correlations.They performed dual bicubic interpolation with edge detection for both oblique and vertical pixels, leading to improved results.Recognizing that high-frequency information has lower polarization direction differences, Li et al. [15] proposed a polarization differential domain interpolation method using Newton polynomial interpolation.Morimatsu et al. [16] presented a method using intensity-guided edge-aware residual interpolation based on channel differences, demonstrating notable improvements despite the requirement for complex computations.In a recent development, Wu et al. [17] introduced a simple yet effective polarization image demosaicking technique, leveraging prior information about polarization channel differences to process MPFA demosaicking.While Xin et al. [18] made further enhancements, neither one nor the other of these methods yielded significantly superior results.
Numerous learning-based methods have also been proposed and applied to the PFA demosaicking in recent years, and some have achieved commendable results.For MPFA Sensors 2024, 24, 3018 3 of 13 demosaicking specifically, there exist dictionary-learning-based [19][20][21][22] and deep-learningbased methods [23][24][25][26], but their efficacy is highly contingent on the datasets employed during training and the demand a high computational costs and runtimes.These limitations underscore the impracticality of their actual applications in image sensors and emphasize the crucial need to further investigate intensity-based interpolation methods.
In this paper, we present a fast and highly effective demosaicking model comprising a three-stage interpolation procedure designed to address this challenging issue.Our model achieves superior results compared with state-of-the-art methods.The first two stages introduce a novel linear interpolation method tailored for MPFAs, incorporating customized weights for neighboring pixel values.For the first stage, we separately interpolate along the 45-degree and 135-degree diagonal directions to recover the missing polarization information that is perpendicular to the polarizer angle formed on the chip.A continuous adaptive edge detection method that relies on the variance difference between these two diagonal directions is then utilized to perform weighted averaging for this bidirectional interpolation.For the second stage, we apply the same method twice in both horizontal and vertical directions, thereby acquiring pixel values in the remaining two absent directions for each pixel location.The complete polarization information for each pixel position is obtained hereto, but we choose to refine the results through a residual interpolation step in the third stage.Instead of using the residual value between two different polarization orientations, we utilize the total intensity value calculated from the results obtained in the previous stages.
Our proposed method has undergone validation on publicly available datasets, demonstrating its superior performance in global evaluation and local details when compared with existing intensity-based interpolation methods.Additionally, it offers relatively fast computational times and holds the potential for further enhancement through parallel acceleration.

Polarization Channel Difference Interpolation
Most CFA interpolation methods are based on color channel difference interpolation [27,28].The same holds true for PFA interpolation, as many studies have discussed the difference interpolation between polarization channels [15,17,18].After various tests and verification, a relatively reliable assumption is made: high-frequency energy is reduced in the polarization channel difference domain [17].By using the channel difference for interpolation, demosaicking accuracy can be improved compared with that of independent interpolation for each polarization channel.
The purpose of basic linear polarization channel difference interpolation is to obtain the average of the polarization channel differences on both sides of the pixel.Figure 2 is a 7 × 7 pixel array arranged in an MPFA pattern.Each angle marked in a small square represents the polarization orientation at a specific pixel position within the MPFA.In Figure 2, I 0 (i, j) is the available information in the central position, while I d (i, j), d ∈ {45, 90, 135}, denotes the missing values that require interpolation.For illustration, only horizontal interpolation is performed in this section, and the polarization information of orientation 135 • can be estimated in position (i, j), as shown below: To avoid confusion with polarization orientations, "⋯" is used in this paper to represent the horizontal interpolating direction.We use "~" to indicate the unknown value that needs to be temporarily represented and does not participate in the actual calculation process."^" is used to represent the estimated result obtained in each interpolation stage.In Equation ( 1),  ⋯ (, ) represents the value estimated through horizontal interpolation.All remaining terms except for  (,  − 1) and  (,  + 1) are known in the original MPFA image.The challenge then lies in estimating these two unknown values.Most methods simply treat these two terms as averages of their respective left and right neighbors [13].
After rearranging Equations ( 1) and ( 2), we obtain the following: where the first term is the first-order mean, and the second term, , which could be written as   ⋯ (, ), represents the second-order gradient.
While this method is effective, it may not provide sufficient precision for MPFA interpolation at a higher level.If we repeatedly estimate  (,  − 1) and  (,  + 1) using the same method as that applied to  ⋯ (, ) in Equation ( 1), it will contain a relatively high proportion of  (,  − 1) and  (,  + 1).This would be consequently subtracted by the first-order mean and potentially impact the primary interpolation results.
Therefore, we propose a compromise solution in which we replace the interpolation of the polarization channel difference,  (,  1) −  (,  1) , with a new approach:  (,  1) −  (,  1).Then, Equation ( 2) is replaced by the following: For example, (i, j) is the location of the center pixel in the pattern.
To avoid confusion with polarization orientations, "• • • " is used in this paper to represent the horizontal interpolating direction.We use "~" to indicate the unknown value that needs to be temporarily represented and does not participate in the actual calculation process."ˆ" is used to represent the estimated result obtained in each interpolation stage.In Equation (1), Î••• 135 (i, j) represents the value estimated through horizontal interpolation.All remaining terms except for ∼ I 0 (i, j − 1) and ∼ I 0 (i, j + 1) are known in the original MPFA image.The challenge then lies in estimating these two unknown values.Most methods simply treat these two terms as averages of their respective left and right neighbors [13].
After rearranging Equations ( 1) and ( 2), we obtain the following: where the first term is the first-order mean, and the second term, , which could be written as , represents the second-order gradient.While this method is effective, it may not provide sufficient precision for MPFA interpolation at a higher level.If we repeatedly estimate ∼ I 0 (i, j − 1) and ∼ I 0 (i, j + 1) using the same method as that applied to Î••• 135 (i, j) in Equation ( 1), it will contain a relatively high proportion of I 135 (i, j − 1) and I 135 (i, j + 1).This would be consequently subtracted by the first-order mean and potentially impact the primary interpolation results.
Therefore, we propose a compromise solution in which we replace the interpolation of the polarization channel difference, ∼ I 0 (i, j ± 1) − I 135 (i, j ± 1), with a new approach: . Then, Equation ( 2) is replaced by the following: Sensors 2024, 24, 3018 5 of 13 where ∼ I 135 (i, j) and ∼ I 135 (i, j ± 2) can be obtained by averaging their neighbors, as in Equation ( 2).After simplifying Equations ( 1) and ( 4), we arrive at the final interpolation procedure in this stage:

Interpolation Sequence and Edge Detection
Section 2.1 exclusively focuses on clarifying linear interpolation in a single direction.However, in many cases, the values within an image can exhibit significant variations, especially near edges.When solely relying on one-directional interpolation for missing values, the resulting outcome may lack reliability.Consequently, there is a need for multi-directional interpolation.In actual practice, bidirectional interpolation is commonly employed, while edge detection is then highly required for directional weighting.
In the process of CFA demosaicking, it is a common practice to prioritize the interpolation of the green channel, both horizontally and vertically, due to it having double the number of pixels compared with other color channels.However, as for MPFA interpolation, polarization information is uniformly distributed in channels across the four directions on the chip, as shown in Figure 1.Bidirectional interpolation for missing polarization information poses a more challenging task.
In our initial stage, we employ diagonal interpolations along the 45-degree and 135-degree directions as a deliberate choice, instead of common interpolation starting from horizontal and vertical directions.This approach allows us to simultaneously estimate information from two perpendicular directions and thus minimizes interference stemming from incorrect directional interpolation near image edges.We can then gather polarization information perpendicular to the polarization orientation of each pixel's polarizer on the MPFA, as illustrated in Figure 3a.
Similarly, we use ". . ." and " . .." to denote the 45-degree and 135-degree interpolation directions, respectively.Î90 (i, j) in bidirectional interpolation can then be expressed by the following equation: where Î . . .90 (i, j) and Î . . .90 (i, j) are the results of interpolation in the respective directions.While ω . . .90 (i, j) and ω . . .90 (i, j) represent the coefficients of directional weights assigned to their corresponding directions, ω . . .90 (i, j) + ω . . .90 (i, j) = 1 can be obtained.Both weights are determined based on the rate of change of pixel values in the respective directions within a specific distance.In many studies [13,14,27], the rate of change, taking the position (i, j) along the 45-degree interpolation direction as an example, is typically defined as follows: When changes along a certain direction are more rapid than those along the perpendicular direction, they imply the potential presence of an edge at the corresponding pixel location.In such cases, it is advisable to assign a smaller coefficient to that direction, given the uneven variation across the edge along that direction.Many edge detection methods adopt a similar coefficient assignment strategy: if one direction undergoes more significant changes than the other, the coefficient for that direction is set to 0, and in the converse scenario, it is set to 1.When the changes in both directions are equal, the coefficients for both directions are set to 0.5 [8,27].However, this common approach may result in discontinuities across the image and improper estimated values due to the rare occur-rence of ν . . .90 (i, j) = ν . . .90 (i, j).This leads to a situation where unidirectional interpolation is predominantly employed in practical applications instead of the expected bidirectional interpolation.
number of pixels compared with other color channels.However, as for MPFA interpolation, polarization information is uniformly distributed in channels across the four directions on the chip, as shown in Figure 1.Bidirectional interpolation for missing polarization information poses a more challenging task.
In our initial stage, we employ diagonal interpolations along the 45-degree and 135degree directions as a deliberate choice, instead of common interpolation starting from horizontal and vertical directions.This approach allows us to simultaneously estimate information from two perpendicular directions and thus minimizes interference stemming from incorrect directional interpolation near image edges.We can then gather polarization information perpendicular to the polarization orientation of each pixel's polarizer on the MPFA, as illustrated in Figure 3a. (b,c) The second stage of obtaining Î135 (i, j) and Î45 (i, j), respectively. (d-f) The third stage of obtaining the final estimation of Î90 (i, j), Î45 (i, j), and Î135 (i, j).
An enhanced edge detection method incorporates the utilization of gradients in both directions and computes coefficients in accordance with the following formulation: where ε is a small coefficient introduced to avoid division by zero.This method is more effective and robust.However, in cases where ν . . .90 (i, j) and ν . . .90 (i, j) manifest substantial disparities, the smaller weight still contributes as a non-negligible part.This often necessitates the additional establishment of complex boundary conditions for constraints [15].These conditions are typically based on their ratio, but may lack universality for various image patterns.Moreover, when ν . . .90 (i, j) and ν . . .90 (i, j) are small values, the ratio between them can easily become imbalanced, leading to an approximation of unidirectional interpolation.In light of this, directional coefficients are suggested to be established based on the difference between the rates of change in this paper.To ensure that the weights change Sensors 2024, 24, 3018 7 of 13 continuously with the difference, ω . . .90 (i, j) and ω . . .90 (i, j) should have the same form and satisfy the following requirements: where f represents the custom function determined to satisfy the equation and maintain the sum of the two directional weights.Similar to the approach in [27], the following function is selected to determine the coefficients of directional weights: and ω . . .
where k is a positive number used to adjust the changes of weights.In our experiments, we set parameter k to 20 when pixel values in the MPFA image ranged from 0 to 1.
Once we acquired the polarization information of each pixel in the orthogonal direction through the first stage of interpolation, we proceed with the same method in the vertical and horizontal directions to obtain the remaining missing information for each pixel.For position (i, j), the second stage aims to obtain Î135 (i, j) and Î45 (i, j), as shown in Figure 3b,c.
To illustrate, the estimation of Î45 (i, j) through bidirectional interpolation is taken as an example, as shown in Equation (11).The vertical result, Î. . .45 (i, j), is obtained through interpolation based on values directly from the original MPFA image and goes in the same direction as the green arrow does.On the other hand, Î••• 45 (i, j) is derived from interpolation along the direction of the yellow arrow, utilizing the results obtained from the initial stage of our procedure to proceed with the second stage:

Third-Stage Interpolation
After completing two stages of bidirectional interpolations, as described in 2.1 and 2.2, the full polarization information was retrieved and the results were already adequate for calculating the polarization parameters, including total light intensity ( Ŝ0 ), degree of linear polarization (DOLP), and angle of linear polarization (AOLP).
However, we decided to go further to enhance estimation accuracy by employing a method similar to that of residual interpolation in CFA [11].For each position, we computed the mean of the four polarization orientations, which accounts for half of the total intensity.For position (i, j) in Figure 2, the mean value is calculated via the following: Next, we calculated the discrepancy between the mean value and the corresponding pixel value in the original MPFA image, resulting in the residual value for each pixel position.The residual map is decomposed into four images, I 0 * , I 90 * , I 45 * , and I 135 * , as illustrated in Figure 4.Each contains residuals of a different orientation, with the remaining positions filled with zeros.Subsequently, we performed linear interpolation of these four images via convolution.After incorporating the mean value, we restored the original MPFA values and updated pixel values at zero positions for these four images, as shown in the following formula: where I maskD MFPA is one of the decomposed images, with D ∈ {0, 45, 90, 135}, and F is the filter of the convolution: Sensors 2024, 24, x FOR PEER REVIEW 8 of 13

Experiment
To evaluate the effectiveness of our proposed method, we employed a publicly available dataset [16], comprising 40 distinct scenes.Each scene consisted of 12 monochrome images, each with a size of 1024 × 768 pixels, representing different color channels and polarization orientations.Notably, this dataset was captured using a 3-CCD camera, enabling the separate recording of each color channel.Unlike images captured with common cameras, sets of four images with the same color channel and scene but taken at different polarizer angles (0°, 45°, 90°, and 135°), excluding the need for CFA demosaicking, prove to be a more suitable choice for evaluating MPFA demosaicking solutions compared with other public datasets.
In our study, we utilized and tested all three color channels separately for experiments.To simulate the polarization images captured by an actual monochrome DoFP camera, each set was down-sampled, creating an MPFA pattern identical to the one depicted in Figure 2. Consequently, each scene, stored in four captured color images, could generate three simulated MPFA images for testing, as illustrated in Figure 5.Besides our proposed method, we also tested other state-of-the-art demosaicking algorithms to demonstrate our method's superiority, including intensity correlation in polarization channel (ICPC, 2016) [14], pseudo-panchromatic image difference (PPID, 2018) [13], Newton's polynomial interpolation and polarization difference (NPPD, 2019) [15], Taking pixel position (i, j) in Figure 2 as an example, the original value, I 0 (i, j), in MPFA remains unchanged in I 0 * , while the additional polarization values Î90 (i, j), Î45 (i, j), and Î135 (i, j) are estimated, respectively, in I 90 * , I 45 * , and I 135 * , as illustrated in Figure 3d-f.
The entire procedure of our proposed method, including the first two stages, is illustrated in Figure 4.This diagram provides a comprehensive overview of the process, and the final estimation of the complete polarization information can be obtained from MPFA.

Experiment
To evaluate the effectiveness of our proposed method, we employed a publicly available dataset [16], comprising 40 distinct scenes.Each scene consisted of 12 monochrome images, each with a size of 1024 × 768 pixels, representing different color channels and polarization orientations.Notably, this dataset was captured using a 3-CCD camera, enabling the separate recording of each color channel.Unlike images captured with common cameras, sets of four images with the same color channel and scene but taken at different polarizer angles (0 • , 45 • , 90 • , and 135 • ), excluding the need for CFA demosaicking, prove to be a more suitable choice for evaluating MPFA demosaicking solutions compared with other public datasets.
In our study, we utilized and tested all three color channels separately for experiments.To simulate the polarization images captured by an actual monochrome DoFP camera, each set was down-sampled, creating an MPFA pattern identical to the one depicted in Figure 2. Consequently, each scene, stored in four captured color images, could generate three simulated MPFA images for testing, as illustrated in Figure 5.
Besides our proposed method, we also tested other state-of-the-art demosaicking algorithms to demonstrate our method's superiority, including intensity correlation in polarization channel (ICPC, 2016) [14], pseudo-panchromatic image difference (PPID, 2018) [13], Newton's polynomial interpolation and polarization difference (NPPD, 2019) [15], edgeaware residual interpolation (EARI, 2019) [16], polarization channel difference prior (PCDP, 2021) [17], and polarization difference based on edge compensation (PDEC, 2023) [18], on the extracted down-sampled images with an MPFA pattern.The codes of all these methods were downloaded from the internet and implemented in MATLAB.We slightly adjusted and standardized these codes for comparative purposes.This ensured the consistency of the input and output formats of various methods for the effective evaluation of performance during the comparison.
In our study, we utilized and tested all three color channels separately for experiments.To simulate the polarization images captured by an actual monochrome DoFP camera, each set was down-sampled, creating an MPFA pattern identical to the one depicted in Figure 2. Consequently, each scene, stored in four captured color images, could generate three simulated MPFA images for testing, as illustrated in Figure 5.Besides our proposed method, we also tested other state-of-the-art demosaicking algorithms to demonstrate our method's superiority, including intensity correlation in polarization channel (ICPC, 2016) [14], pseudo-panchromatic image difference (PPID, 2018) [13], Newton's polynomial interpolation and polarization difference (NPPD, 2019) [15], edge-aware residual interpolation (EARI, 2019) [16], polarization channel difference prior (PCDP, 2021) [17], and polarization difference based on edge compensation (PDEC, 2023) [18], on the extracted down-sampled images with an MPFA pattern.The codes of all these methods were downloaded from the internet and implemented in MATLAB.We slightly adjusted and standardized these codes for comparative purposes.This ensured the consistency of the input and output formats of various methods for the effective evaluation of performance during the comparison.The full polarization information (I0, I45, I90, and I135) of each pixel position could be obtained with these methods.Polarization parameters S0, DoLP and AoLP are also computed separately using different methods for comparison.For quality assessment, we compared the recovered results with ground truth values to evaluate the performance of each method.The root mean square error (RMSE) is selected in this paper to quantify how far predictions are from the ground truth values, and a lower RMSE indicates better performance.

RMSE
where m and n represent the pixel dimensions of the images.In the dataset used in our study, m = 1024, and n = 768.I est p denotes the estimated value of the compared polarization parameter, and I gt p is the ground truth value obtained from the captured images.In some previous studies, PSNR was regarded as a critical indicator of image quality [11,20].However, it is essential to note that the range of PSNR varies at different scales and is represented on a logarithmic scale.It can be readily derived from RMSE, as shown in Equation ( 16): where MAX I is the maximum possible pixel value of the images.During these measurements, we excluded the 12-pixel border around the images to mitigate any boundary effects related to demosaicking.Additionally, we measured the runtime for each demosaicking method.
In order to enhance the visualization and make better comparisons, the error images based on an intensity value of S0 obtained through different methods mentioned above are presented, which quantify the deviation of the estimated value from the ground truth.Considering the difficulty in comparing the DOLP results in Figure 5, we simultaneously represent the estimated values of both DOLP and AOLP within an image, following the structure of HSV color space, which is analogous to that in [16].Specifically, we mapped the range of AOLP angles to the hue channel within the [0, π] range.Given that DOLP values tend to be relatively small, we scaled them up by a factor of 20 and mapped the resulting values to the saturation channel.The V channel was kept constant at 0.5 to maintain the image intensity at a moderate level.These visual comparisons in specific regions allowed us to depict the polarization information in a manner that made it easier to discern and compare the outcomes achieved through various methods.

Overall Performance Comparison
Table 1 presents the average RMSE values of the 40 scenes in a single-color channel within the dataset for different demosaicking methods.One color channel in a particular scenario may exhibit a considerably higher intensity compared with the other two channels, Sensors 2024, 24, 3018 10 of 13 leading to notable disparities in the error of estimated results during MPFA demosaicking.Consequently, comparison results obtained from a single-color channel may vary.The average RMSE of the three color channels are thus all presented, and the results of these methods reveal a nearly consistent difference in deviation among the three color channels, effectively demonstrating and distinguishing the varying performance of these methods.Among these results, our method consistently demonstrates superior performance across most parameters, encompassing the four polarization images (I0, I45, I90, and I135),Stokes image S0, and the AOLP image.The EARI method also exhibits noteworthy efficacy among the existing methods, particularly in DOLP results, as it attains the lowest RMSE values in DOLP images within the green and blue color channels.However, it does not match the overall performance of our method across most parameters.
Moreover, it is important to note that the EARI method necessitates significantly more computational time compared with other methods, whereas our method operates nearly three times faster than EARI does under the same conditions, as shown in Table 2.

Visual Performance in Specific Regions
In the context of images captured using MPFA cameras, it is crucial to note that regions with reflections and glass surfaces may exhibit more pronounced variations in polarization information compared with other areas.Nevertheless, these specific scenarios may not be adequately represented within the dataset.Consequently, it becomes imperative to account for these cases when assessing demosaicking algorithms.We then extracted specific 256 × 256-pixel regions from the captured images in the dataset, where particular objects such as bottles and bulbs occupy a significant portion of the region, as illustrated by the red box region from the original image in Figure 6.These were chosen to assess the algorithm's performance in these particular scenarios.
degree of deviation of each pixel location within the selected region.The vivid visual representations in these specific regions serve as a testament to the effectiveness of our approach in preserving the finer details of polarization information.The pseudo-colored images, which span from blue to red, offer a clear and intuitive means of conveying the extent of deviation in each pixel within the selected regions.Notably, the ground truth is consistently depicted in blue, allowing for the easy differentiation and comparison of outcomes achieved via various methods.In the comparative analysis of the localized regions selected from the two scenarios, it is evident that our method outperforms others.In the S0 domain, the degree of deviation is relatively low, indicating the accuracy and precision of our approach in capturing polarization information.Within the AOLP-DOLP representation, we observe a reduced level of blurriness, further underlining the superior performance of our technique.In conclusion, the localized analysis of these specific regions not only showcases the superior performance of our method in capturing polarization information with precision but also For detailed comparison, we further selected specific local regions from the AOLP-DOLP results.The error images based on intensity are also presented, which quantify the deviation from the ground truth of S0.We use pseudo-colored images to represent the degree of deviation of each pixel location within the selected region.The vivid visual representations in these specific regions serve as a testament to the effectiveness of our approach in preserving the finer details of polarization information.The pseudo-colored images, which span from blue to red, offer a clear and intuitive means of conveying the extent of deviation in each pixel within the selected regions.Notably, the ground truth is consistently depicted in blue, allowing for the easy differentiation and comparison of outcomes achieved via various methods.
In the comparative analysis of the localized regions selected from the two scenarios, it is evident that our method outperforms others.In the S0 domain, the degree of deviation is relatively low, indicating the accuracy and precision of our approach in capturing polarization information.Within the AOLP-DOLP representation, we observe a reduced level of blurriness, further underlining the superior performance of our technique.In conclusion, the localized analysis of these specific regions not only showcases the superior performance of our method in capturing polarization information with precision but also highlights the reduced blurriness within the AOP-DOP domain.These results affirm the effectiveness and potential of our approach in applications involving MPFA cameras and situations with reflections and glass surfaces, which often present challenges in polarization imaging.

Conclusions
In summary, we presented a novel MPFA demosaicking model that demonstrated exceptional performance, surpassing current state-of-the-art techniques in various crucial aspects.Through a well-structured three-stage interpolation procedure, we achieved remarkable success in the more accurate recovery of polarization information, setting a new benchmark in the field.The validation on a public diverse dataset highlights our method's superiority in terms of RMSE values, especially in capturing different polarization parameters including the four polarization images (I0, I45, I90, and I135), Stokes image S0, and the AOLP image.Specific region analysis further underscored the precision and accuracy of our technique, particularly in challenging scenarios involving reflections and glass surfaces.These results collectively underscore the potential and effectiveness of our approach in the realm of MPFA cameras and challenging polarization imaging applications.
Looking ahead, we are committed to exploring algorithm optimization to enhance computational efficiency.Though our approach excels in computational efficiency, operating nearly three times faster than the closest competitor, it still has high potential for further acceleration.Moreover, our work will encompass a broader range of environmental conditions and scene complexities, ensuring the robustness and versatility of our approach.This ongoing research will undoubtedly contribute to the continual evolution of polarization imaging, with far-reaching implications in domains including remote sensing, medical diagnostics, and advanced material characterization.

Figure 1 .
Figure 1.Monochrome polarizing filter array of an industrial image sensor.(a) Sony polarization image sensor IMX250MYR-C.(b) Macro image of the four directional polarizers [9].(c) MPFA layout of four polarization angles.

Figure 1 .
Figure 1.Monochrome polarizing filter array of an industrial image sensor.(a) Sony polarization image sensor IMX250MYR-C.(b) Macro image of the four directional polarizers [9].(c) MPFA layout of four polarization angles.

Figure 2 .
Figure 2. A classic 7 × 7 MPFA pattern with rectangular coordinates for algorithm representation.For example, (, ) is the location of the center pixel in the pattern.

Figure 2 .
Figure 2. A classic 7 × 7 MPFA pattern with rectangular coordinates for algorithm representation.For example, (i, j) is the location of the center pixel in the pattern.

Figure 3 .Figure 3 .
Figure 3. Illustration of interpolation procedure at coordinate (, ) using our three-stage MPFA demosaicking method.The degrees within the grid represent the intensity values of the corresponding Figure 3. Illustration of interpolation procedure at coordinate (i, j) using our three-stage MPFA demosaicking method.The degrees within the grid represent the intensity values of the corresponding polarization orientation at each location.(a) The first stage of obtaining Î90(i, j). (b,c) The second stage of obtaining Î135 (i, j) and Î45 (i, j), respectively. (d-f) The third stage of obtaining the final estimation of Î90 (i, j), Î45 (i, j), and Î135 (i, j).

Figure 4 .
Figure 4. Overall process of our proposed method.

Figure 5 .
Figure 5. Input and output of the demosaicking algorithm.

Figure 4 .
Figure 4. Overall process of our proposed method.

Figure 5 .
Figure 5. Input and output of the demosaicking algorithm.

Figure 5 .
Figure 5. Input and output of the demosaicking algorithm.

Table 1 .
The average RMSE comparisons between different methods for each color channel.

Table 2 .
Running time comparisons between different methods for a single MPFA demosaicking process.