2D Dual Contouring

DualContouring2D.py

bl_info = {
    "name": "2D Implicit Function Contouring",
    "description": "Contour 2D implicit function",
    "author": "TheBusyTypist",
    "location": "View3D > Add > Mesh",
    "category": "Add Mesh"
}

import bpy
import bmesh
import numpy as np
import math

class Config:
    def __init__(self, begin=-2.0, end=2.0, ticks=20):
        self.begin = begin
        self.end = end
        self.ticks = ticks

def Sample(f, config):
    begin = config.begin
    end = config.end
    ticks = config.ticks

    x = np.linspace(begin, end, ticks)
    y = np.linspace(begin, end, ticks)
    X, Y = np.meshgrid(x, y)
    vF = np.vectorize(f)
    s = vF(X, Y)
    return s

def SolveIntersection(f, p0, p1, config):
    EPS = 1e-12
    m = (p0 + p1) * 0.5
    v = f(m[0], m[1])
    while abs(v) >= EPS:
        if v > 0:
            p1 = m
        else:
            p0 = m
        m = (p0 + p1) * 0.5
        v = f(m[0], m[1])
    return m

def HaveOppositeSigns(a, b):
    return (
        a < 0 and b >= 0
        or
        a >= 0 and b < 0
    )

def ComputeHermiteData(f, samples, config):
    def D(f, p):
        h = 1e-8
        fx1 = f(p[0] + h, p[1])
        fx0 = f(p[0] - h, p[1])
        fy1 = f(p[0], p[1] + h)
        fy0 = f(p[0], p[1] - h)
        n = np.array([(fx1 - fx0) * 0.5 / h, (fy1 - fy0) * 0.5 / h])
        return n

    n = config.ticks
    g = [None for i in range(n) for j in range(n)]
    begin = config.begin
    end = config.end
    step = (end - begin) / (n - 1)

    for y in range(n):
        for x in range(n):
            o = np.array([begin + x * step, begin + y * step])
            u = np.array([begin + (x + 1) * step, begin + y * step])
            v = np.array([begin + x * step, begin + (y + 1) * step])

            fo = samples[y, x]
            if x != n - 1:
                fu = samples[y, x + 1]
            else:
                fu = f(u[0], u[1])

            if y != n - 1:
                fv = samples[y + 1, x]
            else:
                fv = f(v[0], v[1])

            a = None
            if HaveOppositeSigns(fo, fu):
                if fo < fu:
                    a = SolveIntersection(f, o, u, config)
                else:
                    a = SolveIntersection(f, u, o, config)

            b = None
            if HaveOppositeSigns(fo, fv):
                if fo < fv:
                    b = SolveIntersection(f, o, v, config)
                else:
                    b = SolveIntersection(f, v, o, config)

            if a is not None:
                na = D(f, a)

            if b is not None:
                nb = D(f, b)

            h = (
                None if a is None else (a, na),
                None if b is None else (b, nb)
            )

            g[y * n + x] = h

    return g

def ConstructQEF(a):
    A = np.array([])
    b = []
    px = []
    py = []
    g = None

    if len(a) >= 2:
        for p, n in a:
            px.append(p[0])
            py.append(p[1])
            b.append(n.dot(p))
            if A.shape[0] == 0:
                A = n
            else:
                A = np.vstack((A, n))

        g = np.array([sum(px) / len(px), sum(py) / len(py)])

    return A, np.array(b), g
    
def SolveQEF(A, b, g):
    ATA = A.T.dot(A)
    ATb = A.T.dot(b)
    d = ATb - ATA.dot(g)
    c = np.linalg.pinv(ATA).dot(d)
    return c + g
    # return np.linalg.pinv(ATA).dot(ATb)

def IsInside(p, lx, ux, ly, uy):
    return p[0] <= ux and p[0] >= lx and p[1] <= uy and p[1] >= ly

def Contour(H, config):
    n = config.ticks - 1
    m = config.ticks
    v = [None for i in range(n) for j in range(n)]
    for y in range(n):
        for x in range(n):
            a = []
            
            h = H[y * m + x]
            h_up = H[(y + 1) * m + x]
            h_right = H[y * m + x + 1]

            if h[0] is not None:
                a.append(h[0])
            if h[1] is not None:
                a.append(h[1])

            if h_up[0] is not None:
                a.append(h_up[0])

            if h_right[1] is not None:
                a.append(h_right[1])

            A, b, g = ConstructQEF(a)

            if len(a) >= 2:
                p = SolveQEF(A, b, g)

                step = (config.end - config.begin) / n
                x0 = config.begin + step * x
                y0 = config.begin + step * y
                if not IsInside(p, x0, x0 + step, y0, y0 + step):
                    p = g
                    
                v[y * n + x] = (p[0], p[1])

    return v

def GenerateMesh(v, config, context):
    n = config.ticks - 1

    indices = [None for i in range(n) for j in range(n)]
    vertices = []
    for i, e in enumerate(v):
        if e is not None:
            indices[i] = len(vertices)
            vertices.append((e[0], e[1], 0))

    edges = []
    for y in range(n):
        for x in range(n):
            ci = y * n + x
            c = indices[ci]

            if x != n - 1:
                ri = y * n + x + 1
                r = indices[ri]

            if y != n - 1:
                ui = (y + 1) * n + x
                u = indices[ui]

            if c is not None:
                if x != n - 1 and r is not None:
                    edges.append((c, r))
                if y != n - 1 and u is not None:
                    edges.append((c, u))

    mesh = bpy.data.meshes.new("Mesh")
    mesh.from_pydata(vertices, edges, [])
    obj = bpy.data.objects.new("MeshObj", mesh)
    context.scene.objects.link(obj)

class DualContouring2D(bpy.types.Operator):
    """Contour 2D implicit function"""
    bl_idname = "object.contour_2d"
    bl_label = "Contour 2D implicit function"
    bl_options = {"REGISTER", "UNDO"}

    Ticks = bpy.props.IntProperty(name="Ticks", default=64)
    FunctionSource = bpy.props.StringProperty(name="Function",
        default="(x ** 2 + y ** 2 - 1) ** 3 - x ** 2 * y ** 3")

    def execute(self, context):
        f = eval("lambda x, y: " + self.FunctionSource)
        config = Config()
        config.ticks = self.Ticks
        samples = Sample(f, config)

        H = ComputeHermiteData(f, samples, config)

        v = Contour(H, config)

        GenerateMesh(v, config, context)
        
        return {"FINISHED"}

def register():
    bpy.utils.register_class(DualContouring2D)

def unregister():
    bpy.utils.unregister_class(DualContouring2D)

if __name__ == "__main__":
    register()