import bpy
import addon_utils
addon_utils.enable("add_curve_extra_objects")
import copy
from math import radians, pi
import mathutils
import torch
import numpy as np
import bmesh
import math
from mathutils import Vector, Quaternion
import itertools

USE_MINIMUM_FACE=True
AVERAGE_EDGE_LENGTH=None


"""
Delete all objects in the scene
"""
def delete_all():
    try:
        bpy.ops.object.select_all(action='SELECT')
        bpy.ops.object.delete()
    except:
        bpy.ops.wm.read_homefile()
        bpy.ops.object.select_all(action='SELECT')
        bpy.ops.object.delete()


"""
Change the center of the object
"""
def change_origin(name, location):
    bpy.context.scene.cursor.location = location
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.origin_set(type='ORIGIN_CURSOR', center='MEDIAN')

"""
Make sure the normal of the mesh is all facing outward or inside
"""
def recalculate_normals(name,inside=False):
    obj = bpy.data.objects[name]
    bpy.ops.object.select_all(action='DESELECT')
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    # go edit mode
    bpy.ops.object.mode_set(mode='EDIT')
    # select al faces
    bpy.ops.mesh.select_all(action='SELECT')
    # recalculate outside normals 
    bpy.ops.mesh.normals_make_consistent(inside=inside)
    # go object mode again
    bpy.ops.object.editmode_toggle()

"""
Extract the mesh points and polygons from the model generated by Blender
"""
def get_faces():
    #print("faces:",bpy.data.objects)
    names = bpy.context.scene.objects.keys()

    def get_one_obj_faces(name):  
        obj = bpy.data.objects[name]
        bpy.ops.object.select_all(action='DESELECT')
        obj.select_set(True)
        modifier = obj.modifiers.new(name="Tri", type='TRIANGULATE')
        bpy.context.view_layer.objects.active = obj
        bpy.ops.object.modifier_apply(modifier="Tri")
        obj_eval = obj.evaluated_get(depsgraph=bpy.context.evaluated_depsgraph_get())
        mesh = obj_eval.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get())
        bm = bmesh.new()
        bm.from_mesh(mesh)
        bm.transform(obj.matrix_world)
        bm.faces.ensure_lookup_table()

        vertices = []
        indices = []

        for f in bm.faces:
            indices.append([v.index for v in f.verts])
        for v in bm.verts:
            vertices.append(v.co)

        indices = np.array(indices)
        vertices = np.array(vertices)

        bm.free()

        return vertices, indices

    verts_all, faces_all= [], []
    cnt=0
    for name in names:
        if bpy.data.objects[name].type=='MESH':
            verts,faces=get_one_obj_faces(name)
            verts_all+=verts.tolist()
            faces+=cnt
            faces_all+=faces.tolist()
            cnt+=len(verts)

    verts_all = torch.tensor(np.array(verts_all),dtype=torch.float32).contiguous()
    faces_all = torch.tensor(np.array(faces_all),dtype=torch.long).contiguous()

    return verts_all,faces_all


"""
fill caps to mesh
"""
def make_caps(name,type="start"):
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.mode_set(mode = 'EDIT')

    object = bpy.data.objects[name]
    mesh = bmesh.from_edit_mesh(object.data) 
    mesh.verts.ensure_lookup_table()
    
    
    if type=="start":
        for i in range(len(mesh.verts)):
            if mesh.verts[i].is_boundary and i<len(mesh.verts)/2:
                mesh.verts[i].select=True
        bpy.ops.mesh.edge_face_add()
    elif type=="end":
        for i in range(len(mesh.verts)):
            if mesh.verts[i].is_boundary and i>len(mesh.verts)/2:
                mesh.verts[i].select=True
        bpy.ops.mesh.edge_face_add()
    
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.mode_set(mode = 'OBJECT')


def split_cone_z(name,n=3):
    # Get the currently active object (make sure it is a cone)
    cone = bpy.data.objects[name]
    if not cone or cone.type != 'MESH':
        raise Exception("Please select a cone object first")

    # Enter edit mode
    bpy.ops.object.mode_set(mode='EDIT')
    
    # Create BMesh instance
    bm = bmesh.from_edit_mesh(cone.data)
    bm.verts.ensure_lookup_table()
    bm.edges.ensure_lookup_table()

    # Find the top vertex (maximum Z coordinate)
    top_vert = max(bm.verts, key=lambda v: v.co.z)
    
    # Select all edges connected to the top vertex (vertical edges)
    vertical_edges = [edge for edge in bm.edges if top_vert in edge.verts]
    
    # Deselect all and select these edges
    bpy.ops.mesh.select_all(action='DESELECT')
    for edge in vertical_edges:
        edge.select = True

    # Perform subdivision (generate n-1 cuts)
    bmesh.ops.subdivide_edges(
        bm,
        edges=vertical_edges,
        cuts=n-1,         # Number of cuts = number of segments - 1
        use_grid_fill=True # Keep ring structure
    )

    # Update mesh and exit edit mode
    bmesh.update_edit_mesh(cone.data)
    bpy.ops.object.mode_set(mode='OBJECT')
    primitive = bpy.data.objects[name]
    
    return primitive
    
def subdivide_primitive(name, resolution, axes):
    # Get the currently active object
    obj = bpy.data.objects[name]
    
    if not obj or obj.type != 'MESH':
        raise Exception("Please select a mesh object")

    # Enter edit mode
    bpy.ops.object.mode_set(mode='EDIT')
    
    # Get BMesh
    mesh = bmesh.from_edit_mesh(obj.data)
    
    def select_edges(axis):
        """Select edges along the specified axis"""
        for edge in mesh.edges:
            v1, v2 = edge.verts
            delta = v2.co - v1.co
            # Calculate absolute values of each component
            dx, dy, dz = (abs(delta.x), abs(delta.y), abs(delta.z))
            
            if axis == 'X' and dy < 1e-6 and dz < 1e-6 and dx > 1e-6:
                edge.select = True
            elif axis == 'Y' and dx < 1e-6 and dz < 1e-6 and dy > 1e-6:
                edge.select = True
            elif axis == 'Z' and dx < 1e-6 and dy < 1e-6 and dz > 1e-6:
                edge.select = True
    
    print(resolution,axes)
    for res,axis in zip(resolution,axes):
        bpy.ops.mesh.select_all(action='DESELECT')
        select_edges(axis)
        bpy.ops.mesh.subdivide(number_cuts=res-1, smoothness=0)
    
    # Update mesh and return to object mode
    bmesh.update_edit_mesh(obj.data)
    bpy.ops.object.mode_set(mode='OBJECT')
    primitive = bpy.data.objects[name]

    return primitive

"""
Add thickness to the face
"""
def solidify(name, thickness=0.01, offset=-1, boundary_mode=None):
    # apply transform to the object
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.mode_set(mode = 'OBJECT')
    bpy.data.objects[name].select_set(True)
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=True)
    bpy.data.objects[name].select_set(False)
    # add solidify modifier
    bpy.data.objects[name].modifiers.new("Solidify", "SOLIDIFY")
    bpy.data.objects[name].modifiers["Solidify"].thickness = thickness
    bpy.data.objects[name].modifiers["Solidify"].offset = offset
    if boundary_mode:
        bpy.data.objects[name].modifiers["Solidify"].solidify_mode = 'NON_MANIFOLD'
        bpy.data.objects[name].modifiers["Solidify"].nonmanifold_boundary_mode = boundary_mode
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.modifier_apply(modifier="Solidify")


"""
Creates a translational object of a line trajectory
"""
def create_curve(name, profile_name=None,control_points=[],points_radius=[],handle_type=[],closed=False, center="POINT",thickness=None, fill_caps="none",flip_normals=False, bevel_width=None, bevel_segments=8,use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH, resolution=24, volumn_origin=True):
    if isinstance(name, str):
        type_dict={0:"AUTO", 1:"VECTOR", 2:"ALIGNED", 3:"FREE"}

        control_points = np.array(control_points).tolist()
        control_points_tmp = copy.deepcopy(control_points)
        num_handle_co = handle_type.count(3)
        num_control_points = len(control_points) - num_handle_co

        curveData = bpy.data.curves.new(name, type='CURVE')
        curveData.dimensions = '3D'

        bezierSpline = curveData.splines.new('BEZIER')
        bezierSpline.bezier_points.add(num_control_points - 1)
        bezierSpline.use_cyclic_u = closed
        

        for i in range(num_control_points):
            bezier_point = bezierSpline.bezier_points[i]
            bezier_point.handle_left_type = type_dict[handle_type[2*i]]
            if type_dict[handle_type[2*i]]=="FREE":
                bezier_point.handle_left = control_points.pop(0)
            bezier_point.co = control_points.pop(0)
            bezier_point.handle_right_type = type_dict[handle_type[2*i+1]]
            if type_dict[handle_type[2*i+1]]=="FREE":
                bezier_point.handle_right = control_points.pop(0)
            bezier_point.radius = points_radius[i] if len(points_radius)!=0 else 1.0

        assert len(control_points)==0, "cannot create curve"
        if use_minimum_face:
            use_resolution = 12
        elif not average_edge_length is None:
            for i in range(len(bezierSpline.bezier_points) - 1):
                p1 = bezierSpline.bezier_points[i].co
                p2 = bezierSpline.bezier_points[i + 1].co
                total_length += (p2 - p1).length
            use_resolution = total_length/average_edge_length
        
        if resolution:
            use_resolution = resolution
        curveData.resolution_u = use_resolution
        
        curveOB = bpy.data.objects.new(name, curveData)

        if profile_name != None:
            curveData.bevel_mode = "OBJECT"
            curveData.splines[0].use_smooth = False
            scn = bpy.context.scene.collection
            scn.objects.link(curveOB)

            if bevel_width!=None:
                bevel(name=name, width=bevel_width, segments=bevel_segments)
                curveData = bpy.data.objects[name].data
                curveData.bevel_mode = 'OBJECT'

            
            curveData.bevel_object = bpy.data.objects[profile_name]
            if fill_caps=="both":
                curveData.use_fill_caps = True
            else:
                curveData.use_fill_caps = False
            
            if use_minimum_face:
                use_resolution = 24
            elif not average_edge_length is None:
                for i in range(len(bezierSpline.bezier_points) - 1):
                    p1 = bezierSpline.bezier_points[i].co
                    p2 = bezierSpline.bezier_points[i + 1].co
                    total_length += (p2 - p1).length
                use_resolution = total_length/average_edge_length
            
            if resolution:
                use_resolution = resolution

            curveData.resolution_u = use_resolution


            bpy.context.view_layer.objects.active = bpy.data.objects[name]
            bpy.ops.object.mode_set(mode = 'OBJECT')
            bpy.data.objects[name].select_set(True)
            bpy.ops.object.convert(target='MESH')
            bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
            if volumn_origin:
                bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')

            if fill_caps in ["start","end"]:
                make_caps(name,fill_caps)

            if flip_normals:
                recalculate_normals(name,inside=True)
            else:
                recalculate_normals(name,inside=False)

            if thickness>1e-10:
                solidify(name,thickness)

            weld(name,1e-5)

            return curveOB
        else:
            scn = bpy.context.scene.collection
            scn.objects.link(curveOB)
            if center=="MEDIAN":
                bpy.data.objects[name].select_set(True)
                bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY', center='MEDIAN')
            points=np.array(control_points_tmp)
            return {"name":name, "points":points, "handle_type":handle_type, "closed":closed, "center":center}
        
    elif profile_name==None:
        if isinstance(profile_name, str):
            profile_name = [profile_name]*len(name)
        if len(points_radius) != 0 and (isinstance(points_radius[0], float) or isinstance(points_radius[0], int)):
            points_radius = [points_radius]*len(name)
        elif len(points_radius) == 0:
            points_radius = [[]] * len(name)
        if isinstance(handle_type[0], int):
            handle_type = [handle_type]*len(name)
        if isinstance(closed, bool):
            closed = [closed]*len(name)
        if isinstance(center, str):
            center = [center]*len(name)
        for i in range(len(name)):
            create_curve(name=name[i], control_points=control_points[i], points_radius=points_radius[i], handle_type=handle_type[i], closed=closed[i], center=center[i], thickness=thickness, fill_caps=fill_caps, flip_normals=flip_normals, resolution=resolution)
        points = np.array(copy.deepcopy(control_points))
        return {"name":name, "points":points, "handle_type":handle_type, "closed":closed, "center":center}

