Modern biomedical research relies on mathematical tools to process and analyze the acquired data. Therefore, it is important to fully understand the data analysis approaches and to provide proper interpretation of each processing step. Here, we shall analyze and improve a recently proposed analysis scheme.

Intensity measurements in noninvasive, noncontact imaging can be curvature biased due to the object’s shape. Image correction is necessary to enable quantitative physiological measurements and several attempts to assess shape and curvature have been proposed.^1,2,3,4 Kainerstorfer et al. developed in Ref. 5, a numerical scheme to remove curvature bias without measuring the object’s shape directly. The approach is based on assumptions on the separation of the physiological signal and the background component. Data correction is proposed through iterative averages of univariate polynomial fits of image rows and columns. This computational scheme is also applied in Refs. 6 and 7. Although polynomial fitting and averaging may have some smoothening effects, the separate application to rows and columns does not allow any direct interpretation of the fitting by means of standard low- or high-pass filtering. The aim of this method is to better understand the proposed scheme in Ref. 5 from a signal and image processing point of view and to perform some modifications that enable significantly wider applications. Indeed, we shall propose a replacement with a simpler, faster, and more intuitive scheme that allows for immediate interpretation of its effects on the data. Our contribution enables modifications of the scheme that allow to significantly weaken the required assumptions on the separation of background and physiological signal.

More precisely, the approach proposed in Ref. 5 fits each horizontal and vertical pixel line to one-dimensional (1-D) polynomials. Outcomes are then averaged pixelwise and the fits are repeated. Iterations are stopped after a fixed number of steps. In this article, we mathematically verify that the scheme converges toward a standard two-dimensional polynomial fit of the underlying image. Thus, we do not need any iterations and have a proper interpretation of the fitted image as the bivariate polynomial fit corresponds to a spatial low-pass filter. Moreover, using irregular grids, we can now additionally perform the scheme when parts of the image are missing or need to be excluded from the analysis. For instance, if regions of interest can be spatially separated from the background parts as in Ref. 8, then the use of irregular grids can greatly reduce inaccuracies. Such modifications are not possible using the scheme in Ref. 5.

We shall introduce a few mathematical concepts helpful to our further analysis. The fitting scheme in Ref. 5 is then formulated using these concepts.

Let $u={[u(x_1),\dots ,u(x_n)]}^{\top }\in {\mathbb{R}}^n$ be a data vector at the locations ${(x_j)}_{j=1}^n\subset \mathbb{R}$. Suppose that we want to fit $u$ by a polynomial of degree $k$. This can be formulated as the least-squares problem $$\underset{a\in {\mathbb{R}}^{k+1}}{\min }{\parallel u-Xa\parallel }^2,$$where $a={(a_j)}_{j=0}^k\in {\mathbb{R}}^{k+1}$ is the coefficients of a polynomial $p(x)=\sum _{j=0}^k{a}_j{x}^j$ and $X$ denotes the Vandermonde matrix of ${(x_j)}_{j=1}^n$, so that

We implicitly suppose that $n$ is bigger than $k$ and that ${(x_j)}_{j=1}^n$ is chosen in sufficient general position such that the columns of $X$ are linearly independent. The minimizer $\widehat{a}$ satisfies the normal equation $X^{\top }X\widehat{a}=X^{\top }u$ and is given by $\widehat{a}=X^{†}u$, where $X^{†}≔{(X^{\top }X)}^{-1}{X}^{\top }$ is the pseudoinverse of $X$. Thus, the term $X{X}^{†}u=X\widehat{a}$ is the polynomial least-squares fit.

To discuss bivariate polynomial fitting, we need to introduce the Kronecker product $A\odot B$ of two matrices $A\in {\mathbb{R}}^{d_1\times d_2}$ and $B\in {\mathbb{R}}^{d_3\times d_4}$, i.e., for $A={({\alpha }_{i,j})}_{i,j}$ we define the block-matrix $A\odot B≔{({\alpha }_{i,j}B)}_{i,j}\in {\mathbb{R}}^{d_1{d}_3\times d_2{d}_4}$. The relations

hold, where $C$ and $D$ are the matrices with compatible dimensions, i.e., $C\in {\mathbb{R}}^{d_2\times d_5}$ and $D\in {\mathbb{R}}^{d_4\times d_6}$.

