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 9fb6a514ee7908901f3ccdc15ce7d950b5171425..8f2d0242ea0112bafa524d7b9ff148d892ec01dc 100644 (file)
 #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);
+/* 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;
@@ -137,6 +150,117 @@ static void bl_debug_draw_point(const float pos[3],const float thickness)
 }
 #endif
 
+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
+
+
+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
  *
  * Avoid duplicate checks, internal logic only,
@@ -543,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;
@@ -4744,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;
@@ -5019,39 +5149,36 @@ static void SCULPT_OT_set_persistent_base(wmOperatorType *ot)
 
 /****************** Topology tools Silhouette ******************/
 
-/*     init data:
- *     Silhouette Data */
-typedef struct SilhouetteStroke {
-       float *points;
-       int totvert;
-       int max_verts;
-}SilhouetteStroke;
+/* Axis-aligned bounding box (Same as in pbvh_intern.h)*/
 
-typedef enum {
-       SIL_INIT = 0,
-       SIL_DRAWING = 1,
-       SIL_OP = 2
-} SilhouetteState;
+/* TODO: import BB functions */
+/* Expand the bounding box to include a new coordinate */
+static void BB_expand(BB *bb, const float co[3])
+{
+       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]);
+       }
+}
 
-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;
+static void BB_reset(BB *bb)
+{
+       bb->bmin[0] = bb->bmin[1] = bb->bmin[2] = FLT_MAX;
+       bb->bmax[0] = bb->bmax[1] = bb->bmax[2] = -FLT_MAX;
+}
 
-       float add_col[3];               /* preview color */
-       float last_mouse_pos[2];
+static bool bb_intersect(BB *bb1, BB *bb2) {
+       int i;
 
-       SilhouetteState state;  /* Operator state */
+       /* 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;
+               }
+       }
 
-       float depth;                    /* Depth or thickness of the generated shape */
-       float smoothness;               /* Smoothness of the generated shape */
-       float anchor[3];                /* Origin point of the reference plane */
-       float z_vec[3];                 /* Orientation of the reference plane */
-} SilhouetteData;
+       return true;
+}
 
 static void silhouette_stroke_free(SilhouetteStroke *stroke)
 {
@@ -5063,23 +5190,28 @@ static void silhouette_stroke_free(SilhouetteStroke *stroke)
        }
 }
 
-static SilhouetteStroke *silhouette_stroke_new(int max_verts)
+static SilhouetteStroke *silhouette_stroke_new()
 {
        SilhouetteStroke *stroke = MEM_callocN(sizeof(SilhouetteStroke), "SilhouetteStroke");
-       stroke->points = 0;
-       stroke->points = MEM_callocN(sizeof(float) * 3 * max_verts,"SilhouetteStrokePoints");/* TODO: Dynamic length */
+       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 = max_verts;
+       stroke->max_verts = SIL_STROKE_STORE_CHUNK;
+       BB_reset(&stroke->bb);
        return stroke;
 }
 
-static SilhouetteData *silhouette_data_new(bContext *C, wmOperator *op, bool rna_full)
+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);
+
        sil->ar = CTX_wm_region(C);
-       sil->current_stroke = silhouette_stroke_new(1024);
+       sil->current_stroke = silhouette_stroke_new();
        view3d_set_viewcontext(C, &sil->vc);
 
        sil->add_col[0] = 1.00; /* add mode color is light red */
@@ -5087,14 +5219,19 @@ static SilhouetteData *silhouette_data_new(bContext *C, wmOperator *op, bool rna
        sil->add_col[2] = 0.39;
 
        /*Load RNA Data if present */
-       sil->smoothness = RNA_float_get(op->ptr, "smoothness");
-       sil->depth = RNA_float_get(op->ptr, "depth");
-       if (rna_full) {
-               RNA_float_get_array(op->ptr, "z_vec", sil->z_vec);
-               RNA_float_get_array(op->ptr, "anchor", sil->anchor);
-               RNA_float_get_array(op->ptr, "points", sil->current_stroke->points);
-               sil->current_stroke->totvert = RNA_int_get(op->ptr, "totvert");
-       }
+       sil->smoothness = sd->silhouette_smoothness / 100.0f;
+       sil->depth = sd->silhouette_depth;
+       sil->resolution = sd->silhouette_resolution;
+
+       copy_v3_v3(sil->anchor, fp);
+
+       /* 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;
@@ -5108,6 +5245,25 @@ static void silhouette_data_free(struct wmOperator *op)
        data = op->customdata;
        if (data) {
                silhouette_stroke_free(data->current_stroke);
+
+               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);
        }
 }
@@ -5118,54 +5274,62 @@ static void silhoute_stroke_point_to_3d(SilhouetteData *sil, int point, float r_
        /*ED_view3d_win_to_3d(sil->vc.v3d, sil->ar, sil->anchor, &sil->current_stroke->points[point], r_v);*/
 }
 
