Added method to remove vertices from basic meshes by remapping the v/e/l/p custom...
[blender.git] / source / blender / editors / sculpt_paint / sculpt.c
index 602fe686e7e9c1c4535ae9ce97be930413a6950c..8f2d0242ea0112bafa524d7b9ff148d892ec01dc 100644 (file)
 #include <stdlib.h>
 #include <string.h>
 
-#ifdef _OPENMP
-#include <omp.h>
+/*Maybe move silhouette GL Drawing into a seperate file*/
+#include "BIF_gl.h"
+#include "BIF_glutil.h"
+#include "ED_space_api.h"
+#include "bmesh.h"
+#include "bmesh_tools.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_editmesh.h"
+#include "BKE_editmesh_bvh.h"
+#include "ED_mesh.h"
+#include "BLI_polyfill2d.h"
+#include "BLI_polyfill2d_beautify.h"
+#include "BLI_memarena.h"
+#include "BLI_heap.h"
+#include "BLI_edgehash.h"
+#include "BLI_linklist.h"
+#include "BLI_alloca.h"
+#include "BLI_array.h"
+
+#define SIL_STROKE_STORE_CHUNK 512
+/* Store bias is used to bias a close estimate since resizing is more expensive than bigger array on first allocate*/
+#define STORE_ESTIMATE_BIAS 0.1f
+/* Fillet Blur determines the fuzziness wether a vert is intersecting or not.
+ * Important for example if two shapes with the same thickness intersect. */
+#define SIL_FILLET_BLUR_MAX 0.3f
+#define SIL_FILLET_BLUR_MIN 0.001f
+
+#define DEBUG_DRAW
+#ifdef DEBUG_DRAW
+/* static void bl_debug_draw(void);*/
+/* add these locally when using these functions for testing */
+extern void bl_debug_draw_quad_clear(void);
+extern void bl_debug_draw_quad_add(const float v0[3], const float v1[3], const float v2[3], const float v3[3]);
+extern void bl_debug_draw_edge_add(const float v0[3], const float v1[3]);
+extern void bl_debug_color_set(const unsigned int col);
+
+static void bl_debug_draw_medge_add(Mesh *me, int e){
+       bl_debug_draw_edge_add(me->mvert[me->medge[e].v1].co, me->mvert[me->medge[e].v2].co);
+}
+
+static void bl_debug_draw_point(const float pos[3],const float thickness)
+{
+       float h = thickness*0.5;
+       float v1[] = {pos[0]-h,pos[1]-h,pos[2]};
+       float v2[] = {pos[0]-h,pos[1]+h,pos[2]};
+       float v3[] = {pos[0]+h,pos[1]-h,pos[2]};
+       float v4[] = {pos[0]+h,pos[1]+h,pos[2]};
+       bl_debug_draw_quad_add(v1,v2,v3,v4);
+}
 #endif
 
