An FPGA-Oriented Algorithm for Real-Time Filtering of Poisson Noise in Video Streams, with Application to X-Ray Fluoroscopy

Abstract

In this paper we propose a new algorithm for real-time filtering of video sequences corrupted by Poisson noise. The algorithm provides effective denoising (in some cases overcoming the filtering performances of state-of-the-art techniques), is ideally suited for hardware implementation, and can be implemented on a small field-programmable gate array using limited hardware resources. The paper describes the proposed algorithm, using X-ray fluoroscopy as a case study. We use IIR filters for time filtering, which largely simplifies hardware cost with respect to previous FIR filter-based implementations. A conditional reset is implemented in the IIR filter, to minimize motion blur, with the help of an adaptive thresholding approach. Spatial filtering performs a conditional mean to further reduce noise and to remove isolated noisy pixels. IIR filter hardware implementation is optimized by using a novel technique, based on Steiglitz–McBride iterative method, to calculate fixed-point filter coefficients with minimal number of nonzero elements. Implementation results using the smallest StratixIV FPGA show that the system uses only, at most, the 22% of the resources of the device, while performing real-time filtering of 1024 × 1024@49fps video stream. For comparison, a previous FIR filter-based implementation, on the same FPGA, in the same conditions and constraints (1024 × 1024@49fps), requires the 80% of the logic resources of the FPGA.

Appendix

1.1 Derivation of Temporal Threshold

As stated in Sect. 2, p(i, j, n) ~ N[µ(i, j, n), σ²(i, j, n)], where, according to (1), σ²(i, j, n) is a function of µ(i, j, n). In order to identify whether p(i, j, n) belongs to a different moving object we can compare p(i, j, n) with the previous output of the T-filter, p_T(i, j, n − 1), which represents the most accurate information about the state of pixels at position (i, j) in previous frames. p_T(i, j, n − 1) is a weighted sum of p(i, j, n′), independent of each other, and belongs to the same object; therefore, p_T(i, j, n − 1) ~ N(µ_T(i, j, n − 1), σ ² _T (i, j, n − 1)). The problem of defining whether p(i, j, n) belongs to the same object of p_T(i, j, n − 1) is the statistical hypothesis testing problem:

$$ {\mathbf{H}}_{{\mathbf{0}}} :\;\mu \left( {i,j,n} \right) = \mu_{T} \left( {i,j,n - 1} \right) $$

(37)

The problem is solved by defining a threshold T_T(i, j, n):

$$ P\left[ {\left| {p(i,j,n) - p_{T} (i,j,n - 1)} \right| > T_{T} (i,j,n)\;\;|\;\;{\mathbf{H}}_{{\mathbf{0}}} } \right] = \alpha $$

(38)

where 1 − α represents the confidence level of the test. The two quantities p(i, j, n) and p_T(i, j, n − 1) are independent of each other; therefore, under hypothesis H₀ we have:

$$ p\left( {i,j,n} \right) - p_{T} \left( {i,j,n - 1} \right)\; \sim \;N\left[ {0,\;\sigma_{D}^{2} (i,j,n - 1)} \right] $$

(39)

where

$$ \sigma_{D}^{2} (i,j,n - 1) = \sigma^{2} (i,j,n) + \sigma_{T}^{2} (i,j,n - 1) $$

(40)

The variance σ ² _T (i, j, n − 1) of the IIR filter response is computed approximating the filter impulse response with a rectangular window of length M:

$$ p_{T} (i,j,n - 1) = \sum\limits_{k = 1}^{\infty } {h_{m} (k - 1)p^{\prime}(i,j,n - k)} \cong \frac{1}{M}\sum\limits_{k = 1}^{M} {p^{\prime}(i,j,n - k)} $$

(41)

Let us now name m(i, j, n − 1) ∈ [1, M] the number of frames from which the pixel p_T(i, j, n − 1) can be considered static (i.e., m(i, j, n − 1) = k in the case in which the last filter reset was at frame n − k). In this hypothesis, we have:

