content: # modified from https://github.com/Profactor/continuous-remeshing
  import torch
  import torch.nn.functional as tfunc
  import torch_scatter
  from typing import Tuple
  
  def prepend_dummies(
          vertices:torch.Tensor, #V,D
          faces:torch.Tensor, #F,3 long
      )->Tuple[torch.Tensor,torch.Tensor]:
      """prepend dummy elements to vertices and faces to enable "masked" scatter operations"""
      V,D = vertices.shape
      vertices = torch.concat((torch.full((1,D),fill_value=torch.nan,device=vertices.device),vertices),dim=0)
      faces = torch.concat((torch.zeros((1,3),dtype=torch.long,device=faces.device),faces+1),dim=0)
      return vertices,faces
  
  def remove_dummies(
          vertices:torch.Tensor, #V,D - first vertex all nan and unreferenced
          faces:torch.Tensor, #F,3 long - first face all zeros
      )->Tuple[torch.Tensor,torch.Tensor]:
      """remove dummy elements added with prepend_dummies()"""
      return vertices[1:],faces[1:]-1
  
  
  def calc_edges(
          faces: torch.Tensor,  # F,3 long - first face may be dummy with all zeros
          with_edge_to_face: bool = False
      ) -> Tuple[torch.Tensor, ...]:
      """
      returns Tuple of
      - edges E,2 long, 0 for unused, lower vertex index first
      - face_to_edge F,3 long
      - (optional) edge_to_face shape=E,[left,right],[face,side]
  
      o-<-----e1     e0,e1...edge, e0<e1
      |      /A      L,R....left and right face
      |  L /  |      both triangles ordered counter clockwise
      |  / R  |      normals pointing out of screen
      V/      |      
      e0---->-o     
      """
  
      F = faces.shape[0]
      
      # make full edges, lower vertex index first
      face_edges = torch.stack((faces,faces.roll(-1,1)),dim=-1) #F*3,3,2
      full_edges = face_edges.reshape(F*3,2)
      sorted_edges,_ = full_edges.sort(dim=-1) #F*3,2
  
      # make unique edges
      edges,full_to_unique = torch.unique(input=sorted_edges,sorted=True,return_inverse=True,dim=0) #(E,2),(F*3)
      E = edges.shape[0]
      face_to_edge = full_to_unique.reshape(F,3) #F,3
  
      if not with_edge_to_face:
          return edges, face_to_edge
  
      is_right = full_edges[:,0]!=sorted_edges[:,0] #F*3
      edge_to_face = torch.zeros((E,2,2),dtype=torch.long,device=faces.device) #E,LR=2,S=2
      scatter_src = torch.cartesian_prod(torch.arange(0,F,device=faces.device),torch.arange(0,3,device=faces.device)) #F*3,2
      edge_to_face.reshape(2*E,2).scatter_(dim=0,index=(2*full_to_unique+is_right)[:,None].expand(F*3,2),src=scatter_src) #E,LR=2,S=2
      edge_to_face[0] = 0
      return edges, face_to_edge, edge_to_face
  
  def calc_edge_length(
          vertices:torch.Tensor, #V,3 first may be dummy
          edges:torch.Tensor, #E,2 long, lower vertex index first, (0,0) for unused
          )->torch.Tensor: #E
  
      full_vertices = vertices[edges] #E,2,3
      a,b = full_vertices.unbind(dim=1) #E,3
      return torch.norm(a-b,p=2,dim=-1)
  
  def calc_face_normals(
          vertices:torch.Tensor, #V,3 first vertex may be unreferenced
          faces:torch.Tensor, #F,3 long, first face may be all zero
          normalize:bool=False,
          )->torch.Tensor: #F,3
      """
           n
           |
           c0     corners ordered counterclockwise when
          / \     looking onto surface (in neg normal direction)
        c1---c2
      """
      full_vertices = vertices[faces] #F,C=3,3
      v0,v1,v2 = full_vertices.unbind(dim=1) #F,3
      face_normals = torch.cross(v1-v0,v2-v0, dim=1) #F,3
      if normalize:
          face_normals = tfunc.normalize(face_normals, eps=1e-6, dim=1) 
      return face_normals #F,3
  
  def calc_vertex_normals(
          vertices:torch.Tensor, #V,3 first vertex may be unreferenced
          faces:torch.Tensor, #F,3 long, first face may be all zero
          face_normals:torch.Tensor=None, #F,3, not normalized
          )->torch.Tensor: #F,3
      F = faces.shape[0]
  
      if face_normals is None:
          face_normals = calc_face_normals(vertices,faces)
      
      vertex_normals = torch.zeros((vertices.shape[0],3,3),dtype=vertices.dtype,device=vertices.device) #V,C=3,3
      vertex_normals.scatter_add_(dim=0,index=faces[:,:,None].expand(F,3,3),src=face_normals[:,None,:].expand(F,3,3))
      vertex_normals = vertex_normals.sum(dim=1) #V,3
      return tfunc.normalize(vertex_normals, eps=1e-6, dim=1)
  
  def calc_face_ref_normals(
          faces:torch.Tensor, #F,3 long, 0 for unused
          vertex_normals:torch.Tensor, #V,3 first unused
          normalize:bool=False,
          )->torch.Tensor: #F,3
      """calculate reference normals for face flip detection"""
      full_normals = vertex_normals[faces] #F,C=3,3
      ref_normals = full_normals.sum(dim=1) #F,3
      if normalize:
          ref_normals = tfunc.normalize(ref_normals, eps=1e-6, dim=1)
      return ref_normals
  
  def pack(
          vertices:torch.Tensor, #V,3 first unused and nan
          faces:torch.Tensor, #F,3 long, 0 for unused
          )->Tuple[torch.Tensor,torch.Tensor]: #(vertices,faces), keeps first vertex unused
      """removes unused elements in vertices and faces"""
      V = vertices.shape[0]
      
      # remove unused faces
      used_faces = faces[:,0]!=0
      used_faces[0] = True
      faces = faces[used_faces] #sync
  
      # remove unused vertices
      used_vertices = torch.zeros(V,3,dtype=torch.bool,device=vertices.device)
      used_vertices.scatter_(dim=0,index=faces,value=True,reduce='add') 
      used_vertices = used_vertices.any(dim=1)
      used_vertices[0] = True
      vertices = vertices[used_vertices] #sync
  
      # update used faces
      ind = torch.zeros(V,dtype=torch.long,device=vertices.device)
      V1 = used_vertices.sum()
      ind[used_vertices] =  torch.arange(0,V1,device=vertices.device) #sync
      faces = ind[faces]
  
      return vertices,faces
  
  def split_edges(
          vertices:torch.Tensor, #V,3 first unused
          faces:torch.Tensor, #F,3 long, 0 for unused
          edges:torch.Tensor, #E,2 long 0 for unused, lower vertex index first
          face_to_edge:torch.Tensor, #F,3 long 0 for unused
          splits, #E bool
          pack_faces:bool=True,
          )->Tuple[torch.Tensor,torch.Tensor]: #(vertices,faces)
  
      #   c2                    c2               c...corners = faces
      #    . .                   . .             s...side_vert, 0 means no split
      #    .   .                 .N2 .           S...shrunk_face
      #    .     .               .     .         Ni...new_faces
      #   s2      s1           s2|c2...s1|c1
      #    .        .            .     .  .
      #    .          .          . S .      .
      #    .            .        . .     N1    .
      #   c0...(s0=0)....c1    s0|c0...........c1
      #
      # pseudo-code:
      #   S = [s0|c0,s1|c1,s2|c2] example:[c0,s1,s2]
      #   split = side_vert!=0 example:[False,True,True]
      #   N0 = split[0]*[c0,s0,s2|c2] example:[0,0,0]
      #   N1 = split[1]*[c1,s1,s0|c0] example:[c1,s1,c0]
      #   N2 = split[2]*[c2,s2,s1|c1] example:[c2,s2,s1]
  
      V = vertices.shape[0]
      F = faces.shape[0]
      S = splits.sum().item() #sync
  
      if S==0:
          return vertices,faces
      
      edge_vert = torch.zeros_like(splits, dtype=torch.long) #E
      edge_vert[splits] = torch.arange(V,V+S,dtype=torch.long,device=vertices.device) #E 0 for no split, sync
      side_vert = edge_vert[face_to_edge] #F,3 long, 0 for no split
      split_edges = edges[splits] #S sync
  
      #vertices
      split_vertices = vertices[split_edges].mean(dim=1) #S,3
      vertices = torch.concat((vertices,split_vertices),dim=0)
  
      #faces
      side_split = side_vert!=0 #F,3
      shrunk_faces = torch.where(side_split,side_vert,faces) #F,3 long, 0 for no split
      new_faces = side_split[:,:,None] * torch.stack((faces,side_vert,shrunk_faces.roll(1,dims=-1)),dim=-1) #F,N=3,C=3
      faces = torch.concat((shrunk_faces,new_faces.reshape(F*3,3))) #4F,3
      if pack_faces:
          mask = faces[:,0]!=0
          mask[0] = True
          faces = faces[mask] #F',3 sync
  
      return vertices,faces
  
  def collapse_edges(
          vertices:torch.Tensor, #V,3 first unused
          faces:torch.Tensor, #F,3 long 0 for unused
          edges:torch.Tensor, #E,2 long 0 for unused, lower vertex index first
          priorities:torch.Tensor, #E float
          stable:bool=False, #only for unit testing
          )->Tuple[torch.Tensor,torch.Tensor]: #(vertices,faces)
          
      V = vertices.shape[0]
      
      # check spacing
      _,order = priorities.sort(stable=stable) #E
      rank = torch.zeros_like(order)
      rank[order] = torch.arange(0,len(rank),device=rank.device)
      vert_rank = torch.zeros(V,dtype=torch.long,device=vertices.device) #V
      edge_rank = rank #E
      for i in range(3):
          torch_scatter.scatter_max(src=edge_rank[:,None].expand(-1,2).reshape(-1),index=edges.reshape(-1),dim=0,out=vert_rank)
          edge_rank,_ = vert_rank[edges].max(dim=-1) #E
      candidates = edges[(edge_rank==rank).logical_and_(priorities>0)] #E',2
  
      # check connectivity
      vert_connections = torch.zeros(V,dtype=torch.long,device=vertices.device) #V
      vert_connections[candidates[:,0]] = 1 #start
      edge_connections = vert_connections[edges].sum(dim=-1) #E, edge connected to start
      vert_connections.scatter_add_(dim=0,index=edges.reshape(-1),src=edge_connections[:,None].expand(-1,2).reshape(-1))# one edge from start
      vert_connections[candidates] = 0 #clear start and end
      edge_connections = vert_connections[edges].sum(dim=-1) #E, one or two edges from start
      vert_connections.scatter_add_(dim=0,index=edges.reshape(-1),src=edge_connections[:,None].expand(-1,2).reshape(-1)) #one or two edges from start
      collapses = candidates[vert_connections[candidates[:,1]] <= 2] # E" not more than two connections between start and end
  
      # mean vertices
      vertices[collapses[:,0]] = vertices[collapses].mean(dim=1) 
  
      # update faces
      dest = torch.arange(0,V,dtype=torch.long,device=vertices.device) #V
      dest[collapses[:,1]] = dest[collapses[:,0]]
      faces = dest[faces] #F,3 
      c0,c1,c2 = faces.unbind(dim=-1)
      collapsed = (c0==c1).logical_or_(c1==c2).logical_or_(c0==c2)
      faces[collapsed] = 0
  
      return vertices,faces
  
  def calc_face_collapses(
          vertices:torch.Tensor, #V,3 first unused
          faces:torch.Tensor, #F,3 long, 0 for unused
          edges:torch.Tensor, #E,2 long 0 for unused, lower vertex index first
          face_to_edge:torch.Tensor, #F,3 long 0 for unused
          edge_length:torch.Tensor, #E
          face_normals:torch.Tensor, #F,3
          vertex_normals:torch.Tensor, #V,3 first unused
          min_edge_length:torch.Tensor=None, #V
          area_ratio = 0.5, #collapse if area < min_edge_length**2 * area_ratio
          shortest_probability = 0.8
          )->torch.Tensor: #E edges to collapse
      
      E = edges.shape[0]
      F = faces.shape[0]
  
      # face flips
      ref_normals = calc_face_ref_normals(faces,vertex_normals,normalize=False) #F,3
      face_collapses = (face_normals*ref_normals).sum(dim=-1)<0 #F
      
      # small faces
      if min_edge_length is not None:
          min_face_length = min_edge_length[faces].mean(dim=-1) #F
          min_area = min_face_length**2 * area_ratio #F
          face_collapses.logical_or_(face_normals.norm(dim=-1) < min_area*2) #F
          face_collapses[0] = False
  
      # faces to edges
      face_length = edge_length[face_to_edge] #F,3
  
      if shortest_probability<1:
          #select shortest edge with shortest_probability chance
          randlim = round(2/(1-shortest_probability))
          rand_ind = torch.randint(0,randlim,size=(F,),device=faces.device).clamp_max_(2) #selected edge local index in face
          sort_ind = torch.argsort(face_length,dim=-1,descending=True) #F,3
          local_ind = sort_ind.gather(dim=-1,index=rand_ind[:,None])
      else:
          local_ind = torch.argmin(face_length,dim=-1)[:,None] #F,1 0...2 shortest edge local index in face
      
      edge_ind = face_to_edge.gather(dim=1,index=local_ind)[:,0] #F 0...E selected edge global index
      edge_collapses = torch.zeros(E,dtype=torch.long,device=vertices.device)
      edge_collapses.scatter_add_(dim=0,index=edge_ind,src=face_collapses.long()) 
  
      return edge_collapses.bool()
  
  def flip_edges(
          vertices:torch.Tensor, #V,3 first unused
          faces:torch.Tensor, #F,3 long, first must be 0, 0 for unused
          edges:torch.Tensor, #E,2 long, first must be 0, 0 for unused, lower vertex index first
          edge_to_face:torch.Tensor, #E,[left,right],[face,side]
          with_border:bool=True, #handle border edges (D=4 instead of D=6)
          with_normal_check:bool=True, #check face normal flips
          stable:bool=False, #only for unit testing
          ):
      V = vertices.shape[0]
      E = edges.shape[0]
      device=vertices.device
      vertex_degree = torch.zeros(V,dtype=torch.long,device=device) #V long
      vertex_degree.scatter_(dim=0,index=edges.reshape(E*2),value=1,reduce='add')
      neighbor_corner = (edge_to_face[:,:,1] + 2) % 3 #go from side to corner
      neighbors = faces[edge_to_face[:,:,0],neighbor_corner] #E,LR=2
      edge_is_inside = neighbors.all(dim=-1) #E
  
      if with_border:
          # inside vertices should have D=6, border edges D=4, so we subtract 2 for all inside vertices
          # need to use float for masks in order to use scatter(reduce='multiply')
          vertex_is_inside = torch.ones(V,2,dtype=torch.float32,device=vertices.device) #V,2 float
          src = edge_is_inside.type(torch.float32)[:,None].expand(E,2) #E,2 float
          vertex_is_inside.scatter_(dim=0,index=edges,src=src,reduce='multiply')
          vertex_is_inside = vertex_is_inside.prod(dim=-1,dtype=torch.long) #V long
          vertex_degree -= 2 * vertex_is_inside #V long
  
      neighbor_degrees = vertex_degree[neighbors] #E,LR=2
      edge_degrees = vertex_degree[edges] #E,2
      #
      # loss = Sum_over_affected_vertices((new_degree-6)**2)
      # loss_change = Sum_over_neighbor_vertices((degree+1-6)**2-(degree-6)**2)
      #                   + Sum_over_edge_vertices((degree-1-6)**2-(degree-6)**2)
      #             = 2 * (2 + Sum_over_neighbor_vertices(degree) - Sum_over_edge_vertices(degree))
      #
      loss_change = 2 + neighbor_degrees.sum(dim=-1) - edge_degrees.sum(dim=-1) #E
      candidates = torch.logical_and(loss_change<0, edge_is_inside) #E
      loss_change = loss_change[candidates] #E'
      if loss_change.shape[0]==0:
          return
  
      edges_neighbors = torch.concat((edges[candidates],neighbors[candidates]),dim=-1) #E',4
      _,order = loss_change.sort(descending=True, stable=stable) #E'
      rank = torch.zeros_like(order)
      rank[order] = torch.arange(0,len(rank),device=rank.device)
      vertex_rank = torch.zeros((V,4),dtype=torch.long,device=device) #V,4
      torch_scatter.scatter_max(src=rank[:,None].expand(-1,4),index=edges_neighbors,dim=0,out=vertex_rank)
      vertex_rank,_ = vertex_rank.max(dim=-1) #V
      neighborhood_rank,_ = vertex_rank[edges_neighbors].max(dim=-1) #E'
      flip = rank==neighborhood_rank #E'
  
      if with_normal_check:
          #  cl-<-----e1     e0,e1...edge, e0<e1
          #   |      /A      L,R....left and right face
          #   |  L /  |      both triangles ordered counter clockwise
          #   |  / R  |      normals pointing out of screen
          #   V/      |      
          #   e0---->-cr    
          v = vertices[edges_neighbors] #E",4,3
          v = v - v[:,0:1] #make relative to e0 
          e1 = v[:,1]
          cl = v[:,2]
          cr = v[:,3]
          n = torch.cross(e1,cl,dim=-1) + torch.cross(cr,e1,dim=-1) #sum of old normal vectors 
          flip.logical_and_(torch.sum(n*torch.cross(cr,cl,dim=-1),dim=-1)>0) #first new face
          flip.logical_and_(torch.sum(n*torch.cross(cl-e1,cr-e1,dim=-1),dim=-1)>0) #second new face
  
      flip_edges_neighbors = edges_neighbors[flip] #E",4
      flip_edge_to_face = edge_to_face[candidates,:,0][flip] #E",2
      flip_faces = flip_edges_neighbors[:,[[0,3,2],[1,2,3]]] #E",2,3
      faces.scatter_(dim=0,index=flip_edge_to_face.reshape(-1,1).expand(-1,3),src=flip_faces.reshape(-1,3))
  