A study on recovering the cloud-top height for the EUSO mission

In this paper we present some preliminary results on an optical-flow based technique aimed at recovering the cloud-top height from infra-red image sequences. This work has been carried out in the context of the development of the ”Extreme Universe Space Observatory” mission (EUSO), an ESA led international mission for the investigation of the nature and origin of Extreme Energy Cosmic Rays. The knowledge of the cloud scenario is critical to measure the primary energy and the composition of EECRs. In this work we explore the feasibility for the cloud-top height recovery, of a technique based on a robust multi-resolution optical-flow algorithm. The robustness is achieved adopting a Least Median of Squares paradigm. The algorithm has been tested on semi-synthetic data (i.e. real data that have been synthetically warped in order to have a reliable ground truth for the motion field), and on real short sequences (pairs of frames) coming from the ATSR2 data set. Since we assumed the same geometry as for the ATSR2 data, the cloud top height could be recovered from the motion field by means of the widely used Prata and Turner equation.


INTRODUCTION
The work presented in this paper has been carried out in the context of the development of the "Extreme Universe Space Observatory" mission (EUSO), that is an ESA led international mission designed to be accommodated onto the Columbus module on the International Space Station.EUSO will investigate the nature and origin of Extreme Energy Cosmic Rays.Looking downward to the Earth atmosphere EUSO will measure the UV (330-400 nm) fluorescence tracks produced by these EECRs, travelling to the ground at the the speed of light, and the diffusely reflected Cerenkov light, signature of the impact point of the shower front on the land, sea or cloud surface.In order to measure the primary energy and composition of EECRs, the knowledge of the cloud scenario is therefore critical.
Reliable estimation of cloud parameters, and in particular cloud-top height retrieval, is one of the main goals of the production of measurements from the space imagery for numerical weather forecasting, atmosphere studies and climate modelling.
Focusing on cloud-top height, different approaches are exploited to evaluate this parameter from satellite data.The most widely used procedures are based on radiative transfer methods that rely on extra information derived from radiosonde measurements or objective analyses, whose uncertainties could provide height estimates with errors as large as 1-3km, as claimed in [1].Techniques like brightness temperature (infrared window method) [2,3] for example, combine cloud Infra Red emissivity estimates, external temperature profiles or lapse rate with the cloudtop temperature measured by the 11-12µm IR-channels of satellite instruments.For optically thin clouds the observed brightness temperature needs to be modified due to the semi-transparency of the cloud.In this case in fact, the cloud temperature is affected by ground /cloud and atmosphere below.In very diffuse clouds like cirrus, the adjusted temperature may correspond closer to the cloud centre than to the cloud-top being the mass of the cloud spread out over a depth of several kilometres even though the optical depth may only be 1 or 2. Whereas none correction is applied for optically thick cloud, because all the vertically propagating IR radiation comes from the cloudtop or very close to it.Finally the cloud brightness temperature is compared to a temperature profile of the atmosphere and the lowest altitude with the same temperature is assigned as cloud height.
The CO 2 slicing method [4] uses the 15m thermal emission region where the CO 2 has a moderate to strong absorption characteristic and it is based on the ratio for two nearby bands of differences of radiation intensities measured in a pixel covered by a cloud and intensities expected on a cloudy free pixel.Finally the Oxygen A-band method [5] analyses the reflected solar radiation.
Stereoscopy is a different approach that has the advantage of depending only on the geometry of the observations.In the past stereo measurements were achieved by means of data coming from different configurations of satellites such as simultaneous geostationary satellite image pairs and also geostationary and polar orbiter image pairs [6,7].
At present there are some single polar satellites having on board instruments with two views that are currently operative: ERS2/ATSR2, ENVISAT/AATSR, EOS Terra-ASTER.Besides the single polar orbiters EOS Terra and Aqua carry an instrument (MISR) with nine views.
The Multi angle Imaging Spectro Radiometer (MISR) collects data from 4 spectral bands (Visible and NearIR) and provides 9 different views of the same scene within a total temporal range of 7 min, with an inter-camera delay of 45-60sec.It provides operational and simultaneous retrieval of cloud-top heights and motion obtained from a stereoscopic approach [8].Feature and area-based stereo-matchers are applied to suitable triplet of properly selected views to provide cloud-motion evaluation on 70.4 km grids with an accuracy of about 3m/s and cloud-height values on 1.1 km grids (due to processing time constraints) with a height resolution of 562m.
The Along-Track Scanning Radiometer-2 (ATSR2) on board the ERS2 satellite (780km altitude), collects data from 7 spectral bands (Visible and IR) and views the surface along the direction of the orbit track at an incidence angle of 55 as it flies toward the scene (forward view).Then, some 120s later, ATSR records a second observation of the scene at an angle close to 0 (nadir view) [9].Some studies on ATSR data have lead to the development of a multi-resolution feature-based matching procedure to retrieve cloud-top height maps with wind correction and results are compared in [10] to those from MISR.
In this paper we present a novel approach to the problem, that is computer vision based and it assumes that EUSO will be provided of an infra-red camera and that pictures of the Field of View will be taken at a suitable inter-frame rate.
The main problem in reconstructing a scene from images is establishing correspondences between the images [11].All the algorithms for recovering the cloud-top height from images that we have found in literature approach the correspondences problem as a stereo matching problem.Here we exploit the fact that, under the assumptions that images can be acquired at relatively high frame rate, the apparent motion in the image plane cannot be very large.Therefore correspondences can be obtained by computing this apparent motion, that is known in literature as optical flow [11].Assuming that the imaging geometry is the same as for the ATSR2 data, the reconstruction is simply obtained applying the Prata and Turner equation [12].The paper is structured as follows.Next section briefly explain the geometry of the imaging system.The algorithm is detailed in Section 3. Experiments are reported and discussed in Section 4. Section 5 is left to final remarks.