+#if 0
 /* TODO: Add dynamic memory allocation */
 static void silhouette_stroke_add_3Dpoint(SilhouetteStroke *stroke, float point[3])
 {
-       if (stroke->totvert < stroke->max_verts) {
-               copy_v3_v3(&stroke->points[stroke->totvert * 3], point);
-               stroke->totvert ++;
-       } else {
-               printf("Stroke reached maximum vert count.\n");
+       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_v3_v3(&stroke->points[stroke->totvert * 3], point);
+       stroke->totvert ++;
 }
+#endif
 
 static void silhouette_stroke_add_point(SilhouetteData *sil, SilhouetteStroke *stroke, float point[2])
 {
-       if (stroke->totvert < stroke->max_verts) {
-               ED_view3d_win_to_3d(sil->vc.v3d, sil->ar, sil->anchor, point, &sil->current_stroke->points[stroke->totvert * 3]);
-               stroke->totvert ++;
-       } else {
-               printf("Stroke reached maximum vert count.\n");
+       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(const bContext *C, SilhouetteData *sil)
+static void silhouette_set_ref_plane(SilhouetteData *sil)
 {
-       Scene *scene = CTX_data_scene(C);
-       View3D *v3d = CTX_wm_view3d(C);
-       const float *fp = ED_view3d_cursor3d_get(scene, v3d);
-
-       ED_view3d_global_to_vector(sil->ar->regiondata, (float[3]){0.0f,0.0f,0.0f}, sil->z_vec);
-       copy_v3_v3(sil->anchor, fp);
+       /*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(bContext *C, float mouse[2], SilhouetteData *sil)
+static void sculpt_silhouette_stroke_update(float mouse[2], SilhouetteData *sil)
 {
-       silhouette_set_ref_plane( C, 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_VERT_GEN = 1,    /* Verts on the branch ends are stored */
-       BRANCH_EDGE_GEN = 2             /* Edges on the ends are stored (used primarly for bridging) */
+       BRANCH_EDGE_GEN = 1     /* Edges on the ends are stored (used primarly for bridging) */
 } BranchState;
 
 
@@ -5188,8 +5352,10 @@ typedef struct SpineBranch{
        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 *gen_verts;                 /* Verts on the branch ends are stored */
        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 */
@@ -5198,6 +5364,31 @@ typedef struct Spine{
        SpineBranch **branches; /* All branches. Can contain Null pointers if branches got removed*/
 }Spine;
 
+typedef struct {
+       SculptSession *ss;
+       BB *bb_target;
+       bool original;
+} SculptSearchBBData;
+
+/* 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;
+
+       BKE_pbvh_node_get_BB(node, bb_min, bb_max);
+
+       /* 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;
+}
+
 static int get_adjacent_faces(BMFace *f, BMFace **ad_f, BMFace *last_f)
 {
        int ad_faces = 0;
@@ -5288,11 +5479,9 @@ static void free_spine_branch(SpineBranch *branch)
        if (branch->terminal_points) {
                MEM_freeN(branch->terminal_points);
        }
-       if (branch->flag & BRANCH_VERT_GEN) {
-               MEM_freeN(branch->gen_verts);
-       }
        if (branch->flag & BRANCH_EDGE_GEN) {
                MEM_freeN(branch->e_start_arr);
+               MEM_freeN(branch->e_flip_side_ends);
        }
        MEM_freeN(branch);
 }
@@ -5344,7 +5533,7 @@ static SpineBranch *new_spine_branch(int idx, int max_alloc, int hull_max)
        /*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) * hull_max * 2, "Spine terminal");
+       branch->terminal_points = MEM_callocN(sizeof(int) * 3 * 2, "Spine terminal");
 
        branch->idx = idx;
        return branch;
@@ -5374,7 +5563,7 @@ static SpineBranch *spine_branchoff(Spine *spine, SpineBranch *current_branch, i
        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;
+       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]);
@@ -5533,6 +5722,95 @@ static Spine *silhouette_generate_spine(SilhouetteData *sil, SilhouetteStroke *s
        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)
 {
@@ -5635,8 +5913,13 @@ static void bridge_loops(Mesh *me, int e_start_a, int e_start_b, int totvert, bo
                me->medge[e_start + i].flag = 0;
 
                if (i < totvert - 1) {
-                       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 + (flip ? -i : i) * b_stride].v1;
+                       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;
@@ -5644,10 +5927,16 @@ static void bridge_loops(Mesh *me, int e_start_a, int e_start_b, int totvert, bo
                        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;
 
-                       me->mloop[l_start + i * 4 + 2].v = me->medge[e_start_b + (flip ? -i : i) * b_stride].v2;
-                       me->mloop[l_start + i * 4 + 2].e = e_start_b + (flip ? -i : i) * b_stride;
+                       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)].v = me->medge[e_start_b + (flip ? -i : 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;
@@ -5656,49 +5945,126 @@ static void bridge_loops(Mesh *me, int e_start_a, int e_start_b, int totvert, bo
                        me->mpoly[p_start + i].flag = 0;
                        me->mpoly[p_start + i].pad = 0;
                } 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 + (flip ? -(i - 1) : (i - 1)) * b_stride].v2;
+                       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;
+                       }
                }
        }
 }
 
-/* TODO: Redundant */
-static void bridge_loop_flip(Mesh *me, int v_start_a, int v_start_b, int totvert, bool flip, int l_start_a, int l_start_b, int l_stride)
+/*
+ * 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;
-       ED_mesh_edges_add(me, NULL, totvert);
-       ED_mesh_loops_add(me, NULL, 4 * totvert - 4);
-       ED_mesh_polys_add(me, NULL, totvert - 1);
+       bool b_flip = false;
 
-       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;
-               me->medge[e_start + i].v1 = v_start_a + i;
-               me->medge[e_start + i].v2 = v_start_b + (flip ? -i : i);
+       e_a = me->medge[edge_a];
+       e_b = me->medge[edge_b];
+       e_c = me->medge[edge_c];
 
-               if (i < totvert - 1) {
-                       me->mloop[l_start + i * 4 + 0].v = v_start_a + i;
-                       me->mloop[l_start + i * 4 + 0].e = l_start_a + i * l_stride;
+       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;
+       }
 
-                       me->mloop[l_start + i * 4 + 1].v = v_start_a + i + 1;
-                       me->mloop[l_start + i * 4 + 1].e = e_start + i + 1;
+       ED_mesh_edges_add(me, NULL, 1);
 
-                       me->mloop[l_start + i * 4 + 2].v = v_start_b + (flip ? -i - 1 : i + 1);
-                       me->mloop[l_start + i * 4 + 2].e = l_start_b + (flip ? -i : i) * l_stride;
+       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;
 
-                       me->mloop[l_start + i * 4 + 3].v = v_start_b + (flip ? -i : i);
-                       me->mloop[l_start + i * 4 + 3].e = e_start + i;
+       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->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;
-               }
+       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?*/
@@ -5708,162 +6074,77 @@ static int cmpfunc (const void * a, const void * 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 ss_steps, int w_steps, float w_fact, int *e_start_sides)
+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)
 {
-       int v_steps = ss_steps;
-       int square_ss_steps = u_steps * v_steps;
        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], z_vec_b[3];
-       int l_u_pos_i = 1, r_u_pos_i = totr - 1;
-       int v_pos_w = 0;
+       float v1[3], v2[3], v3[3], v4[3];
+       int l_u_pos_i = 1, r_u_pos_i = totr - 2;
 
-       int ve_start_a, ve_start_ar, ve_start_b;
        const int v_start = me->totvert;
-       MVert ref_v;
-       ref_v.flag = 0; ref_v.bweight = 0;
-       ED_mesh_vertices_add(me, NULL, square_ss_steps * 2 + (w_steps - 1) * u_steps);
+       const int e_start = me->totedge;
 
        for (int u = 0; u < u_steps; u++) {
-               while (left[l_u_pos_i * 4 + 3] <= step_l * (float)u && l_u_pos_i < totl) {
+               while (l_u_pos_i < totl - 1 && left[l_u_pos_i * 4 + 3] <= step_l * (float)u) {
                        l_u_pos_i ++;
                }
 
-               while (right[r_u_pos_i * 4 + 3] > step_r * (float)(u_steps - u - 1) && r_u_pos_i > 0) {
+               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 */
-               a = left[l_u_pos_i * 4 - 1];
-               b = left[l_u_pos_i * 4 + 3];
-               f = (step_l * (float)u - a) / (b - a);
-               interp_v3_v3v3(v1, &left[l_u_pos_i * 4 - 4], &left[l_u_pos_i * 4], f);
-
-               a = right[r_u_pos_i * 4 + 7];
-               b = right[r_u_pos_i * 4 + 3];
-               f = (step_r * (float)(u_steps - u - 1) - a) / (b - a);
-               interp_v3_v3v3(v2, &right[r_u_pos_i * 4 + 4], &right[r_u_pos_i * 4], f);
-
-               copy_v3_v3(z_vec_b, z_vec);
-
-               /*TODO: Different interpolation, different shapes, smoothness etc*/
-               /* Interpolate point position in z-direction. Smoothing etc is done here. Redundant?*/
-               for (int v = 0; v < ss_steps; v++) {
-                       add_v3_v3(v1, z_vec_b);
-                       add_v3_v3(v2, z_vec_b);
-                       f = (float)v/(float)(v_steps * 2);
-                       f = (sin(M_PI * f - M_PI * 0.5f) * 0.5f + 0.5f) * w_fact;
-                       interp_v3_v3v3(v3, v1, v2, f);
-                       me->mvert[v_start + u * v_steps + v].flag = 0;
-                       me->mvert[v_start + u * v_steps + v].bweight = 0;
-                       copy_v3_v3(me->mvert[v_start + u * v_steps + v].co,v3);
-                       interp_v3_v3v3(v4, v2, v1, f);
-                       me->mvert[v_start + u * v_steps + v + square_ss_steps].flag = 0;
-                       me->mvert[v_start + u * v_steps + v + square_ss_steps].bweight = 0;
-                       copy_v3_v3(me->mvert[v_start + u * v_steps + v + square_ss_steps].co,v4);
-                       mul_v3_fl(z_vec_b, cos(0.5f * M_PI * ((float)v / (float)v_steps)));
-                       if (v == ss_steps - 1) {
-                               for (int w = 1; w < w_steps; w++) {
-                                       f = (float)w / (float)w_steps;
-                                       interp_v3_v3v3(v1, v3, v4, f);
-                                       v_pos_w = v_start + ss_steps * u_steps * 2 + u * (w_steps - 1) + (w - 1);
-                                       me->mvert[v_pos_w].flag = 0;
-                                       me->mvert[v_pos_w].bweight = 0;
-                                       copy_v3_v3(me->mvert[v_pos_w].co,v1);
-                               }
+               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]);
                }
-       }
-
-       /* Generate the Mesh and store border edges for later */
-       ve_start_a = me->totedge + v_steps * 2 - 2;
-       e_start_sides[0] = me->totedge;
-       e_start_sides[1] = v_steps;
-       e_start_sides[2] = 2;
-       e_start_sides[9] = me->totedge + (v_steps * 2 - 1) * (u_steps - 1);
-       e_start_sides[10] = v_steps;
-       e_start_sides[11] = 1;
-       generate_mesh_grid_f_e(me, u_steps, v_steps, v_start, false);
-
-       ve_start_ar = me->totedge + v_steps * 2 - 2;
-       e_start_sides[6] = me->totedge;
-       e_start_sides[7] = v_steps;
-       e_start_sides[8] = 2;
-       e_start_sides[15] = me->totedge + (v_steps * 2 - 1) * (u_steps - 1);
-       e_start_sides[16] = v_steps;
-       e_start_sides[17] = 1;
-       generate_mesh_grid_f_e(me, u_steps, v_steps, v_start + square_ss_steps, true);
-
-       ve_start_b = me->totedge + 1;
-       e_start_sides[3] = me->totedge;
-       e_start_sides[4] = w_steps - 1;
-       e_start_sides[5] = 2;
-       e_start_sides[12] = me->totedge + ((w_steps - 1) * 2 - 1) * (u_steps - 1);//(w_steps * 2 - 1) * (u_steps - 1);
-       e_start_sides[13] = w_steps - 1;
-       e_start_sides[14] = 1;
-       generate_mesh_grid_f_e(me, u_steps, w_steps - 1, v_start + square_ss_steps * 2, false);
-
-       bridge_loops(me,
-                                ve_start_a,
-                                ve_start_b,
-                                u_steps,
-                                false,
-                                (v_steps * 2 - 1),
-                                (w_steps * 2 - 3),
-                                false);
 
-       bridge_loops(me,
-                                ve_start_ar,
-                                ve_start_b + (w_steps * 2 - 3) - 2,
-                                u_steps,
-                                false,
-                                (v_steps * 2 - 1),
-                                (w_steps * 2 - 3),
-                                true);
-}
-
-/* Calculate an Arc between two points with additional interpolation in z direction */
-static void calc_fork_arc(float v1[3], float v2[3], float *out, int offloop, const float z_vec[3], int ss_steps, int w_steps, float w_fact)
-{
-       float v3[3], v4[3], z_vec_b[3];
-       float f;
-
-       copy_v3_v3(z_vec_b, z_vec);
-       for (int v = 0; v < ss_steps; v++) {
-               add_v3_v3(v1, z_vec_b);
-               add_v3_v3(v2, z_vec_b);
-               f = (float)v/(float)(ss_steps * 2);
-               f = (sin(M_PI * f - M_PI * 0.5f) * 0.5f + 0.5f) * w_fact;
-               interp_v3_v3v3(v3, v1, v2, f);
-               copy_v3_v3(&out[offloop + v * 3], v3);
-               interp_v3_v3v3(v4, v2, v1, f);
-               copy_v3_v3(&out[offloop + (ss_steps + w_steps - 1 + v) * 3], v4);
-               mul_v3_fl(z_vec_b, cos(0.5f * M_PI * ((float)v / (float)ss_steps)));
-               if (v == ss_steps - 1) {
-                       for (int w = 1; w < w_steps; w++) {
-                               f = (float)w / (float)w_steps;
-                               interp_v3_v3v3(v1, v3, v4, f);
-                               copy_v3_v3(&out[offloop + (ss_steps + w - 1) * 3], v1);
+               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;
        }
 }
 
-/* Generate the Cap for branches with one ends. */
-static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts, Mesh *me, float z_vec[3], float depth, int ss_steps, int w_steps, float w_fact)
+static int get_cyclic_offset(SpineBranch *branch)
 {
-       float cap_p[branch->tot_hull_points * 4];
-       float v1[3], v2[3], v3[3], z_vec_b[3];
-       float totlength = 0.0f;
-       float step_size = depth / (float)w_steps;
-       float loop[3 * (ss_steps * 2 + w_steps)];
        int cyclic_offset = 0, n_i = 0;
-       gen_verts = NULL; /*TODO: Remove gen_verts */
-
-       /* 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 > 0) {
@@ -5878,10 +6159,36 @@ static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts,
                        }
                }
        }
+       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);
 
-       printf("New Cap:\n");
        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];
@@ -5889,17 +6196,16 @@ static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts,
                        cap_p[i * 4 + 3] = 0.0f;
                }
                copy_v3_v3(v1, &cap_p[i * 4]);
-               printf("Added Point %i; length: %f\n", branch->hull_points[n_i % branch->tot_hull_points], cap_p[i * 4 + 3]);
        }
 
-       totlength = cap_p[branch->tot_hull_points * 4 - 1];
-
-       if (!gen_verts) {
-               copy_v3_v3(v1, &cap_p[0]);
-               copy_v3_v3(v2, &cap_p[branch->tot_hull_points * 4 - 4]);
-               calc_fork_arc(v1, v2, loop, 0, z_vec, ss_steps, w_steps + 1, w_fact);
+       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;
@@ -5907,18 +6213,21 @@ static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts,
        int u_pos_i = 0;
        int v_start = me->totvert;
        int u_steps;
-       int ve_start_a = -1;
-       int e_start_tube[18];
+       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) {
-               printf("Gen Tube before cap\n");
                int side_l = cap_pos;
-               int totl = 0, totr = 0;
-               float left[branch->tot_hull_points * 4], right[branch->tot_hull_points * 4];
+               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];
@@ -5933,6 +6242,7 @@ static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts,
                        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];
@@ -5942,129 +6252,294 @@ static void add_ss_cap(SilhouetteData *sil, SpineBranch *branch, int *gen_verts,
                }
 
                if (totl >= 1 && totr >= 1) {
-                       u_steps = fmax(2.0f, fmax(left[totl * 4 - 1], right[totr * 4 - 1]) / (float)(2 * depth / ss_steps));
-                       ve_start_a = me->totedge + ss_steps * (u_steps - 1) * 2 - 2;
-                       fill_tube(me, left, right, totl, totr, u_steps, z_vec, ss_steps, w_steps + 1, w_fact, e_start_tube);
-                       for (int u = 0; u < w_steps; u++) {
-                               if (!gen_verts) {
-                                       copy_v3_v3(&loop[(ss_steps + u) * 3], me->mvert[me->totvert - w_steps + u].co);
-                               } else {
-                                       gen_verts[ss_steps + u] = me->totvert - w_steps + u;
-                               }
+                       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);
        }
 
-       u_pos_i = 1;
+       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;
 
-       ED_mesh_vertices_add(me, NULL, w_steps * ss_steps);
+       u_pos_i = 1;
+       v_start = me->totvert;
 
-       /*TODO: Corners need two additional loops*/
        for (int u = 0; u < w_steps; u++) {
-               while (cap_p[u_pos_i * 4 + 3] <= cap_pos) {
+               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];
-               f = (cap_pos - a) / (b - a);
-               interp_v3_v3v3(v1, &cap_p[u_pos_i * 4 - 4], &cap_p[u_pos_i * 4], f);
-               copy_v3_v3(z_vec_b, z_vec);
-
-               for (int v = 0; v < ss_steps; v++) {
-                       add_v3_v3(v1, z_vec_b);
-                       f = (float)v / (float)(ss_steps * 2);
-                       f = sin(M_PI * f - M_PI * 0.5f) + 1.0f;
-                       if (!gen_verts) {
-                               interp_v3_v3v3(v3, v1, &loop[(ss_steps + u) * 3], f * w_fact);
-                       } else {
-                               interp_v3_v3v3(v3, v1, me->mvert[gen_verts[ss_steps + u]].co, f * w_fact);
-                       }
-                       me->mvert[v_start + u * ss_steps + v].flag = 0;
-                       me->mvert[v_start + u * ss_steps + v].bweight = 0;
-                       copy_v3_v3(me->mvert[v_start + u * ss_steps + v].co,v3);
-                       mul_v3_fl(z_vec_b, cos(0.5f * M_PI * ((float)v / (float)ss_steps)));
+               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;
 
-       int ve_start_b = me->totedge;
-       generate_mesh_grid_f_e(me, w_steps, ss_steps, v_start, false);
-       if (ve_start_a > 0 && e_start_tube) {
+       if (flip_side) {
+               cap_end_flip_e_start = me->totedge;
                bridge_loops(me,
-                                e_start_tube[9],
-                                ve_start_b,
-                                e_start_tube[10],
+                                        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,
-                                true);
-               bridge_loops(me,
-                                        e_start_tube[15],
-                                        ve_start_b + (ss_steps * 2 - 1) * (w_steps - 1),
-                                        e_start_tube[16],
-                                        false,
-                                        1,
-                                        1,
-                                        false);
-               bridge_loops(me,
-                                        e_start_tube[12],
-                                        ve_start_b + (ss_steps * 2 - 1) - 1,
-                                        e_start_tube[13],
-                                        false,
-                                        1,
-                                        ss_steps * 2 - 1,
-                                        true);
+                                !n_g_flip);
 
-               if (!(branch->flag & BRANCH_EDGE_GEN)) {
-                       branch->e_start_arr = MEM_callocN(sizeof(int) * 9,"edge startposition array");
-                       branch->flag |= BRANCH_EDGE_GEN;
-               }
-               memcpy(branch->e_start_arr, e_start_tube, 9 * sizeof(int));
-               /*Close corners!*/
-       } else {
-               if (!(branch->flag & BRANCH_EDGE_GEN)){
-                       branch->e_start_arr = MEM_callocN(sizeof(int) * 9,"edge startposition array");
-                       branch->flag |= BRANCH_EDGE_GEN;
-               }
+       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;
 
-               if (branch->hull_points[0] == 0) {
-                       branch->e_start_arr[6] = ve_start_b;
-                       branch->e_start_arr[7] = ss_steps;
-                       branch->e_start_arr[8] = 2;
+       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);
 
-                       branch->e_start_arr[0] = ve_start_b + (ss_steps * 2 - 1) * (w_steps - 1);
-                       branch->e_start_arr[1] = ss_steps;
-                       branch->e_start_arr[2] = 1;
+       e_corner_b = me->totedge - 1;
 
-                       /* TODO: Inverse fix?*/
-                       branch->e_start_arr[3] = ve_start_b + (ss_steps * 2 - 1) - 1;
-                       branch->e_start_arr[4] = w_steps;
-                       branch->e_start_arr[5] = ss_steps * 2 - 1;
-               } else {
-                       branch->e_start_arr[0] = ve_start_b;
-                       branch->e_start_arr[1] = ss_steps;
-                       branch->e_start_arr[2] = 2;
+       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];
+       }
 
-                       branch->e_start_arr[6] = ve_start_b + (ss_steps * 2 - 1) * (w_steps - 1);
-                       branch->e_start_arr[7] = ss_steps;
-                       branch->e_start_arr[8] = 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);
 
-                       branch->e_start_arr[3] = ve_start_b + (ss_steps * 2 - 1) - 1;
-                       branch->e_start_arr[4] = w_steps;
-                       branch->e_start_arr[5] = ss_steps * 2 - 1;
+                       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, SpineBranch *branch, Mesh *me, float z_vec[3], int ss_steps, int w_steps, float w_fact)
+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[branch->tot_hull_points * 4 * 3];
+       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);
@@ -6088,6 +6563,7 @@ static void add_ss_tinter(SilhouetteData *sil, SpineBranch *branch, Mesh *me, fl
        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;
@@ -6106,195 +6582,279 @@ static void add_ss_tinter(SilhouetteData *sil, SpineBranch *branch, Mesh *me, fl
                }
        }
 
-       int loop_length = (ss_steps * 2 + w_steps);
-       /* a - b */
-       float loop[3 * loop_length * 3];
-       float v1[3], v2[3], v3[3];
-
-       copy_v3_v3(v1, &sa[b_tot[0] * 4 - 4]);
-       copy_v3_v3(v2, &sa[b_start[1]]);
-       calc_fork_arc(v1, v2, loop, 0, z_vec, ss_steps, w_steps + 1, w_fact);
-
-       /* b - c */
-       copy_v3_v3(v1, &sa[b_start[1] + b_tot[1] * 4 - 4]);
-       copy_v3_v3(v2, &sa[b_start[2]]);
-       calc_fork_arc(v1, v2, loop, loop_length * 3, z_vec, ss_steps, w_steps + 1, w_fact);
-
-       /* c - a */
-       copy_v3_v3(v1, &sa[b_start[2] + b_tot[2] * 4 - 4]);
-       copy_v3_v3(v2, &sa[0]);
-       calc_fork_arc(v1, v2, loop, loop_length * 3 * 2, z_vec, ss_steps, w_steps + 1, w_fact);
-
-       float center[3], ind_center[3 * 3];
-       int center_idx = (ss_steps + (int)floorf(w_steps * 0.5f)) * 3;
+       u_steps = 5;
+       zero_v3(center);
        for (int s = 0; s < 3; s++) {
-               copy_v3_v3(&ind_center[s * 3], &loop[s * loop_length * 3 + center_idx]);
-               add_v3_v3(&ind_center[s * 3], &loop[s * loop_length * 3 + center_idx - 3]);
-               mul_v3_fl(&ind_center[s * 3], 0.5f);
+               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]);
        }
-
-       copy_v3_v3(center, ind_center);
-       add_v3_v3(center, ind_center + 3);
-       add_v3_v3(center, ind_center + 6);
        mul_v3_fl(center, 1.0f/3.0f);
+       add_v3_v3v3(center_up, center, z_vec);
 
-       int v_start = me->totvert;
-       int v_center_guides = v_start;
-       int u_steps = 9; /* TODO: Calc from U_steps and length */
-       int verts_per_side = w_steps * u_steps + w_steps;
-       float f;
-       int e_start = me->totedge;
-       int bridge_in_out_start[3];
+       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));
+       }
 
