User:Shishami/Surface-splatting

The Surface Splatting algorithm can be used to render graphical surfaces from point clouds. It was developed by Matthias Zwicker, Hanspeter Pfister, Jeroen van Baar and Markus Gross in 2000 based on Elliptical Wieghted Average (EWA) filtering which was developed by Ned Greene and Paul Heckbert in 1986. Surface splatting algorithm supports order independent transparent surfaces, anti-aliasing, anisotropic filtering and hidden surface removal.

Algorithm
The pseudo-code algorithm of surface splatting is as follows:

1 SurfaceSplat(points P[k]) 2 {  3      for each point P[k] 4         project P[k] to screen space 5         calculate resampling kernel ρ[k] 6         splat ρ[k]; 7     for each pixel s in frame buffers 8         shade s;  9  }

The input of the algorithm is a set of points P[k] that represent the irregularly spaced points in the 3D space that we wish to apply surface spatting algorithm on (e.g. point cloud). Each point is represented by a 3 dimensional position and the normal vector that determines its front and back. Input points either have color attributes or color attributes are calculately later in the algorithm from 2D textures that we wish to apply on the final result.

Mapping
The second line in the code is to find the projection of these points from the input 3D space into screen space coordinates. This is done through a series of affine transformations to:
 * map objects to camera space (the virtual camera that the final view of the object is based on its position and angle)
 * map perspective projection to screen space (Graphical projection)
 * map viewport to viewport coordinates



Resampling Kernel
The third line of the pseudo-code is to calculate the resampling kernel (ρ) for each pixel. Resampling kernel can be computed as:


 * $$\rho(x) = \frac{1}{\mathbf{J}^{-1}}\mathcal{G}_{\mathbf{J}V_{k}^{r}\mathbf{J}^T + I}(x-m(u_k))$$

Symbol $$\mathbf{J}$$ denotes the Jacobian of the projection that is performed to map objects to viewport. As mentioned previously, this mapping has three steps; hence its Jacobian has three parts. For the first and third part (the first and third step of projection), we only allow uniform scaling and translation in our mappings, so their Jacobian can be represented by uniform scaling factors $$s_{mv}$$ and $$s_{vp}$$; leaving us with the Jacobian of the second step, perspective projection ($$\mathbf{J_{pr}}$$):
 * $$\mathbf{J} = s_{vp}\cdot\mathbf{J_{pr}}\cdots_{mv}$$

