text:
# The WARPS survey – IV: The X-ray luminosity–temperature relation of high redshift galaxy clusters

## 1 Introduction

Clusters of galaxies, as the largest virialised objects in the Universe, allow us a unique insight into the formation and evolution of mass clustering on cosmological time-scales. The X-ray emission from galaxy clusters originates from hot ($`\mathrm{T}10^{78}`$ K) intra-cluster gas which, in relaxed systems, is in hydrostatic equilibrium with the cluster’s total gravitational field. Since the temperature of the gas is proportional to its mass, the observed correlation between X-ray luminosities and temperatures implies a relationship between cluster baryon mass and total mass. Thus, studying the evolution of the L–T relation probes the interrelated evolution of these mass components. One of the first attempts at modelling these properties assumed self-similar evolution of both the gravitational potential and the baryonic intracluster medium (ICM) at varying cosmological epochs (?).

These self-similar predictions are in conflict with observational evidence. The strong positive evolution of the cluster X-ray luminosity function (XLF) originally predicted (assuming realistic cosmological parameters) is in stark contrast to the negative evolution reported by ? and ? for X-ray luminous clusters. The prediction is also in conflict with the results of more recent cluster surveys which find no significant evolution, at any luminosity, out to moderate redshift ($`z0.3`$, ?; ?), and again no significant (or mild negative) evolution at low and intermediate luminosities and high redshifts (?; ?; ?; ?). To account for such discrepancies, modifications were suggested to the self-similar theory. ? and ? modelled a ’pre-heated ICM’, in which an initial injection of energy broke the direct self-similar scaling by introducing an entropy floor in cluster cores. It also brought theoretical predictions of cluster properties into better agreement with observational constraints.

Tests of gravitational self-similar scaling are thus fundamental to the understanding of cluster evolution. One such test used in this analysis is the study of the cluster X-ray luminosity–temperature (L–T) relation, which was predicted by the original self-similar model to be described by an $`\mathrm{L}\mathrm{T}^2`$ law. However, a relationship closer to $`\mathrm{L}\mathrm{T}^3`$, has been found by several cluster surveys and compilations (e.g. ?; ?; ?), whilst even further steepening of the relation is required to agree with galaxy group measurements (?; ?). In addition to predicting the slope of the L–T relation, self-similar theory also predicts significant evolution in the normalisation of this relation at different redshifts. However, the evolution predicted is reduced in low density cosmologies, or by significant entropy injection at early epochs (e.g. ?; ?; ?).

Attempts at constraining evolution in the L–T relation, at redshifts above $`z0.1`$, are complicated by the inherent difficulty of temperature estimation from often poor signal-to-noise data. ?, using $`ASCA`$ data from a sample of 38 clusters at $`z>0.14`$, found no significant evolution out to $`z0.3`$. Prior to this work, the only attempt at quantifying evolution of the L–T relation to $`z0.3`$ was that of ?, whose results were also consistent with no evolution. ? used low signal-to-noise spectra of 67 clusters observed with the $`Einstein`$ observatory, each individually spectrally fitted, and then averaged together in redshift bins. Recently ? have used $`ASCA`$ measurements of clusters drawn from the $`Einstein`$ Medium Sensitivity Survey (EMSS, ?; ?) to extend the conclusion of no evolution to $`z>0.5`$.

In this paper we present an L–T relation based on $`ROSAT`$ spectra of WARPS clusters, in various redshift and luminosity bins. Although this procedure is analogous to the one adopted by ?, a more sophisticated simultaneous fitting method is adopted here. We discuss the data analysis techniques used to obtain temperature estimates from often low signal-to-noise $`ROSAT`$ PSPC data, and present results from simulated spectra to demonstrate their validity. In order to allow self-consistent discussion of any evolutionary effects in our relation, we compare our results with a local L–T relation over a wide range of temperatures and luminosities.

## 2 Sample selection

Traditional studies of galaxy clusters have relied on optically selected cluster samples. Whilst these surveys were often large, they suffered from biasing caused by projection effects and erroneous identifications (?; ?; ?). The recent trend, therefore, has been to base large cluster surveys on X-ray selection techniques. This takes advantage of the X-ray emission from diffuse gas, trapped in a cluster’s gravitational potential well, representing direct evidence of a three-dimensionally bound system. Several recent serendipitous cluster surveys have been compiled, all based on archival $`ROSAT`$ PSPC data (e.g. ?; ?; ?; ?; ?).

The Wide Angle $`ROSAT`$ Pointed Survey (WARPS, ?) is also based on archival $`ROSAT`$ PSPC data. It is a serendipitous cluster survey which uses a sophisticated source detection algorithm, VTP (?), capable of detecting arbitrarily shaped emission, down to very low surface brightness levels. An extensive optical follow-up program then spectroscopically confirms the concordance of galactic redshifts within each candidate cluster (?). When complete, the WARPS survey will have generated a statistically complete, flux-limited catalogue covering a solid angle of more than 73 $`\mathrm{deg}^2`$ and including approximately 150 X-ray selected galaxy clusters.

The sample used in this paper comprises 91 clusters identified in the WARPS survey, with redshifts in the range from $`z=0.11`$ to $`z=0.833`$. The WARPS survey is flux-limited at $`6\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ (0.5–2 keV) in total cluster flux and is restricted to $`ROSAT`$ pointings of more than 8-ks exposure at $`|\mathrm{b}|>20\mathrm{deg}`$. Only the most sensitive area of the PSPC detector, at off-axis angles of $`3<r<15`$ arcmin, is used. In addition to the archival $`ROSAT`$ PSPC data, deep optical imaging of the clusters, optical spectroscopy of all the brightest cluster galaxies (BCGs) and in some cases $`ROSAT`$ HRI data are available (?; ?).

## 3 Data reduction

### 3.1 Spectral fitting

For this work, the archival $`ROSAT`$ PSPC data, containing the 91 WARPS clusters selected, were initially reduced and analysed using the Starlink ASTERIX package (?). The data were screened to minimize particle background rates by restraining the Master Veto Rate housekeeping parameter to a range of 0 to 170. The Aspect Error and Accepted Event Rate parameters were also limited to ranges of 0 to 2 and 0 to 30 respectively. The observation time thus removed was on average 5 per cent. The resulting cleaned data were then used to produce binned spectral cubes, with one spectral and two spatial dimensions, from $`ROSAT`$ channels 11 through 231, with energy ranges of 0.114 to 2.295 keV

Background subtraction for each field was achieved by creating a model background spectral cube. Firstly an annulus was selected between 9 and 15 arcmin off-axis. This area is included in the WARPS analysis and was selected to minimize vignetting extrapolation errors in the background model. The background cube was projected to an image and searched for point sources using the ASTERIX PSS algorithm. Removal of the sources resulted in a new background model. This process was iteratively repeated until no more point sources were detected. The point source free cube was then further edited manually to eliminate any extended sources, such as cluster emission. The background level of this cube was then extrapolated over all off-axis angles, using the known energy-dependent vignetting behaviour of the detector, to produce a final background model.

After background subtraction, a spectrum of the cluster was obtained by projecting a selected area of the spectral cube on to the energy axis, before correcting for energy-dependent vignetting and exposure effects. The circular extraction regions around the clusters had initial radii of 700 and 900 kpc, for $`z<0.5`$ and $`z>0.5`$ respectively (we assumed $`\mathrm{H}_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`q_0=0`$). These radii were manually re-selected, where required, to eliminate contamination from nearby non-cluster sources. A background spectrum from the same area of sky, obtained from the source-free background model, was then associated with the cluster spectrum to allow spectral fitting using the maximum-likelihood C-statistic (?). The number of counts in the cluster spectra ranged from 25 to 1001 (0.11–2.3 keV) with a median value of $`120`$ counts. The signal-to-noise ratios of the data varied dramatically, depending primarily on background level, exposure time and extraction radius. The average signal-to-noise level was $`6.8`$ in the 0.1–2.3 keV band, but was significantly higher in the 0.5–2 keV band used in the initial detection (?). Thus C-statistic fitting was selected to avoid biases inherent in $`\chi ^2`$ fitting for low count-per-bin data (?).

The spectral models fitted to the data were the ? and MEKAL (?; ?) codes for an optically thin thermal plasma. The column density was fixed at the Galactic value determined from radio measurements (?), and a metallicity of 0.3 $`Z_{}`$ was assumed. It was found that both codes gave very good agreement; the results presented here are from the MEKAL model.

Due to the very low number of counts in many of the individual cluster spectra, the spectral fits were often poorly constrained, with several not converging. In addition fits to a number of high redshift, low signal-to-noise spectra converged to low temperatures with small errors (e.g. T=$`1.3\pm 0.1`$ keV). This was of concern as, according to the observed, low redshift, L–T relation, this implied luminosities too low for these clusters to be detected above the flux-limit of the survey.

There are two primary characteristics of the data which dominate the accuracy and reliability of the cluster temperature measurements, and cause the effect discussed in the previous paragraph. Firstly, the signal-to-noise of the cluster above the background, and secondly, the accuracy to which the background level can be determined. To investigate possible systematic biases caused by poor photon statistics and erroneous background estimates, extensive simulations were performed, using the ASTERIX SSIM package. Fig. 1 shows spectrally fitted temperatures from 100 datasets all simulated from a 5keV MEKAL model. Each cluster was simulated as excess emission on top of a PSPC background of 800 counts (0.1–2.3 keV, in the extraction region). This background is brighter than any found in the actual observations; hence the background level can be accurately determined in these simulations. The total number of counts is varied according to Poisson statistics. Spectral fits were performed on the emission after the background was subtracted. It can be seen from the histogram that for the spectra (in this case with an average of $`190`$ counts in the 0.1–2.3 keV band, giving a signal-to-noise of 6.0), three distinct temperature populations were found. Spectral fits that were poorly constrained had the best-fitting temperatures pegged at a maximum value (17 keV in this case). For a different subset of spectra the fit essentially recovered the simulated temperature. However, a third population of low temperature fits with small errors was also found, as observed with some of the $`ROSAT`$ datasets. This is a demonstration of the bias introduced with low signal-to-noise in the cluster spectra, but with the background level reasonably well determined.

The combination of the biases due to low signal-to-noise and poor background statistics, is demonstrated in Fig. 2, which shows the results from 200 simulated spectra, each containing $`120`$ counts (0.1–2.3 keV). The background level in this simulation was tailored to give a signal-to-noise values around 7 (i.e. to have $`175`$ counts, 0.1–2.3 keV). Here the individual fits cannot recover the simulated temperature of 5 keV.

A method of increasing the signal-to-noise of the spectra was thus required to enable accurate temperature estimates to be made. This was achieved via simultaneous multiple fitting of data sets in binned groups.

### 3.2 Multiple fitting