GEOMETRY OF THE SATELLITE VIEWING SYSTEM
Two views of the same scene are acquired at different time.One view is taken by the sensor while the satellite flies towards the scene at time t 1 , and a second one is acquired at a time t 2 such that the most of the same scene is in the field of view of the sensor at the two different times.The angles χ t1 and χ t2 are the satellite zenith angles projected on the sub-satellite track plane [12], due to the curvature of the Earth's surface these angles are different for each point of the image rows (across-track direction).Assuming that the effect of the wind is negligible it is easy to see, from a simple geometrical reasoning on Fig. 1, that the height h of the cloud can be recovered from the two angles χ t1 and χ t2 and the pixel shift u y along the vertical direction, that is the vertical component of the 2D motion vector on the image.By doing some simple algebra we obtain the formula h = u y 1 tan χ t1 − tan χ t2 (1)

DESCRIPTION OF THE METHOD
The scheme of the method that our system prototype follows in order to retrieve the cloud-top height is depicted in Fig. 2. In a nutshell the algorithm first establishes dense correspondences between two consecutive frames in the sequence, and then uses the disparity field together with the angle information, for recovering the height for cloud point.Points on the earth are filtered out from this computation by mean of a precomputed binary cloud mask, where all the pixels are flagged as earth points (i.e.image of point on the heart) or cloud points (images of clouds).
The system is composed by the following main modules: 1. optical flow estimation: estimates the 2D motion between two consecutive frames, say I t1 and I t2 , in the sequence; the motion field is computed both ways, i.e. from frame I t1 to frame I t2 , and from frame I t2 to frame I t1 .2. consistency check and interpolation: the two motion fields computed during the previous step are validated and merged into a single field; the optical flow for pixels for which the motion vector is not considered accurate is computed by means of interpolation.The rest of this section is dedicated to a more detailed description of the three modules of the system.

Optical flow estimation
The optical flow between the two input frames I t1 and I t2 is computed in order to establish a dense pixel correspondences map, necessary for synthesising frame I t .We implemented a non-iterative version of the Lucas-Kanade optical flow algorithm [13], where for each pixel p the motion vector u is obtained solving the system where the pixels p i are all the points in a neighbourhood of p.In order to avoid the contribution of corrupted pixels, situation that is likely to happen with this kind of images, the optical flow is robustly computed: the system is solved using a robust statistical method based on the Least Median of Squares paradigm [14,15].
If the motion vector u has a large magnitude, the optical flow algorithm cannot return accurate results, and in our case this can happen if the time interval between the two consecutive frames is large.A standard way for coping with large motion vectors is to use a multi-resolution approach [16]: the motion is computed for coarsest images and this estimation is used for predicting the solution for the finer level.Our multi-resolution optical flow, similar to the one described in [17], works as follows: 1. build a Gaussian pyramid of depth N for each one of the input images; ; this gives the residual optical flow ∆u N −1 ; 6. set u N −1 = ûN−1 + ∆u N −1 ; 7. iterate the process until the finest level of the pyramid and take u = u 1 .
The estimation of the optical flow is speeded up by computing the motion vector only for the cloud pixels.

Consistency check and interpolation
The two motion fields computed by the module described in the last section can contain corrupted vectors.These vectors should be detected and removed.A method that is straightforward and widely used, especially in stereo matching [18], is the so called consistency check.
In our case the method takes the two motion fields u t1t2 and u t2t1 , and for each pixel p it produces a new optical flow ût1t2 via the following rules where λ is a threshold value that we set to 1, meaning that the error we can tolerate in the motion field is not larger than one pixel.
The pixels that not pass the consistency test (the ones for which the motion field is set to ∞) are assigned a motion vector from the consistent vectors, by an interpolation rule.At this stage of this work we use a simple nearest point interpolation rule [19].

Height estimation
The height estimation is very simple as it is a straightforward application of Eqn.(1).The zenith angles are measured directly from the satellite, and we can assume that this measurement is accurate enough.Therefore the error in the recovered height depends only on the error on the image motion vector; in particular on the vertical component of the optical flow, as it is only u y appearing in Eqn.(1).

