Skip to content

marching_squares

marching_squares(S, GV, nx, ny)

Marching squares algorithm for extracting isocontours from a scalar field. Output is given as a list of (unordered) vertices and edge indices into the vertex list.

Parameters:

Name Type Description Default
S (nx

Scalar field

 required
GV (nx

Grid vertex positions

 required
nx int

Number of grid vertices in x direction

 required
ny int

Number of grid vertices in y direction

 required

Returns:

Name Type Description
V (nv,2) numpy double array

Vertex positions
E (ne,2) numpy int array

Edge indices into V

See Also

marching_cubes

Examples:

import numpy as np
import matplotlib.pyplot as plt
import gpytoolbox as gpt
# Create a scalar field
nx = 100
ny = 100
x = np.linspace(-1,1,nx)
y = np.linspace(-1,1,ny)
X,Y = np.meshgrid(x,y)
S = np.exp(-X**2-Y**2)
# Extract isocontours
V,E = gpt.marching_squares(S,np.c_[X.flatten(),Y.flatten()],nx,ny)
# Plot
plt.figure()
plt.imshow(S.reshape((nx,ny),order='F'),extent=[-1,1,-1,1])
for i in range(E.shape[0]):
        plt.plot([V[E[i,0],0],V[E[i,1],0]],
                 [V[E[i,0],1],V[E[i,1],1]],
                 'k-')
plt.show()
plt.axis('equal')
plt.show()
Source code in src/gpytoolbox/marching_squares.py 
  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
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def marching_squares(S,GV,nx,ny):
    """
    Marching squares algorithm for extracting isocontours from a scalar field. Output is given as a list of (unordered) vertices and edge indices into the vertex list.

    Parameters
    ----------
    S : (nx*ny,) numpy double array
        Scalar field
    GV : (nx*ny,2) numpy double array
        Grid vertex positions
    nx : int
        Number of grid vertices in x direction
    ny : int
        Number of grid vertices in y direction

    Returns
    -------
    V : (nv,2) numpy double array
        Vertex positions
    E : (ne,2) numpy int array
        Edge indices into V

    See Also
    --------
    marching_cubes

    Examples
    --------
    ```python
    import numpy as np
    import matplotlib.pyplot as plt
    import gpytoolbox as gpt
    # Create a scalar field
    nx = 100
    ny = 100
    x = np.linspace(-1,1,nx)
    y = np.linspace(-1,1,ny)
    X,Y = np.meshgrid(x,y)
    S = np.exp(-X**2-Y**2)
    # Extract isocontours
    V,E = gpt.marching_squares(S,np.c_[X.flatten(),Y.flatten()],nx,ny)
    # Plot
    plt.figure()
    plt.imshow(S.reshape((nx,ny),order='F'),extent=[-1,1,-1,1])
    for i in range(E.shape[0]):
            plt.plot([V[E[i,0],0],V[E[i,1],0]],
                     [V[E[i,0],1],V[E[i,1],1]],
                     'k-')
    plt.show()
    plt.axis('equal')
    plt.show()
    ```
    """
    S = np.reshape(S,(nx,ny),order='F')
    # Create empty list for
    verts = []
    edge_list = [] # index of edge vertices
    # Loop over all grid points
    for i in range(nx-1):
        for j in range(ny-1):
            # Get the scalar values at the corners of the grid cell
            a = S[i,j]
            b = S[i+1,j]
            c = S[i+1,j+1]
            d = S[i,j+1]
            # Get the contour index
            k = 0
            if a > 0:
                k += 1
            if b > 0:
                k += 2
            if c > 0:
                k += 4
            if d > 0:
                k += 8
            # Use symmetry
            flip = False
            if k > 7:
                flip = True
                k = 15 - k

            # Get the contour line segments
            if k == 1:
                # x = i
                x = i - a/(b-a)
                y = j
                verts.append([x,y])
                x = i
                y = j - a/(d-a)
                # y = j
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-2,len(verts)-1])
                else:
                    edge_list.append([len(verts)-1,len(verts)-2])
            elif k == 2:
                x = i - a/(b-a)
                # x = i
                y = j
                verts.append([x,y])
                x = i + 1
                y = j - b/(c-b)
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-1,len(verts)-2])
                else:
                    edge_list.append([len(verts)-2,len(verts)-1])
            elif k == 3:
                x = i
                y = j - a/(d-a)
                verts.append([x,y])
                x = i + 1
                y = j - b/(c-b)
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-1,len(verts)-2])
                else:
                    edge_list.append([len(verts)-2,len(verts)-1])
            elif k == 4:
                x = i + 1
                y = j- b/(c-b)
                verts.append([x,y])
                x = i - d/(c-d)
                y = j + 1
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-1,len(verts)-2])
                else:
                    edge_list.append([len(verts)-2,len(verts)-1])
            elif k == 5:
                x = i - a/(b-a)
                y = j
                verts.append([x,y])
                x = i
                y = j - a/(d-a)
                # y = j
                verts.append([x,y])
                x = i + 1
                y = j- b/(c-b)
                verts.append([x,y])
                x = i - d/(c-d)
                y = j + 1
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-4,len(verts)-3])
                else:
                    edge_list.append([len(verts)-3,len(verts)-4])
                if flip:
                    edge_list.append([len(verts)-1,len(verts)-2])
                else:
                    edge_list.append([len(verts)-2,len(verts)-1])
            elif k == 6:
                x = i - a/(b-a)
                y = j
                verts.append([x,y])
                x = i - d/(c-d)
                y = j + 1
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-1,len(verts)-2])
                else:
                    edge_list.append([len(verts)-2,len(verts)-1])
            elif k == 7:
                x = i - d/(c-d)
                y = j + 1
                verts.append([x,y])
                x = i
                y = j - a/(d-a)
                verts.append([x,y])
                if flip:
                    edge_list.append([len(verts)-2,len(verts)-1])
                else:
                    edge_list.append([len(verts)-1,len(verts)-2])
            else:
                pass

    # Convert list to numpy array
    verts = np.array(verts)
    edges = np.array(edge_list)
    verts, SVI, SVJ, edges = remove_duplicate_vertices(verts,faces=edges,
        epsilon=np.sqrt(np.finfo(verts.dtype).eps))

    # Remove trivial edges
    edges = edges[np.not_equal(edges[:,0], edges[:,1]), :]

    # # Remove duplicate edges
    # edges = np.unique(edges, axis=0)

    # Rescale to original grid
    verts[:,0] = verts[:,0]/(nx-1)
    verts[:,1] = verts[:,1]/(ny-1)
    verts = verts*(GV.max(axis=0)-GV.min(axis=0)) + GV.min(axis=0)

    return verts, edges