Considering the reality of image formation as well as the statistical characteristics of contourlet coefficients, we propose a method for statistical image fusion in the contourlet domain. For high-frequency subbands, an image formation model was established that has a continuous-valued blur factor to reflect the actual imaging situation. According to this model, fused contourlet coefficients were estimated by use of an expectation-maximization (EM) algorithm. During the estimation, a contourlet hidden Markov tree model was adopted to grasp all dependencies among coefficients and aim for better estimation results. Because the blur factor is a continuous variable, we exploited an explicit expression to update the factor in the EM algorithm, which contributed to a decline in the number of iterations for convergence. Experimental results indicated that, especially for multifocus images, the proposed method provides more satisfying fusion results in terms of visual effects and objective evaluations than some typical fusion methods based on multiscale decomposition and some statistical fusion methods using a discrete blur factor.