def find_circle(points):
    """
    Calculate the center (cx, cy) and radius r of the circle passing through three points.
    Returns (cx, cy, r).
    """
    x1,y1 = points[0][0],points[0][1]
    x2,y2 = points[1][0],points[1][1]
    x3,y3 = points[2][0],points[2][1]
    D = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
    cx = ((x1**2 + y1**2) * (y2 - y3) + (x2**2 + y2**2) * (y3 - y1) + (x3**2 + y3**2) * (y1 - y2)) / (D+1e-10)
    cy = ((x1**2 + y1**2) * (x3 - x2) + (x2**2 + y2**2) * (x1 - x3) + (x3**2 + y3**2) * (x2 - x1)) / (D+1e-10)
    r = math.sqrt((cx - x1)**2 + (cy - y1)**2)
    return cx, cy, r

def calculate_angle(x1, y1, x2, y2):
    """
    Calculate the angle in radians between the positive x-axis and the line from (x1, y1) to (x2, y2).
    """
    dx = x2 - x1
    dy = y2 - y1
    return math.atan2(dy, dx)

"""
Create primitive object
"""
def create_primitive(name, primitive_type="cube", location=None, scale=None, rotation=None, rotation_mode='QUATERNION', apply=False, x_subdivisions=None, y_subdivisions=None, use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH, resolution=None):
    
    if primitive_type=="uv_sphere":
        if not use_minimum_face:
            if average_edge_length !=None:
                res=int(2*pi//average_edge_length)
                segments, ring_count=res, res    
            else:
                segments, ring_count=resolution[0], resolution[1]
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")(segments=segments,ring_count=ring_count)            
            primitive = bpy.context.object
            primitive.name = name
        else:
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()
            primitive = bpy.context.object
            primitive.name = name
    
    if primitive_type in ["cylinder","cone"]:
        if not use_minimum_face:
            if average_edge_length !=None:
                res_1=int(2*pi//average_edge_length)
                vertices=res_1
                res_2=int(2//average_edge_length)
            else:
                vertices,res_2=resolution[0], resolution[1]
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")(vertices=vertices)
            primitive = bpy.context.object
            primitive.name = name
            if primitive_type=="cylinder":
                primitive=subdivide_primitive(name,[res_2],['Z'])
            else:
                primitive=split_cone_z(name,res_2)           
        else:
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()
            primitive = bpy.context.object
            primitive.name = name
    
    if primitive_type=="torus":
        if not use_minimum_face:
            if average_edge_length !=None:
                res_1=int(2*pi//average_edge_length)
                res_2=int(0.5*pi//average_edge_length)
            else:
                res_1, res_2=resolution[0], resolution[1]
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")(major_segments=res_1,minor_segments=res_2)            
            primitive = bpy.context.object
            primitive.name = name
        else:
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()
            primitive = bpy.context.object
            primitive.name = name
    
    if primitive_type=="cube":
        if not use_minimum_face:
            if average_edge_length !=None:
                res=int(2//average_edge_length)
                if res>1:
                    resolution=[res,res,res]
            else:
                pass
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()  
            primitive = bpy.context.object
            primitive.name = name          
            primitive=subdivide_primitive(name,resolution,['X','Y','Z'])
        else:
            getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()
            primitive = bpy.context.object
            primitive.name = name
    
    if primitive_type=="grid":
        getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")(x_subdivisions=x_subdivisions,y_subdivisions=y_subdivisions)
        primitive = bpy.context.object
        primitive.name = name
    #getattr(bpy.ops.mesh, f"primitive_{primitive_type}_add")()


    if location:
        primitive.location = location
    if scale:
        primitive.scale = scale
    if rotation:
        if rotation_mode=='XYZ':
            primitive.rotation_euler = [angle * pi for angle in rotation]
        elif rotation_mode=='QUATERNION':
            primitive.rotation_mode = 'QUATERNION'
            primitive.rotation_quaternion = rotation
        elif rotation_mode=='MATRIX':
            mat = np.eye(4)
            rotation = np.array(rotation).reshape([3,3])
            mat[:3,:3] = rotation
            bpy.context.view_layer.update()
            world_matrix = torch.tensor(bpy.data.objects[name].matrix_world)
            scale_now = world_matrix.norm(dim=0)[:3]
            scale_matrix = torch.eye(4)
            scale_matrix[0,0],scale_matrix[1,1],scale_matrix[2,2] =scale_now[0],scale_now[1],scale_now[2]
            scale_matrix_inv = scale_matrix.clone()
            for i in range(3):
                # inverse matrix, if the diagonal element is 0, keep it as 0
                if scale_matrix_inv[i,i]>1e-10:
                    scale_matrix_inv[i,i]=1.0 / scale_matrix_inv[i,i]
            mat = scale_matrix_inv@torch.tensor(mat,dtype=torch.float32)@scale_matrix
            mat = mathutils.Matrix(np.array(mat))
            bpy.data.objects[name].matrix_world = bpy.data.objects[name].matrix_world@mat

    if apply:
        bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)

    return primitive

"""
Create primitive object
"""
def create_polygon(name: list[str] | str, sides: list[int] | int=3, location: list[list[float]] | list[float]=[0., 0., 0.], profile_name=None, rotation: list[list[float]]=[1., 0., 0., 0.], radius: list[float] | float=1., scale: list[list[float]] | list[float]=[1., 1., 1.], resolution: list[int] | int=12, center="MEDIAN", thickness=None, fill_caps="none"):
    """
    Create a regular polygon with given sides number and return the object's name and parameters
    Args:
        name (list[str] | str): Name(s) of the polygon object(s)
        sides (list[int] | int): Number of sides for the polygon(s), default is 3
        location (list[list[float]] | list[float]): Location(s) of the polygon(s), default is [0,0,0]
        profile_name (str, optional): Name of the profile object for the bevel. None means no bevel
        rotation (list[list[float]]): Rotation quaternion(s) of the polygon(s), default is [1,0,0,0] 
        radius (list[float] | float): Radius(es) of the polygon(s), default is 1.0
        scale (list[list[float]] | list[float]): Scale factor(s) for the polygon(s), default is [1,1,1]
        center (str): Center type, either "MEDIAN" or "POINT", default is "MEDIAN"
        thickness (float, optional): Thickness of the polygon if solidify is needed
        fill_caps (str): Type of cap filling - "none", "start", "end" or "both", default is "none"
        resolution (int): Resolution of the curve, default is 24
    Returns:
        dict: Dictionary containing polygon parameters (name, radius, center)
    """
    location = location.copy()
    rotation = rotation.copy()
    scale = scale.copy()
    if not isinstance(name, list):
        name = [name]
    if not isinstance(location[0], list):
        location = [location]
    # Name and location must have same length and each polygon needs unique name/location
    assert len(name) == len(location)

    # Convert single values to lists if needed
    if not isinstance(rotation[0], list):
        rotation = [rotation] * len(name)
    if not isinstance(radius, list):
        radius = [radius] * len(name)
    if not isinstance(scale[0], list):
        scale = [scale] * len(name)
    if not isinstance(sides, list):
        sides = [sides] * len(name)

    if profile_name is None:  
        # Create each polygon separately
        for i in range(len(name)):
            bpy.ops.curve.simple(align='WORLD', location=location[i], rotation=[0, 0, 0], Simple_Type='Polygon', shape='3D', Simple_sides=sides[i], Simple_radius=radius[i], outputType='BEZIER', use_cyclic_u=True, edit_mode=False)
            curve = bpy.context.object
            curve.name = name[i]   
            curve.rotation_mode = 'QUATERNION'
            curve.rotation_quaternion = rotation[i]
            if center=="POINT":
                curve.location = [radius[i],0,0]
                bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
            curveData = curve.data
            curveData.dimensions = '3D'
            bpy.data.objects[name[i]].scale = scale[i]
            
        return {"name":name[i], "radius":radius[i], "center":center}
    else:
        assert len(name) == 1
        i = 0
        name = name[i]
        bpy.ops.curve.simple(align='WORLD', location=location[i], rotation=[0,0,0], Simple_Type='Polygon',shape='3D',Simple_sides=sides[i],Simple_radius=radius[i], outputType='BEZIER', use_cyclic_u=True, edit_mode=False)
        curve = bpy.context.object
        curve.rotation_mode = 'QUATERNION'
        curve.rotation_quaternion = rotation[i]
        curve.name = name 
        curveData = bpy.data.objects[name].data

        curveData.resolution_u = 24
        curveData.bevel_mode = "OBJECT"
        if fill_caps == "both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.bevel_object = bpy.data.objects[profile_name]

        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)
        solidify(name,thickness)

"""
Perform a Boolean operation on both objects and delete the second object
"""
def boolean_operation(name1, name2, operation="UNION", flip_normals=False, solver_mode="EXACT"):
    bpy.data.objects[name1].modifiers.new("Boolean", "BOOLEAN")
    bpy.data.objects[name1].modifiers["Boolean"].object = bpy.data.objects[name2]
    bpy.data.objects[name1].modifiers["Boolean"].operation = operation
    bpy.data.objects[name1].modifiers["Boolean"].solver = solver_mode 
    bpy.context.view_layer.objects.active = bpy.data.objects[name1]
    bpy.ops.object.modifier_apply(modifier="Boolean")
    bpy.data.objects.remove(bpy.data.objects[name2], do_unlink=True)
    recalculate_normals(name1, inside=flip_normals)
    

"""
The arc turns to the Bézier control point, and the arc angle is an acute angle
"""
def arc2bezier(start_angle,end_angle,r,cx,cy):
    start_angle = start_angle*pi/180
    end_angle = end_angle*pi/180
    x0 = cx + np.cos(start_angle)*r
    y0 = cy + np.sin(start_angle)*r
    x3 = cx + np.cos(end_angle)*r
    y3 = cy + np.sin(end_angle)*r
    arc_angle = end_angle-start_angle
    a = 4*np.tan(arc_angle/4)/3
    x1 = x0 - a*(y0 - cy)
    y1 = y0 + a*(x0 - cx)
    x2 = x3 + a*(y3 - cy)
    y2 = y3 - a*(x3 - cx)
    
    return [x0,y0,0],[x1,y1,0],[x2,y2,0],[x3,y3,0]
    
"""
Return to the Bézier control point at any angle
"""
def get_bezier_points_by_angles(start_angle,end_angle,r,cx,cy):
    
    if start_angle > end_angle:
        tmp = start_angle
        start_angle = end_angle
        end_angle = tmp
    add_angle = end_angle - start_angle
    points_list = []
    if add_angle>90:
        times = int(np.ceil(add_angle/90))
        angle_delta = add_angle/times
        for i in range(times):
            if i != times-1:
                s_an = start_angle + i*angle_delta
                e_an = start_angle + (i+1)*angle_delta
                points0,points1,points2,points3 = arc2bezier(s_an,e_an,r,cx,cy)
                
                if i==0:
                    points_list+=[points0,points1,points2,points3]
                else:
                    points_list+=[points1,points2,points3]
            else:
                s_an = start_angle + i*angle_delta
                e_an = end_angle
                points0,points1,points2,points3 = arc2bezier(s_an,e_an,r,cx,cy)
                points_list+=[points1,points2,points3]
    else:
        points0,points1,points2,points3 = arc2bezier(start_angle,end_angle,r,cx,cy)
        points_list+=[points0,points1,points2,points3]
                    
    
    return points_list   


"""
Draw a Bezier arc based on the three 3D coordinates on the arc
"""
def create_arc_by_3Dpoints(name,profile_name=None,control_points=[],center="POINT",points_radius=[1.0,1.0],closed=False,thickness=None,fill_caps='none'):
    
    points = copy.deepcopy(control_points)
    p1, p2, p3 = np.array(points[0]), np.array(points[1]), np.array(points[2])
    n = np.cross((p1-p2),(p1-p3))
    x_axis = p1-p2
    y_axis = np.cross((p1-p2),n)
    z_axis = np.cross(x_axis,y_axis)
    if np.linalg.norm(x_axis)<1e-4:
        x_axis*=1e5
    if np.linalg.norm(y_axis)<1e-4:
        y_axis*=1e5
    if np.linalg.norm(z_axis)<1e-4:
        z_axis*=1e5
    x_axis/=np.linalg.norm(x_axis)
    y_axis/=np.linalg.norm(y_axis)
    z_axis/=np.linalg.norm(z_axis)
    
    point3D_to_xyplane=np.array([x_axis,y_axis,z_axis])
    try:
        xyplane_to_point3D = np.linalg.inv(point3D_to_xyplane)
    except:
        raise ValueError("Three-point collinear")

    p1_xy = point3D_to_xyplane@np.array(p1)
    p2_xy = point3D_to_xyplane@np.array(p2)
    p3_xy = point3D_to_xyplane@np.array(p3)
    
    points = np.array([p1_xy,p2_xy,p3_xy])
    points_z = points[0][2]
    points_2D = points[:,:2]
    
    cx,cy,r = find_circle(points_2D)
        
    angle_start = calculate_angle(cx, cy, points[0][0], points[0][1])
    angle_end = calculate_angle(cx, cy, points[2][0], points[2][1])

    
    # Calculate angles in degrees
    angle_start = math.degrees(angle_start) if angle_start>0 else math.degrees(angle_start) + 360
    angle_end = math.degrees(angle_end) if angle_end>0 else math.degrees(angle_end) + 360

    
    control_points = np.array(get_bezier_points_by_angles(angle_start,angle_end,r,cx,cy))
    control_points[:,2] = points_z 
    control_points = xyplane_to_point3D@np.array(control_points).T
    control_points = control_points.T
    if np.isnan(control_points).any() or np.isinf(control_points).any():
        raise ValueError("points can't form an arc")
    
    curveData = bpy.data.curves.new(name, type='CURVE')
    curveData.dimensions = '3D'

    bezierSpline = curveData.splines.new('BEZIER')

    num_cycle = int(np.ceil(len(control_points)/3))
    bezierSpline.bezier_points.add(num_cycle-1)
    
    # Generate all points' radius
    points_radius= [points_radius[0] + i * (points_radius[1] - points_radius[0]) / (num_cycle - 1) for i in range(num_cycle)]

    
    for i in range(num_cycle):
        bezierSpline.bezier_points[i].co = control_points[i*3]
        bezierSpline.bezier_points[i].radius = points_radius[i]
        # The last cycle has only one point
        if i<num_cycle-1:
            bezierSpline.bezier_points[i].handle_right = control_points[i*3+1]
            bezierSpline.bezier_points[i+1].handle_left = control_points[i*3+2]

    bezierSpline.use_cyclic_u = closed
    curveOB = bpy.data.objects.new(name, curveData)
    
    if profile_name==None:

        scn = bpy.context.scene.collection
        scn.objects.link(curveOB)
        if center=="MEDIAN":
            bpy.context.view_layer.objects.active = bpy.data.objects[name]
            bpy.ops.object.mode_set(mode = 'OBJECT')
            bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY', center='MEDIAN')

    else:
        curveData.bevel_mode = "OBJECT"
        if fill_caps=="both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.splines[0].use_smooth = False
        curveData.bevel_object = bpy.data.objects[profile_name]

        scn = bpy.context.scene.collection
        scn.objects.link(curveOB)

        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)
            
        if thickness>1e-10:
            solidify(name,thickness)
        

        return curveOB
"""
Creates an arc and returns three control points
"""
def create_arc_points(name, radius, start_angle, end_angle, rotation_rad=0,location=[0,0,0],rotation=[0,0,0],center="POINT", closed=False):
    # start_angle = start_angle*180
    # end_angle = end_angle*180
    bpy.ops.curve.simple(align='WORLD', location=location, rotation=rotation, Simple_Type='Arc',Simple_sides=3,Simple_radius=radius,Simple_startangle=start_angle, Simple_endangle=end_angle, outputType='BEZIER', use_cyclic_u=closed)
    curve = bpy.context.object
    curve.name = name
    curveData = curve.data
    curveData.dimensions = '3D'
    
    bpy.ops.object.mode_set(mode = 'OBJECT')
    curve.rotation_euler = [0,0,rotation_rad*pi]
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)

    mat = np.array(bpy.data.objects[name].matrix_world)
    front = mat @ np.array(list(curveData.splines[0].bezier_points[0].co) + [1])
    middle = mat @ np.array(list(curveData.splines[0].bezier_points[2].co) + [1])
    end = mat @ np.array(list(curveData.splines[0].bezier_points[4].co) + [1])
    
    points = np.array([front[:3],middle[:3],end[:3]])
    
    return {"name":name, "points":points - points[0],"center":center}



def add_simple_deform_modifier(name, mode='BEND', angle=0., origin=[0, 0, 0], rotation=Quaternion([1, 0, 0, 0]), axis='Z'):
    """
    Add a simple deform modifier to the specified object and apply it, then delete the Empty object used as the deformation origin.

    :param obj_name: The name of the object to add the modifier to
    :param mode: The deformation mode, default is 'BEND' (other options include 'TAPER', 'TWIST', 'SQUEEZE')
    :param angle: The bending angle (in radians, normalized to -1~1)
    :param origin: The coordinates of the bending starting point (requires an empty object as the starting point)
    :param axis: The bending axis, default is 'Z' (options include 'X', 'Y', 'Z')
    """
    origin = origin.copy()
    rotation = rotation.copy()
    if isinstance(rotation[-1], list):
        rotation = Quaternion(rotation[-1])
    else:
        rotation = Quaternion(rotation)
    if isinstance(origin[-1], list):
        origin = origin[-1]
    # Get the object
    obj = bpy.data.objects.get(name)
    if obj is None:
        raise ValueError(f"未找到名为 '{name}' 的对象")

    # Ensure in object mode
    if bpy.context.object and bpy.context.object.mode != 'OBJECT':
        bpy.ops.object.mode_set(mode='OBJECT')

    # Set the target object as the active object and select it
    bpy.context.view_layer.objects.active = obj
    for o in bpy.context.selected_objects:
        o.select_set(False)
    obj.select_set(True)

    # Create a simple deform modifier
    modifier = obj.modifiers.new(name="SimpleDeform", type='SIMPLE_DEFORM')
    modifier.deform_method = mode
    modifier.deform_axis = axis
    modifier.angle = angle * pi

    # Create an empty object as the starting point of deformation
    bpy.ops.object.empty_add(type='PLAIN_AXES', location=origin)
    empty = bpy.context.active_object
    empty.name = f"{name}_DeformOrigin"
    empty.rotation_mode = 'QUATERNION'
    empty.rotation_quaternion = rotation

    # Set the starting point of the modifier
    modifier.origin = empty

    # Ensure the active object is the target object and select it
    bpy.context.view_layer.objects.active = obj
    for o in bpy.context.selected_objects:
        o.select_set(False)
    obj.select_set(True)

    # Apply the modifier
    try:
        bpy.ops.object.modifier_apply(modifier=modifier.name)
    except RuntimeError as e:
        raise RuntimeError(f"应用修改器时出错: {e}")

    # Delete the Empty object
    # Since we named the empty when creating it, we can directly access it here
    for o in bpy.context.selected_objects:
        o.select_set(False)
    empty.select_set(True)
    bpy.context.view_layer.objects.active = empty
    bpy.ops.object.delete()


def create_trapezoid_body(name="trapezoid_body", size=[1, 1, 0.5, 0.5, 1], location=None, rotation=None, rotation_mode='QUATERNION', apply=False):
    """
    Create a trapezoid body, directly operate on Mesh data
    for size:
    :param height: The height of the trapezoid body
    :param top_length: The length of the top of the trapezoid body
    :param top_width: The width of the top of the trapezoid body
    :param bottom_length: The length of the bottom of the trapezoid body
    :param bottom_width: The width of the bottom of the trapezoid body
    """
    # Create a new mesh and object
    mesh = bpy.data.meshes.new("TrapezoidMesh")
    obj = bpy.data.objects.new(name, mesh)
    bpy.context.collection.objects.link(obj)

    top_length, top_width, bottom_length, bottom_width, height = size

    # Define vertices
    vertices = [
        # Top vertices
        (height / 2, top_length / 2, top_width / 2),
        (height / 2, -top_length / 2, top_width / 2),
        (height / 2, -top_length / 2, -top_width / 2),
        (height / 2, top_length / 2, -top_width / 2),
        # Bottom vertices
        (- height / 2,  bottom_length / 2, - bottom_width / 2),
        (- height / 2, - bottom_length / 2, - bottom_width / 2),
        (- height / 2, - bottom_length / 2, bottom_width / 2),
        (- height / 2, bottom_length / 2, bottom_width / 2)
    ]

    # Define faces (through vertex indices)
    faces = [
        # Top face
        (0, 1, 2, 3),
        # Bottom face
        (7, 6, 5, 4),
        # Four side faces
        (0, 3, 4, 7),
        (1, 0, 7, 6),
        (2, 1, 6, 5),
        (3, 2, 5, 4)
    ]

    # Create mesh data
    mesh.from_pydata(vertices, [], faces)
    mesh.update()

    # Switch to object mode
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='OBJECT')

    if location:
        obj.location = location
    if rotation:
        if rotation_mode=='XYZ':
            obj.rotation_euler = [angle * pi for angle in rotation]
        elif rotation_mode=='QUATERNION':
            obj.rotation_mode = 'QUATERNION'
            obj.rotation_quaternion = rotation
        elif rotation_mode=='MATRIX':
            mat = np.eye(4)
            rotation = np.array(rotation).reshape([3,3])
            mat[:3,:3] = rotation
            bpy.context.view_layer.update()
            world_matrix = torch.tensor(bpy.data.objects[name].matrix_world)
            scale_now = world_matrix.norm(dim=0)[:3]
            scale_matrix = torch.eye(4)
            scale_matrix[0,0],scale_matrix[1,1],scale_matrix[2,2] =scale_now[0],scale_now[1],scale_now[2]
            scale_matrix_inv = scale_matrix.clone()
            for i in range(3):
                # Invert the diagonal matrix, if the element on the diagonal is 0, it remains 0 
                if scale_matrix_inv[i,i]>1e-10:
                    scale_matrix_inv[i,i]=1.0 / scale_matrix_inv[i,i]
            mat = scale_matrix_inv@torch.tensor(mat,dtype=torch.float32)@scale_matrix
            mat = mathutils.Matrix(np.array(mat))
            bpy.data.objects[name].matrix_world = bpy.data.objects[name].matrix_world@mat
    
    if apply:
        for o in bpy.context.selected_objects:
            o.select_set(False)
        obj.select_set(True)
        bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)

    return obj

