Parameter-based imaging from passive multispectral polarimetric measurements

: A modiﬁed pBRDF model with a diﬀuse scattering component is applied to estimate the complex refractive index, slope variance roughness, and diﬀuse scattering coeﬃcients of object surfaces from time sequences of polarimetric images. The approach is used for the ﬁrst time to produce parameter-based images from multispectral Stokes imagery of outdoor target scenes collected by the Ground Multiangle Spectro-Polarimetric Imager. The images of the estimated surface parameters show distinctions between diﬀerent objects in the scenes and the parameter values are consistent with reasonable expectations for the object surfaces


Introduction
Optical polarimetry allows the measurement or estimation of surface material properties and supports remote sensing applications [1] such as material classification [2], shape extraction [3], and target recognition [4].Surface parameters of interest include the refractive index, scattering function, and surface roughness.Early work on the utility of optical polarimetry for surface parameter estimation includes Priest et al. [5] who initially derived a polarimetric bidirectional reflectance distribution function (pBRDF) based on a microfacet model by Torrance and Sparrow [6] to describe both the strength and polarization relationships between the incident and scattered Stokes parameters.Thilak et al. [7] and Martin et al. [8] utilized this microfacet pBRDF model to estimate the index of refraction from multiple degree-of-polarization (DOP) measurements of the reflected intensity as the source-target-receiver geometry changes.The microfacet model assumes that a rough surface is composed of a collection of small, randomly oriented, perfect specular facets obeying Fresnel reflection and incident light is assumed to be reflected only once before returning to the observer.Thus, the microfacet pBRDF, in general, cannot describe depolarization phenomena and its application for refractive index estimation is limited to specular materials.
Gartley et al. [9] and Litvinov et al. [10] modified the microfacet pBRDF by introducing several coefficients to account for the depolarization behavior.The coefficient values are determined by fitting the model to measured data but the development of the coefficients was not specifically physics-based.Hyde et al. [11] developed a parameter ρ DHR to describe diffuse scattering that is based on directional hemispherical reflectance (DHR) and energy conservation theories.Following the approach of Hyde et al., Yang et al. [12] jointly estimated the refractive index and surface roughness of several non-specular samples from a series of DOP observations.However, their model and estimation approach are only valid for highly reflective materials since optical absorption inside the material is neglected.In addition, ρ DHR is a complex integral function of the scattering angles, the surface roughness parameters and the complex refractive index over the entire hemisphere above the material surface.The formulation of ρ DHR suggests it is difficult to measure directly and estimating surface parameters using the associated model could be problematic for remote sensing applications [12].We also note that the previous work cited in [7][8][9][10][11][12] addresses a fixed point on the surface of interest and the evaluations are limited to indoor laboratory measurements.
In recent studies [13][14], we describe two independent approaches for estimating the refractive index, roughness and diffuse scattering coefficient of a surface from measurements of the degree of polarization (DOP) and the Stokes parameters, respectively.Theoretical developments were presented that incorporate a diffuse scattering component in the description of the observed DOP and Stokes signals.The parameter estimation results using both approaches are encouraging, and synthetic simulation results of DOP and Stokes parameters generated by the proposed models show close correspondence to measured values.These proposed approaches significantly simplify the surface parameter estimation process without losing accuracy especially for remote sensing where target surfaces are typically rough compared to the incident wavelength.Furthermore, it was recognized that using Stokes parameter measurements rather than the DOP could provide a richer set of data and improve estimation for a limited number of observations or in the presence of measurement noise.
In this paper, we apply the Stokes parameter-based estimation approach to polarimetric image data of outdoor scenes.To our knowledge, this is the first time that images of surface parameters have been generated from polarimetric field data.The polarimetric data for this work was collected with the Ground Multiangle SpectroPolarimetric Imager (Ground-MSPI) [15][16], which is a camera system developed to provide Stokes parameter images at narrowband wavelengths of 470 nm, 660 nm and 865 nm [17][18][19][20][21][22].

