key: cord-0044816-hde9waz3
authors: Kornak, John; Boylan, Ross; Young, Karl; Wolf, Amy; Cobigo, Yann; Rosen, Howard
title: Bayesian Image Analysis in Fourier Space Using Data-Driven Priors (DD-BIFS)
date: 2020-05-16
journal: Information Processing and Management of Uncertainty in Knowledge-Based Systems
DOI: 10.1007/978-3-030-50153-2_29
sha: 6fe96eb60bdab5db7c00f0385fb3341c6ed138d9
doc_id: 44816
cord_uid: hde9waz3

Statistical image analysis is an extensive field that includes problems such as noise-reduction, de-blurring, feature enhancement, and object detection/identification, to name a few. Bayesian image analysis can improve image quality, by balancing a priori expectations of image characteristics, with a model for the noise process via Bayes Theorem. We have previously given a reformulation of the conventional Bayesian image analysis paradigm in Fourier space, i.e. the prior distribution (the prior) and likelihood are given in terms of spatial frequency signals. By specifying the Bayesian model in Fourier space, spatially correlated priors, that are relatively difficult to model and compute in conventional image space, can be efficiently modeled as a set of independent processes across Fourier space. The originally inter-correlated and high-dimensional problem in image space is thereby broken down into a series of (trivially parallelizable) independent one-dimensional problems. In this paper we adapt this Fourier space process into a data-driven framework in which the Fourier space priors are built empirically from a database of images and then used to enhance future images. We will describe the data-driven Bayesian image analysis in Fourier space (DD-BIFS) modeling approach, illustrate it’s computational efficiency and speed. Finally, we give specific applications of DD-BIFS to improve the quality of arterial-spin-labeling (ASL) perfusion images via a database of human brain positron emission tomography (PET) images.

Bayesian image analysis models provide a solution for improving image quality in image reconstruction/enhancement problems by incorporating a priori expectations of image characteristics along with a model for image noise. We have previously presented an approach to reformulating the Bayesian Image Supported by NIH R01 EB022055 and U19 AG063911. analysis in Fourier Space (BIFS) [15] . Spatially correlated prior distributions (priors) that are difficult to model and compute in conventional image space can be modeled as independent across locations in Fourier space. The original high-dimensional problem in image space is thereby broken down into a series of one-dimensional problems, leading to easier specification and implementation, and faster computation.

Consider x to be a true (or idealized) image that we wish to recover from a suboptimal image dataset y. (We use the common shorthand notation of not explicitly distinguishing the random variables and the corresponding image realizations [5] , i.e., we use lower case x and y throughout.) The Bayesian image analysis paradigm incorporates a priori desired spatial characteristics (e.g., smoothness) via a prior distribution for the true image x and defines the noise degradation process via the likelihood.

Instead of the conventional Bayesian image analysis approach of generating prior and likelihood models for the true image x based on image data y directly, we formulate them via their discrete Fourier transforms representations: Fx and Fy. After applying Bayes' Theorem, the posterior becomes

The key component of the BIFS formulation that leads to its useful properties of easy definition and computational speed, is that we define both the prior and likelihood (and therefore the posterior) to consist of a set of independent distributions across Fourier space locations. Desired spatial correlation in image space is induced by allowing the parameters of the distributions to change over Fourier space [18, 20] . This independence definition can be contrasted with conventional Bayesian image analysis using Markov random field (MRF) priors, where Markovian neighborhood structures are used to induce correlation structures across pixels via joint or conditional distributional specifications [4, 5, 9] .

When defining a spatially correlated prior in image space via a set of independent processes across Fourier space, the full conditional posterior at a Fourier space location k = (k x , k y ), or for volumetric data (k x , k y , k z ), now only depends on the prior at that location k, i.e.,

where we use Fx k as shorthand for (Fx) k . The joint posterior density for the image is then

where K is the set of all Fourier space point locations.

The standard process of generating the BIFS prior distribution given in Kornak [15] is based on choosing a pair of distributions to be applied as priors at each location in Fourier space (one for the modulus and the other for the argument of the complex value signal) and a set of parameter functions to define how the parameters of the distributions vary over Fourier space. In contrast, for the datadriven approach, although we again choose the probability distribution forms across Fourier space, the parameters are estimated empirically from a database of transformed images. That is, all of the images in the database are Fourier transformed, the data at each location in Fourier space are extracted, and the distribution parameters for that Fourier space location are estimated from that data. In this way empirical maps of parameter estimates is generated over space. These parameter estimates are then used as the parameters for the prior at each Fourier space location. Separate priors and associated parameter functions are defined for each of the modulus and argument of the complex value at each Fourier space location. Working with the modulus and argument provides a more natural framework for defining prior information at specific Fourier space locations (i.e., specific spatial frequencies) than working with the real and imaginary components, because prior information (e.g., smoothness, edges, or features of interest) most strongly relates to the magnitude of the process involved.