-       if (!(branch->flag & BRANCH_EDGE_GEN)){
-               branch->e_start_arr = MEM_callocN(sizeof(int) * 27,"edge startposition array");
+       /*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;
        }
 
-       ED_mesh_vertices_add(me, NULL, 3 * verts_per_side);
-       int center_face[3 * 2];
+       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);
 
-       /* Generate the topside Geometry. (3 Sides * w_steps / 2 * u_steps + One Tri in the middle)
-        * Interpolate for every side between three points: The two sides and the centerpoint */
        for (int s = 0; s < 3; s++) {
-               for (int w = 0; w < w_steps/2; w++) {
-                       copy_v3_v3(v1, &loop[s * loop_length * 3 + (ss_steps + w_steps - w - 1) * 3]);
-                       copy_v3_v3(v2, &loop[((s + 1) % 3) * loop_length * 3 + (ss_steps + w) * 3]);
-                       for (int u = 0; u < u_steps; u++) {
-                               f = (float)u / (float)(u_steps - 1);
-                               interp_v3_v3v3(v3, v1, v2, f);
-                               f = -2.0 * ((f - 0.5f) * (f - 0.5f)) + 0.5f;
-                               interp_v3_v3v3(v3, v3, center, f / (w_steps/2 - w + 1));
-                               copy_v3_v3(me->mvert[v_start + s * verts_per_side + w * u_steps + u].co, v3);
-
-                               if (u == u_steps / 2 && w == w_steps / 2 - 1) {
-                                       center_face[s * 2] = v_start + s * verts_per_side + w * u_steps + u;
-                               }
-                       }
-               }
-               int e_start_off = me->totedge;
-               generate_mesh_grid_f_e(me, w_steps / 2, u_steps, v_start + s * verts_per_side, true);
-               branch->e_start_arr[s * 9 + 3] = e_start_off + 1;
-               branch->e_start_arr[s * 9 + 4] = w_steps / 2;
-               branch->e_start_arr[s * 9 + 5] = (u_steps - 1) * 2 + 1;
+               step_length = sa[b_start[s] + b_tot[s] * 4 - 1] / (float)u_steps;
 
-               int e_start_tot_ps = (u_steps * (w_steps / 2)) * 2 - u_steps - (w_steps / 2);
+               add_v3_v3v3(v3, &center_s[s * 3], z_vec);
 
-               bridge_in_out_start[s] = e_start + s * e_start_tot_ps;
-       }
+               v_start = me->totvert;
 
-       int edges_per_side = u_steps / 2 + 1;
-       for (int s = 0; s < 3; s++) {
-               /* Todo: redundant */
-               bridge_loop_flip(me,
-                                       v_start + s * verts_per_side + u_steps * (w_steps / 2 - 1),
-                                       v_start + ((s + 2) % 3) * verts_per_side + u_steps * (w_steps / 2) - 1,
-                                       u_steps / 2 + 1,
-                                       true,
-                                       e_start + (s + 1) * edges_per_side - u_steps + 1,
-                                       e_start + ((s + 2) % 3 + 1) * edges_per_side - 1,
-                                       1);
-               center_face[s * 2 + 1] = me->totedge - 1;
-       }
-
-       ED_mesh_loops_add(me, NULL, 3);
-       ED_mesh_polys_add(me, NULL, 1);
-       for (int i = 0; i < 3; i++) {
-               me->mloop[me->totloop - 3 + i].v = center_face[i * 2];
-               me->mloop[me->totloop - 3 + i].e = center_face[i * 2 + 1];
-       }
+               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;
 
-       me->mpoly[me->totpoly - 1].loopstart = me->totloop - 3;
-       me->mpoly[me->totpoly - 1].totloop = 3;
+               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);
 