Model Development
In this section, we summarize the modified polarimetric bidirectional reflectance distribution function (pBRDF) model and the parameter estimation approach used for this work [14].Figure 1 defines the scattering geometry for the pBRDF and the measurements.Incident sunlight is assumed to be unpolarized with a Stokes vector S I =[S 0 I 0 0] T and it illuminates a target with a complex index of refraction η=n-ik where n and k are the real and imaginary part of the complex refractive index and i = (−1) 1/2 .The scattered light, defined by the Stokes vector S = [S 0 S 1 S 2 ] T , is observed by a detector such as the Ground-MSPI camera.The mean surface normal, indicated with z, is relative to the mean plane associated with the surface's macro-surface.θ i and θ r are the incident and reflected zenith angles relative to the mean surface normal; φ i and φ r are the incident and reflected azimuth angles, and a relative azimuth angle is defined as φ=φ r -φ i .Although the direct solar illumination of the surface is unpolarized, the surface also receives downwelling irradiance from the sky that has a polarization preference.For clear sky conditions, the direct illumination generally dominates unless the Sun is low in the sky.Thus, in this work, as well as other passive remote sensing studies [1,[15][16][17][18][19][20][21][22], the scene illumination is considered to be unpolarized.This assumption may not be accurate for certain conditions such as a partly cloudy sky or for relatively smooth surfaces that produce significant specular reflections of the partially polarized sky irradiance.However, if the illumination polarization is known, the extension of the analysis to this case is relatively straightforward.We note that the fourth Stokes component, corresponding to the circular polarization state, is not measured by the Ground-MSPI due to its insignificant contribution in most natural scenes [1][2]4,[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].
For this work we separate scattering phenomena into two general categories: surface and diffuse [9][10][11][12][13][14]. Surface scattering occurs at the interface between air and the target surface.Diffuse scattering is usually considered to be due to subsurface scattering in the material volume that re-emerges into air.For a rough surface it is also possible to have multiple reflections involving a number of surface facets that can introduce an additional depolarization component.However, we only consider surface scattering from a single surface reflection and any multi-surface scattering is presumed incorporated in the diffuse scattering component.We assume the diffuse scattering is completely unpolarized and Lambertian where the light encounters enough multiple scattering events within the bulk of the material that the scattering directions and polarization state are uniformly random.In general, the diffuse scattering may exhibit some polarization and directional preference due to the re-emergence processes [23], but the assumptions of unpolarized and Lambertian scattering are reasonable especially for remote sensing scenarios where the target's surface roughness is typically large compared to the incident wavelength [9][10][11][12][13][14][15][16][17][18][19][20][21][22].In contrast, the surface scattering that develops based on a microfacet model, can produce a directional lobe and polarization preference.
The modified pBRDF, which relates the incident (S I ) and scattered (S) Stokes parameters and includes both surface (f s ) and diffuse (f d ) component matrices, is expressed as [14], where The jth row and lth column are indicated for the pBRDF matrices associated with the surface (M jl s ) and diffuse (M jl d ) scattering Mueller matrices.
The terms in f s that multiply the matrices M jl s R(α i ) represent the scalar BRDF as derived from the microfacet model.The slope variance σ 2 is a surface roughness parameter.θ denotes the angle between the microfacet and the macro-surface (mean) normal and is given by θ=cos −1 [(cosθ i +cosθ r )/(2cosβ)].β is the incident angle relative to the microfacet normal of interest and is defined by cos(2β)=cosθ i cosθ r +sinθ i sinθ r cosφ.The geometric attenuation factor (GAF), which models masking and shadowing effects arising from the surface variations, is included to keep the microfacet pBRDF bounded for large incident and reflection angles and is written as GAF = min[1, 2cosθ i cosθ/cosβ, 2cosθ r cosθ/cosβ] [24].The matrix R(α i ) is responsible for rotating an incident Stokes vector in the incident plane to the scattering plane.α i is the angle between the two planes and is defined by In the M jl matrix, * is the complex conjugate operation; r s and r p are the Fresnel reflection coefficients of the incident s-and p-polarizations, respectively.In terms of the pBRDF coordinate definition, the Fresnel reflection coefficients are given as [25] In the diffuse scattering pBRDF component (f d ), the Lambertian radiance is constant with angle (1/π) although the angular distribution of the observed intensity is cos θ r /π [24,26].M jl d describes the depolarization phenomena of the diffuse scattering, which is defined with no polarization preference in this case [9][10][11][12][13][14][15][16][17][18][19][20][21][22].A convenient reflectance parameter ρ d , based on Kubelka-Munk (KM) theory [27] and the work by Le Hors et al. [26], has been developed to describe the intensity ratio of the diffuse scattering to the scattering from a perfect Lambertian surface.In general, ρ d can be a function of the incident/reflected angles, complex refractive index, material absorptive properties and surface roughness.However, here ρ d is assumed to be independent of scattering angles and surface roughness.Previous studies support this assumption as the change in intensity for typical diffuse scattering over a limited range of illumination angles is small, especially for targets with rough surfaces [13][14]26,28].However, if a target surface is comprised of multiple thin, low-absorptive layers of different materials, then the penetrations and interactions of light among the layers may cause a polarization preference that is not modeled by M jl d in Eq. ( 2) and could also invalidate the assumption that ρ d is constant with angle.
If the Stokes values are measured and the source-target-receiver geometry is known, then Eq. ( 1) is a function of only the four physics-based parameters n, k, σ 2 and ρ d .Collecting multiple Stokes measurements throughout the day provides changing illumination and scattering geometries and leads to a nonlinear system of equations based on Eq. ( 1).If the measurement number is greater than four then the system of equations is overdetermined and the values of n, k, σ 2 and ρ d can be jointly estimated using an unconstrained, non-linear solving algorithm [29].
In essence, the estimation of the refractive index parameters relies on the polarimetric signature due to f s as a function of Sun position.However, the unpolarized signal of f d effectively modifies the result so we jointly estimate ρ d for the diffuse component from the data, which accounts for this effect.
We note that our model and estimation approach is general and other scattering models [30] can be incorporated in the pBRDF expressions.The approach can also be extended to include a partially polarized diffuse scattering component, compensate for atmospheric turbulence effects [31][32], or introduce illumination properties due to cloudy/overcast conditions [33] to improve the parameter estimation accuracy.