The effect of multiply fitting poor signal-to-noise datasets can be seen in the inset panels in both Fig. 1 and Fig. 2. These show temperatures estimated from simultaneously fitting subsets of five of the simulated datasets discussed in the previous section. In most cases, values around the simulated temperature of 5 keV were found by these fits.

As can be seen from Fig. 2, the primary effect of reducing the cluster count level, together with incorporating a low statistical accuracy background measurement, is to produce a number of simultaneous fits that are not constrained and yield the maximum allowed temperature of 17 keV. As the simulated signal-to-noise level is similar to the poorer PSPC spectra used in this work, we might expect to see this effect when using this method on WARPS clusters. No such trend was observed though. This was probably because in the real measurements a large background area was sampled in order to derive an accurate background level. This is in contrast to the simulations, where a much noisier background measurement resulted from using an area $`100`$ times smaller (as required by the SSIM program). This difference also accounts for the increased accuracy of the simulations with higher backgrounds, but similar signal-to-noise. In the analysis of PSPC data, the large area used in the estimation of the individual background level for each cluster serves to minimize the second of the two biases discussed in the previous section.

Numerous simulations were carried out to test the effectiveness of the method. Gaussian distributions in both temperature (e.g. $`\pm 1\mathrm{k}\mathrm{e}\mathrm{V}`$ at 90 per cent confidence) and metallicity ($`\pm 0.1Z_{}`$), around the simulated values, were included. These extended simulations also recovered a correct temperature distribution, though with a larger spread spanning around $`\pm 2\mathrm{k}\mathrm{e}\mathrm{V}`$. The temperature errors obtained from the fits, however, generally encompassed the distribution range.

From all the simulations, we conclude that multiple fits of five spectra, with each individual spectra containing 100–200 counts (i.e. a signal-to-noise ratio of 6–7 for PSPC spectra), with one free parameter, and the background statistically well determined, produce temperatures in which the level of any systematic bias is smaller than the random error. There may be an additional temperature uncertainty at high temperatures due to PSPC calibration uncertainties. ? found evidence of significant differences between PSPC and $`ASCA`$ temperatures for T $`5`$ keV and higher, although the size of the differences are around the same size as our statistical error.

Multiple fitting then, was the method adopted for analysis of the PSPC spectra. In order to produce averaged spectra, we binned the clusters into redshift and luminosity bins. Where possible at least five clusters were included in each bin. Simultaneous multiple fitting was done within Version 10.0 of the XSPEC fitting package (?) in the XANADU X-ray software suite. This enables several sets of data to be fitted to the same model, with certain parameters either frozen or tied to be the same for all datasets. In our case, the clusters were fitted to the same temperature in the plasma code, but with individually varying normalisations. Again, column densities were fixed at Galactic values for each cluster and metallicities were fixed at 0.3 $`Z_{}`$.

The main complication arose from the need to use maximum likelihood statistics. To correctly fit using the Cash statistic, XSPEC requires a non-background subtracted cluster spectrum, with a pre-fitted background component. Thus multiple fitting first involved individually fitting the relevant background for each cluster. The model used was an absorbed combination of a MEKAL plasma, a power law and a Gaussian emission line profile, with an initial energy of 0.54 keV (as used by ? to model the X-ray background spectrum).

### 3.3 Luminosities

Flux estimates in the 0.5–2 keV band were obtained for each cluster from the individual fits to the spectra. These were then k-corrected, and corrected for low surface brightness emission outside the selected area assuming a King profile with $`\beta =\frac{2}{3}`$ and cluster core radius estimated from the VTP data, as described by ?. The fluxes were then converted to luminosities, using values of $`\mathrm{H}_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\mathrm{q}_0=0`$. These luminosity estimates were then used for the binning of the clusters within each redshift range.

Due to the need for an accurate bolometric conversion from 0.5–2 keV luminosities, which relies on temperature estimates, the bolometric luminosities of each cluster used in the final results were obtained after the multiple fit. Each cluster was bolometrically corrected using the characteristic temperature of the bin, before again being k-corrected and convolved through a King profile correction to produce accurate luminosities. The typical bolometric correction resulted in a factor of 2–3 increase from the 0.5–2 keV luminosity. The cluster luminosities within the bin were then averaged to provide a characteristic bin luminosity.

The steep slope and large inherent scatter of the low redshift L–T relation (e.g. ?; ?), heightens the importance of restricting the range of luminosities within each luminosity bin. A large spread in luminosities could result in large temperature differences in the clusters contained within any bin. Which would then invalidate the simultaneous fitting method employed. The flux-limited nature of the WARPS data meant that a large spread in luminosities was found only when attempting to bin together a number of the more luminous clusters within each redshift range. In two cases, for instance, there was no other cluster of similar luminosity to the most luminous cluster, and so this cluster was individually fitted for temperature. Both WARPJ1416.4+2315 at $`z=0.137`$ and WARPJ1418.5+2510 at $`z=0.294`$ had more than 800 counts in the 0.1–2.3 keV band, allowing a robust temperature fit. This also had the effect of reducing the largest spread of signal-to-noise ratios within any bin to less than a factor of two.

## 4 Results

### 4.1 X-ray luminosity–temperature relations

Fig. 3 shows the L–T relation from binned PSPC spectra of WARPS clusters. Different symbols represent the redshift bins used. The luminosity errors are negligible in most cases, and primarily due to temperature uncertainties affecting the bolometric correction.

To put this L–T relation in context, Fig. 4 through Fig. 6 show the WARPS data points superposed on cluster L–T measurements drawn from the literature. Data are included from a number of different sources, including ?, ? and recent results from ?. The ? dataset includes cluster measurements from a number of other works (?; ?; ?; ?). Data from a recent survey of 24 galaxy groups showing diffuse X-ray emission (?), are also included.

To enable us to discuss any possible evolution within the WARPS L–T relation, it is important to have a suitable low redshift relation to compare to. ? have found no evolution in the L–T relation out to $`z0.3`$. So we include all clusters at $`z<0.2`$ in our low redshift sample (Fig. 4). The WARPS points are in good agreement with previous measurements and help fill the gap in the L–T observations at T=1–3 keV.

In order to define a fiducial low-redshift L–T relation, we obtained the best-fitting power law to the data in Fig. 4 at $`\mathrm{T}>1\mathrm{keV}`$. Below T=1 keV, there may be a systematic trend towards lower luminosities, away from the best-fitting line obtained at higher temperatures (?; ?). To allow for the systematic error intrinsic to the use of datasets from three different X-ray observatories and measured over different energy ranges, a minimum error of 5 per cent was assigned to each cluster temperature where the existing error was not already greater than this value. The best fitting slope was then found using $`\chi ^2`$ minimisation on the temperature axis. For $`\mathrm{L}_{\mathrm{Bol}}=\mathrm{CT}^\alpha `$, we find $`\alpha =3.15\pm 0.06`$ and $`\mathrm{C}=6.04\pm 1.47\times 10^{42}\mathrm{ergs}^1`$ (90 per cent confidence regions). Most previous cluster L–T relations had a best-fitting slope of around $`\mathrm{L}\mathrm{T}^{2.8\pm 0.2}`$ (e.g. ?, ?). Our result has a slightly steeper slope primarily due to the inclusion of less luminous clusters not analysed in previous investigations.

### 4.2 Measuring evolution

For the purposes of constraining evolution in the normalisation of the WARPS L–T relation we assume a fixed slope of $`\alpha =3.15`$ (section 4.1). We then further parametrize the L–T relation as $`\mathrm{L}_{\mathrm{Bol}}=\mathrm{C}(1+z)^\eta \mathrm{T}^\alpha `$ where $`\mathrm{C}=6.04\times 10^{42}\mathrm{ergs}^1`$, the best-fitting normalisation of the local L–T relation. The $`\eta `$ term then parametrizes evolution at any redshift $`z`$.

