extfieldsofstudy: Medicine Physics

text:
Analysis of 25% duty-cycle up-conversion passive mixer for UWB transmitter

The performance of differential up-conversion passive mixer operating in the 3–5 GHz band for UWB transmitter driven by 25% duty-cycle clocks is studied and analyzed. A reasonably accurate LTI equivalent circuit accounting for harmonics is derived. We have demonstrated that the conversion gain, input and output impedances of the proposed LTI equivalent circuit matches those of the LTV up-conversion passive mixer. The LTI model can be used to define the reachable design space. The mixer input and output impedance can be tuned by varying resistor at the mixer output allowing for impedance match to connected circuits. We have shown how each design parameters impacts critical performance of the up-conversion passive mixer. Simulations of the proposed up-conversion passive mixer show −3.68 dB of conversion gain, 15.9 dBm of OIP3, 116.7dBm of OIP2 and a NF as low as 5.2 dB while consuming 1.5 pW. Furthermore linearity and ports isolation performances due to voltage threshold mismatch variation of ±60 mV shows an OIP3 greater than 15.85dBm and an LO_RF isolation smaller than −51.9 dB.

Introduction
The standard Wireless Body Area Networks (WBAN) is known to provide efficient, low power and optimized wireless link between implantable and/or wearable devices for continuous and real time physiological signal monitoring [1,2]. In particular WBAN wearable devices can be used to wirelessly monitor Coronavirus patients without exposing medical staff to the virus, therefore limiting the spread of the virus within the hospital.
The Ultra-wideband (UWB) offers robust performance for WBAN with low complexity transmitter and ultra low power operation. The UWB power spectral density emission limit of -41.3 dBm/MHz provide safe power levels for the human body [3][4][5]. The wireless sensor node should be able to operate until ten years under battery energy supply. Passive mixer-based UWB transmitter allows high performance and low power wearable devices implementation, therefore increasing battery life. Unlike passive mixers for receivers, there are only a few works that deal with passive mixers for UWB transmitters. Active mixers have the advantage of having high conversion gain and good baseband-RF ports isolation. Various circuit techniques have been proposed to increase the conversion gain and improve the linearity of active mixers over a wide bandwidth [6,7].
On the other hand, passive mixers have the following performances which do not exist in active mixers [8][9][10][11]: extremely low power consumption because passive mixer dissipates no DC power, robust linearity, good noise figure performance, the passive mixer-based transmitter occupies a very small chip area, better LO-RF and LO-baseband ports isolation performance [12].
The 50% duty cycle CMOS up-conversion mixer in [13] is based on a chain of switched CMOS inverters for 2.45 GHz and 400 MHz bands. The circuit has similarities with the H-Bridge ring mixer. The up-conversion mixer achieves a conversion gain of À5dB and consumes 23.1mW from 3.3 V supply voltage. In [14] authors studied and analyzed conversion gain and noise performances of 25% duty-cycle passive mixer for IQ direct-conversion transmitters. They made Linear Time Variant (LTV) as well as Linear Time Invariant (LTI) analysis to derive the input/output transfer functions at LO and its harmonics. The transmitter proposed in [15] employs a pseudo double-balanced up-conversion passive mixer with more than 22 dB LO leakage rejection. The power consumption is 14.4mW under 1.2 V supply.
The double-balanced IQ passive mixers with 50% duty cycle are widely used in receivers as down converter [9]. To be able to use this structure in transmitter as up-converter it is necessary to https://doi.org/10.1016/j.aeue.2020.153295 1434-8411/Ó 2020 Elsevier GmbH. All rights reserved. use a duty cycle less than or equal to 25%. Indeed for a duty cycle of 50% the two channels I and Q on the RF sides cannot be shorted together because at any given instant there is a switch from I channel and another from Q channel at the ON state. On the other hand, for a duty cycle less than or equal to 25% there is a single switch in the ON state which eliminates any possible interaction between I and Q channels.
In this paper we present an analysis of a 25% duty cycle upconversion passive mixer intended to be used in UWB transmitter for WBAN applications.
Section 2 presents transmitter block diagram, modeling details of the 25% duty-cycle up-conversion passive mixer and derives the LTI equivalent circuit. A comparative study between LTI equivalent circuit accounting for harmonics and the LTV up-conversion passive mixer is carried out in Section 3. Implication of LTI equivalent circuit for design of passive mixer-based transmitter is discussed in Section 4. Section 5 concludes the paper.

Modeling of the 25% duty-cycle up-conversion passive mixer
The proposed transmitter block diagram front-end is shown in Fig. 1. The differential baseband (BB) analog signal from the DAC is first filtered with a Low Pass Filter (LPF) to remove image components and quantization noise. The filtered baseband signal are up-converted to UWB RF using a differential double-balanced upconversion passive mixer. The differential output UWB RF signal is buffered with a high input impedance buffer then converted to single ended signal. The output signal of the driver is typically amplified by a power amplifier (PA). Differential double balanced up-conversion passive mixer can eliminate BB and LO feed through. Voltage mode passive mixer operation was adopted because designing of a voltage buffer is simpler than the current buffer at UWB RF frequency. Current-mode passive mixer is typically used in receivers [14]. Fig. 2 depicts the proposed differential double balanced upconversion passive mixer with 25% LO duty cycle. Transistor switches are biased at V p via the load resistors R RF . C RF is the input impedance seen from the buffer. The series capacitors (C bb ) isolate the DC bias of the LPF and the mixer. R lpf is the LPF output resistance. In the following sections we will study and analyze in time domain the proposed time-variant up-conversion mixer circuit.

RF voltages
To analyze the differential double balanced up-conversion passive mixer we begin by presenting its basic operation in the time domain. We assume that the switches are ideal, except that they have non zero ON resistance R on . Fig. 3 shows the proposed half circuit used to analyze the up-conversion passive mixer. C par models the input parasitic capacitor of the NMOS switch. The parasitic capacitors of the switches can be lumped into the baseband impedance. These capacitors are ignored in the following analysis.
The differential baseband input voltage from LPF, V bb (t) is given by: us define the following periodic functions corresponding to the four LOs.
The output voltage V R1 t ð Þ shown in Fig. 3 is given by: In the ON phases, the charge flow Qc into RF capacitor C RF during T LO 4 can be written as: where In the OFF phases the capacitor C RF discharges through R RF and the voltage V R1 t ð Þ is given by: and the discharge of C RF through R RF correspond to equation (8).
Conservation of charge implies Q ch ¼ Q dis from which V M can be computed by:

Down-converted voltage V d t ð Þ for the up-conversion passive mixer
Up-conversion passives mixers simultaneously down-converts the RF signals to baseband and up-converts those to RF (transparency property) as reported in [11]. Since the LO phases are not overlapping in time, the down-converted voltage V d t ð Þ can be written as (10).
LO signals can be represented by their Fourier Series.
Then we can derive V d t ð Þ (16). For n ¼ Ç1, the voltage V d t ð Þ is given by (17) Replacing V M by its expression (9), the component of the downconverted voltage at the baseband frequency f bb is obtained (18).

LTI equivalent circuit of the 25% duty cycle up-conversion passive mixer
At baseband frequency the total input current of the upconversion passive mixer I bb;fbb (Fig. 2) is given by (20). According to (19) and (20) the input current relation can be determined as shown in the equation (21).
Using (21) the expression between V bb jx bb ð Þ and I bb ðjx bb Þ can be derived as follows: This means that R RF is in series with Z bb jx bb ð Þ. Then we can establish in Fig. 4 the linear time invariant (LTI) equivalent circuit of the differential double balanced up-conversion passive mixer shown in Fig. 2. In this model the weighting factor q accounts for the conversion gain and Z s represents an additional power loss due to harmonics down-conversion from RF side to baseband side. From the equivalent circuit of Fig. 4 the voltage V bb jx bb ð Þ is given by (23) where The factor q and the impedance Z s can be computed as shown in (25).
To get maximum power transfer between the LPF and the upconversion mixer, the input impedance of the mixer should be matched to the output impedance of the LPF. From the equivalent circuit shown in Fig. 4, the input impedance can be determined: We can notice that in equation (26) the input impedance Z in can be tuned by changing R RF . The range is limited as shown in (27).
This analysis will be confirmed by simulation in the next section.

Switching ON resistance optimization
In order to set the optimal size of NMOS switches that set the switching ON resistance R on , we made S-parameter simulation of the circuit of Fig. 4 using ADS (Advanced Design System) tool. Fig. 5(a) shows S22(dB) for swept R on . From this simulation results we can observe that R on should be set greater than 34 O to have a good match for RF frequency between 4 and 5 GHz. To optimize the value of R on we run Harmonic Balance simulation on the upconversion passive mixer to plot output second-order intercept point OIP2 versus V GS . A two-tone generator is applied to the baseband input. One tone at 500.5 MHz and the other at 499.5 MHz. f LO = 4 GHz and the LO duty cycle = 25%. We used 1 ps for the rise and fall time of LO signals. For TSMC 0.18um technology, the dimensions of the NMOS are W = 42.2 mm and L = 0.18 mm. The plot is shown in Fig. 5(b). We can see that the optimum point OIP2 = 116.6dBm is obtained for V GS = 0.76 V, this corresponds to R on = 38.3 O. Fig. 6(a) shows time domain simulation results of the proposed up-conversion passive mixer with f LO set to 4 GHz, duty-cycle equals to 25% and f bb set to 500Mhz. V bb , V d are signals at baseband side and V RF is the up-converted signal at RF side. The mixer output power is set to À41.4dBm for f RF equals to 4.5 GHz. The up-conversion passive mixer output spectrum (VRF_spectrum) and down converted signal spectrum (Vd_spectrum) are given in Fig. 6(b).
The . So the 25% duty-cycle differential double balanced passive mixer did the up-conversion correctly. Each transistor drain is biased at V p = 300 mV and without oscillator circuit the up-conversion passive mixer consumes only 5pA.