"""
create an circle and return radius
"""
def create_circle(name: list[str] | str, location: list[list[float]] | list[float]=[0., 0., 0.], profile_name=None, rotation: list[list[float]]=[1., 0., 0., 0.], radius: list[float] | float=1., scale: list[list[float]] | list[float]=[1., 1., 1.], use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH, resolution: list[int] | int=32, bend_angle: list[float] | float=0., center="MEDIAN", thickness=None, fill_caps="none"):
    """
    Create a circle object or multiple circle objects.
    Args:
        name (list[str] | str): Name or list of names for circle objects.
        location (list[list[float]] | list[float]): Location or list of locations for circle objects. 
        profile_name (str, optional): Name of profile object. If None, generates each circle separately.
        rotation (list[list[float]], optional): List of rotations for circle objects. Defaults to [1., 0., 0., 0.].
        radius (list[float] | float, optional): List of radii or single radius value. Defaults to 1..
        scale (list[list[float]] | list[float], optional): List of scales or single scale value. Defaults to [1., 1., 1.].
        bend_angle (list[float] | float, optional): List of bend angles or single bend angle. Defaults to 0.. Values normalized to -1~1.
        center (str, optional): Center type for circle objects. Defaults to "MEDIAN".
        thickness (float, optional): Thickness of circle objects.
        fill_caps (str, optional): End cap fill type for circle objects. Defaults to "none".
        resolution (list[int] | int): Resolution of the circle. Higher values create smoother circles.
    """
    location = location.copy()
    rotation = rotation.copy()
    scale = scale.copy()
    if not isinstance(name, list):
        name = [name]
    if not isinstance(location[0], list):
        location = [location]
    # The length of name and location must be the same, and each circle must have different name and location
    assert len(name) == len(location)

    # rotation, radius, scale, bend_angle can be the same value
    if not isinstance(rotation[0], list):
        rotation = [rotation] * len(name)
    if not isinstance(radius, list):
        radius = [radius] * len(name)
    if not isinstance(scale[0], list):
        scale = [scale] * len(name)
    if not isinstance(bend_angle, list):
        bend_angle = [bend_angle] * len(name)
    
    # Determine the resolution
    if use_minimum_face:
        resolution = 32
    elif not average_edge_length is None:
        resolution = (np.array(radius) * 2 * pi / average_edge_length).tolist()
    
    if not isinstance(resolution, list):
        resolution = [resolution] * len(name)

    if profile_name is None:  
        # Generate each circle separately
        for i in range(len(name)):
            bpy.ops.curve.simple(align='WORLD', location=location[i], rotation=[0,0,0], Simple_Type='Circle',shape='3D',Simple_sides=4,Simple_radius=radius[i], outputType='BEZIER', use_cyclic_u=True, edit_mode=False)
            curve = bpy.context.object
            curve.name = name[i]   
            curve.rotation_mode = 'QUATERNION'
            curve.rotation_quaternion = rotation[i]
            if center=="POINT":
                curve.location = [radius[i],0,0]
                bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
            curveData = curve.data
            curveData.dimensions = '3D'
            bpy.data.objects[name[i]].scale = scale[i]
            curveData.resolution_u = (resolution[i] + 2) // 4
            # If you need to bend, convert to a mesh object
            if bend_angle[i] != 0.:
                bpy.ops.object.convert(target='MESH')
                add_simple_deform_modifier(name[i], mode='BEND', origin=location[i], angle=bend_angle[i])
        
        return {"name":name[i], "radius":radius[i], "center":center}
    else:
        assert len(name) == 1
        i = 0
        name = name[i]
        bpy.ops.curve.simple(align='WORLD', location=location[i], rotation=[0,0,0], Simple_Type='Circle',shape='3D',Simple_sides=4,Simple_radius=radius[i], outputType='BEZIER', use_cyclic_u=True, edit_mode=False)
        curve = bpy.context.object
        curve.rotation_mode = 'QUATERNION'
        curve.rotation_quaternion = rotation[i]
        curve.name = name 
        curveData = bpy.data.objects[name].data
        curveData.resolution_u = 24
        curveData.bevel_mode = "OBJECT"
        if fill_caps == "both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.bevel_object = bpy.data.objects[profile_name]

        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)
        solidify(name,thickness)

