SnapSystem: begun implementing better state and event handling for the snap system

compilation.

extern/libmv/third_party/ceres/internal/ceres/schur_ordering.cc

@@ -39,16 +39,16 @@
 #include "ceres/residual_block.h"
 #include "ceres/internal/scoped_ptr.h"
 
-CERES_HASH_NAMESPACE_START
-
+//CERES_HASH_NAMESPACE_START
+namespace std {
 // Allow us to hash pointers as if they were int's
 template<> struct hash< ::ceres::internal::ParameterBlock*> {
   size_t operator()(::ceres::internal::ParameterBlock* x) const {
     return reinterpret_cast<size_t>(x);
   }
 };
-
-CERES_HASH_NAMESPACE_END
+}
+//CERES_HASH_NAMESPACE_END
 
 namespace ceres {
 namespace internal {

release/scripts/startup/bl_ui/properties_data_modifier.py

@@ -315,6 +315,24 @@ class DATA_PT_modifiers(ModifierButtonsPanel, Panel):
             row = layout.row()
             row.operator("object.hook_select", text="Select")
             row.operator("object.hook_assign", text="Assign")
+    def LAPLACIANSMOOTH(self, layout, ob, md):
+        #split = layout.split(percentage=0.25)
+
+        #col = split.column()
+
+        #col = split.column()
+        layout.prop(md, "iterations")
+        layout.prop(md, "lamb")
+        layout.prop(md, "lambdaborder")
+        row = layout.row()
+        row.label(text="Axis: ")
+        row.prop(md, "use_x")
+        row.prop(md, "use_y")
+        row.prop(md, "use_z")
+        row = layout.row()
+        row.prop(md, "volume_preservation")
+        layout.label(text="Vertex Group:")
+        layout.prop_search(md, "vertex_group", ob, "vertex_groups", text="")
 
     def LATTICE(self, layout, ob, md):
         split = layout.split()

release/scripts/startup/bl_ui/space_view3d.py

@@ -1699,6 +1699,7 @@ class VIEW3D_MT_edit_mesh_specials(Menu):
         layout.operator("mesh.select_all", text="Select Inverse").action = 'INVERT'
         layout.operator("mesh.flip_normals")
         layout.operator("mesh.vertices_smooth", text="Smooth")
+        layout.operator("mesh.vertices_smooth_laplacian", text="Laplacian Smooth")
         layout.operator("mesh.inset")
         layout.operator("mesh.bevel", text="Bevel")
         layout.operator("mesh.bridge_edge_loops")

source/blender/blenkernel/BKE_bvhutils.h

@@ -109,6 +109,7 @@ void free_bvhtree_from_mesh(struct BVHTreeFromMesh *data);
  */
 float bvhtree_ray_tri_intersection(const BVHTreeRay *ray, const float m_dist, const float v0[3], const float v1[3], const float v2[3]);
 float nearest_point_in_tri_surface(const float v0[3], const float v1[3], const float v2[3], const float p[3], int *v, int *e, float nearest[3]);
+float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, const float m_dist, const float v0[3], const float v1[3], const float v2[3]);
 
 /*
  * BVHCache

source/blender/blenkernel/BKE_snap.h

@@ -0,0 +1,178 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * Contributor(s): Luke Frisken 2012
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+
+#ifndef __BKE_SNAP_H__
+#define __BKE_SNAP_H__
+
+#include "DNA_scene_types.h"
+#include "DNA_object_types.h"
+#include "DNA_view3d_types.h"
+#include "DNA_screen_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_mesh_types.h"
+#include "BKE_context.h"
+#include "WM_types.h"
+
+typedef struct{
+	//perhaps more info later I might guess
+	float location[3];
+	float normal[3];
+	int r_dist;
+	float r_depth;
+}SnapPoint;
+
+typedef struct SnapSystem SnapSystem;
+
+typedef enum{
+	SNAPSYSTEM_MODE_SELECT_ALL = 0,
+	SNAPSYSTEM_MODE_SELECT_NOT_SELECTED = 1,
+	SNAPSYSTEM_MODE_SELECT_NOT_OBEDIT = 2
+}SnapSystem_mode_select;
+
+typedef enum{
+	SNAPSYSTEM_STATE_SNAPPING,
+	SNAPSYSTEM_STATE_INIT_SNAP,
+	SNAPSYSTEM_STATE_IDLE
+}SnapSystem_state;
+
+/*This function creates a new instance of the snapping system*/
+/*callback_data: is the data to be used by the code using the snapping system (e.g transform)
+  when the callbacks are called by snapsystem.
+
+  update_callback: gets called when the snapping system has not yet recieved confirmation from the user
+  of the chosen snap point, but there is a temporary snap_point available for providing realtime feedback
+  to the user of their actions when moving the mouse. In the case of transform, this means the geometry
+  moves towards the snap point without permenantly being applied.
+
+  finish_callback: same as update_callback, but the snap_point is to be applied this time.
+
+  cancel_callback: tells the code using the snapsystem that the current snapsystem snap search has been cancelled.
+
+  nofeedback_callback: will probably be removed/changed/refactored soon.
+  */
+/*The return value pointer to SnapSystem, when finished requires the code using the system to call SnapSystem_free
+  on it to free allocated data within the system*/
+SnapSystem* SnapSystem_create( Scene* scene, View3D* v3d, ARegion* ar, Object* obedit, bContext* C,
+								void* callback_data,
+								void (*update_callback)(SnapSystem *ssystem, void* callback_data, SnapPoint sp),
+							   void (*nofeedback_callback)(SnapSystem *ssystem, void* callback_data),
+								void (*finish_callback)(SnapSystem *ssystem, void* callback_data, SnapPoint sp),
+							   void (*cancel_callback)(SnapSystem *ssystem, void* callback_data));
+
+void SnapSystem_test_run(SnapSystem* ssystem);
+void SnapSystem_reset_ob_list(SnapSystem* ssystem); //reset the list of objects the snapsystem uses for snapping
+void SnapSystem_evaluate_stack(SnapSystem* ssystem);
+void SnapSystem_find_snap(SnapSystem* ssystem);
+void SnapSystem_pick_snap(SnapSystem* ssystem);
+void SnapSystem_add_ob(SnapSystem* ssystem, Object* ob);
+Object* SnapSystem_get_ob(SnapSystem* ssystem, int index);
+void SnapSystem_set_mode_select(SnapSystem* ssystem, SnapSystem_mode_select mode_select);
+int SnapSystem_get_retval(SnapSystem* ssystem);
+bContext* SnapSystem_get_C(SnapSystem* ssystem);
+SnapSystem_state SnapSystem_get_state(SnapSystem* ssystem);
+
+/*This function handles the events of the snapsystem, and returns a 1 if the event has been
+  handled (and consumed) by the snapsystem*/
+int SnapSystem_Event(SnapSystem* ssystem, wmEvent* event);
+/*This function draws all the currently active snaps within the snapsystem*/
+void SnapSystem_draw(SnapSystem* ssystem);
+
+void SnapSystem_reset_snappoint(SnapSystem* ssystem);
+void SnapSystem_clear_pick(SnapSystem* ssystem);
+void SnapSystem_reset(SnapSystem* ssystem);
+void SnapSystem_free(SnapSystem* ssystem);
+
+/*default callbacks*/
+void SnapSystem_default_object_iterator(SnapSystem* ssystem, void* callback_data);
+void SnapSystem_default_object_handler(SnapSystem* ssystem, void* callback_data, Object* ob);
+
+typedef struct Snap Snap;
+
+typedef enum{
+	SNAP_TYPE_MESH,
+	SNAP_TYPE_CURVE,
+	SNAP_TYPE_BONE,
+	SNAP_TYPE_AXIS
+}Snap_type;
+
+typedef enum{
+	SNAPMESH_DATA_TYPE_DerivedMesh,
+	SNAPMESH_DATA_TYPE_BMEditMesh
+}SnapMesh_data_type;
+
+typedef enum{
+	SNAPMESH_TYPE_VERTEX,
+	SNAPMESH_TYPE_FACE,
+	SNAPMESH_TYPE_PLANAR,
+	SNAPMESH_TYPE_EDGE,
+	SNAPMESH_TYPE_EDGE_MIDPOINT,
+	SNAPMESH_TYPE_EDGE_PARALLEL
+}SnapMesh_type;
+
+typedef enum{
+	SNAPMESH_DAT_vert,
+	SNAPMESH_DAT_edge,
+	SNAPMESH_DAT_face
+}SnapMesh_data_array_type;
+
+/*This struct provides a storage of face + verts in one place for using
+  as return geometry in face picking, and gets used for planar snapping*/
+typedef struct{
+	MFace *face;
+	MVert *verts;
+	int nverts;
+	float no[3];
+}MeshData_pickface;
+
+/*Return data type for SnapMesh class*/
+typedef enum{
+	SNAPMESH_RET_DAT_vert_index,
+	SNAPMESH_RET_DAT_edge_index,
+	SNAPMESH_RET_DAT_pface //see MeshData_pickface struct
+}SnapMesh_ret_data_type;
+
+/*Flags for use in VertCheck or EdgeCheck or FaceCheck MeshData functions.*/
+typedef enum{
+	MESHDATA_CHECK_NONE = 0,
+	MESHDATA_CHECK_HIDDEN = 1,
+	MESHDATA_CHECK_SELECTED = 2
+}MeshData_check;
+
+void Snap_run(Snap* s);
+void Snap_free(Snap* s);
+
+int Snap_getretval(Snap* s);
+SnapPoint* Snap_getSnapPoint(Snap* s);
+void Snap_setClosestPoint(Snap* s, SnapPoint* sp);
+void Snap_setpick(Snap* s, Snap* pick);
+
+Snap* SnapMesh_create(	void* mesh_data,
+						SnapMesh_data_type data_type,
+						int free_mesh_data,
+						MeshData_check check,
+						SnapMesh_type sm_type,
+						Scene *scene, Object *ob, View3D *v3d, ARegion *ar, bContext *C, int mval[2]);
+
+void SnapMesh_draw(Snap *sm);
+
+#endif // __BKE_SNAP_H__

source/blender/blenkernel/CMakeLists.txt

@@ -29,6 +29,7 @@ set(INC
 	../blenfont
 	../blenlib
 	../blenloader
+	../editors/include
 	../gpu
 	../ikplugin
 	../imbuf
@@ -44,6 +45,8 @@ set(INC
 	../../../intern/mikktspace
 	../../../intern/opennl/extern
 	../../../intern/raskter
+	../editors/transform
+	
 
 	# XXX - BAD LEVEL CALL WM_api.h
 	../windowmanager
@@ -133,6 +136,7 @@ set(SRC
 	intern/shrinkwrap.c
 	intern/sketch.c
 	intern/smoke.c
+	intern/snap.c
 	intern/softbody.c
 	intern/sound.c
 	intern/speaker.c
@@ -217,6 +221,7 @@ set(SRC
 	BKE_shrinkwrap.h
 	BKE_sketch.h
 	BKE_smoke.h
+	BKE_snap.h
 	BKE_softbody.h
 	BKE_sound.h
 	BKE_speaker.h
@@ -255,7 +260,6 @@ endif()
 if(WITH_MOD_CLOTH_ELTOPO)
 	list(APPEND INC
 		../../../extern/eltopo
-		../../../extern/eltopo/eltopo3d
 	)
 	add_definitions(-DWITH_ELTOPO)
 endif()

source/blender/blenkernel/SConscript

@@ -7,7 +7,7 @@ sources = env.Glob('intern/*.c')
 incs = '. #/intern/guardedalloc #/intern/memutil'
 incs += ' ../blenlib ../blenfont ../makesdna ../windowmanager'
 incs += ' ../render/extern/include #/intern/decimation/extern ../makesrna'
-incs += ' ../imbuf ../ikplugin ../avi #/intern/elbeem/extern ../nodes ../modifiers'
+incs += ' ../imbuf ../ikplugin ../avi #/intern/elbeem/extern ../nodes ../modifiers ../editors/include'
 incs += ' #/intern/iksolver/extern ../blenloader'
 incs += ' #/extern/bullet2/src'
 incs += ' #/intern/opennl/extern #/intern/bsp/extern'
@@ -18,6 +18,7 @@ incs += ' #/intern/mikktspace'
 incs += ' #/intern/audaspace/intern'
 incs += ' #/intern/ffmpeg'
 incs += ' #/intern/raskter'
+incs += ' ../editors/include ../editors/transform'
 
 incs += ' ' + env['BF_OPENGL_INC']
 incs += ' ' + env['BF_ZLIB_INC']

source/blender/blenkernel/intern/bvhutils.c

@@ -57,7 +57,7 @@ float bvhtree_ray_tri_intersection(const BVHTreeRay *ray, const float UNUSED(m_d
 	return FLT_MAX;
 }
 
-static float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, const float m_dist, const float v0[3], const float v1[3], const float v2[3])
+float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, const float m_dist, const float v0[3], const float v1[3], const float v2[3])
 {
 	
 	float idist;

source/blender/bmesh/CMakeLists.txt

@@ -30,6 +30,7 @@ set(INC
 	../blenlib
 	../makesdna
 	../../../intern/guardedalloc
+	../../../intern/opennl/extern
 )
 
 set(INC_SYS
@@ -52,6 +53,7 @@ set(SRC
 	operators/bmo_mirror.c
 	operators/bmo_primitive.c
 	operators/bmo_removedoubles.c
+	operators/bmo_smooth_laplacian.c
 	operators/bmo_subdivide.c
 	operators/bmo_subdivide.h
 	operators/bmo_triangulate.c

source/blender/bmesh/SConscript

@@ -13,6 +13,7 @@ incs = [
 	'../makesdna',
 	'../blenkernel',
 	'#/intern/guardedalloc',
+	'#/intern/opennl/extern',
 	]
 
 defs = []

source/blender/bmesh/intern/bmesh_inline.h

@@ -29,6 +29,9 @@
 #ifndef __BMESH_INLINE_H__
 #define __BMESH_INLINE_H__
 
+#include "bmesh_class.h"
+#include "BLI_utildefines.h"
+
 /* stuff for dealing with header flags */
 #define BM_elem_flag_test(     ele, hflag)      _bm_elem_flag_test     (&(ele)->head, hflag)
 #define BM_elem_flag_test_bool(ele, hflag)      _bm_elem_flag_test_bool(&(ele)->head, hflag)

source/blender/bmesh/intern/bmesh_opdefines.c

@@ -113,6 +113,26 @@ static BMOpDefine bmo_smooth_vert_def = {
 	0
 };
 
+/*
+ * Vertext Smooth Laplacian 
+ * Smooths vertices by using Laplacian smoothing propose by.
+ * Desbrun, et al. Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow
+*/
+static BMOpDefine bmo_smooth_laplacian_vert_def = {
+	"smooth_laplacian_vert",
+	{{BMO_OP_SLOT_ELEMENT_BUF, "verts"}, //input vertices
+	 {BMO_OP_SLOT_FLT, "lambda"}, //lambda param
+	 {BMO_OP_SLOT_FLT, "lambda_border"}, //lambda param in border
+	 {BMO_OP_SLOT_BOOL, "use_x"}, //Smooth object along X axis
+	 {BMO_OP_SLOT_BOOL, "use_y"}, //Smooth object along Y axis
+	 {BMO_OP_SLOT_BOOL, "use_z"}, //Smooth object along Z axis
+	 {BMO_OP_SLOT_BOOL, "volume_preservation"}, //Apply volume preservation after smooth
+	{0} /* null-terminating sentinel */,
+	},
+	bmo_smooth_laplacian_vert_exec,
+	0
+};
+
 /*
  * Right-Hand Faces
  *
@@ -420,8 +440,15 @@ static BMOpDefine bmo_contextual_create_def = {
  */
 static BMOpDefine bmo_bridge_loops_def = {
 	"bridge_loops",
-	{{BMO_OP_SLOT_ELEMENT_BUF, "edges"}, /* input edge */
+    //{{BMO_OP_SLOT_ELEMENT_BUF, "edges"}, /* input edge */
+    {{BMO_OP_SLOT_ELEMENT_BUF, "edgefacein"}, /* input faces and edge */
 	 {BMO_OP_SLOT_ELEMENT_BUF, "faceout"}, /* new face */
+
+     {BMO_OP_SLOT_INT,  "segmentation"}, /* input segmentation param */
+     {BMO_OP_SLOT_INT,  "interpolation"}, /* input interpolation param */
+     {BMO_OP_SLOT_FLT,  "strenght"}, /* input strenght param */
+	 {BMO_OP_SLOT_BOOL, "n_gon"},	 /* flag N-gon optimization */
+
 	 {0, /* null-terminating sentinel */}},
 	bmo_bridge_loops_exec,
 	0,
@@ -1043,6 +1070,9 @@ static BMOpDefine bmo_bevel_def = {
 	                                  * modifier uses this. We could do this as another float setting */
 	 {BMO_OP_SLOT_INT, "lengthlayer"}, /* which PROP_FLT layer to us */
 	 {BMO_OP_SLOT_FLT, "percent"}, /* percentage to expand beveled edge */
+
+	 {BMO_OP_SLOT_FLT, "amount"},
+	 {BMO_OP_SLOT_INT, "segmentation"},
 	 {0} /* null-terminating sentinel */},
 	bmo_bevel_exec,
 	BMO_OP_FLAG_UNTAN_MULTIRES
@@ -1246,6 +1276,7 @@ BMOpDefine *opdefines[] = {
 	&bmo_translate_def,
 	&bmo_triangle_fill_def,
 	&bmo_triangulate_def,
+	&bmo_smooth_laplacian_vert_def,
 	&bmo_weld_verts_def,
 	&bmo_wireframe_def,
 

source/blender/bmesh/intern/bmesh_operator_api.h

@@ -32,6 +32,7 @@ extern "C" {
 #endif
 
 #include "BLI_ghash.h"
+#include "BLI_utildefines.h"
 
 #include <stdarg.h>
 

source/blender/bmesh/intern/bmesh_operators_private.h

@@ -60,10 +60,13 @@ void bmo_dissolve_faces_exec(BMesh *bm, BMOperator *op);
 void bmo_dissolve_limit_exec(BMesh *bm, BMOperator *op);
 void bmo_dissolve_verts_exec(BMesh *bm, BMOperator *op);
 void bmo_duplicate_exec(BMesh *bm, BMOperator *op);
+void bmo_edgebisect_exec(BMesh *bm, BMOperator *op);
 void bmo_edgenet_fill_exec(BMesh *bm, BMOperator *op);
 void bmo_edgenet_prepare(BMesh *bm, BMOperator *op);
+void bmo_edgerotate_exec(BMesh *bm, BMOperator *op);
 void bmo_extrude_discrete_faces_exec(BMesh *bm, BMOperator *op);
 void bmo_extrude_edge_only_exec(BMesh *bm, BMOperator *op);
+void bmo_extrude_face_indiv_exec(BMesh *bm, BMOperator *op);
 void bmo_extrude_face_region_exec(BMesh *bm, BMOperator *op);
 void bmo_extrude_vert_indiv_exec(BMesh *bm, BMOperator *op);
 void bmo_find_doubles_exec(BMesh *bm, BMOperator *op);
@@ -91,6 +94,7 @@ void bmo_similar_faces_exec(BMesh *bm, BMOperator *op);
 void bmo_similar_verts_exec(BMesh *bm, BMOperator *op);
 void bmo_slide_vert_exec(BMesh *bm, BMOperator *op);
 void bmo_smooth_vert_exec(BMesh *bm, BMOperator *op);
+void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op);
 void bmo_solidify_face_region_exec(BMesh *bm, BMOperator *op);
 void bmo_spin_exec(BMesh *bm, BMOperator *op);
 void bmo_split_edges_exec(BMesh *bm, BMOperator *op);

source/blender/bmesh/operators/bmo_connect.c

@@ -26,6 +26,7 @@
 
 #include "MEM_guardedalloc.h"
 
+#include "BLI_listbase.h"
 #include "BLI_math.h"
 #include "BLI_array.h"
 #include "BLI_utildefines.h"
@@ -34,12 +35,38 @@
 
 #include "intern/bmesh_operators_private.h" /* own include */
 
+
 #define VERT_INPUT	1
 #define EDGE_OUT	1
 #define FACE_NEW	2
 #define EDGE_MARK	4
 #define EDGE_DONE	8
 
+#define FACE_MARK       1
+#define EDGE_CONNECTED  7
+#define EDGE_NON_CONNECTED  9
+
+#define LINEAR_INTER  2
+#define CUBIC_INTER   4
+
+#define BRIDGE_EPSILON  1e-6
+
+typedef struct VertexItem{
+    struct VertexItem *next, *prev;
+    BMVert *v;
+}VertexItem;
+
+
+typedef struct BridgeParams{
+    ListBase newVertices; // list of new vertex
+    int inter;  // interpolation
+    int seg; // segmentation param
+    float strenght;
+    float centrod1[3]; // coordinate of centrod input loops
+    float centrod2[3];
+} BridgeParams;
+
+
 void bmo_connect_verts_exec(BMesh *bm, BMOperator *op)
 {
 	BMIter iter, liter;
@@ -178,7 +205,7 @@ static int clamp_index(const int x, const int len)
 			BLI_array_append(arr2, arr_tmp[i]);                               \
 		}                                                                     \
 		BLI_array_free(arr_tmp);                                              \
-	} (void)0
+	}
 
 /* get the 2 loops matching 2 verts.
  * first attempt to get the face corners that use the edge defined by v1 & v2,
@@ -210,7 +237,715 @@ static void bm_vert_loop_pair(BMesh *bm, BMVert *v1, BMVert *v2, BMLoop **l1, BM
 	*l2 = BM_iter_at_index(bm, BM_LOOPS_OF_VERT, v2, 0);
 }
 
-void bmo_bridge_loops_exec(BMesh *bm, BMOperator *op)
+//TODO: rename function
+// return count entry element into list
+int bmo_bridge_array_entry_elem(BMEdge **list, int size,  BMEdge *elem)
+{
+    int i = 0, count = 0;
+
+    for (i = 0; i < size; i++)
+    {
+        if (list[i] == elem)
+            count ++;
+    }
+    return count;
+}
+
+/*
+* return coordinate of middle edge
+*/
+void  get_middle_edge(BMEdge* e, float v[3]){
+    BMVert *v1, *v2;
+    v1 = e->v1;
+    v2 = BM_edge_other_vert(e, v1);
+    mid_v3_v3v3(v, v1->co, v2->co);
+}
+
+/*
+*   this function find centroid of group vertex
+*   return coordinate of centroid
+*/
+void get_centroid(BMVert **vv, int vert_count, float center[3]){
+    int i;
+    center[0] = 0; center[1] = 0; center[2] = 0;
+    for (i = 0; i < vert_count; i++){
+        add_v3_v3(center, vv[i]->co);
+    }
+    mul_v3_fl(center, 1.0f/vert_count);
+}
+
+
+BMVert* bridge_check_existing_vertex(BridgeParams *bp, float co[3])
+{
+    float epsilon = 1e-6;
+    BMVert *vert = NULL;
+    VertexItem *item;
+    for (item = bp->newVertices.first; item; item = item->next)
+    {
+        if (compare_v3v3(item->v->co, co, epsilon))
+            vert = item->v;
+    }
+    return vert;
+}
+/*
+* this function return vertex of linear interpolation between v1 and v2
+*/
+BMVert* get_linear_seg(BMesh *bm, BridgeParams* bp, int n, BMVert *v1, BMVert *v2)
+{
+    float co[3];
+    BMVert *v;
+    VertexItem *item;
+
+    co[0] = v1->co[0] + (v2->co[0]-v1->co[0])*n / bp->seg;
+    co[1] = v1->co[1] + (v2->co[1]-v1->co[1])*n / bp->seg;
+    co[2] = v1->co[2] + (v2->co[2]-v1->co[2])*n / bp->seg;
+
+    v = bridge_check_existing_vertex(bp, co);
+    if (!v)
+    {
+        item = (VertexItem*)MEM_callocN(sizeof(VertexItem), "VertexItem");
+        item->v = BM_vert_create(bm, co, NULL);
+        BLI_addtail(&bp->newVertices, item);
+        v = item->v;
+    }
+    return v;
+}
+
+/*
+* this function return vertex of cubic interpolation between v1 and v2
+* Parametric cubic spline in Hermite form
+* 0 <= t <= 1
+*/
+BMVert* get_cubic_seg(BMesh *bm, BridgeParams *bp, int n, BMVert *v1, BMVert *v2){
+    BMVert *v = NULL;
+    VertexItem *item;
+    float t, co[3], r1[3], r2[3]/*, mn[3]*/;
+    t = (float) (n) / bp->seg;
+
+    //sub_v3_v3v3(r1,  bp->centrod2 , v1->co);
+    //sub_v3_v3v3(r2,  v2->co, bp->centrod1);
+
+    sub_v3_v3v3(r1,  bp->centrod1 , v1->co);
+    sub_v3_v3v3(r2,  v2->co, bp->centrod2);
+    //sub_v3_v3v3(mn, v1->co, v2->co);
+
+    //mul_v3_fl(r1, bp->strenght);
+    //mul_v3_fl(r2, bp->strenght);
+    //add_v3_v3(r1, mn);
+    //add_v3_v3(r2, mn);
+
+    mul_v3_fl(r1, 2*asin(bp->strenght));
+    mul_v3_fl(r2, 2*asin(bp->strenght));
+
+    co[0] = v1->co[0] * (2 * t * t * t - 3 * t * t + 1) +
+            v2->co[0] * (-2 * t * t * t + 3 * t * t ) +
+            r1[0] * (t * t * t - 2 * t * t + t) +
+            r2[0] * (t * t * t - t * t);
+    co[1] = v1->co[1] * (2 * t * t * t - 3 * t * t + 1) +
+            v2->co[1] * (-2 * t * t * t + 3 * t * t ) +
+            r1[1] * (t * t * t - 2 * t * t + t) +
+            r2[1] * (t * t * t - t * t);
+    co[2] = v1->co[2] * (2 * t * t * t - 3 * t * t + 1) +
+            v2->co[2] * (-2 * t * t * t + 3 * t * t ) +
+            r1[2] * (t * t * t - 2 * t * t + t) +
+            r2[2] * (t * t * t - t * t);
+
+    v = bridge_check_existing_vertex(bp, co);
+    if (!v)
+    {
+        item = (VertexItem*)MEM_callocN(sizeof(VertexItem), "VertexItem");
+        item->v = BM_vert_create(bm, co, NULL);
+        BLI_addtail(&bp->newVertices, item);
+        v = item->v;
+    }
+    return v;
+}
+
+
+float get_cos_v3v3(float a[3], float b[3]){
+    float cos_ab = (a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) /
+            (sqrt(a[0]*a[0]+ a[1]*a[1]+ a[2]*a[2]) * sqrt(b[0]*b[0]+ b[1]*b[1]+ b[2]*b[2]));
+    return cos_ab;
+}
+/*
+  create polygons between two edge
+  */
+BMFace* bmo_edge_face_connect(BMesh *bm, BridgeParams* bp, BMEdge *e1, BMEdge *e2, BMFace *f_example)
+{
+    int i;
+    BMVert *vv[4];
+    BMVert *v1, *v2, *v3, *v4;
+    BMVert *vi1, *vi2, *vi3, *vi4; // input vertex
+	BMFace *f = NULL;
+    float vect1[3], vect2[3], vect3[3], vect4[3], normalA[3], normalB[3], normalVector[3];
+    float vp1[3], vp2[3], vp3[3], vp4[3];
+
+    vi1 = e1->v1;
+    vi2 = e2->v1;
+    vi3 = BM_edge_other_vert(e2,vi2);
+    vi4 = BM_edge_other_vert(e1,vi1);
+// calculate vector betwen input vertx
+    sub_v3_v3v3(vect1, vi1->co, vi2->co);  // v1-v2
+    sub_v3_v3v3(vect2, vi2->co, vi3->co);  // v2-v3
+    sub_v3_v3v3(vect3, vi3->co, vi4->co);  // v3-v4
+    sub_v3_v3v3(vect4, vi4->co, vi1->co);  // v4-v1
+
+// cros product
+    cross_v3_v3v3(vp1, vect1, vect2);
+    cross_v3_v3v3(vp2, vect2, vect3);
+    cross_v3_v3v3(vp3, vect3, vect4);
+    cross_v3_v3v3(vp4, vect4, vect1);
+// calculate normal vector
+   mid_v3_v3v3(normalA, bp->centrod1, bp->centrod2);
+   vv[0] = vi1;
+   vv[1] = vi2;
+   vv[2] = vi3;
+   vv[3] = vi4;
+   get_centroid(vv, 4, normalB);
+   sub_v3_v3v3(normalVector, normalB, normalA);
+
+    // check direction cros produc result
+    if (!(get_cos_v3v3(vp1, vp2)>0 &&
+        get_cos_v3v3(vp2, vp3)>0 &&
+        get_cos_v3v3(vp3, vp4)>0 &&
+        get_cos_v3v3(vp4, vp1)>0))
+    {
+        vi1 = e1->v1;
+        vi2 = BM_edge_other_vert(e2, e2->v1);
+        vi3 = BM_edge_other_vert(e2,vi2);
+        vi4 = BM_edge_other_vert(e1,vi1);
+    }
+    // vi1  - v1 - v2 ... vi2
+    // |      |    |       |
+    // vi4 .. v4 - v3 ... vi3
+   v1 = vi1;
+   v2 = vi2;
+   v3 = vi3;
+   v4 = vi4;
+    if (bp->seg > 1){
+        for (i = 1; i < bp->seg; i++){
+            if (bp->inter == LINEAR_INTER)
+            {
+                v2 = get_linear_seg (bm, bp, i, vi1, vi2);
+                v3 = get_linear_seg(bm, bp,  i, vi4, vi3);
+            }
+            if (bp->inter == CUBIC_INTER)
+            {
+                v2 = get_cubic_seg (bm, bp, i, vi1, vi2);
+                v3 = get_cubic_seg(bm, bp,  i, vi4, vi3);
+            }
+           f =  BM_face_create_quad_tri(bm, v1, v2, v3, v4, f_example, TRUE);
+           BM_face_normal_update(f);
+           //check normal BM_face_normal_flip(bm, f);
+           if (get_cos_v3v3(f->no, normalVector)<0)
+               BM_face_normal_flip(bm, f);
+           v1 = v2;
+           v4 = v3;
+        }
+       f =  BM_face_create_quad_tri(bm, v1, vi2, vi3, v4, f_example, TRUE);
+       BM_face_normal_update(f);
+       if (get_cos_v3v3(f->no, normalVector)<0)
+           BM_face_normal_flip(bm, f);
+    }
+    else {
+        f = BM_face_create_quad_tri(bm, v1, v2, v3, v4, f_example, TRUE);
+        BM_face_normal_update(f);
+        if (get_cos_v3v3(f->no, normalVector)<0)
+            BM_face_normal_flip(bm, f);
+    }
+	return f;
+}
+/*
+This function create polygons between edge and vert
+*/
+BMFace* bmo_edge_vert_connect(BMesh *bm, BridgeParams *bp, BMEdge *e, BMVert *v, BMFace *f_example)
+{
+    BMVert *v1, *v2, *v3, *v4, *vv[3];
+	BMFace *f = NULL;
+    float normalA[3], normalB[3], normalVector[3];
+    int i;
+    v1 = e->v1;
+    v2 = v;
+    v3 = NULL;
+    v4 = BM_edge_other_vert(e, e->v1);
+
+    // calculate normal vector
+       mid_v3_v3v3(normalA, bp->centrod1, bp->centrod2);
+       vv[0] = v1;
+       vv[1] = v2;
+       vv[2] = v4;
+       get_centroid(vv, 3, normalB);
+       sub_v3_v3v3(normalVector, normalB, normalA);
+
+    if (bp->seg > 1){
+        for (i = 1; i < bp->seg; i++){
+            if (bp->inter == LINEAR_INTER)
+            {
+                v2 = get_linear_seg(bm, bp, i, v, e->v1);
+                v3 = get_linear_seg(bm, bp, i, v, BM_edge_other_vert(e, e->v1));
+            }
+            if (bp->inter == CUBIC_INTER)
+            {
+                v2 = get_cubic_seg (bm, bp, i, v,  e->v1);
+                v3 = get_cubic_seg (bm, bp, i, v, BM_edge_other_vert(e, e->v1));
+            }
+            if (i == 1)
+                f = BM_face_create_quad_tri(bm, v2, v, v3, NULL, f_example, TRUE);
+            else
+                f = BM_face_create_quad_tri(bm, v1, v2, v3, v4, f_example, TRUE);
+
+            BM_face_normal_update(f);
+            if (get_cos_v3v3(f->no, normalVector)<0)
+                BM_face_normal_flip(bm, f);
+            v1 = v2;
+            v4 = v3;
+        }
+        f = BM_face_create_quad_tri(bm, e->v1, v2, v3, BM_edge_other_vert(e, e->v1), f_example, TRUE);
+        BM_face_normal_update(f);
+        if (get_cos_v3v3(f->no, normalVector)<0)
+            BM_face_normal_flip(bm, f);
+
+    }
+    else{
+        f = BM_face_create_quad_tri(bm, e->v1, BM_edge_other_vert(e, e->v1), v, NULL, f_example, TRUE);
+        BM_face_normal_update(f);
+        if (get_cos_v3v3(f->no, normalVector)<0)
+            BM_face_normal_flip(bm, f);
+    }
+	return f;
+}
+
+
+
+int bridge_check_planar_point(float p1[3], float p2[3], float p3[3], float v[3])
+{
+	float P[3], N[3];
+	float vv1[3], vv2[3];
+	int result = 0;
+
+	/* calculation of the normal to the surface */
+	sub_v3_v3v3(vv1, p1, p2);
+	sub_v3_v3v3(vv2, p3, p2);
+	cross_v3_v3v3(N, vv1, vv2);
+
+	copy_v3_v3(P, p2);
+
+	if (fabs(N[0] * (v[0]-P[0]) + N[1] * (v[1]-P[1]) + N[2] * (v[2]-P[2])) < BRIDGE_EPSILON){
+		/* point located on plane */
+		result = 1;
+	}
+	return result;
+}
+
+int check_planar_faces(BMFace* fa, BMFace *fb)
+{
+	int result = 0;
+	float p1[3], p2[3], p3[3];
+	BMLoop *l = fa->l_first;
+	copy_v3_v3(p1, l->v->co);
+	l = l->next;
+	copy_v3_v3(p2, l->v->co);
+	l = l->next;
+	copy_v3_v3(p3, l->v->co);
+
+	l = fb->l_first;
+	do {
+		if (bridge_check_planar_point(p1, p2, p3, l->v->co))
+			result = 1;
+		else{
+			result = 0;
+			break;
+		}
+		l = l->next;
+	} while  (l != fb->l_first);
+	return result;
+
+}
+
+/*
+* return 1 if face content e
+*/
+int bridge_is_edge_of_face (BMFace *f, BMEdge* e)
+{
+	int result  = 0;
+	BMLoop *l = f->l_first;
+	do {
+		if (l->e == e)
+			result = 1;
+		l = l->next;
+	} while  (l != f->l_first);
+	return result;
+}
+
+int check_content_faces(BMFace *f, BMFace** faces, int count)
+{
+	int i, result = 0;
+	for(i = 0; i< count; i++){
+		if (f == faces[i])
+			result = 1;
+	}
+	return result;
+}
+
+void connect_faces(BMesh *bm, BMFace *fa, BMFace* fb, BMEdge *e)
+{
+	BMVert **vv = NULL;
+	BMLoop *l;
+	int count;
+	BLI_array_declare(vv);
+	BLI_array_append(vv, e->v1);
+	BLI_array_append(vv, BM_edge_other_vert(e, e->v1));
+
+	l = fa->l_first;
+	do {
+		if ((l->v != e->v1) && (l->v != BM_edge_other_vert(e, e->v1)))
+			BLI_array_append(vv, l->v);
+		l = l->next;
+	} while  (l != fa->l_first);
+
+	l = fb->l_first;
+	do {
+		if ((l->v != e->v1) && (l->v != BM_edge_other_vert(e, e->v1)))
+			BLI_array_append(vv, l->v);
+		l = l->next;
+	} while  (l != fb->l_first);
+
+	count = BLI_array_count(vv);
+	if (count > 2){
+		BM_face_create_ngon_vcloud(bm, vv, count, 0);
+	}
+	BLI_array_free(vv);
+}
+
+/*
+  1 - unic element
+  0 - dublicated element
+  */
+int check_unic_edge(BMEdge **ee, int count, BMEdge*e)
+{
+	int i, result = 1;
+	for (i = 0; i < count; i++) {
+		if (ee[i] == e)
+			result = 0;
+	}
+	return result;
+}
+
+void bridge_n_gon_optimization(BMesh *bm, BMFace **faces, int count)
+{
+	int i = 0, j = 0, t_count = 0 ;
+	BMFace *f , **exp_faces = NULL;
+	BMEdge *e = NULL, **ee = NULL;
+	BMLoop *l;
+	BLI_array_declare(exp_faces);
+	BLI_array_declare(ee);
+	for (i = 0; i < count; i++){
+		l = faces[i]->l_first;
+		f = NULL;
+		do {
+			for (j = 0; j < count; j++){
+				if (!check_content_faces(faces[j], exp_faces,  BLI_array_count(exp_faces))){
+					if (faces[i] != faces[j]){
+						if (bridge_is_edge_of_face(faces[j], l->e)){
+							if (check_planar_faces(faces[i], faces[j])){
+								connect_faces(bm, faces[i],faces[j], l->e);
+								f = faces[i];
+								e = l->e;
+								break;
+							}
+						}
+					}
+				}
+			}
+			l = l->next;
+		} while  (l != faces[i]->l_first);
+
+		if ( f != NULL){
+			BLI_array_append(exp_faces, f);
+			if (check_unic_edge(ee, BLI_array_count(ee), e))
+				BLI_array_append(ee, e);
+		}
+	}
+
+	t_count = BLI_array_count(ee);
+	if(t_count > 0 ) {
+		for (i = 0; i < t_count; i++ )
+			BM_edge_kill(bm, ee[i]);
+	}
+
+	BLI_array_free(exp_faces);
+}
+
+void  bmo_bridge_loops_exec(BMesh *bm, BMOperator *op)
+{
+    BMEdge *e;
+    BMOIter siter;
+    BMIter iter;
+    BMEdge **skip_edge = NULL, **all_edge = NULL, **ee1=NULL, **ee2=NULL;
+    BMVert **vv1 = NULL, **vv2 = NULL;
+    BridgeParams bp;
+
+    int i = 0, j = 0, loops_count = 0, cl1 = 0, cl2 = 0;
+	BMFace *f, **face_list = NULL;
+
+    BLI_array_declare(skip_edge);
+    BLI_array_declare(all_edge);
+    BLI_array_declare(ee1);
+    BLI_array_declare(ee2);
+    BLI_array_declare(vv1);
+    BLI_array_declare(vv2);
+	BLI_array_declare(face_list);
+
+    // init bridge param
+    bp.newVertices.first = bp.newVertices.last = NULL;
+    bp.seg = BMO_slot_int_get(op, "segmentation");
+    if (BMO_slot_int_get(op,"interpolation"))
+        bp.inter = CUBIC_INTER;
+    else
+        bp.inter = LINEAR_INTER;
+    bp.strenght = BMO_slot_float_get(op, "strenght");
+
+    // find all edges in faces
+    BMO_ITER (f, &siter, bm, op, "edgefacein", BM_FACE)
+    {
+       BMLoop* loop;
+       loop = f->l_first;
+       do
+       {
+           BLI_array_append(all_edge, loop->e);
+           loop = loop->next;
+       }
+       while (loop != f->l_first);
+    }
+    // count entry edges
+    for (i = 0; i< BLI_array_count(all_edge); i++)
+    {
+        int counter = 0;
+        for (j = 0; j< BLI_array_count(all_edge); j++)
+        {
+            if ((all_edge[i] == all_edge[j]) && (i != j))
+                counter ++;
+        }
+        if ((counter > 0) &&
+            (bmo_bridge_array_entry_elem(skip_edge,
+                                         BLI_array_count(skip_edge),
+                                         all_edge[i]) == 0 ))
+        {
+            BLI_array_append(skip_edge, all_edge[i]);
+            BMO_elem_flag_enable(bm, all_edge[i], EDGE_DONE); // skiped adjacent faces edge
+
+        }
+    }
+    //TODO develop non connected case
+    // create loop according with skiped edges
+    BMO_slot_buffer_flag_enable(bm, op, "edgefacein", BM_EDGE, EDGE_MARK);
+    BMO_ITER (e, &siter, bm, op, "edgefacein", BM_EDGE)
+    {
+        if (!BMO_elem_flag_test(bm, e, EDGE_DONE))
+        {
+        BMVert *v, *ov;
+        BMEdge *e2, *e3;
+
+        if (loops_count > 2)
+        {
+            BMO_error_raise(bm, op, BMERR_INVALID_SELECTION, "Select only two edge loops");
+            goto cleanup_2;
+        }
+            v = e->v1;
+            e2 = e;
+            ov = v;
+            do
+            {
+                if (loops_count == 0)
+                {
+                    BLI_array_append(ee1, e2);
+                    BLI_array_append(vv1, v);
+                }
+                else
+                {
+                    BLI_array_append(ee2, e2);
+                    BLI_array_append(vv2, v);
+                }
+
+                BMO_elem_flag_enable(bm, e2, EDGE_DONE);
+
+                v = BM_edge_other_vert(e2, v);
+                BM_ITER_ELEM (e3, &iter, v, BM_EDGES_OF_VERT)
+                {
+                    if ((e3 != e2)
+                        && (BMO_elem_flag_test(bm, e3, EDGE_MARK))
+                        && (!BMO_elem_flag_test(bm, e3, EDGE_DONE)))
+                    {
+                        break;
+                    }
+                }
+                if (e3)
+                    e2 = e3;
+            }
+            while (e3 && e2 != e);
+            /* test for connected loops, and set cl1 or cl2 if so */
+            if (v == ov) {
+                if (loops_count == 0) {
+                    cl1 = 1;
+                }
+                else {
+                    cl2 = 1;
+                }
+            }
+            loops_count ++;
+        }
+    }
+// -----------------------------------------------
+// ---------------CREATE BRIDGE-------------------
+// -----------------------------------------------
+    if (ee1 && ee2)
+    {
+       float min = 1e32;
+       int min_j;
+       float center1[3], center2[3];
+       BMEdge **non_conected = NULL, **conected = NULL;
+       BLI_array_declare(non_conected);
+       BLI_array_declare(conected);
+
+       if(BLI_array_count(ee1) > BLI_array_count(ee2))
+       {
+           ARRAY_SWAP(BMVert *, vv1, vv2);
+           ARRAY_SWAP(BMEdge *, ee1, ee2);
+       }
+
+       get_centroid(vv1, BLI_array_count(vv1), center1);
+       get_centroid(vv2, BLI_array_count(vv2), center2);
+
+       copy_v3_v3(bp.centrod1, center1);
+       copy_v3_v3(bp.centrod2, center2);
+
+       for (i = 0; i < BLI_array_count(ee1); i++)
+            BMO_elem_flag_enable(bm, ee1[i], EDGE_NON_CONNECTED);
+       for (i = 0; i < BLI_array_count(ee2); i++)
+            BMO_elem_flag_enable(bm, ee2[i], EDGE_NON_CONNECTED);
+
+       for (i = 0; i < BLI_array_count(ee1); i++)
+       {
+           float mid_V_i[3], mid_V_j[3]/*, centerV[3]*/;
+           float co1[3], angel;
+
+           get_middle_edge(ee1[i], mid_V_i);
+           sub_v3_v3v3(co1, mid_V_i, bp.centrod1);
+         /*  sub_v3_v3v3(centerV, bp.centrod2, bp.centrod1);
+           angel = angle_v3v3(centerV, co1);
+           mul_v3_fl(co1, cos(angel - M_PI_2));*/
+           min = 1e32;
+           min_j =  -1;
+           for (j = 0; j < BLI_array_count(ee2); j++)
+           {
+               float co2[3], summ[3]; // new coordinate with center of loops
+               if (BMO_elem_flag_test(bm, ee2[j], EDGE_NON_CONNECTED))
+               {
+                   get_middle_edge(ee2[j],mid_V_j);
+                   sub_v3_v3v3(co2, mid_V_j, bp.centrod2);
+
+                   //---------------------------
+                   /*
+                   sub_v3_v3v3(centerV, bp.centrod1, bp.centrod2);
+                   angel = angle_v3v3(centerV, co2);
+                   mul_v3_fl(co2, cos(angel-M_PI_2));
+                   */
+                   //-----------------
+                   //sub_v3_v3v3(co1, mid_V_i, bp.centrod1);
+                   //sub_v3_v3v3(co2, mid_V_j, bp.centrod2);
+                   //---------------------------------
+                   //sub_v3_v3v3(co1, co1, bp.centrod1);
+                   //sub_v3_v3v3(co2, co2, bp.centrod2);
+                   //---------------------------------
+                   angel = angle_v3v3(co1, co2);
+                   add_v3_v3v3(summ, co1, co2);
+                   //angel = angle_v3v3(co1, summ);
+                   if (len_v3v3(co1, co2) < min)
+                   {
+                       min = len_v3v3(co1, co2);
+                       min_j = j;
+                   }
+               }
+           }
+           if (min_j != -1){
+			   f = bmo_edge_face_connect(bm, &bp, ee1[i], ee2[min_j], NULL);
+			   if (f != NULL)
+				   BLI_array_append(face_list, f);
+
+               BLI_array_append(conected, ee2[min_j]);
+               BMO_elem_flag_enable(bm, ee2[min_j], EDGE_CONNECTED);
+           }
+       }
+
+
+        // create triangle
+       for (i = 0; i < BLI_array_count(ee2); i++)
+       {
+           int flag = 0;
+           for (j = 0; j < BLI_array_count(conected); j++)
+           {
+               if (ee2[i] == conected[j])
+                   flag = 1;
+            }
+           if (flag  == 0)
+               BLI_array_append(non_conected, ee2[i]);
+
+        }
+
+       // найти ближающую точку в vv1; и приконектить...
+        for (i = 0; i< BLI_array_count(non_conected); i++)
+        {
+            //if (BMO_elem_flag_test(bm, ee2[i], EDGE_CONNECTED))
+            //{
+                float  mid[3];
+                float co1[3], co2[3]; // new coordinate with center of loops
+
+                get_middle_edge(non_conected[i], mid);
+                sub_v3_v3v3(co1, mid, center2);
+                min = 1e32;
+                min_j = 0;
+                for (j = 0; j < BLI_array_count(vv1); j++) {
+                    sub_v3_v3v3(co2, vv1[j]->co, center1);
+                    if (len_v3v3(co1, co2) < min)
+                    {
+                        min = len_v3v3(co1, co2);
+                        min_j = j;
+                    }
+                }
+				f = bmo_edge_vert_connect(bm, &bp,  non_conected[i], vv1[min_j], NULL);
+				if (f != NULL)
+					BLI_array_append(face_list, f);
+            //}
+        }
+        BLI_array_free(non_conected);
+        BLI_array_free(conected);
+    }
+
+// ---------------------------------------------
+	if (BMO_slot_bool_get(op,"n_gon"))
+		bridge_n_gon_optimization(bm, face_list, BLI_array_count(face_list));
+// -----kill input data-------------------------
+    BMO_ITER (f, &siter, bm, op, "edgefacein", BM_FACE)
+    {
+        BM_face_kill(bm,f);
+    }
+    for (i=0; i< BLI_array_count(skip_edge); i++){
+        BM_edge_kill(bm, skip_edge[i]);
+    }
+
+cleanup_2:
+    BLI_array_free(ee1);
+    BLI_array_free(ee2);
+    BLI_array_free(vv1);
+    BLI_array_free(vv2);
+    BLI_array_free(skip_edge);
+    BLI_array_free(all_edge);
+    BLI_freelistN(&bp.newVertices);
+
+}
+
+
+void bmo_bridge_loops_exec_old(BMesh *bm, BMOperator *op)
 {
 	BMEdge **ee1 = NULL, **ee2 = NULL;
 	BMVert **vv1 = NULL, **vv2 = NULL;

source/blender/bmesh/operators/bmo_smooth_laplacian.c

@@ -0,0 +1,613 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * Contributor(s): Alexander Pinzon
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/bmesh/operators/bmo_smooth_laplacian.c
+ *  \ingroup bmesh
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_meshdata_types.h"
+
+#include "BLI_array.h"
+#include "BLI_heap.h"
+#include "BLI_math.h"
+#include "BLI_math_geom.h"
+#include "BLI_smallhash.h"
+
+#include "BKE_customdata.h"
+#include "BKE_mesh.h"
+
+#include "bmesh.h"
+
+#include "ONL_opennl.h"
+
+#include "intern/bmesh_operators_private.h" /* own include */
+
+#define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f
+#define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f
+#define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8
+#define SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE 0.15
+
+struct BLaplacianSystem {
+	float *eweights;		/* Length weights per Edge */
+	float (*fweights)[3];   /* Cotangent weights per face */
+	float *ring_areas;		/* Total area per ring*/
+	float *vlengths;		/* Total sum of lengths(edges) per vertice*/
+	float *vweights;		/* Total sum of weights per vertice*/
+	int numEdges;			/* Number of edges*/
+	int numFaces;			/* Number of faces*/
+	int numVerts;			/* Number of verts*/
+	short *zerola;			/* Is zero area or length*/
+
+	/* Pointers to data*/
+	BMesh *bm;
+	BMOperator *op;
+	NLContext *context;
+
+	/*Data*/
+	float min_area;
+};
+typedef struct BLaplacianSystem LaplacianSystem;
+
+static float compute_volume(BMesh *bm, BMOperator *op);
+static float cotan_weight(float *v1, float *v2, float *v3);
+static int vert_is_boundary(BMVert *v);
+static LaplacianSystem * init_laplacian_system( int a_numEdges, int a_numFaces, int a_numVerts);
+static void init_laplacian_matrix(LaplacianSystem * sys);
+static void delete_laplacian_system(LaplacianSystem * sys);
+static void delete_void_pointer(void * data);
+static void fill_laplacian_matrix(LaplacianSystem * sys);
+static void memset_laplacian_system(LaplacianSystem *sys, int val);
+static void validate_solution(LaplacianSystem * sys, int usex, int usey, int usez, float lambda, float lambda_border, int volumepreservation);
+static void volume_preservation(BMesh *bm, BMOperator *op, float vini, float vend, int usex, int usey, int usez);
+
+static void delete_void_pointer(void * data)
+{
+	if (data) {
+		MEM_freeN(data);
+		data = NULL;
+	}
+}
+
+static void delete_laplacian_system(LaplacianSystem * sys) 
+{
+	delete_void_pointer(sys->eweights);
+	delete_void_pointer(sys->fweights);
+	delete_void_pointer(sys->ring_areas);
+	delete_void_pointer(sys->vlengths);
+	delete_void_pointer(sys->vweights);
+	delete_void_pointer(sys->zerola);
+	if (sys->context) {
+		nlDeleteContext(sys->context);
+	}
+	sys->bm = NULL;
+	sys->op = NULL;
+	MEM_freeN(sys);
+}
+
+static void memset_laplacian_system(LaplacianSystem *sys, int val)
+{
+	memset(sys->eweights	, val, sizeof(float) * sys->numEdges);
+	memset(sys->fweights	, val, sizeof(float) * sys->numFaces * 3);
+	memset(sys->ring_areas	, val, sizeof(float) * sys->numVerts);
+	memset(sys->vlengths	, val, sizeof(float) * sys->numVerts);
+	memset(sys->vweights	, val, sizeof(float) * sys->numVerts);
+	memset(sys->zerola		, val, sizeof(short) * sys->numVerts);
+}
+
+static LaplacianSystem * init_laplacian_system( int a_numEdges, int a_numFaces, int a_numVerts) 
+{
+	LaplacianSystem * sys; 
+	sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
+	sys->numEdges = a_numEdges;
+	sys->numFaces = a_numFaces;
+	sys->numVerts = a_numVerts;
+
+	sys->eweights =  MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
+	if (!sys->eweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->fweights =  MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
+	if (!sys->fweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->ring_areas =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
+	if (!sys->ring_areas) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->vlengths =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
+	if (!sys->vlengths) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	sys->vweights =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
+	if (!sys->vweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	sys->zerola =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
+	if (!sys->zerola) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	return sys;
+}
+
+/* Compute weigth between vertice v_i and all your neighbors
+ * weight between v_i and v_neighbor 
+ * Wij = cot(alpha) + cot(beta) / (4.0 * total area of all faces  * sum all weight)
+ *        v_i *
+ *          / | \
+ *         /  |  \
+ *  v_beta*   |   * v_alpha
+ *         \  |  /
+ *          \ | /
+ *            * v_neighbor
+*/
+
+static void init_laplacian_matrix(LaplacianSystem * sys)
+{
+	float areaf;
+	float *v1, *v2, *v3, *v4;
+	float w1, w2, w3, w4;
+	int i, j;
+	int has_4_vert ;
+	unsigned int idv1, idv2, idv3, idv4, idv[4];
+	BMEdge *e;
+	BMFace *f;
+	BMIter eiter;
+	BMIter fiter;
+	BMIter vi;
+	BMVert *vn;
+	BMVert *vf[4];
+
+	BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) {
+		if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) {
+			v1 = e->v1->co;
+			v2 =  e->v2->co;
+			idv1 = BM_elem_index_get(e->v1);
+			idv2 = BM_elem_index_get(e->v2);
+			
+			w1 = len_v3v3(v1, v2);
+			if (w1 > sys->min_area) {
+				w1 = 1.0f / w1;
+				i = BM_elem_index_get(e);
+				sys->eweights[i] = w1;
+				sys->vlengths[idv1] += w1;
+				sys->vlengths[idv2] += w1;
+			}else{
+				sys->zerola[idv1] = 1;
+				sys->zerola[idv2] = 1;
+			}
+		}
+	}
+
+	BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
+		if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+
+			BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+				vf[i] = vn;
+			}
+			has_4_vert = (i == 4) ? 1 : 0;
+			idv1 = BM_elem_index_get(vf[0]);
+			idv2 = BM_elem_index_get(vf[1]);
+			idv3 = BM_elem_index_get(vf[2]);
+			idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0;
+
+			v1 = vf[0]->co;
+			v2 = vf[1]->co;
+			v3 = vf[2]->co;
+			v4 = has_4_vert ? vf[3]->co : 0;
+
+			if (has_4_vert) {
+				areaf = area_quad_v3(v1, v2, v3, v4);
+			} else {
+				areaf = area_tri_v3(v1, v2, v3);
+			}
+
+			if (fabs(areaf) < sys->min_area) { 
+				sys->zerola[idv1] = 1;
+				sys->zerola[idv2] = 1;
+				sys->zerola[idv3] = 1;
+				if (has_4_vert) sys->zerola[idv4] = 1;
+			}
+
+			sys->ring_areas[idv1] += areaf;
+			sys->ring_areas[idv2] += areaf;
+			sys->ring_areas[idv3] += areaf;
+			if (has_4_vert) sys->ring_areas[idv4] += areaf;
+
+			if (has_4_vert) {
+			
+				idv[0] = idv1;
+				idv[1] = idv2;
+				idv[2] = idv3;
+				idv[3] = idv4;
+
+				for (j = 0; j < 4; j++) {
+					idv1 = idv[j];
+					idv2 = idv[(j + 1) % 4];
+					idv3 = idv[(j + 2) % 4];
+					idv4 = idv[(j + 3) % 4];
+
+					v1 = vf[j]->co;
+					v2 = vf[(j + 1) % 4]->co;
+					v3 = vf[(j + 2) % 4]->co;
+					v4 = vf[(j + 3) % 4]->co;
+
+					w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+					w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+					w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+	
+					sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
+				}
+			} else {			
+				i = BM_elem_index_get(f);
+
+				w1 = cotan_weight(v1, v2, v3);
+				w2 = cotan_weight(v2, v3, v1);
+				w3 = cotan_weight(v3, v1, v2);
+
+				sys->fweights[i][0] += w1;
+				sys->fweights[i][1] += w2;
+				sys->fweights[i][2] += w3;
+			
+				sys->vweights[idv1] += w2 + w3;
+				sys->vweights[idv2] += w1 + w3;
+				sys->vweights[idv3] += w1 + w2;
+			}
+		}
+	}
+}
+
+static void fill_laplacian_matrix(LaplacianSystem * sys)
+{
+	float *v1, *v2, *v3, *v4;
+	float w2, w3, w4;
+	int i, j;
+	int has_4_vert ;
+	unsigned int idv1, idv2, idv3, idv4, idv[4];
+
+	BMEdge *e;
+	BMFace *f;
+	BMIter eiter;
+	BMIter fiter;
+	BMIter vi;
+	BMVert *vn;
+	BMVert *vf[4];
+
+	BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
+		if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+			BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+				vf[i] = vn;
+			}
+			has_4_vert = (i == 4) ? 1 : 0;
+			if (has_4_vert) {
+				idv[0] = BM_elem_index_get(vf[0]);
+				idv[1] = BM_elem_index_get(vf[1]);
+				idv[2] = BM_elem_index_get(vf[2]);
+				idv[3] = BM_elem_index_get(vf[3]);
+				for (j = 0; j < 4; j++) {
+					idv1 = idv[j];
+					idv2 = idv[(j + 1) % 4];
+					idv3 = idv[(j + 2) % 4];
+					idv4 = idv[(j + 3) % 4];
+
+					v1 = vf[j]->co;
+					v2 = vf[(j + 1) % 4]->co;
+					v3 = vf[(j + 2) % 4]->co;
+					v4 = vf[(j + 3) % 4]->co;
+
+					w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+					w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+					w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+
+					w2 = w2 / 4.0f;
+					w3 = w3 / 4.0f;
+					w4 = w4 / 4.0f;
+		
+					if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) { 
+						nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
+						nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
+						nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
+					}
+				}
+			} else {
+				idv1 = BM_elem_index_get(vf[0]);
+				idv2 = BM_elem_index_get(vf[1]);
+				idv3 = BM_elem_index_get(vf[2]);
+				/* Is ring if number of faces == number of edges around vertice*/
+				i = BM_elem_index_get(f);
+				if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) { 
+					nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+					nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+				}
+				if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) { 
+					nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+					nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+				}
+				if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) { 
+					nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+					nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+				}
+			}
+		}
+	}
+	BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
+		if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e) ) {
+			v1 = e->v1->co;
+			v2 =  e->v2->co;
+			idv1 = BM_elem_index_get(e->v1);
+			idv2 = BM_elem_index_get(e->v2);
+			if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
+				i = BM_elem_index_get(e);
+				nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+				nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+			}
+		}
+	}
+}
+
+static float cotan_weight(float *v1, float *v2, float *v3)
+{
+	float a[3], b[3], c[3], clen;
+
+	sub_v3_v3v3(a, v2, v1);
+	sub_v3_v3v3(b, v3, v1);
+	cross_v3_v3v3(c, a, b);
+
+	clen = len_v3(c);
+
+	if (clen == 0.0f)
+		return 0.0f;
+	
+	return dot_v3v3(a, b) / clen;
+}
+
+static int vert_is_boundary(BMVert *v)
+{
+	BMEdge *ed;
+	BMFace *f;
+	BMIter ei;
+	BMIter fi;
+	BM_ITER_ELEM(ed, &ei, v, BM_EDGES_OF_VERT) {
+		if (BM_edge_is_boundary(ed)) {
+			return 1;
+		}
+	}
+	BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) {
+		if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+			return 1;
+		}
+	}
+	return 0;
+}
+
+static float compute_volume(BMesh *bm, BMOperator *op)
+{
+	float vol = 0.0f;
+	float x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
+	int i;
+	BMFace *f;
+	BMIter fiter;
+	BMIter vi;
+	BMVert *vn;
+	BMVert *vf[4];
+	
+	BM_ITER_MESH (f, &fiter, bm, BM_FACES_OF_MESH) {
+		BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+			vf[i] = vn;
+		}
+		x1 = vf[0]->co[0];
+		y1 = vf[0]->co[1];
+		z1 = vf[0]->co[2];
+
+		x2 = vf[1]->co[0];
+		y2 = vf[1]->co[1];
+		z2 = vf[1]->co[2];
+
+		x3 = vf[2]->co[0];
+		y3 = vf[2]->co[1];
+		z3 = vf[2]->co[2];
+
+		vol += (1.0 / 6.0) * (0.0 - x3*y2*z1 + x2*y3*z1 + x3*y1*z2 - x1*y3*z2 - x2*y1*z3 + x1*y2*z3);
+
+		if (i == 4) {
+			x4 = vf[3]->co[0];
+			y4 = vf[3]->co[1];
+			z4 = vf[3]->co[2];
+			vol += (1.0 / 6.0) * (x1*y3*z4 - x1*y4*z3 - x3*y1*z4 + x3*z1*y4 + y1*x4*z3 - x4*y3*z1);
+		}
+	}
+	return fabs(vol);
+}
+
+static void volume_preservation(BMesh *bm, BMOperator *op, float vini, float vend, int usex, int usey, int usez)
+{
+	float beta;
+	BMOIter siter;
+	BMVert *v;
+
+	if (vend != 0.0f) {	
+		beta  = pow (vini / vend, 1.0f / 3.0f);
+		BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+			if (usex) {
+				v->co[0] *=  beta;
+			}
+			if (usey) {
+				v->co[1] *= beta;
+			}
+			if (usez) {
+				v->co[2] *= beta;
+			}
+			
+		}
+	}
+}
+
+static void validate_solution(LaplacianSystem * sys, int usex, int usey, int usez, float lambda, float lambda_border, int volumepreservation)
+{
+	int m_vertex_id;
+	float leni, lene;
+	float vini, vend;
+	float *vi1, *vi2, ve1[3], ve2[3];
+	unsigned int idv1, idv2;
+	BMOIter siter;
+	BMVert *v;
+	BMEdge *e;
+	BMIter eiter;
+
+	BM_ITER_MESH  (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
+		idv1 = BM_elem_index_get(e->v1);
+		idv2 = BM_elem_index_get(e->v2);
+		vi1 = e->v1->co;
+		vi2 =  e->v2->co;
+		ve1[0] = nlGetVariable(0, idv1);
+		ve1[1] = nlGetVariable(1, idv1);
+		ve1[2] = nlGetVariable(2, idv1);
+		ve2[0] = nlGetVariable(0, idv2);
+		ve2[1] = nlGetVariable(1, idv2);
+		ve2[2] = nlGetVariable(2, idv2);
+		leni = len_v3v3(vi1, vi2);
+		lene = len_v3v3(ve1, ve2);
+		if ( lene > leni* SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni*SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
+			sys->zerola[idv1] = 1;
+			sys->zerola[idv2] = 1;
+		}
+	}
+
+	if (volumepreservation) {
+		vini = compute_volume(sys->bm, sys->op);
+	}
+	BMO_ITER (v, &siter, sys->bm, sys->op, "verts", BM_VERT) {
+		m_vertex_id = BM_elem_index_get(v);
+		if (sys->zerola[m_vertex_id] == 0) {
+			if (usex) {
+				v->co[0] =  nlGetVariable(0, m_vertex_id);
+			}
+			if (usey) {
+				v->co[1] =  nlGetVariable(1, m_vertex_id);
+			}
+			if (usez) {
+				v->co[2] =  nlGetVariable(2, m_vertex_id);
+			}
+		}
+	}
+	if (volumepreservation) {
+		vend = compute_volume(sys->bm, sys->op);
+		volume_preservation(sys->bm, sys->op, vini, vend, usex, usey, usez);
+	}
+
+}
+
+void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
+{
+	int i;
+	int m_vertex_id;
+	int usex, usey, usez, volumepreservation;
+	float lambda, lambda_border;
+	float w;
+	BMOIter siter;
+	BMVert *v;
+	LaplacianSystem * sys;
+
+	sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert);
+	if (!sys) return;
+	sys->bm = bm;
+	sys->op = op;
+
+	memset_laplacian_system(sys, 0);
+
+	BM_mesh_elem_index_ensure(bm, BM_VERT);
+	lambda = BMO_slot_float_get(op, "lambda");
+	lambda_border = BMO_slot_float_get(op, "lambda_border");
+	sys->min_area = 0.00001f;
+	usex = BMO_slot_bool_get(op, "use_x");
+	usey = BMO_slot_bool_get(op, "use_y");
+	usez = BMO_slot_bool_get(op, "use_z");
+	volumepreservation = BMO_slot_bool_get(op, "volume_preservation");
+
+
+	nlNewContext();
+	sys->context = nlGetCurrent();
+
+	nlSolverParameteri(NL_NB_VARIABLES, bm->totvert);
+	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+	nlSolverParameteri(NL_NB_ROWS, bm->totvert);
+	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
+
+	nlBegin(NL_SYSTEM);	
+	for (i=0; i < bm->totvert; i++) {
+		nlLockVariable(i);
+	}
+	BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+		m_vertex_id = BM_elem_index_get(v);
+		nlUnlockVariable(m_vertex_id);
+		nlSetVariable(0,m_vertex_id, v->co[0]);
+		nlSetVariable(1,m_vertex_id, v->co[1]);
+		nlSetVariable(2,m_vertex_id, v->co[2]);
+	}
+
+	nlBegin(NL_MATRIX);
+	init_laplacian_matrix(sys);
+	BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+		m_vertex_id = BM_elem_index_get(v);
+		nlRightHandSideAdd(0, m_vertex_id, v->co[0]);
+		nlRightHandSideAdd(1, m_vertex_id, v->co[1]);
+		nlRightHandSideAdd(2, m_vertex_id, v->co[2]);
+		i = m_vertex_id;
+		if (sys->zerola[i] == 0) {
+			w = sys->vweights[i] * sys->ring_areas[i];
+			sys->vweights[i] = (w == 0.0f) ? 0.0f : - lambda  / (4.0f * w);
+			w = sys->vlengths[i];
+			sys->vlengths[i] = (w == 0.0f) ? 0.0f : - lambda_border  * 2.0f / w;
+
+			if (!vert_is_boundary(v)) { 
+				nlMatrixAdd(i, i,  1.0f + lambda / (4.0f * sys->ring_areas[i]));
+			} else { 
+				nlMatrixAdd(i, i,  1.0f + lambda_border * 2.0f);
+			}
+		} else {
+			nlMatrixAdd(i, i, 1.0f);
+		}	
+	}
+	fill_laplacian_matrix(sys);
+		
+	nlEnd(NL_MATRIX);
+	nlEnd(NL_SYSTEM);
+
+	if (nlSolveAdvanced(NULL, NL_TRUE) ) {
+		validate_solution(sys, usex, usey, usez, lambda, lambda_border, volumepreservation);
+	}
+		
+	delete_laplacian_system(sys);
+}