Input impedance and output impedance versus R RF , C RF and BB_freq
Now we will compare input and output impedances of the LTI equivalent circuit shown in Fig. 4 and the LTV up-conversion passive mixer at LO duty cycle of 25% shown in Fig. 2. For the LTI equivalent circuit, input and output impedances are denoted by Zinequ and Zoutequ respectively. For the LTV up-conversion passive mixer, input and output impedances are denoted by Zinupmix and Zoutupmix respectively. Fig. 7(a) shows simulation results comparing the magnitude of input and output impedances for swept R RF from 1 O to 40 KO. Fig. 7(b) presents input and output impedances for both circuit versus BB frequency.
We can observe that input and output impedances of the LTI equivalent circuit matches the input and output impedances of the LTV up-conversion passive mixer. This simulation proves that the LTI equivalent circuit can be used to analyze the 25% duty cycle differential up-conversion passive mixer for UWB transmitter. We will study in Section 3 the LTI equivalent circuit accounting for harmonics.

The LTI equivalent circuit accounting for harmonics
In the previous analysis, we defined the down-converted voltage V d only at baseband frequency f bb . However, Fig. 6(b) indicates that V d contains harmonics at (2f LO ± f bb ), (2f LO ± 3f bb ), 3f bb ,. . .as well as its fundamental. Therefore, the input current I bb can be computed: The first term at f bb represents the current coming from the LPF and the other terms represent the return current due to V d . Eq. (28) show that the current I bb decreases compared to the fundamental current used in the previous analysis. Thus, the input impedance Zinupmix will increases. In order to account for the dissipation of power due to these harmonics we multiply the impedances by a factor h as shown in Fig. 8. This model presents a reasonably accurate LTI equivalent circuit of the up-conversion passive mixer accounting for harmonics.
A comparison between conversion gain of the 25% upconversion mixer and its LTI equivalent circuit implies that h = 2 and q 1 ¼ 0:9cosu where cosu-0.

Conversion gain, input and output impedances versus R RF and BB_freq
In the context of the above LTI equivalent model, Fig. 9 shows simulation results comparing the conversion gain of the LTI  equivalent circuit and the LTV up-conversion passive mixer (denoted by C_Gainequ and C_Gainupmix repectively) versus BB_freq and for swept R RF accounting for harmonics. We can see that C_Gainequ follows closely C_Gainupmix. Fig. 10  Thus, in order to match the LPF output impedance as well as the buffer input impedance, the mixer input/output impedance can be tuned by varying R RF .
3.2. Noise factor, S11 and S22 for swept R on accounting for harmonics In order to reduce noise we need to reduce R on and increase R RF as shown in Fig. 11(a) and (b). But this will increasing the size of the transistors and thus will increase the parasitic input capacitor of the NMOS switches. In addition it is recommended to have an impedance matching at the input and output of the upconversion passive mixer. For this we have made S parameters simulation for several values of R on as shown in Fig. 11(c) for S11 and Fig. 11(d) for S22. R RF was set to 3KO for a noise factor of 3.36. From Fig. 11(c, d), we deduce that R on must be less than 46 O to get impedances matching. In order to reduce parasitic input capacitor of the NMOS switches we choose the same dimensions of the NMOS (W = 42.2 mm and L = 0.18 mm) used in the previous simulation taking into account fundamental components only. We can see also whether these dimensions hold or not if we account for harmonics.

Performances of the 25% duty-cycle up-conversion passive mixer
Now we will evaluate the performance of the up-conversion passive mixer using the dimensions of the transistors defined in the previous section. To do this, we start by evaluating OIP2 for different widths of the transistors. In the Fig. 12(a) we can see that the maximum of OIP2 corresponds to W = 41 mm. This value is very close to the value found in the previous section (W = 42.2 mm). This confirms that the LTI equivalent circuit model accounting for harmonics in Fig. 8 can be used to find the up-conversion passive mixer design parameters. Fig. 12(b) shows the plot of mixer output power (Pout) versus baseband input power (Pbb). The mixer inputreferred 1 dB compression point IP 1-dB = 8.818 dBm and outputreferred 1 dB compression point OP 1-dB = 3.912 dBm. To choose the optimized value of the capacitor C RF , the simulation of the conversion gain as a function of C RF was carried out and presented in Fig. 12(c). It is clear that the maximum gain corresponds to C RF equal to 0.18 pF. For this value the up-conversion passive mixer achieves a conversion gain greater than À3.68 dB for any input power lower than 8 dBm.
When we studied the up-conversion passive mixer circuit we assumed that all switches exhibit identical characteristics and bias. In reality NMOS transistors used in this up-conversion mixer are not perfectly matched or balanced. They exhibit voltage threshold (V th ) and W/L ratio mismatch caused by wire bonding variation and process tolerance. V th depends on the doping levels in the channel and the gate, and these levels vary randomly from one transistor to another [16]. V th mismatch lead to DC offset at the circuit nodes. First this may drive the stages of the circuit into nonlinear operation or saturation. Second this will produce some LO to RF, LO to BB and BB to RF feedthrough. Since this up-conversion mixer is based on identical NMOS transistors connected to achieve a double balanced differential structure, it is important to evaluate its inter modulation distortion performance as well as its isolation between ports in case of voltage threshold mismatch. Fig. 13 presents the simulation of OIP2 and output third-order intercept point OIP3 as a function of the variation of Vth (Vth_mismatch) of one of the mixer switches. We can see that OIP2 degrades quickly from 116.7dBm to 45.92dBm (for Vth_mismatch = 0.2 V) and to 43dBm (for Vth_mismatch = À0.2 V). On the other hand OIP3 decreases slowly from 15.9dBm to 15.44dBm (for Vth_mis match = 0.2 V) and to 14dBm (for Vth_mismatch = À0.2 V). Fig. 14 (a) presents the simulation results of leakage between mixer ports without mismatch for an UWB output power of À41.4dBm. The up-conversion passive mixer achieves a BB_RF isolation = À5 6.4 dB, LO_RF isolation = À106 dB and LO_BB isolation = À154.9 d B. With Vth_mismatch = 60 mV Fig. 14(b) shows the isolation between ports as follow: BB_RF isolation = À45.dB, LO_RF  isolation = À51.9 dB and LO_BB isolation = À136.1 dB. The presence of LO component is due to bias mismatch between NMOS gates in the mixer. Therefore the OIP2 and LO_RF isolation of this upconversion passive mixer are sensitive to Vth variation. But for Vth_mismatch variation of ± 60 mV the OIP3 is almost constant around 15.9dBm, OIP2 is higher than 57dBm and LO_RF isolation is less than À51.9 dB. Those values are good enough for linearity and ports isolation performances of the differential up-conversion passive mixer with duty cycle of 25%. However, LO signals should be designed low phase-noise and a careful balanced layout of the up-conversion mixer must be done to maintain balanced operation.

Implication of LTI equivalent circuit for design of passive mixer-based transmitter
The above analysis aims to offer an equivalent LTI circuit as close as possible to a real differential up-conversion passive mixer with LO duty cycle = 25% for UWB transmitter. The accuracy verification of the proposed LTI model was done using the simulation results of input and output impedances as well as the conversion gain shown in Figs. 7, 9 and 10. The resulting curves matches those found with LTV up-conversion passive mixer which clearly shows that the LTI model is equivalent to the actual circuit. The model proposed in Fig. 8 takes into account LO and BB harmonics, mixer input AC coupling capacitors Cbb, filter output impedance and buffer input impedance. This model is the result of analysis based on time domain equations to capture transient effects. In the simulation we used a nonzero LO rise and fall times. The maximum conversion gain is À3.68 dB. Compared to the study carried out in [14], the authors used a purely capacitive output load, they neglected the effect of input AC coupling capacitors and they assumed the LO clocks are ideal with 25% duty-cycles and zero rise and fall times.
They find a maximum conversion gain of À13 dB. One implication of having this LTI equivalent model is that designers will  quickly define the reachable design space from the constraints on the performance of the mixer. Thus allowing to optimize components sizes of the up-conversion passive mixer. Using the simple LTI equivalent circuit, we can see how each design parameters impacts critical performance of the up-conversion passive mixer like input/output matching, linearity, noise and power consumption.
Another implication is that both LTI equivalent circuit of the 25% duty cycle up-conversion passive mixer in Figs. 4 and 8 can be used for different values of BB_freq, R RF , C RF , R lpf and R on . For medical applications using devices worn by patients and powered by battery, the transmitter must be low power. For this, the proposed mixer is the best candidate thanks to its negligible power consumption. Table 1 presents the key performances of the proposed differential 25% duty-cycle up-conversion passive mixer and compare them with published works. The proposed upconversion passive mixer performs higher OIP3 than works in [17][18][19][20][21] while consuming only 1.5pW. Regarding conversion gain this work obtains À3.68 dB. This value is greater than the gains from passive mixers published in [13] and [19]. This upconversion passive mixer achieves the best ports isolation performances.

Conclusion
We have presented analysis and simulation of a differential upconversion passive mixer with 25% duty cycle intended to be used in UWB transmitter for WBAN sensor node. We proposed an equivalent LTI model of this mixer then we studied the implication of this model in the design. After that we analyzed the influence of the different design parameters on linearity, conversion gain and noise figure. We studied the impact of voltage threshold mismatch on linearity and ports isolation performances. Using the results found in this paper, RF circuit designers can quickly optimize the different parameters of the up-conversion passive mixer without going through time consuming transient simulations. This results in a guidance for optimizing low power UWB transmitter performance for WBAN applications. In particular this UWB transmitter can be used to design wireless sensors worn by corona virus patients thereby limiting patient-to-caregiver contact.


extfieldsofstudy:
Computer Science
Medicine

text:
Fusion of Range and Stereo Data for High-Resolution Scene-Modeling

This paper addresses the problem of range-stereo fusion, for the construction of high-resolution depth maps. In particular, we combine low-resolution depth data with high-resolution stereo data, in a maximum a posteriori (MAP) formulation. Unlike existing schemes that build on MRF optimizers, we infer the disparity map from a series of local energy minimization problems that are solved hierarchically, by growing sparse initial disparities obtained from the depth data. The accuracy of the method is not compromised, owing to three properties of the data-term in the energy function. Firstly, it incorporates a new correlation function that is capable of providing refined correlations and disparities, via subpixel correction. Secondly, the correlation scores rely on an adaptive cost aggregation step, based on the depth data. Thirdly, the stereo and depth likelihoods are adaptively fused, based on the scene texture and camera geometry. These properties lead to a more selective growing process which, unlike previous seed-growing methods, avoids the tendency to propagate incorrect disparities. The proposed method gives rise to an intrinsically efficient algorithm, which runs at 3FPS on 2.0MP images on a standard desktop computer. The strong performance of the new method is established both by quantitative comparisons with state-of-the-art methods, and by qualitative comparisons using real depth-stereo data-sets.