-       int u_pos_i, e_tmp_start;
-       float step_l, a, b;
-       float z_vec_b[3];
-       verts_per_side = u_steps * ss_steps;
-       v_start = me->totvert;
-       ED_mesh_vertices_add(me, NULL, 3 * verts_per_side);
+               e_start[s] = me->totedge;
 
-       /* Generate the three lower sides (3 Sides * u_steps * ss_steps)*/
-       for (int s = 0; s < 3; s++) {
-               u_pos_i = fmin(1, b_tot[s] - 1);
-               step_l = sa[b_start[s] + b_tot[s] * 4 - 1] / u_steps;
+               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 = 0; u < u_steps; u++) {
-                       while (u_pos_i < b_tot[s] - 1 && sa[b_start[s] + u_pos_i * 4 + 3] <= step_l * (float)u) {
-                               u_pos_i ++;
+               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 ++;
                        }
 
-                       if (u_pos_i == 0 || step_l == 0) {
+                       /* 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] + u_pos_i * 4 - 1];
-                               b = sa[b_start[s] + u_pos_i * 4 + 3];
-                               BLI_assert(b > a);
-                               f = (step_l * (float)u - a) / (b - a);
-                               interp_v3_v3v3(v1, &sa[b_start[s] + u_pos_i * 4 - 4], &sa[b_start[s] + u_pos_i * 4], f);
+                               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);
                        }
 
-                       copy_v3_v3(z_vec_b, z_vec);
-                       sub_v3_v3v3(v2, center, z_vec);
-                       copy_v3_v3(v2, me->mvert[v_center_guides + ((s + 2) % 3) * (w_steps * (u_steps + 1)) + u].co);
-
-                       for (int v = 0; v < ss_steps; v++) {
-                               add_v3_v3(v1, z_vec_b);
-                               f = (float)v/(float)(ss_steps);
-                               f = sin(0.5f * f * M_PI - 0.5f * M_PI) + 1.0f;
-                               interp_v3_v3v3(v3, v1, v2, f);
-                               me->mvert[v_start + verts_per_side * s + u * ss_steps + v].flag = 0;
-                               me->mvert[v_start + verts_per_side * s + u * ss_steps + v].bweight = 0;
-                               copy_v3_v3(me->mvert[v_start + verts_per_side * s + u * ss_steps + v].co,v3);
-                               mul_v3_fl(z_vec_b, cos(0.5f * M_PI * ((float)v / (float)ss_steps)));
+                       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);
                }
-               e_tmp_start = me->totedge;
-               generate_mesh_grid_f_e(me, u_steps, ss_steps, v_start + s * verts_per_side, false);
-               if (s > 0) {
-                       branch->e_start_arr[s * 9 + 6] = e_tmp_start;
-                       branch->e_start_arr[s * 9 + 7] = ss_steps;
-                       branch->e_start_arr[s * 9 + 8] = 2;
-               } else {
-                       branch->e_start_arr[s * 9 + 0] = e_tmp_start;
-                       branch->e_start_arr[s * 9 + 1] = ss_steps;
-                       branch->e_start_arr[s * 9 + 2] = 2;
+
+               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;
                }
 
-               if ((s + 1) % 3 > 0) {
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 0] = e_tmp_start + (u_steps - 1) * ss_steps * 2 - u_steps + 1;
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 1] = ss_steps;
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 2] = 1;
-               } else {
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 6] = e_tmp_start + (u_steps - 1) * ss_steps * 2 - u_steps + 1;
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 7] = ss_steps;
-                       branch->e_start_arr[(s + 1) % 3 * 9 + 8] = 1;
+               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,
-                                        bridge_in_out_start[(s + 2) % 3],
-                                        e_tmp_start + ss_steps * 2 - 2,
-                                        u_steps,
+                                        e_start_inner[s],
+                                        e_start[s],
+                                        v_steps -(flip_side ? 1 : 0) + w_steps / 2,
                                         false,
                                         2,
-                                        ss_steps * 2 -1,
-                                        true);
+                                        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;
        }
 
-       /*printf("Created a T-Intersection. Edges are %i, %i, %i Total length of: %f, %f, %f\n", b_tot[0], b_tot[1], b_tot[2], sa[b_tot[0] * 4 - 1], sa[b_start[1] + b_tot[1] * 4 - 1], sa[b_start[2] + b_tot[2] * 4 - 1]);*/
+       BLI_array_free(sa);
 }
 
-static void add_ss_tube(SilhouetteData *sil, SpineBranch *branch, Mesh *me, float z_vec[3], float depth, int ss_steps, int w_steps, float w_fact)
+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[branch->tot_hull_points * 4], right[branch->tot_hull_points * 4];
+       float *left = NULL, *right = NULL;
        int totl = 0, totr = 0;
        int u_steps = 0;
-       const int v_start = me->totvert;
-       int square_ss_steps;
        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);
 