"""
create an oval and return param
"""
def create_oval(name, a,b,rotation_rad=0,center="MEDIAN"):
    bpy.ops.curve.simple(align='WORLD', location=[0,0,0], rotation=[0,0,0], Simple_Type='Circle',shape='3D',Simple_sides=4,Simple_radius=1.0, outputType='BEZIER', use_cyclic_u=True,edit_mode=False)
    curve = bpy.context.object
    curve.name = name
    if center=="POINT":
        curve.location = [1,0,0]
        bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
    curve.scale = [a,b,0]
    curve.rotation_euler=[0,0,rotation_rad*pi]
    bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)

    return {"name":name, "a":a,"b":b,"rotation_rad":rotation_rad,"center":center} 


"""
Creates a generic quad and returns four control points
"""
def create_quad(name, control_points,center="POINT", resolution=12):
    if isinstance(name, str):
        control_points = copy.deepcopy(control_points)
        curveData = bpy.data.curves.new(name, type='CURVE')
        curveData.dimensions = '3D'
        curveData.resolution_u = resolution

        bezierSpline = curveData.splines.new('BEZIER')
        bezierSpline.bezier_points.add(len(control_points) - 1) 
        bezierSpline.use_cyclic_u = True

        for i, coord in enumerate(control_points):
            x, y, z = coord
            bezier_point = bezierSpline.bezier_points[i]
            bezier_point.co = (x, y, z)
            bezier_point.handle_left_type = 'VECTOR'
            bezier_point.handle_right_type = 'VECTOR'

        curveOB = bpy.data.objects.new(name, curveData)

        scn = bpy.context.scene.collection
        scn.objects.link(curveOB)
        if center=="MEDIAN":
            bpy.data.objects[name].select_set(True)
            bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY', center='MEDIAN')
        points = np.array(control_points)
        
        return {"name":name,"points":points - points[0],"center":center}
    #return {"name":name,"points":points_sorted - points_sorted[0]}
    else:
        points = []
        for i in range(len(name)):
            points.append(create_quad(name[i], control_points[i],center=center,resolution=resolution)['points'])
        return {"name":name,"points":control_points,"center":center}



"""
create a square
"""
def create_square(name,profile_name=None,length=1.0,rotation_rad=0, location=None, rotation=None, closed=True, center="MEDIAN", thickness=None, fill_caps='none',use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH, resolution=24):
    bpy.ops.curve.simple(align='WORLD',location=(0,0,0), rotation=(0,0,0),Simple_Type='Rectangle', Simple_width=length,shape='3D', Simple_length=length, use_cyclic_u=closed,edit_mode=False)
    curve = bpy.context.object
    curve.name = name
    curveData = curve.data
    curveData.resolution_u = 12
    curveData.dimensions = '3D'

    if profile_name==None:
        if center=="POINT":
            curve.location=[length/2,length/2,0]
            bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
        curve.rotation_euler=[0,0,rotation_rad*pi]
        bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
        
        return {"name":name,"length":length,"rotation_rad":rotation_rad,"center":center}
    
    else:
        curveData.location = location
        curve.rotation_mode = "QUATERNION"
        curve.rotation_quaternion = rotation
        curveData.bevel_mode = "OBJECT"
        if fill_caps=="both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.splines[0].use_smooth = False
        curveData.bevel_object = bpy.data.objects[profile_name]
        
        #bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        #bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)
        solidify(name,thickness)