I. INTRODUCTION
Many computer vision methodologies, including dense 3D reconstruction [1], [2], gesture recognition [3], [4], and object detection [5] have benefited from recently-developed depth sensors. These sensors rely on active-light principles, including modulated-light and pulsed-light cameras, commonly denoted time-of-flight (TOF) [6], [7], or projected-pattern triangulation cameras [8]. Regardless of the working principle, however, these sensors provide low-resolution (LR) or mid-resolution depth maps that are inadequate for a number of applications such as 3DTV and film production. For example, many tasks in the film production industry can greatly benefit from a highresolution (HR) and high-quality depth map [9].
While HR depth maps can be obtained from multiple-view matching and reconstruction using standard color cameras, it is well known that stereo matching is problematic when the scene contains weakly textured areas, repetitive patterns, or occlusions; these situations are very common in both indoor and outdoor environments. Active-light sensors do not suffer from these limitations, although their own depth data are quite noisy in the presence of scattering, non-Lambertian materials, and slanted surfaces. The complementary nature of HR stereo and LR depth sensors leads to the design of mixed camera systems [10]- [18], which seem to be the most promising approach, at present, for high-quality 3D depth maps.
In this context, this paper addresses the problem of HR 3D reconstruction from the combination of a photometric camera pair and an active-light camera, provided that the multiplecamera setup is calibrated [19]. The combination of a stereo matching algorithm and of an active-light sensor raises the central question of devising a matching algorithm with the following features: (i) it considerably increases the resolution of the depth data, e.g., by a factor of ten, (ii) it eliminates depth-sensor errors wherever possible, (iii) it overcomes the limitations of stereo algorithms in textureless areas, and (iv) it is able to compete with a depth sensor in terms of speed. Hence, the availability of an efficient and robust stereo algorithm that takes advantage of LR depth sensors and that provides dense and accurate HR depth maps, possibly with subpixel accuracy, is particularly desirable.
To this end, we propose a 3D reconstruction method that merges depth-sensor measurements with photo-consistency stereo matching. We address the problem from the perspective of seed-growing, starting from a small number of control points whose disparities are then propagated to yield a dense disparity map. We show that this can be cast into maximum a posteriori (MAP) formulation (Sec. III), which leads, in turn, to a series of local optimization problems that are solved hierarchically by a novel region-growing process (Sec. V). While the proposed method may not reach the global optimum, it allows us to devise an intrinsically efficient methodology that bridges the gap between global optimizers based on Markov random fields (MRF) and locally-optimal winner-takeall (WTA) strategies (Sec V-E).
Efficient stereo-only or stereo-depth fusion methods rely on control points, by exploiting either feature correspondences in stereo [20], [21] or depth data in fusion [14], [15], and they assume that these points are of very good quality. A key contribution of this paper is that this requirement is relaxed in order to devise a method tolerant to bad control points. We propose to truly combine LR depth-sensor data with HR rich photometric information, whenever and wherever possible, showing that fusion is helpful, even in the early stage of depth initialization (Sec. IV-A & IV-B). The data term of the proposed MAP formulation benefits from a new cross- Fig. 1: The pipeline of the proposed depth-stereo fusion method. The low-resolution (LR) depth data are projected onto the color data and refined to yield a high-resolution (HR) sparse disparity map. Starting from these disparity seeds, an upsampling process provides an initial HR dense disparity map. Both the HR seeds and the initial dense disparity map are then used by the region-growing depth-stereo fusion to produce the final HR depth map. A prominent feature of our method is that fusion takes place at several data processing stages. correlation function, which provides real-valued disparities via subpixel corrections computed in closed-form (Sec.V-B1), and takes advantage of a depth-guided cross-product aggregation (Sec. V-B2). Moreover, the data-term adaptively merges the stereo-and depth-consistency terms guided by the scene texture and the camera geometry (Sec. V-C). These advantages lead to a more selective growing-of-disparities process, thus preventing the algorithm from propagating erroneous depthsensor data -a phenomenon that is often associated with propagation techniques -as experimentally verified (Sec. VI). It is important to note that the proposed method can be used 'as is' with any depth sensor, as it requires neither sensordependent confidence maps nor a sensor-dependent model. The proposed fusion pipeline is illustrated in Fig. 1.
Supplementary materials, in particular image datasets and Matlab code are available online. 1

II. RELATED WORK
We review pure stereo matching methods with emphasis on local algorithms, owing to their computational suitability for HR images. We also review upsampling methods and depth-stereo fusion methods. More detailed surveys of stereo matching algorithms and depth-stereo fusion methods can be found in [22] and [9], respectively.

A. Pure Stereo Matching
Stereo matching methods can be broadly classified into global and local [22]. Global algorithms [23] typically adopt an MRF formulation and solve a single optimization problem based on a MAP criterion. Despite their superiority over local methods, global algorithms are extremely time-consuming, and hence unattractive for fast fusing of depth-sensor data with high-resolution images. Local algorithms solve per-pixel optimization problems and the state-of-the-art methods build on adaptive cost aggregation [24]- [27]. Most methods, however, must visit the entire cost volume to find an optimal disparity value at each pixel. This volume grows rapidly with respect to the input, as the width of the image typically multiplies the number of disparities. Therefore, although they are able 1 https://team.inria.fr/perception/research/dsfusion/ to provide LR disparity maps in real-time, they remain slow and subject to memory issues in HR stereo. Note that global algorithms need several approximations to obtain LR disparity maps in real-time [28].
More interestingly, algorithms that rely on control-point correspondences [29] are drastically more efficient since they avoid visiting the whole cost volume. Region-growing approaches start from reliable but sparse correspondences (seeds) and propagate them in textured areas [20]. [30] suggests a similar propagation scheme where orientation-consistent disparities are propagated to neighbors at the cost of finding a plane equation per pixel. [21] proposed a generative model, where the prior disparity comes from a 2D triangulated mesh whose vertices (control points) are obtained from matches between low-level features. As with [20], textureless areas remain intractable and the final map is reliable only when the matches are dense and uniformly distributed over the images. Notice that the proposed method relies on the idea of control points that are transferred from a depth sensor, thus avoiding the limitation owing to untextured areas.
To obtain continuous disparities that are required in many scenarios, e.g., 3D reconstruction [31], local stereo algorithms typically employ two strategies: (i) fitting a curve around the correlation peak [14], [22] or (ii) integrating an intensity interpolation kernel into the (dis)similarity function [32] whose optimization leads to subpixel correction. The latter is also the case in the fusion framework of [16], [17] that inherently takes advantage of inter-pixel depth estimations.

B. Depth Upsampling and Depth-Stereo Fusion
Any prior depth information, even at low resolution, is likely to help dense disparity estimation. Apart from a naive interpolation, the bilateral filter [33] can post-process an interpolated map using color HR images for guidance [31]. Alternatively, a joint bilateral filter applies spatial and range kernels to the LR depth map and HR color image respectively, so that upsampling is a by-product of filtering [34]. Upsampling methods, however, are limited to clean depth LR data and cannot reconstruct accurate HR maps when LR data are delivered by depth sensors.
The above limitation of the upsampling methods gives rise to fusion approaches that can merge depth and stereo data in either early or late stages, once the sensors are calibrated. Late fusion suggests merging two depth maps, one obtained with stereo and one from the range sensor, possibly upsampled [18], [35]. The majority of fusion methods, however, merge the stereo and range-sensor data at an earlier processing level. [10] estimates a LR TOF-based disparity map, which initializes a coarse-to-fine stereo scheme [36]. A semi-global dynamic programming approach is followed in [14] with the TOF-based disparities, wherever available, being considered as error-free matches. In a global framework, [11] produces mid-resolution depth maps by merging TOF and stereo data within a graphcut stereo-matching paradigm; each energy term exploits both modalities. Likewise, MRF-based formulations have been proposed [12], [13], [37]. In [37], ground control points reflect an extra regularization term in the MRF energy function. The work described in [12] uses an MRF scheme to merge depth distributions of each sensor alone, but the goal is a LR depth map. [13] extends [12] by means of weighted fusion. A balanced fusion (based on several confidence maps) within a total variation model is also proposed in [17]. A similar variational model that infers the HR map in a coarse-to-fine manner is adopted by [16]. Note that, unlike most fusion methods and similar to the proposed one, [16] is not tuned to a specific depth sensor.
More closely to the present work, [15] fuses the data within the seed-growing method of [20]. In particular, TOFbased disparities constitute seeds while a triangulation-based interpolated TOF map regularizes the seed-growing process. When the TOF data are noisy, however, this approach tends to produce incorrect disparities, and to propagate false positives during the growing process. The proposed method differs considerably from [15] in terms of depth initialization, cost function, and fusion strategy. The proposed initial map is robust to depth discontinuities while it also guides the cost aggregation inside a window. Moreover, our likelihood term integrates functions that are capable of providing sub-pixel corrections [32]. This turns out to be very beneficial, not only for the continuous nature of the final map, but also for the growing process itself, thus propagating high-confidence disparities only. Note that the subpixel disparity correction is obtained from a closed-form solution -an interesting feature for efficiency. Finally, our algorithm benefits from an adaptive fusion scheme that better balances the contribution of each modality (depth or color), and that results in fewer unmatched pixels.

