Skip to content

compactly_supported_normal

compactly_supported_normal(x, n=2, sigma=1.0, center=None, second_derivative=-1)

Approximated compact Gaussian density functiona

This computes a compactly supported approximation of a Gaussian density function by convolving a box filter with itself n times in each dimension and multiplying, as described in "Poisson Surface Reconstruction" by Kazhdan et al. 2006.

Parameters:

Name Type Description Default
x (nx,dim) numpy double array

Coordinates where the function is to be evaluated

required
n int, optional (default 2)

Number of convolutions (between 1 and 4), a higher number will more closely resemble a Gaussian but have broader support

2
sigma double, optional (default 1.0)

Scale / standard deviation function (bigger leades to broader support)

1.0
center

Coordinates where the function is centered (i.e., distribution mean). If None, center at the origin.

None
second_derivative (int, optional(default - 1))

If greater or equal than zero, this is the index of the dimension along which to compute the second partial derivative instead of the pure density.

-1

Returns:

Name Type Description
v (nx,) numpy double array

Values of compacly supported approximate Gaussian function at x

Examples:

from gpytoolbox import compactly_supported_normal
import matplotlib.pyplot as plt
x = np.reshape(np.linspace(-4,4,1000),(-1,1))
plt.plot(x,compactly_supported_normal(x, n=4,center=np.array([0.5])))
plt.plot(x,compactly_supported_normal(x, n=3,center=np.array([0.5])))
plt.plot(x,compactly_supported_normal(x, n=2,center=np.array([0.5])))
plt.plot(x,compactly_supported_normal(x, n=1,center=np.array([0.5])))
plt.show()
Source code in src/gpytoolbox/compactly_supported_normal.py
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def compactly_supported_normal(x,n=2,sigma=1.,center=None,second_derivative=-1):
    """Approximated compact Gaussian density functiona

    This computes a compactly supported approximation of a Gaussian density  function by convolving a box filter with itself n times in each dimension and multiplying, as described in "Poisson Surface Reconstruction" by Kazhdan et al. 2006.

    Parameters
    ----------
    x : (nx,dim) numpy double array
        Coordinates where the function is to be evaluated
    n : int, optional (default 2)
        Number of convolutions (between 1 and 4), a higher number will more closely resemble a Gaussian but have broader support
    sigma : double, optional (default 1.0)
        Scale / standard deviation function (bigger leades to broader support)
    center: (dim,) numpy double array (default None)
        Coordinates where the function is centered (i.e., distribution mean). If None, center at the origin.
    second_derivative : int, optional (default -1)
        If greater or equal than zero, this is the index of the dimension along which to compute the second partial derivative instead of the pure density.

    Returns
    -------
    v : (nx,) numpy double array
        Values of compacly supported approximate Gaussian function at x

    Examples
    --------
    ```python
    from gpytoolbox import compactly_supported_normal
    import matplotlib.pyplot as plt
    x = np.reshape(np.linspace(-4,4,1000),(-1,1))
    plt.plot(x,compactly_supported_normal(x, n=4,center=np.array([0.5])))
    plt.plot(x,compactly_supported_normal(x, n=3,center=np.array([0.5])))
    plt.plot(x,compactly_supported_normal(x, n=2,center=np.array([0.5])))
    plt.plot(x,compactly_supported_normal(x, n=1,center=np.array([0.5])))
    plt.show()
    ```
    """

    ndim = x.ndim
    if ndim==0:
        x = np.array([[x]])
    if ndim==1:
        x = np.reshape(x,(-1,1))
    # print(ndim)
    dim = x.shape[1]

    if (center is None):
        center = np.zeros(dim)
    x_new = (x - np.tile(center[None,:],(x.shape[0],1)))/sigma

    # This is the function centered at 0 and with sigma=1
    def compactly_supported_normal_centered(x,n,second_derivative):
        vals = np.ones(x.shape[0])
        # I got this by writing into Wolfram Alpha:
        # Convolve[rect(y),rect(y)]
        # Convolve[rect(y),triangle(y)]
        # Convolve[rect(y),Piecewise[{{1/6 (-abs(y)^3 + 6 y^2 - 12 abs(y) + 8),1<=abs(y)<=2},{1/6 (3abs(y)^3 - 6y^2 + 4 ),abs(y)<1}}]]
        for i in range(dim):
            xx = x[:,i]
            if second_derivative==i:
                if n==2:
                    xx_val = (np.abs(xx)<0.5)*(-2.0) + (np.abs(xx)<1.5)*(np.abs(xx)>0.5)*(1.0)    
                elif n==1:
                    xx_val = 0.0*np.abs(xx)
                elif n==3:
                    xx_val = (np.abs(xx)<1)*(3*3*2*(np.abs(xx)) - 12)*(1/6)  +  (np.abs(xx)<=2)*(np.abs(xx)>=1)*( -3*2*(np.abs(xx)) + 12)*(1/6.)
                elif n==4:
                    xx_val = (np.abs(xx)<=2.5)*(np.abs(xx)>=1.5)*( 16*4*3*(np.abs(xx)**2.0) - 160*3*2*(np.abs(xx)) + 600*2)*(1/384) + (np.abs(xx)<=1.5)*(np.abs(xx)>=0.5)*( -16*4*3*(np.abs(xx)**2.0) + 80*3*2*(np.abs(xx)) - 120*2)*(1/96) + (np.abs(xx)<=0.5)*( 48*4*3*(np.abs(xx)**2.0) - 120*2)*(1/192)
            else:
                if n==2:
                    xx_val = (np.abs(xx)<0.5)*(-(xx**2.0)+0.75) + (np.abs(xx)<1.5)*(np.abs(xx)>0.5)*(0.5*(xx**2.0) - 1.5*np.abs(xx) + 1.125)
                elif n==1:
                    xx_val = (np.abs(xx)<1.0)*(1 - np.abs(xx))
                elif n==3:
                    xx_val = (np.abs(xx)<1)*(3*(np.abs(xx)**3.0) - 6*(np.abs(xx)**2.0) + 4)*(1/6)  +  (np.abs(xx)<=2)*(np.abs(xx)>=1)*( -(np.abs(xx)**3.0) + 6*(np.abs(xx)**2.0) - 12*np.abs(xx) + 8 )*(1/6.)
                elif n==4:
                    xx_val = (np.abs(xx)<=2.5)*(np.abs(xx)>=1.5)*( 16*(np.abs(xx)**4.0) - 160*(np.abs(xx)**3.0) + 600*(np.abs(xx)**2.0) - 1000*(np.abs(xx)) + 625)*(1/384) + (np.abs(xx)<=1.5)*(np.abs(xx)>=0.5)*( -16*(np.abs(xx)**4.0) + 80*(np.abs(xx)**3.0) - 120*(np.abs(xx)**2.0) + 20*(np.abs(xx)) + 55)*(1/96) + (np.abs(xx)<=0.5)*( 48*(np.abs(xx)**4.0) - 120*(np.abs(xx)**2.0) + 115)*(1/192)

            vals = vals * xx_val
        return vals

    val = compactly_supported_normal_centered(x_new,n,second_derivative)/(sigma**dim)
    if second_derivative>=0:
        val = val/(sigma**2.0)

    # We call it with the transformed x and rescale by sigma
    return val