As for the prior, the BIF likelihood is again modeled separately for the Modulus and Argument of the signal at each Fourier space location. Because we model based on independence across Fourier space points, a range of different noise structures (defined in Fourier space) can readily be incorporated into the likelihood π(Fy k |Fx k ). For example, the combination of independent and identically distributed (i.i.d. ) zero mean Rayleigh noise/Rician likelihood [19] for the modulus with uniform argument on the circle in Fourier space corresponds to the likelihood model of i.i.d. Gaussian noise in image space.

It is at the posterior estimation stage that the computational gains of the independent BIFS formulation are ultimately realized. Posterior estimation in conventional Bayesian image analysis tends to focus on maximum a posteriori (MAP) estimation (minimizing a 0-1 loss function), because it is the most computationally tractable. In the BIFS formulation the MAP estimate can be estimated with added efficiency by maximizing the posterior distribution separately at each Fourier space location and then taking the inverse Fourier transform of the Fourier space MAP estimates, i.e,

This computational efficiency contrasts with conventional Bayesian image analysis, where even the most computationally convenient MAP estimate requires iterative computation methods such as conjugate gradients or expectationmaximization to obtain it.

Implementation of DD-BIFS modeling requires the following steps:

1. Fast Fourier transform (FFT) all images in the database that are to be used to build the DD-BIFS prior, i.e., from image space into Fourier space. 2. Choose the distributional form of the prior at each location in Fourier space 3. Estimate the parameters of the prior at each location in Fourier space using the data from the corresponding Fourier space locations across the subjects in the database. 4. Define the likelihood in Fourier space. 5. FFT the dataset to be reconstructed from image space into Fourier space. 6. Combine the DD-BIFS prior and likelihood for the image at each Fourier space location via Bayes' Theorem to generate the DD-BIFS posterior 7. Generate Fourier space MAP estimate by maximizing the posterior at each Fourier space location. 8. Inverse FFT the Fourier space MAP estimate back to image space and display

In this example we use simulations to drive the formulation of a DD-BIFS prior which is subsequently applied to independent data. We simulated 1000 256 × 256 images representing lesions/tumor patterns. The number of lesions was modeled as a Poisson process and the lesions were simulated as randomly positioned truncated Gaussian probability density functions (resembling bumps) with random intensity, and standard deviation on each axis, and with correlation distributed uniformly between −1 and 0, so that the process was anisotropic (i.e. so the spatial autocorrelation was not uniform in all directions).

We use the following model applied at each Fourier space location:

-Gaussian (normal) prior for the modulus:

where k is the complex noise at Fourier space location k. (Note this model is not Gaussian noise in image space, for which we use a Rayleigh noise model/Rician likelihood, for the modulus.)

The prior for the modulus at each Fourier space location was then generated from the empirical estimates of the mean (μ k ) and standard deviation (τ k ) at the corresponding location across the simulated datasets. The global posterior mode was then be obtained by generating the posterior mode at each Fourier space location based on conjugate Bayes for the Gaussian distribution [8] with

and Arg(x k,MAP ) = Arg(x k,MAP ) (when the prior for the argument is uninformative). Figure 1 shows a single new realization (i.e. that was not included in the simulation set) of the process at the top-left, the same image with added noise at the top-right, a parameter function-based BIFS MAP reconstruction on the bottom-left (approximating a pairwise absolute difference MRF prior Bayesian reconstruction); and the DD-BIFS prior on the bottom-right. The parameter function-based BIFS prior does a reasonable job of enhancing lesions, in particular the blurred lesion furthest to the right. However, the DD-BIFS reconstruction clearly improves the enhancement of the simulated lesions beyond that of the parameter function-based BIFS. The DD-BIFS reconstruction is able to better retain the detail of the lesions, in particular, their non-isotropic elongated form.