III. PROBLEM FORMULATION
The main mathematical notations that are used throughout the paper are summarized in Table I. As discussed, the direct upsampling of LR depth data suffers from limitations, in particular when the upsampling factor is high. 2 Therefore, our goal is to build D by jointly taking advantage of both sensing modalities. Given a proper calibration, (e.g., [19], [38]), the  with D, D 0 and D↓ respectively, with dp ∈ D, d 0 p ∈ D 0 and d p↓ ∈D↓ tp: Subpixel disparity correction of p µp = (dp; tp): Disparity-correction pair referred here to as metadisparity with |tp| < 1 M = {µp}: Set of meta-disparities S = {sp}: Set of observed stereo pixel intensities Ip: Intensity of pixel p u(x, y): a vectorized (zero-mean) form of an intensity window centered at the 2D position (x, y) Np: Neighborhood of pixel p med: 2-d median operator I R , D R : Intensity and disparity maps defined on a sub-region R E S (dp), E D (dp): Stereo-based and depth-based energy of p for given disparity g(ξ; γ) = e −|ξ|/γ : Exponential mapping of ξ with scale γ mapping of the LR depth image onto the rectified HR color images will typically yield a sparse disparity map D↓, which is almost evenly distributed across the HR grid. Note that in [20], [21], the initial matches between the stereo images do not correspond to a uniform sparse version of D, as they are unpredictably distributed, due to the reliability of properly detecting interest points in images. Instead, the sparse map obtained with a depth sensor can be used to guide a stereo algorithm, regardless of the presence or absence of scene texture.
We propose to model the estimation of D, and therefore the map D, as a maximum a posteriori (MAP) problem, based on the available depth and stereo observations. However, instead of immediately using D↓, we first estimate a dense initial map D 0 and its associated set D 0 . Then, we obtain the final disparity map by solving the following optimization problem: where P (D|S, D 0 ) is the posterior distribution of disparities given the observations S and D 0 . The proposed depth initialization method is described in Sec. IV and the proposed solution to the MAP formulation (1) is described in Sec. V. Note that a reliable estimation of D 0 is quite important, since it guides several components of the fusion methodology.
IV. DEPTH INITIALIZATION A two-step approach is proposed in order to obtain the initial disparity map D 0 . First, we refine D↓ to deal with mapping errors. Second, we upsample the refined sparse map in a novel way using color information to obtain the initial dense map D 0 . As shown below, this leads to an initialization robust to depth discontinuities, which in turn helps the growing.

A. Sparse-Depth Refinement
We assume a camera setup with a depth camera mounted in between the two color cameras; other depth-stereo setups are (a) (b) (c) Fig. 2: (Best viewed on screen) (a) The mapping of depth data onto the left image causes artifacts in the presence of depth discontinuities. A cascade of (b) geometry-consistency and (c) color-consistency filters refines the sparse disparity map. Depth values are color-coded from red (close) to blue (far). Fig. 3: (Best viewed on screen) Depth upsampling results using (a) triangulation-based interpolation [15] after cutting big triangles, (b) joint bilateral filter [34] and (c) our method. The depth values are color-coded from red (close) to blue (far), while black areas correspond to non-available values. The white edges in close-ups show the color edges of the image. discussed in [9]. Regardless of the depth sensor technology and type, the projection of the depth map onto the left and right images implies a parallax effect, and hence occlusions. Moreover, this causes gaps as well as areas with overlapping depth data close to depth boundaries [15]. In the case of TOF cameras, these areas are further contaminated from jumpedge errors [39], or 'flying pixels' [40], while a structuredlight camera, e.g., Kinect, leaves more gaps due to the offset between the position of the light projector and the position of the infrared sensor. Fig. 2(a) illustrates the artifacts that we briefly discussed: flying pixels and depth-data overlap in the top and bottom closeups, respectively. In order to eliminate these artifacts we apply a geometry-consistency cascade of two filters: the first one removes isolated pixels (mostly flying pixels) and the second one keeps the foremost pixel inside a window to compensate for the above mentioned overlap. An example of applying this filtering is shown in Fig. 2(b). In practice this does not fully refine the sparse depth map. We still observe mismatches near depth discontinuities, because of depth bias and calibration errors. Therefore, a second filter that imposes color consistency is applied, as described below.
We consider a window centered at p↓ and split into four equally sized sub-windows W i , i = {1 . . . 4}, such that their intersection is only the pixel p↓ (see Fig. 4). The output, d p↓ , of the filter is: with The output d p↓ is the median disparity of the adjacent subwindow whose median intensity is closest to that at p↓. For color images, the term |I p↓ −med(I Wi )| can be replaced by the average deviation from the median, over the color channels. This filter leads to a further refinement near depth discontinuities which are pathological areas for stereo algorithms. The result of this kind of filtering is shown in Fig. 2(c). Note that both refinement filters apply to sparse locations only so that their complexity is negligible. Fig. 4: The window split for the color-consistency filter. The pixel p is linked with the closest (shaded) sub-window in terms of the color consistency (links represent color distances from I p to med(I Wi )).

B. Sparse-Depth Upsampling
While a naive upsampling of the sparse disparity map could be performed, e.g., [15], thereby producing an initial map, strong depth discontinuities are likely to contaminate such an interpolation. Alternatively, cross-bilateral filtering [41] or joint bilateral upsampling [34] may upsample the map using the HR color image as a guide. The filter support of these methods jointly counts on spatial and range (color) kernels. However, both methods apply a linear smoothing once the filter support per pixel is computed. Instead, we propose a non-linear upsampling strategy that builds on the median filter.
Specifically, the depth (or disparity) at a dense pixel location p is initialized by where N c p is a constrained neighborhood of p, that is N c p ⊂ N p , which contains only sparse depth measurements whose color is consistent with I p : with g(ξ; γ) being an exponential mapping (see Table I).
Unlike common bilateral filters, our upsampling process makes a more definitive selection of pixels, thus preserving the depth discontinuities of the scene, while filtering out some of the noise in the depth data. Once the N c p is defined, one can optionally consider a spatial kernel and compute a weighted average instead, in order to better deal with slanted surfaces. Fig. 3 compares our initialization with the upsampling results obtained by [15] and [34]; in the detailed views, the intensity edges are also shown. The proposed method provides more discriminative depth boundaries that coincide with intensity edges. Missing values may be observed in highly textured areas, since N c p may be an empty set. In this example, the radius of N p is 20, γ c = 10, and c = 0.2. The same radius is used for the method of [34] while the scales for spatial and color kernel are 10 and 20 pixels, respectively. The geometryconsistency filter reasonably applies in all cases while our method benefits from our color-consistency refinement as well.
Since the median operator is chosen to account for outliers within a window, the mean operator can be invoked instead when the depth variance almost vanishes (homogeneous areas), thus drastically reducing the computational burden of the upsampling process. Since the vast majority of pixels belong to such areas, the complexity of our filter approaches that of joint upsampling filter [34]. Note that the latter has been extended in [18] by integrating color segmentation results.

V. DEPTH-STEREO FUSION
Let d p ∈ D have N possible discrete states; the goal is to estimate the disparity (state) of each HR image location through the MAP formulation (1). Once D 0 has been initialized, one can assume that S and D 0 are conditionally independent, so that the posterior distribution of (1) can be decomposed as As mentioned, global solutions are prohibitively expensive for high-resolution disparity maps. Therefore, we focus on approximate solutions that allow for the decomposition of the global optimization problem into many local (per-pixel) optimization problems. The proposed method is based on the seeded region-growing framework [20], [42], where the known message of a location (the parent) is propagated to its neighbor (the child). This implies that the estimation of d p is also conditioned by a parent known disparity d pa p , hence dealing with the principal where d pa p is the parent of d p and the probability has a uniform distribution. In other words, the disparity range for each node is a function of the assigned disparity of his parent. We consider a narrow support area, i.e., a constant low value for r, e.g., 1 or 2; node-dependent parametrization of r is left to future work.

A. MAP as Energy Minimization
As is customary, likelihoods are chosen from the exponential family which leads to an energy minimization framework. Our model assumes an energy-dependent distribution (Boltzmann) for P (s p |d p ) and a Laplacian one for P (d 0 p |d p ): Based on (7-10), Fig. 6 shows an example with the distribution of d p being constrained by single or joint observations. Because of the exponential terms, the pixel-wise maximization of the posterior distribution reduces to the minimization of the local energy where is the (stereo) data-term which is defined below. The term E D (d p ) guides the inference in textureless areas, while it penalizes mismatches due to depth discontinuities when the latter are well preserved in D 0 . Notice that a smoothness constraint is implicitly enforced because of the prior term P (d p |d pa p ) owing to the low value of r. In other words, the support area of E(d p ) is truncated as shown in Fig. 6, and the pixel-wise optimization problem becomes The order of visiting pixels and solving (12) obeys a mostconfident first-solved rule, as discussed in Sec V-D. Below, once the data-term is defined, we modify (11) to adaptively combine the data-term with the regularizer.

B. The Data Term
The proposed data-term benefits from two properties: its ability to provide energy for subpixel disparities and to adaptively consider pixel-wise similarities within a window. Apart from the advantage of a locally continuous map [31], it is more important to note that these properties lead to a more selective growing in areas with texture (owing to the former) and varying depth (owing to the latter) because of more accurate and reliable local energies.
1) Subpixel disparity: To obtain energy at interpixel locations, we make use of meta-disparities (d p ; t p ) (see Table I). The estimation of t p builds on [32], which has been shown to be superior to the parabola-fitting around the peak correlation: the pixel p = (x, y) of the left image is matched with the subpixel p = (x + d p + t p , y) of the right image by defining a bounded correlation function C dp (t p ) (|C dp (t p )| < 1) which is maximized with respect to t p for a given d p . This allows us to define our data-term as where Consequently, the total energy (11) is parameterized by t * p as well. Note that if t * p can be estimated via an analytic solution, the extra computational cost of subpixel disparity estimation is negligible. It becomes clear that, since C dp (t * p ) ≥ C dp (0) in (13), more accurate energy values are estimated and more reliable messages are propagated.
Since invariance to photometric distortions is important in stereo matching, we propose to adopt the normalized correlation coefficient, and one of its variants [43], for C dp (t p ). While the former has been already extended to deal with subpixel disparities [32], the latter has never been extended before. Both use zero-mean vectorized forms of the windows around p and p , let u L (x, y) and u R (x + d + t, y), with the latter being written via a first-order Taylor approximation as where ∆ is a difference operator along the x−axis. This is also the case in variational framework [16], [17] where subpixel accuracy is obtained by the early interpolation of the intensity, rather than the late interpolation of correlation around the peak [22].
Enhanced Correlation Coefficient (ECC): The ECC function [32] results from the integration of the above linear interpolation kernel into Pearson coefficient: If the denominator of (15) is non-degenerate, then C P d (t) is a quasi-concave function of t and its maximization results in a closed-form solution [32].
Enhanced Moravec Correlation Coefficient (EMCC): The Moravec coefficient [43] replaces the denominator of (15) with the mean of the variances. This also allows us to introduce a left-right symmetry, thereby estimating a left disparity −t/2 and a right disparity t/2 instead of t (such a modification with ECC leads to a complex optimization problem). The enhanced Moravec correlation coefficient (EMCC) is defined by Note that one can easily show that the integration of an interpolation kernel into the cost function of [15] is equivalent with the EMCC scheme. Although (16) is a rational function of t, the next proposition guarantees that the maximizer has an analytic form. We refer the reader to the appendix for the proof and the exact maximizer.
Proposition I: A rational function of two second-degree polynomials, as in (16), attains at most one global maximum if the denominator is non-degenerate; its maximizer is given by a closed-form solution.
Note that, the estimation of t could be unreliable for strictly homogeneous areas. The value of the window variance or entropy is a good criterion to assess the reliability of subpixel correction, and to enable it accordingly.
2) Adaptive similarity aggregation: The best-performing local stereo algorithms benefit from an adaptive cost aggregation strategy [24]. This strategy is based on the assumption that depth discontinuities are most likely to coincide with color discontinuities, so that each pixel within a window contributes differently to the (dis)similarity cost based on its spatial and color distance from the central pixel. However, only a few color edges correspond to depth edges and the above assumption should be followed only in the absence of any prior information about the depth. Since in our scenario the prior depth information is available, the spatial and color consistency can be replaced by a depth consistency term.
To be specific, we adopt here the exponential g(·) (see Table I) to compute pixel-wise weights w q : with q ∈ N p . The weights apply element-wise to u L , u R , ∆u R and ∆u L in (15) and (16). In other words, we compute the subpixel correction and the optimum local energy after down-weighting pixels in the window that belong to another surface compared to the one of the central pixel. It becomes clear now that not only the term E D (d p ) but also the stereo term E S (d p ) in (11) benefit from an upsampling method that is robust to depth discontinuities. Note that, even if the initial depth map is biased, it is sufficient enough to guide the aggregation step within the window.