source/blender/editors/mesh/editmesh_tools.c

@@ -99,7 +99,7 @@ static int edbm_subdivide_exec(bContext *C, wmOperator *op)
 		RNA_enum_set(op->ptr, "quadcorner", SUBD_INNERVERT);
 	}
 	
-	BM_mesh_esubdivide(em->bm, BM_ELEM_SELECT,
+    BM_mesh_esubdivide(em->bm, BM_ELEM_SELECT,
 	                   smooth, fractal, along_normal,
 	                   cuts,
 	                   SUBDIV_SELECT_ORIG, RNA_enum_get(op->ptr, "quadcorner"),
@@ -1607,6 +1607,90 @@ void MESH_OT_vertices_smooth(wmOperatorType *ot)
 	RNA_def_boolean(ot->srna, "zaxis", 1, "Z-Axis", "Smooth along the Z axis");
 }
 
+static int edbm_do_smooth_laplacian_vertex_exec(bContext *C, wmOperator *op)
+{
+	Object *obedit = CTX_data_edit_object(C);
+	BMEditMesh *em = BMEdit_FromObject(obedit);
+	ModifierData *md;
+	int usex = TRUE, usey = TRUE, usez = TRUE, volume_preservation = TRUE;
+	int i, repeat;
+	float lambda = 0.1f;
+	float lambda_border = 0.1f;
+	BMIter fiter;
+	BMFace *f;
+
+	/* Check if select faces are triangles	*/
+	BM_ITER_MESH (f, &fiter, em->bm, BM_FACES_OF_MESH) {
+		if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+			if(f->len > 4) {
+				BKE_report(op->reports, RPT_WARNING, "Selected faces must be triangles or quads");
+				return OPERATOR_CANCELLED;
+			}	
+		}
+	}
+
+	/* mirror before smooth */
+	if (((Mesh *)obedit->data)->editflag & ME_EDIT_MIRROR_X) {
+		EDBM_verts_mirror_cache_begin(em, TRUE);
+	}
+
+	repeat = RNA_int_get(op->ptr, "repeat");
+	lambda = RNA_float_get(op->ptr, "lambda");
+	lambda_border = RNA_float_get(op->ptr, "lambda_border");
+	usex = RNA_boolean_get(op->ptr, "use_x");
+	usey = RNA_boolean_get(op->ptr, "use_y");
+	usez = RNA_boolean_get(op->ptr, "use_z");
+	volume_preservation = RNA_boolean_get(op->ptr, "volume_preservation");
+	if (!repeat)
+		repeat = 1;
+	
+	for (i = 0; i < repeat; i++) {
+		if (!EDBM_op_callf(em, op,
+		                   "smooth_laplacian_vert verts=%hv lambda=%f lambda_border=%f use_x=%b use_y=%b use_z=%b volume_preservation=%b",
+		                   BM_ELEM_SELECT, lambda, lambda_border, usex, usey, usez, volume_preservation))
+		{
+			return OPERATOR_CANCELLED;
+		}
+	}
+
+	/* apply mirror */
+	if (((Mesh *)obedit->data)->editflag & ME_EDIT_MIRROR_X) {
+		EDBM_verts_mirror_apply(em, BM_ELEM_SELECT, 0);
+		EDBM_verts_mirror_cache_end(em);
+	}
+
+	EDBM_update_generic(C, em, TRUE);
+
+	return OPERATOR_FINISHED;
+}
+
+void MESH_OT_vertices_smooth_laplacian(wmOperatorType *ot)
+{
+	/* identifiers */
+	ot->name = "Laplacian Smooth Vertex";
+	ot->description = "Laplacian smooth of selected vertices";
+	ot->idname = "MESH_OT_vertices_smooth_laplacian";
+	
+	/* api callbacks */
+	ot->exec = edbm_do_smooth_laplacian_vertex_exec;
+	ot->poll = ED_operator_editmesh;
+	
+	/* flags */
+	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
+
+	RNA_def_int(ot->srna, "repeat", 1, 1, 200, 
+					"Number of iterations to smooth the mesh", "", 1, 200);
+	RNA_def_float(ot->srna, "lambda", 0.00005f, 0.0000001f, 1000.0f, 
+					"Lambda factor", "", 0.0000001f, 1000.0f);
+	RNA_def_float(ot->srna, "lambda_border", 0.00005f, 0.0000001f, 1000.0f, 
+					"Lambda factor in border", "", 0.0000001f, 1000.0f);
+	RNA_def_boolean(ot->srna, "use_x", 1, "Smooth X Axis", "Smooth object along	X axis");
+	RNA_def_boolean(ot->srna, "use_y", 1, "Smooth Y Axis", "Smooth object along	Y axis");
+	RNA_def_boolean(ot->srna, "use_z", 1, "Smooth Z Axis", "Smooth object along	Z axis");
+	RNA_def_boolean(ot->srna, "volume_preservation", 1, "Preserve Volume", "Apply volume preservation after smooth");
+
+}
+
 /********************** Smooth/Solid Operators *************************/
 
 static void mesh_set_smooth_faces(BMEditMesh *em, short smooth)
