diff --git a/◯ᴥᗱᗴᗝИNᗱᗴᙁ⚭ⵙ⚭ᙁᗱᗴИNᗝᗱᗴᴥ◯/2.90/SCRIPTS/ADDONS/TISSUE-MASTER/COLORS_GROUPS_EXCHANGER.PY b/◯ᴥᗱᗴᗝИNᗱᗴᙁ⚭ⵙ⚭ᙁᗱᗴИNᗝᗱᗴᴥ◯/2.90/SCRIPTS/ADDONS/TISSUE-MASTER/COLORS_GROUPS_EXCHANGER.PY
deleted file mode 100644
index dee8bd62..00000000
--- a/◯ᴥᗱᗴᗝИNᗱᗴᙁ⚭ⵙ⚭ᙁᗱᗴИNᗝᗱᗴᴥ◯/2.90/SCRIPTS/ADDONS/TISSUE-MASTER/COLORS_GROUPS_EXCHANGER.PY
+++ /dev/null
@@ -1,3235 +0,0 @@
-# ##### BEGIN GPL LICENSE BLOCK #####
-#
-#  This program is free software; you can redistribute it and/or
-#  modify it under the terms of the GNU General Public License
-#  as published by the Free Software Foundation; either version 2
-#  of the License, or (at your option) any later version.
-#
-#  This program is distributed in the hope that it will be useful,
-#  but WITHOUT ANY WARRANTY; without even the implied warranty of
-#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-#  GNU General Public License for more details.
-#
-#  You should have received a copy of the GNU General Public License
-#  along with this program; if not, write to the Free Software Foundation,
-#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# ##### END GPL LICENSE BLOCK #####
-
-#-------------------------- COLORS / GROUPS EXCHANGER -------------------------#
-#                                                                              #
-# Vertex Color to Vertex Group allow you to convert colors channles to weight  #
-# maps.                                                                        #
-# The main purpose is to use vertex colors to store information when importing #
-# files from other softwares. The script works with the active vertex color    #
-# slot.                                                                        #
-# For use the command "Vertex Clors to Vertex Groups" use the search bar       #
-# (space bar).                                                                 #
-#                                                                              #
-#                          (c)  Alessandro Zomparelli                          #
-#                                     (2017)                                   #
-#                                                                              #
-# http://www.co-de-it.com/                                                     #
-#                                                                              #
-################################################################################
-
-import bpy, bmesh
-import numpy as np
-import math, timeit, time
-from math import *#pi, sin
-from statistics import mean, stdev
-from mathutils import Vector
-from numpy import *
-try: from .numba_functions import numba_reaction_diffusion
-except: pass
-try: import numexpr as ne
-except: pass
-
-from bpy.types import (
-        Operator,
-        Panel,
-        PropertyGroup,
-        )
-
-from bpy.props import (
-    BoolProperty,
-    EnumProperty,
-    FloatProperty,
-    IntProperty,
-    StringProperty,
-    FloatVectorProperty,
-    IntVectorProperty
-)
-
-from .utils import *
-
-def reaction_diffusion_add_handler(self, context):
-    # remove existing handlers
-    old_handlers = []
-    for h in bpy.app.handlers.frame_change_post:
-        if "reaction_diffusion" in str(h):
-            old_handlers.append(h)
-    for h in old_handlers: bpy.app.handlers.frame_change_post.remove(h)
-    # add new handler
-    bpy.app.handlers.frame_change_post.append(reaction_diffusion_def)
-
-class formula_prop(PropertyGroup):
-    name : StringProperty()
-    formula : StringProperty()
-    float_var : FloatVectorProperty(name="", description="", default=(0, 0, 0, 0, 0), size=5)
-    int_var : IntVectorProperty(name="", description="", default=(0, 0, 0, 0, 0), size=5)
-
-class reaction_diffusion_prop(PropertyGroup):
-    run : BoolProperty(default=False, update = reaction_diffusion_add_handler,
-        description='Compute a new iteration on frame changes. Currently is not working during  Render Animation')
-
-    time_steps : bpy.props.IntProperty(
-        name="Steps", default=10, min=0, soft_max=50,
-        description="Number of Steps")
-
-    dt : bpy.props.FloatProperty(
-        name="dt", default=1, min=0, soft_max=0.2,
-        description="Time Step")
-
-    diff_a : bpy.props.FloatProperty(
-        name="Diff A", default=0.1, min=0, soft_max=2, precision=3,
-        description="Diffusion A")
-
-    diff_b : bpy.props.FloatProperty(
-        name="Diff B", default=0.05, min=0, soft_max=2, precision=3,
-        description="Diffusion B")
-
-    f : bpy.props.FloatProperty(
-        name="f", default=0.055, min=0, soft_max=0.5, precision=3,
-        description="Feed Rate")
-
-    k : bpy.props.FloatProperty(
-        name="k", default=0.062, min=0, soft_max=0.5, precision=3,
-        description="Kill Rate")
-
-    diff_mult : bpy.props.FloatProperty(
-        name="Scale", default=1, min=0, soft_max=1, max=2, precision=2,
-        description="Multiplier for the diffusion of both substances")
-
-def compute_formula(ob=None, formula="rx", float_var=(0,0,0,0,0), int_var=(0,0,0,0,0)):
-    verts = ob.data.vertices
-    n_verts = len(verts)
-
-    f1,f2,f3,f4,f5 = float_var
-    i1,i2,i3,i4,i5 = int_var
-
-    do_groups = "w[" in formula
-    do_local = "lx" in formula or "ly" in formula or "lz" in formula
-    do_global = "gx" in formula or "gy" in formula or "gz" in formula
-    do_relative = "rx" in formula or "ry" in formula or "rz" in formula
-    do_normal = "nx" in formula or "ny" in formula or "nz" in formula
-    mat = ob.matrix_world
-
-    for i in range(1000):
-        if "w["+str(i)+"]" in formula and i > len(ob.vertex_groups)-1:
-            return "w["+str(i)+"] not found"
-
-    w = []
-    for i in range(len(ob.vertex_groups)):
-        w.append([])
-        if "w["+str(i)+"]" in formula:
-            vg = ob.vertex_groups[i]
-            for v in verts:
-                try:
-                    w[i].append(vg.weight(v.index))
-                except:
-                    w[i].append(0)
-            w[i] = array(w[i])
-
-    start_time = timeit.default_timer()
-    # compute vertex coordinates
-    if do_local or do_relative or do_global:
-        co = [0]*n_verts*3
-        verts.foreach_get('co', co)
-        np_co = array(co).reshape((n_verts, 3))
-        lx, ly, lz = array(np_co).transpose()
-        if do_relative:
-            rx = np.interp(lx, (lx.min(), lx.max()), (0, +1))
-            ry = np.interp(ly, (ly.min(), ly.max()), (0, +1))
-            rz = np.interp(lz, (lz.min(), lz.max()), (0, +1))
-        if do_global:
-            co = [v.co for v in verts]
-            global_co = []
-            for v in co:
-                global_co.append(mat * v)
-            global_co = array(global_co).reshape((n_verts, 3))
-            gx, gy, gz = array(global_co).transpose()
-    # compute vertex normals
-    if do_normal:
-        normal = [0]*n_verts*3
-        verts.foreach_get('normal', normal)
-        normal = array(normal).reshape((n_verts, 3))
-        nx, ny, nz = array(normal).transpose()
-
-    try:
-        weight = eval(formula)
-        return weight
-    except:
-        return "There is something wrong"
-    print("Weight Formula: " + str(timeit.default_timer() - start_time))
-
-class weight_formula_wiki(bpy.types.Operator):
-    bl_idname = "scene.weight_formula_wiki"
-    bl_label = "Online Documentation"
-    bl_options = {'REGISTER', 'UNDO'}
-
-    def execute(self, context):
-        bpy.ops.wm.url_open(url="https://github.com/alessandro-zomparelli/tissue/wiki/Weight-Tools#weight-formula")
-        return {'FINISHED'}
-
-class weight_formula(bpy.types.Operator):
-    bl_idname = "object.weight_formula"
-    bl_label = "Weight Formula"
-    bl_options = {'REGISTER', 'UNDO'}
-
-    ex = [
-        #'cos(arctan(nx/ny)*6 + sin(rz*30)*0.5)/2 + cos(arctan(nx/ny)*6 - sin(rz*30)*0.5 + pi/2)/2 + 0.5',
-        'cos(arctan(nx/ny)*i1*2 + sin(rz*i3))/i2 + cos(arctan(nx/ny)*i1*2 - sin(rz*i3))/i2 + 0.5',
-        'cos(arctan(nx/ny)*i1*2 + sin(rz*i2))/2 + cos(arctan(nx/ny)*i1*2 - sin(rz*i2))/2',
-        '(sin(arctan(nx/ny)*i1)*sin(nz*i1)+1)/2',
-        'cos(arctan(nx/ny)*f1)',
-        'cos(arctan(lx/ly)*f1 + sin(rz*f2)*f3)',
-        'sin(nx*15)<sin(ny*15)',
-        'cos(ny*rz**2*i1)',
-        'sin(rx*30) > 0',
-        'sin(nz*i1)',
-        'w[0]**2',
-        'sqrt((rx-0.5)**2 + (ry-0.5)**2)*2',
-        'abs(0.5-rz)*2',
-        'rx'
-        ]
-    ex_items = list((s,s,"") for s in ex)
-    ex_items.append(('CUSTOM', "User Formula", ""))
-
-    examples : bpy.props.EnumProperty(
-        items = ex_items, default='CUSTOM', name="Examples")
-
-    old_ex = ""
-
-    formula : bpy.props.StringProperty(
-        name="Formula", default="", description="Formula to Evaluate")
-    bl_description = ("Generate a Vertex Group based on the given formula")
-
-    slider_f01 : bpy.props.FloatProperty(
-        name="f1", default=1, description="Slider")
-    bl_description = ("Slider Float 1")
-    slider_f02 : bpy.props.FloatProperty(
-        name="f2", default=1, description="Slider")
-    bl_description = ("Slider Float 2")
-    slider_f03 : bpy.props.FloatProperty(
-        name="f3", default=1, description="Slider")
-    bl_description = ("Slider Float 3")
-    slider_f04 : bpy.props.FloatProperty(
-        name="f4", default=1, description="Slider")
-    bl_description = ("Slider Float 4")
-    slider_f05 : bpy.props.FloatProperty(
-        name="f5", default=1, description="Slider")
-    bl_description = ("Slider Float 5")
-    slider_i01 : bpy.props.IntProperty(
-        name="i1", default=1, description="Slider")
-    bl_description = ("Slider Integer 1")
-    slider_i02 : bpy.props.IntProperty(
-        name="i2", default=1, description="Slider")
-    bl_description = ("Slider Integer 2")
-    slider_i03 : bpy.props.IntProperty(
-        name="i3", default=1, description="Slider")
-    bl_description = ("Slider Integer 3")
-    slider_i04 : bpy.props.IntProperty(
-        name="i4", default=1, description="Slider")
-    bl_description = ("Slider Integer 4")
-    slider_i05 : bpy.props.IntProperty(
-        name="i5", default=1, description="Slider")
-    bl_description = ("Slider Integer 5")
-
-    def invoke(self, context, event):
-        return context.window_manager.invoke_props_dialog(self, width=350)
-
-    def draw(self, context):
-        layout = self.layout
-        #layout.label(text="Examples")
-        layout.prop(self, "examples", text="Examples")
-        #if self.examples == 'CUSTOM':
-        layout.label(text="Formula")
-        layout.prop(self, "formula", text="")
-        #try: self.examples = self.formula
-        #except: pass
-
-        if self.examples != self.old_ex and self.examples != 'CUSTOM':
-            self.formula = self.examples
-            self.old_ex = self.examples
-        elif self.formula != self.examples:
-            self.examples = 'CUSTOM'
-        formula = self.formula
-
-        layout.separator()
-        if "f1" in formula: layout.prop(self, "slider_f01")
-        if "f2" in formula: layout.prop(self, "slider_f02")
-        if "f3" in formula: layout.prop(self, "slider_f03")
-        if "f4" in formula: layout.prop(self, "slider_f04")
-        if "f5" in formula: layout.prop(self, "slider_f05")
-        if "i1" in formula: layout.prop(self, "slider_i01")
-        if "i2" in formula: layout.prop(self, "slider_i02")
-        if "i3" in formula: layout.prop(self, "slider_i03")
-        if "i4" in formula: layout.prop(self, "slider_i04")
-        if "i5" in formula: layout.prop(self, "slider_i05")
-
-        layout.label(text="Variables (for each vertex):")
-        layout.label(text="lx, ly, lz: Local Coordinates", icon='ORIENTATION_LOCAL')
-        layout.label(text="gx, gy, gz: Global Coordinates", icon='WORLD')
-        layout.label(text="rx, ry, rz: Local Coordinates (0 to 1)", icon='NORMALIZE_FCURVES')
-        layout.label(text="nx, ny, nz: Normal Coordinates", icon='SNAP_NORMAL')
-        layout.label(text="w[0], w[1], w[2], ... : Vertex Groups", icon="GROUP_VERTEX")
-        layout.separator()
-        layout.label(text="f1, f2, f3, f4, f5: Float Sliders", icon='MOD_HUE_SATURATION')#PROPERTIES
-        layout.label(text="i1, i2, i3, i4, i5: Integer Sliders", icon='MOD_HUE_SATURATION')
-        layout.separator()
-        #layout.label(text="All mathematical functions are based on Numpy", icon='INFO')
-        #layout.label(text="https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.math.html", icon='INFO')
-        layout.operator("scene.weight_formula_wiki", icon="HELP")
-        #layout.label(text="(where 'i' is the index of the Vertex Group)")
-
-    def execute(self, context):
-        ob = bpy.context.active_object
-        n_verts = len(ob.data.vertices)
-        #if self.examples == 'CUSTOM':
-        #    formula = self.formula
-        #else:
-        #self.formula = self.examples
-        #    formula = self.examples
-
-        #f1, f2, f3, f4, f5 = self.slider_f01, self.slider_f02, self.slider_f03, self.slider_f04, self.slider_f05
-        #i1, i2, i3, i4, i5 = self.slider_i01, self.slider_i02, self.slider_i03, self.slider_i04, self.slider_i05
-        f_sliders = self.slider_f01, self.slider_f02, self.slider_f03, self.slider_f04, self.slider_f05
-        i_sliders = self.slider_i01, self.slider_i02, self.slider_i03, self.slider_i04, self.slider_i05
-
-        if self.examples != self.old_ex and self.examples != 'CUSTOM':
-            self.formula = self.examples
-            self.old_ex = self.examples
-        elif self.formula != self.examples:
-            self.examples = 'CUSTOM'
-        formula = self.formula
-
-        if formula == "": return {'FINISHED'}
-        vertex_group_name = "Formula " + formula
-        ob.vertex_groups.new(name=vertex_group_name)
-
-        weight = compute_formula(ob, formula=formula, float_var=f_sliders, int_var=i_sliders)
-        if type(weight) == str:
-            self.report({'ERROR'}, weight)
-            return {'CANCELLED'}
-
-        #start_time = timeit.default_timer()
-        weight = nan_to_num(weight)
-        if type(weight) == int or type(weight) == float:
-            for i in range(n_verts):
-                ob.vertex_groups[-1].add([i], weight, 'REPLACE')
-        elif type(weight) == ndarray:
-            for i in range(n_verts):
-                ob.vertex_groups[-1].add([i], weight[i], 'REPLACE')
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-
-        # Store formula settings
-        new_formula = ob.formula_settings.add()
-        new_formula.name = ob.vertex_groups[-1].name
-        new_formula.formula = formula
-        new_formula.int_var = i_sliders
-        new_formula.float_var = f_sliders
-
-        #for f in ob.formula_settings:
-        #    print(f.name, f.formula, f.int_var, f.float_var)
-        return {'FINISHED'}
-
-class _weight_laplacian(bpy.types.Operator):
-    bl_idname = "object._weight_laplacian"
-    bl_label = "Weight Laplacian"
-    bl_description = ("Compute the Vertex Group Laplacian")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    bounds : bpy.props.EnumProperty(
-        items=(('MANUAL', "Manual Bounds", ""),
-            ('POSITIVE', "Positive Only", ""),
-            ('NEGATIVE', "Negative Only", ""),
-            ('AUTOMATIC', "Automatic Bounds", "")),
-        default='AUTOMATIC', name="Bounds")
-
-    mode : bpy.props.EnumProperty(
-        items=(('LENGTH', "Length Weight", ""),
-            ('SIMPLE', "Simple", "")),
-        default='SIMPLE', name="Evaluation Mode")
-
-    min_def : bpy.props.FloatProperty(
-        name="Min", default=0, soft_min=-1, soft_max=0,
-        description="Laplacian value with 0 weight")
-
-    max_def : bpy.props.FloatProperty(
-        name="Max", default=0.5, soft_min=0, soft_max=5,
-        description="Laplacian value with 1 weight")
-
-    bounds_string = ""
-
-    frame = None
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def draw(self, context):
-        layout = self.layout
-        col = layout.column(align=True)
-        col.label(text="Evaluation Mode")
-        col.prop(self, "mode", text="")
-        col.label(text="Bounds")
-        col.prop(self, "bounds", text="")
-        if self.bounds == 'MANUAL':
-            col.label(text="Strain Rate \u03B5:")
-            col.prop(self, "min_def")
-            col.prop(self, "max_def")
-        col.label(text="\u03B5" + ": from " + self.bounds_string)
-
-
-    def execute(self, context):
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-
-        group_id = ob.vertex_groups.active_index
-        input_group = ob.vertex_groups[group_id].name
-
-        group_name = "Laplacian"
-        ob.vertex_groups.new(name=group_name)
-        me = ob.data
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        bm.edges.ensure_lookup_table()
-
-        # store weight values
-        weight = []
-        for v in me.vertices:
-            try:
-                weight.append(ob.vertex_groups[input_group].weight(v.index))
-            except:
-                weight.append(0)
-
-        n_verts = len(bm.verts)
-        lap = [0]*n_verts
-        for e in bm.edges:
-            if self.mode == 'LENGTH':
-                length = e.calc_length()
-                if length == 0: continue
-                id0 = e.verts[0].index
-                id1 = e.verts[1].index
-                lap[id0] += weight[id1]/length - weight[id0]/length
-                lap[id1] += weight[id0]/length - weight[id1]/length
-            else:
-                id0 = e.verts[0].index
-                id1 = e.verts[1].index
-                lap[id0] += weight[id1] - weight[id0]
-                lap[id1] += weight[id0] - weight[id1]
-
-        mean_lap = mean(lap)
-        stdev_lap = stdev(lap)
-        filter_lap = [i for i in lap if mean_lap-2*stdev_lap < i < mean_lap+2*stdev_lap]
-        if self.bounds == 'MANUAL':
-            min_def = self.min_def
-            max_def = self.max_def
-        elif self.bounds == 'AUTOMATIC':
-            min_def = min(filter_lap)
-            max_def = max(filter_lap)
-            self.min_def = min_def
-            self.max_def = max_def
-        elif self.bounds == 'NEGATIVE':
-            min_def = 0
-            max_def = min(filter_lap)
-            self.min_def = min_def
-            self.max_def = max_def
-        elif self.bounds == 'POSITIVE':
-            min_def = 0
-            max_def = max(filter_lap)
-            self.min_def = min_def
-            self.max_def = max_def
-        delta_def = max_def - min_def
-
-        # check undeformed errors
-        if delta_def == 0: delta_def = 0.0001
-
-        for i in range(len(lap)):
-            val = (lap[i]-min_def)/delta_def
-            if val > 0.7: print(str(val) + " " + str(lap[i]))
-            #val = weight[i] + 0.2*lap[i]
-            ob.vertex_groups[-1].add([i], val, 'REPLACE')
-        self.bounds_string = str(round(min_def,2)) + " to " + str(round(max_def,2))
-        ob.vertex_groups[-1].name = group_name + " " + self.bounds_string
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-class weight_laplacian(bpy.types.Operator):
-    bl_idname = "object.weight_laplacian"
-    bl_label = "Weight Laplacian"
-    bl_description = ("Compute the Vertex Group Laplacian")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    steps : bpy.props.IntProperty(
-        name="Steps", default=10, min=0, soft_max=50,
-        description="Number of Steps")
-
-    dt : bpy.props.FloatProperty(
-        name="dt", default=0.2, min=0, soft_max=0.2,
-        description="Time Step")
-
-    diff_a : bpy.props.FloatProperty(
-        name="Diff A", default=1, min=0, soft_max=2,
-        description="Diffusion A")
-
-    diff_b : bpy.props.FloatProperty(
-        name="Diff B", default=0.5, min=0, soft_max=2,
-        description="Diffusion B")
-
-    f : bpy.props.FloatProperty(
-        name="f", default=0.055, min=0, soft_max=0.5,
-        description="Feed Rate")
-
-    k : bpy.props.FloatProperty(
-        name="k", default=0.062, min=0, soft_max=0.5,
-        description="Kill Rate")
-
-    diff_mult : bpy.props.FloatProperty(
-        name="Scale", default=1, min=0, soft_max=1, max=2, precision=2,
-        description="Multiplier for the diffusion of both substances")
-
-    bounds_string = ""
-
-    frame = None
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-
-    def execute(self, context):
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-
-        me = ob.data
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        bm.edges.ensure_lookup_table()
-
-        # store weight values
-        a = []
-        b = []
-        for v in me.vertices:
-            try:
-                a.append(ob.vertex_groups["A"].weight(v.index))
-            except:
-                a.append(0)
-            try:
-                b.append(ob.vertex_groups["B"].weight(v.index))
-            except:
-                b.append(0)
-
-        a = array(a)
-        b = array(b)
-        f = self.f
-        k = self.k
-        diff_a = self.diff_a * self.diff_mult
-        diff_b = self.diff_b * self.diff_mult
-        dt = self.dt
-
-        # initialize
-        n_verts = len(bm.verts)
-        # find max number of edges for vertex
-        max_edges = 0
-        n_neighbors = []
-        id_neighbors = []
-        for v in bm.verts:
-            n_edges = len(v.link_edges)
-            max_edges = max(max_edges, n_edges)
-            n_neighbors.append(n_edges)
-            neighbors = []
-            for e in link_edges:
-                for v1 in e.verts:
-                    if v != v1: neighbors.append(v1.index)
-            id_neighbors.append(neighbors)
-        n_neighbors = array(n_neighbors)
-
-
-        a = [[] for i in range(n_verts)]
-        lap_map = []
-
-        for e in bm.edges:
-            id0 = e.verts[0].index
-            id1 = e.verts[1].index
-            lap_map[id0].append(id1)
-            lap_map[id1].append(id0)
-
-        e1 = array(e1)
-        e2 = array(e2)
-        lap_a = a[e1]
-
-        for i in range(self.steps):
-
-            lap_a = zeros((n_verts))#[0]*n_verts
-            lap_b = zeros((n_verts))#[0]*n_verts
-            for e in bm.edges:
-                id0 = e.verts[0].index
-                id1 = e.verts[1].index
-                lap_a[id0] += a[id1] - a[id0]
-                lap_a[id1] += a[id0] - a[id1]
-                lap_b[id0] += b[id1] - b[id0]
-                lap_b[id1] += b[id0] - b[id1]
-            ab2 = a*b**2
-            a += (diff_a*lap_a - ab2 + f*(1-a))*dt
-            b += (diff_b*lap_b + ab2 - (k+f)*b)*dt
-
-        for i in range(n_verts):
-            ob.vertex_groups['A'].add([i], a[i], 'REPLACE')
-            ob.vertex_groups['B'].add([i], b[i], 'REPLACE')
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-
-class reaction_diffusion(bpy.types.Operator):
-    bl_idname = "object.reaction_diffusion"
-    bl_label = "Reaction Diffusion"
-    bl_description = ("Run a Reaction-Diffusion based on existing Vertex Groups: A and B")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    steps : bpy.props.IntProperty(
-        name="Steps", default=10, min=0, soft_max=50,
-        description="Number of Steps")
-
-    dt : bpy.props.FloatProperty(
-        name="dt", default=0.2, min=0, soft_max=0.2,
-        description="Time Step")
-
-    diff_a : bpy.props.FloatProperty(
-        name="Diff A", default=1, min=0, soft_max=2,
-        description="Diffusion A")
-
-    diff_b : bpy.props.FloatProperty(
-        name="Diff B", default=0.5, min=0, soft_max=2,
-        description="Diffusion B")
-
-    f : bpy.props.FloatProperty(
-        name="f", default=0.055, min=0, soft_max=0.5,
-        description="Feed Rate")
-
-    k : bpy.props.FloatProperty(
-        name="k", default=0.062, min=0, soft_max=0.5,
-        description="Kill Rate")
-
-    bounds_string = ""
-
-    frame = None
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-
-    def execute(self, context):
-        #bpy.app.handlers.frame_change_post.remove(reaction_diffusion_def)
-        reaction_diffusion_add_handler(self, context)
-        set_animatable_fix_handler(self, context)
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-
-        me = ob.data
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        bm.edges.ensure_lookup_table()
-
-        # store weight values
-        a = []
-        b = []
-        for v in me.vertices:
-            try:
-                a.append(ob.vertex_groups["A"].weight(v.index))
-            except:
-                a.append(0)
-            try:
-                b.append(ob.vertex_groups["B"].weight(v.index))
-            except:
-                b.append(0)
-
-        a = array(a)
-        b = array(b)
-        f = self.f
-        k = self.k
-        diff_a = self.diff_a
-        diff_b = self.diff_b
-        dt = self.dt
-        n_verts = len(bm.verts)
-
-        for i in range(self.steps):
-
-            lap_a = zeros((n_verts))#[0]*n_verts
-            lap_b = zeros((n_verts))#[0]*n_verts
-            for e in bm.edges:
-                id0 = e.verts[0].index
-                id1 = e.verts[1].index
-                lap_a[id0] += a[id1] - a[id0]
-                lap_a[id1] += a[id0] - a[id1]
-                lap_b[id0] += b[id1] - b[id0]
-                lap_b[id1] += b[id0] - b[id1]
-            ab2 = a*b**2
-            a += (diff_a*lap_a - ab2 + f*(1-a))*dt
-            b += (diff_b*lap_b + ab2 - (k+f)*b)*dt
-
-            for i in range(n_verts):
-                ob.vertex_groups['A'].add([i], a[i], 'REPLACE')
-                ob.vertex_groups['B'].add([i], b[i], 'REPLACE')
-            ob.vertex_groups.update()
-            ob.data.update()
-
-            bpy.ops.wm.redraw_timer(type='DRAW_WIN_SWAP', iterations=1)
-
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-
-class edges_deformation(bpy.types.Operator):
-    bl_idname = "object.edges_deformation"
-    bl_label = "Edges Deformation"
-    bl_description = ("Compute Weight based on the deformation of edges"+
-        "according to visible modifiers.")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    bounds : bpy.props.EnumProperty(
-        items=(('MANUAL', "Manual Bounds", ""),
-            ('COMPRESSION', "Compressed Only", ""),
-            ('TENSION', "Extended Only", ""),
-            ('AUTOMATIC', "Automatic Bounds", "")),
-        default='AUTOMATIC', name="Bounds")
-
-    mode : bpy.props.EnumProperty(
-        items=(('MAX', "Max Deformation", ""),
-            ('MEAN', "Average Deformation", "")),
-        default='MEAN', name="Evaluation Mode")
-
-    min_def : bpy.props.FloatProperty(
-        name="Min", default=0, soft_min=-1, soft_max=0,
-        description="Deformations with 0 weight")
-
-    max_def : bpy.props.FloatProperty(
-        name="Max", default=0.5, soft_min=0, soft_max=5,
-        description="Deformations with 1 weight")
-
-    bounds_string = ""
-
-    frame = None
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.modifiers) > 0
-
-    def draw(self, context):
-        layout = self.layout
-        col = layout.column(align=True)
-        col.label(text="Evaluation Mode")
-        col.prop(self, "mode", text="")
-        col.label(text="Bounds")
-        col.prop(self, "bounds", text="")
-        if self.bounds == 'MANUAL':
-            col.label(text="Strain Rate \u03B5:")
-            col.prop(self, "min_def")
-            col.prop(self, "max_def")
-        col.label(text="\u03B5" + ": from " + self.bounds_string)
-
-    def execute(self, context):
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-
-        # check if the object is Cloth or Softbody
-        physics = False
-        for m in ob.modifiers:
-            if m.type == 'CLOTH' or m.type == 'SOFT_BODY':
-                physics = True
-                if context.scene.frame_current == 1 and self.frame != None:
-                    context.scene.frame_current = self.frame
-                break
-        if not physics: self.frame = None
-
-        if self.mode == 'MEAN': group_name = "Average Deformation"
-        elif self.mode == 'MAX': group_name = "Max Deformation"
-        ob.vertex_groups.new(name=group_name)
-        me0 = ob.data
-
-        me = simple_to_mesh(ob) #ob.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-        if len(me.vertices) != len(me0.vertices) or len(me.edges) != len(me0.edges):
-            self.report({'ERROR'}, "The topology of the object should be" +
-                "unaltered")
-            return {'CANCELLED'}
-
-        bm0 = bmesh.new()
-        bm0.from_mesh(me0)
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        deformations = []
-        for e0, e in zip(bm0.edges, bm.edges):
-            try:
-                l0 = e0.calc_length()
-                l1 = e.calc_length()
-                epsilon = (l1 - l0)/l0
-                deformations.append(epsilon)
-            except: deformations.append(1)
-        v_deformations = []
-        for v in bm.verts:
-            vdef = []
-            for e in v.link_edges:
-                vdef.append(deformations[e.index])
-            if self.mode == 'MEAN': v_deformations.append(mean(vdef))
-            elif self.mode == 'MAX': v_deformations.append(max(vdef, key=abs))
-            #elif self.mode == 'MIN': v_deformations.append(min(vdef, key=abs))
-
-        if self.bounds == 'MANUAL':
-            min_def = self.min_def
-            max_def = self.max_def
-        elif self.bounds == 'AUTOMATIC':
-            min_def = min(v_deformations)
-            max_def = max(v_deformations)
-            self.min_def = min_def
-            self.max_def = max_def
-        elif self.bounds == 'COMPRESSION':
-            min_def = 0
-            max_def = min(v_deformations)
-            self.min_def = min_def
-            self.max_def = max_def
-        elif self.bounds == 'TENSION':
-            min_def = 0
-            max_def = max(v_deformations)
-            self.min_def = min_def
-            self.max_def = max_def
-        delta_def = max_def - min_def
-
-        # check undeformed errors
-        if delta_def == 0:
-            if self.bounds == 'MANUAL':
-                delta_def = 0.0001
-            else:
-                message = "The object doesn't have deformations."
-                if physics:
-                    message = message + ("\nIf you are using Physics try to " +
-                        "save it in the cache before.")
-                self.report({'ERROR'}, message)
-                return {'CANCELLED'}
-        else:
-            if physics:
-                self.frame = context.scene.frame_current
-
-        for i in range(len(v_deformations)):
-            weight = (v_deformations[i] - min_def)/delta_def
-            ob.vertex_groups[-1].add([i], weight, 'REPLACE')
-        self.bounds_string = str(round(min_def,2)) + " to " + str(round(max_def,2))
-        ob.vertex_groups[-1].name = group_name + " " + self.bounds_string
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        bpy.data.meshes.remove(me)
-        return {'FINISHED'}
-
-class edges_bending(bpy.types.Operator):
-    bl_idname = "object.edges_bending"
-    bl_label = "Edges Bending"
-    bl_description = ("Compute Weight based on the bending of edges"+
-        "according to visible modifiers.")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    bounds : bpy.props.EnumProperty(
-        items=(('MANUAL', "Manual Bounds", ""),
-            ('POSITIVE', "Positive Only", ""),
-            ('NEGATIVE', "Negative Only", ""),
-            ('UNSIGNED', "Absolute Bending", ""),
-            ('AUTOMATIC', "Signed Bending", "")),
-        default='AUTOMATIC', name="Bounds")
-
-    min_def : bpy.props.FloatProperty(
-        name="Min", default=-10, soft_min=-45, soft_max=45,
-        description="Deformations with 0 weight")
-
-    max_def : bpy.props.FloatProperty(
-        name="Max", default=10, soft_min=-45, soft_max=45,
-        description="Deformations with 1 weight")
-
-    bounds_string = ""
-    frame = None
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.modifiers) > 0
-
-    def draw(self, context):
-        layout = self.layout
-        layout.label(text="Bounds")
-        layout.prop(self, "bounds", text="")
-        if self.bounds == 'MANUAL':
-            layout.prop(self, "min_def")
-            layout.prop(self, "max_def")
-
-    def execute(self, context):
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-
-        group_name = "Edges Bending"
-        ob.vertex_groups.new(name=group_name)
-
-        # check if the object is Cloth or Softbody
-        physics = False
-        for m in ob.modifiers:
-            if m.type == 'CLOTH' or m.type == 'SOFT_BODY':
-                physics = True
-                if context.scene.frame_current == 1 and self.frame != None:
-                    context.scene.frame_current = self.frame
-                break
-        if not physics: self.frame = None
-
-        #ob.data.update()
-        #context.scene.update()
-        me0 = ob.data
-        me = simple_to_mesh(ob) #ob.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-        if len(me.vertices) != len(me0.vertices) or len(me.edges) != len(me0.edges):
-            self.report({'ERROR'}, "The topology of the object should be" +
-                "unaltered")
-        bm0 = bmesh.new()
-        bm0.from_mesh(me0)
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        deformations = []
-        for e0, e in zip(bm0.edges, bm.edges):
-            try:
-                ang = e.calc_face_angle_signed()
-                ang0 = e0.calc_face_angle_signed()
-                if self.bounds == 'UNSIGNED':
-                    deformations.append(abs(ang-ang0))
-                else:
-                    deformations.append(ang-ang0)
-            except: deformations.append(0)
-        v_deformations = []
-        for v in bm.verts:
-            vdef = []
-            for e in v.link_edges:
-                vdef.append(deformations[e.index])
-            v_deformations.append(mean(vdef))
-        if self.bounds == 'MANUAL':
-            min_def = radians(self.min_def)
-            max_def = radians(self.max_def)
-        elif self.bounds == 'AUTOMATIC':
-            min_def = min(v_deformations)
-            max_def = max(v_deformations)
-        elif self.bounds == 'POSITIVE':
-            min_def = 0
-            max_def = min(v_deformations)
-        elif self.bounds == 'NEGATIVE':
-            min_def = 0
-            max_def = max(v_deformations)
-        elif self.bounds == 'UNSIGNED':
-            min_def = 0
-            max_def = max(v_deformations)
-        delta_def = max_def - min_def
-
-        # check undeformed errors
-        if delta_def == 0:
-            if self.bounds == 'MANUAL':
-                delta_def = 0.0001
-            else:
-                message = "The object doesn't have deformations."
-                if physics:
-                    message = message + ("\nIf you are using Physics try to " +
-                        "save it in the cache before.")
-                self.report({'ERROR'}, message)
-                return {'CANCELLED'}
-        else:
-            if physics:
-                self.frame = context.scene.frame_current
-
-        for i in range(len(v_deformations)):
-            weight = (v_deformations[i] - min_def)/delta_def
-            ob.vertex_groups[-1].add([i], weight, 'REPLACE')
-        self.bounds_string = str(round(min_def,2)) + " to " + str(round(max_def,2))
-        ob.vertex_groups[-1].name = group_name + " " + self.bounds_string
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        bpy.data.meshes.remove(me)
-        return {'FINISHED'}
-
-class weight_contour_displace(bpy.types.Operator):
-    bl_idname = "object.weight_contour_displace"
-    bl_label = "Contour Displace"
-    bl_description = ("")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    use_modifiers : bpy.props.BoolProperty(
-        name="Use Modifiers", default=True,
-        description="Apply all the modifiers")
-    min_iso : bpy.props.FloatProperty(
-        name="Min Iso Value", default=0.49, min=0, max=1,
-        description="Threshold value")
-    max_iso : bpy.props.FloatProperty(
-        name="Max Iso Value", default=0.51, min=0, max=1,
-        description="Threshold value")
-    n_cuts : bpy.props.IntProperty(
-        name="Cuts", default=2, min=1, soft_max=10,
-        description="Number of cuts in the selected range of values")
-    bool_displace : bpy.props.BoolProperty(
-        name="Add Displace", default=True, description="Add Displace Modifier")
-    bool_flip : bpy.props.BoolProperty(
-        name="Flip", default=False, description="Flip Output Weight")
-
-    weight_mode : bpy.props.EnumProperty(
-        items=[('Remapped', 'Remapped', 'Remap values'),
-               ('Alternate', 'Alternate', 'Alternate 0 and 1'),
-               ('Original', 'Original', 'Keep original Vertex Group')],
-        name="Weight", description="Choose how to convert vertex group",
-        default="Remapped", options={'LIBRARY_EDITABLE'})
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def invoke(self, context, event):
-        return context.window_manager.invoke_props_dialog(self, width=350)
-
-    def execute(self, context):
-        start_time = timeit.default_timer()
-        try:
-            check = bpy.context.object.vertex_groups[0]
-        except:
-            self.report({'ERROR'}, "The object doesn't have Vertex Groups")
-            return {'CANCELLED'}
-
-        ob0 = bpy.context.object
-
-        group_id = ob0.vertex_groups.active_index
-        vertex_group_name = ob0.vertex_groups[group_id].name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.mesh.select_all(action='SELECT')
-        bpy.ops.object.mode_set(mode='OBJECT')
-        if self.use_modifiers:
-            #me0 = ob0.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-            me0 = simple_to_mesh(ob0)
-        else:
-            me0 = ob0.data.copy()
-
-        # generate new bmesh
-        bm = bmesh.new()
-        bm.from_mesh(me0)
-        bm.verts.ensure_lookup_table()
-        bm.edges.ensure_lookup_table()
-        bm.faces.ensure_lookup_table()
-
-        # store weight values
-        weight = []
-        ob = bpy.data.objects.new("temp", me0)
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-        for v in me0.vertices:
-            try:
-                weight.append(ob.vertex_groups[vertex_group_name].weight(v.index))
-            except:
-                weight.append(0)
-
-        # define iso values
-        iso_values = []
-        for i_cut in range(self.n_cuts):
-            delta_iso = abs(self.max_iso - self.min_iso)
-            min_iso = min(self.min_iso, self.max_iso)
-            max_iso = max(self.min_iso, self.max_iso)
-            if delta_iso == 0: iso_val = min_iso
-            elif self.n_cuts > 1: iso_val = i_cut/(self.n_cuts-1)*delta_iso + min_iso
-            else: iso_val = (self.max_iso + self.min_iso)/2
-            iso_values.append(iso_val)
-
-        # Start Cuts Iterations
-        filtered_edges = bm.edges
-        for iso_val in iso_values:
-            delete_edges = []
-
-            faces_mask = []
-            for f in bm.faces:
-                w_min = 2
-                w_max = 2
-                for v in f.verts:
-                    w = weight[v.index]
-                    if w_min == 2:
-                        w_max = w_min = w
-                    if w > w_max: w_max = w
-                    if w < w_min: w_min = w
-                    if w_min < iso_val and w_max > iso_val:
-                        faces_mask.append(f)
-                        break
-
-            #link_faces = [[f for f in e.link_faces] for e in bm.edges]
-
-            #faces_todo = [f.select for f in bm.faces]
-            #faces_todo = [True for f in bm.faces]
-            verts = []
-            edges = []
-            edges_id = {}
-            _filtered_edges = []
-            n_verts = len(bm.verts)
-            count = n_verts
-            for e in filtered_edges:
-                #id0 = e.vertices[0]
-                #id1 = e.vertices[1]
-                id0 = e.verts[0].index
-                id1 = e.verts[1].index
-                w0 = weight[id0]
-                w1 = weight[id1]
-
-                if w0 == w1: continue
-                elif w0 > iso_val and w1 > iso_val:
-                    _filtered_edges.append(e)
-                    continue
-                elif w0 < iso_val and w1 < iso_val: continue
-                elif w0 == iso_val or w1 == iso_val:
-                    _filtered_edges.append(e)
-                    continue
-                else:
-                    v0 = bm.verts[id0].co
-                    v1 = bm.verts[id1].co
-                    v = v0.lerp(v1, (iso_val-w0)/(w1-w0))
-                    if e not in delete_edges:
-                        delete_edges.append(e)
-                    verts.append(v)
-                    edges_id[str(id0)+"_"+str(id1)] = count
-                    edges_id[str(id1)+"_"+str(id0)] = count
-                    count += 1
-                    _filtered_edges.append(e)
-            filtered_edges = _filtered_edges
-            splitted_faces = []
-
-            switch = False
-            # splitting faces
-            for f in faces_mask:
-                # create sub-faces slots. Once a new vertex is reached it will
-                # change slot, storing the next vertices for a new face.
-                build_faces = [[],[]]
-                #switch = False
-                verts0 = [v.index for v in f.verts]
-                verts1 = list(verts0)
-                verts1.append(verts1.pop(0)) # shift list
-                for id0, id1 in zip(verts0, verts1):
-
-                    # add first vertex to active slot
-                    build_faces[switch].append(id0)
-
-                    # try to split edge
-                    try:
-                        # check if the edge must be splitted
-                        new_vert = edges_id[str(id0)+"_"+str(id1)]
-                        # add new vertex
-                        build_faces[switch].append(new_vert)
-                        # if there is an open face on the other slot
-                        if len(build_faces[not switch]) > 0:
-                            # store actual face
-                            splitted_faces.append(build_faces[switch])
-                            # reset actual faces and switch
-                            build_faces[switch] = []
-                            # change face slot
-                        switch = not switch
-                        # continue previous face
-                        build_faces[switch].append(new_vert)
-                    except: pass
-                if len(build_faces[not switch]) == 2:
-                    build_faces[not switch].append(id0)
-                if len(build_faces[not switch]) > 2:
-                    splitted_faces.append(build_faces[not switch])
-                # add last face
-                splitted_faces.append(build_faces[switch])
-                #del_faces.append(f.index)
-
-            # adding new vertices
-            for v in verts: new_vert = bm.verts.new(v)
-            bm.verts.index_update()
-            bm.verts.ensure_lookup_table()
-            # adding new faces
-            missed_faces = []
-            added_faces = []
-            for f in splitted_faces:
-                try:
-                    face_verts = [bm.verts[i] for i in f]
-                    new_face = bm.faces.new(face_verts)
-                    for e in new_face.edges:
-                        filtered_edges.append(e)
-                except:
-                    missed_faces.append(f)
-
-            bm.faces.ensure_lookup_table()
-            # updating weight values
-            weight = weight + [iso_val]*len(verts)
-
-            # deleting old edges/faces
-            bm.edges.ensure_lookup_table()
-            for e in delete_edges:
-                bm.edges.remove(e)
-            _filtered_edges = []
-            for e in filtered_edges:
-                if e not in delete_edges: _filtered_edges.append(e)
-            filtered_edges = _filtered_edges
-
-        name = ob0.name + '_ContourDisp'
-        me = bpy.data.meshes.new(name)
-        bm.to_mesh(me)
-        ob = bpy.data.objects.new(name, me)
-
-        # Link object to scene and make active
-        scn = bpy.context.scene
-        bpy.context.collection.objects.link(ob)
-        bpy.context.view_layer.objects.active = ob
-        ob.select_set(True)
-        ob0.select_set(False)
-
-        # generate new vertex group
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-        #ob.vertex_groups.new(name=vertex_group_name)
-
-        all_weight = weight + [iso_val]*len(verts)
-        #mult = 1/(1-iso_val)
-        for id in range(len(all_weight)):
-            #if False: w = (all_weight[id]-iso_val)*mult
-            w = all_weight[id]
-            if self.weight_mode == 'Alternate':
-                direction = self.bool_flip
-                for i in range(len(iso_values)-1):
-                    val0, val1 = iso_values[i], iso_values[i+1]
-                    if val0 < w <= val1:
-                        if direction: w1 = (w-val0)/(val1-val0)
-                        else: w1 = (val1-w)/(val1-val0)
-                    direction = not direction
-                if w < iso_values[0]: w1 = not self.bool_flip
-                if w > iso_values[-1]: w1 = not direction
-            elif self.weight_mode == 'Remapped':
-                if w < min_iso: w1 = 0
-                elif w > max_iso: w1 = 1
-                else: w1 = (w - min_iso)/delta_iso
-            else:
-                if self.bool_flip: w1 = 1-w
-                else: w1 = w
-            ob.vertex_groups[vertex_group_name].add([id], w1, 'REPLACE')
-
-        ob.vertex_groups.active_index = group_id
-
-        # align new object
-        ob.matrix_world = ob0.matrix_world
-
-        # Displace Modifier
-        if self.bool_displace:
-            ob.modifiers.new(type='DISPLACE', name='Displace')
-            ob.modifiers["Displace"].mid_level = 0
-            ob.modifiers["Displace"].strength = 0.1
-            ob.modifiers['Displace'].vertex_group = vertex_group_name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        print("Contour Displace time: " + str(timeit.default_timer() - start_time) + " sec")
-
-        bpy.data.meshes.remove(me0)
-
-        return {'FINISHED'}
-
-class weight_contour_mask(bpy.types.Operator):
-    bl_idname = "object.weight_contour_mask"
-    bl_label = "Contour Mask"
-    bl_description = ("")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    use_modifiers : bpy.props.BoolProperty(
-        name="Use Modifiers", default=True,
-        description="Apply all the modifiers")
-    iso : bpy.props.FloatProperty(
-        name="Iso Value", default=0.5, soft_min=0, soft_max=1,
-        description="Threshold value")
-    bool_solidify : bpy.props.BoolProperty(
-        name="Solidify", default=True, description="Add Solidify Modifier")
-    normalize_weight : bpy.props.BoolProperty(
-        name="Normalize Weight", default=True,
-        description="Normalize weight of remaining vertices")
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def execute(self, context):
-        start_time = timeit.default_timer()
-        try:
-            check = bpy.context.object.vertex_groups[0]
-        except:
-            self.report({'ERROR'}, "The object doesn't have Vertex Groups")
-            return {'CANCELLED'}
-
-        ob0 = bpy.context.object
-
-        iso_val = self.iso
-        group_id = ob0.vertex_groups.active_index
-        vertex_group_name = ob0.vertex_groups[group_id].name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.mesh.select_all(action='SELECT')
-        bpy.ops.object.mode_set(mode='OBJECT')
-        if self.use_modifiers:
-            me0 = simple_to_mesh(ob0)#ob0.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-        else:
-            me0 = ob0.data.copy()
-
-        # generate new bmesh
-        bm = bmesh.new()
-        bm.from_mesh(me0)
-        bm.verts.ensure_lookup_table()
-        bm.edges.ensure_lookup_table()
-        bm.faces.ensure_lookup_table()
-
-        # store weight values
-        weight = []
-        ob = bpy.data.objects.new("temp", me0)
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-        for v in me0.vertices:
-            try:
-                #weight.append(v.groups[vertex_group_name].weight)
-                weight.append(ob.vertex_groups[vertex_group_name].weight(v.index))
-            except:
-                weight.append(0)
-
-        faces_mask = []
-        for f in bm.faces:
-            w_min = 2
-            w_max = 2
-            for v in f.verts:
-                w = weight[v.index]
-                if w_min == 2:
-                    w_max = w_min = w
-                if w > w_max: w_max = w
-                if w < w_min: w_min = w
-                if w_min < iso_val and w_max > iso_val:
-                    faces_mask.append(f)
-                    break
-
-        filtered_edges = bm.edges# me0.edges
-        faces_todo = [f.select for f in bm.faces]
-        verts = []
-        edges = []
-        delete_edges = []
-        edges_id = {}
-        _filtered_edges = []
-        n_verts = len(bm.verts)
-        count = n_verts
-        for e in filtered_edges:
-            id0 = e.verts[0].index
-            id1 = e.verts[1].index
-            w0 = weight[id0]
-            w1 = weight[id1]
-
-            if w0 == w1: continue
-            elif w0 > iso_val and w1 > iso_val:
-                continue
-            elif w0 < iso_val and w1 < iso_val: continue
-            elif w0 == iso_val or w1 == iso_val: continue
-            else:
-                v0 = me0.vertices[id0].co
-                v1 = me0.vertices[id1].co
-                v = v0.lerp(v1, (iso_val-w0)/(w1-w0))
-                delete_edges.append(e)
-                verts.append(v)
-                edges_id[str(id0)+"_"+str(id1)] = count
-                edges_id[str(id1)+"_"+str(id0)] = count
-                count += 1
-
-        splitted_faces = []
-
-        switch = False
-        # splitting faces
-        for f in faces_mask:
-            # create sub-faces slots. Once a new vertex is reached it will
-            # change slot, storing the next vertices for a new face.
-            build_faces = [[],[]]
-            #switch = False
-            verts0 = list(me0.polygons[f.index].vertices)
-            verts1 = list(verts0)
-            verts1.append(verts1.pop(0)) # shift list
-            for id0, id1 in zip(verts0, verts1):
-
-                # add first vertex to active slot
-                build_faces[switch].append(id0)
-
-                # try to split edge
-                try:
-                    # check if the edge must be splitted
-                    new_vert = edges_id[str(id0)+"_"+str(id1)]
-                    # add new vertex
-                    build_faces[switch].append(new_vert)
-                    # if there is an open face on the other slot
-                    if len(build_faces[not switch]) > 0:
-                        # store actual face
-                        splitted_faces.append(build_faces[switch])
-                        # reset actual faces and switch
-                        build_faces[switch] = []
-                        # change face slot
-                    switch = not switch
-                    # continue previous face
-                    build_faces[switch].append(new_vert)
-                except: pass
-            if len(build_faces[not switch]) == 2:
-                build_faces[not switch].append(id0)
-            if len(build_faces[not switch]) > 2:
-                splitted_faces.append(build_faces[not switch])
-            # add last face
-            splitted_faces.append(build_faces[switch])
-
-        # adding new vertices
-        for v in verts: bm.verts.new(v)
-        bm.verts.ensure_lookup_table()
-
-        # deleting old edges/faces
-        bm.edges.ensure_lookup_table()
-        remove_edges = []
-        for e in delete_edges: bm.edges.remove(e)
-
-        bm.verts.ensure_lookup_table()
-        # adding new faces
-        missed_faces = []
-        for f in splitted_faces:
-            try:
-                face_verts = [bm.verts[i] for i in f]
-                bm.faces.new(face_verts)
-            except:
-                missed_faces.append(f)
-
-        # Mask geometry
-        if(True):
-            all_weight = weight + [iso_val+0.0001]*len(verts)
-            weight = []
-            for w, v in zip(all_weight, bm.verts):
-                if w < iso_val: bm.verts.remove(v)
-                else: weight.append(w)
-
-        # Create mesh and object
-        name = ob0.name + '_ContourMask_{:.3f}'.format(iso_val)
-        me = bpy.data.meshes.new(name)
-        bm.to_mesh(me)
-        ob = bpy.data.objects.new(name, me)
-
-        # Link object to scene and make active
-        scn = bpy.context.scene
-        bpy.context.collection.objects.link(ob)
-        bpy.context.view_layer.objects.active = ob
-        ob.select_set(True)
-        ob0.select_set(False)
-
-        # generate new vertex group
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-
-        if iso_val != 1: mult = 1/(1-iso_val)
-        else: mult = 1
-        for id in range(len(weight)):
-            if self.normalize_weight: w = (weight[id]-iso_val)*mult
-            else: w = weight[id]
-            ob.vertex_groups[vertex_group_name].add([id], w, 'REPLACE')
-        ob.vertex_groups.active_index = group_id
-
-        # align new object
-        ob.matrix_world = ob0.matrix_world
-
-        # Add Solidify
-        if self.bool_solidify and True:
-            ob.modifiers.new(type='SOLIDIFY', name='Solidify')
-            ob.modifiers['Solidify'].thickness = 0.05
-            ob.modifiers['Solidify'].offset = 0
-            ob.modifiers['Solidify'].vertex_group = vertex_group_name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        print("Contour Mask time: " + str(timeit.default_timer() - start_time) + " sec")
-
-        bpy.data.meshes.remove(me0)
-
-        return {'FINISHED'}
-
-
-class weight_contour_mask_wip(bpy.types.Operator):
-    bl_idname = "object.weight_contour_mask"
-    bl_label = "Contour Mask"
-    bl_description = ("")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    use_modifiers : bpy.props.BoolProperty(
-        name="Use Modifiers", default=True,
-        description="Apply all the modifiers")
-    iso : bpy.props.FloatProperty(
-        name="Iso Value", default=0.5, soft_min=0, soft_max=1,
-        description="Threshold value")
-    bool_solidify : bpy.props.BoolProperty(
-        name="Solidify", default=True, description="Add Solidify Modifier")
-    normalize_weight : bpy.props.BoolProperty(
-        name="Normalize Weight", default=True,
-        description="Normalize weight of remaining vertices")
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def execute(self, context):
-        start_time = timeit.default_timer()
-        try:
-            check = bpy.context.object.vertex_groups[0]
-        except:
-            self.report({'ERROR'}, "The object doesn't have Vertex Groups")
-            return {'CANCELLED'}
-
-        ob0 = bpy.context.object
-
-        iso_val = self.iso
-        group_id = ob0.vertex_groups.active_index
-        vertex_group_name = ob0.vertex_groups[group_id].name
-
-        #bpy.ops.object.mode_set(mode='EDIT')
-        #bpy.ops.mesh.select_all(action='SELECT')
-        #bpy.ops.object.mode_set(mode='OBJECT')
-        if self.use_modifiers:
-            me0 = simple_to_mesh(ob0)#ob0.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-        else:
-            me0 = ob0.data.copy()
-
-        # generate new bmesh
-        bm = bmesh.new()
-        bm.from_mesh(me0)
-
-        # store weight values
-        weight = []
-        ob = bpy.data.objects.new("temp", me0)
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-        weight = get_weight_numpy(ob.vertex_groups[vertex_group_name], len(me0.vertices))
-
-        me0, bm, weight = contour_bmesh(me0, bm, weight, iso_val)
-
-        # Mask geometry
-        mask = weight >= iso_val
-        weight = weight[mask]
-        mask = np.logical_not(mask)
-        delete_verts = np.array(bm.verts)[mask]
-        #for v in delete_verts: bm.verts.remove(v)
-
-        # Create mesh and object
-        name = ob0.name + '_ContourMask_{:.3f}'.format(iso_val)
-        me = bpy.data.meshes.new(name)
-        bm.to_mesh(me)
-        ob = bpy.data.objects.new(name, me)
-
-        # Link object to scene and make active
-        scn = bpy.context.scene
-        bpy.context.collection.objects.link(ob)
-        bpy.context.view_layer.objects.active = ob
-        ob.select_set(True)
-        ob0.select_set(False)
-
-        # generate new vertex group
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-
-        if iso_val != 1: mult = 1/(1-iso_val)
-        else: mult = 1
-        for id in range(len(weight)):
-            if self.normalize_weight: w = (weight[id]-iso_val)*mult
-            else: w = weight[id]
-            ob.vertex_groups[vertex_group_name].add([id], w, 'REPLACE')
-        ob.vertex_groups.active_index = group_id
-
-        # align new object
-        ob.matrix_world = ob0.matrix_world
-
-        # Add Solidify
-        if self.bool_solidify and True:
-            ob.modifiers.new(type='SOLIDIFY', name='Solidify')
-            ob.modifiers['Solidify'].thickness = 0.05
-            ob.modifiers['Solidify'].offset = 0
-            ob.modifiers['Solidify'].vertex_group = vertex_group_name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        print("Contour Mask time: " + str(timeit.default_timer() - start_time) + " sec")
-
-        bpy.data.meshes.remove(me0)
-
-        return {'FINISHED'}
-
-
-class weight_contour_curves(bpy.types.Operator):
-    bl_idname = "object.weight_contour_curves"
-    bl_label = "Contour Curves"
-    bl_description = ("")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    use_modifiers : bpy.props.BoolProperty(
-        name="Use Modifiers", default=True,
-        description="Apply all the modifiers")
-
-    min_iso : bpy.props.FloatProperty(
-        name="Min Value", default=0., soft_min=0, soft_max=1,
-        description="Minimum weight value")
-    max_iso : bpy.props.FloatProperty(
-        name="Max Value", default=1, soft_min=0, soft_max=1,
-        description="Maximum weight value")
-    n_curves : bpy.props.IntProperty(
-        name="Curves", default=3, soft_min=1, soft_max=10,
-        description="Number of Contour Curves")
-
-    min_rad : bpy.props.FloatProperty(
-        name="Min Radius", default=1, soft_min=0, soft_max=1,
-        description="Change radius according to Iso Value")
-    max_rad : bpy.props.FloatProperty(
-        name="Max Radius", default=1, soft_min=0, soft_max=1,
-        description="Change radius according to Iso Value")
-
-    @classmethod
-    def poll(cls, context):
-        ob = context.object
-        return len(ob.vertex_groups) > 0 or ob.type == 'CURVE'
-
-    def invoke(self, context, event):
-        return context.window_manager.invoke_props_dialog(self, width=350)
-
-    def execute(self, context):
-        start_time = timeit.default_timer()
-        try:
-            check = bpy.context.object.vertex_groups[0]
-        except:
-            self.report({'ERROR'}, "The object doesn't have Vertex Groups")
-            return {'CANCELLED'}
-        ob0 = bpy.context.object
-
-        group_id = ob0.vertex_groups.active_index
-        vertex_group_name = ob0.vertex_groups[group_id].name
-
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.mesh.select_all(action='SELECT')
-        bpy.ops.object.mode_set(mode='OBJECT')
-        if self.use_modifiers:
-            me0 = simple_to_mesh(ob0) #ob0.to_mesh(preserve_all_data_layers=True, depsgraph=bpy.context.evaluated_depsgraph_get()).copy()
-        else:
-            me0 = ob0.data.copy()
-
-        # generate new bmesh
-        bm = bmesh.new()
-        bm.from_mesh(me0)
-        bm.verts.ensure_lookup_table()
-        bm.edges.ensure_lookup_table()
-        bm.faces.ensure_lookup_table()
-
-        # store weight values
-        weight = []
-        ob = bpy.data.objects.new("temp", me0)
-        for g in ob0.vertex_groups:
-            ob.vertex_groups.new(name=g.name)
-        weight = get_weight_numpy(ob.vertex_groups[vertex_group_name], len(bm.verts))
-
-        #filtered_edges = bm.edges
-        total_verts = np.zeros((0,3))
-        total_segments = []
-        radius = []
-
-        # start iterate contours levels
-        vertices = get_vertices_numpy(me0)
-        filtered_edges = get_edges_id_numpy(me0)
-
-        faces_weight = [np.array([weight[v] for v in p.vertices]) for p in me0.polygons]
-        fw_min = np.array([np.min(fw) for fw in faces_weight])
-        fw_max = np.array([np.max(fw) for fw in faces_weight])
-
-        bm_faces = np.array(bm.faces)
-
-        for c in range(self.n_curves):
-            min_iso = min(self.min_iso, self.max_iso)
-            max_iso = max(self.min_iso, self.max_iso)
-            try:
-                iso_val = c*(max_iso-min_iso)/(self.n_curves-1)+min_iso
-                if iso_val < 0: iso_val = (min_iso + max_iso)/2
-            except:
-                iso_val = (min_iso + max_iso)/2
-
-            # remove passed faces
-            bool_mask = iso_val < fw_max
-            bm_faces = bm_faces[bool_mask]
-            fw_min = fw_min[bool_mask]
-            fw_max = fw_max[bool_mask]
-
-            # mask faces
-            bool_mask = fw_min < iso_val
-            faces_mask = bm_faces[bool_mask]
-
-            n_verts = len(bm.verts)
-            count = len(total_verts)
-
-            # vertices indexes
-            id0 = filtered_edges[:,0]
-            id1 = filtered_edges[:,1]
-            # vertices weight
-            w0 = weight[id0]
-            w1 = weight[id1]
-            # weight condition
-            bool_w0 = w0 < iso_val
-            bool_w1 = w1 < iso_val
-
-            # mask all edges that have one weight value below the iso value
-            mask_new_verts = np.logical_xor(bool_w0, bool_w1)
-
-            id0 = id0[mask_new_verts]
-            id1 = id1[mask_new_verts]
-            # filter arrays
-            v0 = vertices[id0]
-            v1 = vertices[id1]
-            w0 = w0[mask_new_verts]
-            w1 = w1[mask_new_verts]
-            param = np.expand_dims((iso_val-w0)/(w1-w0),axis=1)
-            verts = v0 + (v1-v0)*param
-
-            # indexes of edges with new vertices
-            edges_index = filtered_edges[mask_new_verts][:,2]
-            edges_id = {}
-            for i, id in enumerate(edges_index): edges_id[id] = i+len(total_verts)
-
-            # remove all edges completely below the iso value
-            mask_edges = np.logical_not(np.logical_and(bool_w0, bool_w1))
-            filtered_edges = filtered_edges[mask_edges]
-            if len(verts) == 0: continue
-
-            # finding segments
-            segments = []
-            for f in faces_mask:
-                seg = []
-                for e in f.edges:
-                    try:
-                        seg.append(edges_id[e.index])
-                        if len(seg) == 2:
-                            segments.append(seg)
-                            seg = []
-                    except: pass
-
-
-            #curves_points_indexes = find_curves(segments)
-            total_segments = total_segments + segments
-            total_verts = np.concatenate((total_verts,verts))
-
-            if self.min_rad != self.max_rad:
-                try:
-                    iso_rad = c*(self.max_rad-self.min_rad)/(self.n_curves-1)+self.min_rad
-                    if iso_rad < 0: iso_rad = (self.min_rad + self.max_rad)/2
-                except:
-                    iso_rad = (self.min_rad + self.max_rad)/2
-                radius = radius + [iso_rad]*len(verts)
-        print("Contour Curves, computing time: " + str(timeit.default_timer() - start_time) + " sec")
-        bm = bmesh.new()
-        # adding new vertices
-        for v in total_verts: bm.verts.new(v)
-        bm.verts.ensure_lookup_table()
-
-        # adding new edges
-        for s in total_segments:
-            try:
-                pts = [bm.verts[i] for i in s]
-                bm.edges.new(pts)
-            except: pass
-
-
-        try:
-            name = ob0.name + '_ContourCurves'
-            me = bpy.data.meshes.new(name)
-            bm.to_mesh(me)
-            ob = bpy.data.objects.new(name, me)
-            # Link object to scene and make active
-            scn = bpy.context.scene
-            bpy.context.collection.objects.link(ob)
-            bpy.context.view_layer.objects.active = ob
-            ob.select_set(True)
-            ob0.select_set(False)
-
-            print("Contour Curves, bmesh time: " + str(timeit.default_timer() - start_time) + " sec")
-            bpy.ops.object.convert(target='CURVE')
-            ob = context.object
-            if not (self.min_rad == 0 and self.max_rad == 0):
-                if self.min_rad != self.max_rad:
-                    count = 0
-                    for s in ob.data.splines:
-                        for p in s.points:
-                            p.radius = radius[count]
-                            count += 1
-                else:
-                    for s in ob.data.splines:
-                        for p in s.points:
-                            p.radius = self.min_rad
-            ob.data.bevel_depth = 0.01
-            ob.data.fill_mode = 'FULL'
-            ob.data.bevel_resolution = 3
-        except:
-            self.report({'ERROR'}, "There are no values in the chosen range")
-            return {'CANCELLED'}
-
-        # align new object
-        ob.matrix_world = ob0.matrix_world
-        print("Contour Curves time: " + str(timeit.default_timer() - start_time) + " sec")
-
-        bpy.data.meshes.remove(me0)
-        bpy.data.meshes.remove(me)
-
-        return {'FINISHED'}
-
-class tissue_weight_contour_curves_pattern(bpy.types.Operator):
-    bl_idname = "object.tissue_weight_contour_curves_pattern"
-    bl_label = "Contour Curves Pattern"
-    bl_description = ("")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    use_modifiers : bpy.props.BoolProperty(
-        name="Use Modifiers", default=True,
-        description="Apply all the modifiers")
-
-    min_iso : bpy.props.FloatProperty(
-        name="Min Value", default=0., soft_min=0, soft_max=1,
-        description="Minimum weight value")
-    max_iso : bpy.props.FloatProperty(
-        name="Max Value", default=1, soft_min=0, soft_max=1,
-        description="Maximum weight value")
-    n_curves : bpy.props.IntProperty(
-        name="Curves", default=10, soft_min=1, soft_max=100,
-        description="Number of Contour Curves")
-    min_rad = 1
-    max_rad = 1
-
-    in_displace : bpy.props.FloatProperty(
-        name="In Displace", default=0, soft_min=0, soft_max=10,
-        description="Pattern displace strength")
-    out_displace : bpy.props.FloatProperty(
-        name="Out Displace", default=2, soft_min=0, soft_max=10,
-        description="Pattern displace strength")
-
-    in_steps : bpy.props.IntProperty(
-        name="In Steps", default=1, min=0, soft_max=10,
-        description="Number of layers to move inwards")
-    out_steps : bpy.props.IntProperty(
-        name="Out Steps", default=1, min=0, soft_max=10,
-        description="Number of layers to move outwards")
-    limit_z : bpy.props.BoolProperty(
-        name="Limit Z", default=False,
-        description="Limit Pattern in Z")
-
-    merge : bpy.props.BoolProperty(
-        name="Merge Vertices", default=True,
-        description="Merge points")
-    merge_thres : bpy.props.FloatProperty(
-        name="Merge Threshold", default=0.01, min=0, soft_max=1,
-        description="Minimum Curve Radius")
-
-    bevel_depth : bpy.props.FloatProperty(
-        name="Bevel Depth", default=0, min=0, soft_max=1,
-        description="")
-    remove_open_curves : bpy.props.BoolProperty(
-        name="Remove Open Curves", default=False,
-        description="Remove Open Curves")
-
-    vertex_group_pattern : bpy.props.StringProperty(
-        name="Pattern", default='',
-        description="Vertex Group used for pattern displace")
-
-    object_name : bpy.props.StringProperty(
-        name="Active Object", default='',
-        description="")
-
-    try: vg_name = bpy.context.object.vertex_groups.active.name
-    except: vg_name = ''
-
-    vertex_group_contour : bpy.props.StringProperty(
-        name="Contour", default=vg_name,
-        description="Vertex Group used for contouring")
-    clean_distance : bpy.props.FloatProperty(
-        name="Clean Distance", default=0, min=0, soft_max=10,
-        description="Remove short segments")
-
-
-    @classmethod
-    def poll(cls, context):
-        ob = context.object
-        return len(ob.vertex_groups) > 0 or ob.type == 'CURVE'
-
-    def invoke(self, context, event):
-        return context.window_manager.invoke_props_dialog(self, width=350)
-
-    def draw(self, context):
-        if not context.object.type == 'CURVE':
-            self.object_name = context.object.name
-        ob = bpy.data.objects[self.object_name]
-        if self.vertex_group_contour not in [vg.name for vg in ob.vertex_groups]:
-            self.vertex_group_contour = ob.vertex_groups.active.name
-        layout = self.layout
-        col = layout.column(align=True)
-        col.prop(self, "use_modifiers")
-        col.prop_search(self, 'vertex_group_contour', ob, "vertex_groups", text='Contour')
-        row = col.row(align=True)
-        row.prop(self,'min_iso')
-        row.prop(self,'max_iso')
-        col.prop(self,'n_curves')
-        col.separator()
-        col.prop_search(self, 'vertex_group_pattern', ob, "vertex_groups", text='Pattern')
-        if self.vertex_group_pattern != '':
-            row = col.row(align=True)
-            row.prop(self,'in_steps')
-            row.prop(self,'out_steps')
-            row = col.row(align=True)
-            row.prop(self,'in_displace')
-            row.prop(self,'out_displace')
-            col.prop(self,'limit_z')
-        col.separator()
-        col.label(text='Curves:')
-        col.prop(self,'bevel_depth')
-        col.prop(self,'clean_distance')
-        col.prop(self,'remove_open_curves')
-
-
-    def execute(self, context):
-        start_time = timeit.default_timer()
-        try:
-            check = context.object.vertex_groups[0]
-        except:
-            self.report({'ERROR'}, "The object doesn't have Vertex Groups")
-            return {'CANCELLED'}
-        ob0 = bpy.data.objects[self.object_name]
-
-        dg = bpy.context.evaluated_depsgraph_get()
-        ob = ob0.evaluated_get(dg)
-        me0 = ob.data
-
-        # generate new bmesh
-        bm = bmesh.new()
-        bm.from_mesh(me0)
-        n_verts = len(bm.verts)
-
-        # store weight values
-        try:
-            weight = get_weight_numpy(ob.vertex_groups[self.vertex_group_contour], len(me0.vertices))
-        except:
-            self.report({'ERROR'}, "Please select a Vertex Group for contouring")
-            return {'CANCELLED'}
-
-        try:
-            pattern_weight = get_weight_numpy(ob.vertex_groups[self.vertex_group_pattern], len(me0.vertices))
-        except:
-            self.report({'WARNING'}, "There is no Vertex Group assigned to the pattern displace")
-            pattern_weight = np.zeros(len(me0.vertices))
-
-        #filtered_edges = bm.edges
-        total_verts = np.zeros((0,3))
-        total_segments = []# np.array([])
-        radius = []
-
-        # start iterate contours levels
-        vertices, normals = get_vertices_and_normals_numpy(me0)
-        filtered_edges = get_edges_id_numpy(me0)
-
-        faces_weight = [np.array([weight[v] for v in p.vertices]) for p in me0.polygons]
-        fw_min = np.array([np.min(fw) for fw in faces_weight])
-        fw_max = np.array([np.max(fw) for fw in faces_weight])
-
-        bm_faces = np.array(bm.faces)
-        #bm_faces = np.array([np.array([e.index for e in f.edges]) for f in bm.faces ])
-        #all_edges_by_faces = []
-        #for f in bm.faces:
-        #    all_edges_by_faces += [e.index for e in f.edges]
-        #all_edges_by_faces = np.array(all_edges_by_faces)
-
-        print("Contour Curves, data loaded: " + str(timeit.default_timer() - start_time) + " sec")
-        step_time = timeit.default_timer()
-        for c in range(self.n_curves):
-            min_iso = min(self.min_iso, self.max_iso)
-            max_iso = max(self.min_iso, self.max_iso)
-            try:
-                iso_val = c*(max_iso-min_iso)/(self.n_curves-1)+min_iso
-                if iso_val < 0: iso_val = (min_iso + max_iso)/2
-            except:
-                iso_val = (min_iso + max_iso)/2
-
-            # remove passed faces
-            bool_mask = iso_val < fw_max
-            bm_faces = bm_faces[bool_mask]
-            fw_min = fw_min[bool_mask]
-            fw_max = fw_max[bool_mask]
-
-            # mask faces
-            bool_mask = fw_min < iso_val
-            faces_mask = bm_faces[bool_mask]
-
-            count = len(total_verts)
-
-            new_filtered_edges, edges_index, verts = contour_edges_pattern(self, c, len(total_verts), iso_val, vertices, normals, filtered_edges, weight, pattern_weight)
-            #new_filtered_edges, edges_index, verts = contour_edges_pattern(self.in_steps, self.out_steps, self.in_displace, self.out_displace, self.limit_z, c, iso_val, vertices, normals, filtered_edges, weight, pattern_weight)
-            if verts[0,0] == None: continue
-            else: filtered_edges = new_filtered_edges
-            edges_id = {}
-            for i, id in enumerate(edges_index): edges_id[id] = i + count
-
-            if len(verts) == 0: continue
-
-            # finding segments
-            segments = []
-            for f in faces_mask:
-                seg = []
-                for e in f.edges:
-                    try:
-                        #seg.append(new_ids[np.where(edges_index == e.index)[0][0]])
-                        seg.append(edges_id[e.index])
-                        if len(seg) == 2:
-                            segments.append(seg)
-                            seg = []
-                    except: pass
-
-            #segments = np.array([])
-            #keys = np.array(edges_id.keys())
-            #verts_id = np.arange(len(edges_index))+len(total_verts)
-
-            #seg_mask = np.in1d(edges_index, all_edges_by_faces)
-            #segments = verts_id[seg_mask]
-
-            #try: segments = segments.reshape((len(segments)//2,2))
-            #except: continue
-            '''
-            for f in faces_mask:
-                #print(f)
-                #print(edges_index)
-                #print(verts_id)
-                seg_mask = np.in1d(edges_index,f)
-                #if not seg_mask.any(): continue
-                seg = verts_id[seg_mask]
-                seg = seg.reshape((len(seg)//2,2))
-                if len(seg)==0: continue
-                segments = seg if len(segments)==0 else np.concatenate((segments ,seg))
-            '''
-            #if len(segments)>0:
-            #    total_segments = segments if len(total_segments)==0 else np.concatenate((total_segments, segments))
-
-
-            total_segments = total_segments + segments
-            total_verts = np.concatenate((total_verts,verts))
-
-            if self.min_rad != self.max_rad:
-                try:
-                    iso_rad = c*(self.max_rad-self.min_rad)/(self.n_curves-1)+self.min_rad
-                    if iso_rad < 0: iso_rad = (self.min_rad + self.max_rad)/2
-                except:
-                    iso_rad = (self.min_rad + self.max_rad)/2
-                radius = radius + [iso_rad]*len(verts)
-        print("Contour Curves, points computing: " + str(timeit.default_timer() - step_time) + " sec")
-        step_time = timeit.default_timer()
-        if len(total_segments) > 0:
-            step_time = timeit.default_timer()
-            ordered_points = find_curves(total_segments, len(total_verts))
-            print("Contour Curves, point ordered in: " + str(timeit.default_timer() - step_time) + " sec")
-            step_time = timeit.default_timer()
-            crv = curve_from_pydata(total_verts,ordered_points,ob0.name + '_ContourCurves', self.remove_open_curves, merge_distance=self.clean_distance)
-            context.view_layer.objects.active = crv
-            crv.data.bevel_depth = self.bevel_depth
-
-            crv.select_set(True)
-            ob0.select_set(False)
-            crv.matrix_world = ob0.matrix_world
-            print("Contour Curves, curves created in: " + str(timeit.default_timer() - step_time) + " sec")
-        else:
-            self.report({'ERROR'}, "There are no values in the chosen range")
-            return {'CANCELLED'}
-        print("Contour Curves, total time: " + str(timeit.default_timer() - start_time) + " sec")
-        return {'FINISHED'}
-
-class vertex_colors_to_vertex_groups(bpy.types.Operator):
-    bl_idname = "object.vertex_colors_to_vertex_groups"
-    bl_label = "Vertex Color"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Convert the active Vertex Color into a Vertex Group.")
-
-    red : bpy.props.BoolProperty(
-        name="red channel", default=False, description="convert red channel")
-    green : bpy.props.BoolProperty(
-        name="green channel", default=False,
-        description="convert green channel")
-    blue : bpy.props.BoolProperty(
-        name="blue channel", default=False, description="convert blue channel")
-    value : bpy.props.BoolProperty(
-        name="value channel", default=True, description="convert value channel")
-    invert : bpy.props.BoolProperty(
-         name="invert", default=False, description="invert all color channels")
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.data.vertex_colors) > 0
-
-    def execute(self, context):
-        obj = bpy.context.active_object
-        id = len(obj.vertex_groups)
-        id_red = id
-        id_green = id
-        id_blue = id
-        id_value = id
-
-        boolCol = len(obj.data.vertex_colors)
-        if(boolCol): col_name = obj.data.vertex_colors.active.name
-        bpy.ops.object.mode_set(mode='EDIT')
-        bpy.ops.mesh.select_all(action='SELECT')
-
-        if(self.red and boolCol):
-            bpy.ops.object.vertex_group_add()
-            bpy.ops.object.vertex_group_assign()
-            id_red = id
-            obj.vertex_groups[id_red].name = col_name + '_red'
-            id+=1
-        if(self.green and boolCol):
-            bpy.ops.object.vertex_group_add()
-            bpy.ops.object.vertex_group_assign()
-            id_green = id
-            obj.vertex_groups[id_green].name = col_name + '_green'
-            id+=1
-        if(self.blue and boolCol):
-            bpy.ops.object.vertex_group_add()
-            bpy.ops.object.vertex_group_assign()
-            id_blue = id
-            obj.vertex_groups[id_blue].name = col_name + '_blue'
-            id+=1
-        if(self.value and boolCol):
-            bpy.ops.object.vertex_group_add()
-            bpy.ops.object.vertex_group_assign()
-            id_value = id
-            obj.vertex_groups[id_value].name = col_name + '_value'
-            id+=1
-
-        mult = 1
-        if(self.invert): mult = -1
-        bpy.ops.object.mode_set(mode='OBJECT')
-        sub_red = 1 + self.value + self.blue + self.green
-        sub_green = 1 + self.value + self.blue
-        sub_blue = 1 + self.value
-        sub_value = 1
-
-        id = len(obj.vertex_groups)
-        if(id_red <= id and id_green <= id and id_blue <= id and id_value <= \
-                id and boolCol):
-            v_colors = obj.data.vertex_colors.active.data
-            i = 0
-            for f in obj.data.polygons:
-                for v in f.vertices:
-                    gr = obj.data.vertices[v].groups
-                    if(self.red): gr[min(len(gr)-sub_red, id_red)].weight = \
-                        self.invert + mult * v_colors[i].color[0]
-                    if(self.green): gr[min(len(gr)-sub_green, id_green)].weight\
-                        = self.invert + mult * v_colors[i].color[1]
-                    if(self.blue): gr[min(len(gr)-sub_blue, id_blue)].weight = \
-                        self.invert + mult * v_colors[i].color[2]
-                    if(self.value):
-                        r = v_colors[i].color[0]
-                        g = v_colors[i].color[1]
-                        b = v_colors[i].color[2]
-                        gr[min(len(gr)-sub_value, id_value)].weight\
-                        = self.invert + mult * (0.2126*r + 0.7152*g + 0.0722*b)
-                    i+=1
-            bpy.ops.paint.weight_paint_toggle()
-        return {'FINISHED'}
-
-class vertex_group_to_vertex_colors(bpy.types.Operator):
-    bl_idname = "object.vertex_group_to_vertex_colors"
-    bl_label = "Vertex Group"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Convert the active Vertex Group into a Vertex Color.")
-
-    channel : bpy.props.EnumProperty(
-        items=[('Blue', 'Blue Channel', 'Convert to Blue Channel'),
-               ('Green', 'Green Channel', 'Convert to Green Channel'),
-               ('Red', 'Red Channel', 'Convert to Red Channel'),
-               ('Value', 'Value Channel', 'Convert to Grayscale'),
-               ('False Colors', 'False Colors', 'Convert to False Colors')],
-        name="Convert to", description="Choose how to convert vertex group",
-        default="Value", options={'LIBRARY_EDITABLE'})
-
-    invert : bpy.props.BoolProperty(
-        name="invert", default=False, description="invert color channel")
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def execute(self, context):
-        obj = bpy.context.active_object
-        group_id = obj.vertex_groups.active_index
-        if (group_id == -1):
-            return {'FINISHED'}
-
-        bpy.ops.object.mode_set(mode='OBJECT')
-        group_name = obj.vertex_groups[group_id].name
-        bpy.ops.mesh.vertex_color_add()
-        colors_id = obj.data.vertex_colors.active_index
-
-        colors_name = group_name
-        if(self.channel == 'False Colors'): colors_name += "_false_colors"
-        elif(self.channel == 'Value'):  colors_name += "_value"
-        elif(self.channel == 'Red'):  colors_name += "_red"
-        elif(self.channel == 'Green'):  colors_name += "_green"
-        elif(self.channel == 'Blue'):  colors_name += "_blue"
-        bpy.context.object.data.vertex_colors[colors_id].name = colors_name
-
-        v_colors = obj.data.vertex_colors.active.data
-
-        mult = 1
-        if(self.invert): mult = -1
-
-        i = 0
-        for f in obj.data.polygons:
-            for v in f.vertices:
-                gr = obj.data.vertices[v].groups
-
-                if(self.channel == 'False Colors'): v_colors[i].color = (0,0,0.5,1)
-                else: v_colors[i].color = (0,0,0,1)
-
-                for g in gr:
-                    if g.group == group_id:
-                        w = g.weight
-                        if(self.channel == 'False Colors'):
-                            mult = 0.6+0.4*w
-                            if w < 0.25:
-                                v_colors[i].color = (0, w*4*mult, 1*mult,1)
-                            elif w < 0.5:
-                                v_colors[i].color = (0, 1*mult, (1-(w-0.25)*4)*mult,1)
-                            elif w < 0.75:
-                                v_colors[i].color = ((w-0.5)*4*mult,1*mult,0,1)
-                            else:
-                                v_colors[i].color = (1*mult,(1-(w-0.75)*4)*mult,0,1)
-                        elif(self.channel == 'Value'):
-                            v_colors[i].color = (
-                                self.invert + mult * w,
-                                self.invert + mult * w,
-                                self.invert + mult * w,
-                                1)
-                        elif(self.channel == 'Red'):
-                            v_colors[i].color = (
-                                self.invert + mult * w,0,0,1)
-                        elif(self.channel == 'Green'):
-                            v_colors[i].color = (
-                                0, self.invert + mult * w,0,1)
-                        elif(self.channel == 'Blue'):
-                            v_colors[i].color = (
-                                0,0, self.invert + mult * w,1)
-                i+=1
-        bpy.ops.paint.vertex_paint_toggle()
-        bpy.context.object.data.vertex_colors[colors_id].active_render = True
-        return {'FINISHED'}
-
-class curvature_to_vertex_groups(bpy.types.Operator):
-    bl_idname = "object.curvature_to_vertex_groups"
-    bl_label = "Curvature"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Generate a Vertex Group based on the curvature of the"
-                      "mesh. Is based on Dirty Vertex Color.")
-
-    invert : bpy.props.BoolProperty(
-        name="invert", default=False, description="invert values")
-
-    blur_strength : bpy.props.FloatProperty(
-      name="Blur Strength", default=1, min=0.001,
-      max=1, description="Blur strength per iteration")
-
-    blur_iterations : bpy.props.IntProperty(
-      name="Blur Iterations", default=1, min=0,
-      max=40, description="Number of times to blur the values")
-
-    min_angle : bpy.props.FloatProperty(
-      name="Min Angle", default=0, min=0,
-      max=pi/2, subtype='ANGLE', description="Minimum angle")
-
-    max_angle : bpy.props.FloatProperty(
-      name="Max Angle", default=pi, min=pi/2,
-      max=pi, subtype='ANGLE', description="Maximum angle")
-
-    invert : bpy.props.BoolProperty(
-        name="Invert", default=False,
-        description="Invert the curvature map")
-
-    def execute(self, context):
-        bpy.ops.object.mode_set(mode='OBJECT')
-        bpy.ops.mesh.vertex_color_add()
-        vertex_colors = bpy.context.active_object.data.vertex_colors
-        vertex_colors[-1].active = True
-        vertex_colors[-1].active_render = True
-        vertex_colors[-1].name = "Curvature"
-        for c in vertex_colors[-1].data: c.color = (1,1,1,1)
-        bpy.ops.object.mode_set(mode='VERTEX_PAINT')
-        bpy.ops.paint.vertex_color_dirt(
-            blur_strength=self.blur_strength,
-            blur_iterations=self.blur_iterations, clean_angle=self.max_angle,
-            dirt_angle=self.min_angle)
-        bpy.ops.object.vertex_colors_to_vertex_groups(invert=self.invert)
-        bpy.ops.mesh.vertex_color_remove()
-        return {'FINISHED'}
-
-
-class face_area_to_vertex_groups(bpy.types.Operator):
-    bl_idname = "object.face_area_to_vertex_groups"
-    bl_label = "Area"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Generate a Vertex Group based on the area of individual"
-                      "faces.")
-
-    invert : bpy.props.BoolProperty(
-        name="invert", default=False, description="invert values")
-    bounds : bpy.props.EnumProperty(
-        items=(('MANUAL', "Manual Bounds", ""),
-            ('AUTOMATIC', "Automatic Bounds", "")),
-        default='AUTOMATIC', name="Bounds")
-
-    min_area : bpy.props.FloatProperty(
-        name="Min", default=0.01, soft_min=0, soft_max=1,
-        description="Faces with 0 weight")
-
-    max_area : bpy.props.FloatProperty(
-        name="Max", default=0.1, soft_min=0, soft_max=1,
-        description="Faces with 1 weight")
-
-    def draw(self, context):
-        layout = self.layout
-        layout.label(text="Bounds")
-        layout.prop(self, "bounds", text="")
-        if self.bounds == 'MANUAL':
-            layout.prop(self, "min_area")
-            layout.prop(self, "max_area")
-
-    def execute(self, context):
-        try: ob = context.object
-        except:
-            self.report({'ERROR'}, "Please select an Object")
-            return {'CANCELLED'}
-        ob.vertex_groups.new(name="Faces Area")
-
-        areas = [[] for v in ob.data.vertices]
-
-        for p in ob.data.polygons:
-            for v in p.vertices:
-                areas[v].append(p.area)
-
-        for i in range(len(areas)):
-            areas[i] = mean(areas[i])
-        if self.bounds == 'MANUAL':
-            min_area = self.min_area
-            max_area = self.max_area
-        elif self.bounds == 'AUTOMATIC':
-            min_area = min(areas)
-            max_area = max(areas)
-        elif self.bounds == 'COMPRESSION':
-            min_area = 1
-            max_area = min(areas)
-        elif self.bounds == 'TENSION':
-            min_area = 1
-            max_area = max(areas)
-        delta_area = max_area - min_area
-        if delta_area == 0:
-            delta_area = 0.0001
-            if self.bounds == 'MANUAL':
-                delta_area = 0.0001
-            else:
-                self.report({'ERROR'}, "The faces have the same areas")
-                #return {'CANCELLED'}
-        for i in range(len(areas)):
-            weight = (areas[i] - min_area)/delta_area
-            ob.vertex_groups[-1].add([i], weight, 'REPLACE')
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-
-class harmonic_weight(bpy.types.Operator):
-    bl_idname = "object.harmonic_weight"
-    bl_label = "Harmonic"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Create an harmonic variation of the active Vertex Group")
-
-    freq : bpy.props.FloatProperty(
-        name="Frequency", default=20, soft_min=0,
-        soft_max=100, description="Wave frequency")
-
-    amp : bpy.props.FloatProperty(
-        name="Amplitude", default=1, soft_min=0,
-        soft_max=10, description="Wave amplitude")
-
-    midlevel : bpy.props.FloatProperty(
-        name="Midlevel", default=0, min=-1,
-        max=1, description="Midlevel")
-
-    add : bpy.props.FloatProperty(
-        name="Add", default=0, min=-1,
-        max=1, description="Add to the Weight")
-
-    mult : bpy.props.FloatProperty(
-        name="Multiply", default=0, min=0,
-        max=1, description="Multiply for he Weight")
-
-    @classmethod
-    def poll(cls, context):
-        return len(context.object.vertex_groups) > 0
-
-    def execute(self, context):
-        ob = bpy.context.active_object
-        if len(ob.vertex_groups) > 0:
-            group_id = ob.vertex_groups.active_index
-            ob.vertex_groups.new(name="Harmonic")
-            for i in range(len(ob.data.vertices)):
-                try: val = ob.vertex_groups[group_id].weight(i)
-                except: val = 0
-                weight = self.amp*(sin(val*self.freq) - self.midlevel)/2 + 0.5 + self.add*val*(1-(1-val)*self.mult)
-                ob.vertex_groups[-1].add([i], weight, 'REPLACE')
-            ob.data.update()
-        else:
-            self.report({'ERROR'}, "Active object doesn't have vertex groups")
-            return {'CANCELLED'}
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-
-class tissue_weight_distance(bpy.types.Operator):
-    bl_idname = "object.tissue_weight_distance"
-    bl_label = "Weight Distance"
-    bl_options = {'REGISTER', 'UNDO'}
-    bl_description = ("Create a weight map according to the distance from the "
-                    "selected vertices along the mesh surface")
-
-    def fill_neighbors(self,verts,weight):
-        neigh = {}
-        for v0 in verts:
-            for f in v0.link_faces:
-                for v1 in f.verts:
-                    dist = weight[v0.index] + (v0.co-v1.co).length
-                    w1 = weight[v1.index]
-                    if w1 == None or w1 > dist:
-                        weight[v1.index] = dist
-                        neigh[v1] = 0
-        if len(neigh) == 0: return weight
-        else: return self.fill_neighbors(neigh.keys(), weight)
-
-    def execute(self, context):
-        ob = context.object
-        old_mode = ob.mode
-        if old_mode != 'OBJECT':
-            bpy.ops.object.mode_set(mode='OBJECT')
-
-        me = ob.data
-        bm = bmesh.new()
-        bm.from_mesh(me)
-        bm.verts.ensure_lookup_table()
-        bm.edges.ensure_lookup_table()
-        bm.faces.ensure_lookup_table()
-
-        # store weight values
-        weight = [None]*len(bm.verts)
-
-        selected = [v for v in bm.verts if v.select]
-        if len(selected) == 0:
-            bpy.ops.object.mode_set(mode=old_mode)
-            message = "Please, select one or more vertices"
-            self.report({'ERROR'}, message)
-            return {'CANCELLED'}
-        for v in selected: weight[v.index] = 0
-        weight = self.fill_neighbors(selected, weight)
-
-        weight = np.array(weight)
-        max_dist = np.max(weight)
-        if max_dist > 0:
-            weight /= max_dist
-
-        vg = ob.vertex_groups.new(name='Distance: {:.4f}'.format(max_dist))
-        for i, w in enumerate(weight):
-            if w == None: continue
-            vg.add([i], w, 'REPLACE')
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-        return {'FINISHED'}
-
-
-
-class TISSUE_PT_color(bpy.types.Panel):
-    bl_label = "Tissue Tools"
-    bl_category = "Tissue"
-    bl_space_type = "VIEW_3D"
-    bl_region_type = "UI"
-    #bl_options = {'DEFAULT_CLOSED'}
-    bl_context = "vertexpaint"
-
-    def draw(self, context):
-        layout = self.layout
-        col = layout.column(align=True)
-        col.operator("object.vertex_colors_to_vertex_groups",
-            icon="GROUP_VERTEX", text="Convert to Weight")
-
-class TISSUE_PT_weight(bpy.types.Panel):
-    bl_label = "Tissue Tools"
-    bl_category = "Tissue"
-    bl_space_type = "VIEW_3D"
-    bl_region_type = "UI"
-    #bl_options = {'DEFAULT_CLOSED'}
-    bl_context = "weightpaint"
-
-    def draw(self, context):
-        layout = self.layout
-        col = layout.column(align=True)
-        #if context.object.type == 'MESH' and context.mode == 'OBJECT':
-            #col.label(text="Transform:")
-            #col.separator()
-        #elif bpy.context.mode == 'PAINT_WEIGHT':
-        col.label(text="Weight Generate:")
-        #col.operator(
-        #    "object.vertex_colors_to_vertex_groups", icon="GROUP_VCOL")
-        col.operator("object.face_area_to_vertex_groups", icon="FACESEL")
-        col.operator("object.curvature_to_vertex_groups", icon="SMOOTHCURVE")
-        col.operator("object.tissue_weight_distance", icon="TRACKING")
-        try: col.operator("object.weight_formula", icon="CON_TRANSFORM")
-        except: col.operator("object.weight_formula")#, icon="CON_TRANSFORM")
-        #col.label(text="Weight Processing:")
-        col.separator()
-
-        # TO BE FIXED
-        #col.operator("object.weight_laplacian", icon="SMOOTHCURVE")
-
-        col.operator("object.harmonic_weight", icon="IPO_ELASTIC")
-        col.operator("object.vertex_group_to_vertex_colors", icon="GROUP_VCOL",
-            text="Convert to Colors")
-        col.separator()
-        col.label(text="Deformation Analysis:")
-        col.operator("object.edges_deformation", icon="DRIVER_DISTANCE")#FULLSCREEN_ENTER")
-        col.operator("object.edges_bending", icon="DRIVER_ROTATIONAL_DIFFERENCE")#"MOD_SIMPLEDEFORM")
-        col.separator()
-        col.label(text="Weight Contour:")
-        col.operator("object.weight_contour_curves", icon="MOD_CURVE")
-        col.operator("object.tissue_weight_contour_curves_pattern", icon="FORCE_TURBULENCE")
-        col.operator("object.weight_contour_displace", icon="MOD_DISPLACE")
-        col.operator("object.weight_contour_mask", icon="MOD_MASK")
-        col.separator()
-        col.label(text="Simulations:")
-        #col.operator("object.reaction_diffusion", icon="MOD_OCEAN")
-        col.operator("object.start_reaction_diffusion",
-                    icon="EXPERIMENTAL",
-                    text="Reaction-Diffusion")
-
-        #col.prop(context.object, "reaction_diffusion_run", icon="PLAY", text="Run Simulation")
-        ####col.prop(context.object, "reaction_diffusion_run")
-        #col.separator()
-        #col.label(text="Vertex Color from:")
-        #col.operator("object.vertex_group_to_vertex_colors", icon="GROUP_VERTEX")
-
-
-
-
-class start_reaction_diffusion(bpy.types.Operator):
-    bl_idname = "object.start_reaction_diffusion"
-    bl_label = "Start Reaction Diffusion"
-    bl_description = ("Run a Reaction-Diffusion based on existing Vertex Groups: A and B")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    run : bpy.props.BoolProperty(
-        name="Run Reaction-Diffusion", default=True, description="Compute a new iteration on frame changes")
-
-    time_steps : bpy.props.IntProperty(
-        name="Steps", default=10, min=0, soft_max=50,
-        description="Number of Steps")
-
-    dt : bpy.props.FloatProperty(
-        name="dt", default=1, min=0, soft_max=0.2,
-        description="Time Step")
-
-    diff_a : bpy.props.FloatProperty(
-        name="Diff A", default=0.18, min=0, soft_max=2,
-        description="Diffusion A")
-
-    diff_b : bpy.props.FloatProperty(
-        name="Diff B", default=0.09, min=0, soft_max=2,
-        description="Diffusion B")
-
-    f : bpy.props.FloatProperty(
-        name="f", default=0.055, min=0, soft_max=0.5, precision=4,
-        description="Feed Rate")
-
-    k : bpy.props.FloatProperty(
-        name="k", default=0.062, min=0, soft_max=0.5, precision=4,
-        description="Kill Rate")
-
-    @classmethod
-    def poll(cls, context):
-        return context.object.type == 'MESH'
-
-    def execute(self, context):
-        reaction_diffusion_add_handler(self, context)
-        set_animatable_fix_handler(self, context)
-
-        ob = context.object
-
-        ob.reaction_diffusion_settings.run = self.run
-        ob.reaction_diffusion_settings.dt = self.dt
-        ob.reaction_diffusion_settings.time_steps = self.time_steps
-        ob.reaction_diffusion_settings.f = self.f
-        ob.reaction_diffusion_settings.k = self.k
-        ob.reaction_diffusion_settings.diff_a = self.diff_a
-        ob.reaction_diffusion_settings.diff_b = self.diff_b
-
-
-        # check vertex group A
-        try:
-            vg = ob.vertex_groups['A']
-        except:
-            ob.vertex_groups.new(name='A')
-        # check vertex group B
-        try:
-            vg = ob.vertex_groups['B']
-        except:
-            ob.vertex_groups.new(name='B')
-
-        for v in ob.data.vertices:
-            ob.vertex_groups['A'].add([v.index], 1, 'REPLACE')
-            ob.vertex_groups['B'].add([v.index], 0, 'REPLACE')
-
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-
-        return {'FINISHED'}
-
-class reset_reaction_diffusion_weight(bpy.types.Operator):
-    bl_idname = "object.reset_reaction_diffusion_weight"
-    bl_label = "Reset Reaction Diffusion Weight"
-    bl_description = ("Set A and B weight to default values")
-    bl_options = {'REGISTER', 'UNDO'}
-
-    @classmethod
-    def poll(cls, context):
-        return context.object.type == 'MESH'
-
-    def execute(self, context):
-        reaction_diffusion_add_handler(self, context)
-        set_animatable_fix_handler(self, context)
-
-        ob = context.object
-
-        # check vertex group A
-        try:
-            vg = ob.vertex_groups['A']
-        except:
-            ob.vertex_groups.new(name='A')
-        # check vertex group B
-        try:
-            vg = ob.vertex_groups['B']
-        except:
-            ob.vertex_groups.new(name='B')
-
-        for v in ob.data.vertices:
-            ob.vertex_groups['A'].add([v.index], 1, 'REPLACE')
-            ob.vertex_groups['B'].add([v.index], 0, 'REPLACE')
-
-        ob.vertex_groups.update()
-        ob.data.update()
-        bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-
-        return {'FINISHED'}
-
-from bpy.app.handlers import persistent
-
-@persistent
-def reaction_diffusion_def_blur(scene):
-    for ob in scene.objects:
-        if ob.reaction_diffusion_settings.run:
-            #try:
-            me = ob.data
-            bm = bmesh.new()
-            bm.from_mesh(me)
-            bm.edges.ensure_lookup_table()
-
-            # store weight values
-            a = []
-            b = []
-            for v in me.vertices:
-                try:
-                    a.append(ob.vertex_groups["A"].weight(v.index))
-                except:
-                    a.append(0)
-                try:
-                    b.append(ob.vertex_groups["B"].weight(v.index))
-                except:
-                    b.append(0)
-
-            a = array(a)
-            b = array(b)
-            props = ob.reaction_diffusion_settings
-            dt = props.dt
-            time_steps = props.time_steps
-            f = props.f
-            k = props.k
-            diff_a = props.diff_a * props.diff_mult
-            diff_b = props.diff_b * props.diff_mult
-
-            n_verts = len(bm.verts)
-            #bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-            #ob.data.use_paint_mask_vertex = True
-
-            for i in range(time_steps):
-                ab2 = a*b**2
-                ob.vertex_groups.active = ob.vertex_groups['A']
-                bpy.ops.object.vertex_group_smooth(group_select_mode='ACTIVE', factor=diff_a)
-                ob.vertex_groups.active = ob.vertex_groups['B']
-                bpy.ops.object.vertex_group_smooth(group_select_mode='ACTIVE', factor=diff_b)
-
-                a = []
-                b = []
-                for v in me.vertices:
-                    a.append(ob.vertex_groups["A"].weight(v.index))
-                    b.append(ob.vertex_groups["B"].weight(v.index))
-                a = array(a)
-                b = array(b)
-
-                a += - (ab2 + f*(1-a))*dt
-                b += (ab2 - (k+f)*b)*dt
-
-            a = nan_to_num(a)
-            b = nan_to_num(b)
-
-            for i in range(n_verts):
-                ob.vertex_groups['A'].add([i], a[i], 'REPLACE')
-                ob.vertex_groups['B'].add([i], b[i], 'REPLACE')
-            ob.vertex_groups.update()
-            ob.data.update()
-            #bpy.ops.object.mode_set(mode='EDIT')
-            #bpy.ops.object.mode_set(mode='WEIGHT_PAINT
-            #bpy.ops.paint.weight_paint_toggle()
-            #bpy.ops.paint.weight_paint_toggle()
-
-            #bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-            #except:
-            #    pass
-
-def reaction_diffusion_def_(scene):
-    for ob in scene.objects:
-        if ob.reaction_diffusion_settings.run:
-            #try:
-            me = ob.data
-            bm = bmesh.new()
-            bm.from_mesh(me)
-            bm.edges.ensure_lookup_table()
-
-            # store weight values
-            a = []
-            b = []
-            for v in me.vertices:
-                try:
-                    a.append(ob.vertex_groups["A"].weight(v.index))
-                except:
-                    a.append(0)
-                try:
-                    b.append(ob.vertex_groups["B"].weight(v.index))
-                except:
-                    b.append(0)
-
-            a = array(a)
-            b = array(b)
-            props = ob.reaction_diffusion_settings
-            dt = props.dt
-            time_steps = props.time_steps
-            f = props.f
-            k = props.k
-            diff_a = props.diff_a * props.diff_mult
-            diff_b = props.diff_b * props.diff_mult
-
-            n_verts = len(bm.verts)
-            for i in range(time_steps):
-                lap_a = zeros((n_verts))#[0]*n_verts
-                lap_b = zeros((n_verts))#[0]*n_verts
-                if i == 0:
-                    lap_map = [[] for i in range(n_verts)]
-                    lap_mult = []
-                    for e in bm.edges:
-                        id0 = e.verts[0].index
-                        id1 = e.verts[1].index
-                        lap_map[id0].append(id1)
-                        lap_map[id1].append(id0)
-                    for id in range(n_verts):
-                         lap_mult.append(len(lap_map[id]))
-                    lap_mult = array(lap_mult)
-                    lap_map = array(lap_map)
-                for id in range(n_verts):
-                    map = lap_map[id]
-                    lap_a[id] = a[lap_map[id]].sum()
-                    lap_b[id] = b[lap_map[id]].sum()
-                lap_a -= a*lap_mult
-                lap_b -= b*lap_mult
-                ab2 = a*b**2
-
-                a += (diff_a*lap_a - ab2 + f*(1-a))*dt
-                b += (diff_b*lap_b + ab2 - (k+f)*b)*dt
-
-            a = nan_to_num(a)
-            b = nan_to_num(b)
-
-            for i in range(n_verts):
-                ob.vertex_groups['A'].add([i], a[i], 'REPLACE')
-                ob.vertex_groups['B'].add([i], b[i], 'REPLACE')
-            ob.vertex_groups.update()
-            ob.data.update()
-            #bpy.ops.object.mode_set(mode='EDIT')
-            #bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-            bpy.ops.paint.weight_paint_toggle()
-            bpy.ops.paint.weight_paint_toggle()
-
-            #bpy.ops.object.mode_set(mode='WEIGHT_PAINT')
-            #except:
-            #    pass
-
-def reaction_diffusion_def(scene):
-    for ob in scene.objects:
-        if ob.reaction_diffusion_settings.run:
-
-            start = time.time()
-
-            me = ob.data
-            n_edges = len(me.edges)
-            n_verts = len(me.vertices)
-
-            # store weight values
-            a = np.zeros(n_verts)
-            b = np.zeros(n_verts)
-            #a = thread_read_weight(a, ob.vertex_groups["A"])
-            #b = thread_read_weight(b, ob.vertex_groups["B"])
-            #a = read_weight(a, ob.vertex_groups["A"])
-            #b = read_weight(b, ob.vertex_groups["B"])
-
-            for i in range(n_verts):
-                try: a[i] = ob.vertex_groups["A"].weight(i)
-                except: pass
-                try: b[i] = ob.vertex_groups["B"].weight(i)
-                except: pass
-
-            props = ob.reaction_diffusion_settings
-            dt = props.dt
-            time_steps = props.time_steps
-            f = props.f
-            k = props.k
-            diff_a = props.diff_a * props.diff_mult
-            diff_b = props.diff_b * props.diff_mult
-
-            edge_verts = [0]*n_edges*2
-            me.edges.foreach_get("vertices", edge_verts)
-
-            timeElapsed = time.time() - start
-            print('RD - Preparation Time:',timeElapsed)
-            start = time.time()
-
-            try:
-                edge_verts = np.array(edge_verts)
-                a, b = numba_reaction_diffusion(n_verts, n_edges, edge_verts, a, b, diff_a, diff_b, f, k, dt, time_steps)
-                a = nan_to_num(a)
-                b = nan_to_num(b)
-            except:
-                edge_verts = np.array(edge_verts)
-                arr = np.arange(n_edges)*2
-                id0 = edge_verts[arr]     # first vertex indices for each edge
-                id1 = edge_verts[arr+1]   # second vertex indices for each edge
-                for i in range(time_steps):
-                    lap_a = np.zeros(n_verts)
-                    lap_b = np.zeros(n_verts)
-                    lap_a0 =  a[id1] -  a[id0]   # laplacian increment for first vertex of each edge
-                    lap_b0 =  b[id1] -  b[id0]   # laplacian increment for first vertex of each edge
-
-                    for i, j, la0, lb0 in np.nditer([id0,id1,lap_a0,lap_b0]):
-                        lap_a[i] += la0
-                        lap_b[i] += lb0
-                        lap_a[j] -= la0
-                        lap_b[j] -= lb0
-                    ab2 = a*b**2
-                    a += eval("(diff_a*lap_a - ab2 + f*(1-a))*dt")
-                    b += eval("(diff_b*lap_b + ab2 - (k+f)*b)*dt")
-                    #a += (diff_a*lap_a - ab2 + f*(1-a))*dt
-                    #b += (diff_b*lap_b + ab2 - (k+f)*b)*dt
-
-                    a = nan_to_num(a)
-                    b = nan_to_num(b)
-
-            timeElapsed = time.time() - start
-            print('RD - Simulation Time:',timeElapsed)
-            start = time.time()
-
-            for i in range(n_verts):
-                ob.vertex_groups['A'].add([i], a[i], 'REPLACE')
-                ob.vertex_groups['B'].add([i], b[i], 'REPLACE')
-
-            for ps in ob.particle_systems:
-                if ps.vertex_group_density == 'B' or ps.vertex_group_density == 'A':
-                    ps.invert_vertex_group_density = not ps.invert_vertex_group_density
-                    ps.invert_vertex_group_density = not ps.invert_vertex_group_density
-
-            timeElapsed = time.time() - start
-            print('RD - Closing Time:',timeElapsed)
-
-class TISSUE_PT_reaction_diffusion(Panel):
-    bl_space_type = 'PROPERTIES'
-    bl_region_type = 'WINDOW'
-    bl_context = "data"
-    bl_label = "Tissue - Reaction-Diffusion"
-    bl_options = {'DEFAULT_CLOSED'}
-
-    @classmethod
-    def poll(cls, context):
-        return 'A' and 'B' in context.object.vertex_groups
-
-    def draw(self, context):
-        reaction_diffusion_add_handler(self, context)
-
-        ob = context.object
-        props = ob.reaction_diffusion_settings
-        layout = self.layout
-        col = layout.column(align=True)
-        row = col.row(align=True)
-        if not ("A" and "B" in ob.vertex_groups):
-            row.operator("object.start_reaction_diffusion",
-                        icon="EXPERIMENTAL",
-                        text="Reaction-Diffusion")
-        else:
-            row.operator("object.start_reaction_diffusion",
-                        icon="EXPERIMENTAL",
-                        text="Reset Reaction-Diffusion")
-            row = col.row(align=True)
-            row.prop(props, "run", text="Run Reaction-Diffusion")
-            col = layout.column(align=True)
-            row = col.row(align=True)
-            row.prop(props, "time_steps")
-            row.prop(props, "dt")
-            col.separator()
-            row = col.row(align=True)
-            row.prop(props, "diff_a")
-            row.prop(props, "diff_b")
-            row = col.row(align=True)
-            row.prop(props, "diff_mult")
-            #col.separator()
-            row = col.row(align=True)
-            row.prop(props, "f")
-            row.prop(props, "k")
-
-if False:
-    @jit(["float64[:,:](int32, int32, float64, float64, boolean, int32, float64, float64[:,:], float64[:,:], int32[:,:], float64[:], float64[:])"]) #(nopython=True, parallel=True)
-    def contour_edges_pattern(in_steps, out_steps, in_displace, out_displace, limit_z, c, iso_val, vertices, normals, filtered_edges, weight, pattern_weight):
-        # vertices indexes
-        id0 = filtered_edges[:,0]
-        id1 = filtered_edges[:,1]
-        # vertices weight
-        w0 = weight[id0]
-        w1 = weight[id1]
-        # weight condition
-        bool_w0 = w0 < iso_val
-        bool_w1 = w1 < iso_val
-
-        # mask all edges that have one weight value below the iso value
-        mask_new_verts = np.logical_xor(bool_w0, bool_w1)
-        out_array = np.array([[0]]).astype(type='float64', order='A')
-        if not mask_new_verts.any(): return out_array, out_array, out_array
-
-        id0 = id0[mask_new_verts]
-        id1 = id1[mask_new_verts]
-        # filter arrays
-        v0 = vertices[id0]
-        v1 = vertices[id1]
-        n0 = normals[id0]
-        n1 = normals[id1]
-        w0 = w0[mask_new_verts]
-        w1 = w1[mask_new_verts]
-        pattern0 = pattern_weight[id0]
-        pattern1 = pattern_weight[id1]
-        param = (iso_val-w0)/(w1-w0)
-        # pattern displace
-        #mult = 1 if c%2 == 0 else -1
-        if c%(in_steps + out_steps) < in_steps:
-            mult = -in_displace
-        else:
-            mult = out_displace
-        pattern_value = pattern0 + (pattern1-pattern0)*param
-        disp = pattern_value * mult
-        param2 = np.expand_dims(param,axis=1)
-        #param = param2
-        disp2 = np.expand_dims(disp,axis=1)
-        #disp = disp2
-        verts = v0 + (v1-v0)*param2
-        norm = n0 + (n1-n0)*param2
-        if limit_z:
-            norm2 = np.expand_dims(norm[:,2], axis=1)
-            limit_mult = 1-np.absolute(norm2)
-            disp2 *= limit_mult
-        verts = verts + norm*disp2
-
-        # indexes of edges with new vertices
-        edges_index = filtered_edges[mask_new_verts][:,2]
-
-        # remove all edges completely below the iso value
-        mask_edges = np.logical_not(np.logical_and(bool_w0, bool_w1))
-
-        _filtered_edges = filtered_edges[mask_edges]
-        filtered_edges = _filtered_edges.astype(type='float64', order='A')
-
-        _edges_index = np.expand_dims(edges_index,axis=0)
-        _edges_index = _edges_index.astype(type='float64', order='A')
-
-        _verts = verts.astype(type='float64', order='A')
-        return _filtered_edges, _edges_index, _verts
-
-def contour_edges_pattern(operator, c, verts_count, iso_val, vertices, normals, filtered_edges, weight, pattern_weight):
-    # vertices indexes
-    id0 = filtered_edges[:,0]
-    id1 = filtered_edges[:,1]
-    # vertices weight
-    w0 = weight[id0]
-    w1 = weight[id1]
-    # weight condition
-    bool_w0 = w0 < iso_val
-    bool_w1 = w1 < iso_val
-
-    # mask all edges that have one weight value below the iso value
-    mask_new_verts = np.logical_xor(bool_w0, bool_w1)
-    if not mask_new_verts.any(): return np.array([[None]]), {}, np.array([[None]])
-
-    id0 = id0[mask_new_verts]
-    id1 = id1[mask_new_verts]
-    # filter arrays
-    v0 = vertices[id0]
-    v1 = vertices[id1]
-    n0 = normals[id0]
-    n1 = normals[id1]
-    w0 = w0[mask_new_verts]
-    w1 = w1[mask_new_verts]
-    pattern0 = pattern_weight[id0]
-    pattern1 = pattern_weight[id1]
-    param = (iso_val-w0)/(w1-w0)
-    # pattern displace
-    #mult = 1 if c%2 == 0 else -1
-    if c%(operator.in_steps + operator.out_steps) < operator.in_steps:
-        mult = -operator.in_displace
-    else:
-        mult = operator.out_displace
-    pattern_value = pattern0 + (pattern1-pattern0)*param
-    disp = pattern_value * mult
-    param = np.expand_dims(param,axis=1)
-    disp = np.expand_dims(disp,axis=1)
-    verts = v0 + (v1-v0)*param
-    norm = n0 + (n1-n0)*param
-    if operator.limit_z: disp *= 1-abs(np.expand_dims(norm[:,2], axis=1))
-    verts = verts + norm*disp
-
-    # indexes of edges with new vertices
-    edges_index = filtered_edges[mask_new_verts][:,2]
-
-    # remove all edges completely below the iso value
-    #mask_edges = np.logical_not(np.logical_and(bool_w0, bool_w1))
-    #filtered_edges = filtered_edges[mask_edges]
-    return filtered_edges, edges_index, verts
-
-def contour_edges_pattern_eval(operator, c, verts_count, iso_val, vertices, normals, filtered_edges, weight, pattern_weight):
-    # vertices indexes
-    id0 = eval('filtered_edges[:,0]')
-    id1 = eval('filtered_edges[:,1]')
-    # vertices weight
-    w0 = eval('weight[id0]')
-    w1 = eval('weight[id1]')
-    # weight condition
-    bool_w0 = ne.evaluate('w0 < iso_val')
-    bool_w1 = ne.evaluate('w1 < iso_val')
-
-    # mask all edges that have one weight value below the iso value
-    mask_new_verts = eval('np.logical_xor(bool_w0, bool_w1)')
-    if not mask_new_verts.any(): return np.array([[None]]), {}, np.array([[None]])
-
-    id0 = eval('id0[mask_new_verts]')
-    id1 = eval('id1[mask_new_verts]')
-    # filter arrays
-    v0 = eval('vertices[id0]')
-    v1 = eval('vertices[id1]')
-    n0 = eval('normals[id0]')
-    n1 = eval('normals[id1]')
-    w0 = eval('w0[mask_new_verts]')
-    w1 = eval('w1[mask_new_verts]')
-    pattern0 = eval('pattern_weight[id0]')
-    pattern1 = eval('pattern_weight[id1]')
-    param = ne.evaluate('(iso_val-w0)/(w1-w0)')
-    # pattern displace
-    #mult = 1 if c%2 == 0 else -1
-    if c%(operator.in_steps + operator.out_steps) < operator.in_steps:
-        mult = -operator.in_displace
-    else:
-        mult = operator.out_displace
-    pattern_value = eval('pattern0 + (pattern1-pattern0)*param')
-    disp = ne.evaluate('pattern_value * mult')
-    param = eval('np.expand_dims(param,axis=1)')
-    disp = eval('np.expand_dims(disp,axis=1)')
-    verts = ne.evaluate('v0 + (v1-v0)*param')
-    norm = ne.evaluate('n0 + (n1-n0)*param')
-    if operator.limit_z:
-        mult = eval('1-abs(np.expand_dims(norm[:,2], axis=1))')
-        disp = ne.evaluate('disp * mult')
-    verts = ne.evaluate('verts + norm*disp')
-
-    # indexes of edges with new vertices
-    edges_index = eval('filtered_edges[mask_new_verts][:,2]')
-
-    # remove all edges completely below the iso value
-    mask_edges = eval('np.logical_not(np.logical_and(bool_w0, bool_w1))')
-    filtered_edges = eval('filtered_edges[mask_edges]')
-    return filtered_edges, edges_index, verts
-
-
-def contour_bmesh(me, bm, weight, iso_val):
-    bm.verts.ensure_lookup_table()
-    bm.edges.ensure_lookup_table()
-    bm.faces.ensure_lookup_table()
-
-    # store weight values
-
-    vertices = get_vertices_numpy(me)
-    faces_mask = np.array(bm.faces)
-    filtered_edges = get_edges_id_numpy(me)
-    n_verts = len(bm.verts)
-
-    #############################
-
-    # vertices indexes
-    id0 = filtered_edges[:,0]
-    id1 = filtered_edges[:,1]
-    # vertices weight
-    w0 = weight[id0]
-    w1 = weight[id1]
-    # weight condition
-    bool_w0 = w0 < iso_val
-    bool_w1 = w1 < iso_val
-
-    # mask all edges that have one weight value below the iso value
-    mask_new_verts = np.logical_xor(bool_w0, bool_w1)
-    if not mask_new_verts.any(): return np.array([[None]]), {}, np.array([[None]])
-
-    id0 = id0[mask_new_verts]
-    id1 = id1[mask_new_verts]
-    # filter arrays
-    v0 = vertices[id0]
-    v1 = vertices[id1]
-    w0 = w0[mask_new_verts]
-    w1 = w1[mask_new_verts]
-    param = (iso_val-w0)/(w1-w0)
-    param = np.expand_dims(param,axis=1)
-    verts = v0 + (v1-v0)*param
-
-    # indexes of edges with new vertices
-    #edges_index = filtered_edges[mask_new_verts][:,2]
-
-    edges_id = {}
-    for i, e in enumerate(filtered_edges):
-        #edges_id[id] = i + n_verts
-        edges_id['{}_{}'.format(e[0],e[1])] = i + n_verts
-        edges_id['{}_{}'.format(e[1],e[0])] = i + n_verts
-
-
-    '''
-    for e in filtered_edges:
-        id0 = e.verts[0].index
-        id1 = e.verts[1].index
-        w0 = weight[id0]
-        w1 = weight[id1]
-
-        if w0 == w1: continue
-        elif w0 > iso_val and w1 > iso_val:
-            continue
-        elif w0 < iso_val and w1 < iso_val: continue
-        elif w0 == iso_val or w1 == iso_val: continue
-        else:
-            v0 = me0.vertices[id0].co
-            v1 = me0.vertices[id1].co
-            v = v0.lerp(v1, (iso_val-w0)/(w1-w0))
-            delete_edges.append(e)
-            verts.append(v)
-            edges_id[str(id0)+"_"+str(id1)] = count
-            edges_id[str(id1)+"_"+str(id0)] = count
-            count += 1
-    '''
-
-    splitted_faces = []
-
-    switch = False
-    # splitting faces
-    for f in faces_mask:
-        # create sub-faces slots. Once a new vertex is reached it will
-        # change slot, storing the next vertices for a new face.
-        build_faces = [[],[]]
-        #switch = False
-        verts0 = list(me.polygons[f.index].vertices)
-        verts1 = list(verts0)
-        verts1.append(verts1.pop(0)) # shift list
-        for id0, id1 in zip(verts0, verts1):
-
-            # add first vertex to active slot
-            build_faces[switch].append(id0)
-
-            # try to split edge
-            try:
-                # check if the edge must be splitted
-                new_vert = edges_id['{}_{}'.format(id0,id1)]
-                # add new vertex
-                build_faces[switch].append(new_vert)
-                # if there is an open face on the other slot
-                if len(build_faces[not switch]) > 0:
-                    # store actual face
-                    splitted_faces.append(build_faces[switch])
-                    # reset actual faces and switch
-                    build_faces[switch] = []
-                    # change face slot
-                switch = not switch
-                # continue previous face
-                build_faces[switch].append(new_vert)
-            except: pass
-        if len(build_faces[not switch]) == 2:
-            build_faces[not switch].append(id0)
-        if len(build_faces[not switch]) > 2:
-            splitted_faces.append(build_faces[not switch])
-        # add last face
-        splitted_faces.append(build_faces[switch])
-
-    # adding new vertices
-    for v in verts: bm.verts.new(v)
-    bm.verts.ensure_lookup_table()
-
-    # deleting old edges/faces
-    bm.edges.ensure_lookup_table()
-    remove_edges = [bm.edges[i] for i in filtered_edges[:,2]]
-    #for e in remove_edges: bm.edges.remove(e)
-    #for e in delete_edges: bm.edges.remove(e)
-
-    bm.verts.ensure_lookup_table()
-    # adding new faces
-    missed_faces = []
-    for f in splitted_faces:
-        try:
-            face_verts = [bm.verts[i] for i in f]
-            bm.faces.new(face_verts)
-        except:
-            missed_faces.append(f)
-
-    #me = bpy.data.meshes.new('_tissue_tmp_')
-    bm.to_mesh(me)
-    weight = np.concatenate((weight, np.ones(len(verts))*iso_val))
-
-    return me, bm, weight