C. Adaptive Fusion
While a constant fusion may be reasonable for specific types of scenes (e.g. highly-textured scenes), an adaptive balance of the terms in (11) is usually preferred. This implies that the less we count on the data term E S (d p ), the more we should count on the regularization term E D (d p ) during the inference. 4 This suggests a convex combination when the scene point of p is viewed by all cameras.
A summary of methods that perform weighted fusion is discussed in [17]. However, most of them consider TOF-based weights for the regularizer (e.g., [13]) which contradicts our goal of a sensor independent fusion. Moreover, directly using the confidence map of a TOF image is not a good strategy [40]. Therefore, we only rely on stereo data to obtain the mixing coefficients. It is well known in stereo or optical flow that the matching of a point is reliable when its associated image patch contains sufficient texture [44], [45]. Since a good indicator for the texture presence is the image entropy [44], an entropy filter provides us with am adequate reliability factor e p for each window centered at p.
Let us now consider the left image as reference and compute the initial left-to-right disparities. Likewise, we can build a right-to-left disparity map based on the right image, and a cross-checking of these maps can provide an estimation of the major occlusions due to strong depth discontinuities, with respect to the reference image. We refer to these areas as stereo-occlusions and we denote them as Ω SO . Recall now that some points in the left image are not seen by the depth camera, and that this gives rise to gaps in the initial disparity map which can be easily detected. We refer to these areas as depth-occlusions and we denote them as Ω DO . It becomes obvious that the evaluation of E S (d p ) and E D (d p ), in Ω SO and Ω DO respectively, should be avoided. Hence, we propose the following adaptive fusion with the pair (η S p , η D p ) being defined as where Ω defines the whole image area and e p is the normalized output of the entropy filter. We intentionally add the last 4 Here we omit the disparity correction t. case in (19) which shows that the fusion in Ω SO Ω DO is meaningless and a post-filling method should be followed.
It is important to note that the so-called ordering constraint in stereo is valid when large foreground objects appear in the scene, while it is violated when very thin objects are close to the camera (see Fig. 7). The former implies Ω DO Ω SO = Ω DO and the latter implies Ω SO Ω DO = ∅. While the area Ω SO \ Ω DO can always be predetected, the area Ω DO \ Ω SO is safely predetected only when the ordering constraint is not valid. 5 This is because Ω SO is detected from the crosschecking of disparity maps that already suffer from depthocclusions, since they are computed from the depth-to-stereo mapping. Ideally, if the complexity is not an issue, a stereoocclusion detection scheme (e.g. [46]) based on stereo image pair could be enabled beforehand. Therefore, we prefer to not grow the area of depth-occlusions. Optionally, these areas can be filled in a post-processing step. Fig. 7(c) shows the areas for the example of Fig 2. Points that have been removed during the refinement step (Sec. IV-A) are also marked as depthocclusions, while outliers in the initial maps can produce false positives for Ω SO , e.g. the red points on the right arm of the person. Our region-growing method, however, is able to compensate for such errors.

D. The region-growing algorithm
As has been explained, our method solves pixel-based optimization problems in a region-growing manner, based on the seeds contained in D↓. Since the initial disparities may be noisy and biased, they do not reflect true matches as opposed to [15]. This means that we exploit d p↓ to restrict the disparity range of its neighbor, but once a disparity value is assigned to the latter, p↓ is reset to a pixel p with unknown disparity. For our convenience, the set D↓ is augmented to an initial set M of meta-disparities (d p ; 0) with the same cardinality. We also denote with N (µ p ) = {µ j p } 4 j=1 the image-based neighbors of µ p , that is, µ j p = (dp; tp) withp being any of the four immediate neighbors of p. Note that µ j p does not necessarily belong to M.
Algorithm 1 describes the growing process. The algorithm starts by sorting the elements of M based on their energy 5 It would be possible to detect all areas if one knows a priori that the constraint is everywhere valid or invalid. Set f v (µ p ) = true 8: for all µ j p ∈ N (µ p ) with f a (µ j p ) = false do 9: µ j * p = arg min E(µ j p ) (Eq. 20) 10: if E(µ j * p ) < T then 11: Set f a (µ j * p ) = true 12: Push µ j * p in M w.r.t. sorting 13: end if 14: end for 15: until card(M) is fixed 16: Compute the dense disparity map D from M 17: return D.
value, while it initializes to false the visit-and assignmentflag of all candidate meta-disparities. Next, it considers the lowest-energy µ p with false visit-flag and switches this flag to true. Then, it assigns values at each µ j p ∈ N (µ p ) with false assignment-flag based on the following minimization scheme where E(dp; tp) = η Ŝ p E S (dp; tp) + η D p E D (dp) (21) and d * p is the optimum disparity of the parent node ofp. Note that (21) extends (11) by adding the subpixel disparity parameter into the stereo term, and making the fusion pixeldependent. In other words, what we do for each neighborp with false assignment-flag is the following. For each candidate integer value dp ∈ [d * p − r, d * p + r] the optimum t * p that minimizes the term E S (dp; tp) is obtained based on (16) or (15) and the local energy E(d p ; t * p ) is computed from (21). Among the 2r + 1 values, the disparity minimizer d * p is finally chosen and the corresponding subpixel correction is assigned, as it is shown in (20). Recall that r has a low value in contrast to conventional stereo algorithms where it equals the whole disparity range. If E(d * p ; t * p ) < T , then the meta-disparity (d * p ; t * p ) is pushed into M with respect to the sorting, and its assignment-flag becomes true. The above process is repeated with unvisited meta-disparities until the cardinality of M remains fixed. The visiting order depends on local energies, since the lowest-energy meta-disparity is always picked from the stack. The final set M corresponds to the final dense disparity map while a post-filling method can deal with missing disparities; their number heavily depends on threshold T . It is important to note that the algorithm cannot get stuck in a loop, because it propagates disparities in a tree structure, and a true visit-flag can never be reset to false.

E. Comparison with Other Inference Formulations
State-of-the-art stereo or fusion methods adopt an MRF model that provides a straightforward way to integrate multiple sensor data. If we recall equation (6), the term P (D) can be written as a Gibbs distribution whose energy is a sum of potential functions over maximal cliques [47], hence a sum of pairwise potentials when a 4-pixel neighborhood is considered. This in turn offers a smoothness term in the global energy equation that leads to piecewise smooth disparity map. There exist exact solutions for such models, under very specific conditions [23]. Their computational complexity, however, becomes more and more severe as the number of states increases, thus becoming prohibitively expensive in the case of HRstereo. Therefore, one has to focus on feasible approximate inferences that decompose the global optimization problem into a series of local optimization problems.
MRF stereo generalizes the winner-take-all (WTA) approach [22] which can be seen as the simplest inference. Fig. 8(a) and (b) show an 1D-grid example of the MRF and WTA networks. In essence, links between disparity nodes in WTA disappear and each node is connected only to input data (observations). This implies a uniform prior distribution P (D) and the maximization of P (D|S, D 0 ) reduces to a set of independent pixel-wise maximization problems of type P (d p |s p , d 0 p ). As with the proposed method, we tacitly assume, that s p are not intensities but they represent the stereo data of local windows centered at p and its candidate correspondence. This is necessary in WTA since the likelihood cannot count on single pixels only, while it can optionally used in global MRF solutions as well. It is worth noticing that our algorithm switches to a WTA solution when r = N , since the visiting order of nodes and their connectivity becomes meaningless and the uniform distribution P (d pa p |d p ) is not truncated anymore. Another MRF approximation that finds pixel-wise optimizers without breaking the connectivity between nodes is the iterated conditional modes (ICM) method [48]. After a proper initialization of all disparity nodes, ICM-like schemes visit one node (pixel) at a time and assign the disparity value that maximally contributes to the global posterior distribution. On the contrary, our approach requires the initialization of few nodes only, at least one in principle. Moreover, the graph in our case is a set of directed trees, i.e., a forest, so that each non-root node has only one parent. It is important to note, however, that it is the output of our algorithm that defines the final network since nodes may be re-ordered, links may be cancelled and arrows may be reversed during the inference.
To be more specific, let us consider the graph of Fig. 8(c) and let us assume the clamping of d 2 and d 5 to the observation d 0 2 and d 0 5 . Starting from d 2 , we can 'propagate' d 0 2 to its neighbor d 3 (and d 1 ), i.e., to look for the optimum value of d 3 but being strictly conditioned by the value d 0 2 . Note that ICM would look for the best assignment of d 3 value by taking into account both initial values d 0 2 and d 0 4 . Moreover, ICM would search, in principle, among all (here N ) states for the optimum d 3 's assignment, while our scheme looks only around d 0 2 , namely d 3 takes values in {d 0 2 − r, ..., d 0 2 + r}, with r being a small integer. Once d 3 's assignment is done, d 0 2 can be passed to d 1 in a similar way. Next, another node is visited, here one of d 3 and d 5 , and its disparity is propagated to its neighbors. As a result, the principal graph of Fig. 5 is iteratively considered. Note that ICM assigns a new value to d 3 anyway, while our scheme invokes a criterion that validates the assignment. If d 3 's assignment is not valid, our algorithm will possibly assign a new disparity value to d 3 only after d 4 's assignment. Fig 8(d) and (e) show possible final graphs obtained by different realizations of our algorithm. In the example of Fig 8(e), the initial disparity of d 5 was cancelled by the validation process, thus all nodes were filled due to d 2 .
It is now clear that, as opposed to WTA solution, the final inference we obtain depends on the visiting order of nodes. The visiting order of the ICM scheme is either fixed in advance (e.g. raster-scanning), or random. Inspired by [47] and [20], however, we instead adopt a highest confidence first (HCF) scheme that suggests visiting the nodes based on their local evidence (energy). This means that we keep the nodes sorted with respect to their energy, and we visit each time the leastenergy node that has not been visited yet. For instance, in Fig. 8(d), starting with d 2 implies that d 2 is more confident than d 5 . The above validation process relies on thresholding the local energy as explained in Sec. V-D. Table II summarizes some properties of MRF-based solutions (Graph-cut [49] and ICM [48]), WTA, and our approach.