Truth + noise BIFS problem specific prior BIFS simple (~abs.diff. MRF) It should be noted here that primary objective of the study lies with visualization and not direct quantification of image intensities. If the clinical goal is to detect tumors then the objective of reconstruction is to have visually clear tumors that the radiologist can readily identify. Alternatively, if the goal is differentiation of benign vs cancer tumors, the decision may be based on shape characteristics so that accentuating them may help the clinician to visually differentiate between the lesion types.

The overall objective of this study is to process arterial spin labeling (ASL) perfusion MRI to enhance blood flow patterns in the brain associated with frontotemporal lobar dementia (FTLD). ASL perfusion has been shown to be sensitive to FTLD pathology, showing hypoperfusion in frontal regions, potentially providing a cheaper and radiation free alternative to the conventional (FDG)-PET (fluorodeoxyglucose positron emission tomography) [7, 12] . Preliminary data indicates that ASL based perfusion patterns associated with FTLD are coherent with PET, albeit with reduced image quality due to increased noise (Fig. 2) . The objective of our BIFS modeling procedure is to improve ASL perfusion maps toward the quality of PET by using a DD-BIFS prior generated from a database of PET images. 

PET and ASL perfusion pairs for individuals acquired around the same time are taken from the Frontotemporal Lobar Degeneration Neuroimaging Initiative (FTLDNI), a multi-center biomarker trial aimed at identifying promising markers as surrogate endpoints for FTLD disease progression in clinical trials.

Volumetric MPRAGE sequences at UCSF were used to acquire T1-weighted structural images of the entire brain (Sagital slice orientation; slice thickness = 1.0 mm; slices per slab = 160; in-plane resolution = 1.0 × 1.0 mm; matrix = 240 × 256; TR = 2,300 ms; TE = 2.98 ms; TI = 900 ms; flip angle = 9 • ). Pulsed ASL (PASL) imaging was acquired using QUIPSSII with Thin-slice TI1 Periodic Saturation (Q2TIPS) sequence incorporated in a PICORE (Proximal Inversion with Control of Off-Resonance Effects) labeling scheme [16] . The periodic saturation pulses started at the postlabeling delay inversion time TI1 = 700 ms after the in-plan presaturation radio infrequency pulse; the readout started at the postlabeling delay inversion time TI2 = 1800 ms. The repetition and echo time were TR/TE = 2500/11 ms. Sixteen slices were acquired, each 6 mm thick with a 7.2 mm center to center distance and a matrix 64 × 56 of 4 × 4 mm 2 in-plane voxel resolution.

PET Data: PET data were acquired at the Lawrence Berkeley National Laboratory on a Siemens ECAT EXACT HR scanner or a Siemens Biograph PET/CT scanner. FDG was supplied by a radiopharmacy (IBA Molecular). Six emission frames lasting 5 min each were acquired starting 30 min post-injection. Attenuation correction was performed using a 10 min transmission scan on the ECAT scanner or a low-dose CT scan on the Biograph, both being acquired prior to PET acquisition. For both scanners, PET data were reconstructed using an ordered subset expectation maximization algorithm with weighted attenuation and images were smoothed with a 4 mm Gaussian kernel with scatter correction. Final resolution was calculated using Hoffman phantom: 7 × 7 × 7.5 mm for ECAT and 6.5 × 6.5 × 7.25 mm for Biograph.