"""
create a rectangle and return param
"""
def create_rectangle(name, profile_name=None, width=1.0, length=1.0, rotation_rad=0, location=None, rotation=None, closed=True, bevel_width=None, bevel_segments=8, center="MEDIAN", thickness=None, fill_caps='none', use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH, resolution=24):
    bpy.ops.curve.simple(align='WORLD',location=[0,0,0], rotation=[0,0,0],Simple_Type='Rectangle', Simple_width=width,shape='3D', Simple_length=length, use_cyclic_u=closed,edit_mode=False)
    curve = bpy.context.object
    curve.name = name
    curveData = curve.data
    curveData.resolution_u = 12
    curveData.dimensions = '3D'

    if profile_name==None:
        if center=="POINT":
            curve.location=[width/2,length/2,0]
            bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
        curve.rotation_euler=[0,0,rotation_rad*pi]
        bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)    
        return {"name":name, "width":width,"length":length,"rotation_rad":rotation_rad,"center":center,"closed":closed}
    
    else:
        curve.location = location
        curve.rotation_mode = "QUATERNION"
        curve.rotation_quaternion = rotation
        curveData.bevel_mode = "OBJECT"
        if bevel_width!=None:
            for p in curveData.splines[0].bezier_points:
                p.handle_left_type="FREE"
                p.handle_right_type="FREE"
            bevel(name=name, width=bevel_width, segments=bevel_segments)
            curveData = bpy.data.objects[name].data
            curveData.bevel_mode = 'OBJECT'
        if fill_caps=="both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.splines[0].use_smooth = False
        curveData.bevel_object = bpy.data.objects[profile_name]
        
        #bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        #bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)
        solidify(name,thickness)
        
"""
create a spiral curve
"""
def create_spiral(name, profile_name=None,dif_z=0.5,radius=0.5,turns=5,closed=False, center="POINT",thickness=0, handle_type="AUTO",fill_caps="none",flip_normals=False, direction="COUNTER_CLOCKWISE", resolution=8 ,location=[0,0,0], rotation=[0,0,0,0]):
    bpy.ops.curve.spirals(align='WORLD', location=(0, 0, 0), rotation=(0, 0, 0), spiral_type='ARCH', turns=turns, dif_z=dif_z, radius=radius, shape='3D', curve_type='BEZIER', use_cyclic_u=closed,edit_mode=False, spiral_direction=direction, handleType=handle_type,steps=12)
    curve = bpy.context.object
    curve.name = name
    curveData = curve.data
    curveData.resolution_u = resolution
    #curveOB = bpy.data.objects.new(name, curveData)
    if thickness is None:
        thickness = 0
    if profile_name != None:
        curveData.bevel_mode = "OBJECT"
        if fill_caps=="both":
            curveData.use_fill_caps = True
        else:
            curveData.use_fill_caps = False
        curveData.splines[0].use_smooth = False
        curveData.bevel_object = bpy.data.objects[profile_name]

        #scn = bpy.context.scene.collection
        #scn.objects.link(curveOB)
        bpy.context.view_layer.objects.active = bpy.data.objects[name]
        bpy.ops.object.mode_set(mode = 'OBJECT')
        bpy.data.objects[name].select_set(True)
        bpy.ops.object.convert(target='MESH')
        bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
        bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        
        if flip_normals:
            bpy.data.objects[name].select_set(True)
            bpy.ops.object.mode_set(mode = 'EDIT')
            bpy.ops.mesh.select_all(action='SELECT')
            bpy.ops.mesh.flip_normals()
            bpy.ops.object.mode_set(mode = 'OBJECT')
    
        if fill_caps in ["start","end"]:
            make_caps(name,fill_caps)

        recalculate_normals(name,inside=False)
            
        if thickness>1e-10:
            solidify(name,thickness)

        weld(name,1e-5)
    
    else:
        if center=="MEDIAN":
            bpy.data.objects[name].select_set(True)
            bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY', center='MEDIAN')
    obj = bpy.data.objects[name]
    if rotation is not None:
        obj.rotation_mode = 'QUATERNION'
        #print(rotation)
        obj.rotation_quaternion = rotation
    if location is not None:
        obj.location = location
    bpy.ops.object.transform_apply(True)
    return {
        "name": name, "profile_name": profile_name,"dif_z": dif_z,"radius": radius, "turns": turns,"thickness": thickness, "location": location, "rotation": rotation, "fill_caps": fill_caps, "direction": direction
    }



"""
Create a rotation object
"""
def bezier_rotation(name,profile_name,location=[0,0,0], rotation=[0,0,0,0],thickness=0.,use_smooth=False,resolution=12):
    r=1
    bpy.ops.curve.simple(align='WORLD', location=[0,0,0], rotation=[0,0,0], Simple_Type='Circle',shape='3D',Simple_sides=4,Simple_radius=r, outputType='BEZIER', use_cyclic_u=True,edit_mode=False)
    curve = bpy.context.object
    curve.name = name
    curveData = bpy.data.objects[name].data
    curveData.resolution_u = resolution
    curveData.offset=r
    curve.location = location
    curve.rotation_mode = "QUATERNION"
    curve.rotation_quaternion = rotation
    
    curveData.bevel_mode = "OBJECT"
    curveData.splines[0].use_smooth = use_smooth
    curveData.bevel_object = bpy.data.objects[profile_name]
    
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.mode_set(mode = 'OBJECT')
    bpy.data.objects[name].select_set(True)
    bpy.ops.object.convert(target='MESH')
    bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
    bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
    solidify(name, thickness)

def curve2meshgrid(name, thickness, bevel_width, is_seat_subsurf):
    obj = bpy.data.objects[name]
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    curve2mesh(obj)
    bpy.ops.object.mode_set(mode="EDIT")
    bpy.ops.mesh.select_all(action="SELECT")
    bpy.ops.mesh.fill_grid(use_interp_simple=True)
    bpy.ops.object.mode_set(mode="OBJECT")
    
    solidify(obj.name, thickness)
    subsurf(obj.name, 1, not is_seat_subsurf)
    bevel(obj.name, bevel_width, segments=8)
    return obj

def fill_grid(name, thickness=None, use_interp_simple=True):
    obj = bpy.data.objects[name]
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    if obj.type != 'MESH':
        assert obj.type == 'CURVE', f"the type of obj is {obj.type} rather than \'CURVE\' or \'MESH\'"
        curve2mesh(obj)
    bpy.ops.object.mode_set(mode="EDIT")
    bpy.ops.mesh.select_all(action="SELECT")
    bpy.ops.mesh.fill_grid(use_interp_simple=use_interp_simple)
    bpy.ops.object.mode_set(mode="OBJECT")
    if thickness is not None:
        solidify(name, thickness)

def curve2mesh(obj, weld_mod=True, merge_threshold=1e-4):
    bpy.ops.object.mode_set(mode = 'OBJECT')
    obj.select_set(True)
    bpy.ops.object.convert(target="MESH")
    obj = bpy.context.active_object
    if weld_mod:
        weld(obj.name, merge_threshold=merge_threshold)
    return obj

"""
Apply WELD modifier to the object.
"""
def weld(name, merge_threshold=1e-3):
    bpy.data.objects[name].modifiers.new("Weld", "WELD")
    bpy.data.objects[name].modifiers["Weld"].merge_threshold = merge_threshold
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.modifier_apply(modifier="Weld")


def bevel(name, width, segments=8):
    """
    This function aims to add bevel modifier and apply to the object
    """
    if width==0:
        return
    obj = bpy.data.objects.get(name)
    obj.select_set(True)
    bpy.data.objects[name].modifiers.new("Bevel", "BEVEL")
    bpy.data.objects[name].modifiers["Bevel"].width = width
    bpy.data.objects[name].modifiers["Bevel"].segments = segments
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    if bpy.data.objects[name].type == 'MESH':
        bpy.ops.object.modifier_apply(modifier="Bevel")
    elif bpy.data.objects[name].type == 'CURVE':
        bpy.context.object.modifiers["Bevel"].affect = 'VERTICES'
        bpy.ops.object.convert(target='MESH')
        bpy.ops.object.modifier_apply(modifier="Bevel")
        bpy.ops.object.convert(target='CURVE')

def subsurf(name: list[str] | str, levels: list[int] | int, simple: list[bool] | bool=False):
    """ 
    This function aims to add subdivide surface modifier and apply to the object
    TODO: There are some properties that haven't added to this function yet.
    """
    if not isinstance(name, list):
        name = [name]
    if not isinstance(levels, list):
        levels = [levels] * len(name)
    if not isinstance(simple, list):
        simple = [simple] * len(name)
    for i in range(len(name)):
        if levels[i] > 0:
            if bpy.data.objects.get(name[i]).type == 'CURVE':
                bpy.context.view_layer.objects.active = bpy.data.objects[name[i]]
                bpy.data.objects[name[i]].select_set(True)
                bpy.ops.object.convert(target='MESH')
                bpy.data.objects[name[i]].select_set(False)
            bpy.data.objects[name[i]].modifiers.new("Subsurf", "SUBSURF")
            bpy.data.objects[name[i]].modifiers["Subsurf"].levels = levels[i]
            bpy.data.objects[name[i]].modifiers["Subsurf"].render_levels = levels[i]
            bpy.data.objects[name[i]].modifiers["Subsurf"].subdivision_type = "SIMPLE" if simple[i] else "CATMULL_CLARK"
            bpy.context.view_layer.objects.active = bpy.data.objects[name[i]]
            bpy.ops.object.modifier_apply(modifier="Subsurf")

