Using MC Algorithm to Implement 3 D Image Reconstruction for Yunnan Weather Radar Data

3D image reconstruction for weather radar data can not only help the weatherman to improve the forecast efficiency and accuracy, but also help people to understand the weather conditions easily and quickly. Marching Cubes (MC) algorithm in the surface rendering has more excellent applicability in 3D reconstruction for the slice images; it may shorten the time to find and calculate the isosurface from raw volume data, reflect the shape structure more accurately. In this paper, we discuss a method to reconstruct the 3D weather cloud image by using the proposed Cube Weighting Interpolation (CWI) and MC algorithm. Firstly, we detail the steps of CWI, apply it to project the raw radar data into the cubes and obtain the equally spaced cloud slice images, then employ MC algorithm to draw the isosurface. Some experiments show that our method has a good effect and simple operation, which may provide an intuitive and effective reference for realizing the 3D surface reconstruction and meteorological image stereo visualization.


Introduction
In recent years, most of the meteorological radar products are still expression by using two-dimensional data structure to represent three-dimensional cloud state.This process is both abstract and difficult to understand.The three-dimensional image reconstruction for the radar data can not only transform the abstract meteorological data into the visual image information which make people understand and accept easily, but also explore the further potential information from the meteorological targets through spatial three-dimensional modeling and stereoscopic displaying, provide more accurate data basis for the short-term weather forecasts.
Three-dimensional image reconstruction has more extensive and mature application in medical image processing, mechanical manufacturing, computer animation, film and television design and many other fields [1] [2] [3].In the meteorological field, owing to the weather objects in the entire detection space are discrete, irregular and uneven distributed, most of the visualization products are confined into the two-dimensional structure, as shown in Figure 1, the images reconstructed from the radar data with two elevation angles in Kunming Radar Station at 18:23 on July 18, 2013.If we could provide the three-dimensional cloud structure of this moment by using the three-dimensional reconstruction technology, we would better understand more meteorological object feature information and know the weather changes.
There are two kinds of three-dimensional reconstruction technologies.One is known as the surface drawing method, which approached through the geometric unit splicing fitting surface to describe the object three-dimensional structure, the classical algorithms including Marching Cubes (MC), Marching Tetrahedral and Dividing Cubes, etc.; another is called the volume rendering method based on the volume data, it projected directly the voxel into the display plane, the classic algorithms have: light projection, projection imaging and Frequency domain transform method.The commonality of these algorithms is that using many triangular patches to fit the equivalent surfaces, with the advantage that hardware acceleration can be achieved when rendering some 3D reconstructed images [4]- [9].But more intermediate data may be generated from the original data to the final display; the amount of calculation is very large.
Due to the radar detecting elevation angle is small, the data form is also close to the plane data, then we can obtain the equally spaced cloud slice images by Firstly, we need to establish the 3D data field because the weather radar data are not the equidistant slice images data, then use the MC algorithm to reconstruct the 3D cloud image of the weather radar, and realize the stereo visualization of the two-dimensional meteorological data.

