Wavelet-based noise-model driven denoising algorithm for differential phase contrast mammography

: Traditional mammography can be positively complemented by phase contrast and scattering x-ray imaging, because they can detect subtle differences in the electron density of a material and measure the local small-angle scattering power generated by the microscopic density ﬂuctuations in the specimen, respectively. The grating-based x-ray interferometry technique can produce absorption, differential phase contrast (DPC) and scattering signals of the sample, in parallel, and works well with conventional X-ray sources; thus, it constitutes a promising method for more reliable breast cancer screening and diagnosis. Recently, our team proved that this novel technology can provide images superior to conventional mammography. This new technology was used to image whole native breast samples directly after mastectomy. The images acquired show high potential, but the noise level associated to the DPC and scattering signals is signiﬁcant, so it is necessary to remove it in order to improve image quality and visualization. The noise models of the three signals have been investigated and the noise variance can be computed. In this work, a wavelet-based denoising algorithm using these noise models is proposed. It was evaluated with both simulated and experimental mammography data. The outcomes demonstrated that our method offers a good denoising quality, while simultaneously preserving the edges and important structural features. Therefore, it can help improve diagnosis and implement further post-processing techniques such as fusion of the three signals acquired.


Introduction
Breast cancer is the most ubiquitous cancer in women. Its detection in early stages increments the survival rate and helps improve the life quality of the patient. Nowadays, absorption-based mammography is the most frequently used detection technique. However, it is not able to detect some important malignancies, due to small differences in x-ray attenuation between them and healthy tissue [1].
Phase-contrast and scattering-based x-ray imaging can provide complementary information to the conventional absorption-based method, since they are able to detect subtle electron density differences and measure microscopic density variations in the sample. The acquisition of these types of images can be easily carried out using a grating-based x-ray interferometer, which is able to simultaneously generate absorption, differential phase contrast (DPC) and scattering (DCI) images of the sample and, in addition, performs well with conventional x-ray sources. The first research work with native, non-fixed whole breast samples, including regular and cancerous breast tissue, was presented recently by our group [2]. In that work, a differential phase contrast mammography (mammoDPC) demonstrator based on a conventional x-ray tube was designed, constructed and used to image whole native breast samples directly after mastectomy. In practice, the recorded DPC and scattering images show a significant amount of noise, compared to the absorption signals, because not as many signal averages as in the latter case are acquired. This noise must be removed, for better visualization, diagnosis and further post-processing, such as the image fusion of the three signals [2].
In the last two decades, the wavelet transform has gained popularity as a promising denoising technique, due to properties such as sparsity and multi-resolution structure. The discrete wavelet transform generates a scale-space representation of a signal, by successively low-pass and highpass filtering it [3].
The wavelet-based denoising approach operates as follows: 1) The noisy image is wavelet transformed. 2) Thresholding is applied to the wavelet coefficients in absolute value, assuming that the smallest coefficients correspond to noise. The coefficients above the threshold are left untouched, when using a hard-thresholding method, or shrinked by the value of the threshold, when sticking to a soft-thresholding strategy, which in practice produces more visually agreeable images over the former, as the hard-thresholding method is discontinuous and produces abrupt artifacts. 3) An inverse wavelet transform is applied [4].
Several methods, based on the orthogonal wavelet transform, have been proposed for threshold selection. Donoho and Johnstone [5] introduced the universal threshold VisuShrink, which is the optimal threshold in the asymptotic sense and minimizes the L 2 norm of the difference between the unknown noiseless image and its thresholded version. Although this global threshold might be a good starting point when little is known about the signal, its performance is usually insufficient, because it lacks spatial adaptivity, kills many useful signal coefficients, causing image oversmoothing [6], and yields visual artifacts [7].
Later on, Chang, Yu and Vetterli brought in BayesShrink [8]. This threshold is derived in a Bayesian framework, is adaptive to each wavelet decomposition subband and depends on data-driven estimates of the parameters. More recently, Chen, Bui and Krzyzak [9] introduced NeighShrink, which thresholds the wavelet coefficients by considering the magnitude of the sum of the squares of the wavelet coefficients within a pre-defined neighboring window, under the assumption that a large wavelet coefficient will have large coefficients as its neighbors. BayesShrink and NeighShrink certainly outperform VisuShrink, mainly because they are subband-adaptive.
Although not as robust as the wavelet-based denoising methods, there have been some useful spatial-domain based denoising strategies. The first one is the Wiener filter, a method that relies on the assumption that the noisy image f can be expressed as the sum of the noiseless image x and the noise n, which are both stationary Gaussian processes. It aims to find a linear estimate of x that minimizes the mean square error. When x is also a white Gaussian process, the filter has a very simple scalar form that only depends on the noise variance and local estimators of the signal variances and means, which can be easily computed [10]. Chambolle et al [11] proposed a denoising algorithm based on the total variation minimization. They defined the total variation as a weighted sum of a term favoring smoothness, based on the local derivatives, and a term accounting for the fidelity to the input image. Through a fixed point algorithm, they were able to find the optimal weighting factor and compute the noiseless image simultaneously. Even though their approach results certainly interesting from the theoretical point of view, it has shown to yield oversmoothed outcomes.
The denoising methods described above can be classified as "blind" methods, because they assume that the noise distributes in a Gaussian manner along the whole image and use certain estimators to calculate the variance of this noise. However, in mammoDPC, the noise models of the absorption, DPC and scattering signals have been studied and it has been found that the noise variance, instead of being the same for the whole image, varies in a pixelwise manner and can be calculated [12]. Therefore, a denoising method considering these pixelwise noise maps is expected to be more reliable than one assuming uniform Gaussian noise.
One potential choice to solve this denoising problem is the spatial-domain method proposed by Kisilev, Shaked and Hwan Li [13], which claims that the visual perception of noise depends on both the noise and the underlying image content or signal activity. Based on this argument, known as the masking effect, they introduced a pixelwise noise threshold which is proportional to the ratio of the noise standard deviation to the signal activity. The signal activity was calculated by using local directional derivatives. Another option to solve this particular problem is the Wiener filter, violating the assumption of uniform Gaussian noise and using a different noise variance value for each pixel; a more robust strategy is a wavelet-domain based method incorporating the available noise model, because this transform has shown remarkable advantages. For the design of this method, two issues have to be addressed: 1. The noise variance must be translated from the spatial domain to the wavelet domain. 2. A pixelwise threshold in the latter domain has to be defined.
In this work, a pixelwise wavelet-based noise model driven denoising algorithm (WND) for mammoDPC is proposed. The translation of the noise variance from the spatial domain to the wavelet domain is performed through a modified version of the method proposed in [14]. The pixelwise threshold in the wavelet domain is then derived in a Bayesian framework by modeling the so called "context" of each wavelet coefficient, based on the spatially adaptive wavelet denoising method proposed by Chang et al [15]. Along with this innovative proposal, the two spatial-domain noise-model based possibilities described in the previous paragraph are carried out. For comparison purposes, the VisuShrink, BayesShrink and NeighShrink wavelet-based blind denoising methods, as well as the widely employed median filter, are implemented and tested.The available noise model is used to add different levels of noise to DPC and DCI images in a simulated scenario. Additionally, images of a breast sample are denoised for comparison. The Mean Square Error (MSE) and the Complex Wavelet Structural Similarity Index (CWS-SIM) [17] are utilized to assess the outcomes of the denoising methods with the simulated data set. For the mammoDPC images, the background SNR is calculated and some intensity profiles are analyzed.
The results obtained show that the method proposed manages to denoise the images, while preserving the edges and improving the image quality simultaneously. Therefore, this algorithm can certainly be useful for additional post-processing, visualization of structures and diagnosis.