Next, we apply the Kronecker product to describe bivariate polynomial fitting. An arbitrary bivariate polynomial $p$ of total degree $(k,l)$ in the unknowns $x$ and $y$ can be written as $$p(x,y)=\sum _{i=0}^k\sum _{j=0}^l{a}_{i,j}{x}^i{y}^j,$$where $a_{i,j}$ is the coefficients. Given a data matrix $u(x_i,y_j)$ on a grid $(x_i,y_j)$ for $i=1,\dots ,n$ and $j=1,\dots ,m$, the associated bivariate Vandermonde matrix is $Z=X\odot Y$, where $X$ and $Y$ are the underlying univariate Vandermonde matrices of size $n\times (k+1)$ and $m\times (l+1)$, respectively. We flatten the data matrix $u(x_i,y_j)$ into a vector $u={[u(x_1,y_1),\dots ,u(x_1,y_m),u(x_2,y_1),\dots ,u(x_n,y_m)]}^{\top }\in {\mathbb{R}}^{mn}$. Fitting a bivariate polynomial of degree $(k,l)$ to the data matrix $u(x_i,y_j)$ can be expressed in terms of $Z{Z}^{†}u$. If $E_m$ denotes the $m$-dimensional identity matrix, then using Eq. (2) yields

By denoting $u_i={[u(x_i,y_1),\dots ,u(x_i,y_m)]}^{\top }$, we obtain

which means that we fit univariate polynomials of degree $l$ to each row of the data matrix ${[u(x_i,y_j)]}_{i,j}$. The decomposition of the bivariate polynomial fit in Eq. (3) allows us to take the univariate fit of the rows in Eq. (4) as a new data matrix and fit univariate polynomials of degree $k$ to the columns. This yields the bivariate polynomial fit $Z{Z}^{†}u$ of degree $(k,l)$ in Eq. (3).

The iterative curve fitting scheme in Ref. 5 fits polynomials of degree $k$ to the columns of the original data matrix. It then fits polynomials of degree $l$ to the rows of the original data matrix. The column fit and row fit are averaged pixelwise resulting in a new data matrix, in which the scheme is repeated. Few iterations seem sufficient to derive visually acceptable results. To state the iterative fitting scheme in mathematical terms, we shall use the pseudoinverse and the Vandermonde matrices $X\in {\mathbb{R}}^{n\times (k+1)}$ and $Y\in {\mathbb{R}}^{m\times (l+1)}$. Performing one iteration is represented by the application of the matrix

By using the relations $X{X}^{†}X{X}^{†}=X{X}^{†}$ and $Y{Y}^{†}Y{Y}^{†}=Y{Y}^{†}$, which can be deduced from the definition of the pseudoinverse, we obtain for the second and third iterations

After $r$ iterations we end up with Hence, $L^r$ converges towards $Z$ when the number of iterations $r$ tends to infinity. Since the deviation $\Vert L^r-Z\Vert =\frac{1}{2^{r-1}}\Vert L-Z\Vert$ decreases exponentially in $r$ for any matrix norm $\Vert \cdot \Vert$, the term $L^{r}u$ rapidly approximates $Zu$. In other words, the iteration proposed in Ref. 5 essentially leads to a bivariate polynomial fit of the data matrix after few iterations already. Bivariate polynomial fits are associated with the low-pass filtering of the image.

As we have understood that iterative polynomial fitting essentially acts as a low-pass filter, we can now explicitly settle the assumptions needed in Ref. 5 for an accurate image curvature correction. The approach is based on the separation of the physiological foreground signal and the background component. It was mentioned in Ref. 5, that the background is supposed to induce higher intensity values than the physiological features in the foreground, so that the background dominates the polynomial fitting. Moreover, the foreground must correspond to higher frequency components in the image, so that the low-pass filtering captures the background and avoids the foreground. If one of the two assumptions is violated, then inaccuracies can occur.

Often, physiological features may lead to intensity values that are dominating the background, so that the assumptions in Ref. 5 are violated. If features can be separated spatially from the background and cover only a relatively small area of the image, for instance, if features are in the image center only, then our analysis in the previous section becomes useful because it suggests some advantageous modifications. The iterative fitting proposed in Ref. 5 requires a regular pixel grid. In contrast to the iterative fitting, the bivariate polynomial approach can also be performed on irregular grids, i.e., the pixel grid of an image, in which we removed some pixel areas. The use of irregular pixel grids is beneficial when the physiological features are spatially well separated, so that they can be excluded (masked) from the image. Foreground features can have dominating intensity values, but the bivariate polynomial is determined by the irregular grid that avoids those physiological parts. The computed bivariate polynomial can be evaluated at all pixels of the image, including the masked ones. If the background corresponds to low-frequency image components, then it will also be well approximated by this bivariate polynomial at the masked pixels. Therefore, the polynomial fitting on irregular grids can increase accuracy, especially when measurements are needed to provide quantitative physiological information.