We now consider the WARPS measurements in the redshift ranges shown in Fig. 5 and Fig. 6. Assuming the best-fitting slope from the low redshift clusters, we again fit for normalisation by minimising the $`\chi ^2`$-statistic. Again only measurements above 1 keV were included in the fit. Taking the central redshift of these bins allows us to evaluate $`\eta `$ and its associated error. Most theoretical studies discuss evolution in the L–T relation in terms of an $`\mathrm{\Omega }_0=1`$ cosmology. Table 1 displays our evolutionary constraints ($`\eta `$ values) found for the three $`z>0.2`$ redshift ranges assuming $`\mathrm{\Omega }_0=1`$. Wherever a different cosmology is specified, we corrected cluster luminosities accordingly (we assume $`\mathrm{H}_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ throughout).

In order to further constrain evolution in the L–T plot, the number of data points in each redshift bin must be increased. Using datasets from ? and Donahue et al. (1999) in conjunction with the WARPS bins, we were able to place tighter constraints on evolution in the L–T plot. The data used are represented in Fig.s 5 and 6, and the resulting $`\eta `$ values are listed in Table 1.

It can be seen from Table 1 that in general the WARPS measurements are consistent with no evolution of the normalisation of the L–T relation (though one redshift range provides marginal evidence for evolution). When the WARPS points are combined with other existing cluster measurements the results are also consistent with no evolution. To quantify the level of evolution measured we give the averaged $`\eta `$ value of the three separate redshift range measurements. The averages have been weighted by the inverse errors. The averaged $`\eta `$ value of the combined dataset for $`\mathrm{\Omega }_0=1`$ ($`\eta =0.19\pm 0.38`$), is consistent with no evolution to $`z0.8`$. Changing the cosmology used in the fit to an open $`\mathrm{\Omega }_0=0.3`$, $`\lambda _0=0`$ Universe modifies the conclusion only slightly, producing $`\eta =0.60\pm 0.38`$, still consistent with no evolution at around the $`1.5\sigma `$ level.

## 5 Discussion

### 5.1 Theoretical relations

The implications of the evolution parameter, $`\eta `$, can now be discussed in terms of evolution of X-ray temperatures and luminosities in the L–T relation. Here we focus primarily on comparisons with the entropy driven evolution model of ? and its extension to various cosmologies (?). In this model the evolution of the X-ray properties of clusters is split into contributions from the gravitational heating of the ICM and from purely gas related processes, such as shock heating or radiative cooling. Assuming a minimum entropy (or entropy floor) in cluster cores allows deviations from the pure self-similar scaling of the cluster ICM with the underlying dark matter distribution (e.g. ?; ?). The entropy driven model further allows the exploration of the regime between the purely self-similar scaling predictions of ? and the effect of a pre-heated ICM.

As discussed by ?, the heating and cooling of the ICM gas can be described in terms of evolution of the core entropy, if we assume that the slope of the primordial density fluctuation power spectrum, $`n`$, controls the gravitational growth of clusters. Parameterizing this evolution in the form $`s_{\mathrm{min}}=s_{\mathrm{min}}(z=0)+c_vϵ\mathrm{ln}(1+z)`$ describes the level of the entropy floor, $`s_{\mathrm{min}}`$, at any redshift, $`z`$. The term $`ϵ`$ thus describes how this entropy evolves ($`c_v`$ is just the specific heat capacity of the gas at constant volume). A negative value represents continual shock heating of the ICM during cluster collapse and $`ϵ>0`$ describes a cluster core dominated by radiative cooling, whilst $`ϵ=0`$ is the constant entropy model (e.g. ?). It is then possible to derive a functional form for the L–T relation, with evolution of its normalisation parametrized in terms of $`n`$, and $`ϵ`$. If we assume the density profile parameter is equal to the canonical value of, $`\beta =2/3`$, then both the ? and ? scaling relation models give an L–T relation of the form $`\mathrm{L}_{\mathrm{Bol}}(1+z)^{3ϵ/4}\mathrm{T}^{2.75}`$. Although we find a slightly different value of the slope of the L–T relation than this, due to the assumption of weak self-similarity (section 2.4, ?), we are only concerned with the evolution of normalisation of the relation. We therefore use a fixed slope of $`\mathrm{L}_{\mathrm{Bol}}\mathrm{T}^{3.15}`$ when dealing with these scaling relations. The main effect of the difference in slope is to introduce a slight degeneracy between $`n`$ and $`ϵ`$ when calculating the expected normalisation of the L–T relation (e.g. fig. 2, ?).

We can use the scaling relations to predict the evolution in normalisation of the L–T relation in terms of $`n`$ and $`ϵ`$. An example is shown in Fig. 7. Here we have assumed a temperature scaling relation as in Equation 2.11 of ?. Assuming an underlying cosmology ($`\mathrm{\Omega }_0=0.3`$, $`\lambda _0=0`$ in this case) and a value of the initial power spectrum (here $`n=1.6`$, in agreement with the measured APM galaxy power spectrum, ?), we can use the lack of observed evolution in the L–T relation to constrain evolution in the cluster core entropy. This evolution of the entropy floor is represented by the expected normalisation of the L–T relation changing, with redshift, as a function of the value of $`ϵ`$, as shown in Fig. 7. Here the degeneracy with $`n`$ is introduced in the calculation of Equation 2.11 of ? and also with the effect of a changing epoch of cluster formation as in Equations 2.12 and 2.13. The lack of observed evolution (as demonstrated by the value of the $`\eta `$ parameter being consistent with zero), implies almost no evolution of cluster core entropy to $`z0.8`$ given the assumed cosmology and value of $`n`$. Table 2 details our calculated $`ϵ`$ values in both an open and closed cosmology, for three different $`n`$ values. Also tabulated are the corresponding values required for pure self-similar scaling in an $`\mathrm{\Omega }_0=1`$ cosmology (?), given by $`ϵ_{\mathrm{ss}}=\frac{(n+7)}{(n+3)}`$, which is strongly ruled out.

So far we have not attempted to place any cosmological constraints. To do this we need to break the degeneracy that allows evolution in the X-ray properties to be explained via either gravitational or gas-phase phenomena. This requires the introduction of a second dataset. The XLF, for instance, allows the values of $`n`$ and $`ϵ`$ to be constrained in an almost orthogonal manner to the L–T evolution constraints. This is demonstrated in fig. 3 of ?. If the XLF is non-evolving to high redshift, as several recent studies have found (?; ?; ?; ?), then important constraints can be placed on the underlying cosmology. Given that the value of $`ϵ`$ is more tightly constrained in this work than in ?, any fluctuation spectrum of $`n<1.5`$ requires a Universe with a low value of $`\mathrm{\Omega }_0`$ ($`\mathrm{\Omega }_0<0.5`$). Correspondingly, values of $`n1`$ are required to be consistent with the observed XLF and L–T relation evolutionary constraints in an $`\mathrm{\Omega }_0=1`$ universe. Allowing for the possibility of some negative evolution in the XLF at the highest luminosities, as some surveys find (?, ?), may relax these constraints slightly, although the evolution at these high luminosities is not observationally well determined.

### 5.2 Cooling flows

The abundance of cooling flows within the WARPS cluster survey could have a large effect on the estimated cluster luminosity–temperature relation. ? showed that when the most X-ray luminous clusters were modelled isothermally, non-cooling flow clusters had significantly higher temperatures than cooling flow clusters. Furthermore, they found a significant change in the best-fitting slope of the cluster L–T relation when cooling flow clusters were analysed with an isothermal model, as opposed to a model with a central cooling flow component. They claimed a change in slope from $`\mathrm{L}\mathrm{T}^3`$ to approximately $`\mathrm{L}\mathrm{T}^2`$ when cooling flows were taken into account.

Due to the low signal-to-noise of the $`ROSAT`$ X-ray data used in this work, fitting the data for a cooling flow component is not possible. However, we were able to estimate the fraction of clusters in our sample that exhibit evidence for a cooling flow. It has been shown that most clusters ($`70`$ per cent) have or can form cooling flows (e.g. ?; ?), though the cooling flow size is proportional to the X-ray luminosity. Some cooling flows also exhibit other characteristics, most notably optical line emission from the BCGs. The fraction of clusters in X-ray flux limited surveys displaying these properties is around 33 per cent (e.g. ?; ?), though this fraction is highly dependent on the quality of the optical spectra. Less than this level of cooling flow activity in the WARPS sample would be a good indication that our survey does not over-sample cooling flows.

The optical spectra obtained to measure galaxy redshifts allowed analysis of line emission in the BCG of each cluster. In a preliminary analysis, H$`\alpha `$ or \[OII\]$`\lambda 3727`$ emission lines were detected in 4 of the 26 BCG spectra analysed. This implies a fraction of 16 per cent cooling flow clusters in our survey. This might be expected if the frequency of cooling flows exhibiting optical emission lines increased in more X-ray luminous systems.

Cooling flows would produce a lower temperature estimate for the same X-ray luminosity, thus they would, in effect, mimic evolution in the WARPS results in the direction predicted for $`\mathrm{\Omega }_0=1`$, assuming no core entropy evolution (?). Removal of the effects of a fraction of the sample containing cooling flows would hence favour lower values of $`\mathrm{\Omega }_0`$, strengthening our conclusion.

## 6 Conclusions

We have presented a cluster X-ray luminosity–temperature relation out to high redshifts ($`z0.8`$), based on $`ROSAT`$ PSPC spectra of WARPS clusters. Due to the low signal-to-noise of many of our $`ROSAT`$ PSPC spectra, a multiple fitting technique was adopted in order to constrain temperatures. The 91 clusters analysed were binned into redshift and luminosity bins, and spectrally fit to a characteristic temperature.

In order to analyse the evolution of the L–T relation, a large sample of low redshift ($`z<0.2`$) clusters was compiled from the literature. The best-fitting power law to the data at low redshift and at T$`>`$1 keV was found to be a power law of index $`\mathrm{L}\mathrm{T}^{3.15\pm 0.06}`$ (90 per cent confidence). This relation was then compared with the high redshift measurements to constrain possible evolution in the normalisation of the L–T relation.

The binned $`ROSAT`$ measurements of WARPS clusters were consistent with no evolution of the L–T relation to $`z0.8`$. When combined with previous cluster measurements from the literature they agree with the no evolution result of ?, but give stronger constraints; $`\eta =0.19\pm 0.38`$ for $`\mathrm{\Omega }_0=1`$, or $`\eta =0.60\pm 0.38`$ for $`\mathrm{\Omega }_0=0.3`$ where $`\mathrm{L}_{\mathrm{Bol}}(1+z)^\eta \mathrm{T}^{3.15}`$. Using the entropy driven evolution model of ? and ?, a limit to the implied entropy evolution of cluster cores was found. Whilst the level of entropy evolution is slightly dependent on the value of the primordial density fluctuation spectrum, $`n`$, low amounts of evolution were favoured. The pure self-similar scaling expectations of ? were found to be strongly excluded.

The lack of evolution in the WARPS L–T relation, in conjunction with other recent results, such as a non-evolving XLF, implies either an open cosmology ($`\mathrm{\Omega }_0<1`$) or one closed by a cosmological constant, assuming $`n<1`$. The forthcoming X-ray missions Chandra and XMM will allow more accurate temperature estimates of clusters, out to high redshifts, which will greatly improve our understanding of the L–T relation and its cosmological implications.

## 7 Acknowledgments

The authors would like to thank Megan Donahue and collaborators for advance access to their cluster measurements. Thanks to Scott Kay for help with the entropy evolution model. We thank the Birmingham galaxies and cluster group for useful discussions and input. The data used in this work have been obtained from the Leicester database and archive service (LEDAS). Computing facilities provided by the STARLINK project have been used in this work. BWF acknowledges the receipt of a PPARC studentship and LRJ also acknowledges the support of PPARC. HE acknowledges financial support from SAO contract SV4-64008 and NASA LTSA grant NAG 5-8253.


format: text

text:
# Field and intensity correlation in random media

## A

We give the complete expression for the long range contribution $`C_2`$ including absorption and internal reflection, for quasi-one dimensional geometry. We follow the approach of  in which the correlation function was written in the form:

$$C_2(\mathrm{\Delta }\nu )=2\frac{N_a+N_b+N_\alpha }{gBL},$$
(A1)
with an overall factor of 2 to account for the one channel-in, one channel-out experimental setup. We obtain the following expressions, which differ for $`B`$ and $`N_\alpha `$ from :

$`B`$ $`=`$ $`\left[1+2z_b^2(\gamma _+^2\gamma _{}^2)+z_b^4(\gamma _+^2+\gamma _{}^2)^2\right]\left[\mathrm{cosh}(2\gamma _+L)\mathrm{cos}(2\gamma _{}L)\right]`$ (A3)
$`+`$ $`4z_b^2(\gamma _+^2+\gamma _{}^2)\left[\mathrm{cosh}(2\gamma _+L)+\mathrm{cos}(2\gamma _{}L)\right]`$ (A4)
$`+`$ $`\left\{4\gamma _+z_b\left[1+z_b^2(\gamma _+^2+\gamma _{}^2)\right]\mathrm{sinh}(2\gamma _+L)\right\}+\left\{4\gamma _{}z_b\left[1z_b^2(\gamma _+^2+\gamma _{}^2)\right]\mathrm{sin}(2\gamma _{}L)\right\},`$ (A5)
$`N_a`$ $`=`$ $`{\displaystyle \frac{1}{2\gamma _+(\gamma _+^2\alpha ^2)}}`$ (A6)
$`\times `$ $`\{[(2\gamma _+^2\alpha ^2)+z_b^2(8\gamma _+^4+2\gamma _+^2\gamma _{}^2\gamma _{}^2\alpha ^23\gamma _+^2\alpha ^2+\alpha ^4)+z_b^4(\gamma _+^2+\gamma _{}^2)(2\gamma _+^42\gamma _+^2\alpha ^2+\alpha ^4)]`$ (A8)
$`\mathrm{sinh}(2\gamma _+L)`$
$`+`$ $`[4z_b\gamma _+(3\gamma _+^2\alpha ^2)+4z_b^3\gamma _+(3\gamma _+^4+\gamma _+^2\gamma _{}^22\gamma _+^2\alpha ^2+\alpha ^4)]\mathrm{sinh}^2(\gamma _+L)\},`$ (A9)
$`N_b`$ $`=`$ $`{\displaystyle \frac{1}{2\gamma _{}(\gamma _{}^2+\alpha ^2)}}`$ (A10)
$`\times `$ $`\{[(2\gamma _{}^2+\alpha ^2)+z_b^2(8\gamma _{}^4+2\gamma _+^2\gamma _{}^2+\gamma _+^2\alpha ^2+3\gamma _{}^2\alpha ^2+\alpha ^4)z_b^4(\gamma _+^2+\gamma _{}^2)(2\gamma _{}^4+2\gamma _{}^2\alpha ^2+\alpha ^4)]`$ (A12)
$`\mathrm{sin}(2\gamma _{}L)`$
$`+`$ $`[4z_b\gamma _{}(3\gamma _{}^2+\alpha ^2)4z_b^3\gamma _{}(3\gamma _{}^4+\gamma _+^2\gamma _{}^2+2\gamma _{}^2\alpha ^2+\alpha ^4)]\mathrm{sin}^2(\gamma _{}L)\},`$ (A13)
$`N_\alpha `$ $`=`$ $`{\displaystyle \frac{\gamma _+^2+\gamma _{}^2}{2\alpha (\gamma _+^2\alpha ^2)(\gamma _{}^2+\alpha ^2)}}`$ (A14)
$`\times `$ $`\{[\alpha ^2z_b^2(6\gamma _+^2\gamma _{}^2+5\gamma _+^2\alpha ^25\gamma _{}^2\alpha ^2+\alpha ^4)+z_b^4\alpha ^2(2\gamma _+^2\gamma _{}^2\gamma _+^2\alpha ^2+\gamma _{}^2\alpha ^2)]\mathrm{sinh}(2\alpha L)`$ (A15)
$`+`$ $`[4z_b\alpha ^1(\gamma _+^2\gamma _{}^2\gamma _+^2\alpha ^2+\gamma _{}^2\alpha ^2\alpha ^4)+4z_b^3\alpha (3\gamma _+^2\gamma _{}^22\gamma _+^2\alpha ^2+2\gamma _{}^2\alpha ^2)]\mathrm{sinh}^2(\alpha L)\}.`$ (A16)
## B

We show that the time of flight distribution and the field correlation function with frequency shift are Fourier transform pairs. The response to a wave packet $`f(t)`$ is $`E(t)=𝑑\nu \mathrm{exp}[i2\pi \nu t]F(\nu )E(\nu )`$, where $`F(\nu )`$ is the Fourier transform of $`f(t)`$ and $`E(\nu )`$ is the field measured at a given point at frequency $`\nu `$. From this expression, the ensemble average of the temporal intensity variation of the transmitted pulse may be written as,

$$E^2(t)=𝑑\nu _1𝑑\nu _2\mathrm{exp}[i2\pi (\nu _1\nu _2)t]E(\nu _1)E^{}(\nu _2)F(\nu _1)F^{}(\nu _2).$$
(B1)
If $`f(t)`$ is a short pulse with a bandwidth much wider than the correlation frequency $`\delta \nu `$, $`F(\nu _1)`$ and $`F(\nu _2)`$ are nearly identical as long as $`\nu _2\nu _1`$ is much smaller than the pulse bandwidth. With the change of variables $`\nu _1=\nu +\mathrm{\Delta }\nu /2`$ and $`\nu _2=\nu \mathrm{\Delta }\nu /2`$, $`F(\nu _1)F^{}(\nu _2)|F(\nu )|^2`$ and Eq. (B1) becomes

$$E^2(t)=𝑑\nu 𝑑\mathrm{\Delta }\nu \mathrm{exp}[i2\pi \mathrm{\Delta }\nu t]E(\nu +\mathrm{\Delta }\nu /2)E^{}(\nu \mathrm{\Delta }\nu /2)|F(\nu )|^2.$$
(B2)
Within a frequency range in which the variation in the averaged intensity $`I(\nu )`$ is small, we can write $`E(\nu +\mathrm{\Delta }\nu /2)E^{}(\nu \mathrm{\Delta }\nu /2)=I(\nu )C_E(\nu ,\mathrm{\Delta }\nu )`$, where $`C_E(\nu ,\mathrm{\Delta }\nu )`$ is the normalized field correlation function around the frequency $`\nu `$. In a frequency range in which the change in scattering parameters is small, $`C_E(\nu ,\mathrm{\Delta }\nu )=C_E(\mathrm{\Delta }\nu )`$ and we find

$`E^2(t)`$ $`=`$ $`{\displaystyle 𝑑\nu I(\nu )|F(\nu )|^2𝑑\mathrm{\Delta }\nu \mathrm{exp}[i2\pi \mathrm{\Delta }\nu t]C_E(\mathrm{\Delta }\nu )}`$ (B3)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑tI(t){\displaystyle 𝑑\mathrm{\Delta }\nu \mathrm{exp}[i2\pi \mathrm{\Delta }\nu t]C_E(\mathrm{\Delta }\nu )},`$ (B4)
where the last step follows from Parseval’s theorem. Consequently, the time of flight distribution $`(t)=E^2(t)/𝑑tI(t)`$ and the field correlation function are Fourier transform pairs:

$$(t)=𝑑\mathrm{\Delta }\nu C_E(\mathrm{\Delta }\nu )\mathrm{exp}[i2\pi \mathrm{\Delta }\nu t].$$
(B5)


text:
# Untitled Document

The Uniqueness of the Solution to Inverse Problems of

Interpolation of Positive Operators in Banach Lattices 111The research was supported by the Russian Foundation for Basic Research, grant no. 98-01-00103.

O. E. TIKHONOV and L. V. VESELOVA

Research Institute of Mathematics and Mechanics, Kazan State University, Universitetskaya Str. 17, Kazan, Tatarstan, 420008 Russia

E-mail: Oleg.Tikhonov@ksu.ru

Abstract.  We prove that an interpolation pair of Banach lattices is uniquely determined by the collection of intermediate spaces with the property that these are interpolation spaces for positive operators. A correspondent result for exact interpolation is also presented.

Mathematics Subject Classification (1991): 46B42, 46B70.

Key words: Banach lattice, interpolation pair, positive operator.

In  we have proved that an arbitrary Banach couple (= interpolation pair of Banach spaces) is determined, up to the equivalence of the norms, by the collection of all its interpolation spaces. By making use of this, we have proved in  that a Banach couple is determined by the collection of all its exact interpolation spaces up to the proportionality of the norms. (See also , where inverse interpolation problems were examined in a more general setting.) The goal of the present note is to show that similar assertions hold true when studying interpolation of positive operators in Banach lattices.

We first introduce notation and present some facts concerning interpolation of linear operators in Banach lattices. For basic definitions of the theory we refer the reader to .

For a Banach lattice $`E`$, we denote by $`_E`$ and by $`E^+`$ the norm in $`E`$ and the cone of positive elements, respectively; $`|x|`$ stands for the absolute value of $`xE`$. We denote by $`^E`$ the norm in the dual Banach lattice $`E^{}`$ and use the notation $`T_{EE}`$ to denote the norm of a bounded linear operator $`T`$ in $`E`$. A linear operator $`T`$ in $`E`$ is called positive if $`T`$ maps $`E^+`$ into $`E^+`$. Let $`L^+(E)`$ stand for the cone of all such operators. Recall that a positive operator in a Banach lattice is automatically bounded.

Let $`E`$ and $`F`$ be two Banach lattices with $`EF`$. (Here and subsequently, the notation $`EF`$ for two Banach lattices $`E`$ and $`F`$ means that $`E`$ is an ideal in $`F`$ and there exists a constant $`c>0`$ such that $`x_Ecx_F`$ for all $`xE`$.) If $`\phi F^{}`$ then, clearly, $`\phi |_EE^{}`$. We write $`\phi ^E`$ to denote $`\phi |_E^E`$. Note also that $`|\phi |_E|=|\phi ||_E`$.

Two Banach lattices $`X`$ and $`Y`$ are said to form an interpolation pair $`(X,Y)`$ of Banach lattices if they both are ideals in a certain vector lattice $`U`$ which is a Hausdorff topological vector space such that the embeddings of $`X`$ and of $`Y`$ into $`U`$ are continuous. To each interpolation pair $`(X,Y)`$ of Banach lattices can be canonically associated two Banach lattices, the intersection $`XY`$ and the sum $`X+Y`$. As can be easily seen,

$$(X+Y)^+=\{x+yxX^+,yY^+\}.$$
$`(1)`$
We denote by $`L^+((X,Y))`$ the set of all linear operators $`T`$ in $`X+Y`$ with $`T|_XL^+(X)`$, $`T|_YL^+(Y)`$. One can easily deduce from (1) that $`L^+((X,Y))L^+(X+Y)`$. Also, we denote by $`L_1((X,Y))`$ the set of rank one linear operators $`T`$ in $`X+Y`$ with the property that $`T`$ maps $`X`$ into $`X`$ and $`Y`$ into $`Y`$ continuously.

A Banach lattice $`Z`$ is called intermediate for an interpolation pair $`(X,Y)`$ of Banach lattices if $`XYZX+Y`$. For a constant $`c1`$, we denote by $`\text{P-Int}_c(X,Y)`$ the collection of intermediate lattices $`Z`$ satisfying the following ”interpolation property for positive operators”: $`T(Z)Z`$ and

$$T_{ZZ}c\mathrm{max}\{T_{XX},T_{YY}\}$$
$`(2)`$
for all $`TL^+((X,Y))`$. Set $`\text{P-Int}(X,Y)=\underset{c1}{}\text{P-Int}_c(X,Y)`$. Similarly, we denote by $`\text{R1-Int}_c(X,Y)`$ the collection of intermediate lattices $`Z`$ satisfying the following condition: $`T(Z)Z`$ and (2) holds true for all $`TL_1((X,Y))`$. Set $`\text{R1-Int}(X,Y)=\underset{c1}{}\text{R1-Int}_c(X,Y)`$.

LEMMA 1. Let $`(X,Y)`$ be an interpolation pair of Banach lattices and $`c1`$. Then $`\text{P-Int}_c(X,Y)\text{R1-Int}_c(X,Y)`$.

Proof. First, consider an intermediate Banach lattice $`Z`$ and a continuous linear operator $`T`$ in $`X+Y`$ of rank one. Necessary, $`T`$ is of the form $`T()=\phi ()x`$ with $`\phi (X+Y)^{}`$, $`xX+Y`$. Define $`|T|`$ by $`|T|()=|\phi |()|x|`$. We claim that $`|T|`$ maps $`Z`$ into itself if and only if $`T`$ does, and $`|T|_{ZZ}=T_{ZZ}`$ when $`T(Z)Z`$. Really, the following two cases are possible: $`xZ`$ and $`xZ`$. In the first case we have

$$|T|_{ZZ}=|\phi |^Z|x|_Z=\phi ^Zx_Z=T_{ZZ}.$$

In the second one,

$$T(Z)Z\phi |_Z=0|\phi ||_Z=0|T|(Z)Z$$

and

$$T(Z)ZT=|T|=0.$$

Thus, the claim is proved.

Now, let $`Z\text{P-Int}_c(X,Y)`$. Consider an operator $`TL_1((X,Y))`$. $`T`$ is clearly continuous in $`X+Y`$ and we can define $`|T|`$ as above. Then $`|T|L^+((X,Y))`$. Since $`Z\text{P-Int}_c(X,Y)`$, we have $`|T||_ZL^+(Z)`$. The latter implies that $`T`$ maps $`Z`$ into itself. Moreover,

$$\begin{array}{cc}\hfill T_{ZZ}& =|T|_{ZZ}\hfill \\ & c\mathrm{max}\{|T|_{XX},|T|_{YY}\}\hfill \\ & =c\mathrm{max}\{T_{XX},T_{YY}\}.\hfill \end{array}$$

Thus, $`Z\text{R1-Int}_c(X,Y)`$.

For two Banach lattices $`E`$ and $`F`$, let the notation $`EF`$ mean that $`E`$ and $`F`$ coincide as vector lattices and their norms are equivalent, the notation $`EF`$ mean that the norms are proportional.

For two interpolation pairs of Banach lattices $`(X,Y)`$ and $`(V,W)`$, we write $`(X,Y)(V,W)`$ if $`XV`$ and $`YW`$ or $`XW`$ and $`YV`$. We write $`(X,Y)(V,W)`$ if $`XV`$ and $`YW`$ or $`XW`$ and $`YV`$.

THEOREM 1. Let $`(X,Y)`$ and $`(V,W)`$ be two interpolation pairs of Banach lattices with

$$\text{P-Int}(X,Y)=\text{P-Int}(V,W).$$

Then $`(X,Y)(V,W)`$.

Proof. Since $`X,Y\text{P-Int}(X,Y)`$ and $`V,W\text{P-Int}(V,W)`$, we obtain $`X,Y\text{R1-Int}(V,W)`$ and $`V,W\text{R1-Int}(X,Y)`$ by making use of Lemma 1. So it remains to prove

$$X,Y\text{R1-Int}(V,W)\text{ and }V,W\text{R1-Int}(X,Y)(X,Y)(V,W).$$

By analyzing proofs in  one can verify that this implications was actually proved there. Note also that from \[4, Proposition 1.4 and Theorem 2.8\] the latter implication follows immediately.

THEOREM 2. Let $`(X,Y)`$ and $`(V,W)`$ be two interpolation pairs of Banach lattices with

$$\text{P-Int}_1(X,Y)=\text{P-Int}_1(V,W).$$

Then $`(X,Y)(V,W)`$.

Proof. We can proceed similarly to the proof of Theorem 1, but here we need to refer to  and to \[4, Theorem 3.7\].

References

1. J. Linderstrauss and L. Tzafriri, Classical Banach Spaces II. Function Spaces, Springer-Verlag, Berlin-Heidelberg-New York, 1979.

2. O. E. Tikhonov and L. V. Veselova, A Banach couple is determined by the collection of its interpolation spaces, Proc. Amer. Math. Soc. 126 (1998), 1049–1054.

3. O. E. Tikhonov and L. V. Veselova, The uniqueness of the solution to the inverse problem of exact interpolation, Israel Math. Conf. Proc. 13 (1999), 209–215.

4. L. V. Veselova and O. E. Tikhonov, The uniqueness of the solution to inverse interpolation problems, Research Institute of Mathematics and Mechanics, Preprint no. 95-2, Kazan Mathematics Foundation, Kazan, 1995 (Russian); English transl., XXX e-Print Archive, math.FA/9902108.


text:
# Reconstruction of Gray-scale Images

## 1 Introduction

We consider here the problem of image reconstruction. Suppose a multicolor or gray-scale picture is subjected to noise and an observer has access only to this corrupted version. How can she estimate the original picture?

The analysis of this kind of problem has attracted a lot of interest and many approaches have been considered (\[G\]).

One of the methods proposed is the so called MAP estimator. In this method one assumes that the original image is a random realization of a Markov random field that has been corrupted by some site independent noise. One assumes that the distribution of the field (a priori distribution) is known as well as the distribution of the noise, that is, the conditional distribution of the observed image given the original one. The MAP estimator is the image that has the largest probability of have produced the observed one. This is the mode of the posterior distribution.

In multicolor images the literature proposes to use the so called Potts model as a priori distribution. Roughly speaking, this model is a measure on the set of images that gives more weight to images that have neighboring pixels of the same color. This choice has a number of advantages and disadvantages. An important disadvantage is the fact that the algorithms used to compute the MAP estimator operate in exponential time in the number of pixels (\[GJ\]). Our main point in this paper is to propose an alternative a priori measure with the property that the computation time of the exact MAP estimator is polynomial in the number o pixels. Using this approach we have produced an algorithm and a program that reconstruct dirty images in polynomial time in the number of pixels.

In the remaining of this section we explain these ideas in some detail. In the next section we present some experimental results obtained from the implementation of our method to some images. The next section is more technical in nature. It explains why our MAP estimator can be computed efficiently (polynomial time) and provides a proof for the theoretical result presented below. The question of how to obtain estimators for some important parameters associated to the observed image is discussed in the appendix. We close this paper with some final remarks.

In order to motivate the discussion we have to introduce some notation. Assume the image is a point of $`\{1,\mathrm{},c\}^\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is a finite subset of $`𝒵^2`$ (e.g. a square) with $`N`$ sites. The image associates to each site or pixel in $`\mathrm{\Lambda }`$ one of $`c`$ possible colors. Typically values for $`c`$ are $`c=2^8`$ or $`c=2^{24}`$. Assume $`1K=\mathrm{log}_2c𝒵`$ and let $`\underset{¯}{X}=(X_i)_{i=1,\mathrm{},N}`$ for $`X_i=(X_i^1,X_i^2,\mathrm{},X_i^K)\{0,1\}^K,i\mathrm{\Lambda }`$, $`K1`$, represent the true unknown (random) image. Denote the observed image by $`\underset{¯}{Y}=(Y_i)_{i=1,\mathrm{},N}`$. We denote the space of pictures $`\{\{0,1\}^K\}^\mathrm{\Lambda }`$ by $`\mathrm{\Sigma }_K`$. In this notation we have $`c=2^K`$ possible colors in each pixel. Each $`X_i`$ may correspond to a single binary number which gives the intensity of black in pixel $`i`$ as $`_{k=1}^KX_i^k2^k`$ (in the case of gray-scale picture) or it may correspond to three binary numbers, each giving the intensity of one of the three basic colors. In this second case, for instance, we could have $`K=24`$ and $`_{k=1+c8}^{(1+c)8}X_i^k2^k`$, for $`c=0,1`$ and $`2`$ giving the intensity of, respectively, Red, Green and Blue at site $`i`$. To simplify the discussion we assume, without loss, the first case since our approach in the second one is to reconstruct each color separately in order to reconstruct the whole picture.

In the Bayesian setting we assume (a) that the original picture $`\underset{¯}{X}`$ is random with a known distribution which is called the a priori measure and (b) that we know how to model the noise.

To start we consider the noise. We fix an $`ϵ>0`$. Conditional on $`\underset{¯}{X}`$, at each bit $`k`$ of each pixel $`i`$, independently, the observed value in bite $`k`$ of pixel $`i`$, called $`Y_i^k`$, is equal to the true value $`X_i^k`$ with probability $`1ϵ`$ or, with probability $`ϵ`$, it corresponds to the switched value:

$$P_i(Y_i^k|X_i^k)=\{\begin{array}{cc}1ϵ& \text{if}X_i^k=Y_i^k\hfill \\ ϵ& \text{otherwise.}\hfill \end{array}$$
(1)
Hence

$$P(\underset{¯}{Y}|\underset{¯}{X})=\underset{i\mathrm{\Lambda }}{}\underset{k\{1,\mathrm{},K\}}{}P_i(Y_i^k|X_i^k)$$
(2)
$$=\frac{1}{Z_h}\mathrm{exp}\left\{h\underset{i\mathrm{\Lambda }}{}\underset{k\{1,\mathrm{},K\}}{}\text{1}_{\{X_i^k=Y_i^k\}}(\underset{¯}{Y})\right\}$$
(3)
where

$$h=\mathrm{log}\frac{(1ϵ)}{ϵ}.$$
(4)
and $`Z_h`$ is the normalization constant.

Note that this hypothesis is not the same as the one in \[FFG\] where each pixel was either observed correctly with probability $`1ϵ`$ or chosen uniformly among the other $`c1`$ colors.

Now we define the MAP (maximum a posteriori) estimator.

Denote by $`\mu `$ the a priori measure, that is, the distribution of $`\underset{¯}{X}`$. Then the MAP estimator after observing the image $`\underset{¯}{Y}`$, $`\underset{¯}{\overset{^}{X}}=\underset{¯}{\overset{^}{X}}(\underset{¯}{Y})\mathrm{\Sigma }_K`$ is any image which maximizes the posterior distribution $`P(\underset{¯}{X}|\underset{¯}{Y})P(\underset{¯}{Y}|\underset{¯}{X})\mu (\underset{¯}{X})`$, that is

$$P(\underset{¯}{\overset{^}{X}}|\underset{¯}{Y})=\underset{\underset{¯}{X}\mathrm{\Sigma }_K}{\mathrm{max}}P(\underset{¯}{X}|\underset{¯}{Y})$$
(5)
Now we discuss the a priori measure. We want to consider here the situation on which very little is known about the original picture before the observation besides the information that it is some kind of real photo, perhaps taken by a satellite, as opposed to being a photo of something like a geometrical drawing, a cubist oil painting or a cell of a cartoon picture. We suppose that the only prior information available is that the original picture is locally smooth. By this we mean that the measure should be such that neighbor pixels on the picture are more likely to have colors which are near in some sense. More precisely we assume that there exists a real valued function $`H`$ which indicates how smooth  is an image and take the a priori measure to be:

$$\mu (\underset{¯}{X})=\frac{e^{\beta H(\underset{¯}{X})}}{Z}$$
(6)
where $`Z=_{\underset{¯}{X}}e^{\beta H(\underset{¯}{X})}`$, with the sum taken over all possible images, is a normalization constant and $`\beta 0`$ is a real parameter. This parameter measures the tendency to be smooth since if $`\beta `$ is large $`\mu `$ is concentrated on images with small values of $`H`$ while if $`\beta `$ is small $`\mu `$ is close to the uniform distribution on $`\mathrm{\Sigma }_K`$, the set of all images. The motivation for this formula, in particular the minus sign in front of $`\beta `$, comes from statistical mechanics where it would be called Gibbs measure, $`\beta `$ would be the inverse temperature and $`H`$ would be the Hamiltonian (or Energy) function.

Plugging (3) into (5), and taking logarithms, we get that the images which maximize the posterior distribution (5) are those which maximize

$$\beta H(\underset{¯}{X})+h\underset{i\mathrm{\Lambda }}{}\underset{k\{1,\mathrm{},K\}}{}\text{1}_{\{X_i^k=Y_i^k\}}(\underset{¯}{Y})$$
(7)
The image that maximizes this expression makes the best compromise between being globally smooth, regulated by $`\beta H`$, and agreeing as much as possible with the observed image which is regulated by $`h_{i\mathrm{\Lambda }}_{k\{1,\mathrm{},K\}}\text{1}_{\{X_i^k=Y_i^k\}}(\underset{¯}{Y})`$. Of course there exists only one relevant parameter in the maximization problem, say $`h/\beta `$.

Although any strictly positive measure $`\mu `$ on $`\mathrm{\Sigma }_K`$ can be written as above for some  $`H`$, and thus (6) can be safely assumed for a generic a priori  measure, we are interested only in the case on which $`H`$ is both simple and reasonable  as a measurement of smoothness. This simplicity requirement will basically mean the assumption that $`H`$ has only local dependence as follows

$$H(\underset{¯}{X})=\underset{<i,j>}{}d(X_i,X_j)$$
(8)
where $`d(X_i,X_j)`$ is a measure for the distance between the colors at pixel $`i`$ and $`j`$ and the sum is taken over all pairs of neighbor pixels in the picture. That is, $`<i,j>=\{(i,j):|ij|=1\}`$.

The choice of $`d(X_i,X_j)`$ depends on what kind of local properties one expects in the original picture.

One choice which is common in the literature (see, for instance \[FFG\] and its quotations) is

$$H_P(\underset{¯}{X})=\underset{<i,j>}{}1\mathrm{I}_{\{X_iX_j\}}$$
(9)
where the sum is taken over all pairs of neighbor sites in the lattice.

This choice is related to the Potts model in statistical mechanics (see for instance \[M\]) and explains the subscript. It is a model with very interesting properties that corresponds, in the two-color case, to the Ising model \[MW\]. To find a solution to (5) we have to find what is called in the physics literature a ground state for the Ising model with random magnetic field. This field is induced by the observed image.

Most of the interesting properties of those statistical mechanics models, like phase transition, appear in the so called thermodynamic limit, the limit on which the size of the system grows to the whole lattice (\[R\], \[MW\]). Even though in our discussion the lattice size is kept fixed we will need in section 4 some exact results on the thermodynamic limit in the two-color case (Ising model). Also some finite size considerations like the effect of boundary conditions (in our case we chose free boundary conditions) may be important.

The measure defined by (6), with (9) plugged in, would describe images on which neighbor pixels tend to be equal but if two pixels have different colors the cost of this interface does not depend on how different  they are. If two pixels have different colors the most likely is that each one belongs to a single-color region separated by a sharp interface.

We could also try to represent the situation on which this is not necessarily the case choosing a finer notion of distance between $`X_i`$ and $`X_j`$. Two somewhat natural choices would be

$$H_1(\underset{¯}{X})=\underset{<i,j>}{}\left|\underset{k=1}{\overset{K}{}}(X_i^kX_j^k)\frac{1}{2^k}\right|$$
(10)
or

$$H_2(\underset{¯}{X})=\underset{<i,j>}{}\left(\underset{k=1}{\overset{K}{}}(X_i^kX_j^k)\frac{1}{2^k}\right)^2.$$
(11)
Note that, for any choice of $`H`$, the MAP problem is well posed since we could, at least in principle, check the finitely many values in the right hand side of (7) and choose a image that maximizes it. But since $`|\mathrm{\Sigma }_K|=2_{}^{K}{}_{}{}^{N}`$ with something like $`K=24`$ and $`N=400\times 600`$ this approach is not computationally feasible. The problem with all the above choices is that it is not known how to do better than this time-consuming maximization by inspection in the multicolor case ($`K>1`$) and thus the problem is, in practice, not solvable.

We propose another choice for $`H`$, intermediate between the Potts and the other two mentioned above, which respects nicely the notion of smoothness which led to the choice (10) but nevertheless induces to a polynomial time maximization problem. This choice is:

$$H(\underset{¯}{X})=\underset{k=1}{\overset{K}{}}\underset{<i,j>}{}|X_i^kX_j^k|\frac{1}{2^k}=\underset{k=1}{\overset{K}{}}\underset{<i,j>}{}1\mathrm{I}_{\{X_i^kX_j^k\}}\frac{1}{2^k}$$
(12)
With this, the problem of finding the $`2^K`$-color image which maximizes the posterior distribution $`P(\underset{¯}{X}|\underset{¯}{Y})P(\underset{¯}{Y}|\underset{¯}{X})\mu (\underset{¯}{X})`$ is decomposed into $`K`$ binary color maximization problems. Namely one has to solve

$$P(\underset{¯}{\overset{^}{X}}_k|\underset{¯}{Y}_k)=\underset{\underset{¯}{X}\mathrm{\Sigma }_2}{\mathrm{max}}P(\underset{¯}{Y_k}|\underset{¯}{X})\mu _k(\underset{¯}{X})$$
(13)
for each $`k`$, $`1kK`$ where $`\underset{¯}{Y}_k=\{Y_i^k\}_{\{i\mathrm{\Lambda }\}}`$ is the $`k`$-th component of the observed image and $`\mu _k`$ is the Gibbs measure defined on the space of binary color images, $`\mathrm{\Sigma }_2`$, by (6) with

$$H^k(\underset{¯}{X}_k)=\underset{<i,j>}{}1\mathrm{I}_{\{X_i^kX_j^k\}}\frac{1}{2^k}.$$
(14)
The $`K`$ color image given by $`\underset{¯}{Y}`$ is a solution of (5). Each one of these $`K`$ binary images is called a layer. The weighted sum of the MAP for each layer gives our (gray-scale) image estimator:

$$\widehat{X}_i=\underset{k=1}{\overset{K}{}}2^k\widehat{X}_i^k.$$

On the multicolor case one solves a gray-scale problem for each basic color.

Each binary problem can be solved in polynomial time using the results by Greig, Porteous and Seheult who reformulated it as one involving finding a minimum cut on a capacitated network \[FF\] for which there exist fast algorithms. These ideas are presented briefly in the next section. From the statistical point of view the approach is different whether the parameters $`\beta `$ and $`ϵ`$ are known or not. A truly Bayesian approach would associate an a priori measure to each one of the parameters $`\beta `$ and $`ϵ`$. We leave this alternative to future work. Another possibility is to estimate these parameters from the observed image using classical frequentist analysis. This is possible under the hypothesis on the noise and on the original image being a sample of a Gibbs measure for the Potts Hamiltonian (9) (which is always the case after the decomposition in binary colors) using exact results for the two-dimensional Ising model obtained by Frigessi and Piccioni \[FP\]. We have produced an algorithm based on \[FP\] for estimating the parameters $`ϵ`$ and $`\beta `$. This is explained in section 4 below.

Once one has a fast algorithm to reconstruct images it is natural to ask what is the effect of iterating the whole process. More precisely, what happens if one takes the reconstructed image and applies the method again, perhaps updating the value of $`h/\beta `$?

Note that to apply the method again for the (already) reconstructed image is equivalent to assume that the reconstructed image could be thought as obtained from some original picture which was chosen with respect to $`\mu `$ and then subjected to noise. Even though it is not difficult to verify that this assumption is false the question is interesting both from the mathematical and from the applied point of view. On the mathematical side one has a mapping, $`\underset{¯}{Y}\underset{¯}{\overset{^}{X}}(\underset{¯}{Y})`$, from $`\mathrm{\Sigma }_K`$ into itself, a discrete time dynamical system, and it is natural to ask about iteration properties. On the practical minded side one can ask about the effect of iteration on the quality of reconstruction even if this can only be judged subjectively.

If we decide to update $`h/\beta `$ before each iteration the natural thing to do would be to chose a larger value at each step since we assume the procedure did a good job and removed some of the noise, therefore $`\underset{¯}{\overset{^}{X}}(\underset{¯}{Y})`$ would be more reliable than $`\underset{¯}{Y}`$. It would then be natural to increase $`h`$, without changing $`\beta `$. A natural procedure to find each updated value is to use again the estimators given by Frigessi and Piccioni (\[FP\]). In doing this we find, in fact, that this parameter increases, after one iteration, as expected.

What is rather surprising is that once one tries this iteration procedure, either keeping $`h/\beta `$ fixed or increasing it, at each step, one finds that it has no effect et all. The twice reconstructed image is exactly equal to the once reconstructed one. In other words, $`\underset{¯}{\overset{^}{X}}`$ as a function from $`\mathrm{\Sigma }_K`$ into itself has a fixed point in each once-reconstructed image.

More precisely write $`F_{(h/\beta )}(\underset{¯}{X})=\underset{¯}{\overset{^}{X}}`$ for the function from $`\mathrm{\Sigma }_K`$ into itself defined in (5). Then we have the following

Proposition: $`F_s(\underset{¯}{X})=F_t(F_s(\underset{¯}{X}))`$ for all $`ts`$.

We prove this proposition in the next section after reviewing some results on networks and the connection with the maximization problem considered in this note.

Given this result a natural question arises. Is this property not true in general? Suppose one has a random function assuming values in some space $`S`$ with some unknown parameter $`s`$ which itself belongs to $`S`$ chosen according to some a priori measure (in our case, $`s`$ corresponds to the original image). Define the MAP estimator as usual to obtain a function from $`S`$ into itself. Does it always satisfy the this fixed point property? One could argue that since all the information about the original image which is contained in the observed one should be still present in the reconstructed image then this iteration should give no further results and the twice reconstructed image should always be equal to the once-reconstructed one. As it turns out this is not the case and it is not difficult to find counterexamples.

## 2 Experimental Results

We developed a program (map) in *C* to reconstruct images according to the method proposed. This program reads an image in the *portable bit map* format for true-color (*.ppm*) or gray-scale (*.pgm*) pictures. The intensity of each color in each pixel (or the single gray scale intensity there) is written as a (eight bits) binary number. The program then constructs a graph for each one of those layers (as defined after equation (14)), construct the corresponding graphs and finds the minimum cut for each one of them using a standard implementation of the Ford Fulkerson algorithm. The program uses integer arithmetic scaled by a factor of 10000. The values for $`\beta `$ and $`H`$ are supplied by command line switches. It is also possible to supply a bit mask to select which planes to reconstruct. The program is portable and it has been tested on *SunOS*, *Solaris* and *Linux*.

The main problem in developing the program was the large amount of memory to hold the image, the graph and the auxiliary data structures. The solution was to use a compact representation of the graph, namely a matrix where each cell has a pixel and the values for the flow in each of 5 directions (4 neighbors plus $`s`$ or $`t`$). This avoided the explicit representation of the edges and the extra space for the image itself.

A series of tests were made on a 700Mhz *Pentium III* machine, with 128Mb of main memory, running Linux Debian potato 2.2, kernel version 2.4. The running times for restoration of six sample images are presented in table 1. These images are shown in tables 2, 3 and 4.

The images were obtained from a picture shot with a digital camera and converted to 8 bit gray-scale. The final result was saved in the file “fish.pgm”. Two aditional files were obtained by clipping the image to a rectangle with half of the area (“fish2.pgm”) and to another one with a quarter of the total area (“fish3.pgm”).

The noise was added “artificially” by a program which scan every bit in the image and inverts it with a given probability. The probabilies used were 10% and 15%. The modified files were named fish.10.pgm, fish2.10.pgm, fish3.10.pgm and fish.15.pgm, fish2.15.pgm, fish3.15.pgm, respectively.

Each image was restored in all bit planes using 3 values of $`\beta `$, namely 0.1, 0.3, 0.5. An aditional restoration, just on the most significant bit plane were made using the estimated $`\widehat{\beta }`$ as indicated in the appendix. The values are presented in the tables 2, 3 and 4.

We did not try to define a metric to measure the quality of the reconstruction, relying on a subjective analysis. Besides the removal of noise, it was observed an improvement in the shades and (consequently) in the third dimension perception, even if some blurring is introduced.

These and other examples are available directly from the authors or at `http://www.ime.usp.br/~gubi/MAP/index.html`

## 3 MAP estimator and networks

In this section we present some capacitated network ideas which are used to solve the two-color maximization problem described before as discovered by Greg, Porteous and Seheult \[GPS\] and use them to prove Proposition.

As mentioned before with the choice of $`H`$ given by (12) the maximization in (5) is decomposed into $`K`$ two-color problems. Therefore we assume in this section that we are in the binary color case $`K=1`$.

The MAP estimator given $`\underset{¯}{Y}`$ is the image $`\underset{¯}{\overset{^}{X}}\{0,1\}^\mathrm{\Lambda }`$ which gives the maximum of

$$L_\alpha (\underset{¯}{X}|\underset{¯}{Y})=\alpha \underset{i\mathrm{\Lambda }}{}\text{1}_{\{X_i=Y_i\}}(\underset{¯}{X})+\underset{<i,j>}{}\text{1}_{\{X_i=X_j\}}(\underset{¯}{X})$$
(15)
where $`\alpha =h/\beta `$.

A network $`N`$ is a graph $`G=(V,E)`$, where $`V`$ is a finite set of vertices and $`E`$ is a set of couples of vertices, with a capacity $`c(e)0`$ associated with each edge $`eE`$.

We define networks on the set of vertices $`V`$, given by the sites in $`\mathrm{\Lambda }`$ plus two extra ones denoted by $`s`$ (source) and $`t`$ (sink), that is

$$V=\mathrm{\Lambda }\{s\}\{t\}.$$

For the set of arcs we take

$$E=\left(_i\{(s,i)\}\right)\left(_i\{(i,t)\}\right)\left(_{<i,j>}\{(i,j)\}\right)$$

where the first two unions are taken on $`i\mathrm{\Lambda }`$ and the last one is taken over pairs of nearest neighbor sites in $`\mathrm{\Lambda }`$.

Given the observed image $`\underset{¯}{Y}`$ and a real number $`\alpha `$ we define the capacities of the network $`N_\alpha (\underset{¯}{Y})`$ as follows.

If $`Y(i)=1`$ we set $`c_{\underset{¯}{Y},\alpha }(s,i)=\alpha `$ as its capacity, otherwise set $`c_{\underset{¯}{Y},\alpha }(i,t)=\alpha `$; to each arc $`e=(i,j)`$ of neighbor sites in $`\mathrm{\Lambda }`$ we associate $`c_{\underset{¯}{Y},\alpha }(e)=1`$. All other arcs have capacity zero.

For each image $`\underset{¯}{X}`$ let

$$A(\underset{¯}{X})=\{s\}\{i\mathrm{\Lambda }:X(i)=1\}$$

$$B(\underset{¯}{X})=\{t\}\{i\mathrm{\Lambda }:X(i)=0\}.$$

These two sets define a cut of the network

$$𝒞(\underset{¯}{X})=\{(i,j)E:iA(\underset{¯}{X}),jB(\underset{¯}{X})\}.$$

Notice that the cut $`𝒞(\underset{¯}{X})`$ consists of a set of arcs whose removal (cut) makes it impossible to find a path going from $`s`$ to $`t`$ through arcs with non-zero capacity. If we now define the capacity of the cut $`𝒞(\underset{¯}{X})`$ by the quantity

$$C_{\underset{¯}{Y},\alpha }(\underset{¯}{X})=\underset{(i,j)𝒞(\underset{¯}{X})}{}c_{\underset{¯}{Y},\alpha }(i,j)$$
(16)
it is very simple to check that

$$L_\alpha (\underset{¯}{X}|\underset{¯}{Y})=aC_{\underset{¯}{Y},\alpha }(\underset{¯}{X}),$$

where $`a`$ is a constant which does not depend on $`\underset{¯}{X}`$. Therefore to find the MAP estimator one has to find the cut which minimizes (16): the so called minimum cut \[FF\].

Ford and Fulkerson showed that the value of the capacity of the minimum cut is equal to the maximum flow  through the network from source to sink. Recall that a flow $`𝐟`$ in a network $`N`$ (on $`G(V,E)`$ with capacity $`\{c(e)\}_{eE}`$) from $`s`$ to $`t`$ is a collection of real numbers $`\{𝐟(e)\}_{eE}`$, where $`𝐟(e)`$ can be thought as the amount of fluid per unit time going through the pipeline  $`e`$, such that the flow on each arch does not exceed its capacity, $`0𝐟(e)c(e)`$, for all $`eE`$, and such that the flow is conserved

$$\underset{jV}{}𝐟(i,j)\underset{jV}{}𝐟(j,i)=0$$
(17)
for all $`is`$.

Well known fast algorithms to find this maximum flow exist and therefore the problem is solved in practice. The algorithm presented in this paper uses this method to reconstruct each one of the (8 or 24) layers.

Proof of the Proposition. To simplify the notation write $`\underset{¯}{\overset{^}{X}}=F_\alpha (\underset{¯}{Y})`$. Let $`\alpha ^{}\alpha `$. The networks defined by the pair $`(\underset{¯}{Y},\alpha )`$, denoted by $`N_\alpha (\underset{¯}{Y})`$, and by $`(\underset{¯}{\overset{^}{X}},\alpha ^{})`$, $`N_\alpha ^{}(\underset{¯}{\overset{^}{X}})`$, can differ only in the capacities assigned to each arc.

The proposition asserts that the cut defined by $`\underset{¯}{\overset{^}{X}}`$ also has minimum capacity in $`N_\alpha ^{}(\underset{¯}{\overset{^}{X}})`$, that is we want to prove that

$$C_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(\underset{¯}{\overset{^}{X}})C_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(\underset{¯}{Z}),$$
(18)
for any image $`\underset{¯}{Z}\mathrm{\Sigma }_1`$.

Fix $`\underset{¯}{Z}`$ and let $`E_I=𝒞(\underset{¯}{\overset{^}{X}})𝒞(\underset{¯}{Z})`$, $`E_{II}=𝒞(\underset{¯}{\overset{^}{X}})𝒞(\underset{¯}{Z})`$ and $`E_{III}=𝒞(\underset{¯}{Z})𝒞(\underset{¯}{\overset{^}{X}})`$. By definition of MAP estimator we have

$$C_{\underset{¯}{Y},\alpha }(\underset{¯}{\overset{^}{X}})C_{\underset{¯}{Y},\alpha }(\underset{¯}{Z}),$$
(19)
which implies

$$D_{\underset{¯}{Y},\alpha }(E_I)D_{\underset{¯}{Y},\alpha }(E_{III}).$$
(20)
where, for a set of edges $`E`$,

$$D_{\underset{¯}{Y},\alpha }(E)=\underset{eE}{}c_{\underset{¯}{Y},\alpha }(e)$$

Therefore to check (18) it is enough to verify

$$D_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(E_I)D_{\underset{¯}{Y},\alpha }(E_I)$$
(21)
$$D_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(E_{III}))D_{\underset{¯}{Y},\alpha }(E_{III}).$$
(22)
We start with inequality (21). Assume $`eE_I`$. If $`e`$ is an external edge (connecting some site $`i`$ with either the source $`s`$ or the sink $`t`$), then, since $`e𝒞(\underset{¯}{\overset{^}{X}})`$, it must be $`c_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(e)=0`$ which is less or equal than $`c_{\underset{¯}{Y},\alpha }(e)`$. On the other hand, if $`e`$ is an internal edge (i.e. $`e=(i,j)`$ for $`i`$ and $`j`$ nearest neighbors in $`\mathrm{\Lambda }`$) then its capacity equals $`1`$ in both cases. This concludes the proof of (21).