Signal acquisition
A differential phase-contrast mammography prototype [2] was used for signal acquisition in the present work. This prototype consists of a Talbot-Lau interferometer, that is composed of a Seifert ID 3000 x-ray generator, an unfiltered tungsten line focus tube (operated at 40 kVp with mean energy of 28 keV and a current of 25 mA), a 3-grating interferometer (with periods p 0 = 14µm, p 1 = 3.5µm and p 2 = 2µm) and a Hamamatsu C9732DK flat panel CMOS detector with a 12cm × 12cm field of view and a 50µm × 50µm pixel size. The interferometer was tuned at the fifth Talbot distance, yielding an angular sensitivity of 0.15 µrad. The distance between the source and the G1 grating was 140 cm and the geometry could be considered parallel; thus, flat gratings were employed. The resulting field of view, limited by the size of the G1 and G2 gratings, was 5cm × 5cm, so to get an image of the whole breast various acquisitions were stitched together. The exposure time was set to 9 seconds. To separate absorption, phase, and scattering signals, a phase-stepping approach, where G2 is translated perpendicular to the grating lamina by a fractional distance of the grating period, was followed. For this study, minimum 8 phase steps were applied, resulting in a total exposure time of at least 72 seconds for the 5cm × 5cm field of view. This fact explains the very low noise level in the absorption images.