@@ -4470,7 +4554,11 @@ typedef struct {
 
 static void edbm_bevel_update_header(wmOperator *op, bContext *C)
 {
+#ifdef OLDBEV
 	static char str[] = "Confirm: Enter/LClick, Cancel: (Esc/RMB), factor: %s, Use Dist (D): %s: Use Even (E): %s";
+#else
+	static char str[] = "Confirm: Enter/LClick, Cancel: (Esc/RMB), factor: %s, segments: %d";
+#endif
 
 	char msg[HEADER_LENGTH];
 	ScrArea *sa = CTX_wm_area(C);
@@ -4482,11 +4570,18 @@ static void edbm_bevel_update_header(wmOperator *op, bContext *C)
 			outputNumInput(&opdata->num_input, factor_str);
 		else
 			BLI_snprintf(factor_str, NUM_STR_REP_LEN, "%f", RNA_float_get(op->ptr, "percent"));
+#ifdef OLDBEV
 		BLI_snprintf(msg, HEADER_LENGTH, str,
 		             factor_str,
 		             RNA_boolean_get(op->ptr, "use_dist") ? "On" : "Off",
-		             RNA_boolean_get(op->ptr, "use_even") ? "On" : "Off"
-		            );
+
+					 RNA_boolean_get(op->ptr, "use_even") ? "On" : "Off",
+					 RNA_float_get(op->ptr, "Amount"),
+					 RNA_int_get(op->ptr, "Segmentation")
+					 );
+#else
+		BLI_snprintf(msg, HEADER_LENGTH, str, factor_str, RNA_int_get(op->ptr, "segmentation"));
+#endif
 
 		ED_area_headerprint(sa, msg);
 	}