Suppose now $`eE_{III}`$. Inequality (22) follows from the observation that $`c_{\underset{¯}{\overset{^}{X}},\alpha ^{}}(e)`$ can not be zero.

## 4 Final Remarks

The different approaches used in image reconstruction are based in quite different set of theoretical ideas and it is not clear how to compare their results. One possible measure for the quality  of the reconstruction, used in \[FP\] to compare 9 algorithms, is to evaluate the proportion of pixels classified correctly. Since our main goal here was to present a method for the reconstruction of multicolor images we leave the comparison with other methods for future work. In any case, since we are working with exact MAP’s for the chosen Hamiltonian, our method will be as good (and as bad) as the usual two colors MAP estimators, regarding the proportion of bits reconstructed.

For the usual MAP reconstruction problem in the multicolor case no fast algorithm is known (\[FFG\]). For a probabilistic approach via Simulated Annealing in order to get the exact estimator one needs to decrease very slowly some parameter while the computation goes on and thus needs a prohibitively large amount of time \[Gi\].

One possibility is to accept approximated MAP estimators which can be obtained fast enough. One can do this with simulated annealing by updating the parameter fast enough but then we lose control on how close the approximation is to the exact one. An approximate of the MAP estimator with a probabilistic analysis of the error in the three color case was developed in \[FFG\]. Another approach can be found in \[J\].