@@ -6313,13 +6873,11 @@ static void add_ss_tube(SilhouetteData *sil, SpineBranch *branch, Mesh *me, floa
                }
        }
 
-       totl = 0;
-       f_swap = false;
-       n_i = 0;
-
-       for (int i = 0; i < branch->tot_hull_points; i++) {
+       /* 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];
@@ -6328,6 +6886,7 @@ static void add_ss_tube(SilhouetteData *sil, SpineBranch *branch, Mesh *me, floa
                        }
                        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];
@@ -6344,39 +6903,49 @@ static void add_ss_tube(SilhouetteData *sil, SpineBranch *branch, Mesh *me, floa
                }
        }
 
-       int pos_i_t = 0;
-       float tmp_p[4];
-       for (int i = 0; i < totl; i++) {
-               copy_v4_v4(tmp_p, &left[(pos_i_t + cyclic_offset) % totl * 4]);
-               copy_v4_v4(&left[(pos_i_t + cyclic_offset) % totl * 4], &left[pos_i_t * 4]);
-       }
-
        if (totl < 1 && totr < 1) {
                return;
        }
-       /*printf("Creating tube with %i left points and %i right points length: %f, %f\n", totl, totr, left[totl * 4 - 1], right[totr * 4 - 1]);*/
 
-       u_steps = fmax(2.0f, fmax(left[totl * 4 - 1], right[totr * 4 - 1]) / (float)(2 * depth / ss_steps));
-       square_ss_steps = u_steps * ss_steps;
+       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) * 18,"edge startposition array");
+               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;
        }
