EURASIP Member Hayat Khizar

EURASIP Member William Puech
email: william.puech@lirmm.fr

Naveed Islam

Gilles Gesquière

Seamless Heterogeneous 3D Tessellation via DWT Domain Smoothing and Mosaicking

With todays geobrowsers, the tessellations are far from being smooth due to a variety of reasons: the principal being the light difference and resolution heterogeneity. Whilst the former has been extensively dealt with in the literature through classic mosaicking techniques, the latter has got little attention. We focus on this latter aspect and present two DWT domain methods to seamlessly stitch tiles of heterogeneous resolutions. The first method is local in that each of the tiles that constitute the view, is subjected to one of the three context-based smoothing functions proposed for horizontal, vertical, and radial smoothing, depending on its localization in the tessellation. These functions are applied at the DWT subband level and followed by an inverse DWT to give a smoothened tile. In the second method, though we assume the same tessellation scenario, the view field is thought to be of a sliding window which may contain parts of the tiles from the heterogeneous tessellation. The window is refined in the DWT domain through mosaicking and smoothing followed by a global inverse DWT. Rather than the traditional sense, the mosaicking employed over here targets the heterogeneous resolution. Perceptually, this second method has shown better results than the first one. The methods have been successfully applied to practical examples of both the texture and its corresponding DEM for seamless 3D terrain visualization.

Introduction

The contemporary geobrowsers-like Google Earth, NASA's World Wind, or Microsoft's Virtual Earth-give little importance to individual client characteristics, and with the lowend clients, the performance is far from being satisfactory. In addition, area coverage by these browsers is slack, if not judicious, in the case of the underdeveloped and developing parts of the world as opposed to their developed counterparts. There is a need to reform all these browsers to cater for the diversity of clients through some scalable data structuring. One obvious choice is to create multiple levels of detail (LOD) by representing the shape at different levels of approximation. This is essential for the terrain tessellation adjustment as a function of the view parameters [START_REF] Losasso | Geometry clipmaps: terrain rendering using nested regular grids[END_REF]. For LOD, one can rely on the multiresolution nature of the discrete wavelet transform (DWT). It is better to employ some standard state-of-the-art DWT, like the now widely accepted JPEG2000 standard (The ISO/IEC 15444-1 standard.). In fact JPEG2000 is better than the popular pyramidal resolution representation with JPEG or other format closer to normal graphic cards that are of high storage cost: there is a need to store different resolutions. Additional advantages of the JPEG2000, in this context, include capabilities like better compression and progressive data transfer.

Traditionally, tile-based approaches have been favored to render large terrains. The problem, however, arises when the tiles have either been rendered at different resolutions or not photographed in identical luminance or radiometric conditions. To elaborate this we rely on the two snapshot examples of Google Earth, shown in Figure 1. In the first example (Figure 1(a)), where the observer is changing his view to the right, it can be seen that the resolution is not uniform. Usually there are two reasons for this heterogeneity. First is the download latency and second is the nonavailability of proper resolution due to a variety of reasons like area security, data inaccessibility, and the like.

The second snapshot (Figure 1(b)), which pertains to the area surrounding the Pakistani capital of Islamabad, suffers from both heterogeneous resolution and aerial photography in nonidentical light conditions. Obviously the photographs have been taken at different times and then joined together to what we see in the snapshot.

In essence, one can safely deduce that even for today's geobrowsers the tessellations are far from being smooth due to varying resolutions or heterogeneous luminance conditions or even both when the tiles have been acquired. Whilst, the traditional focus of the mosaicking community has been the seamless rendering of tiles that may vary in luminance, little attention has been paid to the heterogeneity due to resolution. That is why traditional mosaicking may correct the light difference phenomenon but the underlying tiles may still be hardly seamless if they vary in resolution.

In this paper we focus on the resolution part, mainly, as we try to smoothen the tile interfaces to avoid the artifacts like those in Figure 1, for example. Our methods are aimed at rendering texture and DEM tessellations, from JPEG2000coded tiles, in a scenario where only visible tiles are rendered and the focused ones, out of these, have the highest quality depending on the viewer's distance. As the viewer moves in, the tiles corresponding to the new view must be rendered and gradually refined in quality as the focus is coming near to it. This refinement is additive, as more and more subbands are added to it, at runtime. The main challenge to this strategy is that the tile must be seamlessly stitched for both the texture and the DEM and the boundaries must be diluted enough to prevent popping artifacts. In this paper, although we are more inclined to the texture, the DEM has also been treated as its integral part. We, therefore, present two different methods to look for a seamless tessellation that is composed of heterogeneous resolution tiles. In fact this work is a continuation of our earlier efforts about the synchronous unification of the DEM/texture pair through data hiding [START_REF] Hayat | Scalable 3D terrain visualization through reversible JPEG2000-based blind data hiding[END_REF] for scalable 3D visualization. That work is our inspiration for the LOD aspect, since its scalable character can result in various tile qualities for both the DEM and its texture.

The first method is local in nature and relies on contextbased smoothing of the tile boundaries. The smoothing is carried out at the subband level in the DWT domain by employing three different smoothing masks. We assume that, at worst, there can be three possible tile resolutions, namely, high, middle, and low. The lowest-resolution tiles of the panorama are not treated while the highest-resolution tile are smoothed by a radial smoothing function. The rest of the tiles are rendered at mid-level resolution. Of these the left/right neighbors of center are smoothed horizontally, while the above/below neighbors are smoothed vertically. The treated tiles can be joined, after each has been subjected to an inverse DWT, to get a smoother tessellation.

In the second method, we take the change in viewpoint focus, analogous to a window-sliding approach. The worst we expect is that the window may come up to be composed of three different tile qualities with the resultant artifacts at tile interfaces. To dilute these artifacts, we treat the tiles at the subband level, in the DWT domain, by employing operations involving suitable subband-sized composite masks, conceived with smoothing and mosaicking in perspective. The resultant composite subbands are subjected to a global inverse DWT to get the final seamless tessellation.

The rest of the paper is arranged as follows. Section 2 gives a brief literature survey on the subject. The presentation of the proposed methods is the theme of Section 3. Application of our methods to a practical worst case example is being demonstrated in Section 4. Section 5 concludes the paper.

State of the Art

The process of mosaicking involves the connection of two or more images to get a new wide-area image with no visible seams. In the case of aerial or satellite photographs, a wide screen-view is synthesized by the proportional combination of a set of images. The mosaic is hence commonly used to increase the view field with pasting successively overlapped images onto the same plane or other structures, like a cylinder or a sphere. In order to be efficient, a mosaicking method should be robust to light change, moving objects, image noise and to some extent rotation and zooming [START_REF] Zagrouba | An efficient imagemosaicing method based on multifeature matching[END_REF]. An integral part of the process is image registration which pertains to the overlaying two or more images of the same scene taken at different times, from different viewpoints or by different sensors. The final information is gained from variety of sources like image fusion, change detection and multichannel image restoration. In this context Zitova and Flusser [START_REF] Zitova | Image registration methods: a survey[END_REF] classify the various approaches according to their nature (area-based and feature-based) and to four basic steps of image registration procedure, namely, feature detection, feature matching, mapping function design and image transformation and resampling. The thrust of the mosaicking methods, in general, is either on the interest point matching or on theoretical corner models. All these, hardly consider the aspect where a resolution image is to be stitched to a considerably lower-resolution image.