"""
Bridge edge loops of more than 2 loops
"""
def bridge_edge_loops(name: str | list[str], profile_name: list[str] | list[list[str]], number_cuts: int | list[int]=None, smoothness: float | list[float]=0., profile_shape_factor: float | list[float]=0., twist: int | list[int]=0, interpolation: str | list[str]='PATH', profile_shape: str | list[str]='SMOOTH', flip_normals: list[bool] | bool=False, fill_caps: list[str] | str="none", use_minimum_face=USE_MINIMUM_FACE, average_edge_length=AVERAGE_EDGE_LENGTH):
    """
    Bridges edge loops between multiple objects in Blender.

    Args:
        name (str | list[str]): The name(s) to assign to the resulting object(s).
        profile_name (list[str] | list[list[str]]): A list of names of the objects to be bridged.
        number_cuts (int | list[int], optional): The number of cuts to make along the bridge. Defaults to 8.
        smoothness (float | list[float], optional): The smoothness of the bridge. Defaults to 0.
        profile_shape_factor (float | list[float], optional): The factor for the profile shape. Defaults to 0.
        twist (int | list[int], optional): The twist offset for the bridge. Defaults to 0.
        interpolation (str | list[str], optional): The interpolation method for the bridge. Must be one of ['LINEAR', 'PATH', 'SURFACE']. Defaults to 'PATH'.
        profile_shape (str | list[str], optional): The profile shape for the bridge. Must be one of ['SMOOTH', 'SPHERE', 'ROOT', 'INVERSE_SQUARE', 'SHARP', 'LINEAR']. Defaults to 'SMOOTH'.
        
    Raises:
        AssertionError: If fewer than two profile names are provided.
        AssertionError: If an invalid interpolation method is provided.
        AssertionError: If an invalid profile shape is provided.

    """
    if isinstance(name, str):
        assert len(profile_name) >= 2
        assert interpolation in ['LINEAR', 'PATH', 'SURFACE']
        assert profile_shape in ['SMOOTH', 'SPHERE', 'ROOT', 'INVERSE_SQUARE', 'SHARP', 'LINEAR']

        # If the profile is a curve, it needs to be converted to a mesh
        for b_name in profile_name:
            if bpy.data.objects.get(b_name).type == 'CURVE':
                bpy.context.view_layer.objects.active = bpy.data.objects[b_name]
                bpy.data.objects[b_name].select_set(True)
                bpy.ops.object.convert(target='MESH')
                bpy.data.objects[b_name].select_set(False)
        
        if fill_caps=='start':
            bpy.ops.object.select_all(action='DESELECT')
            bpy.data.objects[profile_name[0]].select_set(True)
            bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[0]]
            bpy.ops.object.mode_set(mode='EDIT')
            bpy.ops.mesh.select_all(action='SELECT')
            bpy.ops.mesh.edge_face_add()
            bpy.ops.object.mode_set(mode='OBJECT')
        elif fill_caps=='end':
            bpy.ops.object.select_all(action='DESELECT')
            bpy.data.objects[profile_name[-1]].select_set(True)
            bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[-1]]
            bpy.ops.object.mode_set(mode='EDIT')
            bpy.ops.mesh.select_all(action='SELECT')
            bpy.ops.mesh.edge_face_add()
            bpy.ops.object.mode_set(mode='OBJECT')

        scene = bpy.context.scene

        start_obj = bpy.data.objects[profile_name[0]].copy()
        end_obj = bpy.data.objects[profile_name[-1]].copy()
        start_obj.data = bpy.data.objects[profile_name[0]].data.copy()
        end_obj.data = bpy.data.objects[profile_name[-1]].data.copy()
        scene.collection.objects.link(start_obj)
        scene.collection.objects.link(end_obj)
        profile_name[0] += ".001"
        profile_name[-1] += ".001"

        bpy.ops.object.select_all(action='DESELECT')
        for b_name in profile_name:
            bpy.data.objects[b_name].select_set(True)
        if 0. not in bpy.data.objects[profile_name[0]].scale:
            bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[0]]
        else:
            bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[-1]]

        # mininal face
        if number_cuts == None:
            if use_minimum_face:
                number_cuts = 32
            elif not average_edge_length is None:
                avg_len = 0
                for i in range(len(profile_name) - 1):
                    avg_len += np.linalg.norm(bpy.data.objects[profile_name[i]].location - bpy.data.objects[profile_name[i + 1]].location)
                avg_len /= len(profile_name) - 1
                number_cuts = int(avg_len / average_edge_length)
                # clamp number_cuts to [32, 128]
                number_cuts = max(32, min(128, number_cuts))

        bpy.ops.object.join()
        bpy.context.object.name = name
        bpy.ops.object.mode_set(mode='EDIT')
        bpy.ops.mesh.select_all(action='SELECT')
        bpy.ops.mesh.bridge_edge_loops(
            twist_offset=twist,
            number_cuts=number_cuts,
            smoothness=smoothness,
            profile_shape_factor=profile_shape_factor,
            interpolation=interpolation,
            profile_shape=profile_shape
        )
        if fill_caps=='both':
            bpy.ops.mesh.select_all(action='SELECT')
            bpy.ops.mesh.region_to_loop()
            bpy.ops.mesh.edge_face_add()
            bpy.ops.mesh.select_all(action='DESELECT')
        bpy.ops.object.mode_set(mode='OBJECT')
        if flip_normals:
            recalculate_normals(name,inside=True)
        else:
            recalculate_normals(name,inside=False)
        bpy.data.objects[name].select_set(False)

    else:
        if isinstance(number_cuts, int):
            number_cuts = [number_cuts] * len(name)
        if isinstance(smoothness, float):
            smoothness = [smoothness] * len(name)
        if isinstance(profile_shape_factor, float):
            profile_shape_factor = [profile_shape_factor] * len(name)
        if isinstance(twist, int):
            twist = [twist] * len(name)
        if isinstance(interpolation, str):
            interpolation = [interpolation] * len(name)
        if isinstance(profile_shape, str):
            profile_shape = [profile_shape] * len(name)
        if isinstance(flip_normals, bool):
            flip_normals = [flip_normals] * len(name)
        if isinstance(fill_caps, str):
            fill_caps = [fill_caps] * len(name)

        for i in range(len(name)):
            assert len(profile_name[i]) >= 2
            assert interpolation[i] in ['LINEAR', 'PATH', 'SURFACE']
            assert profile_shape[i] in ['SMOOTH', 'SPHERE', 'ROOT', 'INVERSE_SQUARE', 'SHARP', 'LINEAR']

            # If the circle is a curve, it needs to be converted to a mesh
            for b_name in profile_name[i]:
                if bpy.data.objects.get(b_name).type == 'CURVE':
                    bpy.context.view_layer.objects.active = bpy.data.objects[b_name]
                    bpy.data.objects[b_name].select_set(True)
                    bpy.ops.object.convert(target='MESH')
                    bpy.data.objects[b_name].select_set(False)

            if fill_caps[i]=='start':
                bpy.ops.object.select_all(action='DESELECT')
                bpy.data.objects[profile_name[i][0]].select_set(True)
                bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[i][0]]
                bpy.ops.object.mode_set(mode='EDIT')
                bpy.ops.mesh.select_all(action='SELECT')
                bpy.ops.mesh.edge_face_add()
                bpy.ops.object.mode_set(mode='OBJECT')
            elif fill_caps[i]=='end':
                bpy.ops.object.select_all(action='DESELECT')
                bpy.data.objects[profile_name[i][-1]].select_set(True)
                bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[i][-1]]
                bpy.ops.object.mode_set(mode='EDIT')
                bpy.ops.mesh.select_all(action='SELECT')
                bpy.ops.mesh.edge_face_add()
                bpy.ops.object.mode_set(mode='OBJECT')

            scene = bpy.context.scene

            start_obj = bpy.data.objects[profile_name[i][0]].copy()
            end_obj = bpy.data.objects[profile_name[i][-1]].copy()
            start_obj.data = bpy.data.objects[profile_name[i][0]].data.copy()
            end_obj.data = bpy.data.objects[profile_name[i][-1]].data.copy()
            scene.collection.objects.link(start_obj)
            scene.collection.objects.link(end_obj)
            profile_name[i][0] += ".001"
            profile_name[i][-1] += ".001"

            bpy.ops.object.select_all(action='DESELECT')
            for b_name in profile_name[i]:
                bpy.data.objects[b_name].select_set(True)
            if 0. not in bpy.data.objects[profile_name[i][0]].scale:
                bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[i][0]]
            else:
                bpy.context.view_layer.objects.active = bpy.data.objects[profile_name[i][-1]]

            # mininal face
            if number_cuts[i] == None:
                if use_minimum_face:
                    number_cuts[i] = 32
                elif not average_edge_length is None:
                    avg_len = 0
                    for j in range(len(profile_name[i]) - 1):
                        avg_len += np.linalg.norm(bpy.data.objects[profile_name[i][j]].location - bpy.data.objects[profile_name[i][j + 1]].location)
                    avg_len /= len(profile_name[i]) - 1
                    number_cuts[i] = int(avg_len / average_edge_length)
                    # clamp number_cuts to [32, 128]
                    number_cuts[i] = max(32, min(128, number_cuts[i]))
                
            bpy.ops.object.join()
            bpy.context.object.name = name[i]
            bpy.ops.object.mode_set(mode='EDIT')
            bpy.ops.mesh.select_all(action='SELECT')
            bpy.ops.mesh.bridge_edge_loops(
                twist_offset=twist[i],
                number_cuts=number_cuts[i],
                smoothness=smoothness[i],
                profile_shape_factor=profile_shape_factor[i],
                interpolation=interpolation[i],
                profile_shape=profile_shape[i]
            )
            if fill_caps[i]=='both':
                bpy.ops.mesh.select_all(action='SELECT')
                bpy.ops.mesh.region_to_loop()
                bpy.ops.mesh.edge_face_add()
                bpy.ops.mesh.select_all(action='DESELECT')
            bpy.ops.object.mode_set(mode='OBJECT')
            if flip_normals[i]:
                recalculate_normals(name[i],inside=True)
            else:
                recalculate_normals(name[i],inside=False)
            bpy.data.objects[name[i]].select_set(False)

    objects = bpy.data.objects
    profile_name_list = list(itertools.chain(*profile_name)) if isinstance(profile_name[0], list) else profile_name
    profile_name_list = [profile_name.split('.')[0] for profile_name in profile_name_list]
    profile_name_set = set(profile_name_list)
    profile_to_delete = [obj for obj in objects if obj.name in profile_name_set]
    for obj in profile_to_delete:
        obj.select_set(True)
    bpy.ops.object.delete()

    return name


