style cleanup
[blender-staging.git] / source / blender / modifiers / intern / MOD_laplaciansmooth.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software  Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2005 by the Blender Foundation.
19  * All rights reserved.
20  *
21  * Contributor(s): Alexander Pinzon
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  *
25  */
26
27 /** \file blender/modifiers/intern/MOD_laplaciansmooth.c
28  *  \ingroup modifiers
29  */
30
31
32 #include "DNA_meshdata_types.h"
33 #include "DNA_object_types.h"
34
35 #include "BLI_math.h"
36 #include "BLI_utildefines.h"
37 #include "BLI_string.h"
38
39 #include "MEM_guardedalloc.h"
40
41 #include "BKE_cdderivedmesh.h"
42 #include "BKE_deform.h"
43 #include "BKE_displist.h"
44 #include "BKE_mesh.h"
45 #include "BKE_modifier.h"
46 #include "BKE_object.h"
47 #include "BKE_particle.h"
48 #include "BKE_tessmesh.h"
49
50 #include "MOD_modifiertypes.h"
51 #include "MOD_util.h"
52
53 #include "ONL_opennl.h"
54
55 #define MOD_LAPLACIANSMOOTH_MAX_EDGE_PERCENTAGE 1.8f
56 #define MOD_LAPLACIANSMOOTH_MIN_EDGE_PERCENTAGE 0.02f
57
58 struct BLaplacianSystem {
59         float *eweights;        /* Length weights per Edge */
60         float (*fweights)[3];   /* Cotangent weights per face */
61         float *ring_areas;      /* Total area per ring*/
62         float *vlengths;        /* Total sum of lengths(edges) per vertice*/
63         float *vweights;        /* Total sum of weights per vertice*/
64         int numEdges;           /* Number of edges*/
65         int numFaces;           /* Number of faces*/
66         int numVerts;           /* Number of verts*/
67         short *numNeFa;         /* Number of neighboors faces around vertice*/
68         short *numNeEd;         /* Number of neighboors Edges around vertice*/
69         short *zerola;          /* Is zero area or length*/
70
71         /* Pointers to data*/
72         float (*vertexCos)[3];
73         MFace *mfaces;
74         MEdge *medges;
75         NLContext *context;
76
77         /*Data*/
78         float min_area;
79         float vert_centroid[3];
80 };
81 typedef struct BLaplacianSystem LaplacianSystem;
82
83 static CustomDataMask required_data_mask(Object *UNUSED(ob), ModifierData *md);
84 static int is_disabled(ModifierData *md, int UNUSED(useRenderParams));
85 static float average_area_quad_v3(float *v1, float *v2, float *v3, float *v4);
86 static float compute_volume(float (*vertexCos)[3], MFace *mfaces, int numFaces);
87 static float cotan_weight(float *v1, float *v2, float *v3);
88 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts);
89 static void copy_data(ModifierData *md, ModifierData *target);
90 static void delete_laplacian_system(LaplacianSystem *sys);
91 static void delete_void_pointer(void *data);
92 static void fill_laplacian_matrix(LaplacianSystem *sys);
93 static void init_data(ModifierData *md);
94 static void init_laplacian_matrix(LaplacianSystem *sys);
95 static void memset_laplacian_system(LaplacianSystem *sys, int val);
96 static void volume_preservation(LaplacianSystem *sys, float vini, float vend, short flag);
97 static void validate_solution(LaplacianSystem *sys, short flag, float lambda, float lambda_border);
98
99 static void delete_void_pointer(void *data)
100 {
101         if (data) {
102                 MEM_freeN(data);
103                 data = NULL;
104         }
105 }
106
107 static void delete_laplacian_system(LaplacianSystem *sys)
108 {
109         delete_void_pointer(sys->eweights);
110         delete_void_pointer(sys->fweights);
111         delete_void_pointer(sys->numNeEd);
112         delete_void_pointer(sys->numNeFa);
113         delete_void_pointer(sys->ring_areas);
114         delete_void_pointer(sys->vlengths);
115         delete_void_pointer(sys->vweights);
116         delete_void_pointer(sys->zerola);
117         if (sys->context) {
118                 nlDeleteContext(sys->context);
119         }
120         sys->vertexCos = NULL;
121         sys->mfaces = NULL;
122         sys->medges = NULL;
123         MEM_freeN(sys);
124 }
125
126 static void memset_laplacian_system(LaplacianSystem *sys, int val)
127 {
128         memset(sys->eweights,     val, sizeof(float) * sys->numEdges);
129         memset(sys->fweights,     val, sizeof(float) * sys->numFaces * 3);
130         memset(sys->numNeEd,      val, sizeof(short) * sys->numVerts);
131         memset(sys->numNeFa,      val, sizeof(short) * sys->numVerts);
132         memset(sys->ring_areas,   val, sizeof(float) * sys->numVerts);
133         memset(sys->vlengths,     val, sizeof(float) * sys->numVerts);
134         memset(sys->vweights,     val, sizeof(float) * sys->numVerts);
135         memset(sys->zerola,       val, sizeof(short) * sys->numVerts);
136 }
137
138 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts)
139 {
140         LaplacianSystem *sys;
141         sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
142         sys->numEdges = a_numEdges;
143         sys->numFaces = a_numFaces;
144         sys->numVerts = a_numVerts;
145
146         sys->eweights =  MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
147         if (!sys->eweights) {
148                 delete_laplacian_system(sys);
149                 return NULL;
150         }
151
152         sys->fweights =  MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
153         if (!sys->fweights) {
154                 delete_laplacian_system(sys);
155                 return NULL;
156         }
157
158         sys->numNeEd =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothNumNeEd");
159         if (!sys->numNeEd) {
160                 delete_laplacian_system(sys);
161                 return NULL;
162         }
163
164         sys->numNeFa =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothNumNeFa");
165         if (!sys->numNeFa) {
166                 delete_laplacian_system(sys);
167                 return NULL;
168         }
169
170         sys->ring_areas =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
171         if (!sys->ring_areas) {
172                 delete_laplacian_system(sys);
173                 return NULL;
174         }
175
176         sys->vlengths =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
177         if (!sys->vlengths) {
178                 delete_laplacian_system(sys);
179                 return NULL;
180         }
181
182         sys->vweights =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
183         if (!sys->vweights) {
184                 delete_laplacian_system(sys);
185                 return NULL;
186         }
187
188         sys->zerola =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
189         if (!sys->zerola) {
190                 delete_laplacian_system(sys);
191                 return NULL;
192         }
193
194         return sys;
195 }
196
197 static void init_data(ModifierData *md)
198 {
199         LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
200         smd->lambda = 0.01f;
201         smd->lambda_border = 0.01f;
202         smd->repeat = 1;
203         smd->flag = MOD_LAPLACIANSMOOTH_X | MOD_LAPLACIANSMOOTH_Y | MOD_LAPLACIANSMOOTH_Z | MOD_LAPLACIANSMOOTH_PRESERVE_VOLUME | MOD_LAPLACIANSMOOTH_NORMALIZED;
204         smd->defgrp_name[0] = '\0';
205 }
206
207 static void copy_data(ModifierData *md, ModifierData *target)
208 {
209         LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
210         LaplacianSmoothModifierData *tsmd = (LaplacianSmoothModifierData *) target;
211
212         tsmd->lambda = smd->lambda;
213         tsmd->lambda_border = smd->lambda_border;
214         tsmd->repeat = smd->repeat;
215         tsmd->flag = smd->flag;
216         BLI_strncpy(tsmd->defgrp_name, smd->defgrp_name, sizeof(tsmd->defgrp_name));
217 }
218
219 static int is_disabled(ModifierData *md, int UNUSED(useRenderParams))
220 {
221         LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
222         short flag;
223
224         flag = smd->flag & (MOD_LAPLACIANSMOOTH_X | MOD_LAPLACIANSMOOTH_Y | MOD_LAPLACIANSMOOTH_Z);
225
226         /* disable if modifier is off for X, Y and Z or if factor is 0 */
227         if (flag == 0) return 1;
228
229         return 0;
230 }
231
232 static CustomDataMask required_data_mask(Object *UNUSED(ob), ModifierData *md)
233 {
234         LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *)md;
235         CustomDataMask dataMask = 0;
236
237         /* ask for vertexgroups if we need them */
238         if (smd->defgrp_name[0]) dataMask |= CD_MASK_MDEFORMVERT;
239
240         return dataMask;
241 }
242
243 static float average_area_quad_v3(float *v1, float *v2, float *v3, float *v4)
244 {
245         float areaq = 0.0f;
246         areaq = area_tri_v3(v1, v2, v3) + area_tri_v3(v1, v2, v4) + area_tri_v3(v1, v3, v4);
247         return areaq / 2.0f;
248 }
249
250 static float cotan_weight(float *v1, float *v2, float *v3)
251 {
252         float a[3], b[3], c[3], clen;
253
254         sub_v3_v3v3(a, v2, v1);
255         sub_v3_v3v3(b, v3, v1);
256         cross_v3_v3v3(c, a, b);
257
258         clen = len_v3(c);
259
260         if (clen == 0.0f)
261                 return 0.0f;
262
263         return dot_v3v3(a, b) / clen;
264 }
265
266 static float compute_volume(float (*vertexCos)[3], MFace *mfaces, int numFaces)
267 {
268         float vol = 0.0f;
269         float x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
270         int i;
271         float *vf[4];
272         for (i = 0; i < numFaces; i++) {
273                 vf[0] = vertexCos[mfaces[i].v1];
274                 vf[1] = vertexCos[mfaces[i].v2];
275                 vf[2] = vertexCos[mfaces[i].v3];
276
277                 x1 = vf[0][0];
278                 y1 = vf[0][1];
279                 z1 = vf[0][2];
280
281                 x2 = vf[1][0];
282                 y2 = vf[1][1];
283                 z2 = vf[1][2];
284
285                 x3 = vf[2][0];
286                 y3 = vf[2][1];
287                 z3 = vf[2][2];
288
289
290                 vol +=  (1.0f / 6.0f) * (x2 * y3 * z1 + x3 * y1 * z2 - x1 * y3 * z2 - x2 * y1 * z3 + x1 * y2 * z3 - x3 * y2 * z1);
291                 if ((&mfaces[i])->v4) {
292                         vf[3] = vertexCos[mfaces[i].v4];
293                         x4 = vf[3][0];
294                         y4 = vf[3][1];
295                         z4 = vf[3][2];
296                         vol += (1.0f / 6.0f) * (x1 * y3 * z4 - x1 * y4 * z3 - x3 * y1 * z4 + x3 * z1 * y4 + y1 * x4 * z3 - x4 * y3 * z1);
297                 }
298         }
299         return fabsf(vol);
300 }
301
302 static void volume_preservation(LaplacianSystem *sys, float vini, float vend, short flag)
303 {
304         float beta;
305         int i;
306
307         if (vend != 0.0f) {
308                 beta  = pow(vini / vend, 1.0f / 3.0f);
309                 for (i = 0; i < sys->numVerts; i++) {
310                         if (flag & MOD_LAPLACIANSMOOTH_X) {
311                                 sys->vertexCos[i][0] = (sys->vertexCos[i][0] - sys->vert_centroid[0]) * beta + sys->vert_centroid[0];
312                         }
313                         if (flag & MOD_LAPLACIANSMOOTH_Y) {
314                                 sys->vertexCos[i][1] = (sys->vertexCos[i][1] - sys->vert_centroid[1]) * beta + sys->vert_centroid[1];
315                         }
316                         if (flag & MOD_LAPLACIANSMOOTH_Z) {
317                                 sys->vertexCos[i][2] = (sys->vertexCos[i][2] - sys->vert_centroid[2]) * beta + sys->vert_centroid[2];
318                         }
319
320                 }
321         }
322 }
323
324 static void init_laplacian_matrix(LaplacianSystem *sys)
325 {
326         float *v1, *v2, *v3, *v4;
327         float w1, w2, w3, w4;
328         float areaf;
329         int i, j;
330         unsigned int idv1, idv2, idv3, idv4, idv[4];
331         int has_4_vert;
332         for (i = 0; i < sys->numEdges; i++) {
333                 idv1 = sys->medges[i].v1;
334                 idv2 = sys->medges[i].v2;
335
336                 v1 = sys->vertexCos[idv1];
337                 v2 = sys->vertexCos[idv2];
338
339                 sys->numNeEd[idv1] = sys->numNeEd[idv1] + 1;
340                 sys->numNeEd[idv2] = sys->numNeEd[idv2] + 1;
341                 w1 = len_v3v3(v1, v2);
342                 if (w1 < sys->min_area) {
343                         sys->zerola[idv1] = 1;
344                         sys->zerola[idv2] = 1;
345                 }
346                 else {
347                         w1 = 1.0f / w1;
348                 }
349
350                 sys->eweights[i] = w1;
351         }
352         for (i = 0; i < sys->numFaces; i++) {
353                 has_4_vert = ((&sys->mfaces[i])->v4) ? 1 : 0;
354
355                 idv1 = sys->mfaces[i].v1;
356                 idv2 = sys->mfaces[i].v2;
357                 idv3 = sys->mfaces[i].v3;
358                 idv4 = has_4_vert ? sys->mfaces[i].v4 : 0;
359
360                 sys->numNeFa[idv1] += 1;
361                 sys->numNeFa[idv2] += 1;
362                 sys->numNeFa[idv3] += 1;
363                 if (has_4_vert) sys->numNeFa[idv4] += 1;
364
365                 v1 = sys->vertexCos[idv1];
366                 v2 = sys->vertexCos[idv2];
367                 v3 = sys->vertexCos[idv3];
368                 v4 = has_4_vert ? sys->vertexCos[idv4] : 0;
369
370                 if (has_4_vert) {
371                         areaf = area_quad_v3(v1, v2, v3, sys->vertexCos[sys->mfaces[i].v4]);
372                 }
373                 else {
374                         areaf = area_tri_v3(v1, v2, v3);
375                 }
376                 if (fabsf(areaf) < sys->min_area) {
377                         sys->zerola[idv1] = 1;
378                         sys->zerola[idv2] = 1;
379                         sys->zerola[idv3] = 1;
380                         if (has_4_vert) sys->zerola[idv4] = 1;
381                 }
382
383                 if (has_4_vert) {
384                         sys->ring_areas[idv1] += average_area_quad_v3(v1, v2, v3, v4);
385                         sys->ring_areas[idv2] += average_area_quad_v3(v2, v3, v4, v1);
386                         sys->ring_areas[idv3] += average_area_quad_v3(v3, v4, v1, v2);
387                         sys->ring_areas[idv4] += average_area_quad_v3(v4, v1, v2, v3);
388                 }
389                 else {
390                         sys->ring_areas[idv1] += areaf;
391                         sys->ring_areas[idv2] += areaf;
392                         sys->ring_areas[idv3] += areaf;
393                 }
394
395                 if (has_4_vert) {
396
397                         idv[0] = idv1;
398                         idv[1] = idv2;
399                         idv[2] = idv3;
400                         idv[3] = idv4;
401
402                         for (j = 0; j < 4; j++) {
403                                 idv1 = idv[j];
404                                 idv2 = idv[(j + 1) % 4];
405                                 idv3 = idv[(j + 2) % 4];
406                                 idv4 = idv[(j + 3) % 4];
407
408                                 v1 = sys->vertexCos[idv1];
409                                 v2 = sys->vertexCos[idv2];
410                                 v3 = sys->vertexCos[idv3];
411                                 v4 = sys->vertexCos[idv4];
412
413                                 w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
414                                 w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
415                                 w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
416
417                                 sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
418                         }
419                 }
420                 else {
421                         w1 = cotan_weight(v1, v2, v3) / 2.0f;
422                         w2 = cotan_weight(v2, v3, v1) / 2.0f;
423                         w3 = cotan_weight(v3, v1, v2) / 2.0f;
424
425                         sys->fweights[i][0] = sys->fweights[i][0] + w1;
426                         sys->fweights[i][1] = sys->fweights[i][1] + w2;
427                         sys->fweights[i][2] = sys->fweights[i][2] + w3;
428
429                         sys->vweights[idv1] = sys->vweights[idv1] + w2 + w3;
430                         sys->vweights[idv2] = sys->vweights[idv2] + w1 + w3;
431                         sys->vweights[idv3] = sys->vweights[idv3] + w1 + w2;
432                 }
433         }
434         for (i = 0; i < sys->numEdges; i++) {
435                 idv1 = sys->medges[i].v1;
436                 idv2 = sys->medges[i].v2;
437                 /* if is boundary, apply scale-dependent umbrella operator only with neighboors in boundary */
438                 if (sys->numNeEd[idv1] != sys->numNeFa[idv1] && sys->numNeEd[idv2] != sys->numNeFa[idv2]) {
439                         sys->vlengths[idv1] += sys->eweights[i];
440                         sys->vlengths[idv2] += sys->eweights[i];
441                 }
442         }
443
444 }
445
446 static void fill_laplacian_matrix(LaplacianSystem *sys)
447 {
448         float *v1, *v2, *v3, *v4;
449         float w2, w3, w4;
450         int i, j;
451         int has_4_vert;
452         unsigned int idv1, idv2, idv3, idv4, idv[4];
453
454         for (i = 0; i < sys->numFaces; i++) {
455                 idv1 = sys->mfaces[i].v1;
456                 idv2 = sys->mfaces[i].v2;
457                 idv3 = sys->mfaces[i].v3;
458                 has_4_vert = ((&sys->mfaces[i])->v4) ? 1 : 0;
459
460                 if (has_4_vert) {
461                         idv[0] = sys->mfaces[i].v1;
462                         idv[1] = sys->mfaces[i].v2;
463                         idv[2] = sys->mfaces[i].v3;
464                         idv[3] = sys->mfaces[i].v4;
465                         for (j = 0; j < 4; j++) {
466                                 idv1 = idv[j];
467                                 idv2 = idv[(j + 1) % 4];
468                                 idv3 = idv[(j + 2) % 4];
469                                 idv4 = idv[(j + 3) % 4];
470
471                                 v1 = sys->vertexCos[idv1];
472                                 v2 = sys->vertexCos[idv2];
473                                 v3 = sys->vertexCos[idv3];
474                                 v4 = sys->vertexCos[idv4];
475
476                                 w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
477                                 w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
478                                 w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
479
480                                 w2 = w2 / 4.0f;
481                                 w3 = w3 / 4.0f;
482                                 w4 = w4 / 4.0f;
483
484                                 if (sys->numNeEd[idv1] == sys->numNeFa[idv1] && sys->zerola[idv1] == 0) {
485                                         nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
486                                         nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
487                                         nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
488                                 }
489                         }
490                 }
491                 else {
492                         /* Is ring if number of faces == number of edges around vertice*/
493                         if (sys->numNeEd[idv1] == sys->numNeFa[idv1] && sys->zerola[idv1] == 0) {
494                                 nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
495                                 nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
496                         }
497                         if (sys->numNeEd[idv2] == sys->numNeFa[idv2] && sys->zerola[idv2] == 0) {
498                                 nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
499                                 nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
500                         }
501                         if (sys->numNeEd[idv3] == sys->numNeFa[idv3] && sys->zerola[idv3] == 0) {
502                                 nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
503                                 nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
504                         }
505                 }
506         }
507
508         for (i = 0; i < sys->numEdges; i++) {
509                 idv1 = sys->medges[i].v1;
510                 idv2 = sys->medges[i].v2;
511                 /* Is boundary */
512                 if (sys->numNeEd[idv1] != sys->numNeFa[idv1] &&
513                     sys->numNeEd[idv2] != sys->numNeFa[idv2] &&
514                     sys->zerola[idv1] == 0 &&
515                     sys->zerola[idv2] == 0)
516                 {
517                         nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
518                         nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
519                 }
520         }
521 }
522
523 static void validate_solution(LaplacianSystem *sys, short flag, float lambda, float lambda_border)
524 {
525         int i;
526         float lam;
527         float vini, vend;
528
529         if (flag & MOD_LAPLACIANSMOOTH_PRESERVE_VOLUME) {
530                 vini = compute_volume(sys->vertexCos, sys->mfaces, sys->numFaces);
531         }
532         for (i = 0; i < sys->numVerts; i++) {
533                 if (sys->zerola[i] == 0) {
534                         lam = sys->numNeEd[i] == sys->numNeFa[i] ? (lambda >= 0.0f ? 1.0f : -1.0f) : (lambda_border >= 0.0f ? 1.0f : -1.0f);
535                         if (flag & MOD_LAPLACIANSMOOTH_X) {
536                                 sys->vertexCos[i][0] += lam * (nlGetVariable(0, i) - sys->vertexCos[i][0]);
537                         }
538                         if (flag & MOD_LAPLACIANSMOOTH_Y) {
539                                 sys->vertexCos[i][1] += lam * (nlGetVariable(1, i) - sys->vertexCos[i][1]);
540                         }
541                         if (flag & MOD_LAPLACIANSMOOTH_Z) {
542                                 sys->vertexCos[i][2] += lam * (nlGetVariable(2, i) - sys->vertexCos[i][2]);
543                         }
544                 }
545         }
546         if (flag & MOD_LAPLACIANSMOOTH_PRESERVE_VOLUME) {
547                 vend = compute_volume(sys->vertexCos, sys->mfaces, sys->numFaces);
548                 volume_preservation(sys, vini, vend, flag);
549         }
550 }
551
552 static void laplaciansmoothModifier_do(
553         LaplacianSmoothModifierData *smd, Object *ob, DerivedMesh *dm,
554         float (*vertexCos)[3], int numVerts)
555 {
556         LaplacianSystem *sys;
557         MDeformVert *dvert = NULL;
558         MDeformVert *dv = NULL;
559         float w, wpaint;
560         int i, iter;
561         int defgrp_index;
562
563         DM_ensure_tessface(dm);
564
565         sys = init_laplacian_system(dm->getNumEdges(dm), dm->getNumTessFaces(dm), numVerts);
566         if (!sys) {
567                 return;
568         }
569
570         sys->mfaces = dm->getTessFaceArray(dm);
571         sys->medges = dm->getEdgeArray(dm);
572         sys->vertexCos = vertexCos;
573         sys->min_area = 0.00001f;
574         modifier_get_vgroup(ob, dm, smd->defgrp_name, &dvert, &defgrp_index);
575
576         sys->vert_centroid[0] = 0.0f;
577         sys->vert_centroid[1] = 0.0f;
578         sys->vert_centroid[2] = 0.0f;
579         memset_laplacian_system(sys, 0);
580         nlNewContext();
581         sys->context = nlGetCurrent();
582         nlSolverParameteri(NL_NB_VARIABLES, numVerts);
583         nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
584         nlSolverParameteri(NL_NB_ROWS, numVerts);
585         nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
586
587         init_laplacian_matrix(sys);
588
589         for (iter = 0; iter < smd->repeat; iter++) {
590                 nlBegin(NL_SYSTEM);
591                 for (i = 0; i < numVerts; i++) {
592                         nlSetVariable(0, i, vertexCos[i][0]);
593                         nlSetVariable(1, i, vertexCos[i][1]);
594                         nlSetVariable(2, i, vertexCos[i][2]);
595                         if (iter == 0) {
596                                 add_v3_v3(sys->vert_centroid, vertexCos[i]);
597                         }
598                 }
599                 if (iter == 0 && numVerts > 0) {
600                         mul_v3_fl(sys->vert_centroid, 1.0f / (float)numVerts);
601                 }
602
603                 nlBegin(NL_MATRIX);
604                 dv = dvert;
605                 for (i = 0; i < numVerts; i++) {
606                         nlRightHandSideSet(0, i, vertexCos[i][0]);
607                         nlRightHandSideSet(1, i, vertexCos[i][1]);
608                         nlRightHandSideSet(2, i, vertexCos[i][2]);
609                         if (iter == 0) {
610                                 if (dv) {
611                                         wpaint = defvert_find_weight(dv, defgrp_index);
612                                         dv++;
613                                 }
614                                 else {
615                                         wpaint = 1.0f;
616                                 }
617
618                                 if (sys->zerola[i] == 0) {
619                                         if (smd->flag & MOD_LAPLACIANSMOOTH_NORMALIZED) {
620                                                 w = sys->vweights[i];
621                                                 sys->vweights[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda) * wpaint / w;
622                                                 w = sys->vlengths[i];
623                                                 sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
624                                                 if (sys->numNeEd[i] == sys->numNeFa[i]) {
625                                                         nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda) * wpaint);
626                                                 }
627                                                 else {
628                                                         nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
629                                                 }
630                                         }
631                                         else {
632                                                 w = sys->vweights[i] * sys->ring_areas[i];
633                                                 sys->vweights[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda) * wpaint / (4.0f * w);
634                                                 w = sys->vlengths[i];
635                                                 sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
636
637                                                 if (sys->numNeEd[i] == sys->numNeFa[i]) {
638                                                         nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda) * wpaint / (4.0f * sys->ring_areas[i]));
639                                                 }
640                                                 else {
641                                                         nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
642                                                 }
643                                         }
644                                 }
645                                 else {
646                                         nlMatrixAdd(i, i, 1.0f);
647                                 }
648                         }
649                 }
650
651                 if (iter == 0) {
652                         fill_laplacian_matrix(sys);
653                 }
654
655                 nlEnd(NL_MATRIX);
656                 nlEnd(NL_SYSTEM);
657
658                 if (nlSolveAdvanced(NULL, NL_TRUE)) {
659                         validate_solution(sys, smd->flag, smd->lambda, smd->lambda_border);
660                 }
661         }
662         nlDeleteContext(sys->context);
663         sys->context = NULL;
664         delete_laplacian_system(sys);
665
666 }
667
668 static void deformVerts(ModifierData *md, Object *ob, DerivedMesh *derivedData,
669                         float (*vertexCos)[3], int numVerts, ModifierApplyFlag UNUSED(flag))
670 {
671         DerivedMesh *dm = get_dm(ob, NULL, derivedData, NULL, 0);
672
673         laplaciansmoothModifier_do((LaplacianSmoothModifierData *)md, ob, dm,
674                                    vertexCos, numVerts);
675
676         if (dm != derivedData)
677                 dm->release(dm);
678 }
679
680 static void deformVertsEM(
681         ModifierData *md, Object *ob, struct BMEditMesh *editData,
682         DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
683 {
684         DerivedMesh *dm = get_dm(ob, editData, derivedData, NULL, 0);
685
686         laplaciansmoothModifier_do((LaplacianSmoothModifierData *)md, ob, dm,
687                                    vertexCos, numVerts);
688
689         if (dm != derivedData)
690                 dm->release(dm);
691 }
692
693
694 ModifierTypeInfo modifierType_LaplacianSmooth = {
695         /* name */              "Laplacian Smooth",
696         /* structName */        "LaplacianSmoothModifierData",
697         /* structSize */        sizeof(LaplacianSmoothModifierData),
698         /* type */              eModifierTypeType_OnlyDeform,
699         /* flags */             eModifierTypeFlag_AcceptsMesh |
700                                 eModifierTypeFlag_SupportsEditmode,
701
702         /* copy_data */         copy_data,
703         /* deformVerts */       deformVerts,
704         /* deformMatrices */    NULL,
705         /* deformVertsEM */     deformVertsEM,
706         /* deformMatricesEM */  NULL,
707         /* applyModifier */     NULL,
708         /* applyModifierEM */   NULL,
709         /* initData */          init_data,
710         /* requiredDataMask */  required_data_mask,
711         /* freeData */          NULL,
712         /* isDisabled */        is_disabled,
713         /* updateDepgraph */    NULL,
714         /* dependsOnTime */     NULL,
715         /* dependsOnNormals */  NULL,
716         /* foreachObjectLink */ NULL,
717         /* foreachIDLink */     NULL,
718         /* foreachTexLink */    NULL,
719 };