Terrain LOD algorithms can be classified [START_REF] Losasso | Geometry clipmaps: terrain rendering using nested regular grids[END_REF], on the basis of the hierarchical structure, into four groups, namely, the triangulated irregular networks (TINs), bin-tree hierarchies, bin-tree regions and tiled blocks. The tessellation, in the last category, is carried out with the square tiles of different resolutions as has been done by Wagner [START_REF] Wagner | ShaderX2: Shader Programming Tips and Tricks with DirectX 9.0, chapter Terrain Geomorphing in Vertex Shader[END_REF]. But the process requires seamless stitching at the boundaries. According to Deb et al. [START_REF] Deb | Realtime streaming and rendering of terrains[END_REF], Wagner's technique is not suited for the terrains with large height variations and that is why they tweak the technique by having the projection on realtime average height of terrain for a client/server terrain streaming to handle heterogeneous clients. The blending factors are calculated on a per-tile basis because of the use of a regular tile structure, thus, reducing the amount of computation. Special techniques based on mipmaps and clipmaps [START_REF] Tanner | The clipmap: a virtual mipmap[END_REF] can be found in the literature. Losasso and Hoppe [START_REF] Losasso | Geometry clipmaps: terrain rendering using nested regular grids[END_REF] break the terrains into geometric clipmaps of varying metric sizes which can be used as LODs through a view-centered hierarchy resulting in a simplified spatial and temporal interlevel continuity. A thresholding scheme based on calculated ground sample distances is proposed by Tsai and Chiu [START_REF] Tsai | Adaptive level of detail for large terrain visualization[END_REF] with a nested LOD system for efficient seamless visualization of large datasets. In order to remove the T-junctions they have put forward a mesh refinement algorithm. Ueng and Chuang [START_REF] Ueng | Terrain rendering with view-dependent LOD caching[END_REF] propose a dike structure between two adjacent blocks, for a smooth blending between two meshes of different LOD. This approach is combined by Li et al. [START_REF] Li | Visualization for HPC data-large terrain model[END_REF] with heuristics to dynamically stitch the tile meshes together seamlessly. For the texture part they use high-resolution aerial photos subjected to the similar LOD mechanism as described for the meshes. Larsen and Christensen [START_REF] Larsen | Real-time terrain rendering using smooth hardware optimized level of detail[END_REF] try to avoid popping for low-end users by exploiting the low-level hardware programmability in order to maintain interactive frame-rates. They claim their work as pioneering as far as a smooth LOD implementation in commodity hardware is concerned.

Little attention has been paid to the texture images which require their own LOD structure. There have been proposals, like texture clipmaps [START_REF] Tanner | The clipmap: a virtual mipmap[END_REF], but texture tiling is usually preferred. The multiresolution technique, described by Döllner et al. [START_REF] Döllner | Texturing techniques for terrain visualization[END_REF], presents static LOD terrain models, allowing for combination of multiple large-scale textures of different size. The view-dependent texture management technique of Okamoto et al. [START_REF] Okamoto | Texture management in view dependent application for large 3D terrain visualization[END_REF] manages multiresolution textures in multiple caching levels-database, main memory and graphics card memory-which the authors claim to be suitable for hardware with limited resources. The approach of Buchholz and Döllner [START_REF] Buchholz | View-dependent rendering of multiresolution texture-atlases[END_REF] includes dynamic texture loading for memory management and, combine several textures into atlases for avoiding too much texture switches [START_REF] Wloka | ShaderX3: Advanced Rendering with DirectX and OpenGL, chapter Improved Batching via Texture Atlases[END_REF]. Their computationally complex preprocessing part helps them to significantly reduce runtime rendering overhead. Frueh et al. [START_REF] Frueh | Automated texture mapping of 3D city models with oblique aerial imagery[END_REF] described an approach to create a specialized texture atlas for building facades and supports efficient rendering for virtual flythroughs. The created atlas is static, and different texture resolutions are not considered. Based on the minimum error boundary cut idea of the image quilting algorithm as well as Efros and Freeman [START_REF] Efros | Image quilting for texture synthesis and transfer[END_REF], Somol and Haindl [START_REF] Somol | Novel path search algorithm for image stitching and advanced texture tiling[END_REF] propose their minimum error path search algorithm. To avoid visible seams, for combining incompatible pieces of texture, they propose to find (possibly irregular) boundaries between the image pieces to minimize the visual error. Based on eccentric image overlapping, the source image is cropped at the prospective interface along the minimum error path and placed over the target background image. The problem is, however, if no good path exists in the error map, visible artifacts may be inevitable. Integer wavelet transform has been applied to DEM and a simple coding algorithm with high efficiency was introduced by Chen and Li [START_REF] Chen | DEM compression based on integer wavelet transform[END_REF]. Many references can be found on the subject in the works of Danovaro et al. [START_REF] Danovaro | Level-of-detail for data analysis and exploration: a historical overview and some new perspectives[END_REF] and Buchholz and Döllner [START_REF] Buchholz | View-dependent rendering of multiresolution texture-atlases[END_REF]. A survey about 3D interactive rendering by Pajarola and Gobbetti [START_REF] Pajarola | Survey of semi-regular multiresolution models for interactive terrain rendering[END_REF] is worth reading.

The Smoothing Strategies

In this section we explain the two methods for seamlessly stitching heterogeneous resolution tiles. This heterogeneity is traced back to the concept of the level of detail (LOD). For the LOD we are relying on the multiresolution nature of DWT from the JPEG2000 codec. The reasons for choosing JPEG2000 are its scalability, its wide acceptance as a standard, its support to various wavelet forms for multiresolution, progressive transfer capability and last but not the least is its approximation paradigm for LOD. The various approximations do not differ in size but only in quality which is additive in nature. If each tile is DWTed in the JPEG2000 codec at level L then a level-l (≤ L) approximate image is the one that is constructed with (1/4 l )×100 percent of the total coefficients. These coefficients correspond to the available lowest 3(L - l) + 1 subbands and for the coefficients corresponding to the rest of the 3l subbands, zeros are stuffed. Subsequent application of L-level inverse DWT yields what is known as the l-approximate image. For example, level-0 approximate image is constructed from all the coefficients and level-2 approximate image is constructed from 6.12% of the count of the initial coefficients.

3.1. Scenario. We are assuming a simplistic scenario of a panorama of 3 × 3 tiles given in Figure 2. This use case should represent almost all the possibilities we will have in our application. An individual tile is large enough to circumscribe the field of view if the observer is close enough to require highest visual resolution. If the observer is far off, his view is limited to at least 3 × 3 tessellation with all the tiles at their lowest visual resolutions but as soon as he draws closer, further details may be added to the tile of focus and, to some extent, its 4-neighbors to improve their visual resolution. This leads to variations in the tile qualities of the view tessellation. The tessellation is thus heterogeneous in terms of resolution and a given tile can have one of the three qualities, namely, high, middle and low. The highquality resolution (H) corresponds to the center which is the viewpoint focus for a given distance to the observer in a flyby process. We are rendering the 4-neighbors of the center with the middle-resolution (M) tiles and the rest, that is, the corners, with low resolution (L). The terms high, low and middle are not static and their extent depends on the distance of observer from the view as well as the time passed after the viewer changed his position. It also depends on the capacities of the media employed, as we may have M and L or just L for poorer conditions. There may be a tessellation having concentric square rings with each successive ring having a different quality [START_REF] Tanner | The clipmap: a virtual mipmap[END_REF] but we believe that at any instant one would hardly be dealing with a rendering involving more than three resolutions.

3.2. Shape of the Smoothing Function. Figure 3 depicts the shape of an example smoothing function in the simplest of forms. Note that we are affecting the quality of only a small part on the periphery, of the higher resolution tiles in the view, during a given smoothing process. For any given subband, the smoothing function creates a scalar multiplication mask (over [0, 1]) of the size of the subband. This depends on the position of the tile, to which the subband belongs, in the view. for example, for the subbands of a tile having low resolution neighbors on left/right, we may consider the central column of coefficients as the origin in the figure and x may then be the offset from this central column on right or left with w being half the width of the subband. In a similar fashion, for a tile with up/down neighbors in the view, the origin may be the center row with row offset represented by x. Finally for a tile in the 8-neighbor environment, the center of its subband may be the origin and x may be the offset of the polygonal ring of coefficients. The function f (x) has two parameters, a and b. the parameter a decides where to start the smoothing process in the subband whereas b denotes the index where it should fall to zero. To have a full smoothing in a subband we must have a = 0 and b = w, but to dilute the seam, it is preferable to have b < w and to not deteriorate too much the subband it is preferable to have b ≈ w. To not have a strong smoothing between a and b and avoiding the full deterioration of the subband, one can locate a at one quarter of the subband, that is, a = w/2, where, as already stated, w is roughly half the subband dimension (width, length or radius).