To compute $$\mathbf{J}$$ by $$\mathbf{J_{pr}}$$, we have to compute the the local 2-dimensional surface parameterization of the point ($$\tilde{x}_1$$ and $$\tilde{x}_2$$). We can find these vectors by projecting the screen basis vectors ($$x_1$$ and $$x_2$$ along the line between $$P_k$$ and its projection on screen space. We then select the normalized vectors $$u_0 = \frac{\tilde{x}_0}{||\tilde{x}_0||}$$, $$n_k$$ as the normal vector of point $$P_k$$ and $$u_1$$, so that they form a right handed orthogonal coordinate system. Finally we can calculate the inverse of $$\mathbf{J_{pr}}$$ as:


 * $$\mathbf{J_{pr}^{-1}} =

\begin{pmatrix} ||\tilde{x}_0|| & \tilde{x}_1 \cdot u_0 \\ 0 & \tilde{x}_1 \cdot u_0 \end{pmatrix}$$ The remaining part of the formula for calculating kernel function ($$\rho_k$$) is an ellipical Gaussian ($$\mathcal{G}_{\mathbf{J}V_{k}^{r}\mathbf{J}^T + I}^{(x-m(u_k))}$$) with variance matrix $$\mathbf{J}V_{k}^{r}\mathbf{J}^T + I$$ where:
 * $$V_{k}^{r} = \begin{pmatrix}

\frac{1}{h^2} & 0 \\ 0 & \frac{1}{h^2} \end{pmatrix}$$ The value of $$h$$ is specified in different ways, with the simplest one being the maximum distance between points in a small neighberhood.

Splatting
The value for each pixel is calculated as follows and stored in a buffer:
 * $$g(x) = \sum_{k \in N}w_k\frac{\rho_k(x)}{\sum_{j \in N} \rho_j(x)}$$

As the formula suggests, the value for each pixel is the weighted accumulated contribution of all pixels; although this is usually limited to a certain number to increase performance. The $$w_k$$ denotes the value of each pixel either by scalar pixel color values (RGB) or the values derived from texture. For more on textures, refer to.

Shading
Along with the frame buffer that contains color values for pixel, the algorithm uses another frame buffer to store z values of each pixel in camera space. The z value of each pixel is equal to z of tangent plane of the points that contribute to that particular pixel. This value helps the algorithm decide whether the contribution of a point should be added to the pixel or not. In case the buffer is empty (e.g. right after start of the algorithm), the contribution is added and z value is stored inside buffer. Otherwise: For the case of transparent objects (based on Alpha value in RGBA), the algorithm can be extended same as in common transparency rendering techniques: by using a number of different frame buffers, we can store the contribution of a point to the pixel of one of these buffers, decided based on the z values. One can think of these buffers as layers, the contribution of further points are stored in further layers. Finally, the values of one pixel in different buffers are blended together based on their accumulated Alpha. This formula is used to determine the new color ($$c_0$$) and ($${\alpha}_f$$) values of a pixel based on two pixels in the same position (one being behind and the other in front):
 * if the difference of new z value with previously stored z is close to zero (based on a threshold), we can understand that the new point belongs to the same surface as the points that contributed to that pixel earlier, hence the new contribution should also be added.
 * otherwise, if the value of new z value is less than old z, we understand that the new point belongs to a surface which is in front of the surface and points that contributed earlier; hence the new contribution replaces the older contribution.
 * in case the new z is larger than old z, its surface is behind the current surface and its contribution is discarded.


 * $$c_0 = c_f{\alpha}_f + c_b{\alpha}_b(1-{\alpha}_f)$$
 * $${\alpha}_0 = {\alpha}_f + {\alpha}_b(1-{\alpha}_f)$$

This is done repeatedly for the pixels of all buffers from behind to the front to achieve transparent surfaces.

Time and space complexity
The space complexity of the algorithm is 38 bytes per pixel. 16 bytes are used for RGBA values, 12 bytes for XYZ coordinates, 4 bytes for accumulated contributions of points, 2 bytes for material index (used in rendering) and 4 bytes for the z value (if ignoring transparency). The time complexity of the algorithm can be described as $$O(N^2)$$ where $$N$$ is the number of points. This is due to the fact that for every point, the Gaussian function takes all existing points as support (determined by the value of $$h$$). This can be avoided by only considering a neighbourhood of close points. There are time consuming calculations in this algorithm as well. For example at the final stage, the value of each pixel should be evaluated from frame buffers: $$O({screen size} \cdot {number of buffers})$$.

Logic
The algorithm works based on the fact the color of every pixel of the screen is the mixture and blend of the points that fall into and close to that pixel, at a certain camera view. The effect of the color of a point on a certain pixel is determined by the number and density of its neighbours (or all existing points). And also by applying the algorithm for all of the points and not discarding points that behind other points, the algorithm can take into consideration the transparency of surfaces using the alpha value.

Applications
There are many applications for surface rendering techniques based on point clouds. Laser and optical range scanners need these techniques since their output is a number of spaced points and their 3D coordinates, which are not rendered as surface if simply given to a rendering program. This is applicable to rendering 3D geology models, manufacturing and design, home 3D scanning and many more.

Variations
Different variations and modifications of this algorithm exist. Some of these variations have GPU based implementation that allows parallelized computation. There are also variations that are optimized for specific purposes (e.g. facial rendering ). Some variations are based on Phong shading techniques to provide more detailed lighting of surfaces.