Next, we illustrate our findings by providing some numerical results that are, in fact, numerical examples of the theory outlined in the previous section. We shall verify the advantage of our proposed fitting on irregular grids on some simulated intensity measurements. In our example, both assumptions of the iterative polynomial fitting are violated. The intensity values of the foreground dominate the image and the spatial variation is smaller than the background. As the foreground covers only a relatively small part of the image in Fig. 1, our proposed method for irregular grids can be used. We choose the polynomials degree as 7. The mean intensity of the background is $\approx 6.9$ and the signal has constant intensity $\approx 14.6$. The iterative fitting necessarily applied to the entire image yields pixels with errors up to $\approx 6.1$, which is a relative error wrt the maximal intensity of $\approx 41.5\%$. For the bivariate fit, we first mask the signal area, hence compute a bivariate polynomial using only unmarked pixels, and then evaluate this polynomial at the entire pixel grid including masked pixels. This approach yields pixel errors up to $\approx 0.9$, which is a relative error wrt the maximal intensity of $\le 6.1\%$. Hence, as suggested by our theoretical results, if certain signal-to-background separations are not satisfied, then iterative fitting can cause inaccuracies that can be avoided by using irregular grids.

Although several methods to assess the shape and curvature of objects have been proposed in the biomedical field, the indirect curvature correction based on a computational model and the abandonment of additional measurements in Ref. 5 clearly has some benefits.

We have verified that the iterative curve fitting used in Ref. 5 for curvature correction in multispectral image sets converges toward standard bivariate polynomial fitting of the images. Bivariate polynomial fits are associated with the low-pass filtering, so that our analysis enabled a proper interpretation of this scheme. Moreover, no iterations are needed, so that we were able to simplify the computational approach. If iterations run up to $k$, then our simplification is $k$ times faster, but runtime is not an issue in either case, because 1-D polynomial fits can be computed quickly. The gain of our analysis is that it enabled us to modify the scheme and weaken the required assumptions for successful curvature correction. The authors in Ref. 5 needed two assumptions on the separation between foreground and background. It is assumed that background intensities dominate intensities from physiological features, so that the iterative fitting is driven by the background. Moreover, background needs to correspond to the lower-frequency and foreground to the higher-frequency components. Such assumptions may not hold in many types of optical measurements. Often though, the curvature correction based on the iterative fitting still yields a decent qualitative outcome.^5–8 If measurements are needed to provide quantitative information and one of the assumptions is violated, then the proposed iterative polynomial fitting in Ref. 5 causes inaccuracies and needs to be replaced with more advanced and computationally more expensive tools.

By using our new derivations, we can weaken the required assumptions for curvature correction, in particular, the signal-to-background ratio is not required to be small as long as the signal covers only a small spatial area. The iterative curve fitting in Ref. 5 requires a full image, so that pixel locations form a regular grid. The bivariate polynomial fit can be performed on any type of irregular grid. For instance, the assessment of certain types of skin cancer through multispectral imaging deals with the small patches of target signals.⁸ Cancer area may absorb strongly in certain wavelengths matching the required signal-to-background ratio, but wavelengths with low absorbance of the target area and hence higher image intensities do not match the requirements. Then, the cancer area can be excluded from the fitting and the remaining shape can be fitted through the bivariate polynomial approach. Thus, our modification increases the applicability of the polynomial fitting approach.

We have illustrated the potential of bivariate polynomial fits improving accuracy for simulated intensity measurements. We must mention though that bivariate polynomial fitting on irregular grids is computationally more expensive and can even lead to numerical instabilities with growing grid sizes. Nonetheless, as image processing is usually performed off-line, the increase in accuracy may outweigh the growing runtime. Taking care of stability issues for the larger irregular grid sizes goes beyond the scope of the present article, and so we refer to standard textbooks in numerical analysis. Despite the increase in runtime, our results enable a weakening of required assumptions for the curvature correction that is useful for a larger field of applications.