In Figure 3 the given function is of the form,

f (x) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1, if x ∈ [0, a), 1 - x -a b -a , if x ∈ [a, b],
0, elsewhere.

(

) 1
Of special importance is the shape of the fall from a to b which may be streamlined to emulate a more sophisticated function, for example, tail of a Gaussian distribution. One such function could be:

f (x) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1, if x ∈ [0, a), e -(x-a) 2 /2σ 2 , if x ∈ [a, b],
0, elsewhere.

(

Note that beyond a we shall be dealing with b -a consecutive integers as x; therefore σ 2 = ((b-a) 2 -1)/12. One can borrow from sinusoidal functions too and a smoothing function that we tested is given below:

f (x) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1, if x ∈ [0, a), 1 2 1 + cos x -a b -a π , if x ∈ [a, b],
0, elsewhere.

(

There can be many more functions to emulate but due to space limitation we would restrict to the three given in (1)-(3) only.

Let the high, the middle, and the low resolution of our scenario correspond to l h , l m , and l l approximate images, respectively. The quality gap would essentially bring with it the popping artifacts at tile boundaries. To undo that, we are proposing two different DWT domain smoothing methods in the next two sections.

Method I:

Smoothing the Tiles before Stitching. In this section we present a wavelet-domain context-dependent smoothing strategy illustrated in Figure 4. At the corners (labeled L) are the l l -approximate images obtained from lowest 3(L -l l ) + 1 subbands with the rest of 3l l subbands being replaced by 0s. For the rest of the tiles one can partition their subbands into the following three sets. (ii) The highest 3l h subbands of center tile (H) and 3l m subbands of its 4-neighbors (labeled M) are stuffed with 0s. (iii) Last are the remaining 3(l l -l h ) and 3(l l -l m ) subbands of central tile and its 4-neighbors, respectively.

It is the third set (or its subset thereof) which will be subjected to the DWT domain treatment with three different smoothing functions on each subband individually. These smoothing functions are the following.

(1) The medium-resolution left/right neighbors of center require the horizontal smoothing function.

(2) The medium-resolution above/below neighbors of center need the vertical smoothing function.

(3) The high-resolution center is treated through the radial smoothing function.

The horizontal smoothing mask would have a central band of rows of 1's. The thickness of this band depends on the maximum diameter, δ, of the view focus. If the subband belongs to the kth resolution level, then the thickness should be around δ/2 k . Beyond the band of 1's, the rows gradually attenuate on both sides approaching 0 at the periphery, complying anyone of the (1)-( 3). An example smoothing function, for the horizontal case and loosely based on (1), is described in the appendix in the form of Algorithm 1. Along similar lines we can get a vertical smoothing mask by transposing the horizontal mask for a given subband. An example illustration of horizontal smoothing between two tiles (a high-resolution tile H and a medium-resolution tile M on right) is given in Figure 5.

The radial smoothing mask is a bit tricky. We would have a center core of 1's having a diameter of about δ/2 k for a the subband of the k th resolution level. Ideally, this core should be circular with a being its radius, but due to the discrete character of the subband, the shape is at best polygonal. For the sake of simplicity of Algorithm 2 presented in the appendix, we are considering a square core of 1's with dimensions 2a. Each ring of coefficients around the core would gradually attenuate approaching 0 at the corners (or the bth ring to be precise) and around 0.5 in the up/down and left/right directions.

After the smoothing treatment of a subset or all of the third group of subbands and inverse DWT, one can have a seamless visualization of the type given in Figure 4. In our application, we are using the same strategy for the DEM to avoid popping effect on elevation data.

Method II: Sliding Window Analogy to View Field.

For method II we introduce a sliding window analogy, as shown in Figure 6, where the red square/window represent our initial point of focus which had been rendered with highquality resolution, H. Let the focus changes and the window slides to the position represented by the blue square. The new window may have, at the worst, three different resolutions. This situation is illustrated by the 2×2 panorama in the right part of Figure 6. Rendering three different resolutions, side by side, will obviously give rise to nonsmooth transitions at the interface of heterogeneous tiles.

For a 1D viewpoint displacement from H to M or M to L, whether horizontal or vertical, a smooth transition can be realized by following the scheme given in Figure 7. This figure illustrates the DWT domain, H to M seamless smoothing, in the horizontal direction. For H we have all the subbands available but for the M part we have just level-2 subbands available, and six highest-frequency subbands are missing which will be replaced with zero coefficients. This gives rise to two types of combinations. First, the direct concatenation of nonzero subbands with their counterparts from H. Second, the concatenation of each of the zero subbands with the pretreated corresponding subband of H. The pretreatment involves scalar multiplication with a subband-sized mask generated by one of the functions, f (x), outlined Section 3.2, with x being the column index of the 5) for j ← b to a do [START_REF] Deb | Realtime streaming and rendering of terrains[END_REF] min coeff ← min coeff + 2/(h -2a) [START_REF] Tanner | The clipmap: a virtual mipmap[END_REF] for i ← 0 to w do [START_REF] Tsai | Adaptive level of detail for large terrain visualization[END_REF] B[i, j] ← min coeff [START_REF] Ueng | Terrain rendering with view-dependent LOD caching[END_REF] B[i, h -1j] ← min coeff [START_REF] Li | Visualization for HPC data-large terrain model[END_REF] end for [START_REF] Larsen | Real-time terrain rendering using smooth hardware optimized level of detail[END_REF] end for [START_REF] Döllner | Texturing techniques for terrain visualization[END_REF] end Algorithm 1: Example horizontal smoothing function.

Input: dimensions (w, h) of the DWT domain subband and smoothing parameters (a, b) 3) declare an array center row [w] with all elements initialized to 1 (4) for i ← (w/2 -1a) to 0 do [START_REF] Wagner | ShaderX2: Shader Programming Tips and Tricks with DirectX 9.0, chapter Terrain Geomorphing in Vertex Shader[END_REF] center row 7) end for (8) declare a w × h array B (9) for j ← w/2 -1 to 0 do [START_REF] Li | Visualization for HPC data-large terrain model[END_REF] for i ← w -1 to 0 do [START_REF] Larsen | Real-time terrain rendering using smooth hardware optimized level of detail[END_REF] if j ≥ w/2a then (12)