@@ -4573,6 +4668,8 @@ static int edbm_bevel_calc(bContext *C, wmOperator *op)
 	int recursion = 1; /* RNA_int_get(op->ptr, "recursion"); */ /* temp removed, see comment below */
 	const int use_even = RNA_boolean_get(op->ptr, "use_even");
 	const int use_dist = RNA_boolean_get(op->ptr, "use_dist");
+	int seg = RNA_int_get(op->ptr, "segmentation");
+	float amount = RNA_float_get(op->ptr, "amount");
 
 	/* revert to original mesh */
 	if (opdata->is_modal) {
@@ -4584,8 +4681,22 @@ static int edbm_bevel_calc(bContext *C, wmOperator *op)
 
 
 		if (!EDBM_op_init(em, &bmop, op,
-		                  "bevel geom=%hev percent=%f lengthlayer=%i use_lengths=%b use_even=%b use_dist=%b",
-		                  BM_ELEM_SELECT, fac, opdata->li, TRUE, use_even, use_dist))
+						  "bevel geom=%hev"
+						  " percent=%f"
+						  " lengthlayer=%i"
+						  " use_lengths=%b"
+						  " use_even=%b"
+						  " use_dist=%b"
+						  " amount=%f"
+						  " segmentation=%i",
+						  BM_ELEM_SELECT,
+						  fac,
+						  opdata->li,
+						  TRUE,
+						  use_even,
+						  use_dist,
+						  amount,
+						  seg))
 		{
 			return 0;
 		}