Cube Weighting Interpolation Method
Based on the working principle of the weather radar, we know that each layer of the scanning data obtained by the weather radar is a cone scanning surface, all the objects in one cone scanning surface have the different high attitude, we must search each layer of scanning data to obtain the specific attitude firstly.Secondly, we compute contour graphic data in specific attitude by doing some data preprocessing, such as the coordinate transformation, data interconnection and blindness etc.After above processes finished, the object three-dimensional data field is formed.Then, we use spatial interpolation to produce a three-dimensional weather object image.
In order to establish the 3D data field, we propose a projection method called Cube Weighting Interpolation (CWI) to obtain the equally spaced cloud slice images.
CWI includes the following three steps: 1) Standardize the txt data and map it into the 3D space; 2) Build standard 3D grid data; 3) Correct and optimize the grid data.
Step 1: Setting the original data to ( ) Let the mapping data be ( ) where, , , x y z represent the spatial position coordinates separately when the original data mapping to the three-dimensional space, α is the angle between the scanning cone surface and the horizontal plane, i d indicates the echo in- tensity at the coordinate position, each mapping data i d corresponds to the echo intensity data about a space sampling point during the observation process.
Step 2: It must to use spatial interpolation, since the radar scanning sampling data is distributed in 14 tapes, and distribution of the data samples is irregular and uneven.
For any sample point i d , its real spatial position falls into a unit cube M surrounded by an integer axis with a length and width of 1 km, and the distance between the sampling point i d and the cube M is close, as shown in Figure 2.
Assumed that vertex D is center of sphere whose radius is length of the unit cube M in Figure 2, and these curves in Figure 2, which indicate 1 8 spherical surface of sphere.The weight of each sampling point i d can be calculated ac- cording to the magnitude of the distance between the sampling points i d and vertex D in the eight unit cube adjacent to vertex D, and the weighted average of the sampling points i d can be used as estimated echo intensity value of the vertex D, and specific equation is as follows: , , Let the echo intensity be ( ) , , d x y z expresses the observation data that stored in the three-dimensional matrix ( ) , , x y z , and also represents echo intensity at the ( ) , , x y z in real space.as a placeholder, which is helpful to repair data after that.
( ) where: In the Equation (4), ( ) q i means that the weight is assigned to the different sampling points.γ is a parameter related to the weight distribution, which de- termines the relation between the weight assignment and l ∆ , here γ is 1.
When 1 γ = , there is an inverse relationship between the weight size and l ∆ .It is worth noting that ( ) q i is the weight, and its minimum value is 0, and is not negative.So when ( ) q i is negative, we unified with 0 instead it.
In the Equation ( 6), using the commands round up and round down in MATLAB, we can get quickly the coordinate distance between the data sample points i d and the adjacent eight standard cube vertices.The coordinates of the eight vertexes of the standard cube associated with the data sample points can be quickly obtained by combination of x, y, and z with rounding up or rounding down respectively, and the echo intensity of the sample point i d can be as- signed to the eight vertices according to the weight rule.
Through this step, we can establish a three-dimensional space grid with 1 km interval.We call this grid as a standard grid.The data d of grid intersection is the basis of the next three-dimensional display and model reconstruction.
Step 3: In step 2, we have faced two major problems: 1) Data sampling points are sparse where away from the radar, there will be the case of n = 0 in the Equation ( 3), there is no data to estimate the vacuum point; 2) Data samples are dense and the large amount of data is allocated where near the radar.Although it helps to offset the outliers, the influence of the key data is weakened and the estimation results are balanced and cannot be fully characterized.
In order to solve the above two problems, this paper uses a method to change the sphere radius of the estimated grid point data.When constructing a standard grid point, we always use 1km as the sphere radius of the point D. If the sphere radius was increased, we could reduce the risk of problem 1, and if the sphere radius was reduced, we could alleviate conflict of problem 2.
The concrete operation of this method is to do the scaling of the real space in ( ) , , i i i i d x y z , and when the coordinates ( ) , , x y z become half of the original, the grid node spacing will turn into 2 km, and sphere radius of the data point D also turn into 2 km.Then the number of data in the sphere will increase, which will alleviate problem 1.Although the mesh density is reduced by half, the net interpolation can be used to compensate for the mesh density, and we can get a grid with a precision of 2 km.The grid data is recorded as ( ) x y z × , the grid node distance 0.5 km, the grid density becomes Grid density becomes 2 times the original.
Then we restore the original grid density by extracting even bits in grid nodes as the target grid, this grid alleviates the contradiction of problem 2, and its data is recorded as , , x y z d compared to the three, their estimated accuracy sort is from high to low, their vacuum points number is from more to less.Therefore, the final corrected mesh data is filled by ( )  , ,  d x y z .Experiments show that, after filling the grid da- ta, 18 km above sea leave where still does not appear data vacuum point.Clouds above 18 km are scarce according to actual observation, which are almost negligible, thus the results ( ) , , d x y z of the composite structure can meet the re- quirements of observation and analysis.
The details of CWI flow diagram can be shown as Figure 3.
After above three steps, we can get the three-dimensional data field, as shown in Figure 4.