-       fill_tube(me, left, right, totl, totr, u_steps, z_vec, ss_steps, w_steps, w_fact, branch->e_start_arr);
 
-       branch->flag |= BRANCH_VERT_GEN;
-       int verts_per_intersection = (2 * ss_steps + w_steps - 1);
-       branch->gen_verts = MEM_callocN(sizeof(int) * verts_per_intersection * 2, "Branch Bridges");
-       for (int v = 0; v < ss_steps; v ++) {
-               branch->gen_verts[verts_per_intersection + v] = v_start + v;
-               branch->gen_verts[verts_per_intersection + v + ss_steps + w_steps - 1] = v_start + square_ss_steps + v;
-               branch->gen_verts[v] = v_start + v + u_steps * ss_steps - ss_steps;
-               branch->gen_verts[ss_steps + w_steps - 1 + v] = v_start + square_ss_steps + v + u_steps * ss_steps - ss_steps;
-               if (v < w_steps - 1) {
-                       branch->gen_verts[verts_per_intersection + ss_steps + v] = v_start + square_ss_steps * 2 + v;
-                       branch->gen_verts[ss_steps + v] = v_start + square_ss_steps * 2 + v + u_steps * (w_steps - 1) - (w_steps - 1);
-               }
+       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.*/
@@ -6391,10 +6960,12 @@ static int r_branch_count(Spine *spine, SpineBranch *b)
        return r_forks;
 }
 