Pre-processing: Before any prepossessing of the images, all T1-weighted images were visually inspected for quality control. Images with excessive motion or image artifact were excluded. T1-weighted images underwent bias field correction using N3 algorithm, the segmentation was performed using SPM12 (Wellcome Trust Center for Neuroimaging, London, UK, http://www.fil.ion.ucl.ac. uk/spm) unified segmentation [10] . A customized group template was generated from the segmented gray and white matter tissues and cerebrospinal fluid by non-linear registration template generation using Large Deformation Diffeomorphic Metric Mapping framework [2] . Native subjects space gray and white matter were geometrically normalized to the group template, modulated and then smoothed in the group template. The applied smoothing used a Gaussian kernel with 8 mm full width half maximum. All steps of the transformation were carefully inspected from the native space to the group template. Linear and non-linear transformations between the group template space and International Consortium of Brain Mapping (ICBM) [17] was applied.

Frames of the ASL acquisition were corrected for motion, co-registered with the first frame (M0) using FSL [13] . An automatic quality control process removed tagged/untagged pairs of frames when the relative root mean square (RMS) distance value between two consecutive frames was higher than 0.5 mm.

The subject was dropped if this RMS value was higher than 1 mm. Differential perfusion images were created by subtracting unlabeled from adjacent labeled frames and averaging these subtraction images [1] . Susceptibility artifacts along the phase-encoding direction were corrected for the M0 frame and perfusion map using ANTs [3] restricted to the coronal axis of the patient. Cerebral Blood Flow (CBF) was calculated by applying the Buxton kinetic model to the perfusion map [6] . CBF data was processed to obtain partial volume corrected maps of gray matter perfusion, based on the tissue segmentation using transformation matrix from T1 to M0 as previously described [7, 11, 14] . Normalized CBF values were obtained by dividing the voxel CBF value by the mean Calacarin CBF value region of interest. Calcarine was selected based on the observation that FTD variants do not impact this area neither the acquisition field of view.

Analyses on partial volume corrected, non-partial volume corrected and normalized perfusion images was performed in MNI space using the structural set of geometric transformations and smoothed with an isotropic 12 mm Gaussian kernel full width half maximum. All CBF images were visually inspected in the native and MNI spaces. Poor quality images out of the field of view, with large susceptibility or motion artifacts were removed from the study.

Finally, voxels in the ASL image were rescaled to match the dynamic range of the PET images. To do this the voxels in the ASL were first ordered by intensity, ignoring their spatial coordinates. Then a 10% random sample was drawn from all the voxels in all the PET images of the database, and they were ordered by intensity. The n'th brightest voxel in the ASL to be processed was reassigned the intensity of the m'th brightest voxel in the PET sample, where m = round n * np na , and where n p and n a are the sizes of the PET subsample and the ASL image respectively, and n ranges from 0 to n a − 1.

The same DD-BIFS model set as that used for the lesion simulation study. 101 subjects were used to build the PET prior and 1 individual with corresponding ASL-perfusion scan was reserved for evaluation. Figure 3 shows the results of the reconstruction of an individual's ASL image. In the top-left panel is the PET image that we would like to emulate and in the top-right the corresponding original ASL perfusion MRI for the same individual. In the bottom left the DD-BIFS prior reconstruction is displayed (with mask applied) and at the bottom right is DD-BIFS prior reconstruction with additional shrinkage on the prior variance (specifically the standard deviation is multiplied by a factor of 0.01). The reconstructed image does adjust the ASL-perfusion MRI so that it emulates the characteristics of the PET, and shrinking the prior variance (essentially putting more weight on the prior) serves to further move the original ASL to an image that emulates what the clinician might see with PET.

Ongoing work on this project will focus on performing full clinical validation of DD-BIFS to assess clinical applicability as a surrogate for PET. 

The DD-BIFS modeling framework provides a powerful new approach to using information available in large databases to improve reconstruction in individual images. In particular, the independence across Fourier space specification allows for fast and efficient computation which can be further improved with parallelization. Additionally, the ability to efficiently use empirical prior information from a database of images without the need for explicit modeling provides a powerful approach to improving image quality. Our preliminary application to the reconstruction of ASL perfusion MRI shows great promise and large-scale validation work is currently under way to determine its applicability in clinical practice.

Experimental design and the relative sensitivity of BOLD and perfusion fMRI

Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation

Advanced normalization tools (ANTS)

Spatial interaction and the statistical analysis of lattice systems (with discussion)

Digital image processing: towards Bayesian image analysis

Dynamics of blood flow and oxygenation changes during brain activation: the balloon model

Hypoperfusion in frontotemporal dementia and Alzheimer disease by arterial spin labeling MRI

Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images

The neural correlates of emotion regulation by implementation intentions

A non-parametric approach for co-analysis of multi-modal brain imaging data: application to Alzheimer's disease

Distinct cerebral perfusion patterns in FTLD and AD

Pattern of cerebral hypoperfusion in Alzheimer disease and mild cognitive impairment measured with arterial spin-labeling mr imaging: initial experience

Bayesian image analysis in fourier space (BIFS)

QUIPSS II with thin-slice TI 1 periodic saturation: a method for improving accuracy of quantitative perfusion imaging using pulsed arterial spin labeling

A probabilistic approach for mapping the human brain: the international consortium for brain mapping (ICBM)

Central limit theorem for stationary linear processes

Mathematical analysis of random noise

Exploring an ozone spatial time series in the frequency domain

Renaud La Joie and Amelia Strom for providing information on the PET dataset acquisitions.