Data and Analysis
The analysis presented here utilizes Ground-MSPI data collected for two different target scenes where the images were recorded every 5 minutes under mostly clear sky conditions.The acquisition specifications for the datasets, including the approximate viewing direction for the camera, are shown in Table 1.The Ground-MSPI camera measurements are highly accurate with a measured DOP uncertainty less than 0.005 [16].Figures 2 and 3 present polarimetric image panels of the two scenes for the 660 nm wavelength collected around 10:00 in the morning.One scene ("sidewalk") includes part of a pathway with gravel and a concrete barrier and the other ("vehicle") includes a white car and a garbage container lying on its side.Each scene set is comprised of four images that correspond to the first three Stokes parameters (S 0 , S 1 and S 2 ) and the degree of polarization (DOP = [(S 1 ) 2 +(S 2 ) 2 ] 1/2 /S 0 ).The images are displayed in pseudo-color with the associated color scale shown to the right of each image panel.The DOP frames are included because they show high contrast between different objects in the scene.The sidewalk frames show obvious variations in the returned flux (S 0 , Fig. 2(b)) where the gravel walkway is slightly brighter than the ground and a strip on top of the small concrete barrier is more reflective than the rest of the scene.The barrier features also show vertical polarization differences (S 1 , Fig. 2(c)) as well as some 45°and 135°components (S 2 , Fig. 2(d)).For the vehicle scene frames, the reflected flux (S 0 , Fig. 3(b)) appears relatively uniform with a few accents on the car indicating more flux.The S 1 frame (Fig. 3(c)) shows a few features with a preference for vertical polarization (negative values) and the S 2 frame (Fig. 3(d)) shows other features with a small preference for 45°or 135°linear polarization.For each scene, we examined two subareas indicated with the yellow boxes in the S 0 images (second row of panels of Figs. 2 and 3), where A is the gravel walkway, B is part of the concrete barrier, C is the car hood and D corresponds to a flat section of the container that is approximately parallel to the ground.The curves to the right of each image panel show the associated mean Stokes values inside the selected subareas as a function of time for the 470, 660 and 865-nm wavelengths.The S 0 values are presented on a relative scale and the other Stokes values are normalized to S 0 .The peak values of the mean intensity (S 0 ) measurements appear near 13:00 in all cases.Some of the temporal variations of S 0 are due to intermittent clouds that create illumination changes.An example is the variations in curve C near 13:50 in Fig. 3(b).These same variations are also present in curve D in the same plot, although the display scaling does not allow this to be apparent.The mean S 1 curves for subareas B, C and D show a dip that roughly corresponds with the time of the peak S 0 curves.There are no significant trends for the S 1 and S 2 values of the gravel walkway (A).These Stokes measurement results as a function of time are generally consistent with results presented in previous studies [16][17][18][19][20][21].
To estimate the surface scattering parameters for the Ground-MSPI images, it is necessary to first determine the angles θ i , θ r and φ that are associated with each image pixel at each collection time.This requires knowledge of the position of the sun, the camera, the ground location associated with each image pixel and the direction of the mean surface normal at each pixel.The sun's position is calculated from the time of day and the camera ground position (latitude and longitude) is found with a GPS measurement.The latitude and longitude for each  pixel is derived the camera position, the known camera field-of-view and a determination of the viewing direction.
It is impractical to determine the exact mean surface normal direction for each pixel in the image, so it is assumed to be zenith (vertical direction opposite to the gravitational force, consistent with Fig. 1) for all pixels in the scene.We note that the estimation of the n and k values does not depend directly on the mean surface normal value whereas the estimation of σ 2 and ρ d does.However, the vertical mean normal assumption is still useful for estimating σ 2 and ρ d when individual features (e.g., grass blades) are small compared to the region of interest such that they are treated as small variations in the mean surface [1][2][16][17][18][19][20][21][22].In the future, the estimation accuracy of σ 2 and ρ d could be improved by incorporating any available information about the actual orientations of the mean surface normals relative to zenith.