Using MC Algorithm to Reconstruct 3D Image of Cloud Data
Marching Cubes (MC) algorithm is the most representative method among the isosurfacing algorithm with the three-dimensional spatial rule data field.This method is simple and easy to implement, additionally it has been used widely   The essence of the MC algorithm is to view a series of two-dimensional slice data as a three-dimensional data field, from which extract a substance with a certain threshold, and then using a topological form to connect as a triangular.
The core of the MC algorithm is to individually deal with each voxel in the three-dimensional data field, then use the compare data between each peak in each voxel and the isometric surface to determine the intersection of the voxel and the isometric surface.Computing the isosurface normal vector to draw the equivalent surface, among which the determination of the intersection form between the voxel and the isometric surface and the drawing of the isosurface are two key operations.

Voxel Model
First, we evenly sample the data in the 3D data field.Let , , x y z , in the direction of the sampling interval corresponds to , , x y z ∆ ∆ ∆ .We use this sample pitch to sample the three-dimensional data field.The cube region, formed by eight adjacent sampling points, is just a voxel.Simultaneously, these eight sample points referred as the corners of the voxels.Supposing these eight points in the spatial position index are ( ) + , as shown in Figure 5 [15].Assumed ( )

6
, , P x y z is any point within the voxel, it is necessary to convert its spatial coordinates into image coordinates ( ) , , i j k in order to calculate the value of voxel, where respectively the position index of the voxel in directions of X, Y, and Z in the body data.Using the three-linear interpolation can obtain the value of voxel for the point, and its value can be expressed as ) ( )( ) where: ( ) ( ) By the above equation, we can get the equation as follows: , , f x y z a a x a y a z a xy a yz a xz a xyz In Equation ( 9), 0 1 , , f i j k .The volume data is a three-dimensional array, that the sizes is L M N × × , as shown in Figure 6.

Algorithm Steps
The basic assumption of MC algorithm is that the bulk data is locally continuous.According to this assumption, if two adjacent sampling points are greater than and less than the value of isosurface, there be an equivalent point on the edges of them.If we get the equivalent points on each edge, we can regard these points as the vertices, and use a of triangles to approach isosurface.
The main steps of MC algorithm reconstruct the three-dimensional image are as follows: • The three-dimensional discrete rule data field is read hierarchically.
• Scanning two layers of data, to construct voxels one by one.The 8 corners of each voxel take the adjacent two layers.
• Comparing the values of each corner of the voxel with a given isosurface dimension.According to the comparison result, the index table of the voxel is constructed.
• According to the index table, we obtain the boundary voxels, which have an intersection with the isometric surface.
• Using the linear interpolation algorithm to calculate the intersection of voxel edge and isometric surface.
• Using firstly the central difference method to calculate the normal vector at the corners of voxels.Then calculating the normal vector at the vertex of the triangular patches by the linear interpolation method,.
• Finally, drawing isomorphic images according to the coordinates of each triangle points at the coordinates and the normal vector.
After the above steps, we can get a three-dimensional reconstruction map of weather radar data.Firstly, we can get triangle Grid of clouds data as shown in Figure 7.Then, we get isosurface map of cloud data.asshown in Figure 8.

Experiment and Results
We got the cloud map that meteorological department need, by adding a certain realistic rendering to map that have been reconstruction with MC algorithm.
∆ represents the distance between the sampling point and the standard cube vertex D.

Figure 6 .
Figure 6.Volume data based on Voxel.

(
which are uniquely determined by the values of the eight vertices in the voxel.According to Equation (9) we can obtain the property value of any point inside the voxel.Assumed that the three-dimensional data field constructed by the slice sequence of the meteorological image is a rectangular that the length is L VX × , the width is M VY × , and the height is N VZ × .The surface of the rectangular is parallel to the coordinate plane, and for the voxel ( ) value of volume data is( )

Figure 9 (
Figure 9(a) is cloud section map which the Kunming radar station on the ground height of 4 km was obtained at 18:23 on July 18, 2013.Figure 9(b) shows cloud section map on the ground height of 7 km for the same time.Figure 10

Figure 10 shows
the different intensity maps for the same time at 5 km.Weatherman can employ Figure 9 and Figure 10 to grasp and understand state about cloud, such