Most Mueller matrix imaging polarimeters implement sequential acquisition of at least 16 raw images of the same object with different incident and detected light polarizations. When this technique is implemented in vivo, the unavoidable motions of the subject may shift and distort the raw images to an extent such that the final Mueller images cannot be extracted. We describe a registration algorithm which solves this problem for the typical conditions of in vivo imaging, e.g. with spatially inhomogeneous medium to strong depolarization. The algorithm, based on the so called “optical flow”, is validated experimentally by comparing the Mueller images of a pig skin sample taken in static and in dynamic conditions.
©2007 Optical Society of America
As polarimetric imaging can provide contrasts quite different from ordinary intensity imaging, this technique may prove useful in many fields, and in particular for biomedical applications. It is suitable to discrimination of surface and beneath-the-surface contributions to backscattered light , and may provide an accurate delimitation of melanomas [2, 3] and other skin pathologies  as well as an evaluation of radiological burns .
A polarimeter comprises a Polarization State Generator (PSG) which encodes a given set of polarization basis states in the incident beam, and a Polarization State Analyzer (PSA) which analyzes the emerging light polarization by measuring its projections on another set of basis states. The actual number of basis states involved in the PSG and the PSA designs clearly depends on the completeness of the required polarimetric characterization of the sample.
The most complete such characterization is provided by the sample Mueller matrix M defined by the relation S out = MS in, where S in and S out are the Stokes vectors of the incident and emerging light. Experimental determination of M requires at least four linearly independent basis states for both the PSG and PSA, implying a minimal set of 16 measurements. Of course, this number can be reduced if only a partial knowledge of M is actually needed, or if some independent information is available for example, if the sample is known to be nondepolarizing, only 9 independent measurements are required. This is clearly not the case for most biological samples, for which the measurement of their complete Mueller matrix is necessary in a first step to fully assess their essential polarimetric properties. In a second step, the raw data acquisition may be reduced to the most relevant quantities.
Our imaging polarimeter is similar to that described in Ref. . The PSG consists of a linear polarizer followed by two nematic liquid crystal variable retarders, allowing a sequential generation of exactly four linearly independent S in states which form the column vectors of the modulation matrix W. The PSA is made of the same elements as the PSG in reverse order and it sequentially projects S out over exactly four Stokes vectors which constitute the row vectors of the analyzer matrix A. Mathematically, the 16 raw images taken during a complete measurement run can be expressed as a raw matrix image B related to the Mueller matrix image M by
As a result, M is easily extracted from B, provided A and W are known, i.e. the instrument is calibrated. More details on our acquisition and calibration procedures can be found in ref., and examples of Mueller matrix imaging of biopsies and other biological objects or phantoms (scattering media) are described in refs [5, 7, 8].
For polarimetric imaging with sequential acquisition of the raw data a mandatory condition for Eq. (1) to be safely inverted pixelwise is that the object is static during the whole acquisition cycle. For in vivo imaging, the unavoidable motions of the subjects would be “freezed” only if this cycle duration were shorter than 0.05 to 0.1 sec typically. Such an acquisition rate is currently difficult to achieve, due to a series of practical limitations (liquid crystal switching time, CCD integration time related to the CCD speed and available illumination level, data transfer and storage…).
Of course, this issue can be circumvented by implementing non-sequential raw data acquisitions, such as those based on division of amplitude  or spatial encoding of the polarimetric information on a high frequency carrier  with subsequent image analysis by FFT. However, the division of amplitude method may be expensive and lead to bulky setups (more detectors are required), and is thus best adapted to partial polarimetric measurements; the spatial encoding also is limited in practice to partial polarimetry and it reduces the useful spatial resolution to a carrier period.
Thus, sequential polarization modulation and analysis is still a valuable and practical approach for complete Mueller imaging, making it desirable to develop suitable registration algorithms for moving objects. Obviously, such algorithms must be defined according to the polarimetric characteristics of the objects under study and applied to the raw data before extracting M from B.
As shown below, for typical biological samples such as the skin samples presented in this study, this algorithm must meet the two following requirements:
- The motion of the imaged area may be complex, with arbitrary (but reasonably small) deformations;
- Different raw images may exhibit different histograms and contrast levels, with, however, similar small scale patterns.
A registration method, based on evaluation of “the optical flow” by a gradient approximation, is presented in section 2. Section 3 is devoted to the experimental validation of this method. Finally, a discussion and a few words of conclusion constitute the last section 4.
2. The Registration algorithm
2.1 Typical raw polarimetric images of a pig skin sample
Two raw images of the tattooed pig skin sample studied in this investigation are shown in Figure 1. The incident and emerging light polarizations were almost orthogonal for image (a) and almost identical for image (b). These images illustrate a very common behavior of such samples: the contrast, and, to a lesser extent, the average illumination clearly depend on the incident and detected polarizations, but the highly structured intensity patterns are essentially polarization-independent. This behavior directly stems from the absence of significant dichroism or birefringence. The sample polarimetric response is dominated by medium to strong depolarization (its Mueller matrix is essentially diagonal), with, however, a marked spatial dependence related to the surface local orientations which accounts for the easily recognizable patterns present on all raw images.
The intensity histograms of the images in Fig. 1 are shown in Fig. 2. The peak seen at the lowest levels in both histograms is clearly due to the tattooed zone. While many pixels exhibit comparable illumination levels in both images; the shiny pixels appearing in image (b) form an additional significant component on the high intensity side in the corresponding histogram. The images (a) and (b) shown in this example are those exhibiting the largest differences in this acquisition run. For the other 14 raw images the histograms are intermediate between those shown above.
2.1 Image registration by gradient-based optical flow
As the 16 raw images are characterized by different histograms and contrasts, the added information specifically contained in polarization actually constitutes a drawback for registration: the constant brightness assumption suitable to intensity based algorithms is no longer strictly valid. However, the similarity of small scale patterns suggests that the minimization of the intensity difference implemented in optical flow methods might still be adequate for registration, even though typical residual differences are expected to be larger than in cases where brightness is conserved. Moreover, as living subjects are not “rigid” we have to consider a non uniform transformation between different images.
Many well-known search strategies have been used in motion analysis and image registration problems [11–13] Most strategies involve a multiresolution method [14, 15] based on optical flow [16, 17]. Multiresolution allows better results when the displacement between two different images is large, but for our application it does not seem to be necessary.
Several optical flow techniques allow the determination of the 2-D motion between two images. More precisely, there are four main groups of algorithms:
- Gradient-based motion estimation;
- Block-based motion estimation;
- Model-based motion estimation;
- Frequency-domain motion estimation.
We chose a gradient-based motion estimation, as time-space differentiation induces a “stationarization” of the data which removes large scale patterns and enhances sensitivity to the small scale features which exhibit similar shapes in different raw images. In addition, to allow a non rigid transformation of the plane we chose a locally affine motion model.
2.2 The affine motion model
Six parameters are used for the estimation of the object motion:
where (p(x,y,t),q(x,y,t)) is the displacement induced on pixel (x,y) by the motion of the object between images at instants t and t + 1 . As the displacements are assumed small enough, for two consecutive images the intensity variation on the same physical spot of the moving sample can be expanded, to first order, as
where Ix, Iy and It are the partial derivatives of the intensity i(x,y,t). Moreover, with depolarizing objects such as skin samples, images taken with different polarizations typically exhibit variations in overall illumination and contrast, but the local intensity patterns remain quite similar. Thus the displacement is determined by minimizing the root mean squared intensity variation defined by Eq. (3) or, equivalently, the error function
over the region of interest, with the displacement defined by the affine model, which provides the following system of six linear equations for the six unknowns (a,b,c,d,e,f):
This procedure was applied either to the whole image or to a region of interest. Thus, we respectively called these approaches global affine registration or local affine registration.
3. Experimental validation
We first took “reference” polarimetric images of the sample in static conditions, with a pixelwise inversion of Eq. (1). The resulting Mueller matrix is represented on Fig. 3. As expected, this matrix is essentially diagonal. Then, the same acquisition procedure was run while the sample was moved, so as to mimic typical breathing motions. The speed was low enough to keep the raw images reasonably crisp, but sufficient to introduce clearly visible displacements from one raw image to the next. The same pixelwise treatment as for static data now provided the Mueller matrix presented on Fig. 4. This matrix clearly exhibits severe artifacts, especially in M 22 and off-diagonal terms, with “ghost” images of the tattoo.
To remove these artifacts and recover a Mueller matrix as close as possible to that shown in Fig. 3 we applied the registration scheme to readjust all the raw images to the first one (B 11). In a first step, the registration algorithm was applied on the entire images. In a second step, the displacement field was refined locally in a 10×10 pixel slipping window, providing a non uniform registration with improved accuracy. As a typical example, the displacement field for the registration of the two raw images shown in Fig. 1 is displayed in Fig. 5.
Finally, we applied the image displacement fields to the calibration matrices A and W. These matrices do exhibit significant spatial variations with the pixel position (x,y) over the field of view, meaning that the incident and emerging polarization states do depend on this position. Now, if the object to be studied is displaced by (pn,qn) with respect to its reference position (x,y) for the acquisition of the raw image n, the polarization states to be considered for this object and this raw image n are those defined by A and W at the point (x+pn, y+qn) and not (x,y).
After registration of all the raw images and calibration matrices, we eventually obtained the Mueller matrix shown in Fig. 6. For a better visualization of the details we also show the M 22 terms of the static, dynamic and the readjusted Mueller matrices on Fig. 7.
The improvement provided by the registration is clearly visible. More quantitatively, the correlation r, see Eq. (6), between the readjusted, B, and the static, A, images is close to 0.95 while for the “dynamic” one, this correlation factor r drops to 0.88 in this particular case.
We also evaluated the percentage of pixels for which the Mueller matrix is not a physical one, according to the well known criterion of the coherency matrix eigenvalues . For the static case, only about 0.1% of the pixels were not valid, as expected from the accuracy of our experimental setup. In the dynamic case, we found about 5% of unphysical pixels, which is really a lot for this type of depolarizing sample. The readjusted matrix featured about 0.1% non valid pixels, as the static one.
As it is based on the intensity derivatives of the images, the described registration method is well adapted to raw images with different histograms but similar small scale patterns. This point is particularly important in polarimetry where the contrast varies with the polarization of the input-output polarized light. Moreover, the affine displacement approach, implemented in two steps, does not require any additional parameter to determine the best non uniform transformation, which turns out to be very accurate.
A major drawback for practical use is computation time. The most demanding step, the local registration; typically requires 1.5 hours on a Pentium 4 (2 GHz, 1 Go RAM) for each raw image, adding up to one day for a complete Mueller matrix. However, this task could easily be speeded up at moderate cost by simply distributing the computation over several machines.
The limits of this method for other applications in polarimetric imaging are twofold. First, if two images feature no correlated gradients, the registration algorithm may not converge. This situation may occur, for example, with significant diattenuation or retardation, or with too sharp boundaries between zones with large differences in depolarizing power, making the gradient approximation questionable. For most biomedical applications, however, the relevant samples are typically partial depolarizers with diagonal Mueller matrices and “reasonable” spatial inhomogeneities, for which the proposed algorithm is adequate. Second, if the displacements between raw images are too large for the Taylor expansion to be valid, our algorithm must be run in a multi-resolution scheme. This can be easily achieved by a pyramidal approach , the image resolution is first degraded from 256×256 pixels to 16×16 pixels (for example) before applying a first registration, and then the resolution is gradually retrieved, with a registration being performed at each step.
We have presented in this paper a robust raw image registration method which allows in vivo acquisition of Mueller matrix images with virtually no loss in data accuracy with respect to exvivo imaging. This method may open new possibilities for polarimetric imaging, even though the computational burden currently excludes real time applications. Indeed, complete Mueller matrix imaging can be used for itself, or as an intermediate step to optimize a simpler polarimetric technique for each particular application. In both cases, the possibility to use data acquired in vivo is crucial. The interpretation of the data would greatly benefit from investigation of similar samples by polarization-sensitive optical coherence tomography [20,21], which would provide depth resolved investigation of the skin polarimetric properties. In a forthcoming paper we will present our study of the irradiation effect on pig skin.
References and links
2. M. H. Smith, P. Burke, A. Lompado, E. Tanner, and L. W. Hillman, “Mueller matrix imaging polarimetry in dermatology,” Proc. SPIE , 3911, 210–216 (2000). [CrossRef]
5. F. Boulvert, B. Boulbry, G. Le Brun, B. Le Jeune, S. Rivet, and J. Cariou, “Analysis of the depolarizing properties of irradiated pig skin,” J. Opt. A: Pure Appl. Opt. 7, 21–28 (2005). [CrossRef]
7. S. Manhas, M. K. Swami, P. Buddhiwant, N. Ghosh, P. K. Gupta, and J. Singh, “Mueller matrix approach for determination of optical rotation in chiral turbid media in backscattering geometry,” Opt. Express 14, 190–202 (2006). [CrossRef] [PubMed]
11. N. Ritter, R. Owens, J. Cooper, R. H. Eikelboom, and P. P. v. Saarloos, “Registration of stereo and temporal images of the retina,” IEEE Trans. Med. Image 18, 404–418 (1999). [CrossRef]
12. J. A. Bangham, R. Harvey, and P. D. Ling, “Morphological scale-space preserving transforms in many dimensions,” J. Electron. Imaging 5, 283–299 (1996). [CrossRef]
13. M. Kim, J. C. Jeon, J. S. Kwak, M. H. Lee, and C. Ahn, “Moving object segmentation in video sequences by user interaction and automatic object tracking,” Image Vis. Comput. 19, 245–260 (2001). [CrossRef]
14. J. M. Odobez and P. Bouthemy, “Robust multi-resolution estimation of parametric motion models applied to complex scenes,” IRISA, Technical Report , 788 (1994).
15. Z. Zhang and R. S. Blum, “A hybrid image registration technique for a digital camera image fusion application,” Information Fusion 2, 135–149 (2001). [CrossRef]
16. B. K. P. Horn and B. G. Schunk, “Determining optical flow,” Artif. Intell. 17, 185–203 (1981). [CrossRef]
17. J. A. Perrone, “Simple technique for optical flow estimation,” J. Opt. Soc. Am. A 7, 264–278 (1990). [CrossRef]
18. C. Brosseau, Fundamental of polarized light: a statistical optics approach (Wiley, NY, 1998).
19. P. J. Burt and E. H. Adelson, “The Laplacian Pyramid as a Compact Image Code,” IEEE Trans. Commun. , COM-31, 532–540 (1983). [CrossRef]
20. G. Yao and L. V. Wang, “Two-dimensional depth-resolved Mueller matrix characterization of biological tissue by optical coherence tomography,” Opt. Lett. 24, 537–539 (1999). [CrossRef]
21. M. Yamanari, S. Makita, V. D. Madjarova, T. Yatagai, and Y. Yasuno, “fiber-based polarization sensitive optical coherence tomography using B-scan-oriented polarization modulation method,” Opt. Express 14, 6502–6515 (2006) [CrossRef] [PubMed]