@@ -4797,14 +4908,40 @@ void MESH_OT_bevel(wmOperatorType *ot)
 	RNA_def_boolean(ot->srna, "use_even", FALSE, "Even",     "Calculate evenly spaced bevel");
 	RNA_def_boolean(ot->srna, "use_dist", FALSE, "Distance", "Interpret the percent in blender units");
 
+	RNA_def_float(ot->srna, "amount", 0.0f,  -FLT_MAX, FLT_MAX, "Amount", "", 0.0f, 1.0f);
+	RNA_def_int(ot->srna, "segmentation", 1,  1, 10000, "Segmentation", "", 1, 100);
+
 }
 
+enum {
+    BRIDGE_LINEAR,
+    BRIDGE_CUBIC
+};
+
 static int edbm_bridge_edge_loops_exec(bContext *C, wmOperator *op)
 {
 	Object *obedit = CTX_data_edit_object(C);
 	BMEditMesh *em = BMEdit_FromObject(obedit);
-	
-	if (!EDBM_op_callf(em, op, "bridge_loops edges=%he", BM_ELEM_SELECT))
+    int seg = RNA_int_get(op->ptr, "Segmentation");
+    int interpolation_flag = 0;
+	int ngon = RNA_boolean_get(op->ptr, "n_gon");
+    float strenght = RNA_float_get(op->ptr, "Strenght");
+    if (RNA_enum_get(op->ptr, "Interpolation") == BRIDGE_LINEAR)
+        interpolation_flag = 0;
+    if (RNA_enum_get(op->ptr, "Interpolation") == BRIDGE_CUBIC)
+        interpolation_flag = 1;
+	if (!EDBM_op_callf(em, op,
+					   "bridge_loops "
+					   "edgefacein=%hfe "
+                       "segmentation=%i "
+                       "interpolation=%i "
+					   "strenght=%f "
+					   "n_gon=%b",
+                       BM_ELEM_SELECT,
+                       seg,
+                       interpolation_flag,
+					   strenght,
+					   ngon))
 		return OPERATOR_CANCELLED;
 	
 	EDBM_update_generic(C, em, TRUE);
@@ -4812,9 +4949,16 @@ static int edbm_bridge_edge_loops_exec(bContext *C, wmOperator *op)
 	return OPERATOR_FINISHED;
 }
 
+static EnumPropertyItem prop_mesh_bridge_interpolation_types[] = {
+    {BRIDGE_LINEAR, "LINEAR", 0, "Linear", ""},
+    {BRIDGE_CUBIC, "CUBIC", 0, "Cubic", ""},
+    {0, NULL, 0, NULL, NULL}
+};
+
 void MESH_OT_bridge_edge_loops(wmOperatorType *ot)
+
 {
-	/* identifiers */
+    /* identifiers */
 	ot->name = "Bridge Two Edge Loops";
 	ot->description = "Make faces between two edge loops";
 	ot->idname = "MESH_OT_bridge_edge_loops";
@@ -4826,7 +4970,15 @@ void MESH_OT_bridge_edge_loops(wmOperatorType *ot)
 	/* flags */
 	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
 	
-	RNA_def_boolean(ot->srna, "inside", 0, "Inside", "");
+	RNA_def_boolean(ot->srna, "n_gon", 0, "N-Gon optimization", "");
+    RNA_def_int(ot->srna, "Segmentation",1,1,1000,"Segmentation", "number of segments",0,50  );
+    RNA_def_enum(ot->srna,
+                 "Interpolation",
+                 prop_mesh_bridge_interpolation_types,
+                 BRIDGE_LINEAR,
+                 "Interpolation",
+                 "Interpolation");
+    RNA_def_float(ot->srna, "Strenght", 0.5, -1.0f, 1.0f, "strenght", "strenght", -1.0f, 1.0f );
 }
 
 typedef struct {

source/blender/editors/mesh/mesh_intern.h

@@ -136,6 +136,7 @@ void MESH_OT_loop_multi_select(struct wmOperatorType *ot);
 void MESH_OT_mark_seam(struct wmOperatorType *ot);
 void MESH_OT_mark_sharp(struct wmOperatorType *ot);
 void MESH_OT_vertices_smooth(struct wmOperatorType *ot);
+void MESH_OT_vertices_smooth_laplacian(struct wmOperatorType *ot);
 void MESH_OT_noise(struct wmOperatorType *ot);
 void MESH_OT_flip_normals(struct wmOperatorType *ot);
 void MESH_OT_solidify(struct wmOperatorType *ot);

source/blender/editors/mesh/mesh_ops.c

@@ -132,6 +132,7 @@ void ED_operatortypes_mesh(void)
 	WM_operatortype_append(MESH_OT_mark_seam);
 	WM_operatortype_append(MESH_OT_mark_sharp);
 	WM_operatortype_append(MESH_OT_vertices_smooth);
+	WM_operatortype_append(MESH_OT_vertices_smooth_laplacian);
 	WM_operatortype_append(MESH_OT_noise);
 	WM_operatortype_append(MESH_OT_flip_normals);
 	//WM_operatortype_append(MESH_OT_knife_cut);

source/blender/editors/transform/transform.c

@@ -401,7 +401,7 @@ void removeAspectRatio(TransInfo *t, float vec[2])
 	}
 }
 