-/* 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)
+/* 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 = 0;
+       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) {
@@ -6408,22 +6979,65 @@ static void bridge_all_parts_rec(Mesh *me, Spine *spine, SpineBranch *active_bra
                                                b_fork_off ++;
                                        }
                                }
-                               for (int s = 0; s < 3; s++) {
-                                       if (s == 1 && (comp->totforks == 3 || active_branch->totforks == 3)) {
-                                               /* TODO: connect center */
-                                       } else {
-                                               bridge_loops(me,
-                                                                        active_branch->e_start_arr[a_fork_off * 9 + s * 3],
-                                                                        comp->e_start_arr[b_fork_off * 9 + s * 3],
-                                                                        active_branch->e_start_arr[a_fork_off * 9 + s * 3 + 1],
-                                                                        false,
-                                                                        active_branch->e_start_arr[a_fork_off * 9 + s * 3 + 2],
-                                                                        comp->e_start_arr[b_fork_off * 9 + s * 3 + 2],
-                                                                        !(s & 2));
-                                       }
+
+                               /* 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);
+                       bridge_all_parts_rec(me, spine, comp, active_branch, verts_per_loop, n_g_flip);
                }
                if (comp) {
                        a_fork_off ++;
@@ -6431,7 +7045,7 @@ static void bridge_all_parts_rec(Mesh *me, Spine *spine, SpineBranch *active_bra
        }
 }
 
-static void bridge_all_parts(Mesh *me, Spine *spine)
+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++) {
@@ -6445,137 +7059,792 @@ static void bridge_all_parts(Mesh *me, Spine *spine)
                /* No Branches in the spine. Should not happen. */
                return;
        }
-       bridge_all_parts_rec(me, spine, active_branch, NULL);
+       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 = 3;
-       int ss_steps = (1 << ss_level) + 2;
+       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);
 
-       //Generate spine
+       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 / (float)ss_steps);
+       mul_v3_fl(z_vec, depth);
+       mul_v3_v3fl(inv_z_vec, z_vec, -1.0f);
 
-       int w_steps = ss_steps / 2 + 2;
-       float w_fact = sil->smoothness;
+       int w_steps = v_steps / 2 + 2;
+       float smoothness = sil->smoothness;
 
-       /* Generate Tubes first then Caps and T-Intersections later. */
        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);
-                       if (r_forks == 2) {
-                               add_ss_tube(sil, a_branch, me, z_vec, depth, ss_steps, w_steps, w_fact);
-                               int t_off = 0;
-                               for (int sb = 0; sb < a_branch->totforks; sb++) {
-                                       SpineBranch *sub_branch = spine->branches[a_branch->terminal_points[sb * 2 + 1]];
-                                       if (sub_branch && a_branch->flag & BRANCH_VERT_GEN && r_branch_count(spine, sub_branch) == 1) {
-                                               add_ss_cap(sil, sub_branch, a_branch->gen_verts + t_off, me, z_vec, depth, ss_steps, w_steps - 1, w_fact);
-                                               t_off += ss_steps * 2 + (ss_steps / 2) - 1;
+                       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(sub_branch, spine, sil, 0xff00ff);
+                                       /*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, spine, sil, 0x00ff00);
+                                       /*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;
                        }
                }
        }
 
-       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);
+       bridge_all_parts(me, spine, v_steps * 2 + w_steps, n_ori);
 
-                       if (r_forks == 3) {
-                               add_ss_tinter(sil, a_branch, me, z_vec, ss_steps, w_steps - 1, w_fact);
-                               for (int sb = 0; sb < a_branch->totforks; sb++) {
-                                       SpineBranch *sub_branch = spine->branches[a_branch->terminal_points[sb * 2 + 1]];
-                                       if (sub_branch && r_branch_count(spine, sub_branch) == 1) {
-                                               add_ss_cap(sil, sub_branch, NULL, me, z_vec, depth, ss_steps, w_steps - 1, w_fact);
-#ifdef DEBUG_DRAW
-                                               //debug_branch(sub_branch, spine, sil, 0x0000ff);
-#endif
+       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;
                                        }
                                }
-#ifdef DEBUG_DRAW
-                               //debug_branch(a_branch, spine, sil, 0x0000ff);
-#endif
                        }
+               } 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));
                }
        }
+}
 
-       bridge_all_parts(me, spine);
+/* 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;
 
-       free_spine(spine);
+       int n_totvert, n_totedge, n_totloop, n_totpoly;
 
-       ED_mesh_update(me, C, 1, 1);
-}
+       CustomData vdata, edata, ldata, pdata;
 
-#if 0  /* UNUSED */
-static void debug_mesh(Mesh *me)
-{
-       printf("Logging Mesh:\n");
-       printf("Verts in mesh %i\n",me->totvert);
-       printf("Edges in mesh %i\n",me->totedge);
-       printf("Loops in mesh %i\n",me->totloop);
-       printf("Polys in mesh %i\n",me->totpoly);
+       /* 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");
 
-       printf("\nVert log:\n");
+       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++) {
-               printf("\nVert %i (%f,%f,%f)", i, me->mvert[i].co[0], me->mvert[i].co[1], me->mvert[i].co[2]);
+               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 ++;
+               }
        }
 
-       printf("\nEdge log:\n");
-       for (int i = 0; i < me->totedge;i++) {
-               printf("Edge %i, v1v2(%u,%u)\n",i, me->medge[i].v1, me->medge[i].v2);
+       /* 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 ++;
+               }
        }
 
-       printf("\nLoop Log:\n");
+       /* 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++) {
-               printf("Loop %i, v(%u), e(%u)\n", i, me->mloop[i].v, me->mloop[i].e);
+               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 ++;
+               }
        }
 
-       printf("\nPoly log:\n");
+       /* 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++) {
-               printf("Poly %i, start(%i), totloop(%i)\n", i, me->mpoly[i].loopstart, me->mpoly[i].totloop);
+               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);
 }
-#endif
 
-/* Adds additional points to the stroke if start and end are far apart. */
-static void stroke_smooth_cap(SilhouetteStroke *stroke, float max_dist)
+static void calc_ring_bbs(SilhouetteData *sil, Mesh *me)
 {
-       float v1[3], v2[3], dv[3];
-       float length = 0;
-       copy_v3_v3(v1, &stroke->points[0]);
-       copy_v3_v3(v2, &stroke->points[stroke->totvert * 3 - 3]);
+       int len, edge;
+       BB *curr;
+       sil->fillet_ring_bbs = MEM_callocN(sizeof(BB) * sil->num_rings, "ring bb mem");
 
-       sub_v3_v3v3(dv,v1,v2);
-       length = len_v3(dv);
+       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);
+               }
+       }
+}
 
-       if (length > max_dist) {
-               int steps = floorf(length / (float)max_dist);
-               mul_v3_fl(dv, 1.0f / (float)steps);
-               for (int i = 1; i < steps; i++) {
-                       mul_v3_v3fl(v1, dv, i);
-                       add_v3_v3(v1, v2);
-                       silhouette_stroke_add_3Dpoint(stroke, v1);
+/* 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)
@@ -6584,27 +7853,52 @@ static void sculpt_silhouette_calc_mesh(bContext *C, wmOperator *op)
        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(stroke, 0.3f);
+
+       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 
+       /* 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;
 
-#ifdef DEBUG_DRAW
        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);
+               /*bl_debug_draw_edge_add(v1,v2);*/
        }
 #endif
        
@@ -6624,13 +7918,30 @@ static void sculpt_silhouette_clean_draw(bContext *C, wmOperator *op)
 
 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, op, true);
+               sil = silhouette_data_new(C);
+               silhouette_set_ref_plane(sil);
                op->customdata = sil;
        }
