search

dhermes / bezier

Helper for Bézier Curves, Triangles, and Higher Order Objects
Apache License 2.0
262 stars 36 forks source link

Investigate bad surface-surface case(s) #121

Open dhermes opened 6 years ago

dhermes commented 6 years ago

These came up when I was randomly generating surfaces and intersecting them (and solving some PDEs along the way). There are (at least 4 cases).

Case 1

import bezier
import numpy as np
import seaborn; seaborn.set()

CASE1a = np.asfortranarray(
    [
        [
            -0.8214285714285716,
            -0.8396156993080822,
            -0.8576708331786456,
            -0.8217585564509393,
            -0.8398796873259764,
            -0.822088541473307,
        ],
        [
            1.0,
            0.9817980361507281,
            0.9637340006983077,
            0.9817980361507282,
            0.9636650364998821,
            0.963596072301456,
        ],
    ]
)
CASE1b = np.asfortranarray(
    [
        [
            -0.8348214285714287,
            -0.8158482142857144,
            -0.7968750000000001,
            -0.8348214285714287,
            -0.8158482142857143,
            -0.8348214285714286,
        ],
        [
            0.9866071428571429,
            1.0055803571428572,
            1.0245535714285716,
            1.0055803571428572,
            1.0245535714285716,
            1.0245535714285714,
        ],
    ]
)

surface1 = bezier.Surface(CASE1a, degree=2)
surface2 = bezier.Surface(CASE1b, degree=2)
intersection, = surface1.intersect(surface2)
ax = surface1.plot(256)
surface2.plot(256, ax=ax)
intersection.plot(256, ax=ax)
ax.axis('scaled')
ax.figure.show()

figure_1

dhermes commented 6 years ago

Case 2

CASE2a = np.asfortranarray(
    [
        [
            -0.8214285714285716,
            -0.8299520788750556,
            -0.8385565991179567,
            -0.830357142857143,
            -0.8389211567018355,
            -0.8392857142857144,
        ],
        [
            -1.0,
            -0.9912577228478906,
            -0.982478186840489,
            -1.0,
            -0.9912390934202446,
            -1.0,
        ],
    ]
)
CASE2b = np.asfortranarray(
    [
        [
            -0.8158482142857143,
            -0.8253348214285713,
            -0.8348214285714286,
            -0.8253348214285715,
            -0.8348214285714285,
            -0.8348214285714286,
        ],
        [
            -1.0055803571428572,
            -0.9960937500000001,
            -0.986607142857143,
            -1.0055803571428577,
            -0.9960937500000002,
            -1.0055803571428572,
        ],
    ]
)

surface1 = bezier.Surface(CASE2a, degree=2)
surface2 = bezier.Surface(CASE2b, degree=2)
intersection, = surface1.intersect(surface2)
ax = surface1.plot(256)
surface2.plot(256, ax=ax)
intersection.plot(256, ax=ax)
ax.axis('scaled')
ax.figure.show()

figure_1-1

dhermes commented 6 years ago

Case 3

CASE3a = np.asfortranarray(
    [
        [
            -0.9282473773857585,
            -0.9194403251518133,
            -0.9107142857142857,
            -0.919278299558978,
            -0.910511753723242,
            -0.910309221732198,
        ],
        [
            -0.9822918925640269,
            -0.9911645757096599,
            -1.0,
            -0.9823105219916728,
            -0.9911645757096597,
            -0.9823291514193193,
        ],
    ]
)
CASE3b = np.asfortranarray(
    [
        [
            -0.8917410714285715,
            -0.9012276785714285,
            -0.9107142857142858,
            -0.9012276785714286,
            -0.9107142857142858,
            -0.9107142857142859,
        ],
        [
            -1.0055803571428572,
            -0.9960937500000003,
            -0.986607142857143,
            -1.0055803571428572,
            -0.9960937500000002,
            -1.0055803571428572,
        ],
    ]
)

surface1 = bezier.Surface(CASE3a, degree=2)
surface2 = bezier.Surface(CASE3b, degree=2)
intersection, = surface1.intersect(surface2)
ax = surface1.plot(256)
surface2.plot(256, ax=ax)
intersection.plot(256, ax=ax)
ax.axis('scaled')
ax.figure.show()

figure_1-2

dhermes commented 6 years ago

Case 4

