PHLST WITH ADAPTIVE TILING AND ITS APPLICATION TO ANTARCTIC REMOTE SENSING IMAGE APPROXIMATION

. We propose an eﬃcient nonlinear approximation scheme using the Polyharmonic Local Sine Transform (PHLST) of Saito and Remy combined with an algorithm to tile a given image automatically and adaptively according to its local smoothness and singularities. To measure such local smoothness, we introduce the so-called local Besov indices of an image, which is based on the pointwise modulus of smoothness of the image. Such an adaptive tiling of an image is important for image approximation using PHLST because PHLST stores the corner and boundary information of each tile and consequently it is wasteful to divide a smooth region of a given image into a set of smaller tiles. We demonstrate the superiority of the proposed algorithm using Antarctic remote sensing images over the PHLST using the uniform tiling. Analysis of such images including their eﬃcient approximation and compression has gained its importance due to the global climate change.


1.
Introduction. In 2006, Saito and Remy [11] introduced a new tool for image analysis and synthesis, which is named Polyharmonic Local Sine Transform (PHLST). The PHLST resolves several problems occurring in the Local Trigonometric Transforms (LTTs) of Coifman and Meyer [5] and Malvar [9,8], such as the overlapping windows and the slopes of the bell functions. PHLST first segments (or more precisely "tiles") an image into local pieces (i.e., blocks) using the characteristic functions, then decomposes each piece into two components: the polyharmonic component and the residual. The polyharmonic component is obtained by solving the elliptic boundary value problem associated with the so-called polyharmonic equation (e.g., Laplace's equation, biharmonic equation, etc.) given the boundary values (the pixel values along the boundary created by the characteristic function).
Subsequently this component is subtracted from the original local piece to obtain the residual. Since the boundary values of the residual vanish, its Fourier sine series expansion has quickly decaying coefficients. Consequently, PHLST can distinguish intrinsic singularities in the data from the artificial discontinuities created by the local windowing. Combining this ability with the quickly decaying coefficients of the residuals, PHLST is also effective for image approximation, which was demonstrated using both synthetic and real images in [11]. Of critical importance in PHLST is how to tile an input image. As discussed in [11], there is no need to divide a smooth region of a given image into a set of smaller blocks, and in fact, that is wasteful due to the storage of the corner and boundary information in each block. In [11], however, Saito and Remy did not propose any automatic tiling algorithm for PHLST. In this paper, we will propose the PHLST equipped with an adaptive tiling algorithm for efficiently approximating given input images. We will demonstrate the effectiveness of the PHLST with the adaptive tiling algorithm for approximating Antarctic remote sensing images. Since such remote sensing images often consists of a large smooth part and smaller singular regions represented by snow ripples, fractured ice, and coastlines, our proposed PHLST algorithm is more effective to approximate such images then the PHLST with the uniform tiling as was done in [11].
We also would like to mention the importance of efficiently approximating and compressing such Antarctic remote sensing images. The current world is facing a series of unprecedented major global environmental problems caused by global warming. In the 2007 Fourth Assessment Report (AR4) by the Intergovernmental Panel on Climate Change (IPCC) of the United Nations, it is indicated that most of the observed warming over the last 50 years is likely to have been due to the increasing concentrations of greenhouse gases produced by human activities such as deforestation and burning fossil fuel. In polar regions, warming will be expected to be strongest and cause the retreat of glaciers and sea ices, even the melting of ice sheets. Partial deglaciation of the West Antarctic ice sheet could contribute 4-6 meters or more to sea level rise. This will be a big disaster for human being. A good approach to observe and analyze the change of ice structures is to compare remote sensing images of Antarctica taken at difference times. Such an endeavor forces one to store a huge amount of remote sensing data. In order to save storage space, one needs to develop a new image compression algorithm that can efficiently preserve intrinsic features (or singularities) of Antarctic remote sensing images with small storage costs. This paper is organized as follows. In Section 2, we recall the concept of PHLST. In Section 3, in order to measure the local smoothness of an image, we introduce the concept of local Besov indices and discuss the relation between local Besov indices and global Besov indices. Based on these indices, in Section 4, we derive a fundamental principle of adaptive tiling. In Section 5, we obtain a precise estimate of the nonlinear approximation error of the target function using PHLST with adaptive tiling. It is clear that the obtained nonlinear approximation order using PHLST with adaptive tiling is much better than that using PHLST with uniform tiling. In Section 6, we apply our research on nonlinear approximation using PHLST with adaptive tiling to Antarctic remote sensing image approximation.
Ω j and f j = f χ Ωj , and then split each piece f j in two components: the polyharmonic component u j and the residual v j , i.e., f j = u j + v j . Let ∆ be the Laplace operator in R d , i.e., Then the polyharmonic component u j is the solution of the following polyharmonic equation, i.e., ∆ m u j = 0 in Ω j for some m ∈ N with boundary conditions where ∂ 2l ∂ν 2l is the normal derivative at the boundary. Subtracting u j from f j , we get the residual v j . We do odd extension of v j followed by its periodic extension to extend v j from Ω j to R d . Denote the obtained function by v * j . We expand v * j into the Fourier sine series. The Fourier sine expansion coefficients of v * j decay rapidly if there is no intrinsic singularity in f j . This process is called the Polyharmonic Local Sine Transform (PHLST) [11].
In image compression and approximation, as long as the boundary data are stored and the normal derivatives at the boundary are available, the polyharmonic components can be computed quickly by utilizing the FFT-based Laplace solver developed by Averbuch, Braverman, Israeli, and Vozovoi [1,2], which we shall call the ABIV method. Moreover, the FFT-based Discrete Sine Transform is used to generate the Fourier sine coefficients of the residuals. Hence the PHLST is a fast algorithm with its computational cost O(n log n) where n is a number of pixels in an image.
When m = 1, (1) and (2) are reduced to ∆u j = 0 in Ω j and v j = 0 on ∂Ω j , the corresponding algorithm is called the Laplace Local Sine Transform (LLST). In particular, for d = 2, m = 1, f is a bivariate function and each Ω j is a square and ∆ = ∂ 2 each Ω j is a closed interval [a j , b j ], and ∂Ω j are the endpoints a j and b j . Hence u j becomes a polynomial of degree 2m − 1.
3. Besov space and Local Besov indices. We use the Besov space as the measure of the smoothness of a target function. To reflect the smoothness at each point, we introduce the concept of local Besov indices and explain the relation between the local Besov index at each point and the global Besov index. Based on these indices, in Section 4, we derive a fundamental principle of adaptive tiling.
The definition of the Besov space is based on the notion of moduli of smoothness. Let Ω be a domain in R d and a target function f ∈ L p (Ω) (1 ≤ p < ∞). Denote the For α > 0, Besov spaces are defined as where r is the smallest integer larger than α.
Let Ω be a bounded domain of R d and f ∈ L p (Ω) (1 ≤ p ≤ ∞). Now we define the global Besov index on the domain Ω.
We define local Besov indices as follows.
The following theorem explains the relations between global and local Besov indices.

Theorem 3.3.
Let Ω be a domain of R d and f ∈ L p (Ω). Then we have where α f (Ω) and α f (x) are the global Besov index and the local Besov index of f at the point x, respectively.
and λ > 0, there exists a δ(x) > 0 such that where χ E is the characteristic function of the set E. which covers the closed domain Ω.
Let Ω 1 , Ω 2 and Ω 1 ∪ Ω 2 be both domains. By the definition of the Besov space, we know that the following claim: is a domain, using the above claim together with (5), we have f ∈ B α−λ (L p (Ω)). Since λ is arbitrarily small, by Definition 3.1, it follows that On the other hand, by Definitions 3.1 and 3.2, we easily see that . From this and (6), we get 4. Adaptive tiling. In this section, using local Besov indices, we will derive a fundamental principle of adaptive tiling.
Let Ω be a domain in R 2 and the target function f ∈ L p (Ω) (2 ≤ p < ∞). Using the local Besov index of f at every point in Ω, we may adaptively segment the domain Ω and the function f .
Let ε be a given approximation error. Since f ∈ L p (Ω), by the absolute continuity of integral, for the given ε > 0, there exists a η > 0 such that for any set F ⊂ Ω whose measure |F | ≤ η, we have From this, we see that deleting a non-smooth part with small measure does not affect image approximation too much. More precisely, for the above η > 0, we can choose a set E and an index α 0 > 0 such that Below we give an adaptive tiling of the domain Ω based on the above set E We choose two sets of squares A and B step by step, where A is a set of 'good' squares and B is a set of 'bad' squares as follows.
Step 1. We divide the domain Ω into four squares {Ω j } 4 1 : where Ω j are pairwise disjoint except their boundaries. Initially, we set the "good set" A as an empty set and the " bad set " B as a set of all four squares, i.e., Let Ω j be a square in the bad set B. If Ω j does not intersect with E, i.e., Ω j ∩ E = ∅, then we call this square Ω j a "good square". We remove this good square from the bad set B and add it to the good set A. Denote If Ω j intersects with E, i.e., Ω j ∩ E = ∅, then we call Ω j is a "bad square". We retain the "bad square" in the bad set B. Denote Step 3. If the sum of the measures of squares in B is larger than η: where η is stated as above, then we return Steps 1 and 2. We divide Ω * * k into four squares {Ω * * k,l } 4 l=1 . For each square Ω * * k,l , (i) if Ω * * k,l ∩ E = ∅, then we call Ω * * k,l is a " good square ". We remove this " good square " from B and add it to A; (ii) if Ω * * k,l ∩ E = ∅, then we call Ω * * k,l is a " bad square ". We retain this " bad square " in B.
Since, in each step, B ⊃ E and |E| < η, we may continue this procedure until the sum of measures of squares in B is less than η. We denote the final set of good Q k , we have |Ω\Q| = |B| < η. By (7), Since we consider the L p −approximation of bounded functions and the values of the integrals of bounded functions on sets with small measure are small, we may delete the set Ω \ Q with measure ≤ η and the approximation error makes a slight change as in (8). Therefore, is a desired tiling of Ω.
5. L p approximation order. In this section, based on the above adaptive tiling, we derive the L p approximation orders of target functions on the domains by a combination of polyharmonic functions and sine polynomials.
Let a bivariate function f ∈ L p (Ω) (2 ≤ p < ∞) and Ω be a bounded domain in R 2 . Using the adaptive tiling in Section 4, we obtain Q = M k=1 Q k and |Ω \ Q| ≤ η.
Here {Q k } M 1 are disjoint squares. From (8), we know that The approximation of the non-smooth function f on the domain Ω is reduced to the approximation of the smooth functions f on Q.
Denote the global Besov index of f on Q k by α k = α f (Q k ) (k = 1, 2, ..., M ) and We will show that the approximation order of f is determined by the minimal index α k0 . For each Q k , using PHLST, we decompose the bivariate function f as follows where u k is the polyharmonic component satisfying ∆ m k u k = 0 on Q k with boundary conditions on ∂Q k (l = 0, ..., m k − 1).
Here we choose m k as large as possible, i.e., m k is a maximal integer satisfying 2m − 1 ≤ α k , and v k is the residual. Let Then we have Below we discuss approximations of V (x) and U (x), respectively.
Therefore, for an arbitrary small positive number s, Suppose that the center of the square Q k is θ k = (θ k 1 , θ k 2 ) and the length of its side is l k . Let the space τ consist of all sine polynomials t τ which can be expressed as where each β ν is a constant and Λ ⊂ N 2 is a set of the cardinality ≤ τ . Let σ τ (v * k ; x) be the best approximation sine polynomial of v * k in the space τ . We construct piecewise sine polynomials where n = M k=1 n k . We will choose n 1 , ..., n M such that V − V n L p (Q) ≤ ε 4 and n = n 1 + · · · + n M is small as possible. For this purpose, we choose n k such that By a known formula [4] of trigonometric approximation, we have where C k is a constant. From this, we know that (12) holds if i.e., we should choose that From this, we have n k , where n k is stated in (13). Let V (x) and V n (x) be stated in (10) and (11), respectively. Then From the construction of V n , we know that V n (x) is a piecewise sine polynomial which is determined by n coefficients. From (13), we get the approximation order of V (x) by piecewise sine polynomials V n (x). Proposition 2. Let V (x) and V n (x) be stated in (10) and (11). Then the following estimate holds: where α k0 is stated in (9) and s ′ > 0 is an arbitrarily small number.

5.2.
Approximation of U (x). By (10), we see that the approximation of U (x) is reduced to that of each u k on the square Q k . Each u k satisfies that where m k is the maximal value of m ∈ N satisfying 2(m − 1) ≤ α k . Let the four sides of the square Q k are τ 3 , and τ (k) 4 . We first consider the approximation of u k on each side τ Denote the endpoints of the side τ i , we do one-dimensional PHLST decomposition where g By (16) and (17), When we do odd extension of h (k) i followed by its periodic extension, the obtained function belongs to the periodic Besov space B α k −s * (s > 0). Denote the trigono- i , x . By a known result [4], we have We choose n ′ k such that Let u n ′ k satisfy that 1, 2, 3, 4).
(19) Here the polyharmonic function u n ′ k is determined by 8m k coefficients of polynomials and 4n ′ k coefficients of sine polynomials. Now we estimate the difference between two m k −fold polyharmonic functions u k and u n k where A k is the area of Q k . From this and (17)-(19), by the maximum modulus principle, we get From this, we have the following proposition. .
This implies the following proposition.
From this theorem, we see that U n ′ (x) is a combination of a piecewise polyharmonic function and a piecewise sine polynomial which is determined by 4n ′ coefficients of sine polynomials and 8m coefficients of polynomials where m = n k=1 m k , and V (x) is a piecewise sine polynomial which is determined by n coefficients of sine polynomials.
By Propositions 2 and 4, we get the estimates of approximation errors.
6. Numerical experiments. In this section, we will compare the approximation quality of our proposed PHLST with adaptive tiling and that with the uniform tiling using three representative Antarctic remote sensing images. To approximate an image, PHLST first tiles an image into local blocks using the characteristic functions. The minimal size of each block of the PHLST algorithm is set to either 9 × 9 or 17 × 17 [11]. The quality of image approximation is measured by PSNR (peak signal-to-noise ratio) [6,10] PSNR := 20 × log 10 max where RMSE is the absolute ℓ 2 error between the original and the approximation divided by the square root of the total number of pixels in the original image. The unit of PSNR is decibel (dB). The Erebus Ice Tongue (Figure 1, left) is a mountain outlet glacier that projects 11-12 km into McMurdo Sound from the Ross Island coastline near Cape Evans, Antarctica. It is about 10 meters high and is centered upon 77.6 degrees south latitude, 166.75 degrees east longitude. For the adaptive tiling, we first detect the points where the local Besov index attains the local minimum in Erebus Ice Tongue image. Here we do not need to compute exactly the value of local Besov index at each point. What we need is to find the location of points whose local Besov indices attain the local minima. From classic approximation theory [13,7,14], we know that the small local Besov indices correspond to the large pixel-value differences. Hence we use the popular Canny edge detector [3] to obtain these points by the output of the Canny edge detector (Figure 1,right). Using this information, we apply our adaptive tiling algorithm in Section 4 to the Erebus Ice Tongue image ( Figure 2). Finally, we approximate this image by PHLST with uniform tiling and adaptive tiling. In Figure 3, we show the quality of approximation of Erebus Ice Tongue image when we retain 0.5%-5% of the original coefficients, measured by PSNR values. It is clear that we can better approximate images by PHLST with adaptive tiling than with uniform tiling. Finally, we show the reconstructed images   Figure 4 shows reconstructed images using the top 0.5 % of the coefficients by uniform tiling and adaptive tiling.
The next image we want to examine is that of Antarctica's McMurdo Sound ( Figure 5), which is located bout 1,300 km from the South Pole. It was discovered by Captain James Clark Ross and named after Lt. Archibald McMurdo. The iceclogged waters of Antarctica's McMurdo Sound extend about 55 km long and wide. The sound opens into the Ross Sea to the north. Figure 6 shows uniform tiling and adaptive tiling of the McMurdo Sound remote sensing image. Figure 7 shows the quality of approximation when we retain 0.5%-5% of the original coefficients.
The last image we want to examine is that of the Antarctic peninsula, which is experiencing extraordinary warming. Ice mass loss on the peninsula occurred at a rate of 60 billion tonnes in 2006. Several ice shelves along the Antarctic Peninsula have retreated or disintegrated in the last two decades. Figure 8 is a satellite image of a part of the Antarctic peninsula. Using uniform tiling and adaptive tiling in Figure 9, the quality of approximation is shown in Figure 10.   Conclusion. In this paper, we proposed an efficient nonlinear image approximation scheme using PHLST combined with the automatic and adaptive tiling algorithm based on the local Besov indices of an input image. We demonstrated the importance of tiling an input image adaptively according to its regional smoothness and singularities for PHLST-based image approximation by comparing its performance with that of the uniform tiling using Antarctic remote sensing images, which often consists of a large smooth part and smaller singular regions representated by snow ripples, fractured ice, and coastlines. Finally, we would like to mention that it is important to develop an efficient compression algorithm for a practical use of our approximation algorithm in such Antarctic remote sensing image analysis, which will involve efficient quantization  and entropy coding. We will pursue this direction and hope to report our results at a later date.