def add_curve_modifier_to_object(name, curve_name, origin=[0., 0., 0.], rotation=[1, 0, 0, 0], axis='POS_X'):
    """
    Add a curve modifier to the specified mesh object and link it to the specified curve object.
    
    :param name: Name of the mesh object
    :param curve_name: Name of the curve object
    :param axis: Deformation axis, 'POS_X', 'POS_Y', 'POS_Z', 'NEG_X', 'NEG_Y' or 'NEG_Z'
    """
    origin = origin.copy()
    rotation = rotation.copy()
    
    origin = origin[-1] if isinstance(origin[-1], list) else origin
    rotation = rotation[-1] if isinstance(rotation[-1], list) else rotation

    # Get the mesh object
    mesh_obj = bpy.data.objects.get(name)
    assert mesh_obj is not None and mesh_obj.type == 'MESH'
    
    # Get the curve object
    curve_obj = bpy.data.objects.get(curve_name)
    assert curve_obj is not None and curve_obj.type == 'CURVE'

    # If it is not a cardinal pose, first convert the mesh_obj to a cardinal pose and then transform it back
    bpy.context.view_layer.objects.active = curve_obj
    for o in bpy.context.selected_objects:
        o.select_set(False)
    curve_obj.select_set(True)
    bpy.ops.object.convert(target='MESH')
    if origin != [0, 0, 0]:
        curve_obj.location -= Vector(origin)
        bpy.ops.object.transform_apply(location=True, rotation=False, scale=False)
    if rotation != [1, 0, 0, 0]:
        # Rotate the curve_obj
        curve_obj.rotation_mode = 'XYZ'
        curve_obj.rotation_euler = (Quaternion(rotation).conjugated() @ curve_obj.rotation_quaternion).to_euler()
        bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
    bpy.ops.object.convert(target='CURVE')

    # Rotate the curve_obj back
    if rotation != [1, 0, 0, 0]:
        curve_obj.rotation_mode = 'QUATERNION'
        curve_obj.rotation_quaternion = Quaternion(rotation) @ curve_obj.rotation_quaternion
    # Move the mesh_obj back
    if origin != [0, 0, 0]:
        curve_obj.location += Vector(origin)

    # Add the curve modifier
    modifier = mesh_obj.modifiers.new(name="CurveModifier", type='CURVE')
    
    # Set the curve object
    modifier.object = curve_obj
    
    # Set the deformation axis
    assert axis.upper() in ['POS_X', 'POS_Y', 'POS_Z', 'NEG_X', 'NEG_Y', 'NEG_Z']
    modifier.deform_axis = axis.upper()

    # Let the mesh object be the parent of the curve object
    curve_obj.parent = mesh_obj
    curve_obj.matrix_parent_inverse = mesh_obj.matrix_world.inverted()
    
    # Ensure that the active object is the target object and select it
    bpy.context.view_layer.objects.active = mesh_obj
    for o in bpy.context.selected_objects:
        o.select_set(False)
    mesh_obj.select_set(True)

    # Apply the modifier
    try:
        bpy.ops.object.modifier_apply(modifier=modifier.name)
    except RuntimeError as e:
        raise RuntimeError(f"应用修改器时出错: {e}")

    # Delete the curve object
    for o in bpy.context.selected_objects:
        o.select_set(False)
    curve_obj.select_set(True)
    bpy.context.view_layer.objects.active = curve_obj
    bpy.ops.object.delete()


def delete_obj(sub_name: str):
    objects = bpy.data.objects
    circles_to_delete = [obj for obj in objects if sub_name in obj.name]
    for obj in circles_to_delete:
        obj.select_set(True)
    bpy.ops.object.delete()


def join_obj(name: str, seq_name: list[str], weld_threshold=None):
    for n in seq_name:
        bpy.data.objects[n].select_set(True)
    bpy.context.view_layer.objects.active = bpy.data.objects[seq_name[-1]]
    bpy.ops.object.join()
    bpy.data.objects[seq_name[-1]].name = name
    if weld_threshold:
        weld(name=name, merge_threshold=weld_threshold)
    return name

def create_star(name, outer_radius, inner_radius, num_points, edge="smooth"):
    # Create a simple circle curve using the built-in operator
    bpy.ops.curve.simple(
        align='WORLD',
        location=[0, 0, 0],
        rotation=[0, 0, 0],
        Simple_Type='Circle',
        shape='3D',
        Simple_sides=num_points * 2,  # Each point alternates between outer and inner radius
        Simple_radius=outer_radius,
        outputType='BEZIER',
        use_cyclic_u=True,
        edit_mode=False
    )
    
    curve = bpy.context.object
    curve.name = name
    # if center == "POINT":
    #     curve.location = [outer_radius, 0, 0]
    #     bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
    
    curve_data = curve.data
    curve_data.dimensions = '3D'
    
    # Adjust the control points to create a star shape
    angle_step = math.pi / num_points  # Angle step between each point
    spline = curve_data.splines[0]
    points = spline.bezier_points
    
    handle_type = "VECTOR" if edge == "sharp" else "AUTO"
    
    for i, point in enumerate(points):
        angle = i * angle_step
        # Alternate between outer and inner radius
        radius = outer_radius if i % 2 == 0 else inner_radius
        x = radius * math.cos(angle)
        y = radius * math.sin(angle)
        z = 0
        point.co = (x, y, z)
        
        # Set handles to smooth transition VECTOR
        point.handle_left_type = handle_type
        point.handle_right_type = handle_type
    
    return {"name": name, "outer_radius": outer_radius, "inner_radius": inner_radius, "num_points": num_points}

def star_curve_rotation(name,profile_name,num_points=16, r_diff=0.1, location=[0,0,0], rotation=[0,0,0,0],thickness=0.002,use_smooth=False):
    create_star(name, 1 + r_diff, 1, num_points, edge="sharp")
    curve = bpy.context.object
    curveData = bpy.data.objects[name].data
    curve.location = location
    curve.rotation_mode = "QUATERNION"
    curve.rotation_quaternion = rotation
    curveData.offset=(1 + r_diff)*1.05
    curveData.bevel_mode = "OBJECT"
    curveData.splines[0].use_smooth = use_smooth
    curveData.bevel_object = bpy.data.objects[profile_name]
    
    bpy.context.view_layer.objects.active = bpy.data.objects[name]
    bpy.ops.object.mode_set(mode = 'OBJECT')
    bpy.data.objects[name].select_set(True)
    bpy.ops.object.convert(target='MESH')
    bpy.data.objects.remove(bpy.data.objects[profile_name], do_unlink=True)
    bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
    solidify(name, thickness)
   
def create_fork(name, x_tip, thickness, n_cuts, x_anchors, y_anchors, z_anchors, location=(0.0, 0.0, 0.0), rotation=(0.0, 0.0, 0.0,0.0), scale=(1, 1, 1)):
    x_anchors = np.array(x_anchors)
    y_anchors = np.array(y_anchors)
    z_anchors = np.array(z_anchors)
    has_cut =True
    n = int(2 * (n_cuts + 1))
    obj = create_primitive(name=name, primitive_type="grid", x_subdivisions=len(x_anchors) - 1, y_subdivisions=n - 1)
    x = np.concatenate([x_anchors] * n)
    y = np.ravel(y_anchors[np.newaxis, :] * np.linspace(1, -1, n)[:, np.newaxis])
    z = np.concatenate([z_anchors] * n)
    arr = np.stack([x, y, z], -1)
    obj.data.vertices.foreach_set("co", arr.reshape(-1))

    #simply copy make cuts
    if has_cut:
        bpy.context.view_layer.objects.active = obj
        bpy.ops.object.mode_set(mode='EDIT')
        bm = bmesh.from_edit_mesh(obj.data)
        front_verts = []
        for v in bm.verts:
            if abs(v.co[0] - x_tip) < 1e-3:
                front_verts.append(v)
        front_verts = sorted(front_verts, key=lambda v: v.co[1])
        geom = []
        for f in bm.faces:
            vs = list(v for v in f.verts if v in front_verts)
            if len(vs) == 2:
                if min(front_verts.index(vs[0]), front_verts.index(vs[1])) % 2 == 1:
                    geom.append(f)
        bmesh.ops.delete(bm, geom=geom, context="FACES")
        bmesh.update_edit_mesh(obj.data)
        bpy.ops.object.mode_set(mode='OBJECT')
    solidify(name, thickness)
    subsurf(name, 1)
    subsurf(name, 1)
    arr = np.array([v.co for v in obj.data.vertices])
    center = np.array(((arr[:, 0].max() + arr[:, 0].min()) / 2, (arr[:, 1].max() + arr[:, 1].min()) / 2, (arr[:, 2].max() + arr[:, 2].min()) / 2))
    obj.location = -center
    bpy.ops.object.transform_apply(True)
    if scale is not None:
        obj.scale = scale
        bpy.ops.object.transform_apply(True)
    if rotation is not None:
        obj.rotation_mode = 'QUATERNION'
        #print(rotation)
        obj.rotation_quaternion = rotation
    bpy.ops.object.transform_apply(True)
    if location is not None:
        obj.location = location
    bpy.ops.object.transform_apply(True)
    return {"name": name, "x_tip": x_tip, "thickness": thickness,"n_cuts": n_cuts, "x_anchors": x_anchors, "y_anchors": y_anchors, "z_anchors": z_anchors, "location": location, "rotation": rotation, "scale": scale}

def create_spoon(name, z_depth, thickness, x_anchors, y_anchors, z_anchors, location=(0.0, 0.0, 0.0), rotation=(0.0, 0.0, 0.0,0.0), scale=(1, 1, 1)):
    obj = create_primitive(name=name, primitive_type="grid", x_subdivisions=len(x_anchors) - 1, y_subdivisions=2)
    x_anchors = np.array(x_anchors)
    y_anchors = np.array(y_anchors)
    z_anchors = np.array(z_anchors)
    x = np.concatenate([x_anchors] * 3)
    y = np.concatenate([y_anchors, np.zeros_like(y_anchors), -y_anchors])
    z = np.concatenate([z_anchors] * 3)
    x[len(x_anchors)] += 0.02
    z[len(x_anchors) + 1] = -z_depth
    arr = np.stack([x, y, z], -1)
    obj.data.vertices.foreach_set("co", arr.reshape(-1))
    solidify(name, thickness)
    subsurf(name, 1)
    subsurf(name, 2)
    arr = np.array([v.co for v in obj.data.vertices])
    center = np.array(((arr[:, 0].max() + arr[:, 0].min()) / 2, (arr[:, 1].max() + arr[:, 1].min()) / 2, (arr[:, 2].max() + arr[:, 2].min()) / 2))
    obj.location = -center
    bpy.ops.object.transform_apply(True)
    if scale is not None:
        obj.scale = scale
        bpy.ops.object.transform_apply(True)
    if rotation is not None:
        obj.rotation_mode = 'QUATERNION'
        #print(rotation)
        obj.rotation_quaternion = rotation
    bpy.ops.object.transform_apply(True)
    if location is not None:
        obj.location = location
    bpy.ops.object.transform_apply(True)
    return {"name": name, "z_depth": z_depth, "thickness": thickness,"x_anchors": x_anchors, "y_anchors": y_anchors, "z_anchors": z_anchors, "location": location, "rotation": rotation, "scale": scale}

def corrective_smooth(name, factor=0.500, repeat=5, scale=1.):
    if not isinstance(name, list):
        name = [name]
    if not isinstance(factor, list):
        factor = [factor] * len(name)
    if not isinstance(repeat, list):
        repeat = [repeat] * len(name)
    if not isinstance(scale, list):
        scale = [scale] * len(name)
    for i in range(len(name)):
        if bpy.data.objects.get(name[i]).type == 'CURVE':
            bpy.context.view_layer.objects.active = bpy.data.objects[name[i]]
            bpy.data.objects[name[i]].select_set(True)
            bpy.ops.object.convert(target='MESH')
            bpy.data.objects[name[i]].select_set(False)
        bpy.data.objects[name[i]].modifiers.new("CorrectiveSmooth", "CORRECTIVE_SMOOTH")
        bpy.data.objects[name[i]].modifiers["CorrectiveSmooth"].factor = factor[i]
        bpy.data.objects[name[i]].modifiers["CorrectiveSmooth"].iterations = repeat[i]
        bpy.data.objects[name[i]].modifiers["CorrectiveSmooth"].scale = scale[i]
        bpy.context.view_layer.objects.active = bpy.data.objects[name[i]]
        bpy.ops.object.modifier_apply(modifier="CorrectiveSmooth")