CASE4a = np.asfortranarray(
    [
        [
            -1.0,
            -0.9912364210826126,
            -0.9824398436629882,
            -1.0,
            -0.9912199218314941,
            -1.0,
        ],
        [
            0.8214285714285713,
            0.8301847323610781,
            0.838975375392798,
            0.8303571428571429,
            0.8391305448392564,
            0.8392857142857142,
        ],
    ]
)
CASE4b = np.asfortranarray(
    [
        [
            -1.0055803571428572,
            -0.9960937500000001,
            -0.986607142857143,
            -1.0055803571428577,
            -0.9960937500000002,
            -1.0055803571428572,
        ],
        [
            0.8158482142857142,
            0.8253348214285714,
            0.8348214285714285,
            0.825334821428571,
            0.8348214285714286,
            0.8348214285714285,
        ],
    ]
)

surface1 = bezier.Surface(CASE4a, degree=2)
surface2 = bezier.Surface(CASE4b, degree=2)
intersection, = surface1.intersect(surface2)
ax = surface1.plot(256)
surface2.plot(256, ax=ax)
intersection.plot(256, ax=ax)
ax.axis('scaled')
ax.figure.show()

figure_1-3

dhermes commented 6 years ago

Here are 16 more cases:

(Python script to generate image) ```python import bezier import matplotlib.pyplot as plt import numpy as np import seaborn CASES = { 5: ( np.asfortranarray([ [-0.8214285714285716, -0.8384755863215392, -0.8558466524001767, -0.8392857142857144, -0.8564947547715169, -0.8571428571428572], [-1.0, -0.9825154456957814, -0.9648818559703933, -1.0, -0.9824409279851967, -1.0], ]), np.asfortranarray([ [-0.8348214285714286, -0.8158482142857144, -0.7968750000000001, -0.8055676089602298, -0.7865943946745156, -0.7763137893490313], [-0.986607142857143, -1.0055803571428572, -1.0245535714285714, -0.9799395170043632, -0.9989127312900776, -0.9732718911515833], ]), ), 6: ( np.asfortranarray([ [0.9646385713472886, 0.9823633928063411, 1.0, 0.9646826784799853, 0.9823633928063411, 0.9647267856126821], [-0.8565120963510187, -0.8388914887908152, -0.8214285714285716, -0.8385761083948959, -0.8210343459336724, -0.8206401204387732], ]), np.asfortranarray([ [1.0245535714285714, 1.0055803571428574, 0.9866071428571428, 1.0245535714285716, 1.0055803571428572, 1.0245535714285714], [-0.7968750000000001, -0.8158482142857143, -0.8348214285714286, -0.8158482142857144, -0.8348214285714284, -0.8348214285714287], ]), ), 7: ( np.asfortranarray([ [1.0, 0.9823633928063411, 0.9646385713472886, 1.0, 0.9823192856736442, 1.0], [-0.8214285714285716, -0.8388914887908152, -0.8565120963510187, -0.8392857142857144, -0.8568274767469379, -0.8571428571428572], ]), np.asfortranarray([ [0.9866071428571428, 1.0055803571428572, 1.0245535714285714, 0.9836055999687541, 1.0025788142544687, 0.9806040570803654], [-0.8348214285714286, -0.8158482142857144, -0.7968750000000001, -0.8107796918031034, -0.791806477517389, -0.7867379550347781], ]), ), 8: ( np.asfortranarray([ [0.8214285714285713, 0.8395062499491981, 0.8574957142044314, 0.8216491070920552, 0.8396826784799852, 0.8218696427555392], [-1.0, -0.981748631647958, -0.9636549534938755, -0.981748631647958, -0.9635761083948957, -0.963497263295916], ]), np.asfortranarray([ [0.8348214285714285, 0.8158482142857142, 0.7968749999999999, 0.8348214285714286, 0.8158482142857142, 0.8348214285714285], [-0.9866071428571429, -1.0055803571428572, -1.0245535714285714, -1.0055803571428572, -1.0245535714285716, -1.0245535714285714], ]), ), 9: ( np.asfortranarray([ [-0.8385565991179567, -0.8299520788750556, -0.8214285714285716, -0.8295875212911764, -0.8210235074464837, -0.8206184434643962], [-0.982478186840489, -0.9912577228478905, -1.0, -0.9824968162681353, -0.9912577228478907, -0.9825154456957814], ]), np.asfortranarray([ [-0.8158482142857143, -0.8253348214285713, -0.8348214285714286, -0.8253348214285715, -0.8348214285714285, -0.8348214285714286], [-1.0055803571428572, -0.9960937500000001, -0.986607142857143, -1.0055803571428577, -0.9960937500000002, -1.0055803571428572], ]), ), 10: ( np.asfortranarray([ [-1.0, -0.9906663645893408, -0.9814137419750993, -0.990666364589341, -0.9813732355768907, -0.981332729178682], [-0.8214285714285716, -0.830543437133605, -0.8396210439833462, -0.8216148657050335, -0.8307111019824209, -0.8218011599814958], ]), np.asfortranarray([ [-0.9866071428571429, -0.9960937500000004, -1.0055803571428572, -0.9960937500000002, -1.0055803571428572, -1.0055803571428572], [-0.8348214285714287, -0.8253348214285716, -0.8158482142857144, -0.8348214285714288, -0.8253348214285716, -0.8348214285714287], ]), ), 11: ( np.asfortranarray([ [-0.9282473773857585, -0.9194403251518133, -0.9107142857142857, -0.919278299558978, -0.910511753723242, -0.910309221732198], [-0.9822918925640269, -0.9911645757096599, -1.0, -0.9823105219916728, -0.9911645757096597, -0.9823291514193193], ]), np.asfortranarray([ [-0.9296875000000001, -0.920200892857143, -0.9107142857142859, -0.9202008928571429, -0.9107142857142859, -0.9107142857142858], [-0.9866071428571431, -0.9960937500000002, -1.0055803571428572, -0.9866071428571431, -0.9960937500000003, -0.986607142857143], ]), ), 12: ( np.asfortranarray([ [0.8214285714285713, 0.8304674106888847, 0.8394841963828497, 0.8215388392603133, 0.8305666517374526, 0.8216491070920552], [-1.0, -0.990874315823979, -0.9817880541974479, -0.9908743158239788, -0.9817683429227031, -0.981748631647958], ]), np.asfortranarray([ [0.8348214285714285, 0.8253348214285715, 0.8158482142857142, 0.8348214285714286, 0.8253348214285714, 0.8348214285714285], [-0.9866071428571429, -0.9960937500000004, -1.0055803571428572, -0.9960937500000002, -1.0055803571428572, -1.0055803571428572], ]), ), 13: ( np.asfortranarray([ [0.9823413392399927, 0.9911816964031703, 1.0, 0.9823523660231668, 0.9911816964031706, 0.9823633928063411], [-0.8389309113403052, -0.8301600301096934, -0.8214285714285716, -0.8299826286369891, -0.821231458681122, -0.8210343459336724], ]), np.asfortranarray([ [1.0055803571428572, 0.99609375, 0.9866071428571428, 1.0055803571428572, 0.99609375, 1.0055803571428572], [-0.8158482142857143, -0.8253348214285713, -0.8348214285714286, -0.8253348214285715, -0.8348214285714285, -0.8348214285714286], ]), ), 14: ( np.asfortranarray([ [0.9822310714082507, 0.9911265624872994, 1.0, 0.9822420981914253, 0.9911265624872996, 0.9822531249745992], [-0.928413738373469, -0.9195443007691324, -0.9107142857142857, -0.9194654556701528, -0.9106157293405611, -0.9105171729668362], ]), np.asfortranarray([ [1.0055803571428572, 0.99609375, 0.9866071428571428, 1.0055803571428572, 0.99609375, 1.0055803571428572], [-0.8917410714285715, -0.9012276785714285, -0.9107142857142858, -0.9012276785714286, -0.9107142857142858, -0.9107142857142859], ]), ), 15: ( np.asfortranarray([ [0.9822310714082507, 0.9911265624872994, 1.0, 0.9822420981914253, 0.9911265624872996, 0.9822531249745992], [-0.928413738373469, -0.9195443007691324, -0.9107142857142857, -0.9194654556701528, -0.9106157293405611, -0.9105171729668362], ]), np.asfortranarray([ [0.9866071428571428, 0.99609375, 1.0055803571428572, 0.9866071428571428, 0.99609375, 0.9866071428571428], [-0.9296875000000001, -0.920200892857143, -0.9107142857142859, -0.9202008928571429, -0.9107142857142859, -0.9107142857142858], ]), ), 16: ( np.asfortranarray([ [0.9107142857142857, 0.9196979910587284, 0.9286596428368222, 0.9107694196301569, 0.9197420981914252, 0.9108245535460278], [-1.0, -0.9909728721977038, -0.9819851669448975, -0.9909728721977038, -0.9819654556701525, -0.9819457443954076], ]), np.asfortranarray([ [0.9107142857142856, 0.9012276785714286, 0.8917410714285713, 0.9107142857142858, 0.9012276785714286, 0.9107142857142857], [-0.986607142857143, -0.9960937500000001, -1.0055803571428574, -0.9960937500000002, -1.0055803571428574, -1.0055803571428572], ]), ), 17: ( np.asfortranarray([ [-0.8214285714285716, -0.830522135368327, -0.8395827008058454, -0.8215935639397555, -0.8306706286283924, -0.8217585564509393], [1.0, 0.9908990180753638, 0.981832518249941, 0.9908990180753641, 0.9818152772003346, 0.9817980361507281], ]), np.asfortranarray([ [-0.8348214285714287, -0.8253348214285716, -0.8158482142857144, -0.8348214285714288, -0.8253348214285716, -0.8348214285714287], [0.9866071428571429, 0.9960937500000004, 1.0055803571428572, 0.9960937500000002, 1.0055803571428577, 1.0055803571428572], ]), ), 18: ( np.asfortranarray([ [-1.0, -0.9912364210826126, -0.9824398436629882, -1.0, -0.9912199218314941, -1.0], [0.8214285714285713, 0.8301847323610781, 0.838975375392798, 0.8303571428571429, 0.8391305448392564, 0.8392857142857142], ]), np.asfortranarray([ [-0.9859141689170303, -0.9862606558870869, -0.986607142857143, -0.9957472630299438, -0.9960937500000002, -1.0055803571428572], [0.8131424774977621, 0.8239819530345953, 0.8348214285714285, 0.8144953458917382, 0.8253348214285715, 0.8158482142857142], ]), ), 19: ( np.asfortranarray([ [-1.0, -0.9911704240781387, -0.9823078496540409, -1.0, -0.9911539248270205, -1.0], [0.8928571428571428, 0.9016822679880755, 0.910541875218221, 0.9017857142857142, 0.9106280804662537, 0.9107142857142857], ]), np.asfortranarray([ [-0.9866071428571431, -0.9960937500000002, -1.0055803571428572, -0.9866071428571431, -0.9960937500000003, -0.986607142857143], [0.9296875, 0.9202008928571428, 0.9107142857142855, 0.920200892857143, 0.9107142857142854, 0.9107142857142856], ]), ), 20: ( np.asfortranarray([ [-0.9107142857142857, -0.9197253533984493, -0.9287034225803757, -0.9107967819698778, -0.9197913504029228, -0.9108792782254695], [1.0, 0.9909852233233962, 0.9820049287460056, 0.9909852233233966, 0.981987687696399, 0.9819704466467927], ]), np.asfortranarray([ [-0.9107142857142858, -0.9012276785714286, -0.8917410714285715, -0.9107142857142858, -0.9012276785714286, -0.9107142857142859], [0.9866071428571428, 0.99609375, 1.0055803571428572, 0.99609375, 1.0055803571428572, 1.0055803571428572], ]), ), } def main(): fig, axes = plt.subplots(4, 4) axes = axes.flatten() for i in range(5, 21): ax = axes[i - 5] nodes1, nodes2 = CASES[i] surface1 = bezier.Surface.from_nodes(nodes1) surface2 = bezier.Surface.from_nodes(nodes2) intersection, = surface1.intersect(surface2) surface1.plot(256, ax=ax) surface2.plot(256, ax=ax) intersection.plot(256, ax=ax) ax.axis('scaled') ax.set_title('Case {}'.format(i)) plt.show() if __name__ == '__main__': seaborn.set() main() ```

figure_1-1

dhermes commented 6 years ago

It's worth noting, these were "detected" by approximating each surface by a polygon (by just evaluating at 64 points along each edge, i.e. 0, 1/64, ..., 63/64) and then intersecting the polygons using Shapely.

I tried to be careful to figure out what type of intersection there was (i.e. Polygon vs. MultiPolygon vs. LineString etc.) but I'm fairly certain that accessing .area will work in a consistent way:

s_vals = np.linspace(0, 1, 65)[:64]

evaluated1 = np.empty((2, 192), order='F')
index = 0
for edge in surface1.edges:
    evaluated1[:, index:index + 64] = edge.evaluate_multi(s_vals)
    index += 64
polygon1 = shapely.geometry.Polygon(evaluated1.T)

evaluated2 = np.empty((2, 192), order='F')
index = 0
for edge in surface2.edges:
    evaluated2[:, index:index + 64] = edge.evaluate_multi(s_vals)
    index += 64
polygon2 = shapely.geometry.Polygon(evaluated2.T)

intersection = polygon1.intersection(polygon2)
print(intersection.area)