VI. EXPERIMENTS
In this section, we evaluate our algorithm, and quantitatively compare it to the state-of-the-art, based on both simulated and real data-sets. We also test our algorithm and provide qualitative comparison on a new and challenging dataset.

A. Simulated Data
We use the Middlebury database, and focus on a challenging data-set which contains 1.5MP images (1300 × 1100) along with ground-truth maps (GTM) [50]. To simulate an LR disparity map from another viewpoint, we proceed as follows. Given the calibration parameters, we transform the groundtruth disparities into 3D points, as viewed from the midpoint of the baseline. We then apply a 3D rotation to the pointset and we downsample the points by a factor of 10. The rotation is such that the average disparity bias is more than 2 pixels. Finally, we back-translate the points into sparse biased disparities and we add colored noise, that is, a 2D midfrequency sinusoidal signal with peak-to-peak distance equal to 2σ, where σ denotes the noise deviation. Note that [15], [31], [34] only downsample the GTMs in their experimental setup. Algorithm performance is quantified in terms of the so-called bad matching pixels (BMP) percentage in the nonoccluded areas [22], i.e., (1/N o ) p (|D p − G p | > δ) where G is the GTM and N o is the number of the non-occluded pixels. While the threshold level δ = 1 is mostly used for mid-resolution images, HR stereo justifies the value δ = 2 as well [21]; a value δ < 1 is chosen when subpixel accuracy is to be evaluated.
The same parameter settings are used for our method, in all of the experiments. The radius of the upsampling filter is 20 and the values γ c and e c are 10 and 0.2 respectively. Because of the propagation strategy, we choose a relatively small window, i.e., 9 × 9. The local energy in (11) and the weights in (17) are obtained with λ = 0.01 and γ d = 5, respectively. We enable subpixel correction when the normalized entropy in the (left) window is above 0.4. As mentioned, we use a fixed (and strict) search range around the disparity parent, that is r = 1, which leads to the most efficient solution. As for the threshold, we set T = 0.5. Note that the energy validation threshold implies a tradeoff between accuracy and density. We recommend setting a middle threshold value and post-filling sparse missing disparities, e.g., with the upsampling filter, rather than using a high threshold that incorporates erroneous disparities in a fully dense map. The density obtained with this strategy is about 90% in HR images. As with all algorithms, large remaining gaps are filled with a streak-based filling [21]. We refer to our methods as Fusion-ECC (F-ECC) and Fusion-EMCC (F-EMCC).
Before comparing with the state-of-the-art, we show the performance gain in terms of the new modules that are integrated compared to EPC method [15], thus quantifying the contribution of the new data-term and the adaptive fusion. Note that [15] (EPC) follows a seed-growing approach by using a quadratic model for both stereo and depth consistency terms respectively. To better evaluate the contribution of the enhanced correlation coefficients presented in the data-terms, we use LR images (450 × 375 on an average) whose GTM contains subpixel disparities. The noise deviation in the sparse map is σ = 2. We intentionally do not fill any remaining large gaps, in order to assess the net contribution of each module, and we compute the error for the filled area only (85% density). Table III shows the BMP error averaged over eight images. We also evaluate our approach when none of the modules are enabled, i.e. pure correlation is used along with a fixed fusion of the terms; we refer to this method as simple fusion. All of the variants start from the same initial depth map, obtained by our upsampling method. As can be seen, both the data-term and the adaptive fusion process contribute to a better reconstruction, compared to simple fusion. The use of the proposed energy data-term leads to more accurate  [23], [49] results, while the adaptive fusion eliminates large errors (see the error with δ = 2). Even the simple fusion has a lower BMP than EPC, owing to the different stereo-and depth-consistency terms (we use a linear model for the latter). The two proposed criteria behave similarly, with the F-ECC being slightly more accurate, since it achieves higher correlation values (see also Table IV). Note that the results are systematically worse if we use the initialization of [15].
To compare with the state-of-the-art, we implemented the upsampling methods of [34] and [31] (two-view version), as well as the MRF-based fusion [12] by using the MRFstereo toolbox of [23], referred here to as F-MRF. While [12] uses belief propagation, we experimentally found that Graph-Cuts [49] perform better. Specifically, we tried ten different parameter settings and we found that the best performing algorithm is the expansion mode with Birchfield-Tomasi cost [51] truncated at 7, linear disparity differences truncated at 5 and quadratic cost for the smoothness-energy; the weights for the stereo-, depth-and smoothness terms were 1, 1.2 and 10. We refer here to [12] instead of [13] for MRF-based fusion since the latter relies on a TOF-based reliability fusion which cannot be implemented here. Fig. 9 plots the BMP curves of the upsampling and fusion algorithms as a function of the noise deviation for the challenging low-texture HR image Lampshade1. We just add noise here in the down-sampled GTMs. Except for [34], all schemes start from the same HR map, obtained by a naive interpolation, while its BMP curve is plotted as well. As can be seen, F-MRF and EPC are more affected by the noisy prior disparity, in contrast to the proposed algorithm, which is less sensitive to initialization. It is clear that the pure up-sampling methods provide acceptable results only when the initial LR disparity map is very accurate.
We now proceed with a detailed comparison including efficient and well-known stereo algorithms as well. Specifically, we include four recently proposed methods, [21] (ELAS), [20] (CGS), [25] (FastAgg), [27] (NonLocalAgg) and two MRFbased stereo algorithms, Graph Cuts (GC) [52] and constantspace belief propagation (CSBP) [28]. The top-performing lo- cal algorithms, FastAgg and NonLocalAgg, build cost volumes that depend on both the image size and the disparity range. This leads to a huge memory footprint (∼ 3GB) in the case of HR images, and the authors' implementations could not be run as is. In order to be able to run FastAgg, the cost volume has been split into slices and cached on disk. The NonLocalAgg was run using the maximum allowed resolution while the disparity map produced by the algorithm was finally upsampled. Note that both methods invoke a left-right consistency checking, combine color and gradient information and enable refinement steps. Authors' implementations for ELAS, GCS, GraphCut, CSBP (local-minima version+bilateral postprocessing) were used in the comparisons, with the default settings suggested by the authors. We also implemented the ICM algorithm for the MRF-based fusion using the same parameters with GC. For a fair comparison, all fusion schemes merge stereo data with the same initial disparity map, which is obtained here by our upsampling process. Due to the simulated experimental setup, however, the error of the initial map obtained from this process is very close to that of Kopf et al.'s method [34] (the average difference of their BMP error is below 0.5) and is thus omitted. Table IV provides the BMP error for eight HR images with error threshold δ = 1, while the corresponding table for δ = 2 is given in the appendix. Bold and underlined numbers mark the lowest and second lowest errors per column. Weakly textured scenes (Lampshade1, Monopoly) seem to be problematic for conventional stereo algorithms, while cylindrical surfaces (Bowling2, Baby1) present another challenge. It is not surprising, however, that stereo methods outperform fusion methods in highly textured images (e.g. Rocks2), or in images with many thin objects (e.g. Art), since the sparse noisy initial map negatively affects the fusion. ELAS and FastAgg behave better than other stereo algorithms. Similar results would be expected from the NonLocalAgg method, if we were able to run it at full resolution.
Stereoscopic and depth data are better fused in general by the proposed criteria than the other fusion methods. EPC  and up-sampling methods verify the sensitiveness to their initialization, with the former being more effective due to the fusion process. As far as the MRF solution is concerned, the benefit due to the depth data is verified from the results, i.e. the F-MRF (GC) scheme behaves better than pure-stereo GC. Moreover, GC in fusion provides better results than the simple ICM algorithm. Unlike the proposed criteria, F-MRF deals better with thin objects, since it does not aggregate costs in a window, hence the lower error in Art and Reindeer. Recall that [12] refines the TOF depth map based on stereo data, without increasing the resolution. By putting aside the high complexity, it seems that MRF-based solutions need to be reformulated for HR stereo-depth fusion, e.g. high-order connectivity might be more helpful, semi-global solutions could be investigated, and conditional random fields [50] might need to be extended to the fusion framework.
Between the two proposed criteria in our fusion scheme, ECC and EMCC, it is the image content that makes one outperform the other. While F-ECC may be slightly better on average, F-EMCC deals better with images of very low texture, e.g. Monopoly, which supports Moravec's argument for introducing MCC [43]. Moreover, F-EMCC is more affected by the filling, as it provides less dense maps than F-ECC, provided that the threshold is the same. Fig. 10 shows the average performance of all fusion competitors as a function of the error threshold δ, while the stereo baseline of FastAgg is added for reference. In essence, this figure reflects the distribution of errors. Evidently, the contribution of the LR depth prior in the fusion schemes is verified, as opposed to the stereo baseline whose performance is bounded. The proposed schemes are more accurate compared to the fusion baselines. However, F-MRF provides lower errors when the tolerance is not that strict (δ > 2) owing to its global smoothness constraint. Note that the performance of [34] approaches the performance of the fusion schemes as δ is increasing. The other upsampling method of [31] seems to produce large errors.
Computational efficiency is an important feature of any depth-stereo fusion method. Table V shows the execution times of the algorithms for the simulated data, as well as some of their implementation details. A combined Matlab-C version of our fusion algorithm requires 2.0s per image triplet while it takes 0.3s (for 2MP images) when GPU hardware is used for some initializations, and the SSE instruction set accounts for the online computation of correlation values between windows. We also developed a GPU-based implementation of a triangulation-based upsampling (interpolation) that takes less than 50ms. This allows one to envisage real-time execution of the proposed depth-stereo fusion framework, despite the high resolution. Note that it is the Matlab implementation that makes FastAgg and upsampling schemes slow, while the time for the NonLocalAgg method is based on matching at half resolution. We also point out that CSBP attains a solution in a reasonable time, despite its MRF-based formulation. 1.0 GC [52] > 10 3 CSBP [28] 15 FastAgg [25] > 10 3 NonLocalAgg [27] 5.0 Yang et al. [31] > 10 2 Kopf et al [34] 88 EPC [15] 1.5 F-MRF (GC) [12] > 10 3 F-MRF (ICM) [12] > 10 2 Our upsampling 95 F-ECC, F-EMCC 2.0 Upsampling of [15] < 0.1 * F-ECC, F-EMCC 0.3 * * time needed for the 2MP images of our real data-set (see Fig.13)