-static void viewRedrawForce(const bContext *C, TransInfo *t)
+void viewRedrawForce(const bContext *C, TransInfo *t)
 {
 	if (t->spacetype == SPACE_VIEW3D) {
 		/* Do we need more refined tags? */
@@ -776,6 +776,7 @@ int transformEvent(TransInfo *t, wmEvent *event)
 	float mati[3][3] = MAT3_UNITY;
 	char cmode = constraintModeToChar(t);
 	int handled = 1;
+	int snapping_handled = 0;
 
 	t->redraw |= handleMouseInput(t, &t->mouse, event);
 
@@ -799,14 +800,20 @@ int transformEvent(TransInfo *t, wmEvent *event)
 		t->redraw |= handleSnapping(t, event);
 	}
 
+	snapping_handled = SnapSystem_Event(t->tsnap.ssystem, event);
+
 	/* handle modal keymap first */
 	if (event->type == EVT_MODAL_MAP) {
 		switch (event->val) {
 			case TFM_MODAL_CANCEL:
-				t->state = TRANS_CANCEL;
+				if(!snapping_handled){
+					t->state = TRANS_CANCEL;
+				}
 				break;
 			case TFM_MODAL_CONFIRM:
-				t->state = TRANS_CONFIRM;
+				if(!snapping_handled){
+					t->state = TRANS_CONFIRM;
+				}
 				break;
 			case TFM_MODAL_TRANSLATE:
 				/* only switch when... */
@@ -1009,7 +1016,9 @@ int transformEvent(TransInfo *t, wmEvent *event)
 	else if (event->val == KM_PRESS) {
 		switch (event->type) {
 			case RIGHTMOUSE:
-				t->state = TRANS_CANCEL;
+				if(!snapping_handled){
+					t->state = TRANS_CANCEL;
+				}
 				break;
 			/* enforce redraw of transform when modifiers are used */
 			case LEFTSHIFTKEY:
@@ -1924,7 +1933,7 @@ void transformApply(bContext *C, TransInfo *t)
 
 	if ((t->redraw & TREDRAW_HARD) || (t->draw_handle_apply == NULL && (t->redraw & TREDRAW_SOFT))) {
 		selectConstraint(t);
-		if (t->transform) {
+		if (t->transform && (SnapSystem_get_state(t->tsnap.ssystem) != SNAPSYSTEM_STATE_INIT_SNAP)) {
 			t->transform(t, t->mval);  // calls recalcData()
 			viewRedrawForce(C, t);
 		}

source/blender/editors/transform/transform.h

@@ -40,6 +40,7 @@
 
 #include "BLI_smallhash.h"
 #include "BKE_tessmesh.h"
+#include "BKE_snap.h"
 
 /* ************************** Types ***************************** */
 
@@ -73,6 +74,7 @@ typedef struct TransSnapPoint {
 } TransSnapPoint;
 
 typedef struct TransSnap {
+	SnapSystem *ssystem;
 	short	mode;
 	short	target;
 	short	modePoint;
@@ -82,6 +84,8 @@ typedef struct TransSnap {
 	char	snap_self;
 	short	peel;
 	short  	status;
+	int r_dist;
+	float r_depth;
 	float	snapPoint[3]; /* snapping from this point */
 	float	snapTarget[3]; /* to this point */
 	float	snapNormal[3];
@@ -462,6 +466,7 @@ void projectFloatView(TransInfo *t, const float vec[3], float adr[2]);
 
 void applyAspectRatio(TransInfo *t, float *vec);
 void removeAspectRatio(TransInfo *t, float *vec);
+void viewRedrawForce(const bContext *C, TransInfo *t);
 
 void initWarp(TransInfo *t);
 int handleEventWarp(TransInfo *t, struct wmEvent *event);
@@ -552,6 +557,24 @@ void drawPropCircle(const struct bContext *C, TransInfo *t);
 
 struct wmKeyMap *transform_modal_keymap(struct wmKeyConfig *keyconf);
 
+/* NOTE: these defines are saved in keymap files, do not change values but just add new ones */
+#define TFM_MODAL_CANCEL		1
+#define TFM_MODAL_CONFIRM		2
+#define TFM_MODAL_TRANSLATE		3
+#define TFM_MODAL_ROTATE		4
+#define TFM_MODAL_RESIZE		5
+#define TFM_MODAL_SNAP_INV_ON	6
+#define TFM_MODAL_SNAP_INV_OFF	7
+#define TFM_MODAL_SNAP_TOGGLE	8
+#define TFM_MODAL_AXIS_X		9
+#define TFM_MODAL_AXIS_Y		10
+#define TFM_MODAL_AXIS_Z		11
+#define TFM_MODAL_PLANE_X		12
+#define TFM_MODAL_PLANE_Y		13
+#define TFM_MODAL_PLANE_Z		14
+#define TFM_MODAL_CONS_OFF		15
+#define TFM_MODAL_ADD_SNAP		16
+#define TFM_MODAL_REMOVE_SNAP	17
 
 /*********************** transform_conversions.c ********** */
 struct ListBase;
@@ -623,6 +646,8 @@ void snapGridAction(TransInfo *t, float *val, GearsType action);
 int activeSnap(TransInfo *t);
 int validSnap(TransInfo *t);
 
+void initSnappingSystem(TransInfo *t, bContext *C);
+void runSnappingSystem(TransInfo *t, float *UNUSED);
 void initSnapping(struct TransInfo *t, struct wmOperator *op);
 void applyProject(TransInfo *t);
 void applySnapping(TransInfo *t, float *vec);

source/blender/editors/transform/transform_generics.c

@@ -1259,6 +1259,7 @@ int initTransInfo(bContext *C, TransInfo *t, wmOperator *op, wmEvent *event)
 	}
 
 	setTransformViewMatrices(t);
+	initSnappingSystem(t, C);
 	initNumInput(&t->num);
 	
 	return 1;
@@ -1327,6 +1328,8 @@ void postTrans(bContext *C, TransInfo *t)
 	if (t->mouse.data) {
 		MEM_freeN(t->mouse.data);
 	}
+
+	SnapSystem_free(t->tsnap.ssystem);
 }
 
 void applyTransObjects(TransInfo *t)

source/blender/editors/transform/transform_snap.c

@@ -61,6 +61,8 @@
 #include "BKE_context.h"
 #include "BKE_tessmesh.h"
 #include "BKE_mesh.h"
+#include "BKE_snap.h"
+#include "BKE_context.h"
 
 #include "ED_armature.h"
 #include "ED_image.h"
@@ -78,6 +80,9 @@
 
 #include "transform.h"
 
+
+//#include "blendef.h" /* for selection modes */
+
 #define USE_BVH_FACE_SNAP
 
 /********************* PROTOTYPES ***********************/
@@ -154,9 +159,10 @@ void drawSnapping(const struct bContext *C, TransInfo *t)
 			RegionView3D *rv3d = CTX_wm_region_view3d(C);
 			float imat[4][4];
 			float size;
-			
+
+			SnapSystem_draw(t->tsnap.ssystem);
 			glDisable(GL_DEPTH_TEST);
-			
+
 			size = 2.5f * UI_GetThemeValuef(TH_VERTEX_SIZE);
 			
 			invert_m4_m4(imat, rv3d->viewmat);
@@ -273,7 +279,7 @@ int  handleSnapping(TransInfo *t, wmEvent *event)
 	}
 #endif
 	if (event->type == MOUSEMOVE) {
-		status |= updateSelectedSnapPoint(t);
+		//status |= updateSelectedSnapPoint(t);
 	}
 	
 	return status;
@@ -429,9 +435,15 @@ static void initSnappingMode(TransInfo *t)
 			/* Exclude editmesh if using proportional edit */
 			if ((obedit->type == OB_MESH) && (t->flag & T_PROP_EDIT)) {
 				t->tsnap.modeSelect = SNAP_NOT_OBEDIT;
+				SnapSystem_set_mode_select(t->tsnap.ssystem, SNAPSYSTEM_MODE_SELECT_NOT_OBEDIT);
 			}
 			else {
 				t->tsnap.modeSelect = t->tsnap.snap_self ? SNAP_ALL : SNAP_NOT_OBEDIT;
+				if(t->tsnap.snap_self){
+					SnapSystem_set_mode_select(t->tsnap.ssystem, SNAPSYSTEM_MODE_SELECT_ALL);
+				}else{
+					SnapSystem_set_mode_select(t->tsnap.ssystem, SNAPSYSTEM_MODE_SELECT_NOT_OBEDIT);
+				}
 			}
 		}
 		/* Particles edit mode*/
@@ -439,12 +451,14 @@ static void initSnappingMode(TransInfo *t)
 		         (obedit == NULL && BASACT && BASACT->object && BASACT->object->mode & OB_MODE_PARTICLE_EDIT))
 		{
 			t->tsnap.modeSelect = SNAP_ALL;
+			SnapSystem_set_mode_select(t->tsnap.ssystem, SNAPSYSTEM_MODE_SELECT_ALL);
 		}
 		/* Object mode */
 		else if (t->tsnap.applySnap != NULL && // A snapping function actually exist
 		         (obedit == NULL) ) // Object Mode
 		{
 			t->tsnap.modeSelect = SNAP_NOT_SELECTED;
+			SnapSystem_set_mode_select(t->tsnap.ssystem, SNAPSYSTEM_MODE_SELECT_NOT_SELECTED);
 		}
 		else {
 			/* Grid if snap is not possible */
@@ -468,6 +482,85 @@ static void initSnappingMode(TransInfo *t)
 	}
 }
 
+//callback with no feedback, to prevent verts from moving from original position..
+//maybe rename to resetfeedback
+void SnappingNFBCallback(SnapSystem *ssystem, void* trans_info){
+	TransInfo* t = (TransInfo*)trans_info;
+	//TODO: find out where the snap point is applied and move it to the other callback functions
+	//below to prevent things moving after thiupdateSelectedSnapPoints function has been called.
+	restoreTransObjects(t); // calls recalcData()
+	viewRedrawForce(SnapSystem_get_C(ssystem), t);
+}
+
+void SnappingConfirmCallback(SnapSystem *ssystem, void* trans_info, SnapPoint sp){
+	TransInfo* t = (TransInfo*)trans_info;
+	copy_v3_v3(t->tsnap.snapPoint, sp.location);
+	copy_v3_v3(t->tsnap.snapNormal, sp.normal);
+	t->state = TRANS_CONFIRM;
+	t->redraw |= updateSelectedSnapPoint(t);
+
+	//TODO: clean this up and same with update function
+	//applyProject(TransInfo *t)
+	//t->tsnap.applySnap
+
+//	copy_v3_v3(t->tsnap.snapPoint, sp->location);
+//	copy_v3_v3(t->tsnap.snapNormal, sp->normal);
+//	if(retval){
+//		t->tsnap.status |=  POINT_INIT;
+//	}
+//	else {
+//		t->tsnap.status &= ~POINT_INIT;
+//	}
+
+}
+
+void SnappingUpdateCallback(SnapSystem *ssystem, void* trans_info, SnapPoint sp){
+	TransInfo* t = (TransInfo*)trans_info;
+	//applyProject(TransInfo *t)
+	//t->tsnap.applySnap does this actuall fully apply it? (is it used while the mouse is moving?)
+	//applySnapping is another good function to look at, same as its callers Translation, Rotation and Scale
+	copy_v3_v3(t->tsnap.snapPoint, sp.location);
+	copy_v3_v3(t->tsnap.snapNormal, sp.normal);
+	t->redraw |= updateSelectedSnapPoint(t);
+//	if(retval){
+//		t->tsnap.status |=  POINT_INIT;
+//	}
+//	else {
+//		t->tsnap.status &= ~POINT_INIT;
+//	}
+}
+
+void SnappingCancelCallback(SnapSystem* UNUSED(ssystem), void* trans_info){
+	/*TODO: fix this function to work properly in cancelling the snapping and resetting it*/
+	//TransInfo* t = (TransInfo*)trans_info;
+	//t->state = TRANS_CANCEL;
+}
+
+
+void runSnappingSystem(TransInfo *t, float *UNUSED){
+	SnapSystem_reset_snappoint(t->tsnap.ssystem);
+	SnapSystem_evaluate_stack(t->tsnap.ssystem);
+
+
+	if(SnapSystem_get_retval(t->tsnap.ssystem)){
+		t->tsnap.status |=  POINT_INIT;
+	}
+	else {
+		t->tsnap.status &= ~POINT_INIT;
+	}
+//	SnapSystem_free(t->tsnap.ssystem);
+}
+
+void initSnappingSystem(TransInfo *t, bContext *C){
+	t->tsnap.ssystem = SnapSystem_create(t->scene, t->view, t->ar, t->obedit, C, t,
+									SnappingUpdateCallback,
+								   SnappingNFBCallback,
+									SnappingConfirmCallback,
+									SnappingCancelCallback);
+//	runSnappingSystem(t, NULL);
+}
+
+
 void initSnapping(TransInfo *t, wmOperator *op)
 {
 	ToolSettings *ts = t->settings;
@@ -526,7 +619,7 @@ void initSnapping(TransInfo *t, wmOperator *op)
 
 static void setSnappingCallback(TransInfo *t)
 {
-	t->tsnap.calcSnap = CalcSnapGeometry;
+	t->tsnap.calcSnap = runSnappingSystem;//CalcSnapGeometry;
 
 	switch (t->tsnap.target) {
 		case SCE_SNAP_TARGET_CLOSEST:
@@ -541,7 +634,6 @@ static void setSnappingCallback(TransInfo *t)
 		case SCE_SNAP_TARGET_ACTIVE:
 			t->tsnap.targetSnap = TargetSnapActive;
 			break;
-
 	}
 
 	switch (t->mode) {
@@ -805,8 +897,8 @@ static void UNUSED_FUNCTION(CalcSnapGrid) (TransInfo * t, float *UNUSED(vec))
 static void CalcSnapGeometry(TransInfo *t, float *UNUSED(vec))
 {
 	if (t->spacetype == SPACE_VIEW3D) {
-		float loc[3];
-		float no[3];
+		float loc[3]; //temporary snapPoint location
+		float no[3]; //snapNormal
 		float mval[2];
 		int found = 0;
 		int dist = SNAP_MIN_DISTANCE; // Use a user defined value here
@@ -910,8 +1002,8 @@ static void CalcSnapGeometry(TransInfo *t, float *UNUSED(vec))
 				copy_v3_v3(t->tsnap.snapTangent, tangent);
 			}
 			
-			copy_v3_v3(t->tsnap.snapPoint, loc);
-			copy_v3_v3(t->tsnap.snapNormal, no);
+			//copy_v3_v3(t->tsnap.snapPoint, loc);
+//			copy_v3_v3(t->tsnap.snapNormal, no); //disable current snapping system
 
 			t->tsnap.status |=  POINT_INIT;
 		}
@@ -1135,7 +1227,7 @@ static void TargetSnapClosest(TransInfo *t)
 }
 /*================================================================*/
 #ifndef USE_BVH_FACE_SNAP
-static int snapFace(ARegion *ar, float v1co[3], float v2co[3], float v3co[3], float *v4co, float mval[2], float ray_start[3], float ray_start_local[3], float ray_normal_local[3], float obmat[][4], float timat[][3], float loc[3], float no[3], int *dist, float *depth)
+//static int snapFace(ARegion *ar, float v1co[3], float v2co[3], float v3co[3], float *v4co, float mval[2], float ray_start[3], float ray_start_local[3], float ray_normal_local[3], float obmat[][4], float timat[][3], float loc[3], float no[3], int *dist, float *depth)
 {
 	float lambda;
 	int result;
@@ -1269,6 +1361,54 @@ static int snapEdge(ARegion *ar, float v1co[3], short v1no[3], float v2co[3], sh
 	return retval;
 }
 
+static int snapEdgeMiddle(ARegion *ar, float v1co[3], short v1no[3], float v2co[3], short v2no[3], float obmat[][4], float timat[][3],
+					const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3], const float mval[2],
+					float r_loc[3], float r_no[3], int *r_dist, float *r_depth)
+{
+	/*TODO: still need to calculate normal output properly using snapEdge code*/
+	float middle[3];
+	int retval = 0;
+	float dvec[3];
+
+
+	mid_v3_v3v3(middle, v1co, v2co); /*find the spot in the middle of the two edge vertices*/
+
+	sub_v3_v3v3(dvec, middle, ray_start_local);
+
+	if (dot_v3v3(ray_normal_local, dvec) > 0) {
+		float location[3];
+		float new_depth;
+		int screen_loc[2];
+		int new_dist;
+
+		copy_v3_v3(location, middle);
+
+		mul_m4_v3(obmat, location);
+
+		new_depth = len_v3v3(location, ray_start);
+
+		project_int(ar, location, screen_loc);
+		new_dist = abs(screen_loc[0] - (int)mval[0]) + abs(screen_loc[1] - (int)mval[1]);
+
+		if (new_dist <= *r_dist && new_depth < *r_depth) {
+			*r_depth = new_depth;
+			retval = 1;
+
+			copy_v3_v3(r_loc, location);
+
+			if (r_no) {
+				normal_short_to_float_v3(r_no, v1no);
+				mul_m3_v3(timat, r_no);
+				normalize_v3(r_no);
+			}
+
+			*r_dist = new_dist;
+		}
+	}
+
+	return retval;
+}
+
 static int snapVertex(ARegion *ar, float vco[3], short vno[3], float obmat[][4], float timat[][3],
                       const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3], const float mval[2],
                       float r_loc[3], float r_no[3], int *r_dist, float *r_depth)
@@ -1584,7 +1724,7 @@ static int snapDerivedMesh(short snap_mode, ARegion *ar, Object *ob, DerivedMesh
 					
 					for (i = 0; i < totedge; i++) {
 						BMEdge *eed = NULL;
-						MEdge *e = edges + i;
+						MEdge *e = edges + i; //is this pointer arithmatic necessary? why not e = edges[i] ?
 						
 						test = 1; /* reset for every vert */
 					
@@ -1602,6 +1742,8 @@ static int snapDerivedMesh(short snap_mode, ARegion *ar, Object *ob, DerivedMesh
 							else {
 								eed = EDBM_edge_at_index(em, index);
 								
+								/*if edge is hidden, or vertex on either end of the edge is selected, then
+								  don't use it to check for snapping*/
 								if (eed && (BM_elem_flag_test(eed, BM_ELEM_HIDDEN) ||
 								            BM_elem_flag_test(eed->v1, BM_ELEM_SELECT) ||
 								            BM_elem_flag_test(eed->v2, BM_ELEM_SELECT)))
@@ -1616,6 +1758,61 @@ static int snapDerivedMesh(short snap_mode, ARegion *ar, Object *ob, DerivedMesh
 						}
 					}
 
+					if (em != NULL) {
+						EDBM_index_arrays_free(em);
+					}
+					break;
+				}
+				case SCE_SNAP_MODE_EDGE_MIDDLE:
+				{
+					MVert *verts = dm->getVertArray(dm);
+					MEdge *edges = dm->getEdgeArray(dm);
+					int totedge = dm->getNumEdges(dm);
+					int *index_array = NULL;
+					int index = 0;
+					int i;
+
+					if (em != NULL) {
+						index_array = dm->getEdgeDataArray(dm, CD_ORIGINDEX);
+						EDBM_index_arrays_init(em, 0, 1, 0);
+					}
+
+					for ( i = 0; i < totedge; i++) {
+						BMEdge *eed = NULL;
+						MEdge *e = edges + i;
+
+						test = 1; /* reset for every vert */
+
+						if (em != NULL) {
+							if (index_array) {
+								index = index_array[i];
+							}
+							else {
+								index = i;
+							}
+
+							if (index == ORIGINDEX_NONE) {
+								test = 0;
+							}
+							else {
+								eed = EDBM_edge_at_index(em, index);
+
+								/*if edge is hidden, or vertex on either end of the edge is selected, then
+								  don't use it to check for snapping*/
+								if (eed && (BM_elem_flag_test(eed, BM_ELEM_HIDDEN) ||
+									BM_elem_flag_test(eed->v1, BM_ELEM_SELECT) ||
+									BM_elem_flag_test(eed->v2, BM_ELEM_SELECT)))
+								{
+									test = 0;
+								}
+							}
+						}
+
+						if (test) {
+							retval |= snapEdgeMiddle(ar, verts[e->v1].co, verts[e->v1].no, verts[e->v2].co, verts[e->v2].no, obmat, timat, ray_start, ray_start_local, ray_normal_local, mval, r_loc, r_no, r_dist, r_depth);
+						}
+					}
+
 					if (em != NULL) {
 						EDBM_index_arrays_free(em);
 					}
@@ -1628,7 +1825,7 @@ static int snapDerivedMesh(short snap_mode, ARegion *ar, Object *ob, DerivedMesh
 	return retval;
 } 
 
-static int snapObject(Scene *scene, ARegion *ar, Object *ob, int editobject, float obmat[][4],
+static int snapObject(Scene *scene, View3D *v3d, ARegion *ar, Object *ob, int editobject, float obmat[][4],
                       const float ray_start[3], const float ray_normal[3], const float mval[2],
                       float r_loc[3], float r_no[3], int *r_dist, float *r_depth)
 {
@@ -1649,6 +1846,7 @@ static int snapObject(Scene *scene, ARegion *ar, Object *ob, int editobject, flo
 			dm = mesh_get_derived_final(scene, ob, CD_MASK_BAREMESH);
 		}
 		
+
 		retval = snapDerivedMesh(ts->snap_mode, ar, ob, dm, em, obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, r_depth);
 
 		dm->release(dm);
@@ -1673,7 +1871,7 @@ static int snapObjects(Scene *scene, View3D *v3d, ARegion *ar, Object *obedit, c
 	if (mode == SNAP_ALL && obedit) {
 		Object *ob = obedit;
 
-		retval |= snapObject(scene, ar, ob, 1, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
+		retval |= snapObject(scene, v3d, ar, ob, 1, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
 	}
 
 	/* Need an exception for particle edit because the base is flagged with BA_HAS_RECALC_DATA
@@ -1684,7 +1882,7 @@ static int snapObjects(Scene *scene, View3D *v3d, ARegion *ar, Object *obedit, c
 	base = BASACT;
 	if (base && base->object && base->object->mode & OB_MODE_PARTICLE_EDIT) {
 		Object *ob = base->object;
-		retval |= snapObject(scene, ar, ob, 0, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
+		retval |= snapObject(scene, v3d, ar, ob, 0, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
 	}
 
 	for (base = FIRSTBASE; base != NULL; base = base->next) {
@@ -1703,13 +1901,13 @@ static int snapObjects(Scene *scene, View3D *v3d, ARegion *ar, Object *obedit, c
 				for (dupli_ob = lb->first; dupli_ob; dupli_ob = dupli_ob->next) {
 					Object *dob = dupli_ob->ob;
 					
-					retval |= snapObject(scene, ar, dob, 0, dupli_ob->mat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
+					retval |= snapObject(scene, v3d, ar, dob, 0, dupli_ob->mat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
 				}
 				
 				free_object_duplilist(lb);
 			}
 			
-			retval |= snapObject(scene, ar, ob, 0, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
+			retval |= snapObject(scene, v3d, ar, ob, 0, ob->obmat, ray_start, ray_normal, mval, r_loc, r_no, r_dist, &depth);
 		}
 	}
 	

source/blender/makesdna/DNA_modifier_types.h

@@ -78,6 +78,7 @@ typedef enum ModifierType {
 	eModifierType_DynamicPaint      = 40,
 	eModifierType_Remesh            = 41,
 	eModifierType_Skin              = 42,
+	eModifierType_LaplacianSmooth,
 	NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -1093,4 +1094,17 @@ enum {
 	MOD_SKIN_SMOOTH_SHADING = 1
 };
 
+/* Smooth modifier flags */
+#define MOD_LAPLACIANSMOOTH_X (1<<1)
+#define MOD_LAPLACIANSMOOTH_Y (1<<2)
+#define MOD_LAPLACIANSMOOTH_Z (1<<3)
+#define MOD_LAPLACIANSMOOTH_VOLUME_PRESERVATION (1<<4)
+
+typedef struct LaplacianSmoothModifierData {
+	ModifierData modifier;
+	float lambda, lambda_border, pad1;
+	char defgrp_name[64]; /* MAX_VGROUP_NAME */
+	short flag, repeat;
+} LaplacianSmoothModifierData;
+
 #endif

source/blender/makesdna/DNA_scene_types.h

@@ -1374,6 +1374,9 @@ typedef struct Scene {
 #define SCE_SNAP_MODE_NODE_X	5
 #define SCE_SNAP_MODE_NODE_Y	6
 #define SCE_SNAP_MODE_NODE_XY	7
+#define SCE_SNAP_MODE_EDGE_MIDDLE	8
+#define SCE_SNAP_MODE_EDGE_PARALLEL 9
+#define SCE_SNAP_MODE_PLANAR 10
 
 /* toolsettings->selectmode */
 #define SCE_SELECT_VERTEX	1 /* for mesh */

source/blender/makesrna/RNA_access.h

@@ -477,6 +477,7 @@ extern StructRNA RNA_SmokeDomainSettings;
 extern StructRNA RNA_SmokeFlowSettings;
 extern StructRNA RNA_SmokeModifier;
 extern StructRNA RNA_SmoothModifier;
+extern StructRNA RNA_LaplacianSmoothModifier;
 extern StructRNA RNA_SoftBodyModifier;
 extern StructRNA RNA_SoftBodySettings;
 extern StructRNA RNA_SolidifyModifier;

source/blender/makesrna/intern/rna_modifier.c

@@ -84,6 +84,7 @@ EnumPropertyItem modifier_type_items[] = {
 	{eModifierType_Curve, "CURVE", ICON_MOD_CURVE, "Curve", ""},
 	{eModifierType_Displace, "DISPLACE", ICON_MOD_DISPLACE, "Displace", ""},
 	{eModifierType_Hook, "HOOK", ICON_HOOK, "Hook", ""},
+	{eModifierType_LaplacianSmooth, "LAPLACIANSMOOTH", ICON_MOD_SMOOTH, "Laplacian Smooth", ""},
 	{eModifierType_Lattice, "LATTICE", ICON_MOD_LATTICE, "Lattice", ""},
 	{eModifierType_MeshDeform, "MESH_DEFORM", ICON_MOD_MESHDEFORM, "Mesh Deform", ""},
 	{eModifierType_Shrinkwrap, "SHRINKWRAP", ICON_MOD_SHRINKWRAP, "Shrinkwrap", ""},
@@ -210,6 +211,8 @@ static StructRNA *rna_Modifier_refine(struct PointerRNA *ptr)
 			return &RNA_RemeshModifier;
 		case eModifierType_Skin:
 			return &RNA_SkinModifier;
+		case eModifierType_LaplacianSmooth:
+			return &RNA_LaplacianSmoothModifier;
 		default:
 			return &RNA_Modifier;
 	}
@@ -384,6 +387,12 @@ static void rna_SmoothModifier_vgroup_set(PointerRNA *ptr, const char *value)
 	rna_object_vgroup_name_set(ptr, value, lmd->defgrp_name, sizeof(lmd->defgrp_name));
 }
 
+static void rna_LaplacianSmoothModifier_vgroup_set(PointerRNA *ptr, const char *value)
+{
+	LaplacianSmoothModifierData *lmd = (LaplacianSmoothModifierData *)ptr->data;
+	rna_object_vgroup_name_set(ptr, value, lmd->defgrp_name, sizeof(lmd->defgrp_name));
+}
+
 static void rna_WaveModifier_vgroup_set(PointerRNA *ptr, const char *value)
 {
 	WaveModifierData *lmd = (WaveModifierData *)ptr->data;
@@ -1742,6 +1751,64 @@ static void rna_def_modifier_smooth(BlenderRNA *brna)
 	RNA_def_property_update(prop, 0, "rna_Modifier_update");
 }
 
+static void rna_def_modifier_laplaciansmooth(BlenderRNA *brna)
+{
+	StructRNA *srna;
+	PropertyRNA *prop;
+
+	srna = RNA_def_struct(brna, "LaplacianSmoothModifier", "Modifier");
+	RNA_def_struct_ui_text(srna, "Laplacian Smooth Modifier", "Smoothing effect modifier");
+	RNA_def_struct_sdna(srna, "LaplacianSmoothModifierData");
+	RNA_def_struct_ui_icon(srna, ICON_MOD_SMOOTH);
+
+	prop = RNA_def_property(srna, "use_x", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", MOD_LAPLACIANSMOOTH_X);
+	RNA_def_property_ui_text(prop, "X", "Smooth object along X axis");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "use_y", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", MOD_LAPLACIANSMOOTH_Y);
+	RNA_def_property_ui_text(prop, "Y", "Smooth object along Y axis");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "use_z", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", MOD_LAPLACIANSMOOTH_Z);
+	RNA_def_property_ui_text(prop, "Z", "Smooth object along Z axis");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "volume_preservation", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", MOD_LAPLACIANSMOOTH_VOLUME_PRESERVATION);
+	RNA_def_property_ui_text(prop, "Preserve Volume", "Apply volume preservation after smooth");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "lamb", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "lambda");
+	RNA_def_property_range(prop, -FLT_MAX, FLT_MAX);
+	RNA_def_property_ui_range(prop, 0.0000001, 1000.0, 0.0000001, 8);
+	RNA_def_property_ui_text(prop, "Lambda Factor", "Smooth factor effect");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "lambdaborder", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "lambda_border");
+	RNA_def_property_range(prop, -FLT_MAX, FLT_MAX);
+	RNA_def_property_ui_range(prop, 0.0000001, 1000.0, 0.0000001, 8);
+	RNA_def_property_ui_text(prop, "Lambda border", "Lambda factor in border");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "iterations", PROP_INT, PROP_NONE);
+	RNA_def_property_int_sdna(prop, NULL, "repeat");
+	RNA_def_property_ui_range(prop, 0, 200, 1, 0);
+	RNA_def_property_ui_text(prop, "Repeat", "");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+	
+	prop = RNA_def_property(srna, "vertex_group", PROP_STRING, PROP_NONE);
+	RNA_def_property_string_sdna(prop, NULL, "defgrp_name");
+	RNA_def_property_ui_text(prop, "Vertex Group",
+	                         "Name of Vertex Group which determines influence of modifier per point");
+	RNA_def_property_string_funcs(prop, NULL, NULL, "rna_LaplacianSmoothModifier_vgroup_set");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+}
+
 static void rna_def_modifier_cast(BlenderRNA *brna)
 {
 	StructRNA *srna;
@@ -3356,6 +3423,7 @@ void RNA_def_modifier(BlenderRNA *brna)
 	rna_def_modifier_ocean(brna);
 	rna_def_modifier_remesh(brna);
 	rna_def_modifier_skin(brna);
+	rna_def_modifier_laplaciansmooth(brna);
 }
 
 #endif

source/blender/makesrna/intern/rna_scene.c

@@ -132,7 +132,10 @@ EnumPropertyItem snap_element_items[] = {
 	{SCE_SNAP_MODE_INCREMENT, "INCREMENT", ICON_SNAP_INCREMENT, "Increment", "Snap to increments of grid"},
 	{SCE_SNAP_MODE_VERTEX, "VERTEX", ICON_SNAP_VERTEX, "Vertex", "Snap to vertices"},
 	{SCE_SNAP_MODE_EDGE, "EDGE", ICON_SNAP_EDGE, "Edge", "Snap to edges"},
+    {SCE_SNAP_MODE_EDGE_MIDDLE, "EDGE_MIDDLE", ICON_SNAP_EDGE, "Edge Middle", "Snap to middle of edges"},
+	{SCE_SNAP_MODE_EDGE_PARALLEL, "EDGE_PARALLEL", ICON_SNAP_EDGE, "Edge Parallel", "Snap to the parallel of edge"},
 	{SCE_SNAP_MODE_FACE, "FACE", ICON_SNAP_FACE, "Face", "Snap to faces"},
+	{SCE_SNAP_MODE_PLANAR, "PLANAR", ICON_SNAP_FACE, "Planar", "Snap to plane picked from face"},
 	{SCE_SNAP_MODE_VOLUME, "VOLUME", ICON_SNAP_VOLUME, "Volume", "Snap to volume"},
 	{0, NULL, 0, NULL, NULL}
 };

source/blender/modifiers/CMakeLists.txt

@@ -38,6 +38,7 @@ set(INC
 	../render/extern/include
 	../../../intern/elbeem/extern
 	../../../intern/guardedalloc
+	../../../intern/opennl/extern
 )
 
 set(INC_SYS
@@ -64,6 +65,7 @@ set(SRC
 	intern/MOD_fluidsim_util.c
 	intern/MOD_hook.c
 	intern/MOD_lattice.c
+	intern/MOD_laplaciansmooth.c
 	intern/MOD_mask.c
 	intern/MOD_meshdeform.c
 	intern/MOD_mirror.c

source/blender/modifiers/MOD_modifiertypes.h

@@ -75,6 +75,7 @@ extern ModifierTypeInfo modifierType_WeightVGProximity;
 extern ModifierTypeInfo modifierType_DynamicPaint;
 extern ModifierTypeInfo modifierType_Remesh;
 extern ModifierTypeInfo modifierType_Skin;
+extern ModifierTypeInfo modifierType_LaplacianSmooth;
 
 /* MOD_util.c */
 void modifier_type_init(ModifierTypeInfo *types[]);

source/blender/modifiers/SConscript

@@ -4,7 +4,7 @@ Import ('env')
 sources = env.Glob('intern/*.c')
 
 incs = '. ./intern'
-incs += ' #/intern/guardedalloc #/intern/decimation/extern #/intern/bsp/extern #/intern/elbeem/extern #/extern/glew/include'
+incs += ' #/intern/guardedalloc #/intern/decimation/extern #/intern/bsp/extern #/intern/elbeem/extern #/extern/glew/include #/intern/opennl/extern'
 incs += ' ../render/extern/include ../blenloader ../bmesh'
 incs += ' ../include ../blenlib ../blenfont ../makesdna ../makesrna ../blenkernel ../blenkernel/intern'
 incs += ' ../gpu'

source/blender/modifiers/intern/MOD_laplaciansmooth.c

@@ -0,0 +1,697 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software  Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2005 by the Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Alexander Pinzon
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+/** \file blender/modifiers/intern/MOD_laplaciansmooth.c
+ *  \ingroup modifiers
+ */
+
+
+#include "DNA_meshdata_types.h"
+#include "DNA_object_types.h"
+
+#include "BLI_math.h"
+#include "BLI_utildefines.h"
+#include "BLI_string.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_deform.h"
+#include "BKE_displist.h"
+#include "BKE_mesh.h"
+#include "BKE_modifier.h"
+#include "BKE_object.h"
+#include "BKE_particle.h"
+#include "BKE_tessmesh.h"
+
+#include "MOD_modifiertypes.h"
+#include "MOD_util.h"
+
+#include "ONL_opennl.h"
+
+#define MOD_LAPLACIANSMOOTH_MAX_EDGE_PERCENTAGE 1.8
+#define MOD_LAPLACIANSMOOTH_MIN_EDGE_PERCENTAGE 0.02
+
+struct BLaplacianSystem {
+	float *eweights;		/* Length weights per Edge */
+	float (*fweights)[3];   /* Cotangent weights per face */
+	float *ring_areas;		/* Total area per ring*/
+	float *vlengths;		/* Total sum of lengths(edges) per vertice*/
+	float *vweights;		/* Total sum of weights per vertice*/
+	int numEdges;			/* Number of edges*/
+	int numFaces;			/* Number of faces*/
+	int numVerts;			/* Number of verts*/
+	short *numNeFa;			/* Number of neighboors faces around vertice*/
+	short *numNeEd;			/* Number of neighboors Edges around vertice*/
+	short *zerola;			/* Is zero area or length*/
+
+	/* Pointers to data*/
+	float (*vertexCos)[3];
+	MFace *mfaces;	
+	MEdge *medges;
+	NLContext *context;
+
+	/*Data*/
+	float min_area;
+	float vert_centroid[3];
+};
+typedef struct BLaplacianSystem LaplacianSystem;
+
+static CustomDataMask required_data_mask(Object *UNUSED(ob), ModifierData *md);
+static int is_disabled(ModifierData *md, int UNUSED(useRenderParams));
+static float compute_volume(float (*vertexCos)[3], MFace *mfaces, int numFaces);
+static float cotan_weight(float *v1, float *v2, float *v3);
+static LaplacianSystem * init_laplacian_system( int a_numEdges, int a_numFaces, int a_numVerts);
+static void copy_data(ModifierData *md, ModifierData *target);
+static void delete_laplacian_system(LaplacianSystem * sys);
+static void delete_void_pointer(void * data);
+static void fill_laplacian_matrix(LaplacianSystem * sys);
+static void init_data(ModifierData *md);
+static void init_laplacian_matrix(LaplacianSystem * sys);
+static void memset_laplacian_system(LaplacianSystem *sys, int val);
+static void volume_preservation(LaplacianSystem *sys, float vini, float vend, short flag);
+static void validate_solution(LaplacianSystem * sys, short flag, float lambda, float lambda_border);
+
+static void delete_void_pointer(void * data)
+{
+	if (data) {
+		MEM_freeN(data);
+		data = NULL;
+	}
+}
+
+static void delete_laplacian_system(LaplacianSystem * sys) 
+{
+	delete_void_pointer(sys->eweights);
+	delete_void_pointer(sys->fweights);
+	delete_void_pointer(sys->numNeEd);
+	delete_void_pointer(sys->numNeFa);
+	delete_void_pointer(sys->ring_areas);
+	delete_void_pointer(sys->vlengths);
+	delete_void_pointer(sys->vweights);
+	delete_void_pointer(sys->zerola);
+	if (sys->context) {
+		nlDeleteContext(sys->context);
+	}
+	sys->vertexCos = NULL;
+	sys->mfaces = NULL;	
+	sys->medges = NULL;
+	MEM_freeN(sys);
+}
+
+static void memset_laplacian_system(LaplacianSystem *sys, int val)
+{
+	memset(sys->eweights	, val, sizeof(float) * sys->numEdges);
+	memset(sys->fweights	, val, sizeof(float) * sys->numFaces * 3);
+	memset(sys->numNeEd		, val, sizeof(short) * sys->numVerts);
+	memset(sys->numNeFa		, val, sizeof(short) * sys->numVerts);
+	memset(sys->ring_areas	, val, sizeof(float) * sys->numVerts);
+	memset(sys->vlengths	, val, sizeof(float) * sys->numVerts);
+	memset(sys->vweights	, val, sizeof(float) * sys->numVerts);
+	memset(sys->zerola		, val, sizeof(short) * sys->numVerts);
+}
+
+static LaplacianSystem * init_laplacian_system( int a_numEdges, int a_numFaces, int a_numVerts) 
+{
+	LaplacianSystem * sys; 
+	sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
+	sys->numEdges = a_numEdges;
+	sys->numFaces = a_numFaces;
+	sys->numVerts = a_numVerts;
+
+	sys->eweights =  MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
+	if (!sys->eweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->fweights =  MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
+	if (!sys->fweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	sys->numNeEd =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothNumNeEd");
+	if (!sys->numNeEd) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->numNeFa =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothNumNeFa");
+	if (!sys->numNeFa) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->ring_areas =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
+	if (!sys->ring_areas) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+	
+	sys->vlengths =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
+	if (!sys->vlengths) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	sys->vweights =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
+	if (!sys->vweights) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	sys->zerola =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
+	if (!sys->zerola) {
+		delete_laplacian_system(sys);
+		return NULL;
+	}
+
+	return sys;
+}
+
+static void init_data(ModifierData *md)
+{
+	LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
+	smd->lambda = 0.00001f;
+	smd->lambda_border = 0.00005f;
+	smd->repeat = 1;
+	smd->flag = MOD_LAPLACIANSMOOTH_X | MOD_LAPLACIANSMOOTH_Y | MOD_LAPLACIANSMOOTH_Z | MOD_LAPLACIANSMOOTH_VOLUME_PRESERVATION;
+	smd->defgrp_name[0] = '\0';
+}
+
+static void copy_data(ModifierData *md, ModifierData *target)
+{
+	LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
+	LaplacianSmoothModifierData *tsmd = (LaplacianSmoothModifierData *) target;
+
+	tsmd->lambda = smd->lambda;
+	tsmd->lambda_border = smd->lambda_border;
+	tsmd->repeat = smd->repeat;
+	tsmd->flag = smd->flag;
+	BLI_strncpy(tsmd->defgrp_name, smd->defgrp_name, sizeof(tsmd->defgrp_name));
+}
+
+static int is_disabled(ModifierData *md, int UNUSED(useRenderParams))
+{
+	LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
+	short flag;
+
+	flag = smd->flag & (MOD_LAPLACIANSMOOTH_X | MOD_LAPLACIANSMOOTH_Y | MOD_LAPLACIANSMOOTH_Z);
+
+	/* disable if modifier is off for X, Y and Z or if factor is 0 */
+	if ( flag == 0) return 1;
+
+	return 0;
+}
+
+static CustomDataMask required_data_mask(Object *UNUSED(ob), ModifierData *md)
+{
+	LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *)md;
+	CustomDataMask dataMask = 0;
+
+	/* ask for vertexgroups if we need them */
+	if (smd->defgrp_name[0]) dataMask |= CD_MASK_MDEFORMVERT;
+
+	return dataMask;
+}
+
+static float cotan_weight(float *v1, float *v2, float *v3)
+{
+	float a[3], b[3], c[3], clen;
+
+	sub_v3_v3v3(a, v2, v1);
+	sub_v3_v3v3(b, v3, v1);
+	cross_v3_v3v3(c, a, b);
+
+	clen = len_v3(c);
+
+	if (clen == 0.0f)
+		return 0.0f;
+	
+	return dot_v3v3(a, b) / clen;
+}
+
+static float compute_volume(float (*vertexCos)[3], MFace *mfaces, int numFaces)
+{
+	float vol = 0.0f;
+	float x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
+	int i;
+	float *vf[4];
+	for (i = 0; i<numFaces; i++) {
+		vf[0] = vertexCos[mfaces[i].v1];
+		vf[1] = vertexCos[mfaces[i].v2];
+		vf[2] = vertexCos[mfaces[i].v3];
+
+		x1 = vf[0][0];
+		y1 = vf[0][1];
+		z1 = vf[0][2];
+
+		x2 = vf[1][0];
+		y2 = vf[1][1];
+		z2 = vf[1][2];
+
+		x3 = vf[2][0];
+		y3 = vf[2][1];
+		z3 = vf[2][2];
+		
+
+		vol +=  (1.0 / 6.0) * (x2*y3*z1 + x3*y1*z2 - x1*y3*z2 - x2*y1*z3 + x1*y2*z3 - x3*y2*z1);
+		if ((&mfaces[i])->v4) {
+			vf[3] = vertexCos[mfaces[i].v4];
+			x4 = vf[3][0];
+			y4 = vf[3][1];
+			z4 = vf[3][2];
+			vol += (1.0 / 6.0) * (x1*y3*z4 - x1*y4*z3 - x3*y1*z4 + x3*z1*y4 + y1*x4*z3 - x4*y3*z1);
+		}
+	}
+	return fabs(vol);
+}
+
+static void volume_preservation(LaplacianSystem *sys, float vini, float vend, short flag)
+{
+	float beta;
+	int i;
+
+	if (vend != 0.0f) {	
+		beta  = pow (vini / vend, 1.0f / 3.0f);
+		for (i = 0; i < sys->numVerts; i++) {
+			if (flag & MOD_LAPLACIANSMOOTH_X) {
+				sys->vertexCos[i][0] = (sys->vertexCos[i][0] - sys->vert_centroid[0])* beta + sys->vert_centroid[0];
+			}
+			if (flag & MOD_LAPLACIANSMOOTH_Y) {
+				sys->vertexCos[i][1] = (sys->vertexCos[i][1] - sys->vert_centroid[1])* beta + sys->vert_centroid[1];
+			}
+			if (flag & MOD_LAPLACIANSMOOTH_Z) {
+				sys->vertexCos[i][2] = (sys->vertexCos[i][2] - sys->vert_centroid[2])* beta + sys->vert_centroid[2];
+			}
+			
+		}
+	}
+}
+
+static void init_laplacian_matrix(LaplacianSystem * sys)
+{
+	float *v1, *v2, *v3, *v4;
+	float w1, w2, w3, w4;
+	float areaf;
+	int i, j;
+	unsigned int idv1, idv2, idv3, idv4, idv[4];
+	int has_4_vert ;
+	for ( i = 0; i < sys->numEdges; i++) {
+		idv1 = sys->medges[i].v1;
+		idv2 = sys->medges[i].v2;
+
+		v1 = sys->vertexCos[idv1];
+		v2 = sys->vertexCos[idv2];
+
+		sys->numNeEd[idv1] = sys->numNeEd[idv1] + 1;
+		sys->numNeEd[idv2] = sys->numNeEd[idv2] + 1;
+		w1 = len_v3v3(v1, v2);
+		if (w1 < sys->min_area) {
+			sys->zerola[idv1] = 1;
+			sys->zerola[idv2] = 1;
+		} else {
+			w1 = 1.0f / w1;
+		}
+		
+		sys->eweights[i] = w1;
+	}
+	for ( i = 0; i < sys->numFaces; i++) {
+		has_4_vert = ((&sys->mfaces[i])->v4) ? 1 : 0;
+
+		idv1 = sys->mfaces[i].v1;
+		idv2 = sys->mfaces[i].v2;
+		idv3 = sys->mfaces[i].v3;
+		idv4 = has_4_vert ? sys->mfaces[i].v4 : 0;
+		
+		sys->numNeFa[idv1] += 1;
+		sys->numNeFa[idv2] += 1;
+		sys->numNeFa[idv3] += 1;
+		if (has_4_vert) sys->numNeFa[idv4] += 1;
+
+		v1 = sys->vertexCos[idv1];
+		v2 = sys->vertexCos[idv2];
+		v3 = sys->vertexCos[idv3];
+		v4 = has_4_vert ? sys->vertexCos[idv4] : 0;
+		
+		if (has_4_vert) {
+			areaf = area_quad_v3(v1, v2, v3, sys->vertexCos[sys->mfaces[i].v4]);
+		} else {
+			areaf = area_tri_v3(v1, v2, v3);
+		}
+		if (fabs(areaf) < sys->min_area) { 
+			sys->zerola[idv1] = 1;
+			sys->zerola[idv2] = 1;
+			sys->zerola[idv3] = 1;
+			if (has_4_vert) sys->zerola[idv4] = 1;
+		}
+
+		sys->ring_areas[idv1] += areaf;
+		sys->ring_areas[idv2] += areaf;
+		sys->ring_areas[idv3] += areaf;
+		if (has_4_vert) sys->ring_areas[idv4] += areaf;
+
+		if (has_4_vert) {
+			
+			idv[0] = idv1;
+			idv[1] = idv2;
+			idv[2] = idv3;
+			idv[3] = idv4;
+
+			for (j = 0; j < 4; j++) {
+				idv1 = idv[j];
+				idv2 = idv[(j + 1) % 4];
+				idv3 = idv[(j + 2) % 4];
+				idv4 = idv[(j + 3) % 4];
+
+				v1 = sys->vertexCos[idv1];
+				v2 = sys->vertexCos[idv2];
+				v3 = sys->vertexCos[idv3];
+				v4 = sys->vertexCos[idv4];
+
+				w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+				w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+				w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+	
+				sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
+			}
+		} else {			
+			w1 = cotan_weight(v1, v2, v3);
+			w2 = cotan_weight(v2, v3, v1);
+			w3 = cotan_weight(v3, v1, v2);
+
+			sys->fweights[i][0] = sys->fweights[i][0] + w1;
+			sys->fweights[i][1] = sys->fweights[i][1] + w2;
+			sys->fweights[i][2] = sys->fweights[i][2] + w3;
+			
+			sys->vweights[idv1] = sys->vweights[idv1] + w2 + w3;
+			sys->vweights[idv2] = sys->vweights[idv2] + w1 + w3;
+			sys->vweights[idv3] = sys->vweights[idv3] + w1 + w2;
+		}
+	}
+	for ( i = 0; i < sys->numEdges; i++) {
+		idv1 = sys->medges[i].v1;
+		idv2 = sys->medges[i].v2;
+		/* if is boundary, apply scale-dependent umbrella operator only with neighboors in boundary */
+		if (sys->numNeEd[idv1] != sys->numNeFa[idv1] && sys->numNeEd[idv2] != sys->numNeFa[idv2]) { 
+			sys->vlengths[idv1] += sys->eweights[i];
+			sys->vlengths[idv2] += sys->eweights[i];
+		}
+	}
+
+}
+
+static void fill_laplacian_matrix(LaplacianSystem * sys)
+{
+	float *v1, *v2, *v3, *v4;
+	float w2, w3, w4;
+	int i, j;
+	int has_4_vert ;
+	unsigned int idv1, idv2, idv3, idv4, idv[4];
+	
+	for ( i = 0; i < sys->numFaces; i++) {
+		idv1 = sys->mfaces[i].v1;
+		idv2 = sys->mfaces[i].v2;
+		idv3 = sys->mfaces[i].v3;
+		has_4_vert = ((&sys->mfaces[i])->v4) ? 1 : 0;
+
+		if (has_4_vert) {
+			idv[0] = sys->mfaces[i].v1;
+			idv[1] = sys->mfaces[i].v2;
+			idv[2] = sys->mfaces[i].v3;
+			idv[3] = sys->mfaces[i].v4;
+			for (j = 0; j < 4; j++) {
+				idv1 = idv[j];
+				idv2 = idv[(j + 1) % 4];
+				idv3 = idv[(j + 2) % 4];
+				idv4 = idv[(j + 3) % 4];
+
+				v1 = sys->vertexCos[idv1];
+				v2 = sys->vertexCos[idv2];
+				v3 = sys->vertexCos[idv3];
+				v4 = sys->vertexCos[idv4];
+
+				w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+				w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+				w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+
+				w2 = w2 / 4.0f;
+				w3 = w3 / 4.0f;
+				w4 = w4 / 4.0f;
+		
+				if (sys->numNeEd[idv1] == sys->numNeFa[idv1] && sys->zerola[idv1] == 0) { 
+					nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
+					nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
+					nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
+				}
+			}
+		} else {
+			/* Is ring if number of faces == number of edges around vertice*/
+			if (sys->numNeEd[idv1] == sys->numNeFa[idv1] && sys->zerola[idv1] == 0) { 
+				nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+				nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+			}
+			if (sys->numNeEd[idv2] == sys->numNeFa[idv2] && sys->zerola[idv2] == 0) { 
+				nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+				nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+			}
+			if (sys->numNeEd[idv3] == sys->numNeFa[idv3] && sys->zerola[idv3] == 0) { 
+				nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+				nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+			}
+		}
+	}
+
+	for ( i = 0; i < sys->numEdges; i++) {
+		idv1 = sys->medges[i].v1;
+		idv2 = sys->medges[i].v2;
+		/* Is boundary */
+		if (sys->numNeEd[idv1] != sys->numNeFa[idv1] && 
+			sys->numNeEd[idv2] != sys->numNeFa[idv2] && 
+			sys->zerola[idv1] == 0 && 
+			sys->zerola[idv2] == 0) 
+		{
+			nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+			nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+		}	
+	}
+}
+
+static void validate_solution(LaplacianSystem * sys, short flag, float lambda, float lambda_border)
+{
+	int i, idv1, idv2;
+	float leni, lene;
+	float vini, vend;
+	float *vi1, *vi2, ve1[3], ve2[3];
+	if (flag & MOD_LAPLACIANSMOOTH_VOLUME_PRESERVATION) {
+		vini = compute_volume(sys->vertexCos, sys->mfaces, sys->numFaces);
+	}
+	for ( i = 0; i < sys->numEdges; i++) {
+		idv1 = sys->medges[i].v1;
+		idv2 = sys->medges[i].v2;
+		vi1 = sys->vertexCos[idv1];
+		vi2 = sys->vertexCos[idv2];
+		ve1[0] = nlGetVariable(0, idv1);
+		ve1[1] = nlGetVariable(1, idv1);
+		ve1[2] = nlGetVariable(2, idv1);
+		ve2[0] = nlGetVariable(0, idv2);
+		ve2[1] = nlGetVariable(1, idv2);
+		ve2[2] = nlGetVariable(2, idv2);
+		leni = len_v3v3(vi1, vi2);
+		lene = len_v3v3(ve1, ve2);
+		if ( lene > leni*MOD_LAPLACIANSMOOTH_MAX_EDGE_PERCENTAGE || lene < leni*MOD_LAPLACIANSMOOTH_MIN_EDGE_PERCENTAGE ){
+			sys->zerola[idv1] = 1;
+			sys->zerola[idv2] = 1;
+		}
+	}
+	for (i = 0; i < sys->numVerts; i++) {
+		if (sys->zerola[i] == 0) {
+			if (flag & MOD_LAPLACIANSMOOTH_X) {
+				sys->vertexCos[i][0] = nlGetVariable(0, i);
+			}
+			if (flag & MOD_LAPLACIANSMOOTH_Y) {
+				sys->vertexCos[i][1] = nlGetVariable(1, i);
+			}
+			if (flag & MOD_LAPLACIANSMOOTH_Z) {
+				sys->vertexCos[i][2] = nlGetVariable(2, i);
+			}
+		}
+	}
+	if (flag & MOD_LAPLACIANSMOOTH_VOLUME_PRESERVATION) {
+		vend = compute_volume(sys->vertexCos, sys->mfaces, sys->numFaces);
+		volume_preservation(sys, vini, vend, flag);
+	}
+}
+
+static void laplaciansmoothModifier_do(
+        LaplacianSmoothModifierData *smd, Object *ob, DerivedMesh *dm,
+        float (*vertexCos)[3], int numVerts)
+{
+	LaplacianSystem *sys;
+	MDeformVert *dvert = NULL;
+	MDeformVert *dv = NULL;
+	float w, wpaint;
+	int i, iter;
+	int defgrp_index;
+
+	DM_ensure_tessface(dm);
+
+	sys = init_laplacian_system(dm->getNumEdges(dm), dm->getNumTessFaces(dm), numVerts);
+	if(!sys) return;
+
+	sys->mfaces = dm->getTessFaceArray(dm);
+	sys->medges = dm->getEdgeArray(dm);
+	sys->vertexCos = vertexCos;
+	sys->min_area = 0.00001f;
+	modifier_get_vgroup(ob, dm, smd->defgrp_name, &dvert, &defgrp_index);
+
+	sys->vert_centroid[0] = 0.0f;
+	sys->vert_centroid[1] = 0.0f;
+	sys->vert_centroid[2] = 0.0f;
+	for (iter = 0; iter < smd->repeat; iter++) {
+		memset_laplacian_system(sys, 0);
+		nlNewContext();
+		sys->context = nlGetCurrent();
+		nlSolverParameteri(NL_NB_VARIABLES, numVerts);
+		nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+		nlSolverParameteri(NL_NB_ROWS, numVerts);
+		nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
+
+		init_laplacian_matrix(sys);
+
+		nlBegin(NL_SYSTEM);	
+		for (i = 0; i < numVerts; i++) {
+			nlSetVariable(0, i, vertexCos[i][0]);
+			nlSetVariable(1, i, vertexCos[i][1]);
+			nlSetVariable(2, i, vertexCos[i][2]);
+			if (iter == 0) {
+				sys->vert_centroid[0] += vertexCos[i][0];
+				sys->vert_centroid[1] += vertexCos[i][1];
+				sys->vert_centroid[2] += vertexCos[i][2];
+			}
+		}
+		if (iter == 0 && numVerts >0) {
+			sys->vert_centroid[0] = sys->vert_centroid[0] / numVerts;
+			sys->vert_centroid[1] = sys->vert_centroid[1] / numVerts;
+			sys->vert_centroid[2] = sys->vert_centroid[2] / numVerts;
+		}
+
+		nlBegin(NL_MATRIX);
+		dv = dvert;
+		for (i = 0; i < numVerts; i++) {
+			nlRightHandSideAdd(0, i, vertexCos[i][0]);
+			nlRightHandSideAdd(1, i, vertexCos[i][1]);
+			nlRightHandSideAdd(2, i, vertexCos[i][2]);
+			if (dv) {
+				wpaint = defvert_find_weight(dv, defgrp_index);
+				dv++;
+			} else {
+				wpaint = 1.0f;
+			}
+			if (sys->zerola[i] == 0) {
+				w = sys->vweights[i] * sys->ring_areas[i];
+				sys->vweights[i] = (w == 0.0f) ? 0.0f : - smd->lambda * wpaint / (4.0f * w);
+				w = sys->vlengths[i];
+				sys->vlengths[i] = (w == 0.0f) ? 0.0f : - smd->lambda_border * wpaint * 2.0f / w;
+
+				if (sys->numNeEd[i] == sys->numNeFa[i]) { 
+					nlMatrixAdd(i, i,  1.0f + smd->lambda * wpaint / (4.0f * sys->ring_areas[i]));
+				} else { 
+					nlMatrixAdd(i, i,  1.0f + smd->lambda_border * wpaint * 2.0f);
+				}
+			} else {
+				nlMatrixAdd(i, i, 1.0f);
+			}
+		}
+
+		fill_laplacian_matrix(sys);
+
+		nlEnd(NL_MATRIX);
+		nlEnd(NL_SYSTEM);
+
+		if (nlSolveAdvanced(NULL, NL_TRUE)) {
+			validate_solution(sys, smd->flag, smd->lambda, smd->lambda_border);
+		}
+		nlDeleteContext(sys->context);
+		sys->context = NULL;
+	}
+	delete_laplacian_system(sys);
+
+}
+
+static void deformVerts(ModifierData *md, Object *ob, DerivedMesh *derivedData,
+                        float (*vertexCos)[3], int numVerts, ModifierApplyFlag UNUSED(flag))
+{
+	DerivedMesh *dm = get_dm(ob, NULL, derivedData, NULL, 0);
+
+	laplaciansmoothModifier_do((LaplacianSmoothModifierData *)md, ob, dm,
+	                  vertexCos, numVerts);
+
+	if (dm != derivedData)
+		dm->release(dm);
+}
+
+static void deformVertsEM(
+        ModifierData *md, Object *ob, struct BMEditMesh *editData,
+        DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
+{
+	DerivedMesh *dm = get_dm(ob, editData, derivedData, NULL, 0);
+
+	laplaciansmoothModifier_do((LaplacianSmoothModifierData *)md, ob, dm,
+	                  vertexCos, numVerts);
+
+	if (dm != derivedData)
+		dm->release(dm);
+}
+
+
+ModifierTypeInfo modifierType_LaplacianSmooth = {
+	/* name */              "Laplacian Smooth",
+	/* structName */        "LaplacianSmoothModifierData",
+	/* structSize */        sizeof(LaplacianSmoothModifierData),
+	/* type */              eModifierTypeType_OnlyDeform,
+	/* flags */             eModifierTypeFlag_AcceptsMesh |
+	                        eModifierTypeFlag_SupportsEditmode,
+
+	/* copy_data */          copy_data,
+	/* deformVerts */       deformVerts,
+	/* deformMatrices */    NULL,
+	/* deformVertsEM */     deformVertsEM,
+	/* deformMatricesEM */  NULL,
+	/* applyModifier */     NULL,
+	/* applyModifierEM */   NULL,
+	/* init_data */          init_data,
+	/* required_data_mask */  required_data_mask,
+	/* freeData */          NULL,
+	/* is_disabled */        is_disabled,
+	/* updateDepgraph */    NULL,
+	/* dependsOnTime */     NULL,
+	/* dependsOnNormals */	NULL,
+	/* foreachObjectLink */ NULL,
+	/* foreachIDLink */     NULL,
+	/* foreachTexLink */    NULL,
+};

source/blender/modifiers/intern/MOD_util.c

@@ -277,5 +277,6 @@ void modifier_type_init(ModifierTypeInfo *types[])
 	INIT_TYPE(DynamicPaint);
 	INIT_TYPE(Remesh);
 	INIT_TYPE(Skin);
+	INIT_TYPE(LaplacianSmooth);
 #undef INIT_TYPE
 }