Results
For each target scene, we selected eight sets of Stokes images (a total of twenty-four images) collected from 10:00 to 17:00 with a fixed time interval of one hour at each wavelength.These frames were used to estimate the values of parameters ρ d , n, k and σ 2 .The parameter values at each pixel for each wavelength are jointly estimated from the twenty-four images using Eqs.( 1)-( 6).The Levenberg-Marquardt nonlinear least squares algorithm [29] utilized for the parameter estimation is unconstrained.The initial guesses for n, k, ρ d and σ 2 were 1.5, 0, 0.5 and 0.5, respectively.We find that for a reasonable range of initial guesses (such as 0 to 5 for n) the standard deviation of the results changes slightly but there is little effect on the mean results.This behavior is consistent with modeling results from previous studies [7][8][12][13][14].Unreasonable estimate results (e.g.negative values of n, k, ρ d and σ 2 ) are neglected.In addition, some filtering and contrast enhancement processing have been applied (e.g., Filter2 and Wiener2 algorithms from MATLAB [34]) to the resulting parameter images to increase object definition and reduce noise without losing the original estimation characteristics.Finally, the mean and standard deviation values of the parameters within the selected areas (same as in Fig. 2(b) and 3(b)) are calculated as a function of wavelength.

Image Results
Figure 4 shows the images generated with the estimated parameter values of each pixel for sidewalk (left column) and vehicle (right column) at 660 nm.The top of the concrete barrier in sidewalk is distinguishable in the σ 2 image (Fig. 4(d)) and in the ρ d image (Fig. 4(c)).However, the estimated refractive indices of the concrete and gravely areas show little contrast, especially the imaginary parts (Fig. 4(a)-(b)).For the vehicle scene (Fig. 4(e)-(g)), the parameters values for the vehicle and the container clearly have different attributes from each other and from those of the background materials of concrete and gravel.The images generated at 470 nm and 865 nm produced similar results and, for the sake of brevity, are not shown here.In comparison with the Stokes and DOP images in Figs.2-3, the parameter-based images generally distinguish between the same items in the scenes (e.g., barrier, vehicle, and container).However, the pixel values now contain index of refraction, scattering and surface roughness information, which can aid in the analysis of the imagery.