$$ p^{\prime}\left( {i,j,n - k} \right) = \left\{ {\begin{array}{*{20}l} {p(i,j,n - k)} \hfill & {k \le m(i,j,n - 1)} \hfill \\ {p(i,j,n - m(i,j,n - 1))} \hfill & {k > m(i,j,n - 1)} \hfill \\ \end{array} } \right. $$

(42)

In (42) p(i, j, n − k) are independent of each other. In addition, under hypothesis H₀, p(i, j, n − k) ~ N(µ(i, j, n), σ²(i, j, n)). From (41), (42) we obtain:

$$ \sigma_{T}^{2} (i,j,n - 1) = \sigma^{2} (i,j,n) \cdot \left( {\frac{{m(i,j,n - 1)^{2} }}{{M^{2} }} - \frac{m(i,j,n - 1)}{M}\left( {2 + \frac{1}{M}} \right) + \left( {1 + \frac{2}{M}} \right)} \right) $$

(43)

Substituting (43) into (40) we have:

$$ \sigma_{D}^{2} (i,j,n - 1) = \sigma^{2} (i,j,n) \cdot g\left( {m(i,j,n - 1)} \right) $$

(44)

$$ g\left( m \right) = \frac{{m^{2} }}{{M^{2} }} - \frac{m}{M}\left( {2 + \frac{1}{M}} \right) + \left( {2 + \frac{2}{M}} \right) $$

(45)

From (39), the solution of (38) is:

$$ T_{T} (i,j,n) = \text{CDF}_{N}^{ - 1} \left( {1 - \frac{\alpha }{2}} \right) \cdot \sqrt {\sigma^{2} (i,j,n) \cdot g\left( {m(i,j,n - 1)} \right)} $$

(46)

where CDF ⁻¹_N (p) is the inverse normal cumulative distribution function.

1.2 Derivation of Spatial Threshold

By following the same approach used for T-filters, the generic pixel p_T(i′, j′, n) is considered in the mean operation only if we statistically infer that the pixel belongs to the same object of pixel p_T(i, j, n). We therefore impose another statistical hypothesis testing problem:

$$ {\mathbf{H}}_{{\mathbf{0}}} :\;\;\;\mu_{T} \left( {i^{\prime},j^{\prime},n} \right) = \mu_{T} \left( {i,j,n} \right) $$

(47)

The problem is solved defining a threshold T_S(i, j, i′, j′, n):

$$ P\left[ {\left| {p_{T} (i^{\prime},j^{\prime},n) - p_{T} (i,j,n)} \right| > T_{S} (i,j,i^{\prime},j^{\prime},n)\;\;|\;\;{\mathbf{H}}_{{\mathbf{0}}} } \right] = \alpha $$

(48)

Note that, according to (48), the threshold T_S is a function not only of (i, j), but also of (i′, j′). The variance of variable p_T(i, j, n) − p_T(i′, j′, n), in fact, depends on the variance of both p_T(i, j, n) and p_T(i′, j′, n), and these variances depend on m(i, j, n) and m(i′, j′, n), respectively (see (43)). However, with the help of a large number of simulations, we have verified that it is possible to use the same threshold for all pixels p_T(i′, j′, n) considered by S-filter(i, j) by assuming σ ² _T (i′, j′, n) = σ ² _T (i, j, n). We will therefore use a threshold T_S(i, j, n) defined as:

$$ T_{S} \left( {i,j,n} \right) \triangleq \;\left. {T_{S} \left( {i,j,i^{\prime},j^{\prime},n} \right)} \right|_{{\sigma_{T}^{2} \left( {i^{\prime},j^{\prime},n} \right) = \sigma_{T}^{2} \left( {i,j,n} \right)}} $$

(49)

In this case, we have:

$$ T_{S} \left( {i,j,n} \right) = {\text{CDF}}_{\text{N}}^{ - 1} \left( {1 - \frac{\alpha }{2}} \right) \cdot \sqrt {2 \cdot \sigma_{T}^{2} \left( {i,j,n)} \right)} $$

(50)

where σ ² _T (i, j, n) is given by (43). By considering also (44) and (40), we also have:

$$ T_{S} \left( {i,j,n} \right) = {\text{CDF}}_{\text{N}}^{ - 1} \left( {1 - \frac{\alpha }{2}} \right) \cdot \sqrt {2 \cdot \sigma^{2} (i,j,n) \cdot \left( {g\left( {m(i,j,n)} \right) - 1} \right)} . $$

(51)