Noise calculation
The noise model of the absorption, DPC and scattering signals was investigated by Revol et al [12]. The detector quantum noise and the phase-stepping jitter noise are the two main noise sources in differential phase contrast imaging. In our case, the phase stepping jitter noise is insignificant because the mechanical uncertainty of the piezo motor is in the order of nm [16], which yields a noise variance that is around 1% of that caused by quantum noise. Therefore, we take the quantum noise as the main noise source.
The detector quantum noise variance was found to be linearly proportional to the correspondent mean intensity: The slope f 1 has to do with the signal and noise transfer of the incoming x-rays to the output in arbitrary digital units. Due to the beam hardening effect, f 1 is different for the reference (flat field) scan ( f r 1 ) and the sample measurement ( f s 1 ). The variances σ I translate into errors in the transmission(T ), DPC and dark-field (V ) images σ T,det , σ DPC,det and σ V,det , respectively, that can be computed employing the error propagation formula and Eqs. (2)-(4): where: T : Mean intensity of the transmission image. V : Mean intensity of the dark field image. a r 0 : Mean intensity of the reference measurement. N ps : Number of phase steps acquired over one period p 2 . f r 1 : f 1 for the reference measurement. f s 1 : f 1 for the sample measurement. v r : Visibility of the phase stepping curve for the reference measurement. The dark field signal is computed as the ratio of the visibility of the phase stepping curve in the reference measurement to its visibility in the sample measurement.
To calculate f r 1 and f s 1 , 50 images were acquired for a fixed tube voltage and 13 different anode currents, ranging from 1 mA to 25 mA in steps of 2 mA, in the presence and absence of a 1mm Al filter. The mean and variance of the recorded intensity histogram for both the filtered and unfiltered cases and each anode current, were approximated in a pixel-wise manner. Afterwards, the average mean intensity and average variance were computed over a 100 × 100 pixel region. Finally, a linear regression was performed. The values found for f r 1 and f s 1 in our system were 0.33 and 0.27, respectively.

Wavelet based Noise model driven denoising algorithm (WND)
This method translates the noise variance from the spatial to the wavelet domain and uses the latter to define a pixelwise soft threshold.

Noise variance translation from the spatial to the wavelet domain
Let us assume that σ s (t) corresponds to the known noise standard deviation of pixel t in the spatial domain calculated using Eqs. (3) and (4). Therefore, the noise signal n(t), having that the noise can be assumed pixelwise independent in mammoDPC, can be expressed as: where ε(t) is white Gaussian noise. Take A(ω) as the DTFT of σ s (t), so that: Thus, combining Eqs. (5) and (6) and from [19], with F(ω) the DTFT of ε(t). The variance σ 2 s (t) can be computed as [19]: Let A * (ω) denote the complex conjugate of A(ω). Due to the fact that n(t) is real-valued, and since ε(t) is a white noise process, F(ω) has orthogonal increments. Therefore: Taking Eq. (11) into Eq. (10), we have: Now, let H (s,o) (ω) be the frequency response of the wavelet filters at scale s and orientation o, and H * (s,o) (ω) be its complex conjugate. Then, using a derivation similar to that for Eq. (12), the wavelet domain noise variance is approximately given by [14]: In this way, the connection between the spatial noise variance, σ 2 s (t), and the noise variance in the wavelet domain, [σ (s,o) wavelet (t)] 2 is set up by Eqs. (6), (12) and (13). For simplicity, from now on, the index (s, o) will be omitted.