B. Real Data
Dal Mutto et al. [18] provide real TOF-stereo data along with ground truth disparities, shown in Fig. 12. To be consistent with [18], we upsample the depth in a similar way, using a bilateral filter, which benefits from color segmentation. This procedure is used to initialize the depth in all fusion schemes. Table VI shows the mean square error (MSE) of the disparity estimation for several algorithms. The contribution of the stereo data in fusion is unquestionable in scenes A and B, where the depth varies locally. All of the fusion methods obtain a more accurate map than the initial one. However, scene C contains only planar objects, and the upsampling methods provide good results. Although the proposed scheme does not perform best in all examples, it always improves the initial estimate, which demonstrates the advantage of adaptive fusion (very similar numbers are obtained with F-EMCC). As expected, the use of stereo data only (e.g., [25]) performs well only with the textured scene B. Our upsamping filter is more accurate than [34] and less accurate than the filter proposed by [18]. Note that our implementation of [34] achieves better results than those reported in [18].
We also assess the performance of the fusion methods on the HCIbox data-set [17]. The scene shows the interior of a box that contains some objects (Fig. 11). Note that there is no texture, apart from some horizontal lines on the stairs and the ramp, hence stereo methods tend to fail. We follow the experimental setup of [17], thus evaluating the depth estimation based on some statistics of the absolute error, after excluding inter-reflection areas (see Fig. 11). We do not include the results of [17], since the authors provided a different inter-reflection mask with larger support area than the one used in [17]. 6 Table VII shows the error statistics of the algorithms. All of the fusion methods start from the same initial map, obtained by our upsampling method. The proposed fusion method achieves the lowest mean and median error (similar results are obtained with F-EMCC). The variance of F-MRF is increased (a local bias was observed owing to the global smoothness), while its median remains low. Because of the depth discontinuities, [34] yields a less accurate result compared to our upsampling.
We also captured our own challenging TOF-stereo data-set using a synchronized camera setup, developed in collaboration with 4D View Solutions 7 . Two HR (1624 × 1224) color cameras and a MESA Imaging SR4000 TOF (range) camera The three scenes (cropped) of the dataset used in [18].
(176 × 144) are mounted on a rail. The stereo baseline is 50cm approximately while the TOF camera is mounted in the middle. The algorithm of [19] allows us to transform the depth-sensor measurements into a very sparse stereo disparity map (1% sparsity), with an average error of 1.0 pixel. The sparse map is refined and upsampled as discussed in Sec. IV (see also Fig.3). For a fully real-time upsampling, we use our GPU-based implementation of the triangulation-based interpolation [15]. Note that our stereo rectification puts the principal points in the same position, rather than making the optical axes parallel; this maximizes the overlap between the images, given the relatively wide baseline. Our 'MIXCAM' data-set contains challenging cases, e.g., periodicities, weakly textured areas, thin objects, depth discontinuities, and so on. The fusion algorithm merges the stereo and depth data and the outcome is a dense HR disparity map. We reuse our upsampling filter in a post-processing step to fill missing disparities. A streak-based method fills any remaining gaps in all algorithms. The results are shown in Fig. 13. White areas denote unmatched pixels, while black areas mark the detected TOF-occlusions. The left column shows the left image, with the TOF image shown in the bottom-right corner at the true scale. Next columns show the results of ELAS, FastAgg, F-MRF (GC) and F-ECC algorithms; the last column show the disparity maps of F-ECC after post-processing.
ELAS fills local areas, surrounded by textured points, through an interpolation scheme. We intentionally show the results of ELAS before the streak-based filling; as opposed to FastAgg, where missing disparities after the left-right consistency check are filled. Clearly, a pure stereo algorithm cannot deal with large untextured areas, and the post-filling is unreliable. F-MRF provides fully dense results. Note that we run F-MRF with half-resolution images (812 × 612), due  to its tremendous memory requirements. Moreover, we set a fixed value along the disparity range for the data-term of all TOF-occlusion points, so that the global inference becomes independent of this area. F-MRF provides artifacts in stereo occlusions, that are next to TOF-occlusions when the scene contains large foreground objects. As with [13], the results of F-MRF scheme verify the lack of an adaptive fusion of the depth-and stereo-consistency data terms, as opposed to our methods. However, F-MRF seems to deal better with very thin objects (e.g. the branch of the plant), as already discussed above. Note that the biased range measurements of very slanted surfaces (e.g. the table-top) negatively affect the fusion schemes, in particular when the table surface lacks texture (e.g. first example). The proposed scheme provides very good results on average, especially after the post-processing step, which fills the gaps and refines the disparities. We obtain very similar results with the F-EMCC method, while EPC provides results visually close to ours, but with more gaps. The bilateral upsampling of [34] provides visually good results, but with blurred depth discontinuities (see also Fig. 3).
VII. CONCLUSIONS We have presented a high-resolution stereo matching algorithm that is guided by low-resolution depth data, thus helping the algorithm to compensate for its difficulty in estimating disparities over weakly textured areas. We cast the problem into a MAP formulation whose inference is obtained through a series of local optimization problems, solved hierarchically in a seed-growing manner. The latter characteristic yields an intrinsically efficient solution that allows for near realtime matching of 2.0MP images. The data-term of the energy function benefits from a correlation function that is capable of providing scores at subpixel disparities, from an adaptive cost aggregation step inside the window based on the depth data, and from an adaptive fusion of stereo-and depth-consistency terms based on the scene texture and the camera geometry. These properties lead to a more selective growing process that prevents the algorithm from propagating incorrect disparities. As a result, a low-complexity method builds an accurate highresolution disparity map. A quantitative comparison against pure stereo and stereo-depth fusion algorithms, as well as a qualitative assessment on real data, has validated the strong performance of the proposed method. Future research will include the optimum visiting order for seeds in the growing framework, as well as an adaptive window size, based on the local surface orientation.


extfieldsofstudy:
Computer Science
Mathematics

text:
Joint data imputation and mechanistic modelling for simulating heart-brain interactions in incomplete datasets

The use of mechanistic models in clinical studies is limited by the lack of multi-modal patients data representing different anatomical and physiological processes. For example, neuroimaging datasets do not provide a sufficient representation of heart features for the modeling of cardiovascular factors in brain disorders. To tackle this problem we introduce a probabilistic framework for joint cardiac data imputation and personalisation of cardiovascular mechanistic models, with application to brain studies with incomplete heart data. Our approach is based on a variational framework for the joint inference of an imputation model of cardiac information from the available features, along with a Gaussian Process emulator that can faithfully reproduce personalised cardiovascular dynamics. Experimental results on UK Biobank show that our model allows accurate imputation of missing cardiac features in datasets containing minimal heart information, e.g. systolic and diastolic blood pressures only, while jointly estimating the emulated parameters of the lumped model. This allows a novel exploration of the heart-brain joint relationship through simulation of realistic cardiac dynamics corresponding to different conditions of brain anatomy.

Introduction
Heart and brain are characterized by several common physiological and pathophysiological mechanisms [6]. The study of this multi-organ relationship is of great interest, in particular to better understand neurological diseases such as vascular dementia or Alzheimer's disease. The development of computational models simulating heart and brain dynamics is currently limited by the lack of databases containing information for both organs. Neuroimaging datasets often provide a limited number of cardiac-related measurements, usually restricted to brachial diastolic and systolic blood pressure (DBP, SBP) [7]. These quantities provide a limited assessment of the cardiac function, thus compromising the possibility of further analysis of cardiovascular factors in brain disorders.
On the contrary, the availability of a rich set of cardiac information in heart studies allows the use of cardiovascular models to estimate descriptors of the cardiac function that are not possible to measure in-vivo, such as contractility or stiffness of the heart fibers. These models optimize the parameters through dataassimilation procedures to reproduce the observed clinical measurements [11], but usually do not include neurological factors. The ability to jointly account for cardiovascular descriptors and brain information is key to gather novel insights about the relationship between heart dynamics and brain conditions.
Most of current studies relating heart and brain are based on statistical association models, such as based on multivariate regression [5]. While this kind of analysis allows to easily formulate and test association hypotheses, it usually offers a limited interpretation of the complex relationship between organs. This issue is generally addressed by mechanistic modeling, allowing deeper insights on physiological and biomechanical aspects. These models allow for example to describe the brain vasculature, and to quantify physiological aspects such as blood flow auto-regulation effects in the brain [1], up to the simulation of the whole-body circulation with detailed compartmental components [3]. Although these approaches offer a high level of interpretability, they are usually severely illposed and require large data samples and arrays of measurements to opportunely tune their parameters.
To bridge the gap between data-driven and mechanistic approaches to heartbrain analysis, in this work we propose to learn cardiovascular dynamics from brain imaging and clinical data by leveraging on large-scale datasets with missing cardiac information. This is achieved through an inference framework composed of two nested models accounting respectively for the imputation of missing cardiac information conditioned on the available cardiac and brain features, and for a Gaussian Process emulator that mimics the behavior of a lumped cardiovascular model. This setting allow us to formulate a probabilistic end-to-end generative model enabling imputation of missing measurements and estimation of cardiovascular parameters given a subset of observed heart and brain features.
Results on real data from the UK Biobank show that our framework can be used to reliably estimate and simulate cardiac function from datasets in which we have minimal cardiovascular information, such as SBP and DBP only. Moreover, the proposed framework allows novel exploration of the joint heart-brain relationship through the simulation of realistic cardiac dynamics corresponding to different scenarios of brain anatomy and damage.