def array(obj, fit_type, count=None, fit_length=None, fit_curve=None, relative_offset=None, constant_offset=None, offset_object=None, merge_threshold=None, merge_cap=False, start_cap=None, end_cap=None, apply=True):
    modifier = obj.modifiers.new(name="Array", type='ARRAY')
    # Set the parameters of the Array Modifier
    if fit_type == "FIXED_COUNT":
        modifier.fit_type = 'FIXED_COUNT'
        modifier.count = count
    elif fit_type == "FIT_LENGTH":
        modifier.fit_type = 'FIT_LENGTH'
        modifier.fit_length = fit_length
    elif fit_type == "FIT_CURVE":
        modifier.fit_type = 'FIT_CURVE'
        modifier.curve = fit_curve

    if relative_offset:
        modifier.use_relative_offset = True
        modifier.relative_offset_displace = relative_offset
    else:
        modifier.use_relative_offset = False
    
    if constant_offset:
        modifier.use_constant_offset = True
        modifier.constant_offset_displace = constant_offset
    
    if offset_object:
        modifier.use_object_offset = True
        modifier.offset_object = offset_object

    if merge_threshold:
        bpy.context.object.modifiers["Array"].use_merge_vertices = True
        modifier.merge_threshold = merge_threshold
        if merge_cap:
            modifier.use_merge_vertices_cap = merge_cap

    if start_cap:
        modifier.start_cap = start_cap
    if end_cap:
        modifier.end_cap = end_cap

    if apply:
        bpy.context.view_layer.objects.active = obj
        bpy.ops.object.modifier_apply(modifier="Array")

def array_1d(name, fit_type, count=None, fit_length=None, fit_curve=None, relative_offset=None, constant_offset=None, offset_object=None, merge_threshold=None, merge_cap=False, start_cap=None, end_cap=None, apply=True):
    obj = bpy.data.objects[name]
    array(obj, fit_type, count, fit_length, fit_curve, relative_offset, constant_offset, offset_object, merge_threshold, merge_cap, start_cap, end_cap, apply)

def array_2d(name, fit_type_1, fit_type_2, count_1=None, fit_length_1=None, fit_curve_1=None, relative_offset_1=None, constant_offset_1=None, offset_object_1=None, merge_threshold_1=None, merge_cap_1=False, count_2=None, fit_length_2=None, fit_curve_2=None, relative_offset_2=None, constant_offset_2=None, offset_object_2=None, merge_threshold_2=None, merge_cap_2=False, start_cap=None, end_cap=None, apply=True):
    obj = bpy.data.objects[name]
    array(obj, fit_type_1, count_1, fit_length_1, fit_curve_1, relative_offset_1, constant_offset_1, offset_object_1, merge_threshold_1, merge_cap_1, start_cap, end_cap, apply)
    array(obj, fit_type_2, count_2, fit_length_2, fit_curve_2, relative_offset_2, constant_offset_2, offset_object_2, merge_threshold_2, merge_cap_2, start_cap, end_cap, apply)

def array_3d(name, fit_type_1, fit_type_2, fit_type_3, count_1=None, fit_length_1=None, fit_curve_1=None, relative_offset_1=None, constant_offset_1=None, offset_object_1=None, merge_threshold_1=None, merge_cap_1=False, count_2=None, fit_length_2=None, fit_curve_2=None, relative_offset_2=None, constant_offset_2=None, offset_object_2=None, merge_threshold_2=None, merge_cap_2=False, count_3=None, fit_length_3=None, fit_curve_3=None, relative_offset_3=None, constant_offset_3=None, offset_object_3=None, merge_threshold_3=None, merge_cap_3=False, start_cap=None, end_cap=None, apply=True):
    obj = bpy.data.objects[name]
    array(obj, fit_type_1, count_1, fit_length_1, fit_curve_1, relative_offset_1, constant_offset_1, offset_object_1, merge_threshold_1, merge_cap_1, start_cap, end_cap, apply)
    array(obj, fit_type_2, count_2, fit_length_2, fit_curve_2, relative_offset_2, constant_offset_2, offset_object_2, merge_threshold_2, merge_cap_2, start_cap, end_cap, apply)
    array(obj, fit_type_3, count_3, fit_length_3, fit_curve_3, relative_offset_3, constant_offset_3, offset_object_3, merge_threshold_3, merge_cap_3, start_cap, end_cap, apply)

def array_rotate(name, count=None, rotation_angle=0, rotation_origin=(0, 0, 0), z_offset=0, rotation_axis=None, rotation_mode='XYZ', scale=(1, 1, 1), fit_type="FIXED_COUNT", merge_threshold=None, merge_cap=False, start_cap=None, end_cap=None, apply=True, delete_empty=True):
    obj = bpy.data.objects[name]
    # Recover the original rotation, apply and then rotate as the same as rotation_axis
    for o in bpy.context.selected_objects:
        o.select_set(False)
    obj.select_set(True)
    bpy.ops.object.transform_apply(location=False, rotation=False, scale=True)
    obj.select_set(False)

    if rotation_axis is None or rotation_axis == obj.rotation_quaternion:
        rotation_axis = obj.rotation_quaternion
    else:
        bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY', center='MEDIAN')
        # Assume A and B are known quaternions
        A = obj.rotation_quaternion
        B = Quaternion(rotation_axis)

        # Calculate C = B_inv * A
        C = B.inverted() @ A

        obj.rotation_quaternion = C
        for o in bpy.context.selected_objects:
            o.select_set(False)
        obj.select_set(True)
        bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
        obj.select_set(False)
        obj.rotation_quaternion = B


    # Calculate intersection point
    A = np.array(obj.location)
    B = np.array(rotation_origin)
    # Calculate BA
    BA = A - B

    # Calculate the direction vector after rotation (initially Z axis)
    initial_dir = Vector((0, 0, 1))
    axis_vector = Quaternion(rotation_axis) @ initial_dir  # Quaternion acts on the initial direction

    # Calculate the projection coefficient t
    t = np.dot(BA, axis_vector) / np.dot(axis_vector, axis_vector)

    # Calculate the intersection point C
    C = B + t * axis_vector

    rotation_origin = Vector(C)

    if len(rotation_axis) == 4:
        rotation_axis = Quaternion(rotation_axis).to_euler(rotation_mode)
    bpy.ops.object.transform_apply(location=True, rotation=False, scale=True)
    
    bpy.context.scene.cursor.location = rotation_origin

    for o in bpy.context.selected_objects:
        o.select_set(False)
    obj.select_set(True)
    bpy.ops.object.origin_set(type='ORIGIN_CURSOR', center='MEDIAN')
    obj.select_set(False)

    bpy.ops.object.empty_add(type='PLAIN_AXES', align='WORLD', location=rotation_origin, rotation=rotation_axis, scale=scale)
    
    empty = bpy.data.objects["Empty"]
    empty.select_set(True)
    bpy.ops.object.origin_set(type='ORIGIN_CURSOR', center='MEDIAN')

    empty.select_set(True)
    bpy.context.view_layer.objects.active = empty

    euler = mathutils.Euler(rotation_axis, rotation_mode)
    orient_matrix = euler.to_matrix()

    if isinstance(z_offset, (int, float)):
        z_offset = (0, 0, z_offset) 

    bpy.ops.transform.translate(value=z_offset, orient_type='LOCAL', orient_matrix=orient_matrix, orient_matrix_type='LOCAL', constraint_axis=(False, False, True))

    rotation_angle = rotation_angle * pi
    bpy.ops.transform.rotate(value=rotation_angle, orient_axis='Z', orient_type='LOCAL', orient_matrix=orient_matrix, orient_matrix_type='LOCAL', constraint_axis=(False, False, True))

    array(obj, fit_type, count, offset_object=empty, merge_threshold=merge_threshold, merge_cap=merge_cap, start_cap=start_cap, end_cap=end_cap, apply=True)

    if delete_empty:
        for o in bpy.context.selected_objects:
            o.select_set(False)
        empty.select_set(True)
        bpy.context.view_layer.objects.active = empty
        bpy.ops.object.delete()

    for o in bpy.context.selected_objects:
        o.select_set(False)
    obj.select_set(True)
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
    obj.select_set(False)
    bpy.context.scene.cursor.location = [0, 0, 0]



def zoom_upper_face(name, factor):
    obj = bpy.data.objects[name]
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='EDIT')  # Switch to edit mode

    # Switch to face selection mode
    bpy.ops.mesh.select_mode(type="FACE")
    bpy.ops.mesh.select_all(action='DESELECT')  # Clear all selection

    # Select the upper surface
    bpy.ops.object.mode_set(mode='OBJECT')  # Switch to object mode
    mesh = obj.data
    for face in mesh.polygons:
        if face.center.z > 0:  # Determine the upper surface
            face.select = True
    bpy.ops.object.mode_set(mode='EDIT')  # Switch back to edit mode

    # Resize the upper surface
    resize_factor_sub= factor
    # resize_factor_sub= random.uniform(1.3, 2.0)
    # resize_factor_sub2= random.uniform(1.5, 2.0)
    resize_factor = (resize_factor_sub,resize_factor_sub,resize_factor_sub)
    bpy.ops.transform.resize(
        value=resize_factor,  # Resize the upper surface proportionally
        orient_type='GLOBAL',
        orient_matrix=((1, 0, 0), (0, 1, 0), (0, 0, 1)),
        orient_matrix_type='GLOBAL',
        mirror=True,
        use_proportional_edit=False,
        proportional_edit_falloff='SMOOTH',
        proportional_size=1
    )

    # Exit edit mode
    bpy.ops.object.mode_set(mode='OBJECT')


def create_wave(name, shift=0, frequency=2, altitude=1, segs=200, length=1, rotation_rad=0.0):
    x = np.linspace(shift, 1+shift+length, segs)
    y = np.sin(2 * np.pi * frequency * x) * altitude
    x = np.linspace(-0.5*length, 0.5*length, segs)
    pts = [[x_, y_, 0] for x_, y_ in zip(x, y)]
    pts_ = []
    c_pts = [0] * (len(pts) * 2)
    create_curve(name, None, pts, handle_type=c_pts)
    bpy.data.objects[name].rotation_euler[2] = rotation_rad * np.pi
    bpy.data.objects[name].select_set(True)
    bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
    return {'name':name, 'shift': shift, 'frequency':frequency, 'altitude':altitude, 'segs':segs, 'length': length, 'rotation_rad': rotation_rad}