Our approach is not completely Bayesian as we also consider a situation on which the a priori measure has some unknown parameters but do not assume any prior knowledge about them. In the Appendix we describe a method to estimate these parameters from the picture itself using classical statistical methods, as proposed by Frigessi and Piccioni \[FP\]. In those cases we verify that plugging those estimated values into the formulas used to get the MAP estimator provides a reconstruction which appears to be the best. As mentioned before we do not try to quantify this.

Ricardo Maronna proposed that instead of estimating $`\beta `$ and $`h`$ one could look for the $`\alpha =\beta /h`$ that maximizes (7). The estimator for the true image will then be the $`\underset{¯}{X}`$ which realizes this maximum with the best $`\alpha `$. It is not clear yet how to justify this theoretically. Another alternative would be to choose the uniform distribution on an appropriate range as a priori distribution for $`ϵ`$ and $`\beta `$. However the computation of the MAP estimator in this case seems prohibitive.

As an experimental observation we remark that reconstructing only the first layer of a dirty image (and leaving the others as they are) gives a quite good visual result. A possible explanation of this fact is that, as a consequence of the binary decomposition, each layer is “half” as important as the previous one.

Looking for algorithms which give smooth solutions we propose the following hierarchical procedure. First consider layer $`1`$ and find $`\widehat{Z}^1`$, as before for $`\widehat{X}^1`$. Then given layer $`\widehat{Z}^{\mathrm{}}`$, for $`\mathrm{}=1,\mathrm{},k1`$, define $`\widehat{Z}^k`$ as the (binary) image that maximizes