Problem statement
We denote by ν the vector representing brain image-derived phenotypes (IDPs) and clinical information such as age or body surface area (BSA), and by x the vector of cardiac IDPs and blood pressure measurements. The vector x can be represented as x = {x, x obs } wherex represents the unobserved information we wish to impute and x obs the observed one. Moreover, we assume that for each observation x a corresponding set of parameters y of the associated lumped model is available. We would like to learn a generative model in which we assume that the unobserved measurements are generated by a latent random variable z conditioned on the variables ν. Hence, our generative process can be seen as sampling from a distribution p(z|ν) and then obtainingx with probability p(x|z, ν). Due to the association between cardiac IDPs, x, and lumped model parameters y, we also assume that the latter are dependent from ν through z. Our graphical model is shown in Figure 1, while the evidence lower bound (ELBO) of the joint data marginal p(y,x|x obs , ν) writes as: x y xobs Fig. 1. Graphical model of our framework. From a latent variable z we generate the unobserved featuresx conditioned on the variables ν and we estimate y via Gaussian process regression. During inference x obs and ν are used to estimate the approximated posterior q φ (z|x obs , ν).
The approximation of the posterior distribution by q φ (z|x obs , ν) defines the optimization of the ELBO through variational inference. The variational distributions q φ (z|x obs , ν) and p θ (x|ν, z) are parametrized respectively with parameters φ and θ. The term log p ω (y|x obs , z) in the ELBO denotes the emulator which approximates the mechanistic behavior of the lumped model via Gaussian Process regression parametrized by ω. The choice of a GP as emulator is motivated by the uncertainty of the data, hence it is desirable to obtain a distribution of interpolating functions rather than a single deterministic function. Moreover, GPs have already been proved to be valid emulators of 1D mechanistic vascular models [10]. The GP allow us to sample functions f (x) from a given prior parameterized by a mean µ(x) and covariance Σ(x), i.e: f (x) ∼ N (µ(x), Σ(x)) to obtain the marginal y ∼ N (µ(x), Σ(x) + σ 2 I). The prior mean µ(x) is here set to 0 and we use a radial-basis function (RBF) kernel for the covariance: Since our data is multi-dimensional, k j (x i , x ′ i ) is the kernel for the j th target and the i th predictor. The hyper-parameters of the kernel ω = {α, β} represent the output amplitude α and length scale β of the sampled functions. The goal during training is to learn the hyper-parameters that maximize the marginal likelihood of the observed data y .
The second term of the ELBO is related to the imputation ofx. The term log p θ (x|ν, z) denotes the log-likelihood of the imputed features, which can be seen as the reconstruction error. The last term KL(q φ (z|x obs , ν)||p θ (z|ν)) is the Kullback-Leibler divergence between variational approximation and prior for z, that can be expressed in a closed form given that both distributions are Gaussians. The imputation scheme is equivalent to a conditional variational autoencoder (CVAE) which has become a popular approach to feature imputation [8]. Fast and efficient optimization in our model is possible by means of stochastic gradient descent thanks to the closed form for data fit and KL terms, the use of the reparametrisation trick and Monte Carlo sampling.

Data processing and cardiovascular model
From UK Biobank we selected a subset of 3445 subjects for which T1, T2 FLAIR magnetic resonance images (MRI) and several brain and cardiac IDPs were available. Among the brain IDPs we used total grey matter (GM), total white matter (WM) and ventricles volumes. We used T1 and T2 FLAIR images to obtain the number of white matter hyper-intensities (WMHs) and their total volume relying on the lesion prediction algorithm (LPA), available from the lesion segmentation toolbox (LST) [14] of SPM 3 . WMHs are a common indicator of brain damage of presumably vascular origin [15]. We combined WM and GM volumes into a single measurement that we denoted as brain volume. The WMHs total volume and number of lesions presented a skewed distribution, and were Box-Cox transformed prior to the analysis. Regarding the cardiac IDPs we selected stroke volume (SV), ejection fraction (EF) and end-diastolic volume (EDV) for the left ventricle. All brain-related volumes were normalized by head size. Besides IDPs, we had access to blood pressure measurements (DBP, SBP) and sociodemographic features such as age and body surface area (BSA). We used DBP and SBP to compute MBP as M BP = DP + (SP − DP )/3. Next, we used the lumped cardiovascular model derived in [4] to obtain additional indicators of the cardiac function. In particular we estimated the contractility of the main systemic arteries (τ ), peripheral resistance (R p ) , the radius of the left ventricle (R 0 ), contractility of the cardiac fibers (σ 0 ) and their stiffness (C 1 ). The parameters of the model were selected based on the available clinical data. To obtain the target values for the emulator, the data-assimilation procedure was carried out according to the approach presented in [2,12].

Experiments
The data was split in two sets: one containing the full-information (2309 subjects), and one in which cardiac IDPs (x) and the estimated model parameters (y) were removed (1136 subjects, Table 1). The quality of imputation was compared to conventional methods such as mean, median and k-nearest neighbors (KNN) imputation. We subsequently assessed the relationship learnt by our model between brain features and cardiovascular parameters through simulation. Starting from the mean values of each parameter we sampled along the dimension of the different conditional variables ν that we used to parameterize the prior. This procedure allowed us to assess their influence in the inferred simulation parameters. Table 1. Variables used in our framework. Mean blood pressure (MBP), diastolic blood pressure (DBP), stroke volume (SV), end-diastolic volume (EDV), ejection fraction (EF), heart fibers contractility (σ0), ventricle size (R0), heart fibers stiffness (c1), peripheral resistance (Rp) and aortic compliance (τ )

Results
Data imputation and regression. We assessed the performance of our model by measuring the mean squared error (MSE) on the testing data, for both the imputation of missing cardiac information and the emulation of the lumped model parameters. The results in Figure 2a show that our method gives significantly better estimates than mean and median imputation, and comparable predictions to KNN, for which the optimal number of neighbors was optimized through 10fold cross-validation and corresponds to K = 10. At the same time the emulator consistently gives low errors for the parameters' estimation ( Figure 2b). In supp. Figure 3 we present the most relevant predictors for each emulated feature based on their β values, while a qualitative comparison between the distributions of imputed and emulated features compared to the ground truth data is available in in supp. Cardiovascular dynamics simulation and model plausibility. In Figure  3a we observe the change in the predicted parameters of the cardiovascular model as we sample along the range of values of the different conditional variables ν. Figure 3b shows the pressure-volume (PV) loops generated by the lumped model using the inferred parameters. The simulated PV loops highlight meaningful relationships: -An increase in the volume of WMHs is associated to decreased SV and EDV, together with a smaller reduction of ESV, leading to a decrease of EF which is related to reduced contractility. -Similar dynamics are associated to brain volume loss.
-The number of WMHs exhibits different dynamics than WMHs volume, associated to the increase in afterload and the increase of EDV.
From the plots showing the estimation of the parameters we can observe the ones driving the observed dynamics. For example, WMHs and brain volumes changes are mainly driven by the joint evolution of peripheral resistance (R p ), contractility (σ 0 ) and size of the left ventricle (R 0 ). The changes in the number of WMHs are related to heart-remodelling changes, driven by R 0 and by a decrease of σ 0 and the ventricular stiffness c 1 . In supp. Figure 2 we provide analysis and discussion of the remaining conditional features (age, BSA and brain ventricles volume). Overall, in the simulated dynamics we can identify several physiological responses in line with the clinical literature. Moreover, our model may be useful to give insights in currently controversial topics such as concerning the pathogenesis of WMHs. Our results suggest that the evolution of cardiac function with respect to brain and WMHs volumes is similar to the one due to aging (see supp. similar to the one related to ventricles enlargement. In both cases the changes in R 0 , σ 0 and C 1 suggest that the increase in the number of WMHs and the enlargement of the brain ventricles are related to heart-remodelling processes, such as loss of contractility or decrease in compliance. These findings are in line with clinical observations [9] relating lower cardiac output with higher burden of WMHs and reduced brain volume.

Conclusion
We presented a generative model that enables the analysis of complex physiological relationships between heart and brain in datasets where we have minimal available features. The framework allow us to emulate a lumped cardiovascular model through data-driven inference of mechanistic parameters, and provides us a generative model to explore hypothetical scenarios of heart and brain relationships. In the future, the model will allow to potentially transfer the knowledge learnt in UK Biobank to datasets where we have minimal cardiac information, to explore the relationship between brain conditions and cardiovascular factors in specific clinical contexts, such as in neurodegeneration. Our approach could also be extended to account for deep learning architectures, and the framework could be further improved by jointly accounting for multiple outputs, which are currently modelled independently, or by including spatial information from imaging data, beyond the modelling of scalar volumetric features. Furthermore, while the cardiovascular features considered in this study are rather general, more complex features will allow to study more realistic cardiovascular models. For example, while the mechanistic model used in this study does not simulate cerebral blood flow, previous studies suggested that WMHs may be due to local vascular impairment [13]. Hence, by selecting appropriate clinical features we could constrain the imputation by means of any biophysical model representing the desired aspect of systems biology. This could represent an innovative tool in real world scenarios, for which multi-modal patient data is often limited or not available.