EXPERIMENTAL ASSESSMENT
As the only source of errors for the reconstruction of the cloud-top height is the estimated motion field, we mainly focused the experimental assessment for our method to the optical flow estimation.Moreover it must be noted that having ground truth data it is a not easy task, that we preferred to postpone at this preliminary stage of the work.

Experiments with synthetic motion
In order to evaluate the performance of the optical flow algorithm on infrared satellite images we produced synthetic sequences (with known motion field) starting from real satellite images taken from the ATSR2 data set.In practise we selected some images from the database, and from each one of them we created a synthetic sequence by applying a chain of affine transformation to the original image.More formally starting from a real image I we create a sequence I i as where H i is the 3 × 3 matrix of an affine transformation, and H 0 is the identity matrix.The images I i are then corrupted by additive zero-mean Gaussian noise.
Given a synthetic sequence we computed the optical flow between pairs of consecutive frames.As measures of goodness for the estimation of the optical flow we used the following: a) simple statistics for the error on the motion vector (mean error); b) peak signal to noise ratio between the original image and reconstructed image.
This last measure needs to be explained more in detail.Given the optical flow u i,i+1 between the frames I i and I i+1 , it is possible to warp back I i+1 using the optical flow, obtaining an image Îi , that in the noiseless case and assuming a perfect optical flow should be an exact copy of I i .Therefore we can use as measure the widely used Peak Signal to Noise Ratio (PSNR) that is defined as  Here we present results on two synthetic sequences generated from the first frame in Fig. 5.For the first one the norm of the real motion vector is one pixel for each image point.For the second sequence the true motion is three pixels for each image point.We show for each sequence both the error measures.The PSNR is shown for all the three optical flow computed by our algorithm (the two ways optical flow and the one obtained from the consistency and interpolation step).
For both the sequences we show the evolution of the error measures across the sequence.In Fig. 3 we show the results for the first sequence.The results for the second sequence are shown in Fig. 4. The experiments where run with a three levels pyramid, and the size of the neighbourhood for solving for the motion vector is 9 × 9.
The results reported here in the figures show that the error in recovering the motion field is, on average, always smaller than one pixel, that is a very conforting result.The graphs on the PSNR show that the values relative to the optical flow after the consistency test are, in general, between the values relative to the other two motion fields: this means that the optical flow computed by our routine can be considered as a mean between the two motion fields that can be obtained directly from the images.

Experiments on ATSR2 data
We also run some experiments on pairs of real images, that we extracted from the ATSR2 data set.For these pairs we could run a complete recovery of the cloud-top height.However we must admit that this can of data are not very suitable for our algorithm, as the time interval between the two images is long (120 secs), and therefore the motion vector can sometimes be very large.Moreover the two images are originally acquired at different resolutions, and then the smallest one is rescaled to the same size of the largest.This is a problem for stereo matching algorithm, but is particularly serious for a differential approach like ours.
In Fig. 5 we show a pair of images from the dataset, and in Fig. 6 the relative height map is shown.Darkest pixels represents points closest to the ground level.

CONCLUSIONS
In this paper we presented some preliminary results of an algorithm for recovering the cloud-top height from infrared satellite images.The algorithm assumes that the images can be acquired at a relatively high frame rate (i.e., norm of the motion field no more that 3-4 pixels).The experiments run on synthetic sequences returned promising results, as shown and discussed in Section 4. The results on real data are less exciting, but this is due more to the nature of the particular data set that we used, that makes them not very suitable for our algorithm.
More work is planned.First we need to evaluate the performance of the whole algorithm from real data, getting a ground truth from other kind of sensors.A comparison between the results of our algorithm and the ones obtained from already existing techniques needs to be done.Moreover we would like to explore the use of state of the art stereo matching algorithm in this context.

FIGURE 1 :
FIGURE 1:The geometry we assume for the imaging system

FIGURE 2 :
FIGURE 2: Schematic representation of the algorithm

2 . 1 t1and I N − 1 t2; 4 .; 5 .
compute the optical flow u N between the coarsest images I N t1 and I N t2 ; 3. scale and interpolate u N to obtain a prediction ûN−1 of the flow between I N −warp I N −1 t1 with the flow ûN−1 obtaining the image ÎN−1 t1 compute the optical flow between images ûN−1 and I N −1 t2

2 ,FIGURE 3 :
FIGURE 3: Results of the experiments on the first synthetic sequence.Left: mean error on the OF.Right: PSNR for the reconstructed images

FIGURE 4 :
FIGURE 4: Results of the experiments on the second synthetic sequence.Left: mean error on the OF.Right: PSNR for the reconstructed images

FIGURE 5 :
FIGURE 5: Real pair from the ATSR2 data set

FIGURE 6 :
FIGURE 6: Left: cloudiness binary mask for the real pair in Fig. 5. Right: height map for the same real pair, brightness is proportional to the height