$$L_\alpha (\underset{¯}{X}_k|\underset{¯}{Y})=\alpha \underset{i\mathrm{\Lambda }}{}\text{1}_{\{X_i^k=Y_i^k\}}(\underset{¯}{X})+\underset{<i,j>}{}\underset{\mathrm{}=1}{\overset{k}{}}\text{1}_{\{X_i^{\mathrm{}}=X_i^{\mathrm{}}\}}(\underset{¯}{X}).$$
(23)
In words, this algorithms will try to get the same value for two neighbors in layer $`k`$ if these neighbors have the same value for all previous layers. If not they are not coupled. Visual realizations of this algorithm give also good results. Although there is no Hamiltonian for this model, the algorithm is well defined. Moreover each layer corresponds to the so called diluted Ising model.

Assume that we get $`N`$ independent samples of the same image $`\underset{¯}{Y}_1,\mathrm{},\underset{¯}{Y}_N`$. This means that the original realization $`\underset{¯}{X}`$ of the image is the same but the noises are independent. A generalization of our approach deals with this problem in the following way. We modify the capacities associated to the graph $`G`$. To each arc $`(s,i)`$ connected to the source we associate the capacity

$$\alpha \left(\underset{n=1}{\overset{N}{}}(2Y_{n,i}^k1)\right)^+.$$

