diff --git a/docs/nodes/surface/curvatures.rst b/docs/nodes/surface/curvatures.rst new file mode 100644 index 0000000000000000000000000000000000000000..d5cc00c51a0d5234f7bd92a7f09c25a51dbd6228 --- /dev/null +++ b/docs/nodes/surface/curvatures.rst @@ -0,0 +1,127 @@ +Surface Curvatures +================== + +Functionality +------------- + +This node calculates several types of information about surface curvature: + +* Principal curvature values +* Principal curvature directions +* Gauss curvature +* Mean curvature +* Matrix based on principal curvature directions and surface normal. + +You can refer to Wikipedia_ for more detailed information about these terms. + +.. _Wkikpedia: https://en.wikipedia.org/wiki/Differential_geometry_of_surfaces + +If you need only Gaussian curvature value, you can as well use simpler "Surface Gauss Curvature" node. + +Inputs +------ + +This node has the following inputs: + +* **Surface**. The surface to analyze. This input is mandatory. +* **U**, **V**. Values of U and V surface parameters. These inputs are + available only when **Input mode** parameter isset to **Separate**. The + default value is 0.5. +* **UVPoints**. Points at which the surface is to be analyzed. Only two of + three coordinates will be used; the coordinates used are defined by the + **Orientation** parameter. This input is available and mandatory if the + **Input mode** parameter is set to **Vertices**. + +Parameters +---------- + +This node has the following parameters: + +* **Input mode**. The available options are: + + * **Separate**. The values of U and V surface parameters will be provided in + **U** and **V** inputs, correspondingly. + * **Vertices**. The values of U and V surface parameters will be provided in + **Vertices** input; only two of three coordinates of the input vertices + will be used. + + The default mode is **Separate**. + +* **Input orientation**. This parameter is available only when **Input mode** + parameter is set to **Vertices**. This defines which coordinates of vertices + provided in the **Vertices** input will be used. The available options are + XY, YZ and XZ. For example, if this is set to XY, then the X coordinate of + vertices will be used as surface U parameter, and Y coordinate will be used + as V parameter. The default value is XY. +* **Clamp**. This defines how the node will process the values of + surface U and V parameters which are out of the surface's domain. The + available options are: + + * **As is**. Do not do anything special, just pass the parameters to the + surface calculation algorithm as they are. The behaviour of the surface + when the values of parameters provided are out of domain depends on + specific surface: some will just calculate points by the same formula, + others will give an error. + * **Clamp**. Restrict the parameter values to the surface domain: replace + values that are greater than the higher bound with higher bound value, + replace values that are smaller than the lower bound with the lower bound + value. For example, if the surface domain along U direction is `[0 .. 1]`, + and the value of U parameter is 1.05, calculate the point of the surface + at U = 1.0. + * **Wrap**. Wrap the parameter values to be within the surface domain, i.e. + take the values modulo domain. For example, if the surface domain along U + direction is `[0 .. 1]`, and the value of U parameter is 1.05, evaluate + the surface at U = 0.05. + + The default mode is **As is**. + +* **Sort curvatures**. If checked, then the node will always make sure that + curvature values are calculated so that **Curvature1** value is less than + **Curvature2**. If not checked, this will not be guaranteed. For many simple + surfaces, the order of curvature values calculated without this parameter is + always the same. Checked by default. + +Outputs +------- + +This node has the following outputs: + +* **Curvature1**. The first principal curvature value. If **Sort curvatures** + parameter is checked, then this will be always the smaller principal + curvature value. +* **Curvature2**. The second principal curvature value. If **Sort curvatures** + parameter is checked, then this will be always the bigger principal curvature + value. +* **Dir1**. The first principal curvature direction - one which corresponds to **Curvature1**. +* **Dir2**. The second principal curvature direction - one which corresponds to **Curvature2**. +* **Gauss*. Gauss curvature value. +* **Mean**. Mean curvature value. +* **Matrix**. A matrix composed from principal curvature directions. It's X + axis is looking along **Dir1**, Y axis is looking along **Dir2** and Z axis + is looking along the surface's normal. + +Examples of usage +----------------- + +Use first principal curvature value for vertex colors: + +.. image:: https://user-images.githubusercontent.com/284644/80917635-d9089900-8d79-11ea-982e-ccde3742ffc6.png + +The same with second principal curvature value: + +.. image:: https://user-images.githubusercontent.com/284644/80917636-d9a12f80-8d79-11ea-9563-a77c3447abf5.png + +Gaussian curvature: + +.. image:: https://user-images.githubusercontent.com/284644/80917638-da39c600-8d79-11ea-8243-d74ab1f7b7b5.png + +Use Matrix output to place Suzannes on the surface: + +.. image:: https://user-images.githubusercontent.com/284644/80917634-d8700280-8d79-11ea-99e3-0d4d065d639a.png + +Here Suzannes are looking along the second principal curvature direction. + +Use Matrix output to place cubes on the surface: + +.. image:: https://user-images.githubusercontent.com/284644/80917640-da39c600-8d79-11ea-8e5e-2cfd3e7a0806.png + diff --git a/docs/nodes/surface/gauss_curvature.rst b/docs/nodes/surface/gauss_curvature.rst index 6cb9c9775da3226b03e21360dbba8f173b686509..533b16f953bad3ab543497f865b390261937fd09 100644 --- a/docs/nodes/surface/gauss_curvature.rst +++ b/docs/nodes/surface/gauss_curvature.rst @@ -17,6 +17,8 @@ Here are the main facts about Gaussian curvature that are usually useful: Note that the calculation is done by numerical differentiation, so it may be not very precise in some cases. +This node is a simpler version of the "Surface Curvatures" node. + Inputs ------ diff --git a/docs/nodes/surface/surface_index.rst b/docs/nodes/surface/surface_index.rst index 76f2ed431932dfed1a971e7b3f87245074ec55e1..c4c6f02c11a53c048c7100749837ab2fce8fa680 100644 --- a/docs/nodes/surface/surface_index.rst +++ b/docs/nodes/surface/surface_index.rst @@ -21,6 +21,7 @@ Surface subdomain flip swap + curvatures gauss_curvature tessellate_trim evaluate_surface diff --git a/index.md b/index.md index 347e7ff3a4b712d26b092f95a375c65ef5136ea6..8de9b56db9106b257fbb99563d34d5acd5873e62 100644 --- a/index.md +++ b/index.md @@ -96,6 +96,7 @@ SvFlipSurfaceNode SvSwapSurfaceNode SvSurfaceGaussCurvatureNode + SvSurfaceCurvaturesNode SvExApplyFieldToSurfaceNode SvExTessellateTrimSurfaceNode SvExEvalSurfaceNode diff --git a/nodes/surface/curvatures.py b/nodes/surface/curvatures.py new file mode 100644 index 0000000000000000000000000000000000000000..adca0b2c560e3fc7f8a8d3d2759e883a67f633e5 --- /dev/null +++ b/nodes/surface/curvatures.py @@ -0,0 +1,198 @@ + +import numpy as np + +import bpy +from bpy.props import FloatProperty, EnumProperty, BoolProperty, IntProperty + +import sverchok +from sverchok.node_tree import SverchCustomTreeNode, throttled +from sverchok.data_structure import updateNode, zip_long_repeat, ensure_nesting_level, get_data_nesting_level, repeat_last_for_length +from sverchok.utils.logging import info, exception +from sverchok.utils.surface import SvSurface + +class SvSurfaceCurvaturesNode(bpy.types.Node, SverchCustomTreeNode): + """ + Triggers: Surface Curvature + Tooltip: Calculate surface curvature values and directions + """ + bl_idname = 'SvSurfaceCurvaturesNode' + bl_label = 'Surface Curvature' + bl_icon = 'OUTLINER_OB_EMPTY' + sv_icon = 'SV_EVAL_SURFACE' + + @throttled + def update_sockets(self, context): + self.inputs['U'].hide_safe = self.input_mode == 'VERTICES' + self.inputs['V'].hide_safe = self.input_mode == 'VERTICES' + self.inputs['UVPoints'].hide_safe = self.input_mode == 'PAIRS' + + input_modes = [ + ('PAIRS', "Separate", "Separate U V (or X Y) sockets", 0), + ('VERTICES', "Vertices", "Single socket for vertices", 1) + ] + + input_mode : EnumProperty( + name = "Input mode", + items = input_modes, + default = 'PAIRS', + update = update_sockets) + + clamp_modes = [ + ('NO', "As is", "Do not clamp input values - try to process them as is (you will get either error or extrapolation on out-of-bounds values, depending on specific surface type", 0), + ('CLAMP', "Clamp", "Clamp input values into bounds - for example, turn -0.1 into 0", 1), + ('WRAP', "Wrap", "Wrap input values into bounds - for example, turn -0.1 into 0.9", 2) + ] + + clamp_mode : EnumProperty( + name = "Clamp", + items = clamp_modes, + default = 'NO', + update = updateNode) + + axes = [ + ('XY', "X Y", "XOY plane", 0), + ('YZ', "Y Z", "YOZ plane", 1), + ('XZ', "X Z", "XOZ plane", 2) + ] + + orientation : EnumProperty( + name = "Orientation", + items = axes, + default = 'XY', + update = updateNode) + + u_value : FloatProperty( + name = "U", + description = "Surface U parameter", + default = 0.5, + update = updateNode) + + v_value : FloatProperty( + name = "V", + description = "Surface V parameter", + default = 0.5, + update = updateNode) + + order : BoolProperty( + name = "Sort curvatures", + description = "If enabled, make sure that Curvature1 is always less than Curvature2", + default = True, + update = updateNode) + + def draw_buttons(self, context, layout): + layout.label(text="Input mode:") + layout.prop(self, "input_mode", expand=True) + if self.input_mode == 'VERTICES': + layout.label(text="Input orientation:") + layout.prop(self, "orientation", expand=True) + layout.prop(self, 'clamp_mode', expand=True) + layout.prop(self, 'order') + + def sv_init(self, context): + self.inputs.new('SvSurfaceSocket', "Surface") + self.inputs.new('SvVerticesSocket', "UVPoints") + self.inputs.new('SvStringsSocket', "U").prop_name = 'u_value' + self.inputs.new('SvStringsSocket', "V").prop_name = 'v_value' + self.outputs.new('SvStringsSocket', "Curvature1") + self.outputs.new('SvStringsSocket', "Curvature2") + self.outputs.new('SvVerticesSocket', "Dir1") + self.outputs.new('SvVerticesSocket', "Dir2") + self.outputs.new('SvStringsSocket', "Gauss") + self.outputs.new('SvStringsSocket', "Mean") + self.outputs.new('SvMatrixSocket', "Matrix") + self.update_sockets(context) + + def parse_input(self, verts): + verts = np.array(verts) + if self.orientation == 'XY': + us, vs = verts[:,0], verts[:,1] + elif self.orientation == 'YZ': + us, vs = verts[:,1], verts[:,2] + else: # XZ + us, vs = verts[:,0], verts[:,2] + return us, vs + + def _clamp(self, surface, us, vs): + u_min = surface.get_u_min() + u_max = surface.get_u_max() + v_min = surface.get_v_min() + v_max = surface.get_v_max() + us = np.clip(us, u_min, u_max) + vs = np.clip(vs, v_min, v_max) + return us, vs + + def _wrap(self, surface, us, vs): + u_min = surface.get_u_min() + u_max = surface.get_u_max() + v_min = surface.get_v_min() + v_max = surface.get_v_max() + u_len = u_max - u_min + v_len = v_max - u_min + us = (us - u_min) % u_len + u_min + vs = (vs - v_min) % v_len + v_min + return us, vs + + def process(self): + if not any(socket.is_linked for socket in self.outputs): + return + + surfaces_s = self.inputs['Surface'].sv_get() + surfaces_s = ensure_nesting_level(surfaces_s, 2, data_types=(SvSurface,)) + src_point_s = self.inputs['UVPoints'].sv_get(default=[[]]) + src_point_s = ensure_nesting_level(src_point_s, 4) + src_u_s = self.inputs['U'].sv_get() + src_u_s = ensure_nesting_level(src_u_s, 3) + src_v_s = self.inputs['V'].sv_get() + src_v_s = ensure_nesting_level(src_v_s, 3) + + need_values = self.outputs['Curvature1'].is_linked or self.outputs['Curvature2'].is_linked + need_directions = self.outputs['Dir1'].is_linked or self.outputs['Dir2'].is_linked + need_gauss = self.outputs['Gauss'].is_linked + need_mean = self.outputs['Mean'].is_linked + need_matrix = self.outputs['Matrix'].is_linked + + c1_out = [] + c2_out = [] + dir1_out = [] + dir2_out = [] + mean_out = [] + gauss_out = [] + matrix_out = [] + for surfaces, src_points_i, src_u_i, src_v_i in zip_long_repeat(surfaces_s, src_point_s, src_u_s, src_v_s): + new_curvatures = [] + for surface, src_points, src_us, src_vs in zip_long_repeat(surfaces, src_points_i, src_u_i, src_v_i): + if self.input_mode == 'VERTICES': + us, vs = self.parse_input(src_points) + else: + maxlen = max(len(src_us), len(src_vs)) + src_us = repeat_last_for_length(src_us, maxlen) + src_vs = repeat_last_for_length(src_vs, maxlen) + us, vs = np.array(src_us), np.array(src_vs) + data = surface.curvature_calculator(us, vs, order=self.order).calc(need_values, need_directions, need_gauss, need_mean, need_matrix) + if need_values: + c1_out.append(data.principal_value_1.tolist()) + c2_out.append(data.principal_value_2.tolist()) + if need_directions: + dir1_out.append(data.principal_direction_1.tolist()) + dir2_out.append(data.principal_direction_2.tolist()) + if need_mean: + mean_out.append(data.mean.tolist()) + if need_gauss: + gauss_out.append(data.gauss.tolist()) + if need_matrix: + matrix_out.extend(data.matrix) + + self.outputs['Curvature1'].sv_set(c1_out) + self.outputs['Curvature2'].sv_set(c2_out) + self.outputs['Dir1'].sv_set(dir1_out) + self.outputs['Dir2'].sv_set(dir2_out) + self.outputs['Mean'].sv_set(mean_out) + self.outputs['Gauss'].sv_set(gauss_out) + self.outputs['Matrix'].sv_set(matrix_out) + +def register(): + bpy.utils.register_class(SvSurfaceCurvaturesNode) + +def unregister(): + bpy.utils.unregister_class(SvSurfaceCurvaturesNode) + diff --git a/utils/surface.py b/utils/surface.py index dec6c275a2042c85e3351a62653a7e4a2d5645d1..b91dee451fdbf1334e300900d270c18ab1068214 100644 --- a/utils/surface.py +++ b/utils/surface.py @@ -43,6 +43,203 @@ def rotate_vector_around_vector_np(v, k, theta): s3 = p1 * p2 * k return s1 + s2 + s3 +class SurfaceCurvatureData(object): + """Container class for calculated curvature values""" + def __init__(self): + self.principal_value_1 = self.principal_value_2 = None + self.principal_direction_1 = self.principal_direction_2 = None + self.mean = self.gauss = None + self.matrix = None + +class SurfaceCurvatureCalculator(object): + """ + This class contains pre-calculated first and second surface derivatives, + and calculates any curvature information from them. + """ + def __init__(self, us, vs, order=True): + self.us = us + self.vs = vs + self.order = order + self.fu = self.fv = None + self.duu = self.dvv = self.duv = None + self.nuu = self.nvv = self.nuv = None + self.points = None + self.normals = None + + def set(self, points, normals, fu, fv, duu, dvv, duv, nuu, nvv, nuv): + """Set derivatives information""" + self.points = points + self.normals = normals + self.fu = fu # df/du + self.fv = fv # df/dv + self.duu = duu # (fu, fv), a.k.a. E + self.dvv = dvv # (fv, fv), a.k.a. G + self.duv = duv # (fu, fv), a.k.a F + self.nuu = nuu # (fuu, normal), a.k.a l + self.nvv = nvv # (fvv, normal), a.k.a n + self.nuv = nuv # (fuv, normal), a.k.a m + + def mean(self): + """Calculate mean curvature""" + duu, dvv, duv, nuu, nvv, nuv = self.duu, self.dvv, self.duv, self.nuu, self.nvv, self.nuv + A = duu*dvv - duv*duv + B = duu*nvv - 2*duv*nuv + dvv*nuu + return -B / (2*A) + + def gauss(self): + """Calculate Gaussian curvature""" + duu, dvv, duv, nuu, nvv, nuv = self.duu, self.dvv, self.duv, self.nuu, self.nvv, self.nuv + numerator = nuu * nvv - nuv*nuv + denominator = duu * dvv - duv*duv + return numerator / denominator + + def values(self): + """ + Calculate two principal curvature values. + If "order" parameter is set to True, then it will be guaranteed, + that C1 value is always less than C2. + """ + # It is possible to calculate principal curvature values + # as solutions of quadratic equation, without calculating + # corresponding principal curvature directions. + + # lambda^2 (E G - F^2) - lambda (E N - 2 F M + G L) + (L N - M^2) = 0 + + duu, dvv, duv, nuu, nvv, nuv = self.duu, self.dvv, self.duv, self.nuu, self.nvv, self.nuv + A = duu*dvv - duv*duv + B = duu*nvv - 2*duv*nuv + dvv*nuu + C = nuu*nvv - nuv*nuv + D = B*B - 4*A*C + c1 = (-B - np.sqrt(D))/(2*A) + c2 = (-B + np.sqrt(D))/(2*A) + + c1mask = (c1 < c2) + c2mask = np.logical_not(c1mask) + + c1_r = np.where(c1mask, c1, c2) + c2_r = np.where(c2mask, c1, c2) + + return c1_r, c2_r + + def values_and_directions(self): + """ + Calculate principal curvature values together with principal curvature directions. + If "order" parameter is set to True, then it will be guaranteed, that C1 value + is always less than C2. Curvature directions are always output correspondingly, + i.e. principal_direction_1 corresponds to principal_value_1 and principal_direction_2 + corresponds to principal_value_2. + """ + # If we need not only curvature values, but principal curvature directions as well, + # we have to solve an eigenvalue problem to find values and directions at once. + + # L p = lambda G p + + fu, fv = self.fu, self.fv + duu, dvv, duv, nuu, nvv, nuv = self.duu, self.dvv, self.duv, self.nuu, self.nvv, self.nuv + n = len(self.us) + + L = np.empty((n,2,2)) + L[:,0,0] = nuu + L[:,0,1] = nuv + L[:,1,0] = nuv + L[:,1,1] = nvv + + G = np.empty((n,2,2)) + G[:,0,0] = duu + G[:,0,1] = duv + G[:,1,0] = duv + G[:,1,1] = dvv + + M = np.matmul(np.linalg.inv(G), L) + eigvals, eigvecs = np.linalg.eig(M) + # Values of first and second principal curvatures + c1 = eigvals[:,0] + c2 = eigvals[:,1] + + if self.order: + c1mask = (c1 < c2) + c2mask = np.logical_not(c1mask) + c1_r = np.where(c1mask, c1, c2) + c2_r = np.where(c2mask, c1, c2) + else: + c1_r = c1 + c2_r = c2 + + # dir_1 corresponds to c1, dir_2 corresponds to c2 + dir_1_x = eigvecs[:,0,0][np.newaxis].T + dir_2_x = eigvecs[:,0,1][np.newaxis].T + dir_1_y = eigvecs[:,1,0][np.newaxis].T + dir_2_y = eigvecs[:,1,1][np.newaxis].T + + # another possible approach +# A = duv * nvv - dvv*nuv +# B = duu * nvv - dvv*nuu +# C = duu*nuv - duv*nuu +# D = B*B - 4*A*C +# t1 = (-B - np.sqrt(D)) / (2*A) +# t2 = (-B + np.sqrt(D)) / (2*A) + + dir_1 = dir_1_x * fu + dir_1_y * fv + dir_2 = dir_2_x * fu + dir_2_y * fv + + dir_1 = dir_1 / np.linalg.norm(dir_1, axis=1, keepdims=True) + dir_2 = dir_2 / np.linalg.norm(dir_2, axis=1, keepdims=True) + + if self.order: + c1mask = c1mask[np.newaxis].T + c2mask = c2mask[np.newaxis].T + dir_1_r = np.where(c1mask, dir_1, -dir_2) + dir_2_r = np.where(c2mask, dir_1, dir_2) + else: + dir_1_r = dir_1 + dir_2_r = dir_2 + #r = (np.cross(dir_1_r, dir_2_r) * self.normals).sum(axis=1) + #print(r) + + return c1_r, c2_r, dir_1_r, dir_2_r + + def calc(self, need_values=True, need_directions=True, need_gauss=True, need_mean=True, need_matrix = True): + """ + Calculate curvature information. + Return value: SurfaceCurvatureData instance. + """ + # We try to do as less calculations as possible, + # by not doing complex computations if not required + # and reusing results of other computations if possible. + data = SurfaceCurvatureData() + if need_matrix: + need_directions = True + if need_directions: + # If we need principal curvature directions, then the method + # being used will calculate us curvature values for free. + c1, c2, dir1, dir2 = self.values_and_directions() + data.principal_value_1, data.principal_value_2 = c1, c2 + data.principal_direction_1, data.principal_direction_2 = dir1, dir2 + if need_gauss: + data.gauss = c1 * c2 + if need_mean: + data.mean = (c1 + c2)/2.0 + if need_matrix: + matrices_np = np.dstack((data.principal_direction_2, data.principal_direction_1, self.normals)) + matrices_np = np.transpose(matrices_np, axes=(0,2,1)) + matrices_np = np.linalg.inv(matrices_np) + matrices = [Matrix(m.tolist()).to_4x4() for m in matrices_np] + for matrix, point in zip(matrices, self.points): + matrix.translation = Vector(point) + data.matrix = matrices + if need_values and not need_directions: + c1, c2 = self.values() + data.principal_value_1, data.principal_value_2 = c1, c2 + if need_gauss: + data.gauss = c1 * c2 + if need_mean: + data.mean = (c1 + c2)/2.0 + if need_gauss and not need_directions and not need_values: + data.gauss = self.gauss() + if need_mean and not need_directions and not need_values: + data.mean = self.mean() + return data + class SvSurface(object): def __repr__(self): if hasattr(self, '__description__'): @@ -84,7 +281,7 @@ class SvSurface(object): #self.info("Normals: %s", normal) return normal - def gauss_curvature_array(self, us, vs): + def curvature_calculator(self, us, vs, order=True): if hasattr(self, 'normal_delta'): h = self.normal_delta else: @@ -117,9 +314,25 @@ class SvSurface(object): dvv = np.linalg.norm(fv, axis=1) **2 duv = (fu * fv).sum(axis=1) - numerator = nuu * nvv - nuv*nuv - denominator = duu * dvv - duv*duv - return numerator / denominator + calc = SurfaceCurvatureCalculator(us, vs, order=order) + calc.set(surf_vertices, normal, fu, fv, duu, dvv, duv, nuu, nvv, nuv) + return calc + + def gauss_curvature_array(self, us, vs): + calc = self.curvature_calculator(us, vs) + return calc.gauss() + + def mean_curvature_array(self, us, vs): + calc = self.curvature_calculator(us, vs) + return calc.mean() + + def principal_curvature_values_array(self, us, vs): + calc = self.curvature_calculator(us, vs) + return calc.values() + + def principal_curvatures_array(self, us, vs): + calc = self.curvature_calculator(us, vs) + return calc.values_and_directions() def get_coord_mode(self): return 'UV'