Thresholding
Each wavelet coefficient is assumed to follow a Generalized Gaussian Distribution (GGD) corrupted with noise; this GGD has a standard deviation σ x (t). Therefore, the following soft threshold is a good approximation to the optimal threshold [15]: For the calculation of σ x (t), context modeling, a widespread procedure employed to adapt the coder to varying image characteristics in the field of image compression, is used [15]. The context Z [t] of each wavelet coefficient Y (t) is defined as a weighted average of the absolute value of its p neighbors. If these p neighbors are placed in a vector u t , then the context becomes: The weighting factor w is found through least squares: where U is a N × p matrix with each row being u T t for all t and Y is the N ×1 vector containing all the wavelet coefficients Y (t) belonging to a specific wavelet sub-band (N is the number of coefficients per sub-band) In this case, nine neighbors were considered, eight of them taken from the same subband and the other one from the parent subband [15].
Afterwards, the variance of the random variable is estimated from other coefficients whose context variable are close in value to Z [t] , as follows: where B t comprises the coefficients correspondent to the L closest points above Z [t] , the coefficients correspondent to the L closest points below and Y (t). The term σ 2 wavelet (t) shall be subtracted, because {Y (t)} are noisy coefficients and the variance is independent of the noise.
The proposed method can be used either with the orthogonal or the undecimated wavelet transform. However, the former has been found to yield images with evident ring artifacts. Therefore, in this work, the latter was used. Since the wavelet coefficients produced by this kind of transform are correlated, it is necessary to decorrelate them beforehand. For the sth level of decomposition, the coefficients can be separated into 2 2s sets of uncorrelated coefficients. These sets are given by Y [2 s t + k 1 ] t , k 1 = 0, 1, ..., 2 s − 1 [15].

Data sets
To evaluate the proposed denoising methods, two data sets were used: a. Simulated data set: To generate the simulated images (Figs. 1, 2 and 3), a 1200 × 1200 Matlab Shepp-Logan phantom was used. The noiseless dark field and tranmission images, V and T respectively, were produced by giving different intensity values to the distinct phantom regions. The noiseless DPC image, DPC, was computed by taking the unidirectional gradient of a phase contrast image which was created by giving a variant phase (between 0 and 2π) to the phantom. Spatially variant noise was added to each pixel t of images V and DPC as follows: where randn(t) represents a pseudorandom value drawn from the standard normal distribution. To compute σ DPC,det and σ V,det , Eqs. (3) and (4) were used. The visibility value (v r ) was set to 10 % for the pixels inside a mask corresponding to the center of the image, and to 8 % for the pixels located in the periphery; the number of phase steps (N ps ) was set to 4, to get an appreciable level of noise; the values found for f r 1 and f s 1 in our system were 0.33 and 0.27, respectively; V and T were generated as explained in the previous paragraph; the lowest intensity value of the mean intensity of the reference measurement, a r 0 , was set to 1000 and increased in 1000-steps, in order to obtain five DPC and DCI images, each with a different σ DPC,det and σ V,det , and a 50 % decrease of the noise level per step with respect to the previous one. Although it could seem more natural to add noise to both the reference and sample measurements, the fact that Eqs. (3) and (4) already take into account the noise propagation, allows us to generate the noisy images as presented above.
b. Mammography images: Images of a human breast specimen in the Craneo-caudal (CC) position were acquired. The sample preparation and the imaging protocol are described in [2] and section 2.1 of this paper, respectively.

Performance metrics
To evaluate the performance of the implemented denoising methods in the simulations, the Mean Squared Error(MSE) and the Complex Wavelet Structural Similarity (CWSSIM) index were used.
The MSE between two images x and y with the same number of pixels is defined as: where M is the number of pixels of the images. The CWSSIM index is an extension of the Structural Similarity index (SSIM) to the complex wavelet domain [17]. The principle behind the SSIM is that the Human Visual System (HVS) is mainly adapted to extract the structure of the objects from the images. Consequently, a measure of the structural similarity constitutes a good approximation of perceptual image quality. SSIM tries to ignore non structural distorsions, such as local intensity patterns and is defined in the spatial domain as: where C 1 and C 2 are two small positive constants;µ x and µ y are the mean of images x and y, respectively; σ x and σ y are the standard deviation of images x and y, correspondingly; and σ xy is simply: The maximum value SSIM can achieve is 1, and it is only achieved when x and y are identical. The original motivation to extend the SSIM to the complex wavelet domain was its high sensitivity to rotation, translation and scaling of the images, in the spatial domain [17]. Additionally, the fact that the HVS is organized in a multirresolution manner and that phase contains more structural information than magnitude in natural images, inspired this wavelet-based index. The CWSSIM is defined as: with c x,i and c y,i being the wavelet coefficients of x and y, respectively, at the decomposition sub-band i (N is the number of coefficients per sub-band), and c * denoting the complex conjugate. K is a small positive constant, whose task is to improve the robustness of the method in regions where the Signal-to-Noise Ratio (SNR) is poor [17].
CW SSIM(c x , c y ) can be re-written as [17]: The first component of this product is determined by the magnitude of the wavelet coefficients and, therefore, is equivalent to SSIM. The second component depends on the consistency of the phase changes between c x and c y and achieves its maximum value 1 when their difference in phase is a constant for all i. It is assumed that the structural information of local image features is principally located in the relative phase changes and a constant phase shift does not cause any structural distorsion [17]. Thus, the more similar x and y are, the higher the value of CWSSIM. The maximum value achievable is 1.
When having real data, it is not possible to calculate any of the metrics presented above, because the "ground truth" is unknown. Therefore, metrics such as the SNR can be employed to assess the denoising algorithms. For the experiment data, the SNR was calculated in the background of the image, as a reference. The SNR was computed as: µ represents the mean and σ the standard deviation. Due to the characteristics of the breast anatomy, it is hard to find constant intensity regions inside it and, therefore, not possible to calculate the object SNR. In this case, the visual appearance of the images and some intensity profiles were employed as the main indicator.