-#if defined(__APPLE__) && defined _OPENMP
-#include <sys/sysctl.h>
-
-/* Query how many cores not counting HT aka physical cores we've got. */
-static int system_physical_thread_count(void)
-{
-       int pcount;
-       size_t pcount_len = sizeof(pcount);
-       sysctlbyname("hw.physicalcpu", &pcount, &pcount_len, NULL, 0);
-       return pcount;
+typedef struct {
+       float bmin[3], bmax[3];
+} BB;
+
+/*     init data:
+ *     Silhouette Data */
+typedef struct SilhouetteStroke {
+       float *points;
+       float *points_v2;
+       int totvert;
+       int max_verts;
+       BB bb;
+} SilhouetteStroke;
+
+#ifdef DEBUG_DRAW
+void bl_debug_draw_BB_add(BB *bb,const unsigned int col){
+       float v1[3],v2[3],v3[3],v4[3];
+       float xd[3], yd[3], zd[3];
+
+       bl_debug_color_set(col);
+
+       xd[0] = bb->bmax[0]-bb->bmin[0];
+       xd[1] = 0.0f;
+       xd[2] = 0.0f;
+
+       yd[0] = 0.0f;
+       yd[1] = bb->bmax[1]-bb->bmin[1];
+       yd[2] = 0.0f;
+
+       zd[0] = 0.0f;
+       zd[1] = 0.0f;
+       zd[2] = bb->bmax[2]-bb->bmin[2];
+
+       copy_v3_v3(v1,bb->bmin);
+       copy_v3_v3(v2,bb->bmin);
+       add_v3_v3(v2,xd);
+       add_v3_v3v3(v3,v1,yd);
+       add_v3_v3v3(v4,v2,yd);
+
+       bl_debug_draw_edge_add(v1,v2);
+       bl_debug_draw_edge_add(v1,v3);
+       bl_debug_draw_edge_add(v2,v4);
+
+       copy_v3_v3(v1,v3);
+       copy_v3_v3(v2,v4);
+       add_v3_v3v3(v3,v1,zd);
+       add_v3_v3v3(v4,v2,zd);
+
+       bl_debug_draw_edge_add(v1,v2);
+       bl_debug_draw_edge_add(v1,v3);
+       bl_debug_draw_edge_add(v2,v4);
+
+       copy_v3_v3(v1,v3);
+       copy_v3_v3(v2,v4);
+       sub_v3_v3v3(v3,v1,yd);
+       sub_v3_v3v3(v4,v2,yd);
+
+       bl_debug_draw_edge_add(v1,v2);
+       bl_debug_draw_edge_add(v1,v3);
+       bl_debug_draw_edge_add(v2,v4);
+
+       copy_v3_v3(v1,v3);
+       copy_v3_v3(v2,v4);
+       sub_v3_v3v3(v3,v1,zd);
+       sub_v3_v3v3(v4,v2,zd);
+
+       bl_debug_draw_edge_add(v1,v2);
+       bl_debug_draw_edge_add(v1,v3);
+       bl_debug_draw_edge_add(v2,v4);
+       
 }
-#endif  /* __APPLE__ */
+#endif
+
+
+typedef enum {
+       SIL_INIT = 0,
+       SIL_DRAWING = 1,
+       SIL_OP = 2
+} SilhouetteState;
+
+typedef struct SilhouetteData {
+       ARegion *ar;            /* region that Silhouette started drawn in */
+       void *draw_handle;      /* for drawing preview loop */
+       ViewContext vc;
+       SilhouetteStroke *current_stroke;
+       Object *ob;
+       BMEditMesh *em;                 /*Triangulated stroke for spine generation*/
+       Scene *scene;
+
+       float add_col[3];               /* preview color */
+       float last_mouse_pos[2];
+
+       SilhouetteState state;  /* Operator state */
+
+       float depth;                                    /* Depth or thickness of the generated shape */
+       float smoothness;                               /* Smoothness of the generated shape */
+       int resolution;                                 /* Subdivision of the shape*/
+       float anchor[3];                                /* Origin point of the reference plane */
+       float z_vec[3];                                 /* Orientation of the reference plane */
+       MeshElemMap *emap;                              /* Original Mesh vert -> edges map */
+       GHash *i_edges;                                 /* Edges crossing the both shapes. (only orig mesh)*/
+       int *fillet_ring_orig;                  /* ring_edges to connect to in the orig mesh */
+       int *fillet_ring_orig_start;    /* start positions to each individual ring */
+       int *fillet_ring_new;                   /* ring_edges to connect to in the new mesh */
+       int *fillet_ring_new_start;             /* start positions to each individual ring */
+       int num_rings, fillet_ring_tot;
+       int *inter_edges;                               /* edges crossing the two shapes */
+       int num_inter_edges;                    /* number of edges crossing */
+       BB *fillet_ring_bbs;                            /* every ring gets a Bounding box to check intersection with branches */
+} SilhouetteData;
 
 /** \name Tool Capabilities
  *
@@ -159,10 +304,10 @@ static bool sculpt_tool_is_proxy_used(const char sculpt_tool)
 /**
  * Test whether the #StrokeCache.sculpt_normal needs update in #do_brush_action
  */
-static int sculpt_brush_needs_normal(const Brush *brush)
+static int sculpt_brush_needs_normal(const Brush *brush, float normal_weight)
 {
        return ((SCULPT_TOOL_HAS_NORMAL_WEIGHT(brush->sculpt_tool) &&
-                (brush->normal_weight > 0)) ||
+                (normal_weight > 0.0f)) ||
 
                ELEM(brush->sculpt_tool,
                     SCULPT_TOOL_BLOB,
@@ -175,9 +320,21 @@ static int sculpt_brush_needs_normal(const Brush *brush)
 
                (brush->mtex.brush_map_mode == MTEX_MAP_MODE_AREA));
 }
-
 /** \} */
 
+static bool sculpt_brush_needs_rake_rotation(const Brush *brush)
+{
+       return SCULPT_TOOL_HAS_RAKE(brush->sculpt_tool) && (brush->rake_factor != 0.0f);
+}
+
+/* Factor of brush to have rake point following behind
+ * (could be configurable but this is reasonable default). */
+#define SCULPT_RAKE_BRUSH_FACTOR 0.25f
+
+struct SculptRakeData {
+       float follow_dist;
+       float follow_co[3];
+} SculptRakeData;
 
 typedef enum StrokeFlags {
        CLIP_X = 1,
@@ -209,6 +366,7 @@ typedef struct StrokeCache {
        float pressure;
        float mouse[2];
        float bstrength;
+       float normal_weight;  /* from brush (with optional override) */
 
        /* The rest is temporary storage that isn't saved as a property */
 
@@ -218,13 +376,18 @@ typedef struct StrokeCache {
        float projection_mat[4][4];
 
        /* Clean this up! */
-       ViewContext *vc;
-       Brush *brush;
+       struct ViewContext *vc;
+       struct Brush *brush;
 
        float special_rotation;
        float grab_delta[3], grab_delta_symmetry[3];
        float old_grab_location[3], orig_grab_location[3];
 
+       /* screen-space rotation defined by mouse motion */
+       float   rake_rotation[4], rake_rotation_symmetry[4];
+       bool is_rake_rotation_valid;
+       struct SculptRakeData rake_data;
+
        int symmetry; /* Symmetry index between 0 and 7 bit combo 0 is Brush only;
                       * 1 is X mirror; 2 is Y mirror; 3 is XY; 4 is Z; 5 is XZ; 6 is YZ; 7 is XYZ */
        int mirror_symmetry_pass; /* the symmetry pass we are currently on between 0 and 7*/
@@ -252,7 +415,7 @@ typedef struct StrokeCache {
        float anchored_location[3];
 
        float vertex_rotation; /* amount to rotate the vertices when using rotate brush */
-       Dial *dial;
+       struct Dial *dial;
        
        char saved_active_brush_name[MAX_ID_NAME];
        char saved_mask_brush_tool;
@@ -269,6 +432,19 @@ typedef struct StrokeCache {
        rcti current_r; /* current redraw rectangle */
 } StrokeCache;
 
+/* reduce brush spacing step size when the geometry curves away from the view */
+static void set_adaptive_space_factor(Sculpt *sd)
+{
+       Brush *brush = BKE_paint_brush(&(sd->paint));
+
+       /*TODO: Reasonable 2D View 3D conversion
+        * Currently somewhere about 1bu / 200px
+        */
+
+       brush->adaptive_space_factor = 1.0f/200.0f;
+}
+
+
 /************** Access to original unmodified vertex data *************/
 
 typedef struct {
@@ -345,6 +521,67 @@ static void sculpt_orig_vert_data_update(SculptOrigVertData *orig_data,
        }
 }
 
+static void sculpt_rake_data_update(struct SculptRakeData *srd, const float co[3])
+{
+       float rake_dist = len_v3v3(srd->follow_co, co);
+       if (rake_dist > srd->follow_dist) {
+               interp_v3_v3v3(srd->follow_co, srd->follow_co, co, rake_dist - srd->follow_dist);
+       }
+}
+
+
+static void sculpt_rake_rotate(
+        const SculptSession *ss, const float sculpt_co[3], const float v_co[3], float factor, float r_delta[3])
+{
+       float vec_rot[3];
+
+#if 0
+       /* lerp */
+       sub_v3_v3v3(vec_rot, v_co, sculpt_co);
+       mul_qt_v3(ss->cache->rake_rotation_symmetry, vec_rot);
+       add_v3_v3(vec_rot, sculpt_co);
+       sub_v3_v3v3(r_delta, vec_rot, v_co);
+       mul_v3_fl(r_delta, factor);
+#else
+       /* slerp */
+       float q_interp[4];
+       sub_v3_v3v3(vec_rot, v_co, sculpt_co);
+
+       copy_qt_qt(q_interp, ss->cache->rake_rotation_symmetry);
+       mul_fac_qt_fl(q_interp, factor);
+       mul_qt_v3(q_interp, vec_rot);
+
+       add_v3_v3(vec_rot, sculpt_co);
+       sub_v3_v3v3(r_delta, vec_rot, v_co);
+#endif
+
+}
+
+/**
+ * Align the grab delta to the brush normal.
+ *
+ * \param grab_delta: Typically from `ss->cache->grab_delta_symmetry`.
+ */
+static void sculpt_project_v3_normal_align(SculptSession *ss, const float normal_weight, float grab_delta[3])
+{
+       /* signed to support grabbing in (to make a hole) as well as out. */
+       const float len_signed = dot_v3v3(ss->cache->sculpt_normal_symm, grab_delta);
+
+       /* this scale effectively projects the offset so dragging follows the cursor,
+        * as the normal points towards the view, the scale increases. */
+       float len_view_scale;
+       {
+               float view_aligned_normal[3];
+               project_plane_v3_v3v3(view_aligned_normal, ss->cache->sculpt_normal_symm, ss->cache->view_normal);
+               len_view_scale = fabsf(dot_v3v3(view_aligned_normal, ss->cache->sculpt_normal_symm));
+               len_view_scale = (len_view_scale > FLT_EPSILON) ? 1.0f / len_view_scale : 1.0f;
+       }
+
+       mul_v3_fl(grab_delta, 1.0f - normal_weight);
+       madd_v3_v3fl(grab_delta, ss->cache->sculpt_normal_symm, (len_signed * normal_weight) * len_view_scale);
+}
+
+
 /** \name SculptProjectVector
  *
  * Fast-path for #project_plane_v3_v3v3
@@ -419,7 +656,7 @@ typedef struct SculptThreadedTaskData {
        Sculpt *sd;
        Object *ob;
        Brush *brush;
-    PBVHNode **nodes;
+       PBVHNode **nodes;
        int totnode;
 
        /* Data specific to some callbacks. */
@@ -430,6 +667,9 @@ typedef struct SculptThreadedTaskData {
        float strength;
        bool smooth_mask;
        bool has_bm_orco;
+       SilhouetteData *sil;
+       int *v_to_rm; /* Shared array handle access with mutex! */
+       int num_v_to_rm;
 
        SculptProjectVector *spvc;
        float *offset;
@@ -513,7 +753,7 @@ static void paint_mesh_restore_co(Sculpt *sd, Object *ob)
 
        BLI_task_parallel_range(
                    0, totnode, &data, paint_mesh_restore_co_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && !ss->bm && totnode > SCULPT_OMP_LIMIT));
+                   ((sd->flags & SCULPT_USE_OPENMP) && !ss->bm && totnode > SCULPT_THREADED_LIMIT));
 
        if (nodes)
                MEM_freeN(nodes);
@@ -592,6 +832,7 @@ typedef struct SculptBrushTest {
        float radius_squared;
        float location[3];
        float dist;
+       int mirror_symmetry_pass;
 
        /* View3d clipping - only set rv3d for clipping */
        RegionView3D *clip_rv3d;
@@ -605,6 +846,7 @@ static void sculpt_brush_test_init(SculptSession *ss, SculptBrushTest *test)
        copy_v3_v3(test->location, ss->cache->location);
        test->dist = 0.0f;   /* just for initialize */
 
+       test->mirror_symmetry_pass = ss->cache->mirror_symmetry_pass;
 
        if (rv3d->rflag & RV3D_CLIPPING) {
                test->clip_rv3d = rv3d;
@@ -617,7 +859,12 @@ static void sculpt_brush_test_init(SculptSession *ss, SculptBrushTest *test)
 BLI_INLINE bool sculpt_brush_test_clipping(const SculptBrushTest *test, const float co[3])
 {
        RegionView3D *rv3d = test->clip_rv3d;
-       return (rv3d && (ED_view3d_clipping_test(rv3d, co, true)));
+       if (!rv3d) {
+               return false;
+       }
+       float symm_co[3];
+       flip_v3_v3(symm_co, co, test->mirror_symmetry_pass);
+       return ED_view3d_clipping_test(rv3d, symm_co, true);
 }
 
 static bool sculpt_brush_test(SculptBrushTest *test, const float co[3])
@@ -752,10 +999,9 @@ static float calc_overlap(StrokeCache *cache, const char symm, const char axis,
        flip_v3_v3(mirror, cache->true_location, symm);
 
        if (axis != 0) {
-               float mat[4][4];
-               unit_m4(mat);
-               rotate_m4(mat, axis, angle);
-               mul_m4_v3(mat, mirror);
+               float mat[3][3];
+               axis_angle_to_mat3_single(mat, axis, angle);
+               mul_m3_v3(mat, mirror);
        }
 
        /* distsq = len_squared_v3v3(mirror, cache->traced_location); */
@@ -963,13 +1209,13 @@ static void calc_area_center(
 
        SculptThreadedTaskData data = {
                .sd = sd, .ob = ob, .nodes = nodes, .totnode = totnode,
-           .has_bm_orco = has_bm_orco, .area_cos = area_cos, .area_nos = NULL, .count = count,
+               .has_bm_orco = has_bm_orco, .area_cos = area_cos, .area_nos = NULL, .count = count,
        };
        BLI_mutex_init(&data.mutex);
 
        BLI_task_parallel_range(
                    0, totnode, &data, calc_area_normal_and_center_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
 
        BLI_mutex_end(&data.mutex);
 
@@ -1003,13 +1249,13 @@ static void calc_area_normal(
 
        SculptThreadedTaskData data = {
                .sd = sd, .ob = ob, .nodes = nodes, .totnode = totnode,
-           .has_bm_orco = has_bm_orco, .area_cos = NULL, .area_nos = area_nos, .count = count,
+               .has_bm_orco = has_bm_orco, .area_cos = NULL, .area_nos = area_nos, .count = count,
        };
        BLI_mutex_init(&data.mutex);
 
        BLI_task_parallel_range(
                    0, totnode, &data, calc_area_normal_and_center_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
 
        BLI_mutex_end(&data.mutex);
 
@@ -1041,13 +1287,13 @@ static void calc_area_normal_and_center(
 
        SculptThreadedTaskData data = {
                .sd = sd, .ob = ob, .nodes = nodes, .totnode = totnode,
-           .has_bm_orco = has_bm_orco, .area_cos = area_cos, .area_nos = area_nos, .count = count,
+               .has_bm_orco = has_bm_orco, .area_cos = area_cos, .area_nos = area_nos, .count = count,
        };
        BLI_mutex_init(&data.mutex);
 
        BLI_task_parallel_range(
                    0, totnode, &data, calc_area_normal_and_center_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
 
        BLI_mutex_end(&data.mutex);
 
@@ -1076,7 +1322,9 @@ static void calc_area_normal_and_center(
 /* Return modified brush strength. Includes the direction of the brush, positive
  * values pull vertices, negative values push. Uses tablet pressure and a
  * special multiplier found experimentally to scale the strength factor. */
-static float brush_strength(const Sculpt *sd, const StrokeCache *cache, const float feather, const UnifiedPaintSettings *ups)
+static float brush_strength(
+        const Sculpt *sd, const StrokeCache *cache,
+        const float feather, const UnifiedPaintSettings *ups)
 {
        const Scene *scene = cache->vc->scene;
        const Brush *brush = BKE_paint_brush((Paint *)&sd->paint);
@@ -1168,18 +1416,21 @@ static float brush_strength(const Sculpt *sd, const StrokeCache *cache, const fl
 
 /* Return a multiplier for brush strength on a particular vertex. */
 static float tex_strength(SculptSession *ss, Brush *br,
-                          const float point[3],
+                          const float brush_point[3],
                           const float len,
                           const short vno[3],
                           const float fno[3],
-                          const float mask)
+                          const float mask,
+                          const int thread_id)
 {
        StrokeCache *cache = ss->cache;
        const Scene *scene = cache->vc->scene;
        MTex *mtex = &br->mtex;
        float avg = 1;
        float rgba[4];
-       int thread_num;
+       float point[3];
+
+       sub_v3_v3v3(point, brush_point, cache->plane_offset);
 
        if (!mtex->tex) {
                avg = 1;
@@ -1222,12 +1473,7 @@ static float tex_strength(SculptSession *ss, Brush *br,
                        x += br->mtex.ofs[0];
                        y += br->mtex.ofs[1];
 
-#ifdef _OPENMP
-                       thread_num = omp_get_thread_num();
-#else
-                       thread_num = 0;
-#endif
-                       avg = paint_get_tex_pixel(&br->mtex, x, y, ss->tex_pool, thread_num);
+                       avg = paint_get_tex_pixel(&br->mtex, x, y, ss->tex_pool, thread_id);
 
                        avg += br->texture_sample_bias;
                }
@@ -1564,7 +1810,8 @@ typedef struct SculptDoBrushSmoothGridDataChunk {
        size_t tmpgrid_size;
 } SculptDoBrushSmoothGridDataChunk;
 
-static void do_smooth_brush_mesh_task_cb(void *userdata, int n)
+static void do_smooth_brush_mesh_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -1585,7 +1832,8 @@ static void do_smooth_brush_mesh_task_cb(void *userdata, int n)
                if (sculpt_brush_test(&test, vd.co)) {
                        const float fade = bstrength * tex_strength(
                                               ss, brush, vd.co, test.dist, vd.no, vd.fno,
-                                              smooth_mask ? 0.0f : (vd.mask ? *vd.mask : 0.0f));
+                                              smooth_mask ? 0.0f : (vd.mask ? *vd.mask : 0.0f),
+                                              thread_id);
                        if (smooth_mask) {
                                float val = neighbor_average_mask(ss, vd.vert_indices[vd.i]) - *vd.mask;
                                val *= fade * bstrength;
@@ -1610,7 +1858,8 @@ static void do_smooth_brush_mesh_task_cb(void *userdata, int n)
        BKE_pbvh_vertex_iter_end;
 }
 
-static void do_smooth_brush_bmesh_task_cb(void *userdata, int n)
+static void do_smooth_brush_bmesh_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -1630,7 +1879,8 @@ static void do_smooth_brush_bmesh_task_cb(void *userdata, int n)
        {
                if (sculpt_brush_test(&test, vd.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, smooth_mask ? 0.0f : *vd.mask);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, smooth_mask ? 0.0f : *vd.mask,
+                                              thread_id);
                        if (smooth_mask) {
                                float val = bmesh_neighbor_average_mask(vd.bm_vert, vd.cd_vert_mask_offset) - *vd.mask;
                                val *= fade * bstrength;
@@ -1656,7 +1906,7 @@ static void do_smooth_brush_bmesh_task_cb(void *userdata, int n)
 }
 
 static void do_smooth_brush_multires_task_cb_ex(
-        void *userdata, void *userdata_chunk, const int n, const int UNUSED(thread_id))
+        void *userdata, void *userdata_chunk, const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptDoBrushSmoothGridDataChunk *data_chunk = userdata_chunk;
@@ -1768,7 +2018,8 @@ static void do_smooth_brush_multires_task_cb_ex(
 
                                if (sculpt_brush_test(&test, co)) {
                                        const float strength_mask = (smooth_mask ? 0.0f : *mask);
-                                       const float fade = bstrength * tex_strength(ss, brush, co, test.dist, NULL, fno, strength_mask);
+                                       const float fade = bstrength * tex_strength(
+                                                              ss, brush, co, test.dist, NULL, fno, strength_mask, thread_id);
                                        float f = 1.0f / 16.0f;
 
                                        if (x == 0 || x == gridsize - 1)
@@ -1818,6 +2069,8 @@ static void smooth(
                return;
        }
 
+       set_adaptive_space_factor(sd);
+
        for (iteration = 0; iteration <= count; ++iteration) {
                const float strength = (iteration != count) ? 1.0f : last;
 
@@ -1842,20 +2095,20 @@ static void smooth(
 
                                BLI_task_parallel_range_ex(
                                            0, totnode, &data, data_chunk, size, do_smooth_brush_multires_task_cb_ex,
-                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT), false);
+                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 
                                MEM_freeN(data_chunk);
                                break;
                        }
                        case PBVH_FACES:
-                               BLI_task_parallel_range(
-                                           0, totnode, &data, do_smooth_brush_mesh_task_cb,
-                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                               BLI_task_parallel_range_ex(
+                                           0, totnode, &data, NULL, 0, do_smooth_brush_mesh_task_cb_ex,
+                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
                                break;
                        case PBVH_BMESH:
-                               BLI_task_parallel_range(
-                                           0, totnode, &data, do_smooth_brush_bmesh_task_cb,
-                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                               BLI_task_parallel_range_ex(
+                                           0, totnode, &data, NULL, 0, do_smooth_brush_bmesh_task_cb_ex,
+                                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
                                break;
                }
 
@@ -1870,7 +2123,8 @@ static void do_smooth_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnod
        smooth(sd, ob, nodes, totnode, ss->cache->bstrength, false);
 }
 
-static void do_mask_brush_draw_task_cb(void *userdata, const int n)
+static void do_mask_brush_draw_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -1885,7 +2139,7 @@ static void do_mask_brush_draw_task_cb(void *userdata, const int n)
        BKE_pbvh_vertex_iter_begin(ss->pbvh, data->nodes[n], vd, PBVH_ITER_UNIQUE)
        {
                if (sculpt_brush_test(&test, vd.co)) {
-                       const float fade = tex_strength(ss, brush, vd.co, test.dist, vd.no, vd.fno, 0.0f);
+                       const float fade = tex_strength(ss, brush, vd.co, test.dist, vd.no, vd.fno, 0.0f, thread_id);
 
                        (*vd.mask) += fade * bstrength;
                        CLAMP(*vd.mask, 0, 1);
@@ -1901,14 +2155,16 @@ static void do_mask_brush_draw(Sculpt *sd, Object *ob, PBVHNode **nodes, int tot
 {
        Brush *brush = BKE_paint_brush(&sd->paint);
 
+       set_adaptive_space_factor(sd);
+
        /* threaded loop over nodes */
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_mask_brush_draw_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_mask_brush_draw_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
 static void do_mask_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
@@ -1926,7 +2182,8 @@ static void do_mask_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
        }
 }
 
-static void do_draw_brush_task_cb(void *userdata, const int n)
+static void do_draw_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -1945,7 +2202,8 @@ static void do_draw_brush_task_cb(void *userdata, const int n)
        {
                if (sculpt_brush_test(&test, vd.co)) {
                        /* offset vertex */
-                       const float fade = tex_strength(ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                       const float fade = tex_strength(
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], offset, fade);
 
@@ -1968,18 +2226,21 @@ static void do_draw_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
        mul_v3_v3(offset, ss->cache->scale);
        mul_v3_fl(offset, bstrength);
 
+       set_adaptive_space_factor(sd);
+
        /* threaded loop over nodes */
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .offset = offset,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_draw_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_draw_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_crease_brush_task_cb(void *userdata, const int n)
+static void do_crease_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2001,7 +2262,7 @@ static void do_crease_brush_task_cb(void *userdata, const int n)
                if (sculpt_brush_test(&test, vd.co)) {
                        /* offset vertex */
                        const float fade = tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
                        float val1[3];
                        float val2[3];
 
@@ -2055,18 +2316,21 @@ static void do_crease_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnod
         * Without this we get a 'flat' surface surrounding the pinch */
        sculpt_project_v3_cache_init(&spvc, ss->cache->sculpt_normal_symm);
 
+       set_adaptive_space_factor(sd);
+
        /* threaded loop over nodes */
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .spvc = &spvc, .offset = offset, .flippedbstrength = flippedbstrength,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_crease_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_crease_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_pinch_brush_task_cb(void *userdata, const int n)
+static void do_pinch_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2085,7 +2349,7 @@ static void do_pinch_brush_task_cb(void *userdata, const int n)
        {
                if (sculpt_brush_test(&test, vd.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
                        float val[3];
 
                        sub_v3_v3v3(val, test.location, vd.co);
@@ -2102,16 +2366,19 @@ static void do_pinch_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode
 {
        Brush *brush = BKE_paint_brush(&sd->paint);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_pinch_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_pinch_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_grab_brush_task_cb(void *userdata, const int n)
+static void do_grab_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2136,7 +2403,8 @@ static void do_grab_brush_task_cb(void *userdata, const int n)
 
                if (sculpt_brush_test(&test, orig_data.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f,
+                                              thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], grab_delta, fade);
 
@@ -2152,29 +2420,27 @@ static void do_grab_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
        SculptSession *ss = ob->sculpt;
        Brush *brush = BKE_paint_brush(&sd->paint);
        float grab_delta[3];
-       float len;
 
        copy_v3_v3(grab_delta, ss->cache->grab_delta_symmetry);
 
-       len = len_v3(grab_delta);
-
-       if (brush->normal_weight > 0) {
-               mul_v3_fl(ss->cache->sculpt_normal_symm, len * brush->normal_weight);
-               mul_v3_fl(grab_delta, 1.0f - brush->normal_weight);
-               add_v3_v3(grab_delta, ss->cache->sculpt_normal_symm);
+       if (ss->cache->normal_weight > 0.0f) {
+               sculpt_project_v3_normal_align(ss, ss->cache->normal_weight, grab_delta);
        }
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .grab_delta = grab_delta,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_grab_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_grab_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_nudge_brush_task_cb(void *userdata, const int n)
+static void do_nudge_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2193,8 +2459,8 @@ static void do_nudge_brush_task_cb(void *userdata, const int n)
        BKE_pbvh_vertex_iter_begin(ss->pbvh, data->nodes[n], vd, PBVH_ITER_UNIQUE)
        {
                if (sculpt_brush_test(&test, vd.co)) {
-                       const float fade = bstrength * tex_strength(ss, brush, vd.co, test.dist,
-                                                                   vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                       const float fade = bstrength * tex_strength(
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], cono, fade);
 
@@ -2217,27 +2483,35 @@ static void do_nudge_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode
        cross_v3_v3v3(tmp, ss->cache->sculpt_normal_symm, grab_delta);
        cross_v3_v3v3(cono, tmp, ss->cache->sculpt_normal_symm);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .cono = cono,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_nudge_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_nudge_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_snake_hook_brush_task_cb(void *userdata, const int n)
+static void do_snake_hook_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
        Brush *brush = data->brush;
+       SculptProjectVector *spvc = data->spvc;
        const float *grab_delta = data->grab_delta;
 
        PBVHVertexIter vd;
        SculptBrushTest test;
        float (*proxy)[3];
        const float bstrength = ss->cache->bstrength;
+       const bool do_rake_rotation = ss->cache->is_rake_rotation_valid;
+       const bool do_pinch = (brush->crease_pinch_factor != 0.5f);
+       const float pinch = do_pinch ?
+               (2.0f * (0.5f - brush->crease_pinch_factor) * (len_v3(grab_delta) / ss->cache->radius)) : 0.0f;
 
        proxy = BKE_pbvh_node_add_proxy(ss->pbvh, data->nodes[n])->co;
 
@@ -2247,10 +2521,41 @@ static void do_snake_hook_brush_task_cb(void *userdata, const int n)
        {
                if (sculpt_brush_test(&test, vd.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], grab_delta, fade);
 
+                       /* negative pinch will inflate, helps maintain volume */
+                       if (do_pinch) {
+                               float delta_pinch_init[3], delta_pinch[3];
+
+                               sub_v3_v3v3(delta_pinch, vd.co, test.location);
+
+                               /* important to calculate based on the grabbed location (intentionally ignore fade here). */
+                               add_v3_v3(delta_pinch, grab_delta);
+
+                               sculpt_project_v3(spvc, delta_pinch, delta_pinch);
+
+                               copy_v3_v3(delta_pinch_init, delta_pinch);
+
+                               float pinch_fade = pinch * fade;
+                               /* when reducing, scale reduction back by how close to the center we are,
+                                * so we don't pinch into nothingness */
+                               if (pinch > 0.0f) {
+                                       /* square to have even less impact for close vertices */
+                                       pinch_fade *= pow2f(min_ff(1.0f, len_v3(delta_pinch) / ss->cache->radius));
+                               }
+                               mul_v3_fl(delta_pinch, 1.0f + pinch_fade);
+                               sub_v3_v3v3(delta_pinch, delta_pinch_init, delta_pinch);
+                               add_v3_v3(proxy[vd.i], delta_pinch);
+                       }
+
+                       if (do_rake_rotation) {
+                               float delta_rotate[3];
+                               sculpt_rake_rotate(ss, test.location, vd.co, fade, delta_rotate);
+                               add_v3_v3(proxy[vd.i], delta_rotate);
+                       }
+
                        if (vd.mvert)
                                vd.mvert->flag |= ME_VERT_PBVH_UPDATE;
                }
@@ -2264,32 +2569,38 @@ static void do_snake_hook_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int to
        Brush *brush = BKE_paint_brush(&sd->paint);
        const float bstrength = ss->cache->bstrength;
        float grab_delta[3];
-       float len;
 
-       copy_v3_v3(grab_delta, ss->cache->grab_delta_symmetry);
+       SculptProjectVector spvc;
 
-       len = len_v3(grab_delta);
+       copy_v3_v3(grab_delta, ss->cache->grab_delta_symmetry);
 
        if (bstrength < 0)
                negate_v3(grab_delta);
 
-       if (brush->normal_weight > 0) {
-               mul_v3_fl(ss->cache->sculpt_normal_symm, len * brush->normal_weight);
-               mul_v3_fl(grab_delta, 1.0f - brush->normal_weight);
-               add_v3_v3(grab_delta, ss->cache->sculpt_normal_symm);
+       if (ss->cache->normal_weight > 0.0f) {
+               sculpt_project_v3_normal_align(ss, ss->cache->normal_weight, grab_delta);
+       }
+
+       /* optionally pinch while painting */
+       if (brush->crease_pinch_factor != 0.5f) {
+               sculpt_project_v3_cache_init(&spvc, grab_delta);
        }
 
+
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
-           .grab_delta = grab_delta,
+           .spvc = &spvc, .grab_delta = grab_delta,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_snake_hook_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_snake_hook_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_thumb_brush_task_cb(void *userdata, const int n)
+static void do_thumb_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2314,7 +2625,8 @@ static void do_thumb_brush_task_cb(void *userdata, const int n)
 
                if (sculpt_brush_test(&test, orig_data.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f,
+                                              thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], cono, fade);
 
@@ -2337,17 +2649,20 @@ static void do_thumb_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode
        cross_v3_v3v3(tmp, ss->cache->sculpt_normal_symm, grab_delta);
        cross_v3_v3v3(cono, tmp, ss->cache->sculpt_normal_symm);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .cono = cono,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_thumb_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_thumb_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_rotate_brush_task_cb(void *userdata, const int n)
+static void do_rotate_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2373,7 +2688,8 @@ static void do_rotate_brush_task_cb(void *userdata, const int n)
                if (sculpt_brush_test(&test, orig_data.co)) {
                        float vec[3], rot[3][3];
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, orig_data.co, test.dist, orig_data.no, NULL, vd.mask ? *vd.mask : 0.0f,
+                                              thread_id);
 
                        sub_v3_v3v3(vec, orig_data.co, ss->cache->location);
                        axis_angle_normalized_to_mat3(rot, ss->cache->sculpt_normal_symm, angle * fade);
@@ -2396,17 +2712,20 @@ static void do_rotate_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnod
        static const int flip[8] = { 1, -1, -1, 1, -1, 1, 1, -1 };
        const float angle = ss->cache->vertex_rotation * flip[ss->cache->mirror_symmetry_pass];
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .angle = angle,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_rotate_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_rotate_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_layer_brush_task_cb(void *userdata, const int n)
+static void do_layer_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2438,7 +2757,7 @@ static void do_layer_brush_task_cb(void *userdata, const int n)
 
                if (sculpt_brush_test(&test, orig_data.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
                        float *disp = &layer_disp[vd.i];
                        float val[3];
 
@@ -2477,20 +2796,23 @@ static void do_layer_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode
 
        mul_v3_v3v3(offset, ss->cache->scale, ss->cache->sculpt_normal_symm);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .offset = offset,
        };
        BLI_mutex_init(&data.mutex);
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_layer_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_layer_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 
        BLI_mutex_end(&data.mutex);
 }
 
-static void do_inflate_brush_task_cb(void *userdata, const int n)
+static void do_inflate_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2509,7 +2831,7 @@ static void do_inflate_brush_task_cb(void *userdata, const int n)
        {
                if (sculpt_brush_test(&test, vd.co)) {
                        const float fade = bstrength * tex_strength(
-                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, test.dist, vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
                        float val[3];
 
                        if (vd.fno)
@@ -2531,13 +2853,15 @@ static void do_inflate_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totno
 {
        Brush *brush = BKE_paint_brush(&sd->paint);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_inflate_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_inflate_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
 static void calc_sculpt_plane(
@@ -2665,7 +2989,8 @@ static float get_offset(Sculpt *sd, SculptSession *ss)
        return rv;
 }
 
-static void do_flatten_brush_task_cb(void *userdata, const int n)
+static void do_flatten_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2694,7 +3019,8 @@ static void do_flatten_brush_task_cb(void *userdata, const int n)
 
                        if (plane_trim(ss->cache, brush, val)) {
                                const float fade = bstrength * tex_strength(
-                                                      ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                                      ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f,
+                                                      thread_id);
 
                                mul_v3_v3fl(proxy[vd.i], val, fade);
 
@@ -2728,17 +3054,20 @@ static void do_flatten_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totno
        mul_v3_fl(temp, displace);
        add_v3_v3(area_co, temp);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .area_no = area_no, .area_co = area_co,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_flatten_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_flatten_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_clay_brush_task_cb(void *userdata, const int n)
+static void do_clay_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2771,7 +3100,8 @@ static void do_clay_brush_task_cb(void *userdata, const int n)
                                        /* note, the normal from the vertices is ignored,
                                         * causes glitch with planes, see: T44390 */
                                        const float fade = bstrength * tex_strength(
-                                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f,
+                                                              thread_id);
 
                                        mul_v3_v3fl(proxy[vd.i], val, fade);
 
@@ -2809,17 +3139,20 @@ static void do_clay_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
 
        /* add_v3_v3v3(p, ss->cache->location, area_no); */
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .area_no = area_no, .area_co = area_co,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_clay_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_clay_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_clay_strips_brush_task_cb(void *userdata, const int n)
+static void do_clay_strips_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2854,7 +3187,7 @@ static void do_clay_strips_brush_task_cb(void *userdata, const int n)
                                         * causes glitch with planes, see: T44390 */
                                        const float fade = bstrength * tex_strength(
                                                               ss, brush, vd.co, ss->cache->radius * test.dist,
-                                                              vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                                              vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
 
                                        mul_v3_v3fl(proxy[vd.i], val, fade);
 
@@ -2875,7 +3208,7 @@ static void do_clay_strips_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int t
        const bool flip = (ss->cache->bstrength < 0);
        const float radius    = flip ? -ss->cache->radius : ss->cache->radius;
        const float offset    = get_offset(sd, ss);
-       const float displace  = radius * (0.25f + offset);;
+       const float displace  = radius * (0.25f + offset);
 
        float area_no_sp[3];  /* the sculpt-plane normal (whatever its set to) */
        float area_no[3];     /* geometry normal */
@@ -2917,17 +3250,20 @@ static void do_clay_strips_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int t
        mul_m4_m4m4(tmat, mat, scale);
        invert_m4_m4(mat, tmat);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .area_no_sp = area_no_sp, .area_co = area_co, .mat = mat,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_clay_strips_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_clay_strips_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_fill_brush_task_cb(void *userdata, const int n)
+static void do_fill_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -2958,7 +3294,7 @@ static void do_fill_brush_task_cb(void *userdata, const int n)
                                if (plane_trim(ss->cache, brush, val)) {
                                        const float fade = bstrength * tex_strength(
                                                               ss, brush, vd.co, sqrtf(test.dist),
-                                                              vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                                              vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f, thread_id);
 
                                        mul_v3_v3fl(proxy[vd.i], val, fade);
 
@@ -2994,17 +3330,20 @@ static void do_fill_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode)
        mul_v3_fl(temp, displace);
        add_v3_v3(area_co, temp);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .area_no = area_no, .area_co = area_co,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_fill_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_fill_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_scrape_brush_task_cb(void *userdata, const int n)
+static void do_scrape_brush_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -3034,7 +3373,8 @@ static void do_scrape_brush_task_cb(void *userdata, const int n)
 
                                if (plane_trim(ss->cache, brush, val)) {
                                        const float fade = bstrength * tex_strength(
-                                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f,
+                                                              thread_id);
 
                                        mul_v3_v3fl(proxy[vd.i], val, fade);
 
@@ -3070,17 +3410,20 @@ static void do_scrape_brush(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnod
        mul_v3_fl(temp, displace);
        add_v3_v3(area_co, temp);
 
+       set_adaptive_space_factor(sd);
+
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .area_no = area_no, .area_co = area_co,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_scrape_brush_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_scrape_brush_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
-static void do_gravity_task_cb(void *userdata, const int n)
+static void do_gravity_task_cb_ex(
+        void *userdata, void *UNUSED(userdata_chunk), const int n, const int thread_id)
 {
        SculptThreadedTaskData *data = userdata;
        SculptSession *ss = data->ob->sculpt;
@@ -3098,7 +3441,8 @@ static void do_gravity_task_cb(void *userdata, const int n)
        BKE_pbvh_vertex_iter_begin(ss->pbvh, data->nodes[n], vd, PBVH_ITER_UNIQUE) {
                if (sculpt_brush_test_sq(&test, vd.co)) {
                        const float fade = tex_strength(
-                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f);
+                                              ss, brush, vd.co, sqrtf(test.dist), vd.no, vd.fno, vd.mask ? *vd.mask : 0.0f,
+                                              thread_id);
 
                        mul_v3_v3fl(proxy[vd.i], offset, fade);
 
@@ -3123,15 +3467,17 @@ static void do_gravity(Sculpt *sd, Object *ob, PBVHNode **nodes, int totnode, fl
        mul_v3_v3v3(offset, gravity_vector, ss->cache->scale);
        mul_v3_fl(offset, bstrength);
 
+       set_adaptive_space_factor(sd);
+
        /* threaded loop over nodes */
        SculptThreadedTaskData data = {
            .sd = sd, .ob = ob, .brush = brush, .nodes = nodes,
            .offset = offset,
        };
 
-       BLI_task_parallel_range(
-                   0, totnode, &data, do_gravity_task_cb,
-                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+       BLI_task_parallel_range_ex(
+                   0, totnode, &data, NULL, 0, do_gravity_task_cb_ex,
+                   ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT), false);
 }
 
 
@@ -3274,9 +3620,9 @@ static void do_brush_action(Sculpt *sd, Object *ob, Brush *brush, UnifiedPaintSe
 
                BLI_task_parallel_range(
                            0, totnode, &task_data, do_brush_action_task_cb,
-                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
 
-               if (sculpt_brush_needs_normal(brush))
+               if (sculpt_brush_needs_normal(brush, ss->cache->normal_weight))
                        update_sculpt_normal(sd, ob, nodes, totnode);
 
                if (brush->mtex.brush_map_mode == MTEX_MAP_MODE_AREA)
@@ -3455,7 +3801,7 @@ static void sculpt_combine_proxies(Sculpt *sd, Object *ob)
 
                BLI_task_parallel_range(
                            0, totnode, &data, sculpt_combine_proxies_task_cb,
-                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
        }
 
        if (nodes)
@@ -3536,7 +3882,7 @@ static void sculpt_flush_stroke_deform(Sculpt *sd, Object *ob)
 
                BLI_task_parallel_range(
                            0, totnode, &data, sculpt_flush_stroke_deform_task_cb,
-                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_OMP_LIMIT));
+                           ((sd->flags & SCULPT_USE_OPENMP) && totnode > SCULPT_THREADED_LIMIT));
 
                if (vertCos) {
                        sculpt_vertcos_to_key(ob, ss->kb, vertCos);
@@ -3597,6 +3943,10 @@ static void calc_brushdata_symm(Sculpt *sd, StrokeCache *cache, const char symm,
                flip_v3_v3(cache->gravity_direction, cache->true_gravity_direction, symm);
                mul_m4_v3(cache->symm_rot_mat, cache->gravity_direction);
        }
+
+       if (cache->is_rake_rotation_valid) {
+               flip_qt_qt(cache->rake_rotation_symmetry, cache->rake_rotation, symm);
+       }
 }
 
 typedef void (*BrushActionFunc)(Sculpt *sd, Object *ob, Brush *brush, UnifiedPaintSettings *ups);
@@ -3854,7 +4204,9 @@ static void sculpt_init_mirror_clipping(Object *ob, SculptSession *ss)
 }
 
 /* Initialize the stroke cache invariants from operator properties */
-static void sculpt_update_cache_invariants(bContext *C, Sculpt *sd, SculptSession *ss, wmOperator *op, const float mouse[2])
+static void sculpt_update_cache_invariants(
+        bContext *C, Sculpt *sd, SculptSession *ss,
+        wmOperator *op, const float mouse[2])
 {
        StrokeCache *cache = MEM_callocN(sizeof(StrokeCache), "stroke cache");
        Scene *scene = CTX_data_scene(C);
@@ -3899,6 +4251,15 @@ static void sculpt_update_cache_invariants(bContext *C, Sculpt *sd, SculptSessio
        mode = RNA_enum_get(op->ptr, "mode");
        cache->invert = mode == BRUSH_STROKE_INVERT;
        cache->alt_smooth = mode == BRUSH_STROKE_SMOOTH;
+       cache->normal_weight = brush->normal_weight;
+
+       /* interpret invert as following normal, for grab brushes */
+       if (SCULPT_TOOL_HAS_NORMAL_WEIGHT(brush->sculpt_tool)) {
+               if (cache->invert) {
+                       cache->invert = false;
+                       cache->normal_weight = (cache->normal_weight == 0.0f);
+               }
+       }
 
        /* not very nice, but with current events system implementation
         * we can't handle brush appearance inversion hotkey separately (sergey) */
@@ -4041,7 +4402,7 @@ static void sculpt_update_brush_delta(UnifiedPaintSettings *ups, Object *ob, Bru
 
                /* compute 3d coordinate at same z from original location + mouse */
                mul_v3_m4v3(loc, ob->obmat, cache->orig_grab_location);
-               ED_view3d_win_to_3d(cache->vc->ar, loc, mouse, grab_location);
+               ED_view3d_win_to_3d(cache->vc->v3d, cache->vc->ar, loc, mouse, grab_location);
 
                /* compute delta to move verts by */
                if (!cache->first_time) {
@@ -4090,6 +4451,54 @@ static void sculpt_update_brush_delta(UnifiedPaintSettings *ups, Object *ob, Bru
                        copy_v2_v2(ups->anchored_initial_mouse, cache->initial_mouse);
                        ups->anchored_size = ups->pixel_radius;
                }
+
+
+               /* handle 'rake' */
+               cache->is_rake_rotation_valid = false;
+
+               if (cache->first_time) {
+                       copy_v3_v3(cache->rake_data.follow_co, grab_location);
+               }
+
+               if (sculpt_brush_needs_rake_rotation(brush)) {
+                       cache->rake_data.follow_dist = cache->radius * SCULPT_RAKE_BRUSH_FACTOR;
+
+                       if (!is_zero_v3(cache->grab_delta)) {
+                               const float eps = 0.00001f;
+
+                               float v1[3], v2[3];
+
+                               copy_v3_v3(v1, cache->rake_data.follow_co);
+                               copy_v3_v3(v2, cache->rake_data.follow_co);
+                               sub_v3_v3(v2, cache->grab_delta);
+
+                               sub_v3_v3(v1, grab_location);
+                               sub_v3_v3(v2, grab_location);
+
+                               if ((normalize_v3(v2) > eps) &&
+                                   (normalize_v3(v1) > eps) &&
+                                   (len_squared_v3v3(v1, v2) > eps))
+                               {
+                                       const float rake_dist_sq = len_squared_v3v3(cache->rake_data.follow_co, grab_location);
+                                       const float rake_fade = (rake_dist_sq > SQUARE(cache->rake_data.follow_dist)) ?
+                                               1.0f : sqrtf(rake_dist_sq) / cache->rake_data.follow_dist;
+
+                                       float axis[3], angle;
+                                       float tquat[4];
+
+                                       rotation_between_vecs_to_quat(tquat, v1, v2);
+
+                                       /* use axis-angle to scale rotation since the factor may be above 1 */
+                                       quat_to_axis_angle(axis, &angle, tquat);
+                                       normalize_v3(axis);
+
+                                       angle *= brush->rake_factor * rake_fade;
+                                       axis_angle_normalized_to_quat(cache->rake_rotation, axis, angle);
+                                       cache->is_rake_rotation_valid = true;
+                               }
+                       }
+                       sculpt_rake_data_update(&cache->rake_data, grab_location);
+               }
        }
 }
 
@@ -4260,7 +4669,9 @@ static void sculpt_raycast_detail_cb(PBVHNode *node, void *data_v, float *tmin)
        }
 }
 
-static float sculpt_raycast_init(ViewContext *vc, const float mouse[2], float ray_start[3], float ray_end[3], float ray_normal[3], bool original)
+static float sculpt_raycast_init(
+        ViewContext *vc, const float mouse[2],
+        float ray_start[3], float ray_end[3], float ray_normal[3], bool original)
 {
        float obimat[4][4];
        float dist;
@@ -4460,8 +4871,11 @@ static bool sculpt_stroke_test_start(bContext *C, struct wmOperator *op,
                                     const float mouse[2])
 {
        /* Don't start the stroke until mouse goes over the mesh.
-        * note: mouse will only be null when re-executing the saved stroke. */
-       if (!mouse || over_mesh(C, op, mouse[0], mouse[1])) {
+        * note: mouse will only be null when re-executing the saved stroke.
+        * We have exception for 'exec' strokes since they may not set 'mouse', only 'location', see: T52195. */
+       if (((op->flag & OP_IS_INVOKE) == 0) ||
+           (mouse == NULL) || over_mesh(C, op, mouse[0], mouse[1]))
+       {
                Object *ob = CTX_data_active_object(C);
                SculptSession *ss = ob->sculpt;
                Sculpt *sd = CTX_data_tool_settings(C)->sculpt;
@@ -4491,7 +4905,7 @@ static void sculpt_stroke_update_step(bContext *C, struct PaintStroke *UNUSED(st
        sculpt_restore_mesh(sd, ob);
 
        if (sd->flags & SCULPT_DYNTOPO_DETAIL_CONSTANT) {
-               BKE_pbvh_bmesh_detail_size_set(ss->pbvh, sd->constant_detail / 100.0f);
+               BKE_pbvh_bmesh_detail_size_set(ss->pbvh, 1.0f / sd->constant_detail);
        }
        else if (sd->flags & SCULPT_DYNTOPO_DETAIL_BRUSH) {
                BKE_pbvh_bmesh_detail_size_set(ss->pbvh, ss->cache->radius * sd->detail_percent / 100.0f);
@@ -4660,8 +5074,11 @@ static void sculpt_brush_stroke_cancel(bContext *C, wmOperator *op)
        Object *ob = CTX_data_active_object(C);
        SculptSession *ss = ob->sculpt;
        Sculpt *sd = CTX_data_tool_settings(C)->sculpt;
+       const Brush *brush = BKE_paint_brush(&sd->paint);
 
-       if (ss->cache) {
+       /* XXX Canceling strokes that way does not work with dynamic topology, user will have to do real undo for now.
+        *     See T46456. */
+       if (ss->cache && !sculpt_stroke_is_dynamic_topology(ss, brush)) {
                paint_mesh_restore_co(sd, ob);
        }
 
@@ -4730,117 +5147,3026 @@ static void SCULPT_OT_set_persistent_base(wmOperatorType *ot)
        ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
 }
 
-/************************** Dynamic Topology **************************/
+/****************** Topology tools Silhouette ******************/
 
-static void sculpt_dynamic_topology_triangulate(BMesh *bm)
+/* Axis-aligned bounding box (Same as in pbvh_intern.h)*/
+
+/* TODO: import BB functions */
+/* Expand the bounding box to include a new coordinate */
+static void BB_expand(BB *bb, const float co[3])
 {
-       if (bm->totloop != bm->totface * 3) {
-               BM_mesh_triangulate(bm, MOD_TRIANGULATE_QUAD_BEAUTY, MOD_TRIANGULATE_NGON_EARCLIP, false, NULL, NULL, NULL);
+       for (int i = 0; i < 3; ++i) {
+               bb->bmin[i] = min_ff(bb->bmin[i], co[i]);
+               bb->bmax[i] = max_ff(bb->bmax[i], co[i]);
        }
 }
 
-void sculpt_pbvh_clear(Object *ob)
+static void BB_reset(BB *bb)
 {
-       SculptSession *ss = ob->sculpt;
-       DerivedMesh *dm = ob->derivedFinal;
+       bb->bmin[0] = bb->bmin[1] = bb->bmin[2] = FLT_MAX;
+       bb->bmax[0] = bb->bmax[1] = bb->bmax[2] = -FLT_MAX;
+}
 
-       /* Clear out any existing DM and PBVH */
-       if (ss->pbvh)
-               BKE_pbvh_free(ss->pbvh);
-       ss->pbvh = NULL;
-       if (dm)
-               dm->getPBVH(NULL, dm);
-       BKE_object_free_derived_caches(ob);
+static bool bb_intersect(BB *bb1, BB *bb2) {
+       int i;
+
+       /* min is inclusive max is exclusive? BB*/
+       for (i = 0; i < 3; ++i) {
+               if(bb1->bmin[i] >= bb2->bmax[i] || bb1->bmax[i] < bb2->bmin[i]){
+                       return false;
+               }
+       }
+
+       return true;
 }
 
-void sculpt_dyntopo_node_layers_add(SculptSession *ss)
+static void silhouette_stroke_free(SilhouetteStroke *stroke)
 {
-       int cd_node_layer_index;
+       if (stroke) {
+               if (stroke->points) {
+                       MEM_SAFE_FREE(stroke->points);
+               }
+               MEM_SAFE_FREE(stroke);
+       }
+}
 
-       char layer_id[] = "_dyntopo_node_id";
+static SilhouetteStroke *silhouette_stroke_new()
+{
+       SilhouetteStroke *stroke = MEM_callocN(sizeof(SilhouetteStroke), "SilhouetteStroke");
+       stroke->points = MEM_callocN(sizeof(float) * 3 * SIL_STROKE_STORE_CHUNK,"SilhouetteStrokePoints");
+       stroke->points_v2 = MEM_callocN(sizeof(float) * 2 * SIL_STROKE_STORE_CHUNK,"SilhouetteStrokePoints");
+       stroke->totvert = 0;
+       stroke->max_verts = SIL_STROKE_STORE_CHUNK;
+       BB_reset(&stroke->bb);
+       return stroke;
+}
 
-       cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->vdata, CD_PROP_INT, layer_id);
-       if (cd_node_layer_index == -1) {
-               BM_data_layer_add_named(ss->bm, &ss->bm->vdata, CD_PROP_INT, layer_id);
-               cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->vdata, CD_PROP_INT, layer_id);
-       }
+static SilhouetteData *silhouette_data_new(bContext *C)
+{
+       SilhouetteData *sil = MEM_callocN(sizeof(SilhouetteData), "SilhouetteData");
+       Object *obedit = CTX_data_edit_object(C);
+       Scene *scene = CTX_data_scene(C);
+       Sculpt *sd = CTX_data_tool_settings(C)->sculpt;
+       View3D *v3d = CTX_wm_view3d(C);
+       const float *fp = ED_view3d_cursor3d_get(scene, v3d);
 
-       ss->cd_vert_node_offset = CustomData_get_n_offset(&ss->bm->vdata, CD_PROP_INT,
-                                                         cd_node_layer_index - CustomData_get_layer_index(&ss->bm->vdata, CD_PROP_INT));
+       sil->ar = CTX_wm_region(C);
+       sil->current_stroke = silhouette_stroke_new();
+       view3d_set_viewcontext(C, &sil->vc);
 
-       ss->bm->vdata.layers[cd_node_layer_index].flag |= CD_FLAG_TEMPORARY;
+       sil->add_col[0] = 1.00; /* add mode color is light red */
+       sil->add_col[1] = 0.39;
+       sil->add_col[2] = 0.39;
 
-       cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->pdata, CD_PROP_INT, layer_id);
-       if (cd_node_layer_index == -1) {
-               BM_data_layer_add_named(ss->bm, &ss->bm->pdata, CD_PROP_INT, layer_id);
-               cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->pdata, CD_PROP_INT, layer_id);
-       }
+       /*Load RNA Data if present */
+       sil->smoothness = sd->silhouette_smoothness / 100.0f;
+       sil->depth = sd->silhouette_depth;
+       sil->resolution = sd->silhouette_resolution;
 
-       ss->cd_face_node_offset = CustomData_get_n_offset(&ss->bm->pdata, CD_PROP_INT,
-                                                         cd_node_layer_index - CustomData_get_layer_index(&ss->bm->pdata, CD_PROP_INT));
+       copy_v3_v3(sil->anchor, fp);
 
-       ss->bm->pdata.layers[cd_node_layer_index].flag |= CD_FLAG_TEMPORARY;
-}
+       /* Intersection variables */
+       sil->fillet_ring_orig = NULL;
+       sil->fillet_ring_orig_start = NULL;
+       sil->fillet_ring_new = NULL;
+       sil->fillet_ring_new_start = NULL;
+       sil->inter_edges = NULL;
+       sil->fillet_ring_bbs = NULL;
 
+       sil->scene = scene;
+       sil->ob = obedit;
+       sil->state = SIL_INIT;
+       return sil;
+}
 
-void sculpt_update_after_dynamic_topology_toggle(bContext *C)
+static void silhouette_data_free(struct wmOperator *op)
 {
-       Scene *scene = CTX_data_scene(C);
-       Object *ob = CTX_data_active_object(C);
-       Sculpt *sd = scene->toolsettings->sculpt;
+       SilhouetteData *data;
+       data = op->customdata;
+       if (data) {
+               silhouette_stroke_free(data->current_stroke);
 
-       /* Create the PBVH */
-       BKE_sculpt_update_mesh_elements(scene, sd, ob, false, false);
-       WM_event_add_notifier(C, NC_OBJECT | ND_DRAW, ob);
+               if (data->inter_edges) {
+                       MEM_freeN(data->inter_edges);
+               }
+               if (data->fillet_ring_orig) {
+                       MEM_freeN(data->fillet_ring_orig);
+               }
+               if (data->fillet_ring_orig_start) {
+                       MEM_freeN(data->fillet_ring_orig_start);
+               }
+               if (data->fillet_ring_new) {
+                       MEM_freeN(data->fillet_ring_new);
+               }
+               if (data->fillet_ring_new_start) {
+                       MEM_freeN(data->fillet_ring_new_start);
+               }
+               if (data->fillet_ring_bbs) {
+                       MEM_freeN(data->fillet_ring_bbs);
+               }
+               MEM_SAFE_FREE(data);
+       }
 }
 
-void sculpt_dynamic_topology_enable(bContext *C)
+static void silhoute_stroke_point_to_3d(SilhouetteData *sil, int point, float r_v[3])
 {
-       Scene *scene = CTX_data_scene(C);
-       Object *ob = CTX_data_active_object(C);
-       SculptSession *ss = ob->sculpt;
-       Mesh *me = ob->data;
-       const BMAllocTemplate allocsize = BMALLOC_TEMPLATE_FROM_ME(me);
+       copy_v3_v3(r_v, &sil->current_stroke->points[point]);
+       /*ED_view3d_win_to_3d(sil->vc.v3d, sil->ar, sil->anchor, &sil->current_stroke->points[point], r_v);*/
+}
 
-       sculpt_pbvh_clear(ob);
+#if 0
+/* TODO: Add dynamic memory allocation */
+static void silhouette_stroke_add_3Dpoint(SilhouetteStroke *stroke, float point[3])
+{
+       if (stroke->totvert >= stroke->max_verts) {
+               stroke->max_verts += SIL_STROKE_STORE_CHUNK;
+               stroke->points = MEM_reallocN(stroke->points, sizeof(float) * 3 * stroke->max_verts);
+               stroke->points_v2 = MEM_reallocN(stroke->points_v2, sizeof(float) * 2 * stroke->max_verts);
+       }
 
-       ss->bm_smooth_shading = (scene->toolsettings->sculpt->flags & SCULPT_DYNTOPO_SMOOTH_SHADING) != 0;
+       copy_v3_v3(&stroke->points[stroke->totvert * 3], point);
+       stroke->totvert ++;
+}
+#endif
 
-       /* Dynamic topology doesn't ensure selection state is valid, so remove [#36280] */
-       BKE_mesh_mselect_clear(me);
+static void silhouette_stroke_add_point(SilhouetteData *sil, SilhouetteStroke *stroke, float point[2])
+{
+       if (stroke->totvert >= stroke->max_verts) {
+               stroke->max_verts += SIL_STROKE_STORE_CHUNK;
+               stroke->points = MEM_reallocN(stroke->points, sizeof(float) * 3 * stroke->max_verts);
+               stroke->points_v2 = MEM_reallocN(stroke->points_v2, sizeof(float) * 2 * stroke->max_verts);
+       }
+
+       copy_v2_v2(&stroke->points_v2[stroke->totvert * 2], point);
+       ED_view3d_win_to_3d(sil->vc.v3d, sil->ar, sil->anchor, point, &sil->current_stroke->points[stroke->totvert * 3]);
+       BB_expand(&stroke->bb, &sil->current_stroke->points[stroke->totvert * 3]);
+       stroke->totvert ++;
+}
+
+/* Set reference plane, 3D plane which is drawn on in 2D */
+static void silhouette_set_ref_plane(SilhouetteData *sil)
+{
+       /*Get the view vector. Same as ED_view3d_global_to_vector from the center of screen or orthographic*/
+       negate_v3_v3(sil->z_vec, sil->vc.rv3d->viewinv[2]);
+       normalize_v3(sil->z_vec);
+       /*ED_view3d_global_to_vector(sil->ar->regiondata, sil->anchor, sil->z_vec);*/
+}
+
+static void sculpt_silhouette_stroke_update(float mouse[2], SilhouetteData *sil)
+{
+       float anchor[3];
+       silhouette_set_ref_plane(sil);
+       SilhouetteStroke *stroke = sil->current_stroke;
+
+       sil->last_mouse_pos[0] = mouse[0];
+       sil->last_mouse_pos[1] = mouse[1];
+       silhouette_stroke_add_point(sil, stroke, sil->last_mouse_pos);
+       interp_v3_v3v3(anchor, sil->anchor, &stroke->points[stroke->totvert * 3 - 3], 1 / stroke->totvert);
+       copy_v3_v3(sil->anchor, anchor);
+       ED_region_tag_redraw(sil->ar);
+       copy_v2_v2(sil->last_mouse_pos,mouse);
+}
+
+typedef enum {
+       BRANCH_INIT = 0,
+       BRANCH_EDGE_GEN = 1     /* Edges on the ends are stored (used primarly for bridging) */
+} BranchState;
+
+
+/* Spine Memory 
+ * The Spine is a treelike representation of the drawn shape.
+ * It is used to determine the topology of the shape.
+ * The Spine is generated in 4 steps
+ * 1. Triangulate the hull drawn by the sculptor. This uses a BMesh triangulation method
+ * 2. Traverse all triangles, starting from the first, marching from one tri to its neighbouring tris.
+ * 3. Every step gets stored into a branch. Two new branches are created when a tri with 3 neighbours is hit. Traversing ends when a tri with one neighbour is hit.
+ * 4. Small, stubby, branches are dissolved into the parent branches.
+ * This results in a tree with 3 types of branches. One end, the Caps. Two ends, the Tubes. Three ends, the t-intersections.*/
+
+typedef struct SpineBranch{
+       int totpoints;                  /* Number of points of the spine. */
+       int totforks;                   /* Number of connected adjacent branches */
+       float *points;                  /* Points of the spine. Generated from tri-centers */
+       int tot_hull_points;
+       int *hull_points;               /* Hullpoints, pointer to stroke points */
+       int idx;                                /* Index in the Spine branches array */
+       int *terminal_points;   /* Description of the connected branches. Per fork 2 ints (point,branch_idx) */
+       BranchState flag;
+       int *e_start_arr;               /* Edges on the ends are stored (used primarly for bridging) */
+       int fs_bs_offset;               /* Frontside edge offset to backside*/
+       int *e_flip_side_ends;  /* Front and backside connecting edges of each part*/
+       bool intersecting;
+}SpineBranch;
+
+/* Main Tree Container */
+typedef struct Spine{
+       int totbranches;
+       SpineBranch **branches; /* All branches. Can contain Null pointers if branches got removed*/
+}Spine;
 
-       /* Create triangles-only BMesh */
-       ss->bm = BM_mesh_create(&allocsize);
+typedef struct {
+       SculptSession *ss;
+       BB *bb_target;
+       bool original;
+} SculptSearchBBData;
 
-       BM_mesh_bm_from_me(ss->bm, me, true, true, ob->shapenr);
-       sculpt_dynamic_topology_triangulate(ss->bm);
-       BM_data_layer_add(ss->bm, &ss->bm->vdata, CD_PAINT_MASK);
-       sculpt_dyntopo_node_layers_add(ss);
-       /* make sure the data for existing faces are initialized */
-       BM_mesh_normals_update(ss->bm);
+/* Search nodes to integrate another BB into (used for silhouette)*/
+static bool sculpt_search_BB_cb(PBVHNode *node, void *data_v)
+{
+       SculptSearchBBData *data = data_v;
+       BB *bb_target = data->bb_target;
+       float bb_min[3], bb_max[3];
+       int i;
 
-       /* Enable dynamic topology */
-       me->flag |= ME_SCULPT_DYNAMIC_TOPOLOGY;
-       
-       /* Enable logging for undo/redo */
-       ss->bm_log = BM_log_create(ss->bm);
+       BKE_pbvh_node_get_BB(node, bb_min, bb_max);
 
-       /* Refresh */
-       sculpt_update_after_dynamic_topology_toggle(C);
+       /* min is inclusive max is exclusive? BB*/
+       for (i = 0; i < 3; ++i) {
+               if(bb_target->bmin[i] >= bb_max[i] || bb_target->bmax[i] < bb_min[i]){
+                       return false;
+               }
+       }
+       return true;
 }
 
-/* Free the sculpt BMesh and BMLog
- *
- * If 'unode' is given, the BMesh's data is copied out to the unode
- * before the BMesh is deleted so that it can be restored from */
-void sculpt_dynamic_topology_disable(bContext *C,
-                                     SculptUndoNode *unode)
+static int get_adjacent_faces(BMFace *f, BMFace **ad_f, BMFace *last_f)
 {
-       Object *ob = CTX_data_active_object(C);
-       SculptSession *ss = ob->sculpt;
-       Mesh *me = ob->data;
+       int ad_faces = 0;
+       /* There should only be tris in a triangulated mesh. */
+       BLI_assert(f->len == 3);
 
-       sculpt_pbvh_clear(ob);
+       /* Loop edges in faces */
+       BMLoop *f_t_l = f->l_first;
+       for (int i = 0; i < 3; i++) {
+               /* Loop faces connected to this edge */
+               BMLoop *t_l = f_t_l;
+               while (t_l != f_t_l->radial_prev && ad_faces < 3) {
+                       if (t_l->f != f && t_l->f != last_f) {
+                               ad_f[ad_faces] = t_l->f;
+                               ad_faces++;
+                       }
+                       t_l = t_l->radial_next;
+               }
+               if (t_l->f != f && t_l->f != last_f) {
+                       ad_f[ad_faces] = t_l->f;
+                       ad_faces++;
+               }
+               f_t_l = f_t_l->next;
+       }
+       return ad_faces;
+}
+
+#ifdef DEBUG_DRAW
+static void debug_branch(SpineBranch *sb, unsigned int color)
+{
+       float v1[3], v2[3];
+
+       if (sb->totpoints > 1) {
+               bl_debug_color_set(color);
+               for (int j = 1; j < sb->totpoints; j++) {
+                       copy_v3_v3(v1,&sb->points[j*3-3]);
+                       copy_v3_v3(v2,&sb->points[j*3  ]);
+                       bl_debug_draw_edge_add(v1,v2);
+               }
+       }
+       bl_debug_color_set(color);
+
+       /*for(int j = 0; j < sb->tot_hull_points; j++){
+               silhoute_stroke_point_to_3d(sil, sb->hull_points[j]*2, v1);
+               bl_debug_draw_point(v1, 0.2f);
+        }*/
+
+       /*SpineBranch *rsb;*/
+       /*unsigned int forkcol = 0xFFFF00;
+        for(int j = 1; j < sb->totforks; j++){
+               printf("Branch fork %i, stroke ids: %i, %i\n", j, sb->terminal_points[j*4], sb->terminal_points[j*4+3]);
+               if (spine->branches[sb->terminal_points[j*4+2]] && spine->branches[sb->terminal_points[j*4+2]]->totforks > 0) {
+        bl_debug_color_set(forkcol);
+        silhoute_stroke_point_to_3d(sil, sb->terminal_points[j*4]*2, v1);
+        bl_debug_draw_point(v1, 0.2f);
+        silhoute_stroke_point_to_3d(sil, sb->terminal_points[j*4+3]*2, v1);
+        bl_debug_draw_point(v1, 0.2f);
+        bl_debug_color_set(color);
+               }
+        }*/
+}
+
+static void debug_spine(Spine *spine)
+{
+       char color = 0;
+
+       for (int i = 0; i < spine->totbranches; i++) {
+               if (spine->branches[i]) {
+                       bl_debug_color_set(color);
+                       SpineBranch *sb = spine->branches[i];
+                       debug_branch(sb, color);
+                       color += 90;
+               }
+       }
+
+       bl_debug_color_set(0x000000);
+}
+#endif
+
+static void free_spine_branch(SpineBranch *branch)
+{
+       if (branch->points) {
+               MEM_freeN(branch->points);
+       }
+       if (branch->hull_points) {
+               MEM_freeN(branch->hull_points);
+       }
+       if (branch->terminal_points) {
+               MEM_freeN(branch->terminal_points);
+       }
+       if (branch->flag & BRANCH_EDGE_GEN) {
+               MEM_freeN(branch->e_start_arr);
+               MEM_freeN(branch->e_flip_side_ends);
+       }
+       MEM_freeN(branch);
+}
+
+static void detach_branch(SpineBranch *b, SpineBranch *db)
+{
+       bool clear = false;
+       /* Find the branch to be disconnected in the fork array and shift all following forward */
+       for (int i = 0; i < b->totforks; i++) {
+               if (b->terminal_points[i * 2 + 1] == db->idx) {
+                       clear = true;
+               } else if (clear) {
+                       b->terminal_points[i * 2 - 2] = b->terminal_points[i * 2    ];
+                       b->terminal_points[i * 2 - 1] = b->terminal_points[i * 2 + 1];
+               }
+       }
+       if (clear) {
+               b->totforks --;
+       }
+}
+
+static void dissolve_branch(Spine *spine, SpineBranch *branch, SpineBranch *t_branch)
+{
+       /* Dissolve all connected branches recursively which arent from the target subbranch */
+       for (int i = 0; i < branch->totforks; i++) {
+               if (branch->terminal_points[i * 2 + 1] != t_branch->idx) {
+                       dissolve_branch(spine, spine->branches[branch->terminal_points[i * 2 + 1]], branch);
+               }
+       }
+
+       detach_branch(t_branch, branch);
+
+       /* Copy Hullpoints from the dissolving branch to the target branch */
+       for (int i = 0; i < branch->tot_hull_points; i++) {
+               t_branch->hull_points[t_branch->tot_hull_points] = branch->hull_points[i];
+               t_branch->tot_hull_points ++;
+       }
+
+       spine->branches[branch->idx] = NULL;
+       free_spine_branch(branch);
+}
+
+static SpineBranch *new_spine_branch(int idx, int max_alloc, int hull_max)
+{
+       SpineBranch *branch = MEM_callocN(sizeof(SpineBranch), "Spine Branch");
+       branch->totpoints = 0;
+       branch->totforks = 0;
+
+       /*TODO: way to big, maybe shrink if done creating or dynamic arrays?*/
+       branch->points = MEM_callocN(sizeof(float) * 3 * max_alloc, __func__);
+       branch->hull_points = MEM_callocN(sizeof(int) * hull_max * 2 * 3, "Spine Hull");
+       branch->terminal_points = MEM_callocN(sizeof(int) * 3 * 2, "Spine terminal");
+
+       branch->idx = idx;
+       return branch;
+}
+
+static Spine *new_spine(int max_alloc, int hull_max)
+{
+       Spine *s = MEM_callocN(sizeof(Spine),"Silhouette Spine");
+       s->branches = MEM_callocN(sizeof(SpineBranch *) * max_alloc,"branches");
+       s->totbranches = 1;
+       s->branches[0] = new_spine_branch(0, max_alloc, hull_max);
+       return s;
+}
+
+static void free_spine(Spine *spine)
+{
+       for (int i = 0; i < spine->totbranches; i++) {
+               if (spine->branches[i]) {
+                       free_spine_branch(spine->branches[i]);
+               }
+       }
+       MEM_freeN(spine->branches);
+}
+
+static SpineBranch *spine_branchoff(Spine *spine, SpineBranch *current_branch, int max_alloc, int hull_max)
+{
+       SpineBranch *new_branch = new_spine_branch(spine->totbranches, max_alloc, hull_max);
+       spine->branches[spine->totbranches] = new_branch;
+
+       current_branch->terminal_points[current_branch->totforks * 2    ] = current_branch->totpoints - 1;
+       current_branch->terminal_points[current_branch->totforks * 2 + 1] = spine->totbranches;
+
+       copy_v3_v3(&new_branch->points[0], &current_branch->points[current_branch->totpoints * 3 - 3]);
+       new_branch->totpoints = 1;
+
+       new_branch->terminal_points[0] = 0;
+       new_branch->terminal_points[1] = current_branch->idx;
+
+       current_branch->totforks ++;
+       new_branch->totforks = 1;
+
+       spine->totbranches ++;
+
+       return new_branch;
+}
+
+static void add_face_to_branch(SpineBranch *branch, BMFace *f)
+{
+       float center[3];
+       BM_face_calc_center_mean_weighted(f, center);
+       copy_v3_v3(branch->points + branch->totpoints * 3, center);
+       branch->totpoints++;
+
+       branch->hull_points[branch->tot_hull_points    ] = BM_elem_index_get(f->l_first->v);
+       branch->hull_points[branch->tot_hull_points + 1] = BM_elem_index_get(f->l_first->next->v);
+       branch->hull_points[branch->tot_hull_points + 2] = BM_elem_index_get(f->l_first->prev->v);
+       branch->tot_hull_points += 3;
+}
+
+static int calc_mid_spine_rec(BMFace *f, Spine *spine, SpineBranch *active_branch, BMFace *last_f, int max_alloc, int hull_max)
+{
+       BMFace **ad_f = MEM_callocN(sizeof(BMFace *) * 3,__func__);
+       int adjacent_faces = get_adjacent_faces(f, ad_f, last_f);
+       int added_points = 0;           /* Points added by the current branch */
+       int sub_added_points = 0;       /* Points added by child branches */
+
+       /* TODO Maybe not duplicates? */
+       add_face_to_branch(active_branch, f);
+       added_points ++;
+
+       SpineBranch *new_branch;
+
+       if (adjacent_faces == 1) {
+               added_points += calc_mid_spine_rec(ad_f[0], spine, active_branch, f, max_alloc, hull_max);
+       } else if (adjacent_faces == 2) {
+               /*
+                Longer Branches?
+
+                if(last_f != NULL){
+                       float v_tmp[3], v_origin[3], v_branch1[3], v_branch2[3];
+                       BM_face_calc_center_mean_weighted(f, v_tmp);
+                       BM_face_calc_center_mean_weighted(last_f, v_origin);
+                       sub_v3_v3(v_origin, v_tmp);
+                       BM_face_calc_center_mean_weighted(ad_f[0], v_branch1);
+                       sub_v3_v3(v_branch1, v_tmp);
+                       normalize_v3(v_branch1);
+                       BM_face_calc_center_mean_weighted(ad_f[1], v_branch2);
+                       sub_v3_v3(v_branch2, v_tmp);
+                       normalize_v3(v_branch2);
+
+                       //switch branches depending on angle
+                       if(dot_v3v3(v_origin, v_branch1) < dot_v3v3(v_origin, v_branch2)){
+                               ad_f[2] = ad_f[0];
+                               ad_f[0] = ad_f[1];
+                               ad_f[1] = ad_f[2];
+                       }
+               }*/
+
+               new_branch = spine_branchoff(spine, active_branch, max_alloc, hull_max);
+               sub_added_points = calc_mid_spine_rec(ad_f[0], spine, new_branch, f, max_alloc, hull_max);
+               /* Controls when to remove small/stubby branches */
+               if (sub_added_points < 20 && new_branch->totforks < 3) {
+                       dissolve_branch(spine, new_branch, active_branch);
+               }
+               added_points += sub_added_points;
+
+               new_branch = spine_branchoff(spine, active_branch, max_alloc, hull_max);
+               sub_added_points = calc_mid_spine_rec(ad_f[1], spine, new_branch, f, max_alloc, hull_max);
+               if (sub_added_points < 20 && new_branch->totforks < 3) {
+                       dissolve_branch(spine, new_branch, active_branch);
+               }
+               added_points += sub_added_points;
+       }
+
+       MEM_freeN(ad_f);
+       return added_points;
+}
+
+static Spine *silhouette_generate_spine(SilhouetteData *sil, SilhouetteStroke *stroke)
+{
+       BMesh *bm;
+       const struct BMeshCreateParams bm_create_params = {0};
+       BMVert **vert_arr = MEM_mallocN(sizeof(BMVert *) * stroke->totvert, __func__);
+       BMFace *f;
+       BMVert *v;
+       float v_co[3];
+
+       /* Generate a BMesh from the drawn hull */
+       bm = BM_mesh_create(&bm_mesh_allocsize_default,
+                                               &bm_create_params);
+
+       for (int i = 0; i < stroke->totvert; i++) {
+               silhoute_stroke_point_to_3d(sil, i * 3, v_co);
+               v = BM_vert_create( bm, v_co, NULL, BM_CREATE_NOP);
+               BM_elem_index_set(v,i);
+               vert_arr[i] = v;
+       }
+
+       f = BM_face_create_ngon_verts(bm, vert_arr, stroke->totvert,
+                                                                 NULL, BM_CREATE_NO_DOUBLE,
+                                                                 true, true);
+       BM_face_normal_update(f);
+
+       int save_max = (2 * stroke->totvert - 2) * 2;
+       int faces_array_tot = save_max;/* Upper limit -convexhull verts for precise calc */
+       BMFace **faces_array = BLI_array_alloca(faces_array, faces_array_tot);
+       struct MemArena *pf_arena;
+       struct Heap *pf_heap;
+       struct EdgeHash *pf_ehash;
+
+       pf_arena = BLI_memarena_new(BLI_POLYFILL_ARENA_SIZE, __func__);
+
+       pf_heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
+       pf_ehash = BLI_edgehash_new_ex(__func__, BLI_POLYFILL_ALLOC_NGON_RESERVE);
+
+       /* Triangulate to traverse faces */
+       BM_face_triangulate(bm, f,
+                                               faces_array,
+                                               &faces_array_tot,
+                                               NULL,
+                                               NULL,
+                                               NULL,
+                                               MOD_TRIANGULATE_QUAD_BEAUTY,
+                                               MOD_TRIANGULATE_NGON_BEAUTY,
+                                               false, pf_arena, pf_heap, pf_ehash);
+
+       /*
+        Debug triangulated face:
+        for (int i = 0; i < faces_array_tot; i++){
+               f = faces_array[i];
+               bl_debug_draw_edge_add(f->l_first->e->v1->co,f->l_first->e->v2->co);
+               bl_debug_draw_edge_add(f->l_first->next->e->v1->co,f->l_first->next->e->v2->co);
+               bl_debug_draw_edge_add(f->l_first->next->next->e->v1->co,f->l_first->next->next->e->v2->co);
+       }*/
+
+       /* start traversing at the first face */
+       f = faces_array[0];
+       Spine *spine = new_spine(save_max, stroke->totvert);
+       printf("Verts in stroke: %i\n", stroke->totvert);
+       /* traverse recursively */
+       calc_mid_spine_rec(f, spine, spine->branches[0], NULL, save_max, stroke->totvert);
+       printf("Spine generated with %i Branches.\n", spine->totbranches);
+#ifdef DEBUG_DRAW
+       /*debug_spine(spine);*/
+#endif
+       return spine;
+}
+
+/*
+ *  d      ** | ******c
+ *      *     |       |
+ *    *       |       |
+ *   *        |       |
+ *  *         |       |
+ *  ---------- --------
+ *  *         |       |
+ *  *      smooth     |
+ * a _________|______ b
+ * Interpolate between the three points resulting in a vertex line between a and c.
+ * Smoothness regulates the cutoff to start a circular interpolation
+ */
+static void calc_vert_quarter(Mesh *me, float a[3], float b[3], float c[3], int v_steps, int w_h_steps, float smoothness, bool flip, bool flip_side){
+       int v_start = me->totvert;
+       int v_pos = flip ? v_start + v_steps + w_h_steps - 1 - (flip_side ? 1 : 0): v_start;
+       float inv_smooth = 1.0f - smoothness;
+       float v1[3], v2[3], v3[3], v4[3], v5[3], up[3], side[3], d[3];
+       float f_sin, f_cos;
+       int s_steps_w = inv_smooth * w_h_steps;
+       int s_steps_v = fmax(1, inv_smooth * v_steps);
+       int s_steps_c = v_steps - s_steps_v + w_h_steps - s_steps_w;
+
+       ED_mesh_vertices_add(me, NULL, v_steps + w_h_steps - (flip_side ? 1 : 0));
+
+       sub_v3_v3v3(up, c, b);
+       add_v3_v3v3(d, a, up);
+       mul_v3_fl(up, 1.0f / (float)v_steps);
+       sub_v3_v3v3(side, a, b);
+       if (w_h_steps > 0) {
+               mul_v3_fl(side, 1.0f / (float)(w_h_steps));
+       }
+       mul_v3_v3fl(v2, side, s_steps_w);
+       add_v3_v3(v2, c);
+
+       copy_v3_v3(v1, a);
+       for (int v = 0; v < s_steps_v; v++) {
+               if (!flip_side || v > 0) {
+                       copy_v3_v3(me->mvert[v_pos].co, v1);
+                       me->mvert[v_pos].flag = 0;
+                       me->mvert[v_pos].bweight = 0;
+                       v_pos += flip ? - 1 : 1;
+               }
+               if (v < s_steps_v - 1) {
+                       add_v3_v3(v1, up);
+               }
+       }
+
+       sub_v3_v3v3(v3, d, v1);
+       sub_v3_v3v3(v4, d, v2);
+       sub_v3_v3(v1, v4);
+
+       for (int c = 1; c < s_steps_c + 1; c++) {
+               f_sin = sin(((float)c / (float)(s_steps_c + 1)) * M_PI / 2);
+               f_cos = cos(((float)c / (float)(s_steps_c + 1)) * M_PI / 2);
+
+               mul_v3_v3fl(v2, v4, f_cos);
+               add_v3_v3v3(v5, v2, v1);
+               mul_v3_v3fl(v2, v3, f_sin);
+               add_v3_v3(v5, v2);
+
+               copy_v3_v3(me->mvert[v_pos].co, v5);
+               me->mvert[v_pos].flag = 0;
+               me->mvert[v_pos].bweight = 0;
+               v_pos += flip ? - 1 : 1;
+       }
+
+       add_v3_v3(v1, v3);
+       for (int w = 0; w < s_steps_w; w++) {
+               copy_v3_v3(me->mvert[v_pos].co, v1);
+               me->mvert[v_pos].flag = 0;
+               me->mvert[v_pos].bweight = 0;
+               sub_v3_v3(v1, side);
+               v_pos += flip ? - 1 : 1;
+       }
+}
+
+/* W_steps needs to be uneven! */
+static void calc_vert_half(Mesh *me, float left[3], float right[3], float center[3], float upper_center[3], int v_steps, int w_steps, float smoothness, bool flip_side)
+{
+       int half_w = w_steps / 2;
+       calc_vert_quarter(me, left, center, upper_center, v_steps, half_w, smoothness, false, flip_side);
+       ED_mesh_vertices_add(me, NULL, 1);
+       copy_v3_v3(me->mvert[me->totvert - 1].co, upper_center);
+       me->mvert[me->totvert - 1].flag = 0;
+       me->mvert[me->totvert - 1].bweight = 0;
+       calc_vert_quarter(me, right, center, upper_center, v_steps, half_w, smoothness, true, flip_side);
+}
+
+/* Create standardised Mesh. Simple UxV rectangular grid. (Edges, Loops, Polys):*/
+static void generate_mesh_grid_f_e(Mesh *me, int u_steps, int v_steps, int v_start, bool n_flip)
+{
+       int face_count = (u_steps - 1) * (v_steps - 1);
+       int e_start = me->totedge;
+       int l_start = me->totloop;
+       int p_start = me->totpoly;
+       int p_pos = 0, l_pos = 0, e_pos = 0, v_pos = 0;
+       int edges_per_side = face_count * 2 + u_steps + v_steps - 2;
+       /*TODO: drag out of branch loop*/
+       ED_mesh_edges_add(me, NULL, edges_per_side);
+       ED_mesh_loops_add(me, NULL, face_count * 4);
+       ED_mesh_polys_add(me, NULL, face_count);
+
+       for (int u = 0; u < u_steps; u++) {
+               for (int v = 0; v < v_steps; v++) {
+                       /* add edges */
+                       if (v < v_steps - 1 && u < u_steps - 1) {
+                               e_pos = e_start + u * (v_steps * 2 - 1) + v * 2;
+                               me->medge[e_pos].crease = 0;
+                               me->medge[e_pos].bweight = 0;
+                               me->medge[e_pos].flag = 0;
+                               me->medge[e_pos].v1 = v_start + u * v_steps + v;
+                               me->medge[e_pos].v2 = v_start + u * v_steps + v + 1;
+
+                               me->medge[e_pos + 1].crease = 0;
+                               me->medge[e_pos + 1].bweight = 0;
+                               me->medge[e_pos + 1].flag = 0;
+                               me->medge[e_pos + 1].v1 = v_start + u * v_steps + v;
+                               me->medge[e_pos + 1].v2 = v_start + (u + 1) * v_steps + v;
+
+                               /* add loops */
+                               l_pos = l_start + u * 4 * (v_steps - 1) + v * 4;
+                               v_pos = v_start + u * v_steps + v;
+                               me->mloop[l_pos].v = v_pos;
+                               me->mloop[l_pos].e = e_pos;
+
+                               if (n_flip) {
+                                       /* clockwise */
+                                       me->mloop[l_pos + 1].v = v_pos + 1;
+                                       me->mloop[l_pos + 1].e = e_pos + 3;
+
+                                       me->mloop[l_pos + 2].v = v_pos + v_steps + 1;
+                                       me->mloop[l_pos + 2].e = e_pos + v_steps * 2 - 1;
+
+                                       me->mloop[l_pos + 3].v = v_pos + v_steps;
+                                       me->mloop[l_pos + 3].e = e_pos + 1;
+                               } else {
+                                       /* anti clockwise */
+                                       me->mloop[l_pos + 1].v = v_pos + v_steps;
+                                       me->mloop[l_pos + 1].e = e_pos + 1;
+
+                                       me->mloop[l_pos + 2].v = v_pos + v_steps + 1;
+                                       me->mloop[l_pos + 2].e = e_pos + v_steps * 2 - 1;
+
+                                       me->mloop[l_pos + 3].v = v_pos + 1;
+                                       me->mloop[l_pos + 3].e = e_pos + 3;
+                               }
+
+                               /* add Polys */
+                               p_pos = p_start + u * (v_steps - 1) + v;
+                               me->mpoly[p_pos].totloop = 4;
+                               me->mpoly[p_pos].loopstart = l_pos;
+                               me->mpoly[p_pos].mat_nr = 0;
+                               me->mpoly[p_pos].flag = 0;
+                               me->mpoly[p_pos].pad = 0;
+
+
+                       } else if (v == v_steps - 1 && u != u_steps - 1){
+                               e_pos = e_start + u * (v_steps * 2 - 1) + v * 2;
+                               me->medge[e_pos].crease = 0;
+                               me->medge[e_pos].bweight = 0;
+                               me->medge[e_pos].flag = 0;
+                               me->medge[e_pos].v1 = v_start + u * v_steps + v;
+                               me->medge[e_pos].v2 = v_start + (u + 1) * v_steps + v;
+                       } else if (u == u_steps - 1 && v != v_steps - 1) {
+                               e_pos = e_start + u * (v_steps * 2 - 1) + v;
+                               me->medge[e_pos].crease = 0;
+                               me->medge[e_pos].bweight = 0;
+                               me->medge[e_pos].flag = 0;
+                               me->medge[e_pos].v1 = v_start + u * v_steps + v;
+                               me->medge[e_pos].v2 = v_start + u * v_steps + v + 1;
+                       }
+               }
+       }
+}
+
+static void bridge_loops(Mesh *me, int e_start_a, int e_start_b, int totvert, bool flip, int a_stride, int b_stride, bool n_flip)
+{
+       int e_start = me->totedge;
+       int l_start = me->totloop;
+       int p_start = me->totpoly;
+       ED_mesh_edges_add(me, NULL, totvert);
+       ED_mesh_loops_add(me, NULL, 4 * totvert - 4);
+       ED_mesh_polys_add(me, NULL, totvert - 1);
+
+       for (int i = 0; i < totvert; i++) {
+               me->medge[e_start + i].crease = 0;
+               me->medge[e_start + i].bweight = 0;
+               me->medge[e_start + i].flag = 0;
+
+               if (i < totvert - 1) {
+                       if (flip) {
+                               me->medge[e_start + i].v1 = me->medge[e_start_a + i * a_stride].v1;
+                               me->medge[e_start + i].v2 = me->medge[e_start_b - i * b_stride].v2;
+                       } else {
+                               me->medge[e_start + i].v1 = me->medge[e_start_a + i * a_stride].v1;
+                               me->medge[e_start + i].v2 = me->medge[e_start_b + i * b_stride].v1;
+                       }
+
+                       me->mloop[l_start + i * 4 + 0].v = me->medge[e_start_a + i * a_stride].v1;
+                       me->mloop[l_start + i * 4 + 0].e = e_start_a + i * a_stride;
+
+                       me->mloop[l_start + i * 4 + (n_flip? 3 : 1)].v = me->medge[e_start_a + i * a_stride].v2;
+                       me->mloop[l_start + i * 4 + (n_flip? 3 : 1)].e = e_start + i + 1;
+
+                       if (flip) {
+                               me->mloop[l_start + i * 4 + 2].v = me->medge[e_start_b - i * b_stride].v1;
+                               me->mloop[l_start + i * 4 + 2].e = e_start_b - i * b_stride;
+                               me->mloop[l_start + i * 4 + (n_flip? 1 : 3)].v = me->medge[e_start_b - i * b_stride].v2;
+                       } else {
+                               me->mloop[l_start + i * 4 + 2].v = me->medge[e_start_b + i * b_stride].v2;
+                               me->mloop[l_start + i * 4 + 2].e = e_start_b + i * b_stride;
+                               me->mloop[l_start + i * 4 + (n_flip? 1 : 3)].v = me->medge[e_start_b + i * b_stride].v1;
+                       }
+
+                       me->mloop[l_start + i * 4 + (n_flip? 1 : 3)].e = e_start + i;
+
+                       me->mpoly[p_start + i].loopstart = l_start + i * 4;
+                       me->mpoly[p_start + i].totloop = 4;
+                       me->mpoly[p_start + i].mat_nr = 0;
+                       me->mpoly[p_start + i].flag = 0;
+                       me->mpoly[p_start + i].pad = 0;
+               } else {
+                       if (flip) {
+                               me->medge[e_start + i].v1 = me->medge[e_start_a + (i - 1) * a_stride].v2;
+                               me->medge[e_start + i].v2 = me->medge[e_start_b - (i - 1) * b_stride].v1;
+                       } else {
+                               me->medge[e_start + i].v1 = me->medge[e_start_a + (i - 1) * a_stride].v2;
+                               me->medge[e_start + i].v2 = me->medge[e_start_b + (i - 1) * b_stride].v2;
+                       }
+               }
+       }
+}
+
+/*
+ * Generate a quad from three edges. Returning the newly created edge.
+ *  ___a___
+ *  |      |
+ *  b     new
+ *  |      |
+ *  ___c___
+ */
+static int add_quad(Mesh *me, int edge_b, int edge_a, int edge_c, bool flip)
+{
+       int v_a, v_c;
+       MEdge e_a, e_b, e_c;
+       int e_start = me->totedge;
+       int l_start = me->totloop;
+       int p_start = me->totpoly;
+       bool b_flip = false;
+
+       e_a = me->medge[edge_a];
+       e_b = me->medge[edge_b];
+       e_c = me->medge[edge_c];
+
+       if (e_a.v1 == e_b.v1 || e_a.v1 == e_b.v2){
+               b_flip = (e_a.v1 == e_b.v1);
+               v_a = e_a.v2;
+       } else {
+               b_flip = (e_a.v2 == e_b.v1);
+               v_a = e_a.v1;
+       }
+       if (e_c.v1 == e_b.v1 || e_c.v1 == e_b.v2){
+               v_c = e_c.v2;
+       } else {
+               v_c = e_c.v1;
+       }
+
+       ED_mesh_edges_add(me, NULL, 1);
+
+       me->medge[e_start].v2 = v_a;
+       me->medge[e_start].v1 = v_c;
+       me->medge[e_start].crease = 0;
+       me->medge[e_start].bweight = 0;
+       me->medge[e_start].flag = 0;
+
+       ED_mesh_loops_add(me, NULL, 4);
+       me->mloop[l_start].v = b_flip ? e_b.v1 : e_b.v2;
+       me->mloop[l_start].e = edge_a;
+
+       me->mloop[l_start + (flip ? 1 : 3)].v = v_a;
+       me->mloop[l_start + (flip ? 1 : 3)].e = e_start;
+
+       me->mloop[l_start + 2].v = v_c;
+       me->mloop[l_start + 2].e = edge_c;
+
+       me->mloop[l_start + (flip ? 3 : 1)].v = b_flip ? e_b.v2 : e_b.v1;
+       me->mloop[l_start + (flip ? 3 : 1)].e = edge_b;
+
+       ED_mesh_polys_add(me, NULL, 1);
+       me->mpoly[p_start].loopstart = l_start;
+       me->mpoly[p_start].totloop = 4;
+       me->mpoly[p_start].mat_nr = 0;
+       me->mpoly[p_start].flag = 0;
+       me->mpoly[p_start].pad = 0;
+
+       return e_start;
+}
+
+static void add_face(Mesh *me, int e1, int e2, int e3, int e4, int l_start, int p_start, bool flip){
+       int comp_v1;
+
+       if (me->medge[e1].v1 == me->medge[e2].v1 || me->medge[e1].v1 == me->medge[e2].v2) {
+               comp_v1 = me->medge[e1].v1;
+               me->mloop[l_start].v = me->medge[e1].v2;
+               me->mloop[l_start].e = e1;
+       } else {
+               comp_v1 = me->medge[e1].v2;
+               me->mloop[l_start].v = me->medge[e1].v1;
+               me->mloop[l_start].e = e1;
+       }
+       if (me->medge[e2].v1 == comp_v1) {
+               me->mloop[l_start + (flip ? 3 : 1)].v = me->medge[e2].v1;
+               comp_v1 = me->medge[e2].v2;
+       } else {
+               me->mloop[l_start + (flip ? 3 : 1)].v = me->medge[e2].v2;
+               comp_v1 = me->medge[e2].v1;
+       }
+       me->mloop[l_start + (flip ? 3 : 1)].e = e2;
+
+       if (me->medge[e3].v1 == comp_v1) {
+               me->mloop[l_start + 2].v = me->medge[e3].v1;
+               comp_v1 = me->medge[e3].v2;
+       } else {
+               me->mloop[l_start + 2].v = me->medge[e3].v2;
+               comp_v1 = me->medge[e3].v1;
+       }
+       me->mloop[l_start + 2].e = e3;
+
+       if (me->medge[e4].v1 == comp_v1) {
+               me->mloop[l_start + (flip ? 1 : 3)].v = me->medge[e4].v1;
+               comp_v1 = me->medge[e4].v2;
+       } else {
+               me->mloop[l_start + (flip ? 1 : 3)].v = me->medge[e4].v2;
+               comp_v1 = me->medge[e4].v1;
+       }
+       me->mloop[l_start + (flip ? 1 : 3)].e = e4;
+
+       me->mpoly[p_start].totloop = 4;
+       me->mpoly[p_start].loopstart = l_start;
+       me->mpoly[p_start].mat_nr = 0;
+       me->mpoly[p_start].flag = 0;
+       me->mpoly[p_start].pad = 0;
+}
+
+/* TODO: is there a sort function already?*/
+static int cmpfunc (const void * a, const void * b)
+{
+       return ( *(int*)a - *(int*)b );
+}
+
+/* Generate the Tube shape for branches with two ends. */
+static void fill_tube(Mesh *me, float *left, float *right, int totl, int totr, int u_steps, float z_vec[3], int v_steps, int w_steps, float smoothness, int *r_edge_loop_ends, bool n_g_flip, bool flip_side)
+{
+       float step_l = left[totl * 4 - 1] / (float)u_steps;
+       float step_r = right[totr * 4 - 1] / (float)u_steps;
+       float a, b, f;
+       float v1[3], v2[3], v3[3], v4[3];
+       int l_u_pos_i = 1, r_u_pos_i = totr - 2;
+
+       const int v_start = me->totvert;
+       const int e_start = me->totedge;
+
+       for (int u = 0; u < u_steps; u++) {
+               while (l_u_pos_i < totl - 1 && left[l_u_pos_i * 4 + 3] <= step_l * (float)u) {
+                       l_u_pos_i ++;
+               }
+
+               while (r_u_pos_i > 0 && right[r_u_pos_i * 4 + 3] > step_r * (float)(u_steps - u - 1)) {
+                       r_u_pos_i --;
+               }
+
+               /* Interpolate over the points of each side. Interpolate an even point distribution along the line */
+               if (totl > 1) {
+                       a = left[l_u_pos_i * 4 - 1];
+                       b = left[l_u_pos_i * 4 + 3];
+                       if (a != b) {
+                               f = (step_l * (float)u - a) / (b - a);
+                       } else {
+                               f = 0.0f;
+                       }
+                       interp_v3_v3v3(v1, &left[l_u_pos_i * 4 - 4], &left[l_u_pos_i * 4], f);
+               } else {
+                       copy_v3_v3(v1, &left[0]);
+               }
+
+               if (totr > 1) {
+                       a = right[r_u_pos_i * 4 + 7];
+                       b = right[r_u_pos_i * 4 + 3];
+                       if (a != b) {
+                               f = (step_r * (float)(u_steps - u - 1) - a) / (b - a);
+                       } else {
+                               f = 0.0f;
+                       }
+                       interp_v3_v3v3(v2, &right[r_u_pos_i * 4 + 4], &right[r_u_pos_i * 4], f);
+               } else {
+                       copy_v3_v3(v2, &right[0]);
+               }
+
+               add_v3_v3v3(v3, v1, v2);
+               mul_v3_fl(v3, 0.5f);
+               add_v3_v3v3(v4, v3, z_vec);
+
+               /* v1 left, v2 right, v3 center bottom, v4 center top */
+               calc_vert_half(me, v1, v2, v3, v4, v_steps, w_steps, smoothness, flip_side);
+       }
+       generate_mesh_grid_f_e(me, u_steps, v_steps * 2 + w_steps - (flip_side ? 2 : 0), v_start, n_g_flip);
+       if (!flip_side) {
+               r_edge_loop_ends[0] = e_start;
+               r_edge_loop_ends[1] = 2;
+               r_edge_loop_ends[2] = me->totedge - v_steps * 2 - w_steps + 1;
+               r_edge_loop_ends[3] = 1;
+       } else {
+               r_edge_loop_ends[4] = e_start;
+               r_edge_loop_ends[5] = 2;
+               r_edge_loop_ends[6] = me->totedge - v_steps * 2 - w_steps + 3;
+               r_edge_loop_ends[7] = 1;
+       }
+}
+
+static int get_cyclic_offset(SpineBranch *branch)
+{
+       int cyclic_offset = 0, n_i = 0;
+       if (branch->hull_points[0] == 0) {
+               for (int i = 0; i < branch->tot_hull_points; i++) {
+                       if (n_i > 0) {
+                               cyclic_offset ++;
+                       }
+                       if (branch->hull_points[i] + 1 != branch->hull_points[i + 1] &&
+                               branch->hull_points[i] != branch->hull_points[i + 1] &&
+                               i < branch->tot_hull_points - 1)
+                       {
+                               n_i ++;
+                               cyclic_offset = 0;
+                       }
+               }
+       }
+       return cyclic_offset;
+}
+
+/* Generate the Cap for branches with one ends. */
+static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, Mesh *me, float z_vec[3], float depth, int v_steps, int w_steps, float smoothness, bool n_g_flip, bool flip_side)
+{
+       float *cap_p = NULL;
+       float v1[3], m_center[3], m_center_up[3], left_ref[3], right_ref[3];
+       float totlength = 0.0f;
+       float step_size = depth / (float)w_steps;
+       int cyclic_offset = 0, n_i = 0;
+       int e_cap_start_a, e_cap_start_b, e_cap_start_c;
+       int e_corner_a, e_corner_b;
+       int cap_end_flip_e_start, cap_end_flip_start_a_l, cap_end_flip_start_b_l;
+       int cap_end_flip_start_a_r, cap_end_flip_start_b_r, e_cap_tube_start, e_cap_tube_end;
+       int e_flip_tube_end[2];
+       BLI_array_declare(cap_p);
+
+       if (!flip_side) {
+               branch->fs_bs_offset = me->totedge;
+       }
+
+       /* calc and sort hullpoints for the three sides */
+       qsort (branch->hull_points, branch->tot_hull_points, sizeof(int), cmpfunc);
+
+       cyclic_offset = get_cyclic_offset(branch);
+
+       for (int i = 0; i < branch->tot_hull_points; i++) {
+               n_i = branch->tot_hull_points + i - cyclic_offset;
+               BLI_array_grow_items(cap_p, 4);
+               silhoute_stroke_point_to_3d(sil, branch->hull_points[n_i % branch->tot_hull_points] * 3, &cap_p[i * 4]);
+               if (i > 0) {
+                       cap_p[i * 4 + 3] = len_v3v3(&cap_p[i * 4], v1) + cap_p[i * 4 - 1];
+               } else {
+                       cap_p[i * 4 + 3] = 0.0f;
+               }
+               copy_v3_v3(v1, &cap_p[i * 4]);
+       }
+
+       if (!(branch->flag & BRANCH_EDGE_GEN)) {
+               branch->e_start_arr = MEM_callocN(sizeof(int) * 4, "edge startposition array");
+               branch->e_flip_side_ends = MEM_callocN(sizeof(int) * 2, "fs bs edges");
+               branch->flag |= BRANCH_EDGE_GEN;
+       }
+
+       totlength = cap_p[branch->tot_hull_points * 4 - 1];
+
+       float cap_length = fmin(totlength, (w_steps * step_size));
+       step_size = cap_length / (float)w_steps;
+       float cap_pos = (totlength - cap_length) * 0.5f;
+       float a, b, f;
+       int u_pos_i = 0;
+       int v_start = me->totvert;
+       int u_steps;
+       int e_start_tube[8];
+       int totl = 0, totr = 0;
+       int e_flip_start;
+       int e_flip_offset = 0;  /* carry edgecount difference in both sides to refrence opposing edges by subtracting totedge - flip offset. Only valid if flip_side = true*/
+
+       /* If the cap is big enough a tube is added between the cap and the last branch. */
+       if (totlength > step_size * w_steps) {
+               int side_l = cap_pos;
+               float *left = NULL, *right = NULL;
+               float n_off_right;
+               BLI_array_declare(left);
+               BLI_array_declare(right);
+
+               while (cap_p[u_pos_i * 4 + 3] < side_l) {
+                       BLI_array_grow_items(left, 4);
+                       left[totl * 4    ] = cap_p[u_pos_i * 4    ];
+                       left[totl * 4 + 1] = cap_p[u_pos_i * 4 + 1];
+                       left[totl * 4 + 2] = cap_p[u_pos_i * 4 + 2];
+                       left[totl * 4 + 3] = cap_p[u_pos_i * 4 + 3];
+                       totl ++;
+                       u_pos_i ++;
+               }
+               while (cap_p[u_pos_i * 4 + 3] < totlength - side_l) {
+                       u_pos_i ++;
+               }
+               while (u_pos_i < branch->tot_hull_points) {
+                       if (totr == 0) {
+                               n_off_right = cap_p[u_pos_i * 4 + 3];
+                       }
+                       BLI_array_grow_items(right, 4);
+                       right[totr * 4    ] = cap_p[u_pos_i * 4    ];
+                       right[totr * 4 + 1] = cap_p[u_pos_i * 4 + 1];
+                       right[totr * 4 + 2] = cap_p[u_pos_i * 4 + 2];
+                       right[totr * 4 + 3] = cap_p[u_pos_i * 4 + 3] - n_off_right;
+                       totr ++;
+                       u_pos_i ++;
+               }
+
+               if (totl >= 1 && totr >= 1) {
+                       u_steps = fmax(2.0f, fmax(left[totl * 4 - 1], right[totr * 4 - 1]) / (float)(2 * depth / v_steps));
+                       e_flip_start = me->totedge;
+                       fill_tube(me, left, right, totl, totr, u_steps, z_vec, v_steps, w_steps, smoothness, e_start_tube, n_g_flip, flip_side);
+                       copy_v3_v3(left_ref, &left[totl * 4 - 4]);
+                       copy_v3_v3(right_ref, &right[0]);
+                       cap_length = totlength - left[totl * 4 - 1] - right[totr * 4 - 1];
+
+                       if (flip_side) {
+                               branch->e_flip_side_ends[0] = me->totedge;
+                               bridge_loops(me,
+                                                        e_flip_start + 1,
+                                                        e_flip_start - branch->fs_bs_offset + 1,
+                                                        u_steps,
+                                                        false,
+                                                        ((v_steps - 1) * 2 + w_steps) * 2 - 1,
+                                                        (v_steps * 2 + w_steps) * 2 - 1,
+                                                        !n_g_flip);
+                               branch->e_flip_side_ends[1] = me->totedge;
+                               e_flip_tube_end[0] = me->totedge - 1;
+                               bridge_loops(me,
+                                                        e_flip_start + ((v_steps - 1) * 2 + w_steps) * 2 - 2,
+                                                        e_flip_start - branch->fs_bs_offset + (v_steps * 2 + w_steps) * 2 - 2,
+                                                        u_steps,
+                                                        false,
+                                                        ((v_steps - 1) * 2 + w_steps) * 2 - 1,
+                                                        (v_steps * 2 + w_steps) * 2 - 1,
+                                                        n_g_flip);
+                               e_flip_tube_end[1] = me->totedge - 1;
+                               e_flip_offset += u_steps * 2 - 2;
+                       }
+               }
+               BLI_array_free(left);
+               BLI_array_free(right);
+       }
+
+       if (totlength <= step_size * w_steps || totl == 0 || totr == 0) {
+               copy_v3_v3(left_ref, &cap_p[0]);
+               copy_v3_v3(right_ref, &cap_p[branch->tot_hull_points * 4 - 4]);
+       }
+
+       cap_pos = (totlength - cap_length) / 2.0f;
+       step_size = cap_length / (w_steps + 2);
+
+       interp_v3_v3v3(m_center, left_ref, right_ref, 0.5f);
+       add_v3_v3v3(m_center_up, m_center, z_vec);
+
+       /* Add connecting edge */
+       v_start = me->totvert;
+       e_cap_start_a = me->totedge;
+       calc_vert_half(me,
+                                  left_ref,
+                                  right_ref,
+                                  m_center,
+                                  m_center_up,
+                                  v_steps,
+                                  w_steps,
+                                  smoothness,
+                                  flip_side);
+
+       /*TODO connect to flipside */
+       ED_mesh_edges_add(me, NULL, v_steps * 2 + w_steps - 1 - (flip_side ? 2 : 0));
+       for(int v = 0; v < v_steps * 2 + w_steps - 1 - (flip_side ? 2 : 0); v++){
+               me->medge[e_cap_start_a + v].v1 = v_start + v;
+               me->medge[e_cap_start_a + v].v2 = v_start + v + 1;
+               me->medge[e_cap_start_a + v].crease = 0;
+               me->medge[e_cap_start_a + v].bweight = 0;
+               me->medge[e_cap_start_a + v].flag = 0;
+       }
+
+       e_flip_offset += 2;
+
+       cap_pos += step_size;
+
+       u_pos_i = 1;
+       v_start = me->totvert;
+
+       for (int u = 0; u < w_steps; u++) {
+               while (u_pos_i < branch->tot_hull_points && cap_p[u_pos_i * 4 + 3] <= cap_pos) {
+                       u_pos_i ++;
+               }
+
+               interp_v3_v3v3(m_center, left_ref, right_ref, smoothness * 0.5f + (1.0f - smoothness) * ((float)u / w_steps));
+               add_v3_v3v3(m_center_up, m_center, z_vec);
+
+               a = cap_p[u_pos_i * 4 - 1];
+               b = cap_p[u_pos_i * 4 + 3];
+               if (a != b) {
+                       f = (cap_pos - a) / (b - a);
+               } else {
+                       f = 0.0f;
+               }
+               interp_v3_v3v3(v1, &cap_p[u_pos_i * 4 - 4], &cap_p[u_pos_i * 4], f);
+
+               calc_vert_quarter(me,
+                                                 v1,
+                                                 m_center,
+                                                 m_center_up,
+                                                 v_steps,
+                                                 0,
+                                                 smoothness,
+                                                 false,
+                                                 flip_side);
+
+               cap_pos += step_size;
+       }
+       e_cap_start_b = me->totedge;
+       generate_mesh_grid_f_e(me, w_steps, v_steps - (flip_side ? 1 : 0), v_start, n_g_flip);
+       e_cap_start_c = me->totedge;
+
+       if (flip_side) {
+               cap_end_flip_e_start = me->totedge;
+               bridge_loops(me,
+                                        e_cap_start_b + 1,
+                                        e_cap_start_b - branch->fs_bs_offset + e_flip_offset + 1,
+                                        w_steps,
+                                        false,
+                                        v_steps * 2 - 3,
+                                        v_steps * 2 - 1,
+                                        !n_g_flip);
+               e_flip_offset += w_steps - 1;
+       }
+
+       bridge_loops(me,
+                                e_cap_start_a,
+                                e_cap_start_b,
+                                v_steps - (flip_side ? 1 : 0),
+                                false,
+                                1,
+                                2,
+                                !n_g_flip);
+
+       e_corner_a = me->totedge - 1;
+       cap_end_flip_start_a_l = me->totedge - v_steps + 1;
+       cap_end_flip_start_b_l = me->totedge - v_steps + 1 - branch->fs_bs_offset + e_flip_offset;
+
+       bridge_loops(me,
+                                e_cap_start_a + v_steps - (flip_side ? 1 : 0),
+                                e_cap_start_b + 2 * v_steps - 2 - (flip_side ? 2 : 0),
+                                w_steps,
+                                false,
+                                1,
+                                2 * v_steps - 1 - (flip_side ? 2 : 0),
+                                !n_g_flip);
+
+       e_corner_b = me->totedge - 1;
+
+       bridge_loops(me,
+                                e_cap_start_a + v_steps - (flip_side ? 1 : 0) + w_steps,
+                                e_cap_start_c - 1,
+                                v_steps - (flip_side ? 1 : 0),
+                                true,
+                                1,
+                                1,
+                                !n_g_flip);
+       cap_end_flip_start_a_r = me->totedge - 1;
+       cap_end_flip_start_b_r = me->totedge - branch->fs_bs_offset + e_flip_offset + 1;
+
+       ED_mesh_loops_add(me, NULL, 6);
+       ED_mesh_polys_add(me, NULL, 2);
+       /* corner a */
+       me->mloop[me->totloop - 6].v = me->medge[e_cap_start_a + v_steps -(flip_side ? 1 : 0) - 1].v1;
+       me->mloop[me->totloop - 6].e = e_cap_start_a + v_steps -(flip_side ? 1 : 0) - 1;
+
+       me->mloop[me->totloop - (n_g_flip ? 5 : 4)].v = me->medge[e_corner_a + 1].v1;
+       me->mloop[me->totloop - (n_g_flip ? 5 : 4)].e = e_corner_a + 1;
+
+       me->mloop[me->totloop - (n_g_flip ? 4 : 5)].v = me->medge[e_corner_a].v2;
+       me->mloop[me->totloop - (n_g_flip ? 4 : 5)].e = e_corner_a;
+
+       me->mpoly[me->totpoly - 2].loopstart = me->totloop - 6;
+       me->mpoly[me->totpoly - 2].totloop = 3;
+       me->mpoly[me->totpoly - 2].mat_nr = 0;
+       me->mpoly[me->totpoly - 2].flag = 0;
+       me->mpoly[me->totpoly - 2].pad = 0;
+
+       /* corner b*/
+       me->mloop[me->totloop - 3].v = me->medge[e_cap_start_a + v_steps -(flip_side ? 1 : 0) + w_steps - 1].v1;
+       me->mloop[me->totloop - 3].e = e_cap_start_a + v_steps -(flip_side ? 1 : 0) + w_steps - 1;
+
+       me->mloop[me->totloop - (n_g_flip ? 2 : 1)].v = me->medge[e_corner_b + 1].v1;
+       me->mloop[me->totloop - (n_g_flip ? 2 : 1)].e = e_corner_b + 1;
+
+       me->mloop[me->totloop - (n_g_flip ? 1 : 2)].v = me->medge[e_corner_b].v2;
+       me->mloop[me->totloop - (n_g_flip ? 1 : 2)].e = e_corner_b;
+
+       me->mpoly[me->totpoly - 1].loopstart = me->totloop - 3;
+       me->mpoly[me->totpoly - 1].totloop = 3;
+       me->mpoly[me->totpoly - 1].mat_nr = 0;
+       me->mpoly[me->totpoly - 1].flag = 0;
+       me->mpoly[me->totpoly - 1].pad = 0;
+
+       if (totlength <= step_size * w_steps || totl == 0 || totr == 0) {
+               branch->e_start_arr[flip_side ? 2 : 0] = e_cap_start_a;
+               branch->e_start_arr[flip_side ? 3 : 1] = 1;
+       } else {
+               e_cap_tube_start = me->totedge;
+               bridge_loops(me,
+                                        e_cap_start_a,
+                                        e_start_tube[flip_side ? 6 : 2],
+                                        v_steps * 2 + w_steps - (flip_side ? 2 : 0),
+                                        false,
+                                        1,
+                                        1,
+                                        n_g_flip);
+               e_cap_tube_end = me->totedge - 1;
+               e_flip_offset += 2;
+               branch->e_start_arr[flip_side ? 2 : 0] = e_start_tube[flip_side ? 4 : 0];
+               branch->e_start_arr[flip_side ? 3 : 1] = e_start_tube[flip_side ? 5 : 1];
+       }
+
+       if (!flip_side) {
+               branch->fs_bs_offset = me->totedge - branch->fs_bs_offset;
+       } else {
+               add_quad(me, cap_end_flip_e_start, cap_end_flip_start_a_l, cap_end_flip_start_b_l, !n_g_flip);
+               add_quad(me, cap_end_flip_e_start + w_steps - 1, cap_end_flip_start_a_r, cap_end_flip_start_b_r, n_g_flip);
+               if (totlength <= step_size * w_steps || totl == 0 || totr == 0) {
+                       branch->e_flip_side_ends[0] = me->totedge - 2;
+                       branch->e_flip_side_ends[1] = me->totedge - 1;
+               } else {
+                       ED_mesh_loops_add(me, NULL, 8);
+                       ED_mesh_polys_add(me, NULL, 2);
+                       add_face(me,
+                                        e_flip_tube_end[0],
+                                        e_cap_tube_start - branch->fs_bs_offset + e_flip_offset,
+                                        me->totedge - 2,
+                                        e_cap_tube_start,
+                                        me->totloop - 8,
+                                        me->totpoly - 2,
+                                        n_g_flip);
+
+                       add_face(me,
+                                        e_flip_tube_end[1],
+                                        e_cap_tube_end - branch->fs_bs_offset + e_flip_offset + 2,
+                                        me->totedge - 1,
+                                        e_cap_tube_end,
+                                        me->totloop - 4,
+                                        me->totpoly - 1,
+                                        !n_g_flip);
+               }
+       }
+
+       BLI_array_free(cap_p);
+}
+
+static int calc_branch_orientation(Spine *spine, SpineBranch *branch, float point[3], int f_s)
+{
+       int s = 0;
+       float last_dist = len_v3v3(&spine->branches[branch->terminal_points[1]]->points[spine->branches[branch->terminal_points[1]]->totpoints * 3 - 3], point);
+       if (spine->branches[branch->terminal_points[3]]->totpoints < 2 || spine->branches[branch->terminal_points[5]]->totpoints < 2) {
+               return f_s;
+       }
+       float curr_dist = len_v3v3(&spine->branches[branch->terminal_points[3]]->points[3], point);
+       if(last_dist > curr_dist){
+               last_dist = curr_dist;
+               s = 1;
+       }
+       curr_dist = len_v3v3(&spine->branches[branch->terminal_points[5]]->points[3], point);
+       if(last_dist > curr_dist){
+               last_dist = curr_dist;
+               s = 2;
+       }
+       return s;
+}
+
+/* Generate a T-Intersection for branches with three ends. */
+static void add_ss_tinter(SilhouetteData *sil, Spine *spine, SpineBranch *branch, Mesh *me, float z_vec[3], float depth, int v_steps, int w_steps, float smoothness, bool n_g_flip, bool flip_side)
+{
+       float *sa = NULL;
+       int b_start[3] = {0,0,0}, b_tot[3] = {0,0,0};
+       int filler = 0;
+       int cyclic_offset = 0, n_i = 0;
+       int pos_i_sa, u_steps;
+       float step_length, a, b, f;
+       float v1[3], v2[3], v3[3], v4[3], center[3], center_up[3];
+       float center_s[3 * 3];
+       int w_h_steps = w_steps / 2;
+       int v_start, v_start_center;
+       int e_start[3], e_start_center, e_start_inner[3], e_t_sign[6];
+       int stride_le;
+       int ori[3];
+       int e_flip_offset = 0;  /* carry edgecount difference in both sides to refrence opposing edges by subtracting totedge - flip offset. Only valid if flip_side = true*/
+       int e_flip_start[3];
+       int e_flip_q_l[3], e_flip_q_r[3];
+       BLI_array_declare(sa);
+
+       if (!flip_side) {
+               branch->fs_bs_offset = me->totedge;
+       }
+
+       /* calc and sort hullpoints for the three sides */
+       qsort (branch->hull_points, branch->tot_hull_points, sizeof(int), cmpfunc);
+
+       if (branch->hull_points[0] == 0) {
+               for (int i = 0; i < branch->tot_hull_points; i++) {
+                       if (n_i > 2) {
+                               cyclic_offset ++;
+                       }
+                       if (branch->hull_points[i] + 1 != branch->hull_points[i + 1] &&
+                               branch->hull_points[i] != branch->hull_points[i + 1] &&
+                               i < branch->tot_hull_points - 1)
+                       {
+                               n_i ++;
+                               cyclic_offset = 0;
+                       }
+               }
+       }
+
+       b_start[0] = 0;
+       filler = 0;
+       for (int i = 0; i < branch->tot_hull_points; i++) {
+               n_i = branch->tot_hull_points + i - cyclic_offset;
+               BLI_array_grow_items(sa, 4);
+               silhoute_stroke_point_to_3d(sil, branch->hull_points[n_i % branch->tot_hull_points] * 3, &sa[b_start[filler] + b_tot[filler] * 4]);
+               if(b_tot[filler] == 0){
+                       sa[b_start[filler] + b_tot[filler] * 4 + 3] = 0.0f;
+               }else{
+                       sa[b_start[filler] + b_tot[filler] * 4 + 3] = len_v3v3(&sa[b_start[filler] + b_tot[filler] * 4 - 4], &sa[b_start[filler] + b_tot[filler] * 4]) + sa[b_start[filler] + b_tot[filler] * 4 - 1];
+               }
+               b_tot[filler] ++;
+               if (branch->hull_points[n_i % branch->tot_hull_points] + 1 != branch->hull_points[(n_i + 1) % branch->tot_hull_points] &&
+                  branch->hull_points[n_i % branch->tot_hull_points] != branch->hull_points[(n_i + 1) % branch->tot_hull_points] &&
+                  !(cyclic_offset != 0 && b_tot[0] <= cyclic_offset))
+               {
+                       if ((filler + 1) % 3 > 0){
+                               b_start[(filler + 1) % 3] = b_start[filler] + b_tot[filler] * 4;
+                       }
+                       filler = (filler + 1) % 3;
+               }
+       }
+
+       u_steps = 5;
+       zero_v3(center);
+       for (int s = 0; s < 3; s++) {
+               copy_v3_v3(&center_s[s * 3], &sa[b_start[s]]);
+               add_v3_v3(&center_s[s * 3], &sa[b_start[(s + 2) % 3] + b_tot[(s + 2) % 3] * 4 - 4]);
+               mul_v3_fl(&center_s[s * 3], 0.5f);
+               add_v3_v3(center, &center_s[s * 3]);
+       }
+       mul_v3_fl(center, 1.0f/3.0f);
+       add_v3_v3v3(center_up, center, z_vec);
+
+       for (int s = 0; s < 3; s++) {
+               /*TODO: Better max funktion (int)*/
+               u_steps = fmax(u_steps, len_v3v3(&center_s[(s + 1) % 3 * 3], &center_s[s * 3]) / (float)(2 * depth / v_steps));
+       }
+
+       /*needs to be uneven*/
+       u_steps |= 1;
+
+       if (!(branch->flag & BRANCH_EDGE_GEN)) {
+               branch->e_start_arr = MEM_callocN(sizeof(int) * 12,"edge startposition array");
+               branch->e_flip_side_ends = MEM_callocN(sizeof(int) * 6, "fs bs edges");
+               branch->flag |= BRANCH_EDGE_GEN;
+       }
+
+       v_start_center = me->totvert;
+       ED_mesh_vertices_add(me, NULL, u_steps + u_steps / 2);
+       e_start_center = me->totedge;
+       ED_mesh_edges_add(me, NULL, (u_steps / 2) * 3);
+
+       for (int s = 0; s < 3; s++) {
+               step_length = sa[b_start[s] + b_tot[s] * 4 - 1] / (float)u_steps;
+
+               add_v3_v3v3(v3, &center_s[s * 3], z_vec);
+
+               v_start = me->totvert;
+
+               ori[s] = calc_branch_orientation(spine, branch, &center_s[s * 3], s);
+               branch->e_start_arr[(flip_side ? 6 : 0) + ori[s] * 2] = me->totedge;
+               branch->e_start_arr[(flip_side ? 6 : 0) + ori[s] * 2 + 1] = 1;
+
+               calc_vert_half(me,
+                                          &sa[b_start[s]],
+                                          &sa[b_start[(s + 2) % 3] + b_tot[(s + 2) % 3] * 4 - 4],
+                                          &center_s[s * 3],
+                                          v3,
+                                          v_steps,
+                                          w_steps,
+                                          smoothness,
+                                          flip_side);
+
+               e_start[s] = me->totedge;
+
+               ED_mesh_edges_add(me, NULL, v_steps * 2 + w_steps - 1 - (flip_side ? 2 : 0));
+               for(int v = 0; v < v_steps * 2 + w_steps - 1 - (flip_side ? 2 : 0); v++){
+                       me->medge[e_start[s] + v].v1 = v_start + v;
+                       me->medge[e_start[s] + v].v2 = v_start + v + 1;
+                       me->medge[e_start[s] + v].crease = 0;
+                       me->medge[e_start[s] + v].bweight = 0;
+                       me->medge[e_start[s] + v].flag = 0;
+               }
+
+               e_flip_offset += 2;
+
+               v_start = me->totvert;
+
+               for (int u = 1; u < u_steps - 1; u++) {
+                       pos_i_sa = 1;
+                       while (pos_i_sa < b_tot[s] && sa[b_start[s] + pos_i_sa * 4 + 3] <= step_length * (float)u) {
+                               pos_i_sa ++;
+                       }
+
+                       /* Interpolate over the points of each side. Interpolate an even point distribution along the line */
+                       if (b_tot[s] < 2) {
+                               copy_v3_v3(v1, &sa[b_start[s]]);
+                       } else {
+                               a = sa[b_start[s] + pos_i_sa * 4 - 1];
+                               b = sa[b_start[s] + pos_i_sa * 4 + 3];
+                               if (a != b) {
+                                       f = (step_length * (float)u - a) / (b - a);
+                               } else {
+                                       f = 0.0f;
+                               }
+                               interp_v3_v3v3(v1, &sa[b_start[s] + pos_i_sa * 4 - 4], &sa[b_start[s] + pos_i_sa * 4], f);
+                       }
+
+                       f = fabs((float)(u_steps / 2 - u) / ((float)u_steps / 2.0f));
+                       if (u < u_steps / 2){
+                               interp_v3_v3v3(v4, center, &center_s[s * 3], f);
+                               add_v3_v3v3(v2, v4, z_vec);
+
+                               copy_v3_v3(me->mvert[v_start_center + s * (u_steps / 2) + u].co, v2);
+                               me->mvert[v_start_center + s * (u_steps / 2) + u].flag = 0;
+                               me->mvert[v_start_center + s * (u_steps / 2) + u].bweight = 0;
+                       } else if (u == u_steps / 2) {
+                               copy_v3_v3(v4, center);
+                               add_v3_v3v3(v2, v4, z_vec);
+                               if (s == 0) {
+                                       /* Add center at v2 */
+                                       copy_v3_v3(me->mvert[v_start_center].co, v2);
+                                       me->mvert[v_start_center].flag = 0;
+                                       me->mvert[v_start_center].bweight = 0;
+                               }
+                       } else {
+                               interp_v3_v3v3(v4, center, &center_s[(s + 1) % 3 * 3], f);
+                               add_v3_v3v3(v2, v4, z_vec);
+                       }
+
+                       calc_vert_quarter(me, v1, v4, v2, v_steps, w_h_steps, smoothness, false, flip_side);
+               }
+
+               me->medge[e_start_center + s * (u_steps / 2)].v1 = me->medge[e_start[s] + v_steps - (flip_side ? 1 : 0) + w_steps / 2].v1;
+               me->medge[e_start_center + s * (u_steps / 2)].v2 = v_start_center + s * (u_steps / 2) + 1;
+               me->medge[e_start_center + s * (u_steps / 2)].crease = 0;
+               me->medge[e_start_center + s * (u_steps / 2)].bweight = 0;
+               me->medge[e_start_center + s * (u_steps / 2)].flag = 0;
+
+               for (int u = 1; u < u_steps / 2 - 1; u++) {
+                       me->medge[e_start_center + s * (u_steps / 2) + u].v1 = v_start_center + s * (u_steps / 2) + u;
+                       me->medge[e_start_center + s * (u_steps / 2) + u].v2 = v_start_center + s * (u_steps / 2) + 1 + u;
+                       me->medge[e_start_center + s * (u_steps / 2) + u].crease = 0;
+                       me->medge[e_start_center + s * (u_steps / 2) + u].bweight = 0;
+                       me->medge[e_start_center + s * (u_steps / 2) + u].flag = 0;
+               }
+
+               me->medge[e_start_center + (s + 1) * (u_steps / 2) - 1].v1 = v_start_center + s * (u_steps / 2) + 1 + (u_steps / 2 - 2);
+               me->medge[e_start_center + (s + 1) * (u_steps / 2) - 1].v2 = v_start_center;
+               me->medge[e_start_center + (s + 1) * (u_steps / 2) - 1].crease = 0;
+               me->medge[e_start_center + (s + 1) * (u_steps / 2) - 1].bweight = 0;
+               me->medge[e_start_center + (s + 1) * (u_steps / 2) - 1].flag = 0;
+
+               e_start_inner[s] = me->totedge;
+               generate_mesh_grid_f_e(me, u_steps - 2, v_steps - (flip_side ? 1 : 0) + w_steps / 2, v_start, n_g_flip);
+
+               if (flip_side) {
+                       e_flip_start[s] = me->totedge;
+                       bridge_loops(me,
+                                                e_start_inner[s] + 1,
+                                                e_start_inner[s] - branch->fs_bs_offset + e_flip_offset + 1,
+                                                u_steps - 2,
+                                                false,
+                                                (v_steps - 1 + w_steps / 2) * 2 - 1,
+                                                (v_steps + w_steps / 2) * 2 - 1,
+                                                !n_g_flip);
+
+                       e_flip_offset += u_steps - 3;
+               }
+       }
+
+       for(int s = 0; s < 3; s++){
+               int e_end_a, e_end_b;
+               bridge_loops(me,
+                                        e_start_inner[s],
+                                        e_start[s],
+                                        v_steps -(flip_side ? 1 : 0) + w_steps / 2,
+                                        false,
+                                        2,
+                                        1,
+                                        n_g_flip);
+               e_flip_offset += 1;
+
+               if (flip_side) {
+                       e_flip_q_l[0] = e_flip_start[s];
+                       e_flip_q_l[1] = me->totedge - (v_steps + w_steps / 2) + 1;
+                       e_flip_q_l[2] = me->totedge - (v_steps + w_steps / 2) - branch->fs_bs_offset + e_flip_offset;
+               }
+
+               e_end_a = me->totedge;
+               bridge_loops(me,
+                                        e_start_inner[(s + 2) % 3] + (u_steps - 3) * (2 * (v_steps -(flip_side ? 1 : 0) + w_steps / 2) - 1),
+                                        e_start[s] + v_steps * 2 -(flip_side ? 2 : 0) + w_steps - 2,
+                                        v_steps -(flip_side ? 1 : 0) + w_steps / 2,
+                                        true,
+                                        1,
+                                        1,
+                                        !n_g_flip);
+
+               if (flip_side) {
+                       e_flip_q_r[0] = e_flip_start[(s + 2) % 3] + u_steps - 3;
+                       e_flip_q_r[1] = me->totedge - (v_steps + w_steps / 2) + 1;
+                       e_flip_q_r[2] = me->totedge - (v_steps + w_steps / 2) + 1 - branch->fs_bs_offset + e_flip_offset;
+               }
+
+               e_end_b = me->totedge;
+
+               int e_side_a, e_side_b;
+               e_side_a = add_quad(me, e_start[s] + v_steps -(flip_side ? 1 : 0) + w_steps / 2, e_start_center + s * (u_steps / 2), e_end_b - 1, !n_g_flip);
+               e_side_b = add_quad(me, e_start[s] + v_steps -(flip_side ? 1 : 0) + w_steps / 2 - 1, e_end_a - 1, e_start_center + s * (u_steps / 2), !n_g_flip);
+               stride_le = (2 * (v_steps -(flip_side ? 1 : 0) + w_steps / 2) - 1);
+
+               for (int u = 1; u < u_steps / 2 - 1; u++) {
+                       e_side_b = add_quad(me,
+                                                               e_side_b,
+                                                               e_start_inner[s] - 1 + stride_le * u,
+                                                               e_start_center + s * (u_steps / 2) + u,
+                                                               !n_g_flip);
+                       e_side_a = add_quad(me,
+                                                               e_side_a,
+                                                               e_start_inner[(s + 2) % 3] - 1 + stride_le * (u_steps - 2 - u),
+                                                               e_start_center + s * (u_steps / 2) + u,
+                                                               n_g_flip);
+               }
+
+               e_t_sign[s * 2] = add_quad(me,
+                                e_side_b,
+                                e_start_inner[s] - 1 + stride_le * ((u_steps / 2) - 1),
+                                e_start_center + s * (u_steps / 2) + u_steps / 2 - 1,
+                                !n_g_flip);
+
+               e_t_sign[s * 2 + 1] = e_side_a;
+
+               if (flip_side) {
+                       branch->e_flip_side_ends[ori[s] * 2] = me->totedge;
+                       add_quad(me,
+                                        e_flip_q_l[0],
+                                        e_flip_q_l[1],
+                                        e_flip_q_l[2],
+                                        !n_g_flip);
+                       branch->e_flip_side_ends[ori[s] * 2 + 1] = me->totedge;
+                       add_quad(me,
+                                        e_flip_q_r[0],
+                                        e_flip_q_r[1],
+                                        e_flip_q_r[2],
+                                        n_g_flip);
+                       e_flip_offset -= 1;
+               }
+       }
+
+       for(int s = 0; s < 3; s++){
+               ED_mesh_loops_add(me, NULL, 4);
+               me->mloop[me->totloop - 4].v = me->medge[e_t_sign[(s + 2) % 3 * 2]].v1;
+               me->mloop[me->totloop - 4].e = e_t_sign[(s + 2) % 3 * 2];
+
+               me->mloop[me->totloop - (n_g_flip ? 3 : 1)].v = me->medge[e_start_center + s * (u_steps / 2) + u_steps / 2 - 1].v1;
+               me->mloop[me->totloop - (n_g_flip ? 3 : 1)].e = e_start_center + s * (u_steps / 2) + u_steps / 2 - 1;
+
+               me->mloop[me->totloop - 2].v = me->medge[e_start_inner[(s + 2) % 3] - 1 + stride_le * ((u_steps / 2))].v2;
+               me->mloop[me->totloop - 2].e = e_t_sign[s * 2 + 1];
+
+               me->mloop[me->totloop - (n_g_flip ? 1 : 3)].v = me->medge[e_start_inner[(s + 2) % 3] - 1 + stride_le * ((u_steps / 2))].v1;
+               me->mloop[me->totloop - (n_g_flip ? 1 : 3)].e = e_start_inner[(s + 2) % 3] - 1 + stride_le * ((u_steps / 2));
+
+               ED_mesh_polys_add(me, NULL, 1);
+               me->mpoly[me->totpoly - 1].loopstart = me->totloop - 4;
+               me->mpoly[me->totpoly - 1].totloop = 4;
+               me->mpoly[me->totpoly - 1].mat_nr = 0;
+               me->mpoly[me->totpoly - 1].flag = 0;
+               me->mpoly[me->totpoly - 1].pad = 0;
+
+       }
+
+       if (!flip_side) {
+               branch->fs_bs_offset = me->totedge - branch->fs_bs_offset;
+       }
+
+       BLI_array_free(sa);
+}
+
+static void add_ss_tube(SilhouetteData *sil, SpineBranch *branch, Mesh *me, float z_vec[3], float depth, int v_steps, int w_steps, float w_fact, bool n_g_flip, bool flip_side)
+{
+       /* x y z l (accumulative length)*/
+       float *left = NULL, *right = NULL;
+       int totl = 0, totr = 0;
+       int u_steps = 0;
+       bool f_swap = false;
+       int cyclic_offset = 0, n_i = 0;
+       int e_start = 0;
+       BLI_array_declare(left);
+       BLI_array_declare(right);
+
+       if (!flip_side) {
+               branch->fs_bs_offset = me->totedge;
+       }
+       /* Calc and sort Hullpoints to left and right side */
+       qsort (branch->hull_points, branch->tot_hull_points, sizeof(int), cmpfunc);
+
+       if (branch->hull_points[0] == 0) {
+               for (int i = 0; i < branch->tot_hull_points; i++) {
+                       if (n_i > 1) {
+                               cyclic_offset ++;
+                       }
+                       if (branch->hull_points[i] + 1 != branch->hull_points[i + 1] &&
+                               branch->hull_points[i] != branch->hull_points[i + 1] &&
+                               i < branch->tot_hull_points - 1)
+                       {
+                               n_i ++;
+                               cyclic_offset = 0;
+                       }
+               }
+       }
+
+       /* fill left and right arrays with hull points. */
+       for (int i = 0; i < branch->tot_hull_points; i++){
+               n_i = branch->tot_hull_points + (i - cyclic_offset);
+               if (!f_swap) {
+                       BLI_array_grow_items(left, 4);
+                       silhoute_stroke_point_to_3d(sil, branch->hull_points[n_i % branch->tot_hull_points] * 3, &left[totl * 4]);
+                       if (totl > 0) {
+                               left[totl * 4 + 3] = len_v3v3(&left[totl * 4],&left[totl * 4 - 4]) + left[totl * 4 - 1];
+                       } else {
+                               left[totl * 4 + 3] = 0.0f;
+                       }
+                       totl ++;
+               } else {
+                       BLI_array_grow_items(right, 4);
+                       silhoute_stroke_point_to_3d(sil, branch->hull_points[n_i % branch->tot_hull_points] * 3, &right[totr*4]);
+                       if (totr > 0) {
+                               right[totr * 4 + 3] = len_v3v3(&right[totr * 4],&right[totr * 4 - 4]) + right[totr * 4 - 1];
+                       } else {
+                               right[totr * 4 + 3] = 0.0f;
+                       }
+                       totr ++;
+               }
+               if (branch->hull_points[n_i % branch->tot_hull_points] + 1 != branch->hull_points[(n_i + 1) % branch->tot_hull_points] &&
+                       branch->hull_points[n_i % branch->tot_hull_points] != branch->hull_points[(n_i + 1) % branch->tot_hull_points] &&
+                       !(cyclic_offset != 0 && totl <= cyclic_offset))
+               {
+                       f_swap = !f_swap;
+               }
+       }
+
+       if (totl < 1 && totr < 1) {
+               return;
+       }
+
+       u_steps = fmax(2.0f, fmax(left[totl * 4 - 1], right[totr * 4 - 1]) / (float)(2 * depth / v_steps));
+
+       if (!(branch->flag & BRANCH_EDGE_GEN)) {
+               branch->e_start_arr = MEM_callocN(sizeof(int) * 8, "edge startposition array");
+               branch->e_flip_side_ends = MEM_callocN(sizeof(int) * 4, "fs bs edges");
+               branch->flag |= BRANCH_EDGE_GEN;
+       }
+
+       e_start = me->totedge;
+
+       fill_tube(me, left, right, totl, totr, u_steps, z_vec, v_steps, w_steps, w_fact, branch->e_start_arr, n_g_flip, flip_side);
+
+       if (flip_side) {
+               branch->e_flip_side_ends[0] = me->totedge;
+               bridge_loops(me,
+                                        e_start + 1,
+                                        e_start - branch->fs_bs_offset + 1,
+                                        u_steps,
+                                        false,
+                                        ((v_steps - 1) * 2 + w_steps) * 2 - 1,
+                                        (v_steps * 2 + w_steps) * 2 - 1,
+                                        !n_g_flip);
+               branch->e_flip_side_ends[1] = me->totedge;
+               branch->e_flip_side_ends[2] = me->totedge - 1;
+               bridge_loops(me,
+                                        e_start + ((v_steps - 1) * 2 + w_steps) * 2 - 2,
+                                        e_start - branch->fs_bs_offset + (v_steps * 2 + w_steps) * 2 - 2,
+                                        u_steps,
+                                        false,
+                                        ((v_steps - 1) * 2 + w_steps) * 2 - 1,
+                                        (v_steps * 2 + w_steps) * 2 - 1,
+                                        n_g_flip);
+               branch->e_flip_side_ends[3] = me->totedge - 1;
+       } else {
+               branch->fs_bs_offset = me->totedge - branch->fs_bs_offset;
+       }
+
+       BLI_array_free(left);
+       BLI_array_free(right);
+}
+
+/* TODO: Fix for incomplete Spine Generation. Delete? It calculates the real amount of adjacent branches not counting empty ones.*/
+static int r_branch_count(Spine *spine, SpineBranch *b)
+{
+       int r_forks = 0;
+       for (int i = 0; i < b->totforks; i++) {
+               if (spine->branches[b->terminal_points[i * 2 + 1]]) {
+                       r_forks ++;
+               }
+       }
+       return r_forks;
+}
+
+/* TODO: T-Intersections are sometimes misordered! Connects the different Branches if they have the BRANCH_EDGE_GEN flag set. */
+static void bridge_all_parts_rec(Mesh *me, Spine *spine, SpineBranch *active_branch, SpineBranch *prev_branch, int verts_per_loop, bool n_g_flip)
+{
+       int b_fork_off, a_fork_off, a_fork_off_inv, b_fork_off_inv, comp_e_start, comp_e_start_flip;
+       float dist_a, dist_b;
+       a_fork_off = 0;
+       for (int i = 0; i < active_branch->totforks; i++) {
+               SpineBranch *comp = spine->branches[active_branch->terminal_points[i * 2 + 1]];
+               if (comp && comp != prev_branch) {
+                       if (active_branch->flag & BRANCH_EDGE_GEN && comp->flag & BRANCH_EDGE_GEN) {
+                               b_fork_off = 0;
+                               for (int sb = 0; sb < comp->totforks; sb++) {
+                                       if (spine->branches[comp->terminal_points[sb * 2 + 1]]) {
+                                               if (comp->terminal_points[sb * 2 + 1] == active_branch->idx) {
+                                                       break;
+                                               }
+                                               b_fork_off ++;
+                                       }
+                               }
+
+                               /* TODO: Might fail with thin geometry */
+                               dist_a = len_v3v3(me->mvert[me->medge[active_branch->e_start_arr[a_fork_off * 2]].v1].co, me->mvert[me->medge[comp->e_start_arr[b_fork_off * 2]].v1].co);
+                               dist_b = len_v3v3(me->mvert[me->medge[active_branch->e_start_arr[a_fork_off * 2]].v1].co, me->mvert[me->medge[comp->e_start_arr[b_fork_off * 2] + (verts_per_loop - 2 - 0) * comp->e_start_arr[b_fork_off * 2 + 1]].v2].co);
+
+                               /* TODO: Tube bug inverse?
+                               if (comp->hull_points[0] == 0) {
+                                       bl_debug_color_set(0x0000ff);
+                                       bl_debug_draw_medge_add(me, comp->e_start_arr[b_fork_off * 2] + comp->e_start_arr[b_fork_off * 2 + 1] * (verts_per_loop - 2));
+                                       bl_debug_color_set(0x000000);
+                               }*/
+                               a_fork_off_inv = a_fork_off * 2 + active_branch->totforks * 2;
+                               b_fork_off_inv = b_fork_off * 2 + comp->totforks * 2;
+                               if (dist_a > dist_b) {
+                                       comp_e_start = comp->e_start_arr[b_fork_off * 2] + comp->e_start_arr[b_fork_off * 2 + 1] * (verts_per_loop - 2);
+                                       comp_e_start_flip = comp->e_start_arr[b_fork_off_inv] + comp->e_start_arr[b_fork_off_inv + 1] * (verts_per_loop - 4);
+                               } else {
+                                       comp_e_start = comp->e_start_arr[b_fork_off * 2];
+                                       comp_e_start_flip = comp->e_start_arr[b_fork_off_inv];
+                               }
+                               bridge_loops(me,
+                                                        active_branch->e_start_arr[a_fork_off * 2],
+                                                        comp_e_start,
+                                                        verts_per_loop,
+                                                        (dist_a > dist_b),
+                                                        active_branch->e_start_arr[a_fork_off * 2 + 1],
+                                                        comp->e_start_arr[b_fork_off * 2 + 1],
+                                                        !n_g_flip ^ (dist_a > dist_b));
+
+                               bridge_loops(me,
+                                                        active_branch->e_start_arr[a_fork_off_inv],
+                                                        comp_e_start_flip,
+                                                        verts_per_loop - 2,
+                                                        (dist_a > dist_b),
+                                                        active_branch->e_start_arr[a_fork_off_inv + 1],
+                                                        comp->e_start_arr[b_fork_off_inv + 1],
+                                                        n_g_flip ^ (dist_a > dist_b));
+
+                               ED_mesh_loops_add(me, NULL, 8);
+                               ED_mesh_polys_add(me, NULL, 2);
+                               add_face(me,
+                                                me->totedge - 1,
+                                                active_branch->e_flip_side_ends[a_fork_off * 2 + 1],
+                                                me->totedge - verts_per_loop + 1,
+                                                comp->e_flip_side_ends[b_fork_off * 2 + ((dist_a > dist_b) ? 0 : 1)],
+                                                me->totloop - 8,
+                                                me->totpoly - 2,
+                                                n_g_flip ^ (dist_a > dist_b));
+
+                               add_face(me,
+                                                me->totedge - verts_per_loop + 2,
+                                                active_branch->e_flip_side_ends[a_fork_off * 2],
+                                                me->totedge - verts_per_loop * 2 + 2,
+                                                comp->e_flip_side_ends[b_fork_off * 2 + ((dist_a > dist_b) ? 1 : 0)],
+                                                me->totloop - 4,
+                                                me->totpoly - 1,
+                                                !n_g_flip ^ (dist_a > dist_b));
+                       }
+                       bridge_all_parts_rec(me, spine, comp, active_branch, verts_per_loop, n_g_flip);
+               }
+               if (comp) {
+                       a_fork_off ++;
+               }
+       }
+}
+
+static void bridge_all_parts(Mesh *me, Spine *spine, int verts_per_loop, bool n_g_flip)
+{
+       SpineBranch *active_branch = NULL;
+       for (int i = 0; i < spine->totbranches; i++) {
+               if (!active_branch) {
+                       active_branch = spine->branches[i];
+               } else {
+                       break;
+               }
+       }
+       if (!active_branch) {
+               /* No Branches in the spine. Should not happen. */
+               return;
+       }
+       bridge_all_parts_rec(me, spine, active_branch, NULL, verts_per_loop, n_g_flip);
+}
+
+static bool calc_stroke_normal_ori(SilhouetteStroke *stroke, float z_vec[3]) {
+       float n[3];
+       /*TODO: stroke points to multidimensional array*/
+       cross_poly_v3(n, (float (*)[3])stroke->points, stroke->totvert);
+       return dot_v3v3(n, z_vec) <= 0;
+}
+
+/* Generates a 3D shape from a stroke. */
+static void silhouette_create_shape_mesh(bContext *C, Mesh *me, SilhouetteData *sil, SilhouetteStroke *stroke)
+{
+       float z_vec[3] = {0.0f,0.0f,1.0f};
+       float inv_z_vec[3];
+       float depth = sil->depth;
+       int ss_level = sil->resolution;
+       int v_steps = (1 << ss_level) + 2;
+       bool n_ori = false;
+       /* TODO: RNA Init*/
+
+       copy_v3_v3(z_vec, sil->z_vec);
+       normalize_v3(z_vec);
+
+       n_ori = calc_stroke_normal_ori(stroke, z_vec);
+
+       /* Generate spine */
+       Spine *spine = silhouette_generate_spine(sil, stroke);
+       SpineBranch *a_branch;
+
+       mul_v3_fl(z_vec, depth);
+       mul_v3_v3fl(inv_z_vec, z_vec, -1.0f);
+
+       int w_steps = v_steps / 2 + 2;
+       float smoothness = sil->smoothness;
+
+       for (int i = 0; i < spine->totbranches; i++) {
+               a_branch = spine->branches[i];
+               if (a_branch) {
+                       int r_forks = r_branch_count(spine, a_branch);
+                       switch (r_forks) {
+                               case 1:
+                                       add_ss_cap(sil, a_branch, me, z_vec, depth, v_steps, w_steps, smoothness, n_ori, false);
+                                       add_ss_cap(sil, a_branch, me, inv_z_vec, depth, v_steps, w_steps, smoothness, !n_ori, true);
+#ifdef DEBUG_DRAW
+                                       /*debug_branch(a_branch, 0x00ff00);*/
+#endif
+                                       break;
+                               case 2:
+                                       add_ss_tube(sil, a_branch, me, z_vec, depth, v_steps, w_steps, smoothness, n_ori, false);
+                                       add_ss_tube(sil, a_branch, me, inv_z_vec, depth, v_steps, w_steps, smoothness, !n_ori, true);
+#ifdef DEBUG_DRAW
+                                       /*debug_branch(a_branch, 0xff0000);*/
+#endif
+                                       break;
+                               case 3:
+                                       add_ss_tinter(sil, spine, a_branch, me, z_vec, depth, v_steps, w_steps, smoothness, n_ori, false);
+                                       add_ss_tinter(sil, spine, a_branch, me, inv_z_vec, depth, v_steps, w_steps, smoothness, !n_ori, true);
+#ifdef DEBUG_DRAW
+                                       /*debug_branch(a_branch, 0x0000ff);*/
+#endif
+                                       break;
+                       }
+               }
+       }
+
+       bridge_all_parts(me, spine, v_steps * 2 + w_steps, n_ori);
+
+       free_spine(spine);
+
+       ED_mesh_update(me, C, 1, 1);
+}
+
+/* Adds additional points to the stroke if start and end are far apart. */
+static void stroke_smooth_cap(SilhouetteData *sil, SilhouetteStroke *stroke, float max_dist)
+{
+       float v1[2], v2[2], dv[2];
+       float length = 0;
+       copy_v2_v2(v1, &stroke->points_v2[0]);
+       copy_v2_v2(v2, &stroke->points_v2[stroke->totvert * 2 - 2]);
+
+       sub_v2_v2v2(dv,v1,v2);
+       length = len_v2(dv);
+
+       if (length > max_dist) {
+               int steps = floorf(length / (float)max_dist);
+               mul_v2_fl(dv, 1.0f / (float)steps);
+               for (int i = 1; i < steps; i++) {
+                       mul_v2_v2fl(v1, dv, i);
+                       add_v2_v2(v1, v2);
+                       silhouette_stroke_add_point(sil, stroke, v1);
+               }
+       }
+}
+
+static void remove_connected_from_edgehash(MeshElemMap *emap, GHash *edge_hash, int v) {
+       for (int e = 0; e < emap[v].count; e++) {
+               BLI_ghash_remove(edge_hash, emap[v].indices[e], NULL, NULL);
+       }
+}
+
+static bool has_cross_border_neighbour(Mesh *me, GHash *vert_hash, GHash *edge_hash, MeshElemMap *emap, int edge, int l_v_edge, int depth) {
+       int v_edge;
+
+       v_edge = me->medge[edge].v1 == l_v_edge ? me->medge[edge].v2 : me->medge[edge].v1;
+
+       if (depth == 0) {
+               return BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(v_edge));
+       } else {
+               if(!BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(v_edge))){
+                       for (int e = 0; e < emap[v_edge].count; e++) {
+                               if(emap[v_edge].indices[e] != edge) {
+                                       if(has_cross_border_neighbour(me, vert_hash, edge_hash, emap, emap[v_edge].indices[e], v_edge, depth - 1)){
+                                               return true;
+                                       }
+                               }
+                       }
+               } else {
+                       BLI_ghash_remove(edge_hash, edge, NULL, NULL);
+               }
+       }
+       return false;
+}
+
+/* Get the adjacent edge which connects the edges within the edge_hash. Used to create multiple ordered loops
+ * v_edge is the endpoint off curr_edge from which to branch off
+ * TODO: One wide strips might get cutoff */
+static int get_adjacent_edge(Mesh *me, MeshElemMap *emap, int curr_edge, int v_edge, GHash *edge_hash, GHash *vert_hash)
+{
+       for (int e = 0; e < emap[v_edge].count; e++) {
+               if(emap[v_edge].indices[e] != curr_edge && has_cross_border_neighbour(me, vert_hash, edge_hash, emap, emap[v_edge].indices[e], v_edge, 1)) {
+                       return emap[v_edge].indices[e];
+               }
+       }
+       for (int e = 0; e < emap[v_edge].count; e++) {
+               if(emap[v_edge].indices[e] != curr_edge && has_cross_border_neighbour(me, vert_hash, edge_hash, emap, emap[v_edge].indices[e], v_edge, 2)) {
+                       return emap[v_edge].indices[e];
+               }
+       }
+       /*End Of Loop. Shouldn't happen with two manifold meshes*/
+       return -1;
+}
+
+/*TODO:Remove Temp debug function*/
+static char *cd_type_name(int t)
+{
+       char *name = "";
+       switch(t) {
+               case -1 : name = "CD_AUTO_FROM_NAME      ";break;
+               case  0 : name = "CD_MVERT               ";break;
+               case  1 : name = "CD_MSTICKY             ";break;  /* DEPRECATED */
+               case  2 : name = "CD_MDEFORMVERT         ";break;
+               case  3 : name = "CD_MEDGE               ";break;
+               case  4 : name = "CD_MFACE               ";break;
+               case  5 : name = "CD_MTFACE              ";break;
+               case  6 : name = "CD_MCOL                ";break;
+               case  7 : name = "CD_ORIGINDEX           ";break;
+               case  8 : name = "CD_NORMAL          ";break;
+                       /*      CD_POLYINDEX        = 9, */
+               case 10 : name = "CD_PROP_FLT        ";break;
+               case 11 : name = "CD_PROP_INT        ";break;
+               case 12 : name = "CD_PROP_STR        ";break;
+               case 13 : name = "CD_ORIGSPACE       ";break;  /* for modifier stack face location mapping */
+               case 14 : name = "CD_ORCO            ";break;
+               case 15 : name = "CD_MTEXPOLY        ";break;
+               case 16 : name = "CD_MLOOPUV         ";break;
+               case 17 : name = "CD_MLOOPCOL        ";break;
+               case 18 : name = "CD_TANGENT         ";break;
+               case 19 : name = "CD_MDISPS          ";break;
+               case 20 : name = "CD_PREVIEW_MCOL    ";break;  /* for displaying weightpaint colors */
+                       /*      CD_ID_MCOL          = 21, */
+               case 22 : name = "CD_TEXTURE_MLOOPCOL";break;
+               case 23 : name = "CD_CLOTH_ORCO      ";break;
+               case 24 : name = "CD_RECAST          ";break;
+
+                       /* BMESH ONLY START */
+               case 25 : name = "CD_MPOLY           ";break;
+               case 26 : name = "CD_MLOOP           ";break;
+               case 27 : name = "CD_SHAPE_KEYINDEX  ";break;
+               case 28 : name = "CD_SHAPEKEY        ";break;
+               case 29 : name = "CD_BWEIGHT         ";break;
+               case 30 : name = "CD_CREASE          ";break;
+               case 31 : name = "CD_ORIGSPACE_MLOOP ";break;
+               case 32 : name = "CD_PREVIEW_MLOOPCOL";break;
+               case 33 : name = "CD_BM_ELEM_PYPTR   ";break;
+                       /* BMESH ONLY END */
+               case 34 : name = "CD_PAINT_MASK      ";break;
+               case 35 : name = "CD_GRID_PAINT_MASK ";break;
+               case 36 : name = "CD_MVERT_SKIN      ";break;
+               case 37 : name = "CD_FREESTYLE_EDGE  ";break;
+               case 38 : name = "CD_FREESTYLE_FACE  ";break;
+               case 39 : name = "CD_MLOOPTANGENT    ";break;
+               case 40 : name = "CD_TESSLOOPNORMAL  ";break;
+               case 41 : name = "CD_CUSTOMLOOPNORMAL";break;
+
+               case 42 : name = "CD_NUMTYPES        ";break;
+               default: name = "No Name";
+       }
+       return name;
+}
+
+#if 0
+static void debug_cd(Mesh *me)
+{
+       char *name = "";
+       printf("Debugging Custom Data:\n\n");
+       printf("%i Customdata Layers in vdata\n", CustomData_number_of_layers_typemask(&me->vdata, CD_MASK_EVERYTHING));
+       printf("%i Customdata Layers in edata\n", CustomData_number_of_layers_typemask(&me->edata, CD_MASK_EVERYTHING));
+       printf("%i Customdata Layers in ldata\n", CustomData_number_of_layers_typemask(&me->ldata, CD_MASK_EVERYTHING));
+       printf("%i Customdata Layers in pdata\n", CustomData_number_of_layers_typemask(&me->pdata, CD_MASK_EVERYTHING));
+       for (int i = 0; i < CustomData_number_of_layers_typemask(&me->vdata, CD_MASK_EVERYTHING); i++) {
+               name = cd_type_name(me->vdata.layers[i].type);
+               printf("Layer found with name %s\n", name);
+       }
+       for (int i = 0; i < CustomData_number_of_layers_typemask(&me->edata, CD_MASK_EVERYTHING); i++) {
+               name = cd_type_name(me->edata.layers[i].type);
+               printf("Layer found with name %s\n", name);
+       }
+       for (int i = 0; i < CustomData_number_of_layers_typemask(&me->ldata, CD_MASK_EVERYTHING); i++) {
+               name = cd_type_name(me->ldata.layers[i].type);
+               printf("Layer found with name %s\n", name);
+       }
+       for (int i = 0; i < CustomData_number_of_layers_typemask(&me->pdata, CD_MASK_EVERYTHING); i++) {
+               name = cd_type_name(me->pdata.layers[i].type);
+               printf("Layer found with name %s\n", name);
+       }
+}
+#endif
+
+/* Copy customdata from source to destination. Layers need to be matching!
+ * redirect_map needs to be the new positions ofelements in the dest array. -1 marks elements which do not get copied over. */
+static void CustomData_copy_partial(const struct CustomData *source, struct CustomData *dest, CustomDataMask mask, int *redirect_map, int tot_elem)
+{
+       int num_layers_to_copy;
+       int d_size;
+       CustomDataLayer *curr_l, *dest_l;
+       char *name;
+
+       num_layers_to_copy = CustomData_number_of_layers_typemask(source, mask);
+
+       for (int l = 0; l < num_layers_to_copy; l++) {
+               curr_l = &source->layers[l];
+               dest_l = &dest->layers[l];
+               if (CD_TYPE_AS_MASK(curr_l->type) & mask) {
+                       for (int i = 0; i < tot_elem; i++) {
+                               if(redirect_map[i] != -1) {
+                                       d_size = CustomData_sizeof(curr_l->type);
+                                       memcpy(dest_l->data + redirect_map[i] * d_size, curr_l->data + i * d_size, d_size);
+                               }
+                       }
+                       name = cd_type_name(curr_l->type);
+                       printf("Layer %s copied %i values of size %i.\n", name, tot_elem, CustomData_sizeof(curr_l->type));
+               }
+       }
+}
+
+/* Function used to remove vertices from a basic mesh
+ * Can be optimised easily be introucing multithreading
+ * Parts and ideas how multithreading can be introduced are marked with "MT:"
+ * TODO: Add Multithreading */
+static void remove_verts_from_mesh(Mesh *me, int *v_to_rm, int num_v_to_rm){
+       int *v_rd_table, *e_rd_table, *l_rd_table, *p_rd_table;
+       int next_del = 0, next_del_pos = 0;
+       MEdge e;
+       MLoop l;
+       MPoly p;
+
+       int n_totvert, n_totedge, n_totloop, n_totpoly;
+
+       CustomData vdata, edata, ldata, pdata;
+
+       /* MT: Prefix Sum / Scan to calculate new positions for vertices.
+        * Calculating the new positions with the vertices removed
+        */
+       v_rd_table = MEM_callocN(sizeof(int) * me->totvert, "Vertex redirect table");
+
+       qsort(v_to_rm, num_v_to_rm, sizeof(int), cmpfunc);
+
+       n_totvert = 0;
+       next_del = v_to_rm[0];
+       for (int i = 0; i < me->totvert; i++) {
+               if (i == next_del) {
+                       if (next_del_pos + 1 < num_v_to_rm) {
+                               next_del_pos ++;
+                               next_del = v_to_rm[next_del_pos];
+                       }
+                       v_rd_table[i] = -1;
+               } else {
+                       v_rd_table[i] = n_totvert;
+                       n_totvert ++;
+               }
+       }
+
+       /* MT: parallelize loop, !shared access to sum value
+        * Calculate new edge positions + remap vertice pointer
+        */
+       e_rd_table = MEM_callocN(sizeof(int) * me->totedge, "Edge redirect table");
+       n_totedge = 0;
+       for (int i = 0; i < me->totedge; i++) {
+               e = me->medge[i];
+               if(v_rd_table[e.v1] == -1 || v_rd_table[e.v2] == -1) {
+                       e_rd_table[i] = -1;
+               } else {
+                       e.v1 = v_rd_table[e.v1];
+                       e.v2 = v_rd_table[e.v2];
+                       me->medge[i] = e;
+                       e_rd_table[i] = n_totedge;
+                       n_totedge ++;
+               }
+       }
+
+       /* MT: same as above
+        * Calculate new loop positions + remap edge pointers */
+       l_rd_table = MEM_callocN(sizeof(int) * me->totloop, "Loop redirect table");
+       n_totloop = 0;
+       for (int i = 0; i < me->totloop; i++) {
+               l = me->mloop[i];
+               if(v_rd_table[l.v] == -1 || e_rd_table[l.e] == -1) {
+                       l_rd_table[i] = -1;
+               } else {
+                       l.v = v_rd_table[l.v];
+                       l.e = e_rd_table[l.e];
+                       me->mloop[i] = l;
+                       l_rd_table[i] = n_totloop;
+                       n_totloop ++;
+               }
+       }
+
+       /* MT: same as above
+        * Calculate new poly positions + remap pointers */
+       p_rd_table = MEM_callocN(sizeof(int) * me->totpoly, "Poly redirect table");
+       n_totpoly = 0;
+       for (int i = 0; i < me->totpoly; i++) {
+               p = me->mpoly[i];
+               for(int l = p.loopstart; l < p.loopstart + p.totloop; l++){
+                       if(l_rd_table[l] == -1) {
+                               p_rd_table[i] = -1;
+                               /* TODO: Bad practise? easily solved other way*/
+                               goto skip_poly;
+                       }
+               }
+               me->mpoly[i].loopstart = l_rd_table[me->mpoly[i].loopstart];
+               p_rd_table[i] = n_totpoly;
+               n_totpoly ++;
+               skip_poly:;
+       }
+
+       /*Redirection tables are done. Continue to copy and allocate new Customdata blocks*/
+       CustomData_copy(&me->vdata, &vdata, CD_MASK_EVERYTHING, CD_CALLOC, n_totvert);
+       CustomData_copy(&me->edata, &edata, CD_MASK_EVERYTHING, CD_CALLOC, n_totedge);
+       CustomData_copy(&me->ldata, &ldata, CD_MASK_EVERYTHING, CD_CALLOC, n_totloop);
+       CustomData_copy(&me->pdata, &pdata, CD_MASK_EVERYTHING, CD_CALLOC, n_totpoly);
+
+       CustomData_copy_partial(&me->vdata, &vdata, CD_MASK_EVERYTHING, v_rd_table, me->totvert);
+       CustomData_copy_partial(&me->edata, &edata, CD_MASK_EVERYTHING, e_rd_table, me->totedge);
+       CustomData_copy_partial(&me->ldata, &ldata, CD_MASK_EVERYTHING, l_rd_table, me->totloop);
+       CustomData_copy_partial(&me->pdata, &pdata, CD_MASK_EVERYTHING, p_rd_table, me->totpoly);
+
+       CustomData_free(&me->vdata, me->totvert);
+       CustomData_free(&me->edata, me->totedge);
+       CustomData_free(&me->ldata, me->totloop);
+       CustomData_free(&me->pdata, me->totpoly);
+
+       me->vdata = vdata;
+       me->edata = edata;
+       me->ldata = ldata;
+       me->pdata = pdata;
+
+       me->totvert = n_totvert;
+       me->totedge = n_totedge;
+       me->totloop = n_totloop;
+       me->totpoly = n_totpoly;
+
+
+       BKE_mesh_update_customdata_pointers(me, true);
+}
+
+static void calc_ring_bbs(SilhouetteData *sil, Mesh *me)
+{
+       int len, edge;
+       BB *curr;
+       sil->fillet_ring_bbs = MEM_callocN(sizeof(BB) * sil->num_rings, "ring bb mem");
+
+       for (int r = 0; r < sil->num_rings; r++) {
+               len = r + 1 < sil->num_rings ? sil->fillet_ring_new_start[r + 1] - sil->fillet_ring_new_start[r] : sil->fillet_ring_tot - sil->fillet_ring_new_start[r];
+               curr = &sil->fillet_ring_bbs[r];
+               BB_reset(curr);
+               for (int e = 0; e < len; e++) {
+                       edge = sil->fillet_ring_new[sil->fillet_ring_new_start[r] + e];
+                       BB_expand(curr, me->mvert[me->medge[edge].v1].co);
+                       BB_expand(curr, me->mvert[me->medge[edge].v2].co);
+               }
+       }
+}
+
+/* We calculate a start and an endpoint of the two node ends intersecting. Now we need to determine the right side on which the intersection happens 
+ * Returning the points in correct order (positive looping results in inner side);
+ */
+static void order_positive_is_inside(Mesh *me, SilhouetteData *sil, MeshElemMap *emap, int *r11, int *r12, int r1_start, int r1_tot, int r2_start, int r2_tot)
+{
+       int dist_1;
+       int center;
+       int comp_v1, comp_v2;
+       int comp_e;
+       int tmp_swap;
+       if (*r11 > *r12) {
+               dist_1 = r1_tot - *r11 + *r12;
+       } else {
+               dist_1 = *r12 - *r11;
+       }
+
+       center = r1_start + (*r11 + (dist_1 / 2)) % r1_tot;
+       comp_v1 = me->medge[sil->fillet_ring_new[center]].v2;
+
+       for (int e = 0; e < emap[comp_v1].count; e++) {
+               comp_e = emap[comp_v1].indices[e];
+               if (comp_e != sil->fillet_ring_new[center] && comp_e != sil->fillet_ring_new[center + 1] && comp_e != sil->fillet_ring_new[center - 1]) {
+                       comp_v2 = me->medge[comp_e].v1 == comp_v1 ? me->medge[comp_e].v2 : me->medge[comp_e].v1;
+                       for (int e_l = 0; e_l < r2_tot; e_l ++) {
+                               if (me->medge[sil->fillet_ring_new[r2_start + e_l]].v2 == comp_v2) {
+                                       return;
+                               }
+                       }
+               }
+       }
+       tmp_swap = *r11;
+       *r11 = *r12;
+       *r12 = tmp_swap;
+       return;
+}
+
+typedef enum MergeRingFlag {
+       ADDED_TO_MERGE = 1,
+} MergeRingFlag;
+
+
+typedef struct MergeRingInfo {
+       int r1;
+       int r2;
+       int r1_start;
+       int r2_start;
+       int r1_tot;
+       int r2_tot;
+       int r1_e_a;
+       int r1_e_b;
+       int r2_e_a;
+       int r2_e_b;
+       int flag;
+} MergeRingInfo;
+
+static void join_node_separated_rings(SilhouetteData *sil, Mesh *me, MeshElemMap *emap)
+{
+       MergeRingInfo *merge_info = NULL;
+       int *merged_to_one = NULL;
+       MEdge e1_c, e2_c;
+       MergeRingInfo t_m_info;
+
+       BLI_array_declare(merge_info);
+
+       printf("Joining rings. In total %i rings to check.\n", sil->num_rings);
+       for (int r1 = 0; r1 < sil->num_rings - 1; r1++) {
+               for (int r2 = r1 + 1; r2 < sil->num_rings; r2++) {
+                       if (bb_intersect(&sil->fillet_ring_bbs[r1], &sil->fillet_ring_bbs[r2])) {
+                               t_m_info.r1_start = sil->fillet_ring_new_start[r1];
+                               t_m_info.r2_start = sil->fillet_ring_new_start[r2];
+                               t_m_info.r1_tot = sil->fillet_ring_new_start[r1 + 1] - t_m_info.r1_start;
+                               t_m_info.r2_tot = r2 + 1 < sil->num_rings ? sil->fillet_ring_new_start[r2 + 1] - t_m_info.r2_start : sil->fillet_ring_tot - t_m_info.r2_start;
+                               t_m_info.r1_e_a = -1, t_m_info.r1_e_b = -1, t_m_info.r2_e_a = -1, t_m_info.r2_e_b = -1;
+                               for (int e1 = 0; e1 < t_m_info.r1_tot; e1++) {
+                                       e1_c = me->medge[sil->fillet_ring_new[t_m_info.r1_start + e1]];
+                                       for (int e2 = 0; e2 < t_m_info.r2_tot; e2++) {
+                                               e2_c = me->medge[sil->fillet_ring_new[t_m_info.r2_start + e2]];
+                                               if (e1_c.v1 == e2_c.v1 || e1_c.v1 == e2_c.v2 || e1_c.v2 == e2_c.v1 || e1_c.v2 == e2_c.v2) {
+                                                       if (t_m_info.r1_e_a == -1) {
+                                                               t_m_info.r1_e_a = e1;
+                                                               t_m_info.r2_e_a = e2;
+                                                       } else {
+                                                               if (abs(t_m_info.r1_e_a - e1) > 3) {
+                                                                       t_m_info.r1_e_b = e1;
+                                                                       t_m_info.r2_e_b = e2;
+                                                                       /* Found start and endpoint of the two ring intersections */
+                                                                       order_positive_is_inside(me, sil, emap, &t_m_info.r1_e_a, &t_m_info.r1_e_b, t_m_info.r1_start, t_m_info.r1_tot, t_m_info.r2_start, t_m_info.r2_tot);
+                                                                       order_positive_is_inside(me, sil, emap, &t_m_info.r2_e_a, &t_m_info.r2_e_b, t_m_info.r2_start, t_m_info.r2_tot, t_m_info.r1_start, t_m_info.r1_tot);
+
+                                                                       BLI_array_append(merge_info, t_m_info);
+#ifdef DEBUG_DRAW
+                                                                       bl_debug_color_set(0xffffff);
+                                                                       bl_debug_draw_point(me->mvert[me->medge[sil->fillet_ring_new[t_m_info.r1_start + t_m_info.r1_e_a]].v1].co, 0.2f);
+                                                                       bl_debug_color_set(0x000000);
+                                                                       bl_debug_draw_point(me->mvert[me->medge[sil->fillet_ring_new[t_m_info.r1_start + t_m_info.r1_e_b]].v1].co, 0.3f);
+                                                                       bl_debug_color_set(0x000000);
+
+                                                                       bl_debug_color_set(0x00ff00);
+                                                                       for (int e_ins = 0; e_ins < t_m_info.r1_tot; e_ins ++) {
+                                                                               if((t_m_info.r1_e_a + e_ins) % t_m_info.r1_tot == t_m_info.r1_e_b) {
+                                                                                       bl_debug_color_set(0x000000);
+                                                                                       break;
+                                                                               }
+                                                                               bl_debug_draw_medge_add(me, sil->fillet_ring_new[t_m_info.r1_start + (t_m_info.r1_e_a + e_ins) % t_m_info.r1_tot]);
+                                                                       }
+                                                                       bl_debug_color_set(0x000000);
+#endif 
+                                                                       /* TODO: Is this a bad coding practise?
+                                                                        * Maybe:
+                                                                        * e1 = r1_tot;
+                                                                        * e2 = r2_tot;
+                                                                        * r2++;
+                                                                        */
+                                                                       goto next_ring;
+                                                               }
+                                                       }
+                                               }
+                                       }
+                               }
+                       }
+                       /* Continue with the next ring */
+                       next_ring:;
+               }
+       }
+
+       BLI_array_declare(merged_to_one);
+       for (int i = 0; i < BLI_array_count(merge_info); i++) {
+               if (!merge_info[i].flag & ADDED_TO_MERGE) {
+                       BLI_array_append(merged_to_one, merge_info[0].r1);
+                       merge_info[0].flag |= ADDED_TO_MERGE;
+                       for (int j = i; j < BLI_array_count(merge_info); j++) {
+                               for (int k = 0; k < BLI_array_count(merged_to_one); k++) {
+
+                               }
+                       }
+               }
+       }
+       BLI_array_free(merged_to_one);
+       BLI_array_free(merge_info);
+}
+
+static void prep_int_shared_mem(int **mem, int *r_num, int *r_start, int len, const char *str)
+{
+       if (!mem) {
+               *r_start = 0;
+               *mem = MEM_callocN(sizeof(int) * len, str);
+               *r_num = len;
+       } else {
+               *r_start = *r_num;
+               *r_num = *r_num + len;
+               *mem = MEM_reallocN(*mem, sizeof(int) * (*r_num));
+       }
+}
+
+/* Calculate the intersections with existing geometry. Used to connect the new Silhouette to the old mesh
+ * Fillets are the smooth transition from one part to another.
+ */
+static void do_calc_fillet_line_task_cb_ex(void *userdata, void *UNUSED(userdata_chunk), const int n, const int UNUSED(thread_id))
+{
+       SculptThreadedTaskData *data = userdata;
+       SculptSession *ss = data->ob->sculpt;
+       SilhouetteData *sil = data->sil;
+       Mesh *me = data->ob->data;
+       PBVHNode *curr_node = data->nodes[n];
+
+       PBVHVertexIter vd;
+       float point[2];
+       float sil_plane[4];
+       float fuzz;
+       MEdge e_comp;
+       int v_i, vd_i;
+       GHash *vert_hash = BLI_ghash_int_new("vertices within intersection");
+       GHash *edge_hash = BLI_ghash_int_new("edges on the intersection");
+       int *edge_ring_fillet = NULL;
+       int *ring_start = NULL;
+       int v_rm_start_in_shared_arr = 0, int_e_start_in_shared_arr = 0, fillet_edge_ring_start_shared_arr = 0, fillet_ring_start_start_shared_arr = 0;
+       int r_size;
+       BLI_array_declare(edge_ring_fillet);
+       BLI_array_declare(ring_start);
+       int comp_v, idx;
+
+       /*GHashIterState state;*/
+       GHashIterator gh_iter;
+       int curr_edge = -1, last_edge = -1, start_edge = -1, tmp_curr_edge = -1;
+
+       plane_from_point_normal_v3(sil_plane, sil->anchor, sil->z_vec);
+
+       /* TODO: Is this a unique vert iteration? If not edges might be added twice by different threads */
+       BKE_pbvh_vertex_iter_begin(ss->pbvh, curr_node, vd, PBVH_ITER_UNIQUE)
+       {
+               /* get the interior vertices of the 2d drawn silhouette and all relevant vertices 
+                * Ignores smoothness, assuming the smoothness blures the fillets anyways it should be ok. */
+               fuzz = SIL_FILLET_BLUR_MIN + sil->smoothness * 0.01f * SIL_FILLET_BLUR_MAX;
+               if (dist_squared_to_plane_v3(vd.co, sil_plane) <= sil->depth + fuzz) {
+                       if (!BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(vd.vert_indices[vd.i]))) {
+                               ED_view3d_project_float_v2_m4(sil->ar, vd.co, point, data->mat);
+                               if (isect_point_poly_v2(point, (float(*)[2])sil->current_stroke->points_v2, sil->current_stroke->totvert, false)){
+                                       if (!BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(vd.vert_indices[vd.i]))) {
+                                               BLI_ghash_insert(vert_hash, SET_INT_IN_POINTER(vd.vert_indices[vd.i]), SET_INT_IN_POINTER(vd.vert_indices[vd.i]));
+                                               vd_i = vd.vert_indices[vd.i];
+                                               for (int e = 0; e < sil->emap[vd_i].count; e++){
+                                                       e_comp = me->medge[sil->emap[vd_i].indices[e]];
+                                                       v_i = e_comp.v1 == vd_i ? e_comp.v2 : e_comp.v1;
+                                                       if (!BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(v_i)))
+                                                       {
+                                                               ED_view3d_project_float_v2_m4(sil->ar, me->mvert[v_i].co, point, data->mat);
+                                                               if (!isect_point_poly_v2(point, (float(*)[2])sil->current_stroke->points_v2, sil->current_stroke->totvert, false) || dist_squared_to_plane_v3(me->mvert[v_i].co, sil_plane) > sil->depth + fuzz) {
+                                                                       BLI_ghash_insert(edge_hash, SET_INT_IN_POINTER(sil->emap[vd_i].indices[e]), SET_INT_IN_POINTER(sil->emap[vd_i].indices[e]));
+                                                               }
+                                                       }
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+       BKE_pbvh_vertex_iter_end;
+
+       /* Finished writing all vertices which are within the intersection and need to be removed.
+        * write them to the shared array. Lock the mutex to avoid collisions */
+       BLI_mutex_lock(&data->mutex);
+       prep_int_shared_mem(&data->v_to_rm, &data->num_v_to_rm, &v_rm_start_in_shared_arr, BLI_ghash_size(vert_hash), "verts to remove");
+       prep_int_shared_mem(&sil->inter_edges, &sil->num_inter_edges, &int_e_start_in_shared_arr, BLI_ghash_size(edge_hash), "edges on transition");
+
+       /* Copy vertice data over.*/
+       GHASH_ITER_INDEX (gh_iter, vert_hash, idx) {
+               data->v_to_rm[v_rm_start_in_shared_arr + idx] = BLI_ghashIterator_getKey(&gh_iter);
+       }
+
+       /* Copy edge data over. */
+       GHASH_ITER_INDEX (gh_iter, edge_hash, idx) {
+               sil->inter_edges[int_e_start_in_shared_arr + idx] = BLI_ghashIterator_getKey(&gh_iter);
+       }
+       BLI_mutex_unlock(&data->mutex);
+
+       /* TODO: Workaround, do a BLI_ghash_pop style while loop
+        * TODO: A adjacency search might fail if there is not a single path to be searched, shouldn't be a problem on first thought though.
+        * Breaker is a anti crash method in case the algorithm gets caught in an endless loop. Shouldn't happen!*/
+       int breaker;
+       int debug_test_v;
+       do {
+               breaker = me->totedge;
+               BLI_ghashIterator_init(&gh_iter, edge_hash);
+               if(!BLI_ghashIterator_done(&gh_iter)){
+                       start_edge = BLI_ghashIterator_getValue(&gh_iter);
+                       BLI_ghash_remove(edge_hash, start_edge, NULL, NULL);
+                       comp_v = BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(me->medge[start_edge].v1)) ? me->medge[start_edge].v2 : me->medge[start_edge].v1;
+                       debug_test_v = comp_v;
+                       BLI_assert(!BLI_ghash_haskey(vert_hash, SET_INT_IN_POINTER(comp_v)));
+                       start_edge = get_adjacent_edge(me, sil->emap, start_edge, comp_v, edge_hash, vert_hash);
+                       if(start_edge >= 0) {
+                               BLI_array_append(ring_start, BLI_array_count(edge_ring_fillet));
+                               curr_edge = start_edge;
+                               last_edge = -1;
+
+                               while(!(curr_edge == start_edge && last_edge != -1) && curr_edge != -1 && breaker > 0) {
+                                       BLI_array_append(edge_ring_fillet, curr_edge);
+                                       if(last_edge == -1) {
+                                               comp_v = me->medge[start_edge].v1;
+                                       } else {
+                                               if (me->medge[curr_edge].v1 == me->medge[last_edge].v1 || me->medge[curr_edge].v1 == me->medge[last_edge].v2) {
+                                                       comp_v = me->medge[curr_edge].v2;
+                                               } else {
+                                                       comp_v = me->medge[curr_edge].v1;
+                                               }
+                                       }
+                                       remove_connected_from_edgehash(sil->emap, edge_hash, comp_v);
+                                       tmp_curr_edge = get_adjacent_edge(me, sil->emap, curr_edge, comp_v, edge_hash, vert_hash);
+                                       last_edge = curr_edge;
+                                       curr_edge = tmp_curr_edge;
+                                       breaker --;
+                               }
+                               printf("Found a cut loop!\n");
+                               BLI_assert(breaker > 0);
+                               /* TODO: Bug shouldn't reach but does on some occasion.*/
+                               if (breaker == 0) {
+                                       BLI_array_empty(edge_ring_fillet);
+#ifdef DEBUG_DRAW
+
+                                       bl_debug_color_set(0x00ff00);
+                                       bl_debug_draw_point(me->mvert[debug_test_v].co, 0.2f);
+                                       bl_debug_draw_medge_add(me, start_edge);
+                                       bl_debug_color_set(0x000000);
+#endif
+                               }
+                       }
+               }
+       } while (!BLI_ghashIterator_done(&gh_iter));
+
+       /* Prep ring memory*/
+       BLI_mutex_lock(&data->mutex);
+       prep_int_shared_mem(&sil->fillet_ring_new, &sil->fillet_ring_tot, &fillet_edge_ring_start_shared_arr, BLI_array_count(edge_ring_fillet), "edges on transition");
+       prep_int_shared_mem(&sil->fillet_ring_new_start, &sil->num_rings, &fillet_ring_start_start_shared_arr, BLI_array_count(ring_start), "start of individual rings");
+
+       /* Copy ring memory */
+       for (int i = 0; i < BLI_array_count(edge_ring_fillet); i++) {
+               sil->fillet_ring_new[fillet_edge_ring_start_shared_arr + i] = edge_ring_fillet[i];
+       }
+       /* Offset start pointers to account chunks beforehand */
+       for (int i = 0; i < BLI_array_count(ring_start); i++) {
+               sil->fillet_ring_new_start[fillet_ring_start_start_shared_arr + i] = ring_start[i] + fillet_edge_ring_start_shared_arr;
+       }
+       BLI_mutex_unlock(&data->mutex);
+
+       /*TODO: merge rings from multiple threads / nodes*/
+#ifdef DEBUG_DRAW
+       for(int r = 0; r < BLI_array_count(ring_start); r++) {
+               r_size = r < BLI_array_count(ring_start) - 1 ? ring_start[r + 1] - ring_start[r] : BLI_array_count(edge_ring_fillet) - ring_start[r];
+               for(int i = 0; i < r_size; i++) {
+                       if(i == 0){
+                               bl_debug_color_set(0x00ffff);
+                       } else {
+                               bl_debug_color_set(0xff00ff);
+                       }
+                       bl_debug_draw_medge_add(me, edge_ring_fillet[ring_start[r] + i]);
+                       bl_debug_color_set(0x000000);
+               }
+       }
+#endif
+       BLI_array_free(ring_start);
+       BLI_array_free(edge_ring_fillet);
+       BLI_ghash_free(vert_hash, NULL, NULL);
+       BLI_ghash_free(edge_hash, NULL, NULL);
+}
+
+static void do_calc_fillet_line(Object *ob, SilhouetteData *silhouette, PBVHNode **nodes, int totnode)
+{
+       Mesh *me = ob->data;
+       float projmat[4][4];
+       MeshElemMap *emap;
+       int *emap_mem;
+       int *v_remove = NULL;
+       RegionView3D *rv3d;
+       View3D *v3d;
+       /*rctf viewplane;
+       float clipend;
+       float clipnear;*/
+
+       rv3d = silhouette->ar->regiondata;
+       v3d = silhouette->vc.v3d;
+
+       BKE_mesh_vert_edge_map_create(&emap, &emap_mem, me->medge, me->totvert, me->totedge);
+       silhouette->emap = emap;
+       /* calc the projection matrix used to convert 3d vertice in 2d space */
+       mul_m4_m4m4(projmat, (float (*)[4])rv3d->persmat, ob->obmat);
+       /* TODO: Orthographic projection:
+       ED_view3d_ob_project_mat_get(silhouette->ar->regiondata, ob, projmat);
+
+       float vmat[4][4], winmat[4][4];;
+
+
+       mul_m4_m4m4(vmat, (float (*)[4])rv3d->viewmat, ob->obmat);
+       mul_m4_m4m4(projmat, (float (*)[4])rv3d->persmat, ob->obmat);
+
+       ED_view3d_viewplane_get(v3d, rv3d, silhouette->ar->winx, silhouette->ar->winy, &viewplane, &clipnear, &clipend, NULL);
+       orthographic_m4(winmat, viewplane.xmin, viewplane.xmax, viewplane.ymin, viewplane.ymax, -clipend, clipend);
+
+       mul_m4_m4m4(vmat, (float (*)[4])rv3d->viewmat, ob->obmat);
+       mul_m4_m4m4(projmat, winmat, vmat);*/
+
+       /* threaded loop over nodes */
+       SculptThreadedTaskData data = {
+               .ob = ob, .nodes = nodes,
+               .sil = silhouette, .mat = projmat, .v_to_rm = v_remove
+       };
+
+       BLI_task_parallel_range_ex(
+                                                          0, totnode, &data, NULL, 0, do_calc_fillet_line_task_cb_ex,
+                                                          (totnode > SCULPT_THREADED_LIMIT), false);
+
+       v_remove = data.v_to_rm;
+
+       calc_ring_bbs(silhouette, me);
+
+#ifdef DEBUG_DRAW
+       for (int r = 0; r < silhouette->num_rings; r ++) {
+               bl_debug_draw_BB_add(&silhouette->fillet_ring_bbs[r], 0x0000ff);
+       }
+#endif
+
+       /*TODO: Join multiple parts together when totnode > 1.*/
+       join_node_separated_rings(silhouette, me, emap);
+
+       if (v_remove) {
+               printf("Removing vertices/edges/loops/polys from mesh.\n");
+               remove_verts_from_mesh(me, v_remove, data.num_v_to_rm);
+               MEM_freeN(v_remove);
+       }
+
+       MEM_freeN(emap);
+       MEM_freeN(emap_mem);
+}
+
+static void sculpt_silhouette_calc_mesh(bContext *C, wmOperator *op)
+{
+       /*finalize stroke*/
+       Object *ob = CTX_data_active_object(C);
+       SilhouetteData *sil = op->customdata;
+       Mesh *me = ob->data;
+       /*Sculpt *sd = CTX_data_tool_settings(C)->sculpt;*/
+       SculptSession *ss = ob->sculpt;
+       SculptSearchBBData data;
+       PBVHNode **nodes = NULL;
+
+       SilhouetteStroke *stroke = sil->current_stroke;
+
+       stroke_smooth_cap(sil, stroke, 10.0f);
+
+       int totnode;
+
+       data.ss = ss;
+       data.bb_target = &stroke->bb;
+       /* get the pbvh nodes intersecting the silhouette BB */
+       BKE_pbvh_search_gather(ss->pbvh, sculpt_search_BB_cb, &data, &nodes, &totnode);
+
+       /* Only act if some verts are inside the silhouette drawn */
+       if (totnode) {
+               printf("Connect to geometry\n");
+               do_calc_fillet_line(ob, sil, nodes, totnode);
+       }
+
+       silhouette_create_shape_mesh(C, me, sil, stroke);
+
+       /* Rebuild mesh caches
+        * TODO: Proper PBVH etc. */
+
+       BKE_mesh_calc_edges(me, false, true);
+       BKE_mesh_tessface_clear(me);
+       BKE_mesh_calc_normals(me);
+       DAG_id_tag_update(&me->id, 0);
+       WM_event_add_notifier(C, NC_GEOM | ND_DATA, me);
+
+       BKE_object_free_derived_caches(ob);
+}
+
+static void sculpt_silhouette_stroke_done(bContext *UNUSED(C), wmOperator *op)
+{
+#ifdef DEBUG_DRAW
+       SilhouetteData *sil = op->customdata;
+
+       for (int i = 1; i < sil->current_stroke->totvert; i ++) {
+               float v1[3], v2[3];
+               copy_v3_v3(v1, &sil->current_stroke->points[i * 3 - 3]);
+               copy_v3_v3(v2, &sil->current_stroke->points[i * 3]);
+               /*bl_debug_draw_edge_add(v1,v2);*/
+       }
+#endif
+       
+       /*cleanup*/
+       silhouette_data_free(op);
+       op->customdata = NULL;
+}
+
+/* Stop drawing the UI*/
+static void sculpt_silhouette_clean_draw(bContext *C, wmOperator *op)
+{
+       SilhouetteData *sil = op->customdata;
+       WM_cursor_modal_restore(CTX_wm_window(C));
+       ED_region_draw_cb_exit(sil->ar->type,sil->draw_handle);
+       ED_region_tag_redraw(sil->ar);
+}
+
+static int sculpt_silhouette_exec(bContext *C, wmOperator *op)
+{
+       int v_start, e_start, l_start, p_start;
+       SilhouetteData *sil = op->customdata;
+       Object *ob = CTX_data_active_object(C);
+       Mesh *me = ob->data;
+
+       if (!sil) {
+               sil = silhouette_data_new(C);
+               silhouette_set_ref_plane(sil);
+               op->customdata = sil;
+       }
+
+       /*TODO: Add undo for fillets etc*/
+       if (sil->current_stroke->totvert > 3) {
+               /*sculpt_undo_push_begin("draw Silhouette");*/
+               v_start = me->totvert;
+               e_start = me->totedge;
+               l_start = me->totloop;
+               p_start = me->totpoly;
+               sculpt_silhouette_calc_mesh(C, op);
+               sculpt_silhouette_stroke_done(C, op);
+
+               /*sculpt_undo_silhouette_push(ob, v_start, e_start, l_start, p_start);
+               sculpt_undo_push_end(C);*/
+       }
+
+       return OPERATOR_FINISHED;
+}
+
+static int sculpt_silhouette_modal(bContext *C, wmOperator *op, const wmEvent *event)
+{
+       float mouse[2];
+       SilhouetteData *sil = op->customdata;
+       copy_v2_fl2(mouse, event->mval[0], event->mval[1]);
+       printf(".");
+       if (event->val == KM_RELEASE) {
+               sculpt_silhouette_clean_draw(C, op);
+               if (sil->state == SIL_DRAWING) {
+                       silhouette_set_ref_plane(sil);
+                       return sculpt_silhouette_exec(C, op);
+               }
+               return OPERATOR_FINISHED;
+       } else {
+               if (sil->state == SIL_DRAWING) {
+                       sculpt_silhouette_stroke_update(mouse, op->customdata);
+               }
+               return OPERATOR_RUNNING_MODAL;
+       }
+}
+
+/* UI Indicator for the tool. TODO: beautify */
+static void sculpt_silhouette_draw(const bContext *UNUSED(C),struct ARegion *UNUSED(ar), void *arg)
+{
+       SilhouetteData *sil = arg;
+
+       if (!sil) {
+               return;
+       }
+       const SilhouetteStroke *stroke = sil->current_stroke;
+       if (!stroke) {
+               return;
+       }
+
+       glLineWidth(1.0f);
+       glEnable(GL_BLEND);
+       glEnable(GL_LINE_SMOOTH);
+
+       /* set brush color */
+       glColor4f(sil->add_col[0],sil->add_col[1],sil->add_col[2], 0.2f);
+
+       if (stroke->points) {
+               glBegin(GL_POLYGON);
+               for (int i = 0; i < stroke->totvert; i++) {
+                       glVertex3f(stroke->points[3 * i], stroke->points[3 * i + 1], stroke->points[3 * i + 2]);
+               }
+               glEnd();
+       }
+
+       /* restore GL state */
+       glDisable(GL_BLEND);
+       glDisable(GL_LINE_SMOOTH);
+}
+
+static int sculpt_silhouette_invoke(bContext *C, wmOperator *op, const wmEvent *UNUSED(event))
+{
+       printf("Drawing Silhouette\n");
+
+       SilhouetteData *sil_data;
+
+       sil_data = silhouette_data_new(C);
+
+       op->customdata = sil_data;
+
+       /* Tag for UI to be drawn */
+       ED_region_tag_redraw(sil_data->ar);
+
+       sil_data->draw_handle = ED_region_draw_cb_activate(sil_data->ar->type, sculpt_silhouette_draw, op->customdata, REGION_DRAW_PRE_VIEW);
+       sil_data->state = SIL_DRAWING;
+
+       /* add modal handler */
+       WM_event_add_modal_handler(C, op);
+
+       return OPERATOR_RUNNING_MODAL;
+}
+
+static int sculpt_silhouette_poll(bContext *UNUSED(C))
+{
+       /* TODO; */
+       return 1;
+}
+
+static void SCULPT_OT_silhouette_draw(wmOperatorType *ot)
+{
+       /* identifiers */
+       ot->name = "Draw Silhouette";
+       ot->idname = "SCULPT_OT_silhouette_draw";
+       ot->description = "Draw a new silhoutte for the sculpt";
+
+       /* api callbacks */
+       ot->invoke = sculpt_silhouette_invoke;
+       ot->modal = sculpt_silhouette_modal;
+       ot->exec = sculpt_silhouette_exec;
+       ot->poll = sculpt_silhouette_poll;
+       ot->cancel = sculpt_silhouette_stroke_done;
+
+       /* flags */
+       ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
+}
+/* end Silhouette */
+
+/************************** Dynamic Topology **************************/
+
+static void sculpt_dynamic_topology_triangulate(BMesh *bm)
+{
+       if (bm->totloop != bm->totface * 3) {
+               BM_mesh_triangulate(bm, MOD_TRIANGULATE_QUAD_BEAUTY, MOD_TRIANGULATE_NGON_EARCLIP, false, NULL, NULL, NULL);
+       }
+}
+
+void sculpt_pbvh_clear(Object *ob)
+{
+       SculptSession *ss = ob->sculpt;
+       DerivedMesh *dm = ob->derivedFinal;
+
+       /* Clear out any existing DM and PBVH */
+       if (ss->pbvh)
+               BKE_pbvh_free(ss->pbvh);
+       ss->pbvh = NULL;
+       if (dm)
+               dm->getPBVH(NULL, dm);
+       BKE_object_free_derived_caches(ob);
+}
+
+void sculpt_dyntopo_node_layers_add(SculptSession *ss)
+{
+       int cd_node_layer_index;
+
+       char layer_id[] = "_dyntopo_node_id";
+
+       cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->vdata, CD_PROP_INT, layer_id);
+       if (cd_node_layer_index == -1) {
+               BM_data_layer_add_named(ss->bm, &ss->bm->vdata, CD_PROP_INT, layer_id);
+               cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->vdata, CD_PROP_INT, layer_id);
+       }
+
+       ss->cd_vert_node_offset = CustomData_get_n_offset(
+               &ss->bm->vdata, CD_PROP_INT,
+               cd_node_layer_index - CustomData_get_layer_index(&ss->bm->vdata, CD_PROP_INT));
+
+       ss->bm->vdata.layers[cd_node_layer_index].flag |= CD_FLAG_TEMPORARY;
+
+       cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->pdata, CD_PROP_INT, layer_id);
+       if (cd_node_layer_index == -1) {
+               BM_data_layer_add_named(ss->bm, &ss->bm->pdata, CD_PROP_INT, layer_id);
+               cd_node_layer_index = CustomData_get_named_layer_index(&ss->bm->pdata, CD_PROP_INT, layer_id);
+       }
+
+       ss->cd_face_node_offset = CustomData_get_n_offset(
+               &ss->bm->pdata, CD_PROP_INT,
+               cd_node_layer_index - CustomData_get_layer_index(&ss->bm->pdata, CD_PROP_INT));
+
+       ss->bm->pdata.layers[cd_node_layer_index].flag |= CD_FLAG_TEMPORARY;
+}
+
+
+void sculpt_update_after_dynamic_topology_toggle(bContext *C)
+{
+       Scene *scene = CTX_data_scene(C);
+       Object *ob = CTX_data_active_object(C);
+       Sculpt *sd = scene->toolsettings->sculpt;
+
+       /* Create the PBVH */
+       BKE_sculpt_update_mesh_elements(scene, sd, ob, false, false);
+       WM_event_add_notifier(C, NC_OBJECT | ND_DRAW, ob);
+}
+
+void sculpt_dynamic_topology_enable(bContext *C)
+{
+       Scene *scene = CTX_data_scene(C);
+       Object *ob = CTX_data_active_object(C);
+       SculptSession *ss = ob->sculpt;
+       Mesh *me = ob->data;
+       const BMAllocTemplate allocsize = BMALLOC_TEMPLATE_FROM_ME(me);
+
+       sculpt_pbvh_clear(ob);
+
+       ss->bm_smooth_shading = (scene->toolsettings->sculpt->flags & SCULPT_DYNTOPO_SMOOTH_SHADING) != 0;
+
+       /* Dynamic topology doesn't ensure selection state is valid, so remove [#36280] */
+       BKE_mesh_mselect_clear(me);
+
+       /* Create triangles-only BMesh */
+       ss->bm = BM_mesh_create(
+               &allocsize,
+               &((struct BMeshCreateParams){.use_toolflags = false,}));
+
+       BM_mesh_bm_from_me(
+               ss->bm, me, (&(struct BMeshFromMeshParams){
+                   .calc_face_normal = true, .use_shapekey = true, .active_shapekey = ob->shapenr,
+               }));
+       sculpt_dynamic_topology_triangulate(ss->bm);
+       BM_data_layer_add(ss->bm, &ss->bm->vdata, CD_PAINT_MASK);
+       sculpt_dyntopo_node_layers_add(ss);
+       /* make sure the data for existing faces are initialized */
+       if (me->totpoly != ss->bm->totface) {
+               BM_mesh_normals_update(ss->bm);
+       }
+
+       /* Enable dynamic topology */
+       me->flag |= ME_SCULPT_DYNAMIC_TOPOLOGY;
+       
+       /* Enable logging for undo/redo */
+       ss->bm_log = BM_log_create(ss->bm);
+
+       /* Refresh */
+       sculpt_update_after_dynamic_topology_toggle(C);
+}
+
+/* Free the sculpt BMesh and BMLog
+ *
+ * If 'unode' is given, the BMesh's data is copied out to the unode
+ * before the BMesh is deleted so that it can be restored from */
+void sculpt_dynamic_topology_disable(bContext *C,
+                                     SculptUndoNode *unode)
+{
+       Object *ob = CTX_data_active_object(C);
+       SculptSession *ss = ob->sculpt;
+       Mesh *me = ob->data;
+
+       sculpt_pbvh_clear(ob);
 
        if (unode) {
                /* Free all existing custom data */
@@ -4894,6 +8220,8 @@ static int sculpt_dynamic_topology_toggle_exec(bContext *C, wmOperator *UNUSED(o
        Object *ob = CTX_data_active_object(C);
        SculptSession *ss = ob->sculpt;
 
+       WM_cursor_wait(1);
+
        if (ss->bm) {
                sculpt_undo_push_begin("Dynamic topology disable");
                sculpt_undo_push_node(ob, NULL, SCULPT_UNDO_DYNTOPO_END);
@@ -4906,16 +8234,24 @@ static int sculpt_dynamic_topology_toggle_exec(bContext *C, wmOperator *UNUSED(o
        }
        sculpt_undo_push_end(C);
 
+       WM_cursor_wait(0);
+
        return OPERATOR_FINISHED;
 }
 
+enum eDynTopoWarnFlag {
+       DYNTOPO_WARN_VDATA = (1 << 0),
+       DYNTOPO_WARN_EDATA = (1 << 1),
+       DYNTOPO_WARN_LDATA = (1 << 2),
+       DYNTOPO_WARN_MODIFIER = (1 << 3),
+};
 
-static int dyntopo_warning_popup(bContext *C, wmOperatorType *ot, bool vdata, bool modifiers)
+static int dyntopo_warning_popup(bContext *C, wmOperatorType *ot, enum eDynTopoWarnFlag flag)
 {
        uiPopupMenu *pup = UI_popup_menu_begin(C, IFACE_("Warning!"), ICON_ERROR);
        uiLayout *layout = UI_popup_menu_layout(pup);
 
-       if (vdata) {
+       if (flag & (DYNTOPO_WARN_VDATA | DYNTOPO_WARN_EDATA | DYNTOPO_WARN_LDATA)) {
                const char *msg_error = TIP_("Vertex Data Detected!");
                const char *msg = TIP_("Dyntopo will not preserve vertex colors, UVs, or other customdata");
                uiItemL(layout, msg_error, ICON_INFO);
@@ -4923,7 +8259,7 @@ static int dyntopo_warning_popup(bContext *C, wmOperatorType *ot, bool vdata, bo
                uiItemS(layout);
        }
 
-       if (modifiers) {
+       if (flag & DYNTOPO_WARN_MODIFIER) {
                const char *msg_error = TIP_("Generative Modifiers Detected!");
                const char *msg = TIP_("Keeping the modifiers will increase polycount when returning to object mode");
 
@@ -4939,33 +8275,35 @@ static int dyntopo_warning_popup(bContext *C, wmOperatorType *ot, bool vdata, bo
        return OPERATOR_INTERFACE;
 }
 
-
-static int sculpt_dynamic_topology_toggle_invoke(bContext *C, wmOperator *op, const wmEvent *UNUSED(event))
+static enum eDynTopoWarnFlag sculpt_dynamic_topology_check(bContext *C)
 {
        Object *ob = CTX_data_active_object(C);
        Mesh *me = ob->data;
        SculptSession *ss = ob->sculpt;
 
-       if (!ss->bm) {
-               Scene *scene = CTX_data_scene(C);
-               ModifierData *md;
-               VirtualModifierData virtualModifierData;
-               int i;
-               bool vdata = false;
-               bool modifiers = false;
-
-               for (i = 0; i < CD_NUMTYPES; i++) {
-                       if (!ELEM(i, CD_MVERT, CD_MEDGE, CD_MFACE, CD_MLOOP, CD_MPOLY, CD_PAINT_MASK, CD_ORIGINDEX) &&
-                           (CustomData_has_layer(&me->vdata, i) ||
-                            CustomData_has_layer(&me->edata, i) ||
-                            CustomData_has_layer(&me->ldata, i)))
-                       {
-                               vdata = true;
-                               break;
+       Scene *scene = CTX_data_scene(C);
+       enum eDynTopoWarnFlag flag = 0;
+
+       BLI_assert(ss->bm == NULL);
+       UNUSED_VARS_NDEBUG(ss);
+
+       for (int i = 0; i < CD_NUMTYPES; i++) {
+               if (!ELEM(i, CD_MVERT, CD_MEDGE, CD_MFACE, CD_MLOOP, CD_MPOLY, CD_PAINT_MASK, CD_ORIGINDEX)) {
+                       if (CustomData_has_layer(&me->vdata, i)) {
+                               flag |= DYNTOPO_WARN_VDATA;
+                       }
+                       if (CustomData_has_layer(&me->edata, i)) {
+                               flag |= DYNTOPO_WARN_EDATA;
+                       }
+                       if (CustomData_has_layer(&me->ldata, i)) {
+                               flag |= DYNTOPO_WARN_LDATA;
                        }
                }
+       }
 
-               md = modifiers_getVirtualModifierList(ob, &virtualModifierData);
+       {
+               VirtualModifierData virtualModifierData;
+               ModifierData *md = modifiers_getVirtualModifierList(ob, &virtualModifierData);
 
                /* exception for shape keys because we can edit those */
                for (; md; md = md->next) {
@@ -4973,14 +8311,26 @@ static int sculpt_dynamic_topology_toggle_invoke(bContext *C, wmOperator *op, co
                        if (!modifier_isEnabled(scene, md, eModifierMode_Realtime)) continue;
 
                        if (mti->type == eModifierTypeType_Constructive) {
-                               modifiers = true;
+                               flag |= DYNTOPO_WARN_MODIFIER;
                                break;
                        }
                }
+       }
 
-               if (vdata || modifiers) {
+       return flag;
+}
+
+static int sculpt_dynamic_topology_toggle_invoke(bContext *C, wmOperator *op, const wmEvent *UNUSED(event))
+{
+       Object *ob = CTX_data_active_object(C);
+       SculptSession *ss = ob->sculpt;
+
+       if (!ss->bm) {
+               enum eDynTopoWarnFlag flag = sculpt_dynamic_topology_check(C);
+
+               if (flag) {
                        /* The mesh has customdata that will be lost, let the user confirm this is OK */
-                       return dyntopo_warning_popup(C, op->type, vdata, modifiers);
+                       return dyntopo_warning_popup(C, op->type, flag);
                }
        }
 
@@ -5055,6 +8405,8 @@ static int sculpt_symmetrize_exec(bContext *C, wmOperator *UNUSED(op))
        sculpt_undo_push_node(ob, NULL, SCULPT_UNDO_DYNTOPO_SYMMETRIZE);
        BM_log_before_all_removed(ss->bm, ss->bm_log);
 
+       BM_mesh_toolflags_set(ss->bm, true);
+
        /* Symmetrize and re-triangulate */
        BMO_op_callf(ss->bm, BMO_FLAG_DEFAULTS,
                     "symmetrize input=%avef direction=%i  dist=%f",
@@ -5064,6 +8416,8 @@ static int sculpt_symmetrize_exec(bContext *C, wmOperator *UNUSED(op))
        /* bisect operator flags edges (keep tags clean for edge queue) */
        BM_mesh_elem_hflag_disable_all(ss->bm, BM_EDGE, BM_ELEM_TAG, false);
 
+       BM_mesh_toolflags_set(ss->bm, false);
+
        /* Finish undo */
        BM_log_all_added(ss->bm, ss->bm_log);
        sculpt_undo_push_end(C);
@@ -5125,14 +8479,21 @@ static int sculpt_mode_toggle_exec(bContext *C, wmOperator *op)
                if (mmd)
                        multires_force_update(ob);
 
-               if (flush_recalc || (ob->sculpt && ob->sculpt->bm))
+               /* Always for now, so leaving sculpt mode always ensures scene is in
+                * a consistent state.
+                */
+               if (true || flush_recalc || (ob->sculpt && ob->sculpt->bm)) {
                        DAG_id_tag_update(&ob->id, OB_RECALC_DATA);
+               }
 
                if (me->flag & ME_SCULPT_DYNAMIC_TOPOLOGY) {
                        /* Dynamic topology must be disabled before exiting sculpt
                         * mode to ensure the undo stack stays in a consistent
                         * state */
                        sculpt_dynamic_topology_toggle_exec(C, NULL);
+
+                       /* store so we know to re-enable when entering sculpt mode */
+                       me->flag |= ME_SCULPT_DYNAMIC_TOPOLOGY;
                }
 
                /* Leave sculptmode */
@@ -5146,12 +8507,6 @@ static int sculpt_mode_toggle_exec(bContext *C, wmOperator *op)
                /* Enter sculptmode */
                ob->mode |= mode_flag;
 
-               /* Remove dynamic-topology flag; this will be enabled if the
-                * file was saved with dynamic topology on, but we don't
-                * automatically re-enter dynamic-topology mode when loading a
-                * file. */
-               me->flag &= ~ME_SCULPT_DYNAMIC_TOPOLOGY;
-
                if (flush_recalc)
                        DAG_id_tag_update(&ob->id, OB_RECALC_DATA);
 
@@ -5172,7 +8527,7 @@ static int sculpt_mode_toggle_exec(bContext *C, wmOperator *op)
                if (!ts->sculpt->detail_percent)
                        ts->sculpt->detail_percent = 25;
                if (ts->sculpt->constant_detail == 0.0f)
-                       ts->sculpt->constant_detail = 30.0f;
+                       ts->sculpt->constant_detail = 3.0f;
 
                /* Set sane default tiling offsets */
                if (!ts->sculpt->paint.tile_offset[0]) ts->sculpt->paint.tile_offset[0] = 1.0f;
@@ -5205,6 +8560,52 @@ static int sculpt_mode_toggle_exec(bContext *C, wmOperator *op)
                BKE_paint_init(scene, ePaintSculpt, PAINT_CURSOR_SCULPT);
 
                paint_cursor_start(C, sculpt_poll_view3d);
+
+               /* Check dynamic-topology flag; re-enter dynamic-topology mode when changing modes,
+                * As long as no data was added that is not supported. */
+               if (me->flag & ME_SCULPT_DYNAMIC_TOPOLOGY) {
+                       const char *message_unsupported = NULL;
+                       if (me->totloop != me->totpoly * 3) {
+                               message_unsupported = TIP_("non-triangle face");
+                       }
+                       else if (mmd != NULL) {
+                               message_unsupported = TIP_("multi-res modifier");
+                       }
+                       else {
+                               enum eDynTopoWarnFlag flag = sculpt_dynamic_topology_check(C);
+                               if (flag == 0) {
+                                       /* pass */
+                               }
+                               else if (flag & DYNTOPO_WARN_VDATA) {
+                                       message_unsupported = TIP_("vertex data");
+                               }
+                               else if (flag & DYNTOPO_WARN_EDATA) {
+                                       message_unsupported = TIP_("edge data");
+                               }
+                               else if (flag & DYNTOPO_WARN_LDATA) {
+                                       message_unsupported = TIP_("face data");
+                               }
+                               else if (flag & DYNTOPO_WARN_MODIFIER) {
+                                       message_unsupported = TIP_("constructive modifier");
+                               }
+                               else {
+                                       BLI_assert(0);
+                               }
+                       }
+
+                       if (message_unsupported == NULL) {
+                               /* undo push is needed to prevent memory leak */
+                               sculpt_undo_push_begin("Dynamic topology enable");
+                               sculpt_dynamic_topology_enable(C);
+                               sculpt_undo_push_node(ob, NULL, SCULPT_UNDO_DYNTOPO_BEGIN);
+                       }
+                       else {
+                               BKE_reportf(op->reports, RPT_WARNING,
+                                           "Dynamic Topology found: %s, disabled",
+                                           message_unsupported);
+                               me->flag &= ~ME_SCULPT_DYNAMIC_TOPOLOGY;
+                       }
+               }
        }
 
        if (ob->derivedFinal) /* VBO no longer valid */
@@ -5263,7 +8664,7 @@ static int sculpt_detail_flood_fill_exec(bContext *C, wmOperator *UNUSED(op))
        size = max_fff(bb_max[0], bb_max[1], bb_max[2]);
 
        /* update topology size */
-       BKE_pbvh_bmesh_detail_size_set(ss->pbvh, sd->constant_detail / 100.0f);
+       BKE_pbvh_bmesh_detail_size_set(ss->pbvh, 1.0f / sd->constant_detail);
 
        sculpt_undo_push_begin("Dynamic topology flood fill");
        sculpt_undo_push_node(ob, NULL, SCULPT_UNDO_COORDS);
@@ -5328,7 +8729,8 @@ static void sample_detail(bContext *C, int ss_co[2])
                         ray_start, ray_normal, false);
 
        if (srd.hit) {
-               sd->constant_detail = srd.detail * 100.0f;
+               /* convert edge length to detail resolution */
+               sd->constant_detail = 1.0f / srd.detail;
        }
 }
 
@@ -5398,7 +8800,8 @@ static void SCULPT_OT_sample_detail_size(wmOperatorType *ot)
 
        ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
 
-       RNA_def_int_array(ot->srna, "location", 2, NULL, 0, SHRT_MAX, "Location", "Screen Coordinates of sampling", 0, SHRT_MAX);
+       RNA_def_int_array(ot->srna, "location", 2, NULL, 0, SHRT_MAX,
+                         "Location", "Screen Coordinates of sampling", 0, SHRT_MAX);
 }
 
 
@@ -5412,11 +8815,11 @@ static int sculpt_set_detail_size_exec(bContext *C, wmOperator *UNUSED(op))
        WM_operator_properties_create_ptr(&props_ptr, ot);
 
        if (sd->flags & SCULPT_DYNTOPO_DETAIL_CONSTANT) {
-               set_brush_rc_props(&props_ptr, "sculpt", "constant_detail", NULL, 0);
-               RNA_string_set(&props_ptr, "data_path_primary", "tool_settings.sculpt.constant_detail");
+               set_brush_rc_props(&props_ptr, "sculpt", "constant_detail_resolution", NULL, 0);
+               RNA_string_set(&props_ptr, "data_path_primary", "tool_settings.sculpt.constant_detail_resolution");
        }
        else if (sd->flags & SCULPT_DYNTOPO_DETAIL_BRUSH) {
-               set_brush_rc_props(&props_ptr, "sculpt", "constant_detail", NULL, 0);
+               set_brush_rc_props(&props_ptr, "sculpt", "constant_detail_resolution", NULL, 0);
                RNA_string_set(&props_ptr, "data_path_primary", "tool_settings.sculpt.detail_percent");
        }
        else {
@@ -5450,6 +8853,7 @@ void ED_operatortypes_sculpt(void)
        WM_operatortype_append(SCULPT_OT_brush_stroke);
        WM_operatortype_append(SCULPT_OT_sculptmode_toggle);
        WM_operatortype_append(SCULPT_OT_set_persistent_base);
+       WM_operatortype_append(SCULPT_OT_silhouette_draw);
        WM_operatortype_append(SCULPT_OT_dynamic_topology_toggle);
        WM_operatortype_append(SCULPT_OT_optimize);
        WM_operatortype_append(SCULPT_OT_symmetrize);