and then to each arc $`(i,t)`$ connected to the sink we associate the capacity

$$\alpha \left(\underset{n=1}{\overset{N}{}}(2Y_{n,i}^k1)\right)^{}.$$

In the same vein one could use the remaining two colors to get information about the color being reconstructed.

## 5 Appendix:  Estimation of $`\beta `$ and $`ϵ`$

In this section we apply ideas from Frigessi and Piccioni (\[FP\]) which exploit well known (but highly nontrivial) results from statistical mechanics (\[R\], \[MW\]) to obtain estimators for the parameters $`\beta `$ and $`ϵ=(1+e^h)^1`$.

The result is a program called estima that has as input a multicolor image and returns one estimated $`\beta `$ and one estimated $`ϵ`$ for each layer.

The measure defined in (6) with $`H`$ given by (9) in the two-color case (Ising model) in the thermodynamic limit gives rise to a translation invariant measure which may be ergodic or not according to the value of $`\beta `$.

This choice of $`H`$ does not favor zeros or ones. This symmetry is perhaps easier to see if we represent a configuration as an element $`\underset{¯}{S}`$ in $`\{1,+1\}^{𝒵^2}`$ with $`1`$ replacing zeros. This is the usual Ising notation while the one using zeros and ones is known in statistical mechanics as lattice gas representation. In the situation considered here they are equivalent.