Simulations
The denoised simulated images corresponding to the highest noise level (worst case scenario) are shown in Figs. 2 and 3. The metrics for the whole simulated data set are reported in Table  1, Table 2, Table 3 and Table 4. VisuShrink 9.7e-4 6.4e-4 5.0e-4 4.1e-4 3.5e-4 BayesShrink 8.2e-4 5.9e-4 4.6e-4 3.8e-4 3.3e-4 NeighShrink 1.5e-4 8.7e-5 6.2e-5 4.9e-5 4.2e-5 Median filter 6.3e-3 6.3e-3 6.3e-3 6.3e-3 6.3e-3 SAM 9.8e-4 9.1e-4 8.8e-4 8.7e-4 8.5e-4 Wiener filter 1.9e-4 1.1e-4 7.6e-5 5.9e-5 4.9e-5 WND 1.2e-3 7.9e-4 6.0e-4 4.9e-4 4.2e-4  2(h) and 3(h)], although the MSE value associated to the latter is 1.24 times higher than that associated to the former. Even though, for more than 50 years, the MSE has been the dominant quantitative performance metric in the signal processing field, is simple to calculate and has a clear physical meaning, it has proven to fail to predict human perception of signal quality and fidelity, because it does not take into account any dependency between pixels, such as ordering, textures and patterns. Additionally, it completely ignores the signs of the error and the effects of these signs are significant from a perceptual point of view [18]. Thus, the optimal denoising method can not be selected solely based on the MSE and more reliable metrics should be used.   A safer metric would be the CWSSIM, because it focuses on the distortions that are meant to cause important perceptually alterations [17]. Correspondingly, for the DPC signals, our proposed method yields the highest CWSSIM [ Fig. 2(h)]. However, in the DCI case it indicates that VisuShrink is the best denoising algorithm for the four highest noise levels, although [ Fig.  3(b)] shows that this method causes oversmoothing, edge loss and ring artifacts. This result could be explained by the fact that the SNR of the DCI images is significantly lower than that of the DPC images, and therefore the CWSSIM works better for the latter kind of signals (see Table 5). In order to make the CWSSIM safer for DCI signals, an SNR threshold could be established, so we could make sure that above certain value the calculated CWSSIM would be trustable.
For the DCI signal, WND has the largest CWSSIM for the lowest noise level, and the second largest for the rest; this result is more coherent, because as [ Fig. 3(h)] proves, this technique causes no structural distortions at all and has a notable edge preservation ability.