-       sculpt_silhouette_calc_mesh(C, op);
-       sculpt_silhouette_stroke_done(C, op);
+
+       /*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;
 }
@@ -6638,24 +7949,19 @@ static int sculpt_silhouette_exec(bContext *C, wmOperator *op)
 static int sculpt_silhouette_modal(bContext *C, wmOperator *op, const wmEvent *event)
 {
        float mouse[2];
-       float z_vec[3];
        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) {
-                       RNA_float_set_array(op->ptr, "points", sil->current_stroke->points);
-                       RNA_int_set(op->ptr,"totvert", sil->current_stroke->totvert);
-                       ED_view3d_global_to_vector(sil->ar->regiondata, (float[3]){0.0f,0.0f,0.0f}, z_vec);
-                       RNA_float_set_array(op->ptr, "z_vec", z_vec);
-                       RNA_float_set_array(op->ptr, "anchor", sil->anchor);
+                       silhouette_set_ref_plane(sil);
                        return sculpt_silhouette_exec(C, op);
                }
                return OPERATOR_FINISHED;
        } else {
                if (sil->state == SIL_DRAWING) {
-                       sculpt_silhouette_stroke_update(C,mouse,op->customdata);
+                       sculpt_silhouette_stroke_update(mouse, op->customdata);
                }
                return OPERATOR_RUNNING_MODAL;
        }
@@ -6700,7 +8006,7 @@ static int sculpt_silhouette_invoke(bContext *C, wmOperator *op, const wmEvent *
 
        SilhouetteData *sil_data;
 
-       sil_data = silhouette_data_new(C, op, false);
+       sil_data = silhouette_data_new(C);
 
        op->customdata = sil_data;
 
@@ -6724,8 +8030,6 @@ static int sculpt_silhouette_poll(bContext *UNUSED(C))
 
 static void SCULPT_OT_silhouette_draw(wmOperatorType *ot)
 {
-       PropertyRNA *prop;
-
        /* identifiers */
        ot->name = "Draw Silhouette";
        ot->idname = "SCULPT_OT_silhouette_draw";
@@ -6739,32 +8043,7 @@ static void SCULPT_OT_silhouette_draw(wmOperatorType *ot)
        ot->cancel = sculpt_silhouette_stroke_done;
 
        /* flags */
-       ot->flag = OPTYPE_BLOCKING | OPTYPE_REGISTER | OPTYPE_UNDO;
-
-       /* properties */
-       prop = RNA_def_int(ot->srna, "ss_level", 3, 1, 16, "Subsurface Level", "", 1, 8);
-       RNA_def_property_flag(prop, PROP_SKIP_SAVE);
-
-       prop = RNA_def_float(ot->srna, "depth", 1.5f, 0.0f, 100000000.0f,
-                                                "Depth", "Shape depth", 0.0f, 100.0f);
-       RNA_def_property_flag(prop, PROP_SKIP_SAVE);
-
-       prop = RNA_def_float(ot->srna, "smoothness", 0.75f, 0.0f, 1.0f,
-                                                "Smoothness", "Smoothness factor", 0.0f, 1.0f);
-       RNA_def_property_flag(prop, PROP_SKIP_SAVE);
-
-       /* TODO: RNA points max size 30 */
-       prop = RNA_def_float_vector(ot->srna, "points", 3 * 1024, NULL, -FLT_MAX, FLT_MAX, "3D Points", "Stroke Mouse locations", -FLT_MAX, FLT_MAX);
-       RNA_def_property_flag(prop, PROP_HIDDEN);
-
-       prop = RNA_def_int(ot->srna, "totvert", 0, 0, INT_MAX, "stroke point count", "", 0, INT_MAX);
-       RNA_def_property_flag(prop, PROP_HIDDEN);
-
-       prop = RNA_def_float_vector(ot->srna, "z_vec", 3, NULL, -FLT_MAX, FLT_MAX, "Z Vector", "Silhouette orientation", -FLT_MAX, FLT_MAX);
-       RNA_def_property_flag(prop, PROP_HIDDEN);
-
-       prop = RNA_def_float_vector(ot->srna, "anchor", 3, NULL, -FLT_MAX, FLT_MAX, "Anchor", "Silhouette position", -FLT_MAX, FLT_MAX);
-       RNA_def_property_flag(prop, PROP_HIDDEN);
+       ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
 }
 /* end Silhouette */