More precisely in if $`\underset{¯}{Y}\{0,1\}^\mathrm{\Lambda }`$ let $`\underset{¯}{S}`$ be the corresponding Ising configuration given by

$$\underset{¯}{S}=\{S_i,i\mathrm{\Lambda }:S_i=2Y_i1\}$$
(24)
And instead of Hamiltonian (9) we use

$$H_P(\underset{¯}{S})=\underset{<i,j>}{}S_iS_j$$
(25)
The corresponding infinite volume limit measure is not ergodic if $`\beta `$ is smaller than some (known) critical value $`\beta _c`$. In this case the limit measure is a mixture with equal weights of two measures $`\mu _\beta ^+`$ and $`\mu _\beta ^{}`$, the first favoring configurations with more $`+1`$’s than $`1`$’s and the other favoring configurations with more $`1`$’s than $`+1`$’s. If $`\beta \beta _c`$ then $`\mu _\beta ^+=\mu _\beta ^{}`$. If we denote by $`E_\beta ^+(f)`$ ( $`E_\beta ^{}(f)`$ ) the expected value of a function $`f`$ defined on $`\{1,+1\}^{𝒵^2}`$ the symmetry between $`\mu _\beta ^+`$ and $`\mu _\beta ^{}`$ imply

$$E_\beta ^+(S_iS_j)=E_\beta ^{}(S_iS_j)=r_\beta (|ij|)$$
(26)
$`r_\beta (|ij|)`$ is called the two point correlation function for the infinite volume system at inverse temperature $`\beta `$. We will need two of those correlation functions in what follows: $`r_\beta (1)`$, for nearest neighbors, and $`r_\beta (\sqrt{2})`$ for neighbors along the diagonal on the lattice.

Suppose $`\underset{¯}{S}`$ corresponds (in the Ising notation) to the original image and $`\underset{¯}{R}`$ corresponds to the observed image after noise. Under the hypothesis on the noise

$$E_\beta ^+(R_iR_j)=E_\beta ^+(S_iS_j)(12ϵ)^2$$
(27)
Therefore the ratio between $`r_\beta (1)`$ and $`r_\beta (\sqrt{2})`$ computed for the observed image depends only on $`\beta `$. Call this ratio $`\varphi (\beta )`$.

If $`\underset{¯}{R}\{1,1\}^\mathrm{\Lambda }`$ corresponds to an observed two-color picture \[FP\] found a sequence $`(\widehat{\beta },\widehat{ϵ})`$ of consistent estimators for $`(\beta ,ϵ)`$ given by

$$\widehat{\beta }=\varphi ^1\left(\frac{G^1(\underset{¯}{R})}{G^2(\underset{¯}{R})}\right)$$
(28)
and

$$\widehat{ϵ}=\frac{1}{2}\left\{1\left(\frac{G^1(\underset{¯}{R})}{r_{\widehat{\beta }}(1))}\right)^{\frac{1}{2}}\right\}$$
(29)
with

$$G^a(\underset{¯}{R})=\frac{_{(i,j)\mathrm{\Lambda }^0}^aR_iR_j}{4|\mathrm{\Lambda }^0|}$$
(30)
where for a =1 the sum is taken over pairs of nearest neighbor sites along the lattice directions in $`\mathrm{\Lambda }^0`$, the interior of $`\mathrm{\Lambda }`$, and for a = 2 the sum is over neighbor sites along the two lattice diagonals again in the interior of $`\mathrm{\Lambda }`$.

Expressions for $`r_\beta (1)`$ and $`r_\beta (\sqrt{2})`$ are complicated involving elliptic integrals with different formulas for $`\beta `$ smaller, equal and larger than $`\beta _c`$ and the function $`\varphi ^1`$ must be computed numerically. Our program estima finds these estimators from a observed image.

Acknowledgments.

We thank Ricardo Maronna for comments on a draft of this paper. We thank Eric Soares da Costa and Lucas Meyer dos Santos for developing and debugging the first version of the program. This research is part of FAPESP “Projeto Temático” Grant number 90/3918-5. Partially supported by CNPq.

References

\[AK\] E. Aarts, J. Korst, Simulated Annealing and Boltzmann Machines, Wiley, New York, 1989.

\[B\] F. Barahona, On the computational complexity of Ising spin-glass models, J. of Physics A 15 (1982), 3241.

\[FFG\] P. A. Ferrari, A. Frigessi and P. Gonzaga de Sá, Fast approximate maximum a posteriori restoration of multicolor images, J. R. Statisti. Soc. B (1995), 485.

\[FP\] T. A. Ferryman and S. J. Press, A comparison of nine pixel classification algorithms, technical report, University of California, (1997).

\[FF\] L. R. Ford and D. R. Fulkerson, Flows in networks, Princeton: Princeton University Press. (1962).

\[FP\] A. Frigessi and M. Piccioni, Parameter estimation for two-dimensional Ising fields corrupted by noise, Stoch. Proc. Appl. 34 (1990) 297-311.

\[GJ\] M. R. Garey and D. S. Johnson, Computers and intractability: a guide to the theory of NP-completeness, San Francisco, Freeman (1979).

\[G\] D. Geman, Random Fields and Inverse Problems in Imaging, Lecture Notes in Mathematics 1470 (1990), Springer Verlag, New York.

\[Gi\] B. Gidas Metropolis type Monte Carlo simulation algorithm and simulated annealing, Topics in contemporary probability and its applications, 159–232, Probab. Stochastics Ser., CRC, Boca Raton, FL (1995).

\[GPS\] D. M. Greig, B. T. Porteous, A. M. Seheult, Exact maximum a posteriori estimation for binary images, J. Royal Statistical Society, Series B 51 (1989), 271-279.

\[J\] M. D. Jubb, Ph.D Thesis, unpublished, University of Bath, UK (1989).

\[M\] P. Martin, Potts models and related problems in statistical mechanics. Series on Advances in Statistical Mechanics, 5. World Scientific Publishing Co., Inc., Teaneck, NJ, 1991.

\[MW\] B. M. McCoy and T. T. Wu, The two-dimensional Ising model, Harvard Univ. Press (1973).

\[R\] D. Ruelle, Statistical mechanics: Rigorous results. W. A. Benjamin, Inc., New York-Amsterdam 1969.

\[T\] A. Trouvé Cycle decompositions and simulated annealing.  SIAM J. Control Optim. 34 (1996), no. 3, 966–986.

Instituto de Matemática e Estatística — Universidade de São Paulo Cx. Postal 66.281 — 05315-970 São Paulo SP — Brasil \<pablo@ime.usp.br\> \<gubi@ime.usp.br\> \<neves@ime.usp.br\>