Mammography data
A ROI surrounding a tumor was selected, to illustrate better the denoising effects (Figs. 4 and 5), and an intensity profile of one of the rows was selected to make more clear these effects (Figs. 6 and 7). The mammogram's background SNR values computed are reported in Table 5. As can be grasped from the tables, VisuShrink yields the highest SNR value for both signals. This is understandable, since these values were computed in constant intensity regions with no edges, so the oversmoothing caused by this technique boosts the SNR. Additionally, it is worth to point out that the WND method yielded the highest SNR among the three noise model-based methods in the DPC mammograms. Nevertheless, a disadvantage of the SNR is that it does not assess the preservation of edges and small features. Since these characteristics are important in mammography, this metric is not enough, although it can be a good starting evaluation point, e.g. to confirm that the noise level was indeed reduced. Consequently, to be able to perform a more reliable evaluation, we plotted the intensity profile of one of the rows traversing the tumor ( Figs. 6 and 7). These plots confirm that VisuShrink kills many important signal components and attenuates the edges. They also evidence that the median filter causes the largest reduction to the edge strength. Conversely, the WND method manages to denoise the signal while simultaneously preserving the borders. The latter effect is easily visualized, but cannot be quantified precisely. Regarding all the results obtained, it can be stated that the blind denoising techniques are able to provide an acceptable solution as long as they take into account statistical characteristics of the image and do not rely exclusively on the noise standard deviation. For instance, VisuShrink is the least optimal threshold of this kind, because it ignores these characteristics and strongly depends on the estimated variance value. On the contrary, BayesShrink and NeighShrink offer better solutions in general, because they consider this information.
Another important conclusion is that it is difficult to provide a reliable assessment of a denoising method, relying solely on the existent automatic performance metrics. As it was proved, there is still a long way to go over to develop a robust and exact metric.
A relevant and unique characteristic, and also advantage, of this denoising task is that the variance of the noise associated to each pixel can be calculated. Even so, it is fundamental to keep in mind that in image acquisition systems, there is always random noise that is not possible to calculate. In fact, due to this additional noise, the denoising process yielded much better outcomes in the simulated scenario. Therefore, albeit the noise a-priori knowledge driven denoising methods are more trustable, they are still not perfect. A good way to reduce the uncertainty of these methods is to improve the acquisition hardware, e.g. gratings and detectors.

Conclusion
A wavelet based noise model driven denoising algorithm (WND) for Differential phase contrast mammography data was introduced. The most innovative element of WND is the translation of the noise variance from the spatial to the wavelet domain. Although there have been several related approaches [14], none of them dealt with the fact of having a noise variance different for each image pixel that could be calculated.
Several other blind as well as model-based existing denoising methods were also imple-mented. All these denoising techniques were tested with simulated and real mammography DPC and DCI data. Two performance metrics, the Mean-Squared error (MSE) and the Complex Wavelet Structural Similarity Index (CWSSIM), were used for evaluation in the simulated scenario. Since these metrics demand knowing the ground truth, they can not be employed when dealing with real data. Therefore, the SNR and was used to assess the real data denoising results.
The results yielded by WND were satisfactory. Its outcomes demonstrated a reduction of the noise while showing preservation of the edges in both the simulated and real scenarios. The CWSSIM and SNR metrics also gave a good grade to this method. However, as was largely discussed, these metrics are not equivalent to the human visual system and, therefore, it would be ideal to carry out a reader study with experts, who can provide a more reliable assessment of the proposed method. Another important issue is the noise calculation. The noise model for the three signals has been studied and the noise variance associated to each pixel can be estimated. However, this model only considers the detector and phase-jitter noise and, in this work, the latter was not included because, according to the manual of the motor used to move the grating G2, it was insignificant [16]. Even so, it would be worth to get to measure this uncertainty and try to include other possible sources of noise, disregarded by this model, in the noise calculation, so it can be more reliable and the denoising can yield better outcomes.
As stated in the discussion section, the blind denoising methods that take into account some image information can be acceptable (e.g. Bayes Shrink and NeighShrink). However, since the noise model is known, it is better to have a denoising stage based on it and not on general assumptions. The widely used median filter evidenced more disadvantages than advantages, so it is not recommended to employ it in mammoDPC.
As to the beam hardening effect, it depends on both the sample composition and thickness, and the x-ray energy or tube voltage. However, since the grating-based x-ray interferometry system is designed for a fixed tube voltage, in our case only the sample properties have an effect on beam hardening. To image different objects, in principle, f s 1 should be measured and calibrated with different materials and sample thicknesses in advance. Nonetheless, our major interest here is mammography, and for imaging breast samples, the beam hardening effect is not pronounced, as presented in [20] and as we can grasp from the similarity between f r 1 and f s 1 . Breast samples are usually assumed to have a simple composition (fat and glandular tissue, which have close attenuation properties), therefore f s 1 is estimated with only one material in our case. For situations with severe beam hardening and complex sample compositions with high atomic number materials, we suggest to do the calibration for f s 1 as described above, in order to get a good estimation of the noise. Since the noise models we are using take into account the difference in f 1 between the reference and sample measurements, we believe that our denoising approach would still be useful in these situations.
In conclusion, a working noise model based denoising algorithm for Differential phase contrast mammography was designed, implemented and evaluated, yielding good results and showing potential in this emerging field.