Output: Radial smoothing square mask B of dimensions (w, h) (1) begin (2) set interval ← |b -2/(w + h -4a)| (
[i] ← center row[i + 1] -interval (6) center row[w -1 -i] ← center row[w -2 -i] -interval (
B[i, j] ← center row[i] (13) B[i, w -1 -j] ← center row[i] (14) else (15) B[i, j] ← B[i, j + 1] -interval (16) B[i, w -1 -j] ← B[i, w -j -2] -interval (17)
end if [START_REF] Somol | Novel path search algorithm for image stitching and advanced texture tiling[END_REF] end for [START_REF] Chen | DEM compression based on integer wavelet transform[END_REF] subband coefficient. Note that for a vertical transition, column index will have to be replaced with the row index in the function f (x). The resultant concatenated H/M subbands, upon inverse DWT, results the required seamlessly smoothed image.

Note that we are affecting the quality of only a small part on the periphery of the tile during a given smoothing process. For any given subband, the smoothing function creates a scalar multiplication mask (over [0, 1]) of the size of the subband. For example, the H to M horizontal smoothing mask would have a band of columns of 1's on the left. The thickness of this band depends on the maximum diameter of the view focus. Beyond the band of 1's the columns gradually attenuate toward right, approaching 0 at the periphery. Along similar lines we can get a vertical smoothing mask for a given subband.

In two dimensions, the diagonal transition from H to L is also to be taken into account. One can think of the 2D version as an operation involving either two horizontal 1D smoothings followed by a vertical 1D smoothing, as shown in Figure 8, or two vertical 1D smoothings followed by a horizontal smoothing. Keeping this in view, these operations can be easily carried out in parallel to avoid time delay, thanks to the involvement of subbands and the independence of blocks as envisioned in the JPEG2000 standard. Taking the argument further, since we have three different resolution tiles, for concatenation we may face three subband combinations as follows.

(1) The corresponding subbands of all of the H, M, and L are nonzero implying direct concatenation according to their localization in space. (2) With zero subbands of L corresponding to nonzero subbands of M and H, we would follow the smoothing strategy described in Figure 8.

(3) For nonzero subbands of H, corresponding to zero subbands of L and M, we would concatenate the former in the pretreated form with the zero subbands keeping in view their localization in space.

But we believe that if one think at the outset about the diagonal transition by starting the smoothing process from those subbands of M and H whose counterparts are zero in L, the overall smoothing may still be commendable. In addition the closest part of H to L can be programmed to be degraded more than the rest of its interfaces during smoothing.

Simulation Results

We have applied our smoothing strategies to practical examples provided by IGN (http://www.ign.fr/) France and the results have been encouraging. For the reader we are presenting a representative example of 3072 × 3072 texture tessellation composed of 9 tiles of size 1024 × 1024 pixels. The corresponding DEM tessellation at our disposal has a size of 96 × 96 with each of the nine tiles having 32 × 32 coefficients. The texture tiles were transformed at level-4 reversible DWT in a standard JPEG2000 encoder whereas the DEM tiles were subjected to reversible level-4 DWT simply. In our example visualization, we are going for a case where the corners are being rendered with level-3 approximate tiles which would mean a rendering with only 1.56% of the coefficients. The 4-neighbors of center have the middle quality of level-2 approximation. For the centre we are taking highest possible quality, that is, level-0. The high-quality gap between three qualities is likely to make the tile boundaries conspicuous. For the purpose of clarity we are taking subset 2048 × 2048 tessellation comprising of the tile H and its right, below and base diagonal neighbors-eastern, southern and south eastern, to be exact. A texture tessellation with these tiles without any smoothing treatment is shown in Figure 9(a). To observe the boundaries, we have magnified a small 384×384 area at coordinates (768, 768), corresponding to the small white square in Figure 9(a), from the tessellation in Figure 9(b). We have deliberately selected a region that contains each of the three (H, M, and L) resolutions and smooth quality. This magnification is only for the purpose of comparison and the magnified portion does not represent the visual effect, we are up to, since the distance, it represents, warrants improved qualities at corners in our strategy. For the treated examples, we will not be showing the whole tessellation but rather displaying the corresponding square region in each tessellation for comparison.

Let d imply that 3d + 1 lowest subbands of level-4 DWT decomposition have been avoided during smoothing, for example, d = 1 implies that the lowest 4 subbands have not been touched. The tessellations for various smoothing possibilities are illustrated in Figure 10. Based on d we are defining d h which corresponds to the center tile and d m that corresponds to the 4-neighbors. Note that since the 4-neighbors are level-2 approximates, the highest six subbands are already zero. This implies that d m cannot exceed 1. If we do not tamper the lowest four subbands during the horizontal smoothing of left/right neighbors and vertical smoothing of above/below neighbors, then d m = 1, otherwise d m = 0 where the only lowest level-4 subband is avoided from smoothing. d h , which corresponds to the center H, can assume a value of up to 3 and it was observed that the border between the four tiles gradually smoothen as we decrease the value of d h since more and more subbands are being involved in the radial smoothing of center. Since we are working on a worst case scenario and we are closer than we should be, the borders do not disappear at all but in the original environment the transition is seamless. When subjected to method II, the resultant tessellation was considerably smooth as can be seen in Figure 11 which represents the detailed square area in the tessellation corresponding to Figure 9(a). The difference is obvious as it is very hard now to differentiate the tile interfaces, even though they have different qualities. The better results of method II are also supported by the PSNR values. The original texture detail (Figure 9(b)) has a PSNR of 23.44 dB with respect to the original full-resolution panorama. With method I this value comes down to 21.40 dB for the case where d m = 1, d h = 1 (Figure 10(c)) as against 22.09 dB for method II (Figure 11). The corresponding difference images are shown in Figure 12.

We now apply, to our example, the two other smoothing functions (Equations ( 2) and ( 3)) in liaison with method II. The results are shown in Figure 13. It is evident that the blurring is a bit in excess and the seam is hardly evident. But this has degraded the quality in the vicinity the tile interfaces. This can be judged from the PSNR of these images (20.73 dB and 20.75 dB) as well as their respective differences shown in Figure 14. These results imply that although seam elimination is marginally, the accompanying blur is excessive and the linear form behaves well in this regard. This does not mean that one should entirely ignore the functions. In fact, they are equally effective but further detailed analysis is needed to make them or similar function indispensable. The purpose of this work is to show the way to eliminate the seams using some suitable function tailored to the view environment. The quest for alternatives is open to the user. Therefore, from this point onwards we will be employing the simplest of the potential functions (i.e., the linear one of ( 1)) for smoothing.

The DEM data corresponding to the above example is too small and using the same resolutions will not be a good idea, for example, a level-3 approximate tile would require only four coefficients. Still, for the purpose of comparison we are illustrating a 96 × 96 tessellation along the similar lines as the texture in Figure 15. The original, with all tiles at highest-resolution, is given in Figure 15(a). For the sake of comparison, Figure 15(b) shows a tessellation obtained by the simple replication of available coefficients to fill in for the missing coefficients to get the original size. For example, if only level-2 subband (1/16th of the original in size) of M is available, each of its coefficient has to be quadruplicated horizontally as well as vertically to get the original size and fill One observes that the difference is not that evident. The same is not the case with method II, since the results are far better, as can be seen in Figure 17. Method II has thus an obvious advantage, as far as the DEM stitching is concerned and the results are better than those of method I.

When it comes to 3D, the DEM triangulation results support the arguments given above. For comparison we are showing the 3D triangulation with and without the overlaying of the corresponding texture in Figure 18. As an additional reference, the replication example of Figure 15(b) is shown in the triangulated form in Figure 19. The popping artifacts are conspicuous and obviously this is not we are up to. Things are a bit improved when the level-0, level-2 and level-3 approximation images are tessellated without any treatment, as shown in Figure 20. The artifacts are still visible due to the tile seams, however. Figure 21 shows the 3D view rendered from the heterogeneous tessellation after treatment with method I. The results are not that impressive, as far as the DEM triangulation is concerned. This may probably due to the use of very small tiles, as the 32×32 tiles at our disposal may not allow for effective smoothing. But even with this situation, method II gives interesting results, as illustrated in Figure 22. It is evident that the artifacts have been eliminated to a great extent and in fact there are hardly any feeling of unsmooth transition.

Conclusion

The results, presented above, have been interesting in the face of the fact that we presented a worst case example tessellation in terms of quality difference as well as DEM size. Whilst the two methods were at par in the case of texture stitching, method II had the edge in the case of DEM. The problem of synchronization may arise if the DEM is different enough to warrant its own LOD. One probable answer to this is the fact that since the LOD is wavelet driven, one can always synchronously unify the DEM texture pairs as explained by Hayat et al. [START_REF] Hayat | Scalable 3D terrain visualization through reversible JPEG2000-based blind data hiding[END_REF] or even adapt the synchronization [START_REF] Hayat | Scalable data hiding for online textured 3D terrain visualization[END_REF] to exclude some subbands from embedding. One can still argue that the seam elimination has been realized at the expense of blurring the better resolution tile. The counter argument to this is the fact that the resultant blur is temporary and as soon as further resolution information for a given tile arrives, one can employ complement of the mask to restore the lost information.

In our method we are using a regular grid. To increase the performance, we would like to employ an irregular grid depending on the DEM resolution. The well-known Tcracks problem will be avoided by the application of the junction technique of Pajarola and Gobbetti [START_REF] Pajarola | Survey of semi-regular multiresolution models for interactive terrain rendering[END_REF]. In the continuation of our work we would like to extensively analyze and optimize the parameters of the functions employed for smoothing. Alongside, a better function needs to be discovered to further improve the stitching. The overhead in terms of the DWT and its inverse has to be quantified and optimization must be proposed for practical applications. The ultimate is its extension to real-time environment, with the emphasis on scalability.

Figure 1 :

 1 Figure 1: Google Earth examples.

Figure 2 :Figure 3 :

 23 Figure 2: The tessellation scenario.

Figure 4 :

 4 Figure 4: Method I as applied to the tessellation scenario.

Figure 5 :Figure 6 :

 56 Figure 5: Example horizontal smoothing strategy in one direction (right).

Figure 7 :

 7 Figure 7: The proposed smoothing strategy in 1D.

Figure 8 :

 8 Figure 8: The proposed smoothing strategy as applied in 2D.

Figure 9 : 1 Figure 10 :

 9110 Figure 9: Example texture tessellation: l h = 0, l m = 2, l l = 3.

Algorithm 2 :

 2 Example radial smoothing function.

Figure 11 :

 11 Figure 11: Application of method II to the example: l h = 0, l m = 2, l l = 3.

Figure 12 :

 12 Figure 12: Difference images with the original full-resolution tessellation detail.

Figure 13 :Figure 14 :

 1314 Figure 13: Method II as applied special smoothing functions: l h = 0, l m = 2, l l = 3.

2 Figure 15 :

 215 Figure 15: Example DEM tessellation without smoothing (l h = 0).

Figure 16 :

 16 Figure 16: Method I: example DEM tessellation with l h = 0.

2 Figure 17 :

 217 Figure 17: Method II: Example DEM tessellation with l h = 0.

Figure 18 :

 18 Figure 18: 3D views obtained with the original (l h = l m = l l = 0).

Figure 19 :

 19 Figure 19: 3D triangulation of the replication example from Figure 15(b) (l h = 0, l m = 2, l l = 3).

Figure 20 :

 20 Figure 20: 3D views after tessellation without smoothing (l h = 0, l m = 2, l l = 3).

  (a) 3D triangulation (b) overlaying its texture

Figure 21 :

 21 Figure 21: 3D views with the tessellation after smoothing with Method I: d m = 2, d h = 1 (l h = 0, l m = 2, l l = 3).

Figure 22 :

 22 Figure 22: 3D views with the tessellation after smoothing with Method II: d m = 2, d h = 1 (l h = 0, l m = 2, l l = 3).

(a) The original (2048 × 2048 pixels) (b) zoomed square of (a)

Acknowledgments

This work is in part supported by the Higher Education Commission (HEC), Pakistan, as well as by VOODDO (2008VOODDO ( -2011)), a project of ANR and the region of Languedoc Roussillon, France.

Preliminary call for papers

The 2011 European Signal Processing Conference (EUSIPCO 2011) is the nineteenth in a series of conferences promoted by the European Association for Signal Processing (EURASIP, www.eurasip.org). This year edition will take place in Barcelona, capital city of Catalonia (Spain), and will be jointly organized by the Centre Tecnològic de Telecomunicacions de Catalunya (CTTC) and the Universitat Politècnica de Catalunya (UPC). applications as listed below. Acceptance of submissions will be based on quality, relevance and originality. Accepted papers will be published in the EUSIPCO proceedings and presented during the conference. Paper submissions, proposals for tutorials and proposals for special sessions are invited in, but not limited to, the following areas of interest. • Multimedia signal processing and coding.

Areas of Interest

• Image and multidimensional signal processing.

• Signal detection and estimation.

• Sensor array and multi channel signal processing.

• Sensor fusion in networked systems.

• Signal processing for communications.

• Medical imaging and image analysis.

• Non stationary, non linear and non Gaussian signal processing.

Submissions

Montserrat Nájar (UPC)

Tutorials

Daniel P. Palomar (Hong Kong UST) Beatrice Pesquet Popescu (ENST)

Publicity

Stephan Pfletschinger (CTTC) Mònica Navarro (CTTC)

Publications

Antonio Pascual (UPC) Carles Fernández (CTTC)

I d i l Li i & E hibi

Submissions

Procedures to submit a paper and proposals for special sessions and tutorials will be detailed at www.eusipco2011.org. Submitted papers must be camera ready, no more than 5 pages long, and conforming to the standard specified on the EUSIPCO 2011 web site. First authors who are registered students can participate in the best student paper competition.

Important Deadlines:


D Tsirkov

V Komarov

T Azaryan

D Chiladze

S Dymov

A Dzyuba

M Hartmann

A Kacharava

A Khoukaz

A Kulikov

V Kurbatov

G Macharashvili

S Merzliakov

M Mielke

S Mikirtychiants

M Nekipelov

F Rathmann

V Serdyuk

H Ströher

Yu Uzikov

Yu Valdau

C Wilkin

Energy dependence of hard bremsstrahlung production in proton-proton collisions in the ∆(1232) region

Keywords: numbers: 25.40.Ep, 13.60.-r, 25.20.Dc

Hard bremsstrahlung production in proton-proton collisions has been studied with the ANKE spectrometer at COSY-Jülich in the energy range of 353-800 MeV by detecting the final proton pair {pp} s from the pp → {pp} s γ reaction with very low excitation energy. Differential cross sections were measured at small diproton c.m. angles from 0 • to 20 • and the average over this angular interval reveals a broad peak at a beam energy around 650 MeV with a FWHM ≈ 220 MeV, suggesting the influence of ∆(1232)N intermediate states. Comparison with deuteron photodisintegration shows that the cross section for diproton production is up to two orders of magnitude smaller, due largely to differences in the selection rules.

Introduction

The photodisintegration of the deuteron

γd → pn, (1)
is the simplest reaction of its type and has therefore attracted much attention, both theoretical and experimental, to learn more about how to treat the interaction of the electromagnetic field with nuclei. Significant advances in the understanding of the dynamics at several hundred MeV have been achieved within the framework of models that are either perturbative, taking into account a set of relevant Feynman diagrams [1], or employ a coupled NN, N ∆, ∆∆ and NNγ channel formalism [2,3]. Nevertheless, extensive work on the reaction continues and a full interpretation of the data has still to be achieved [START_REF] Glister | Polarization Observables in Deuteron Photodisintegration below 360 MeV[END_REF].

A kinematically similar reaction involves the photodisintegration of a 1 S 0 proton pair (diproton), here denoted by {pp} s , γ{pp} s → pp .

(

Reactions ( 1) and ( 2) share the common feature that, in the ∆(1232) range and above, there is a large energy transfer from the photon to the target nucleon pair, transforming it into one of high invariant mass and driving the final state deep into the nucleon resonance region [START_REF] Gilman | [END_REF]. Thus, in contrast to other electromagnetic processes, such as electron-deuteron elastic scattering, the nucleon resonance excitation channels are explicitly open. As a consequence, the influence of isobar or mesonic exchange currents is likely to be very strong for both reactions. Since the diproton is the spin-isospin partner of the deuteron, reactions (1) and (2) involve different transitions. Whereas one of the main driving terms for γd → pn is an M1 excitation of an S-wave ∆(1232)N pair that de-excites into pn, an analogous transition is forbidden by angular momentum and parity conservation for the γ{pp} s → pp reaction [1]. There are likely to be cancellations among the large amplitudes as one approaches the pure S-wave diproton limit and this means that extra insight may be gained into the field through a combined study of the photodisintegration of the deuteron and diproton.

However, the free 1 S 0 proton pair is not bound and the diproton investigation has generally been approached through the photodisintegration of a pp pair embedded in a light nucleus, in particular 3 He. By studying events where there were two fast protons emerging from the γ 3 He → ppn reaction with a slow (reconstructed) neutron, this part of phase space was primarily interpreted in terms of an interaction on a diproton, with the neutron merely appearing as a spectator particle, whose influence is only felt through the kinematics [6][7][8][9][10][11][12][13][14]. These data show little indication for the excitation of an intermediate ∆(1232)N state, certainly much less than for those for fast pn pairs [11,12].

As pointed out in several of these publications, the corrections to the simple spectator model picture can be quite large, owing mainly to the fact that the photoabsorption on pn pairs is much stronger than on pp. This makes it harder to separate cleanly these two terms, especially at the lower energies where the Dalitz plot is not sufficiently wide. Final state interactions, where the neutron from a fast pn pair undergoes a charge-exchange reaction on a spectator proton, and other three-body mechanisms have also to be considered.

In view of these potential drawbacks, it is well worth seeing if the information derived from these experiments could be complemented by a direct measurement of the cross section for pp → {pp} s γ .

(

) 3
In order to ensure that the final diproton is almost exclusively in the 1 S 0 state, it is crucial that the pp excitation energy E pp be very small, and this is an important constraint on any measurement. For small E pp , the proton beam energy T p is essentially twice the photon energy E γ for the inverse reaction.

There have been many measurements of bremsstrahlung in proton-proton scattering in the few hundred MeV range [15][16][17][18][19][20]. These were generally accomplished through the use of pairs of counters placed on either side of the primary beam direction. As a consequence they were not sensitive to the behaviour at low E pp , i.e. to the hardest part of the bremsstrahlung spectrum. The one early exception benefitted from the wide acceptance offered by the COSY-TOF spectrometer, where data were taken at 293 MeV [21]. However, the statistics in the low E pp region were severely limited and a cross section for this region was not given.

The first data on hard bremsstrahlung production in reaction (3) at intermediate energies were reported recently from an experiment carried out at the ANKE facility at the Cooler Synchrotron COSY-Jülich [22]. The pp → {pp} s γ differential cross section was measured at energies T p = 353, 500 and 550 MeV. Events were here selected with a final excitation energy E pp < 3 MeV, where it is expected that the contribution from P -waves in the diproton should be minimised. However, due to the limited angular coverage of the forward detector employed in these studies, only small diproton c.m. angles θ pp were covered; 0 • < θ pp < 20 • . Since the cross section is symmetric about 90 • , this effectively means a similar cut on the photon angle θ γ ; 0 • < θ γ < 20 • . In contrast, the full pp → {pp} s γ angular domain was measured in a high statistics experiment at CELSIUS with the same E pp cut [23], though only at the single (lower) energy of 310 MeV.

Taken together, the ANKE and CELSIUS data show a rise in the near-forward cross section from 310 to 550 MeV, but data at higher energy are needed to see if there is a maximum in the region of the ∆(1232)N threshold. To clarify the experimental situation, we have supplemented the earlier COSY-ANKE data [22] by measurements of the small angle pp → {pp} s γ differential cross section at T p = 625, 700 and 800 MeV. Combining the data sets allows us to study the energy dependence of the near-forward cross section throughout the ∆(1232) resonance region. Unlike some of the 3 He photodisintegration results [11,12], our data reveal a significant maximum in a region where ∆(1232)N intermediate states might be expected to play a role and this should provide a useful guide for further theoretical work.

In view of earlier publications, the description of the experiment and its analysis in section 2 could be made more brief. The results and their significance are discussed in section 3, with our conclusions and suggestions for further work being presented in section 4.

Measurement and Analysis

The experiment was carried out using the magnetic spectrometer ANKE [24], which is installed at an internal target station of COSY. A hydrogen cluster-jet target was positioned in the proton beam and the secondary particles were detected with wire chambers and a scintillation hodoscope. The three-momenta and trajectories of the particles were reconstructed on the basis of the known field map of the analysing magnet, assuming that these particles originated from a point-like source situated at the centre of the target-beam interaction volume.

The first step in the identification of the ppγ final state was the selection of two coincident protons from among all the detected pairs of positively charged particles. The scintillation hodoscope allowed the measurement of the difference between the times of flight from the target to the detector for the two recorded particles. The comparison of this value with that calculated from the measured particle momenta and trajectories led to a very good identification of proton pairs, as illustrated in figure 1. The background from misidentified pairs was at the few percent level for all beam energies, as detailed in Table 1. Having identified two protons and determined their momenta, the complete kinematics of the pp → ppX process could be reconstructed. Events with zero missing mass, within the experimental resolution, were accepted as candidates for the pp → ppγ reaction. Those where the pair's kinetic energy E pp in their rest frame is small, specifically E pp < 3 MeV, were classified as belonging to reaction (3).

The experimental techniques, conditions of the measurements, and the datahandling procedure have been described in an earlier publication [22] and initial results given there at 353, 500, and 550 MeV. The new results at energies of 625, 700 and 800 MeV were obtained in three separate beam periods. As a consequence, there were some differences in the adjustment of the magnetic system setup and the number of sensitive planes in the tracking detector. The 625 and 800 MeV data were collected as by-products of the study of other reactions so that the measurement conditions were not optimised for the study of the pp → ppγ reaction. In particular, the larger missing-mass widths found at these two energies were mainly the result of a lower field in the analysing magnet and thus a smaller deflection angle. The data at 700 MeV were obtained with a polarised proton beam and a weighted average was then taken over the two polarisation states.

In order to measure the excitation energy of the proton pair and study the angular dependence of the cross section of the pp → ppγ reaction, sufficient resolution in the Identification of the proton pairs from the pp → ppX reaction on the basis of time-of-flight information at 700 MeV. ∆t meas is the directly measured difference of the times of flight; ∆t calc is the time-of-flight difference calculated using the measured particle momenta and trajectories. The central peak corresponds to proton pairs, whereas those around ±15 ns are associated with pπ + pairs. The full line is a linear approximation to the background arising from accidental coincidences. Events inside the indicated ±5 ns interval were retained for the analysis.

Table 1. Characteristics of the measurements at different energies: integral luminosity L with systematic (first) and normalisation (second) errors; background/signal ratio N bg /N pp for proton pair identification; width FWHM(M 2 π ) of the π 0 peak in the missing-mass-squared distribution from the pp → {pp} s X reaction. The reasons for the broader widths at 625 and 800 MeV are discussed in the text. corresponding variables is needed. The uncertainty σ(θ pp ) in the polar angle θ pp of the diproton in the overall c.m. system ranged from 0.5 • to 2.3 • , depending on the beam energy and the value of θ pp . The uncertainty in the excitation energy E pp , which generally increased with E pp , was between 0.08 to 0.6 MeV for E pp < 3 MeV, depending also on the beam energy. The E pp resolution was thus sufficient for the measurement of the excitation energy spectra up to 3 MeV. The E pp spectra for pp → {pp} s γ, as well as for the pp → {pp} s π 0 reaction, are satisfactorily reproduced by Monte Carlo simulations, where events were generated according to phase space modified by the S-wave pp final state interaction [22,25]. Extra evidence for the S-wave nature of the proton pairs was provided by the isotropy of the acceptance-corrected angular distributions in the diproton rest frame. The distributions in missing mass squared M 2 x for events where θ pp < 20 • are shown in figure 2 for the six beam energies. At 353 MeV there is a clear γ peak that is well separated from the pion peak associated with the pp → {pp} s π 0 reaction. The low background arising from accidental proton coincidences changes weakly with M 2

T p L N bg /N pp FWHM(M 2 π )×10 -3 (MeV) (10 31 cm -2 ) ( ( G e V / c
x and could be taken as linear in the interval -0.02 < M 2

x < +0.06 (GeV/c 2 ) 2 .

›0.02 0 0.02 The shaded area corresponds to the γ peak, the dashed line to the π 0 peak, the dotted to the linear accidental background, and the solid to the sum of these three contributions. The data at the three lower energies were reported in [22], whereas the others are from this work.

The widths of the peaks in figure 2 are governed by the accuracy of the measurement of the proton three-momenta in the relevant experimental runs and this generally gets worse with increasing beam energy. The corresponding growth in the widths reported in Table 1 results at higher energies in the merging of the γ signal with that of the pion. The principal difficulty consists therefore in selecting the small number of γ events from the total distribution.

The M 2 x distributions were fitted by a sum of peaks, corresponding to γ and π 0 production, plus a linear background. In order to estimate the shapes of the peaks as reliably as possible, a detailed Monte-Carlo simulation was undertaken at each energy, taking into account all the known features of the setup. These include, in particular, the smearing caused by the radial distribution of the proton beam at the target region, the multiple scattering in the exit window of the vacuum chamber and detector materials, and the actual clustering of the wires fired in the proportional chambers. The procedure of track reconstruction used in the simulation was the same as that in the data handling. The free parameters of interest used to fit the missing-mass spectra were the number of events in the γ peak, the number of events in the pion peak, and the two constants of the linear background. However, in order to compensate for the lack of knowledge of the beam spatial distribution, additional parameters were inserted into the fits. These were the shift of the pion peak position and correction factors for the γ and pion peak widths.

In order to investigate the systematics, several fitting methods were employed. Either both peaks were fitted together or the pion peak with the background was first fitted separately, excluding the M 2

x range where the γ signal was expected. The correction for the setup acceptance was also accounted for in two different ways. This was either introduced for each event and the weighted M 2

x spectrum fitted, or the uncorrected distribution was fitted and the number of γ events corrected for the average acceptance factor afterwards. Finally, two fitting approaches were tried, namely χ 2 minimisation and logarithmic likelihood, though the latter was used only for fitting uncorrected spectra. Although these techniques were developed primarily in order to identify the γ signal in the more difficult high energy cases, they were also used on the published data [22], but any changes there are of little significance.

At each energy the results obtained using the different methods were completely compatible and their average was taken, with the variation in the deduced cross section being considered as a systematic uncertainty. As a further check on the reliability of the identification of the γ signal in the most complicated cases of higher T p , the M 2 x distributions were also fitted assuming that there were no γ events. In this case the π 0 peak was first fitted, excluding the γ range around M 2 x = 0, and the χ 2 of the difference between the resulting function and the histogram evaluated over the previously excluded γ range. The results of fitting the data at 625, 700 and 800 MeV with and without the γ contribution are shown in figure 3, where it is seen that, unless the γ is included, there is a significant deterioration in the values of χ 2 evaluated for the γ peak region even in the most severe case of 800 MeV.

The luminosity L was determined from the number of the elastically scattered protons detected in parallel. In order to correct for the acceptance, a simulation was carried out where the setup geometry, the efficiency of the proton detection by the multiwire proportional chambers, and the track reconstruction algorithm were taken into account. These measurements were then compared to differential cross sections predicted by the SAID analysis program [26]. Although this does not yield errors, study of experimental data in the literature suggests that the uncertainties range from about 3% at 353 MeV up to 5% at 800 MeV. Small variations of the luminosity (those exceeding the statistical fluctuations), derived at different proton angles, reflect uncertainty in the acceptance correction and were taken as systematic errors in the luminosity.

The total uncertainty was estimated as the quadratic sum of the statistical, systematic, and normalisation errors. As a general rule, statistical errors had a negligible x distributions at 625, 700 and 800 MeV, assuming the presence (upper panels) or absence (lower panels) of a γ contribution, as explained in the text. The notation is as in figure 2. The values χ 2 per degree of freedom in the γ peak region increase significantly for fits that do not include the possibility of a γ signal. At 625 MeV the χ 2 /ndf goes from 20.8/9 to 69.8/10, at 700 MeV from 9.7/9 to 314/10, and at 800 MeV from 9.7/9 to 36/10. effect so that they are not shown in Table 1. Also negligible is the contribution arising from the systematic error in the proton momentum.

The angular dependences at 353, 500 and 550 MeV were reported earlier [22]. The poorer conditions did not allow similar studies to be made at 625 and 800 MeV. However, to estimate the angular dependence of the differential cross section at 700 MeV, we have divided the events into three θ pp intervals, 0 • -7 • , 7 • -12 • , 12 • -20 • , and made separate fits for each of these ranges.

Results and Discussion

In order to study the energy dependence of the pp → {pp} s γ reaction from 353 to 800 MeV, we have evaluated the differential cross section dσ/dΩ (0-20) averaged over the interval 0 • -20 • in diproton c.m. angle. Our measurements show that the angular dependence is rather smooth in this region [22] so that the average cross section should reflect reasonably well the energy dependence at fixed angle.

The numbers of selected γ events, N raw , and the corresponding acceptance-corrected figure, N corr , are given in Table 2 for the six energies. Also to be found there are the average cross sections dσ/dΩ (0-20) , together with the corresponding partial and total errors. The uncertainties are particularly large at 625 and 800 MeV, due to the low luminosity and non-optimal conditions that hampered the identification of the γ signal.

Table 2. Numbers of pp → {pp} s γ events and the corresponding differential cross sections at different beam energies. Data at 353, 500, and 550 MeV were already presented in [22]. Here N raw is the number of registered {pp} s γ events, N corr the number of events corrected for acceptance, and dσ/dΩ (0-20) the differential cross section averaged over 0 • -20 • . The statistical error is denoted by σ stat , that coming from the systematics in the γ event acquisition by σ γ , and that arising from the luminosity uncertainty by σ lum . Adding these contributions quadratically gives a total error of σ tot . As previously remarked, it was also possible to extract information on the angular dependence of the pp → {pp} s γ reaction at 700 MeV by dividing the ANKE range into three intervals and the results are presented in Table 3. These data are compared in figure 4 with those obtained at lower energies [22]. Although the new data show some tendency for a forward dip, this is less strong than that measured at 500 and 550 MeV. In all cases, over the limited angular interval reported here, the cross section is consistent with a linear variation in cos 2 θ pp . Table 3. Angular dependence of the pp → {pp} s γ reaction at 700 MeV.

T p (
θ pp 0 • -7 • 7 • -12 • 12 • -20 • dσ/dΩ (nb/sr) 34±5 37±4 42±6
The measured values of the average cross section dσ/dΩ (0-20) are shown in figure 5a as a function of the proton beam energy. Similar data are also available with the identical 3 MeV cut in E pp at 310 MeV from CELSIUS [23]. Although this experiment had a much wider angular coverage, the value shown in the figure was obtained by taking the 0 • -20 • average. This point joins very smoothly onto the ANKE data.

The energy dependence of the differential cross section of the pp → {pp} s γ reaction shown in figure 5a reveals a broad peak, with a maximum near 650 MeV and FWHM ∼ 220 MeV, and it is tempting to suggest that this might be associated in some way with the excitation of the ∆(1232) isobar. It is well known that this resonance plays a crucial role for the analogous pn → dγ reaction, whose energy dependence at a deuteron production angle of θ d = 20 • with respect to the proton direction is also shown in the same figure, scaled down by a factor of twenty. These values were obtained from the MAMI deuteron photodisintegration γd → np results [27] by using detailed balance.

It is important to note that the maximum in the pp → {pp} s γ cross section appears at about 80 MeV higher than that for pn → dγ. Direct excitation of the ∆-isobar dominates the pn → dγ reaction in this energy range through an M1 transition to an intermediate ∆N state in an S-wave. Such an M1 transition is forbidden in the diproton case [1] but the strength could come from a P -wave ∆N configuration, which requires extra energy in order to overcome the centrifugal barrier. Similar arguments have been used to explain why the ∆ peak is shifted to higher energy also for the pp → {pp} s π 0 reaction [28]. In addition to the displacement of the pp → {pp} s γ peak position to higher energies, it is narrowed significantly compared to pn → dγ.

The ratio of the two differential cross sections

R = dσ/dΩ(pn → dγ) dσ/dΩ(pp → {pp} s γ) (4)
is presented in figure 5b. Where necessary, the pn → dγ data [27] were splineinterpolated in energy. The smallest c.m. angle in these measurements was 20 • and these are compared in the figure to the diproton results both at this angle and averaged between 0 • and 20

• .
What is striking about the ratio shown in figure 5b is the smallness of the diproton cross section compared to that for the deuteron. The factor is strongly energy dependent, dropping smoothly from over 100 at low energies to about 20 at 800 MeV. It should, however, be borne in mind that some of this suppression is to be associated with the difference in the phase space volume for producing the deuteron bound state and diproton continuum state. The theoretical treatment of the pp → {pp} s γ reaction has been far less developed than that for pn → dγ. The first thing to note is that, in the diproton case, the E1 transition is suppressed by the vanishing of an electric dipole operator for the proton pair. The absence of charged pion exchange currents also reduces some of the meson exchange effects. Finally, because the S-wave J P = 1 + ∆N intermediate state cannot couple to an initial pp system [1], the M1 transition that dominates the pn → dγ reaction in the few hundred MeV region is also not present. One might therefore expect the first effects of the ∆ degrees of freedom to show up in an S-wave J P = 2 + ∆N contribution, which would lead to an E2 transition, or through a ∆N configuration in a P or higher wave. The coupled-channel calculations of [29] suggest that at intermediate energies the E2 multipole should dominate and this would give rise to a sin 2 θ γ cos 2 θ γ dependence. In addition, this model showed no sign for any apparent ∆N signal in the energy variation, in contrast to our data. Furthermore, the CELSIUS data at 310 MeV show a linear variation with cos 2 θ γ so that any E2 contribution at this energy should be very small.

In a recent publication [30] the high precision KVI pp bremsstrahlung data at 190 MeV [19] have been successfully described for the first time through the introduction of a phenomenological contact interaction current that explicitly satisfies the generalised Ward-Takahishi identity. However, these data did not sample the hard bremsstrahlung limit. A generalisation of the theoretical model includes some ∆ contributions, but it still neglects entirely the interaction current due to the ∆ that gives rise to the five-point contact current [31]. It is not clear if this is the reason why it did not reproduce the shape of the angular distribution for small E pp from CELSIUS at 310 MeV [23].

A peak in the energy dependence of the pp → {pp} s γ differential cross section due to the ∆N configuration was obtained in calculations [START_REF] Yu | Proc. XIX Int. Baldin Seminar on High Energy Physics Problems[END_REF] within the framework of a simplified one-pion-exchange model. However, gauge invariance was not imposed and other terms neglected.

Summary and conclusions

The differential cross section for hard bremsstrahlung production in proton-proton collisions has been measured in the near-forward direction, 0 • -20 • , for beam energies between 353 and 800 MeV. The energy dependence of this cross section reveals a broad peak around 650 MeV, which is roughly where one might expect a contribution from a ∆(1232)N intermediate state. The peak is shifted to higher energies compared to that in the analogous pn → dγ reaction but this is not surprising because the S-wave J p = 1 + ∆N state that drives the M1 transition in the deuteron case does not couple to an initial pp system.

The suppression of the isobar M1 term leads to a large but energy-dependent factor between the cross sections for pn → dγ and pp → {pp} s γ. It also makes the theoretical description harder to realise because there is then no longer an obviously dominant term to consider. Further theoretical work is clearly needed.

Even with a 3 MeV cut in E pp there might be some small contamination of Pwaves in the pp system but there are similar concerns for the 3 He(γ, pp)n data because variational Monte Carlo calculations suggest that the pp pair in 3 He is not in a pure 1 S 0 state and might also contain a few per cent of higher partial waves [START_REF] Wiringa | [END_REF].

The combined study of the pn → dγ and pp → {pp} s γ reactions has assumed a greater importance in recent years because of the interest in the investigation of short-range pn and pp correlations in nuclei [14,34,35]. These are studied through photoabsorption that leads to the emission of nucleon pairs with large back-to-back momenta in the pair's rest frame. On the experimental side, the separation of the various multipoles in the pp → {pp} s γ reaction would clearly require the measurement of the cross section (and analysing power) over a wider angular region and this will be possible at ANKE through the use of a positive side detector in combination with the forward detector [36] that provided the data reported here. Of great use in the separation would be data on the proton-proton spin correlation parameters and it might be possible to study these at COSY-ANKE through the use of a polarised gas cell target [START_REF] Kacharava | Spin Physics from COSY to FAIR[END_REF].

New experimental data [14] has allowed the comparison of the cross sections for d(γ, p)n and 3 He(γ, pp)n at θ cm = 90 • . These data come from the scaling energy region, E γ > 2.4 GeV, where the influence of the additional nucleon in the 3 He target is of less importance. The ratio of the cross sections found there was ≈ 40, but to compare with the R ratio of equation 4, one needs to transform data with a bound diproton to ones with a scattering pp system, which requires a model calculation.

An extension of the experimental study of the pp → {pp} s γ reaction from the ∆(1232) excitation region to GeV energies could provide an alternative way to investigate the transition to a situation where hadronic internal degrees of freedom determine the interaction. The onset of the QCD scaling regime might take place in pp → {pp} s γ at a beam energy T p ≈ 2 GeV, corresponding to the E γ ≈ 1 GeV which, it is suggested, is the start of the domain for deuteron photodisintegration [START_REF] Rossi | [END_REF].

Although the maximum proton beam energy at COSY is nearly 3 GeV, the identification of the pp → {pp} s γ reaction through the pp missing mass becomes progressively harder as the energy increases, as is well illustrated by the data in figure 2. Therefore any attempt to use this reaction for the investigation of the γNN dynamics and its relation to the quark degrees of freedom would certainly necessitate the detection of the γ in coincidence. This is a challenge for the future.

Figure 1 .

 1 Figure1. Identification of the proton pairs from the pp → ppX reaction on the basis of time-of-flight information at 700 MeV. ∆t meas is the directly measured difference of the times of flight; ∆t calc is the time-of-flight difference calculated using the measured particle momenta and trajectories. The central peak corresponds to proton pairs, whereas those around ±15 ns are associated with pπ + pairs. The full line is a linear approximation to the background arising from accidental coincidences. Events inside the indicated ±5 ns interval were retained for the analysis.

Figure 2 .

 2 Figure 2. Distributions of the missing mass squared in the pp → {pp} s X reaction for proton pairs with θ pp < 20 • compared to fits with different contributions. The expected π 0 position is indicated by the arrow. The shaded area corresponds to the γ peak, the dashed line to the π 0 peak, the dotted to the linear accidental background, and the solid to the sum of these three contributions. The data at the three lower energies were reported in[22], whereas the others are from this work.

Figure 3 .

 3 Figure 3. Descriptions of the M 2x distributions at 625, 700 and 800 MeV, assuming the presence (upper panels) or absence (lower panels) of a γ contribution, as explained in the text. The notation is as in figure2. The values χ 2 per degree of freedom in the γ peak region increase significantly for fits that do not include the possibility of a γ signal. At 625 MeV the χ 2 /ndf goes from 20.8/9 to 69.8/10, at 700 MeV from 9.7/9 to 314/10, and at 800 MeV from 9.7/9 to 36/10.

Figure 4 .

 4 Figure 4. Angular dependence of the differential cross section for the pp → {pp} s γ reaction at four beam energies. The results of the present measurement at 700 MeV are shown by stars. Also shown are the data from [22] at 353 MeV (circles), 500 MeV (squares), and 550 MeV (triangles). The errors coming from the statistics and background subtraction are shown but not the overall ones arising from the luminosity uncertainty. The lines represent linear fits to the four data sets.

Figure 5 .

 5 Figure5. (a) Energy dependence of the pp → {pp} s γ and pn → dγ differential cross sections. Full circles represent COSY-ANKE data for the average differential cross section dσ/dΩ (0-20) for the pp → {pp} s γ reaction whereas the triangle is taken from the CELSIUS results[23] for the same angular and E pp conditions. The open circles show the values of the differential cross section for the pn → dγ reaction at a deuteron c.m. angle of θ d = 20 • with respect to the proton direction, evaluated from the deuteron photodisintegration data of[27] and scaled down by a factor of twenty. (b) Ratio of the differential cross sections for pn → dγ to pp → {pp} s γ, as defined by equation (4). The stars were obtained for the same production angle θ d = θ pp = 20 • while for the circles the angular average of the pp → {pp} s γ cross section over 0 • -20 • was used.

Acknowledgments

The authors wish to thank other members of the ANKE collaboration for their help and assistance in the running of the experiment. We are grateful also to the COSY crew for providing good working conditions. Correspondence with Kanzo Nakayama has been very helpful. This work has been partially supported by the BMBF (grant ANKE COSY-JINR), RFBR (09-02-91332), DFG (436 RUS 113/965/0-1), the JCHP FFE, and the HGF-VIQCD.