Evaluation of Estimated Parameter Values
The curves to the right of each image panel in Fig. 4 show the estimated mean and standard deviation (error bars) of n, k, ρ d and σ 2 as a function of wavelengths for the selected areas in the two scenes.The estimated n of the concrete and gravel are close to 1.5, which is reasonable because their main composition is presumably silicon dioxide [35].The refractive index values for the container and car hood are in the 1.40 to 1.45 range and indicate a small imaginary part (k).The typical materials that would make up the container and the car hood are plastic and aluminum (or maybe carbon fiber) but it is likely that paint on these objects is actually being observed.In this case, the estimated values show some similarity to the paints studied previously [7,[13][14]23,36].
The σ 2 values in the sidewalk scene show a clear distinction between the concrete barrier components and the background, and the ρ d and σ 2 values show a distinction between the objects in the vehicle scene.The σ 2 values within the sample areas lie in the range consistent with a rough surface and satisfy the requirements of the microfacet assumptions [5][6].The estimated roughness value of the gravel is the largest and the value for the car hood is the smallest.The σ 2 values generally decrease with wavelength whereas the ρ d values increase.Note that the measured DOP in Figs.2(a) and 3(a) for the selected areas descend in value in the following order: (D) container, (B) concrete, (C) car hood and (A) gravel.The estimated results of ρ d and σ 2 are consistent with the DOP measurements since a smaller diffuse scattering component and surface roughness value should be associated with reduced depolarization behavior.Although the results are generally consistent with our expectations for the various materials and surfaces found in the scenes, future research on the accuracy of the approach would be aided by placing surfaces with known characteristics in the scenes.In this study, we were limited to the existing targets within the scenes.
To verify the consistency of the estimation process, the estimated parameter mean values for the subareas were inserted back into Eq.( 1) to model the Stokes parameter values as a function of time.Figure 5 shows the model points along with the corresponding average Stokes values measured by Ground-MSPI.The model values show a consistent fit with the measurement curves.In remote sensing applications, the estimation accuracy of the model can be affected by others factors not considered here such as non-Gaussian surface roughness statistics, atmospheric scattering, multiple scattering, and other illumination effects [1,28,[31][32][33].These issues are subjects for study in the future.

Conclusions
In conclusion, we propose and demonstrate the application of a modified pBRDF model with a diffuse scattering component to estimate the refractive indexes and scattering characteristics of materials from temporal sequences of Stokes polarimetric images of outdoor scenes.The images generated from the estimated surface parameters show distinctions between the same objects in the scenes as the DOP images.However, the parameter-based pixel values also provide index of refraction, scattering and surface roughness information that aids in the further analysis of the imagery.Estimated refractive index values for selected sample areas are consistent with expectations (concrete/gravel about 1.

Fig. 2 .
Fig. 2. Stokes images (660 nm) of the scene "sidewalk" alongside temporal plots of the mean Stokes values inside the selected subareas: subarea A (gravel walkway, -x) and B (concrete barrier, -circle)).Subarea locations are indicated in the second row S 0 images.

Fig. 3 .
Fig. 3. Stokes images (660 nm) of the scene "vehicle" alongside temporal plots of the mean Stokes values inside the selected subareas: subarea C (car hood, -plus) and D (container, -square).Subarea locations are indicated in the second row S 0 images.
5, and paints in the 1.40 to 1.45 range).Both the scattering parameters ρ d and σ 2 show variations depending on the object surface and the results are generally consistent with the expectation that a smaller diffuse scattering component and surface roughness value is associated with reduced depolarization behavior.The σ 2 values are generally found to decrease with wavelength whereas the ρ d values increase.The consistency of the estimation process was verified by inserting the estimated parameter values back into the model equations and demonstrating a close fit with the measured Stokes parameter curves as a function of time.