Dynamic Paint:
[blender.git] / source / blender / blenkernel / intern / dynamicpaint.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  * Contributor(s): Miika Hämäläinen
10  *
11  * ***** END GPL LICENSE BLOCK *****
12  */
13
14
15 #include "MEM_guardedalloc.h"
16
17 #include <math.h>
18 #include <stdio.h>
19
20 #include "BLI_blenlib.h"
21 #include "BLI_math.h"
22 #include "BLI_kdtree.h"
23 #include "BLI_utildefines.h"
24
25 /* Platform independend time    */
26 #include "PIL_time.h"
27
28 #include "DNA_anim_types.h"
29 #include "DNA_dynamicpaint_types.h"
30 #include "DNA_group_types.h" /*GroupObject*/
31 #include "DNA_mesh_types.h"
32 #include "DNA_meshdata_types.h"
33 #include "DNA_modifier_types.h"
34 #include "DNA_object_types.h"
35 #include "DNA_scene_types.h"
36 #include "DNA_userdef_types.h"  /* to get temp file path        */
37
38 #include "BKE_animsys.h"
39 #include "BKE_bvhutils.h"       /* bvh tree     */
40 #include "BKE_blender.h"
41 #include "BKE_cdderivedmesh.h"
42 #include "BKE_context.h"
43 #include "BKE_customdata.h"
44 #include "BKE_colortools.h"
45 #include "BKE_deform.h"
46 #include "BKE_depsgraph.h"
47 #include "BKE_DerivedMesh.h"
48 #include "BKE_dynamicpaint.h"
49 #include "BKE_effect.h"
50 #include "BKE_global.h"
51 #include "BKE_main.h"
52 #include "BKE_material.h"
53 #include "BKE_modifier.h"
54 #include "BKE_object.h"
55 #include "BKE_particle.h"
56 #include "BKE_pointcache.h"
57 #include "BKE_report.h"
58 #include "BKE_scene.h"
59 #include "BKE_texture.h"
60
61 #include "RNA_access.h"
62 #include "RNA_define.h"
63 #include "RNA_enum_types.h"
64
65 #include "ED_screen.h"
66 #include "WM_api.h"
67
68 /* for image output     */
69 #include "IMB_imbuf_types.h"
70 #include "IMB_imbuf.h"
71 #include "BKE_image.h"
72 #include "intern/IMB_filetype.h"
73 #ifdef WITH_OPENEXR
74 #include "intern/openexr/openexr_api.h"
75 #endif
76
77 /* uv validate  */
78 #include "intern/MOD_util.h"
79
80 /* to read object material color        */
81 #include "DNA_texture_types.h"
82 #include "../render/intern/include/pointdensity.h"
83 #include "../render/intern/include/render_types.h"
84 #include "../render/intern/include/voxeldata.h"
85 #include "../render/intern/include/texture.h"
86 #include "DNA_material_types.h"
87 #include "RE_render_ext.h"
88
89 #ifdef _OPENMP
90 #include <omp.h>
91 #endif
92
93 #define DPOUTPUT_JPEG 0
94 #define DPOUTPUT_PNG 1
95 #define DPOUTPUT_OPENEXR 2
96
97 struct Object;
98 struct Scene;
99 struct DerivedMesh;
100
101 /* precalculated gaussian factors for 5x super sampling */
102 float gaussianFactors[5] = {    0.996849f,
103                                                                 0.596145f,
104                                                                 0.596145f,
105                                                                 0.596145f,
106                                                                 0.524141f};
107 float gaussianTotal = 3.309425f;
108
109 /*
110 *       UV Image neighbouring pixel table x and y list
111 */
112 int neighX[8] = {1,1,0,-1,-1,-1, 0, 1};
113 int neighY[8] = {0,1,1, 1, 0,-1,-1,-1};
114
115 static int dynamicPaint_doStep(Scene *scene, Object *ob, DynamicPaintSurface *surface, float timescale, float subframe);
116 static int dynamicPaint_calculateFrame(DynamicPaintSurface *surface, Scene *scene, Object *cObject, int frame);
117
118 /***************************** Internal Structs ***************************/
119
120 typedef struct Bounds2D {
121         float min[2], max[2];
122 } Bounds2D;
123
124 typedef struct Bounds3D {
125         int valid;
126         float min[3], max[3];
127 } Bounds3D;
128
129 typedef struct VolumeGrid {
130         int x,y,z;
131         Bounds3D grid_bounds; /* whole grid bounds */
132
133         Bounds3D *bounds;       /* (x*y*z) precalculated grid cell bounds */
134         unsigned int *s_pos; /* (x*y*z) t_index begin id */
135         unsigned int *s_num; /* (x*y*z) number of t_index points */
136         unsigned int *t_index; /* actual surface point index,
137                                                    access: (s_pos+s_num) */
138 } VolumeGrid;
139
140 typedef struct Vec3f {
141         float v[3];
142 } Vec3f;
143
144 typedef struct BakeNeighPoint {
145         float dir[3];   /* vector pointing towards this neighbour */
146         float dist;             /* distance to */
147 } BakeNeighPoint;
148
149 /* Surface data used while processing a frame   */
150 typedef struct PaintBakeNormal {
151         float invNorm[3];  /* current pixel world-space inverted normal */
152         float normal_scale; /* normal directional scale for displace mapping */
153 } PaintBakeNormal;
154
155 /* Temp surface data used to process a frame */
156 typedef struct PaintBakeData {
157         /* point space data */
158         PaintBakeNormal *bNormal;
159         unsigned int *s_pos;    /* index to start reading point sample realCoord */
160         unsigned int *s_num;    /* num of realCoord samples */
161         Vec3f *realCoord;  /* current pixel center world-space coordinates for each sample
162                                            *  ordered as (s_pos+s_num)*/
163
164         /* adjacency info */
165         BakeNeighPoint *bNeighs; /* current global neighbour distances and directions, if required */
166         double average_dist;
167         /* space partitioning */
168         VolumeGrid *grid;               /* space partitioning grid to optimize brush checks */
169
170         /* velocity and movement */
171         Vec3f *velocity;                /* speed vector in global space movement per frame, if required */
172         Vec3f *prev_velocity;
173         float *brush_velocity;  /* special temp data for post-p velocity based brushes like smudge
174                                                         *  3 float dir vec + 1 float str */
175         MVert *prev_verts;              /* copy of previous frame vertices. used to observe surface movement */
176         float prev_obmat[4][4]; /* previous frame object matrix */
177         int clear;                              /* flag to check if surface was cleared/reset -> have to redo velocity etc. */
178
179 } PaintBakeData;
180
181 /* UV Image sequence format point       */
182 typedef struct PaintUVPoint {
183         /* Pixel / mesh data */
184         unsigned int face_index, pixel_index;   /* face index on domain derived mesh */
185         unsigned int v1, v2, v3;                                /* vertex indexes */
186
187         unsigned int neighbour_pixel;   /* If this pixel isn't uv mapped to any face,
188                                                                            but it's neighbouring pixel is */
189         short quad;
190 } PaintUVPoint;
191
192 typedef struct ImgSeqFormatData {
193         PaintUVPoint *uv_p;
194         Vec3f *barycentricWeights;              /* b-weights for all pixel samples */
195 } ImgSeqFormatData;
196
197 typedef struct EffVelPoint {
198         float previous_pos[3];
199         float previous_vel[3];
200 } EffVelPoint;
201
202
203 /* adjacency data flags */
204 #define ADJ_ON_MESH_EDGE (1<<0)
205
206 typedef struct PaintAdjData {
207         unsigned int *n_target;         /* array of neighbouring point indexes,
208                                                                for single sample use (n_index+neigh_num) */
209         unsigned int *n_index;          /* index to start reading n_target for each point */
210         unsigned int *n_num;            /* num of neighs for each point */
211         unsigned int *flags;            /* vertex adjacency flags */
212         unsigned int total_targets; /* size of n_target */
213 } PaintAdjData;
214
215 /***************************** General Utils ******************************/
216
217 /*
218 *       Output error message to both ui and console
219 */
220 static int printError(DynamicPaintCanvasSettings *canvas, char *string)
221 {
222         if (strlen(string)>64) string[63] = '\0';
223
224         /* Add error to canvas ui info label */
225         sprintf(canvas->error, string);
226
227         /* Print console output */
228         printf("DynamicPaint bake failed: %s\n", canvas->error);
229
230         return 0;
231 }
232
233 /* Get number of surface points for cached types */
234 static int dynamicPaint_surfaceNumOfPoints(DynamicPaintSurface *surface)
235 {
236         if (surface->format == MOD_DPAINT_SURFACE_F_PTEX) {
237                 return 0; /* not supported atm */
238         }
239         else if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
240                 if (!surface->canvas->dm) return 0; /* invalid derived mesh */
241                 return surface->canvas->dm->getNumVerts(surface->canvas->dm);
242         }
243         else
244                 return 0;
245 }
246
247 /* checks whether surface's format/type has realtime preview */
248 int dynamicPaint_surfaceHasColorPreview(DynamicPaintSurface *surface)
249 {
250         if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) return 0;
251         else if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
252                 if (surface->type == MOD_DPAINT_SURFACE_T_DISPLACE ||
253                         surface->type == MOD_DPAINT_SURFACE_T_WAVE) return 0;
254                 else return 1;
255         }
256         else return 1;
257 }
258
259 /* get currently active surface (in user interface) */
260 struct DynamicPaintSurface *get_activeSurface(DynamicPaintCanvasSettings *canvas)
261 {
262         DynamicPaintSurface *surface = canvas->surfaces.first;
263         int i;
264
265         for(i=0; surface; surface=surface->next) {
266                 if(i == canvas->active_sur)
267                         return surface;
268                 i++;
269         }
270         return NULL;
271 }
272
273 /* set preview to first previewable surface */
274 void dynamicPaint_resetPreview(DynamicPaintCanvasSettings *canvas)
275 {
276         DynamicPaintSurface *surface = canvas->surfaces.first;
277         int done=0;
278
279         for(; surface; surface=surface->next) {
280                 if (!done && dynamicPaint_surfaceHasColorPreview(surface)) {
281                         surface->flags |= MOD_DPAINT_PREVIEW;
282                         done=1;
283                 }
284                 else
285                         surface->flags &= ~MOD_DPAINT_PREVIEW;
286         }
287 }
288
289 /* set preview to defined surface */
290 static void dynamicPaint_setPreview(DynamicPaintSurface *t_surface)
291 {
292         DynamicPaintSurface *surface = t_surface->canvas->surfaces.first;
293         for(; surface; surface=surface->next) {
294                 if (surface == t_surface)
295                         surface->flags |= MOD_DPAINT_PREVIEW;
296                 else
297                         surface->flags &= ~MOD_DPAINT_PREVIEW;
298         }
299 }
300
301 int dynamicPaint_outputLayerExists(struct DynamicPaintSurface *surface, Object *ob, int index)
302 {
303         char *name;
304
305         if (index == 0)
306                 name = surface->output_name;
307         else if (index == 1)
308                 name = surface->output_name2;
309         else
310                 return 0;
311
312         if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
313                 if (surface->type == MOD_DPAINT_SURFACE_T_PAINT) {
314                         Mesh *me = ob->data;
315                         return (CustomData_get_named_layer_index(&me->fdata, CD_MCOL, name) != -1);
316                 }
317                 else if (surface->type == MOD_DPAINT_SURFACE_T_WEIGHT)
318                         return (defgroup_name_index(ob, surface->output_name) != -1);
319         }
320
321         return 0;
322 }
323
324 static int surface_duplicateOutputExists(void *arg, const char *name)
325 {
326         DynamicPaintSurface *t_surface = (DynamicPaintSurface*)arg;
327         DynamicPaintSurface *surface = t_surface->canvas->surfaces.first;
328
329         for(; surface; surface=surface->next) {
330                 if (surface!=t_surface && surface->type==t_surface->type &&
331                         surface->format==t_surface->format) {
332                         if (surface->output_name[0]!='\0' && !strcmp(name, surface->output_name)) return 1;
333                         if (surface->output_name2[0]!='\0' && !strcmp(name, surface->output_name2)) return 1;
334                 }
335         }
336         return 0;
337 }
338
339 void surface_setUniqueOutputName(DynamicPaintSurface *surface, char *basename, int output)
340 {
341         char name[64];
342         strncpy(name, basename, 62); /* in case basename is surface->name use a copy */
343         if (!output)
344                 BLI_uniquename_cb(surface_duplicateOutputExists, surface, name, '.', surface->output_name, sizeof(surface->output_name));
345         if (output)
346                 BLI_uniquename_cb(surface_duplicateOutputExists, surface, name, '.', surface->output_name2, sizeof(surface->output_name2));
347 }
348
349
350 static int surface_duplicateNameExists(void *arg, const char *name)
351 {
352         DynamicPaintSurface *t_surface = (DynamicPaintSurface*)arg;
353         DynamicPaintSurface *surface = t_surface->canvas->surfaces.first;
354
355         for(; surface; surface=surface->next) {
356                 if (surface!=t_surface && !strcmp(name, surface->name)) return 1;
357         }
358         return 0;
359 }
360
361 void dynamicPaintSurface_setUniqueName(DynamicPaintSurface *surface, char *basename)
362 {
363         char name[64];
364         strncpy(name, basename, 62); /* in case basename is surface->name use a copy */
365         BLI_uniquename_cb(surface_duplicateNameExists, surface, name, '.', surface->name, sizeof(surface->name));
366 }
367
368
369 /* change surface data to defaults on new type */
370 void dynamicPaintSurface_updateType(struct DynamicPaintSurface *surface)
371 {
372         if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) {
373                 surface->output_name[0]='\0';
374                 surface->output_name2[0]='\0';
375                 surface->flags |= MOD_DPAINT_ANTIALIAS;
376                 surface->disp_clamp = 1.0f;
377         }
378         else {
379                 sprintf(surface->output_name, "dp_");
380                 strcpy(surface->output_name2,surface->output_name);
381                 surface->flags &= ~MOD_DPAINT_ANTIALIAS;
382                 surface->disp_clamp = 0.0f;
383         }
384
385         if (surface->type == MOD_DPAINT_SURFACE_T_PAINT) {
386                 strcat(surface->output_name,"paintmap");
387                 strcat(surface->output_name2,"wetmap");
388                 surface_setUniqueOutputName(surface, surface->output_name2, 1);
389         }
390         else if (surface->type == MOD_DPAINT_SURFACE_T_DISPLACE) {
391                 strcat(surface->output_name,"displace");
392         }
393         else if (surface->type == MOD_DPAINT_SURFACE_T_WEIGHT) {
394                 strcat(surface->output_name,"weight");
395         }
396         else if (surface->type == MOD_DPAINT_SURFACE_T_WAVE) {
397                 strcat(surface->output_name,"wave");
398         }
399
400         surface_setUniqueOutputName(surface, surface->output_name, 0);
401
402         /* update preview */
403         if (dynamicPaint_surfaceHasColorPreview(surface))
404                 dynamicPaint_setPreview(surface);
405         else
406                 dynamicPaint_resetPreview(surface->canvas);
407 }
408
409 static int surface_totalSamples(DynamicPaintSurface *surface)
410 {
411         if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ &&
412                 surface->flags & MOD_DPAINT_ANTIALIAS)
413                 return (surface->data->total_points*5);
414         if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX &&
415                 surface->flags & MOD_DPAINT_ANTIALIAS && surface->data->adj_data)
416                 return (surface->data->total_points+surface->data->adj_data->total_targets);
417
418         return surface->data->total_points;
419 }
420
421 /* assumes source alpha > 0.0f or results NaN colors */
422 static void mixColors(float *t_color, float t_alpha, float *s_color, float s_alpha)
423 {
424         float invFact = (s_alpha<t_alpha) ? 1.0f : t_alpha/s_alpha;
425         float factor = 1.0f - invFact;
426
427         /* set initial color depending on existing alpha */
428         t_color[0] = t_color[0]*invFact + s_color[0]*factor;
429         t_color[1] = t_color[1]*invFact + s_color[1]*factor;
430         t_color[2] = t_color[2]*invFact + s_color[2]*factor;
431
432         /* mix final color */
433         factor = s_alpha;
434         invFact = 1.0f - factor;
435         t_color[0] = t_color[0]*invFact + s_color[0]*factor;
436         t_color[1] = t_color[1]*invFact + s_color[1]*factor;
437         t_color[2] = t_color[2]*invFact + s_color[2]*factor;
438 }
439
440 /* set "ignore cache" flag for all caches on this object */
441 static void object_cacheIgnoreClear(Object *ob, int state)
442 {
443         ListBase pidlist;
444         PTCacheID *pid;
445         BKE_ptcache_ids_from_object(&pidlist, ob, NULL, 0);
446
447         for(pid=pidlist.first; pid; pid=pid->next) {
448                 if(pid->cache) {
449                         if (state)
450                                 pid->cache->flag |= PTCACHE_IGNORE_CLEAR;
451                         else
452                                 pid->cache->flag &= ~PTCACHE_IGNORE_CLEAR;
453                 }
454         }
455
456         BLI_freelistN(&pidlist);
457 }
458
459 #define UPDATE_PARENTS (1<<0)
460 #define UPDATE_MESH (1<<1)
461 #define UPDATE_EVERYTHING (UPDATE_PARENTS|UPDATE_MESH)
462
463 static void subframe_updateObject(Scene *scene, Object *ob, int flags, float frame)
464 {
465         int oflags;
466         DynamicPaintModifierData *pmd = (DynamicPaintModifierData *)modifiers_findByType(ob, eModifierType_DynamicPaint);
467
468         /* if other is dynamic paint canvas, dont update */
469         if (pmd && pmd->canvas)
470                 return;
471
472         /* if object has parent, update it too */
473         if ((flags & UPDATE_PARENTS) && ob->parent) subframe_updateObject(scene, ob->parent, 0, frame);
474         if ((flags & UPDATE_PARENTS) && ob->track) subframe_updateObject(scene, ob->track, 0, frame);
475
476         /* for curve */
477         if(ob->type==OB_CURVE) {
478                 Curve *cu= ob->data;
479                 BKE_animsys_evaluate_animdata(scene, &cu->id, cu->adt, frame, ADT_RECALC_ANIM);
480         }
481         
482         /* backup object flags */
483         oflags = ob->recalc;
484
485         ob->recalc |= OB_RECALC_ALL;
486         ob->recalc |= OB_RECALC_DATA;
487         BKE_animsys_evaluate_animdata(scene, &ob->id, ob->adt, frame, ADT_RECALC_ANIM);
488         if (flags & UPDATE_MESH) {
489                 /* ignore cache clear during subframe updates
490                 *  to not mess up cache validity */
491                 object_cacheIgnoreClear(ob, 1);
492                 object_handle_update(scene, ob);
493                 object_cacheIgnoreClear(ob, 0);
494         }
495         else
496                 where_is_object_time(scene, ob, frame);
497
498         /* restore object flags */
499         ob->recalc = oflags;
500 }
501
502 static void scene_setSubframe(Scene *scene, float subframe)
503 {
504         /* dynamic paint subframes must be done on previous frame */
505         scene->r.cfra -= 1;
506         scene->r.subframe = subframe;
507 }
508
509 #define BRUSH_USES_VELOCITY (1<<0)
510
511 static int surface_getBrushFlags(DynamicPaintSurface *surface, Scene *scene, Object *ob)
512 {
513         Base *base = NULL;
514         GroupObject *go = NULL; 
515         Object *brushObj = NULL;
516         ModifierData *md = NULL;
517
518         int flags = 0;
519
520         if(surface->brush_group)
521                 go = surface->brush_group->gobject.first;
522         else
523                 base = scene->base.first;
524
525         while (base || go)
526         {
527                 brushObj = NULL;
528
529                 /* select object */
530                 if(surface->brush_group) {                                              
531                         if(go->ob)      brushObj = go->ob;                                      
532                 }                                       
533                 else                                            
534                         brushObj = base->object;
535
536                 if(!brushObj)                                   
537                 {
538                         if(surface->brush_group) go = go->next;
539                         else base= base->next;                                  
540                         continue;                       
541                 }
542
543                 if(surface->brush_group)
544                         go = go->next;
545                 else
546                         base= base->next;
547
548                 md = modifiers_findByType(brushObj, eModifierType_DynamicPaint);
549                 if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))                                    
550                 {
551                         DynamicPaintModifierData *pmd2 = (DynamicPaintModifierData *)md;
552
553                         if (pmd2->brush)
554                         {
555                                 DynamicPaintBrushSettings *brush = pmd2->brush;
556
557                                 if (brush->flags & MOD_DPAINT_USES_VELOCITY)
558                                         flags |= BRUSH_USES_VELOCITY;
559                         }
560                 }
561         }
562
563         return flags;
564 }
565
566 /* check whether two bounds intersect */
567 static int boundsIntersect(Bounds3D *b1, Bounds3D *b2)
568 {
569         int i=2;
570         if (!b1->valid || !b2->valid) return 0;
571         for (; i>=0; i-=1)
572                 if (!(b1->min[i] <= b2->max[i] && b1->max[i] >= b2->min[i])) return 0;
573         return 1;
574 }
575
576 /* check whether two bounds intersect inside defined proximity */
577 static int boundsIntersectDist(Bounds3D *b1, Bounds3D *b2, float dist)
578 {
579         int i=2;
580         if (!b1->valid || !b2->valid) return 0;
581         for (; i>=0; i-=1)
582                 if (!(b1->min[i] <= (b2->max[i]+dist) && b1->max[i] >= (b2->min[i]-dist))) return 0;
583         return 1;
584 }
585
586 /* check whether bounds intersects a point with given radius */
587 static int boundIntersectPoint(Bounds3D *b, float point[3], float radius)
588 {
589         int i=2;
590         if (!b->valid) return 0;
591         for (; i>=0; i-=1)
592                 if (!(b->min[i] <= (point[i]+radius) && b->max[i] >= (point[i]-radius))) return 0;
593         return 1;
594 }
595
596 /* expand bounds by a new point */
597 static void boundInsert(Bounds3D *b, float point[3])
598 {
599         int i=2;
600         if (!b->valid) {
601                 VECCOPY(b->min, point);
602                 VECCOPY(b->max, point);
603                 b->valid = 1;
604         }
605         else {
606                 for (; i>=0; i-=1) {
607                         if (point[i] < b->min[i]) b->min[i]=point[i];
608                         if (point[i] > b->max[i]) b->max[i]=point[i];
609                 }
610         }
611 }
612
613 static void freeGrid(PaintSurfaceData *data)
614 {
615         PaintBakeData *bData = data->bData;
616         VolumeGrid *grid = bData->grid;
617
618         if (grid->bounds) MEM_freeN(grid->bounds);
619         if (grid->s_pos) MEM_freeN(grid->s_pos);
620         if (grid->s_num) MEM_freeN(grid->s_num);
621         if (grid->t_index) MEM_freeN(grid->t_index);
622
623         MEM_freeN(bData->grid);
624         bData->grid = NULL;
625 }
626
627 static void surfaceGenerateGrid(struct DynamicPaintSurface *surface)
628 {
629         PaintSurfaceData *sData = surface->data;
630         PaintBakeData *bData = sData->bData;
631         Bounds3D *grid_bounds;
632         VolumeGrid *grid;
633         int grid_cells, axis = 3;
634         int *temp_t_index = NULL;
635         int *temp_s_num = NULL;
636
637 #ifdef _OPENMP
638         int num_of_threads = omp_get_max_threads();
639 #else
640         int num_of_threads = 1;
641 #endif
642
643         if (bData->grid)
644                 freeGrid(sData);
645
646         /* allocate separate bounds for each thread */
647         grid_bounds = MEM_callocN(sizeof(Bounds3D)*num_of_threads, "Grid Bounds");
648
649         bData->grid = MEM_callocN(sizeof(VolumeGrid), "Surface Grid");
650         grid = bData->grid;
651
652         if (grid && grid_bounds) {
653                 int index, error = 0;
654                 float dim_factor, volume, dim[3];
655                 float tx,ty,tz;
656                 float min_dim;
657
658                 /* calculate canvas dimensions */
659                 #ifdef _OPENMP
660                 #pragma omp parallel for schedule(static)
661                 #endif
662                 for (index = 0; index < sData->total_points; index++) {
663 #ifdef _OPENMP
664                         int id = omp_get_thread_num();
665                         boundInsert(&grid_bounds[id], (bData->realCoord[bData->s_pos[index]].v));
666 #else
667                         boundInsert(&grid_bounds[0], (bData->realCoord[bData->s_pos[index]].v));
668 #endif
669                 }
670
671                 /* get final dimensions */
672                 for (index = 0; index<num_of_threads; index++) {
673                         boundInsert(&grid->grid_bounds, grid_bounds[index].min);
674                         boundInsert(&grid->grid_bounds, grid_bounds[index].max);
675                 }
676
677                 dim[0] = grid->grid_bounds.max[0]-grid->grid_bounds.min[0];
678                 dim[1] = grid->grid_bounds.max[1]-grid->grid_bounds.min[1];
679                 dim[2] = grid->grid_bounds.max[2]-grid->grid_bounds.min[2];
680
681                 tx = dim[0];
682                 ty = dim[1];
683                 tz = dim[2];
684                 min_dim = MAX3(tx,ty,tz) / 1000.f;
685
686                 /* deactivate zero axises */
687                 if (tx<min_dim) {tx=1.0f; axis-=1;}
688                 if (ty<min_dim) {ty=1.0f; axis-=1;}
689                 if (tz<min_dim) {tz=1.0f; axis-=1;}
690
691                 if (axis == 0 || MAX3(tx,ty,tz) < 0.0001f)
692                         return;
693
694                 /* now calculate grid volume/area/width depending on num of active axis */
695                 volume = tx*ty*tz;
696
697                 /* determine final grid size by trying to fit average 10.000 points per grid cell */
698                 dim_factor = pow(volume / ((double)sData->total_points / 10000.f), 1.0f/axis);
699
700                 /* define final grid size using dim_factor, use min 3 for active axises */
701                 grid->x = (int)floor(tx / dim_factor);
702                 CLAMP(grid->x, (dim[0]>=min_dim) ? 3 : 1, 100);
703                 grid->y = (int)floor(ty / dim_factor);
704                 CLAMP(grid->y, (dim[1]>=min_dim) ? 3 : 1, 100);
705                 grid->z = (int)floor(tz / dim_factor);
706                 CLAMP(grid->z, (dim[2]>=min_dim) ? 3 : 1, 100);
707
708                 grid_cells = grid->x*grid->y*grid->z;
709
710                 //printf("final grid size %i,%i,%i\n", grid->x, grid->y, grid->z);
711
712                 /* allocate memory for grids */
713
714                 grid->bounds = MEM_callocN(sizeof(Bounds3D) * grid_cells, "Surface Grid Bounds");
715                 grid->s_pos = MEM_callocN(sizeof(int) * grid_cells, "Surface Grid Position");
716                 grid->s_num = MEM_callocN(sizeof(int) * grid_cells*num_of_threads, "Surface Grid Points");
717                 temp_s_num = MEM_callocN(sizeof(int) * grid_cells, "Temp Surface Grid Points");
718                 grid->t_index = MEM_callocN(sizeof(int) * sData->total_points, "Surface Grid Target Ids");
719                 temp_t_index = MEM_callocN(sizeof(int) * sData->total_points, "Temp Surface Grid Target Ids");
720
721                 if (!grid->bounds || !grid->s_pos || !grid->s_num || !grid->t_index || !temp_s_num || !temp_t_index)
722                         error = 1;
723
724                 if (!error) {
725                         /* calculate number of points withing each cell */
726                         #ifdef _OPENMP
727                         #pragma omp parallel for schedule(static)
728                         #endif
729                         for (index = 0; index < sData->total_points; index++) {
730                                 int x,y,z;
731                                 x = floor((bData->realCoord[bData->s_pos[index]].v[0] - grid->grid_bounds.min[0])/dim[0]*grid->x);
732                                 CLAMP(x, 0, grid->x-1);
733                                 y = floor((bData->realCoord[bData->s_pos[index]].v[1] - grid->grid_bounds.min[1])/dim[1]*grid->y);
734                                 CLAMP(y, 0, grid->y-1);
735                                 z = floor((bData->realCoord[bData->s_pos[index]].v[2] - grid->grid_bounds.min[2])/dim[2]*grid->z);
736                                 CLAMP(z, 0, grid->z-1);
737
738                                 temp_t_index[index] = x + y * grid->x + z * grid->x*grid->y;
739 #ifdef _OPENMP
740                                 grid->s_num[temp_t_index[index]+omp_get_thread_num()*grid_cells]++;
741 #else
742                                 grid->s_num[temp_t_index[index]]++;
743 #endif
744                         }
745
746                         /* for first cell only calc s_num */
747                         for (index = 1; index<num_of_threads; index++) {
748                                 grid->s_num[0] += grid->s_num[index*grid_cells];
749                         }
750
751                         /* calculate grid indexes */
752                         for (index = 1; index < grid_cells; index++) {
753                                 int id;
754                                 for (id = 1; id<num_of_threads; id++) {
755                                         grid->s_num[index] += grid->s_num[index+id*grid_cells];
756                                 }
757                                 grid->s_pos[index] = grid->s_pos[index-1] + grid->s_num[index-1];
758                         }
759
760                         /* save point indexes to final array */
761                         for (index = 0; index < sData->total_points; index++) {
762                                 int pos = grid->s_pos[temp_t_index[index]] + temp_s_num[temp_t_index[index]];
763                                 grid->t_index[pos] = index;
764
765                                 temp_s_num[temp_t_index[index]]++;
766                         }
767
768                         /* calculate cell bounds */
769                         {
770                                 int x;
771                                 #ifdef _OPENMP
772                                 #pragma omp parallel for schedule(static)
773                                 #endif
774                                 for (x=0; x<grid->x; x++) {
775                                         int y;
776                                         for (y=0; y<grid->y; y++) {
777                                                 int z;
778                                                 for (z=0; z<grid->z; z++) {
779                                                         int b_index = x + y * grid->x + z * grid->x*grid->y;
780
781                                                         /* set bounds */
782                                                         grid->bounds[b_index].min[0] = grid->grid_bounds.min[0] + dim[0]/grid->x*x;
783                                                         grid->bounds[b_index].min[1] = grid->grid_bounds.min[1] + dim[1]/grid->y*y;
784                                                         grid->bounds[b_index].min[2] = grid->grid_bounds.min[2] + dim[2]/grid->z*z;
785
786                                                         grid->bounds[b_index].max[0] = grid->grid_bounds.min[0] + dim[0]/grid->x*(x+1);
787                                                         grid->bounds[b_index].max[1] = grid->grid_bounds.min[1] + dim[1]/grid->y*(y+1);
788                                                         grid->bounds[b_index].max[2] = grid->grid_bounds.min[2] + dim[2]/grid->z*(z+1);
789
790                                                         grid->bounds[b_index].valid = 1;
791                                                 }
792                                         }
793                                 }
794                         }
795                 }
796
797                 if (temp_s_num) MEM_freeN(temp_s_num);
798                 if (temp_t_index) MEM_freeN(temp_t_index);
799
800                 /* free per thread s_num values */
801                 grid->s_num = MEM_reallocN(grid->s_num, sizeof(int) * grid_cells);
802
803                 if (error || !grid->s_num)
804                         freeGrid(sData);
805         }
806
807         if (grid_bounds) MEM_freeN(grid_bounds);
808 }
809
810 /***************************** Freeing data ******************************/
811
812 /* Free brush data */
813 void dynamicPaint_freeBrush(struct DynamicPaintModifierData *pmd)
814 {
815         if(pmd->brush) {
816                 if(pmd->brush->dm)
817                         pmd->brush->dm->release(pmd->brush->dm);
818                 pmd->brush->dm = NULL;
819
820                 if(pmd->brush->paint_ramp)
821                          MEM_freeN(pmd->brush->paint_ramp);
822                 pmd->brush->paint_ramp = NULL;
823                 if(pmd->brush->vel_ramp)
824                          MEM_freeN(pmd->brush->vel_ramp);
825                 pmd->brush->vel_ramp = NULL;
826
827                 MEM_freeN(pmd->brush);
828                 pmd->brush = NULL;
829         }
830 }
831
832 static void dynamicPaint_freeAdjData(PaintSurfaceData *data)
833 {
834         if (data->adj_data) {
835                 if (data->adj_data->n_index) MEM_freeN(data->adj_data->n_index);
836                 if (data->adj_data->n_num) MEM_freeN(data->adj_data->n_num);
837                 if (data->adj_data->n_target) MEM_freeN(data->adj_data->n_target);
838                 if (data->adj_data->flags) MEM_freeN(data->adj_data->flags);
839                 MEM_freeN(data->adj_data);
840                 data->adj_data = NULL;
841         }
842 }
843
844 static void free_bakeData(PaintSurfaceData *data)
845 {
846         PaintBakeData *bData = data->bData;
847         if (bData) {
848                 if (bData->bNormal) MEM_freeN(bData->bNormal);
849                 if (bData->s_pos) MEM_freeN(bData->s_pos);
850                 if (bData->s_num) MEM_freeN(bData->s_num);
851                 if (bData->realCoord) MEM_freeN(bData->realCoord);
852                 if (bData->bNeighs) MEM_freeN(bData->bNeighs);
853                 if (bData->grid) freeGrid(data);
854                 if (bData->prev_verts) MEM_freeN(bData->prev_verts);
855                 if (bData->velocity) MEM_freeN(bData->velocity);
856                 if (bData->prev_velocity) MEM_freeN(bData->prev_velocity);
857
858                 MEM_freeN(data->bData);
859                 data->bData = NULL;
860         }
861 }
862
863 /* free surface data if it's not used anymore */
864 void surface_freeUnusedData(DynamicPaintSurface *surface)
865 {
866         if (!surface->data) return;
867
868         /* free bakedata if not active or surface is baked */
869         if (!(surface->flags & MOD_DPAINT_ACTIVE) ||
870                 (surface->pointcache && surface->pointcache->flag & PTCACHE_BAKED))
871                 free_bakeData(surface->data);
872 }
873
874 static void dynamicPaint_freeSurfaceData(DynamicPaintSurface *surface)
875 {
876         PaintSurfaceData *data = surface->data;
877         if (!data) return;
878         if (data->format_data) {
879                 /* format specific free */
880                 if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) {
881                         ImgSeqFormatData *format_data = (ImgSeqFormatData*)data->format_data;
882                         if (format_data->uv_p)
883                                 MEM_freeN(format_data->uv_p);
884                         if (format_data->barycentricWeights)
885                                 MEM_freeN(format_data->barycentricWeights);
886                 }
887                 MEM_freeN(data->format_data);
888         }
889         /* type data */
890         if (data->type_data) MEM_freeN(data->type_data);
891         dynamicPaint_freeAdjData(data);
892         /* bake data */
893         free_bakeData(data);
894
895         MEM_freeN(surface->data);
896         surface->data = NULL;
897 }
898
899 void dynamicPaint_freeSurface(DynamicPaintSurface *surface)
900 {
901         if (!surface) return;
902
903         /* point cache */
904         BKE_ptcache_free_list(&(surface->ptcaches));
905         surface->pointcache = NULL;
906
907         if(surface->effector_weights)
908                 MEM_freeN(surface->effector_weights);
909         surface->effector_weights = NULL;
910
911         BLI_remlink(&(surface->canvas->surfaces), surface);
912         dynamicPaint_freeSurfaceData(surface);
913         MEM_freeN(surface);
914 }
915
916 /* Free canvas data */
917 void dynamicPaint_freeCanvas(DynamicPaintModifierData *pmd)
918 {
919         if(pmd->canvas) {
920                 /* Free surface data */
921                 DynamicPaintSurface *surface = pmd->canvas->surfaces.first;
922                 DynamicPaintSurface *next_surface = NULL;
923
924                 while (surface) {
925                         next_surface = surface->next;
926                         dynamicPaint_freeSurface(surface);
927                         surface = next_surface;
928                 }
929
930                 /* free dm copy */
931                 if (pmd->canvas->dm)
932                         pmd->canvas->dm->release(pmd->canvas->dm);
933                 pmd->canvas->dm = NULL;
934
935                 MEM_freeN(pmd->canvas);
936                 pmd->canvas = NULL;
937         }
938 }
939
940 /* Free whole dp modifier */
941 void dynamicPaint_Modifier_free(struct DynamicPaintModifierData *pmd)
942 {
943         if(pmd) {
944                 dynamicPaint_freeCanvas(pmd);
945                 dynamicPaint_freeBrush(pmd);
946         }
947 }
948
949
950 /***************************** Initialize and reset ******************************/
951
952 /*
953 *       Creates a new surface and adds it to the list
954 *       A pointer to this surface is returned
955 */
956 struct DynamicPaintSurface *dynamicPaint_createNewSurface(DynamicPaintCanvasSettings *canvas, Scene *scene)
957 {
958         DynamicPaintSurface *surface= MEM_callocN(sizeof(DynamicPaintSurface), "DynamicPaintSurface");
959         if (!surface) return NULL;
960
961         surface->canvas = canvas;
962         surface->format = MOD_DPAINT_SURFACE_F_VERTEX;
963         surface->type = MOD_DPAINT_SURFACE_T_PAINT;
964
965         /* cache */
966         surface->pointcache = BKE_ptcache_add(&(surface->ptcaches));
967         surface->pointcache->flag |= PTCACHE_DISK_CACHE;
968         surface->pointcache->step = 1;
969
970         /* Set initial values */
971         surface->flags = MOD_DPAINT_ANTIALIAS | MOD_DPAINT_MULALPHA | MOD_DPAINT_DRY_LOG | MOD_DPAINT_DISSOLVE_LOG |
972                                          MOD_DPAINT_ACTIVE | MOD_DPAINT_PREVIEW | MOD_DPAINT_OUT1;
973         surface->effect = 0;
974         surface->effect_ui = 1;
975
976         surface->diss_speed = 300;
977         surface->dry_speed = 300;
978         surface->disp_clamp = 0.0f;
979         surface->disp_type = MOD_DPAINT_DISP_DISPLACE;
980         surface->image_fileformat = MOD_DPAINT_IMGFORMAT_PNG;
981
982         surface->init_color[0] = 1.0f;
983         surface->init_color[1] = 1.0f;
984         surface->init_color[2] = 1.0f;
985         surface->init_color[3] = 1.0f;
986
987         surface->image_resolution = 256;
988         surface->substeps = 0;
989
990         if (scene) {
991                 surface->start_frame = scene->r.sfra;
992                 surface->end_frame = scene->r.efra;
993         }
994         else {
995                 surface->start_frame = 1;
996                 surface->end_frame = 250;
997         }
998
999         surface->spread_speed = 1.0f;
1000         surface->color_spread_speed = 1.0f;
1001         surface->shrink_speed = 1.0f;
1002
1003         surface->wave_damping = 0.05f;
1004         surface->wave_speed = 0.8f;
1005         surface->wave_timescale = 1.0f;
1006         surface->wave_spring = 0.20;
1007
1008         sprintf(surface->image_output_path, "%sdynamicpaint/", "/tmp/");
1009         dynamicPaintSurface_setUniqueName(surface, "Surface");
1010
1011         surface->effector_weights = BKE_add_effector_weights(NULL);
1012
1013         dynamicPaintSurface_updateType(surface);
1014
1015         BLI_addtail(&canvas->surfaces, surface);
1016
1017         return surface;
1018 }
1019
1020 /*
1021 *       Initialize modifier data
1022 */
1023 int dynamicPaint_createType(struct DynamicPaintModifierData *pmd, int type, struct Scene *scene)
1024 {
1025         if(pmd)
1026         {
1027                 if(type == MOD_DYNAMICPAINT_TYPE_CANVAS)
1028                 {
1029                         if(pmd->canvas)
1030                                 dynamicPaint_freeCanvas(pmd);
1031
1032                         pmd->canvas = MEM_callocN(sizeof(DynamicPaintCanvasSettings), "DynamicPaint Canvas");
1033                         if (!pmd->canvas)
1034                                 return 0;
1035                         pmd->canvas->pmd = pmd;
1036                         pmd->canvas->dm = NULL;
1037
1038                         /* Create one surface */
1039                         if (!dynamicPaint_createNewSurface(pmd->canvas, scene))
1040                                 return 0;
1041
1042                         pmd->canvas->ui_info[0] = '\0';
1043
1044                 }
1045                 else if(type == MOD_DYNAMICPAINT_TYPE_BRUSH)
1046                 {
1047                         if(pmd->brush)
1048                                 dynamicPaint_freeBrush(pmd);
1049
1050                         pmd->brush = MEM_callocN(sizeof(DynamicPaintBrushSettings), "DynamicPaint Paint");
1051                         if (!pmd->brush)
1052                                 return 0;
1053                         pmd->brush->pmd = pmd;
1054
1055                         pmd->brush->psys = NULL;
1056
1057                         pmd->brush->flags = MOD_DPAINT_ABS_ALPHA | MOD_DPAINT_RAMP_ALPHA;
1058                         pmd->brush->collision = MOD_DPAINT_COL_VOLUME;
1059                         
1060                         pmd->brush->mat = NULL;
1061                         pmd->brush->r = 0.0f;
1062                         pmd->brush->g = 0.0f;
1063                         pmd->brush->b = 1.0f;
1064                         pmd->brush->alpha = 1.0f;
1065                         pmd->brush->wetness = 1.0f;
1066
1067                         pmd->brush->paint_distance = 1.0f;
1068                         pmd->brush->proximity_falloff = MOD_DPAINT_PRFALL_SMOOTH;
1069
1070                         pmd->brush->particle_radius = 0.2f;
1071                         pmd->brush->particle_smooth = 0.05f;
1072
1073                         pmd->brush->wave_factor = 1.0f;
1074                         pmd->brush->wave_clamp = 0.0f;
1075                         pmd->brush->smudge_strength = 0.3f;
1076                         pmd->brush->max_velocity = 1.0f;
1077
1078                         pmd->brush->dm = NULL;
1079
1080                         /* Paint proximity falloff colorramp. */
1081                         {
1082                                 CBData *ramp;
1083
1084                                 pmd->brush->paint_ramp = add_colorband(0);
1085                                 if (!pmd->brush->paint_ramp)
1086                                         return 0;
1087                                 ramp = pmd->brush->paint_ramp->data;
1088                                 /* Add default smooth-falloff ramp.     */
1089                                 ramp[0].r = ramp[0].g = ramp[0].b = ramp[0].a = 1.0f;
1090                                 ramp[0].pos = 0.0f;
1091                                 ramp[1].r = ramp[1].g = ramp[1].b = ramp[1].pos = 1.0f;
1092                                 ramp[1].a = 0.0f;
1093                                 pmd->brush->paint_ramp->tot = 2;
1094                         }
1095
1096                         /* Brush velocity ramp. */
1097                         {
1098                                 CBData *ramp;
1099
1100                                 pmd->brush->vel_ramp = add_colorband(0);
1101                                 if (!pmd->brush->vel_ramp)
1102                                         return 0;
1103                                 ramp = pmd->brush->vel_ramp->data;
1104                                 ramp[0].r = ramp[0].g = ramp[0].b = ramp[0].a = ramp[0].pos = 0.0f;
1105                                 ramp[1].r = ramp[1].g = ramp[1].b = ramp[1].a = ramp[1].pos = 1.0f;
1106                                 pmd->brush->paint_ramp->tot = 2;
1107                         }
1108                 }
1109         }
1110         else
1111                 return 0;
1112
1113         return 1;
1114 }
1115
1116 void dynamicPaint_Modifier_copy(struct DynamicPaintModifierData *pmd, struct DynamicPaintModifierData *tpmd)
1117 {
1118         /* Init modifier        */
1119         tpmd->type = pmd->type;
1120         if (pmd->canvas)
1121                 dynamicPaint_createType(tpmd, MOD_DYNAMICPAINT_TYPE_CANVAS, NULL);
1122         if (pmd->brush)
1123                 dynamicPaint_createType(tpmd, MOD_DYNAMICPAINT_TYPE_BRUSH, NULL);
1124
1125         /* Copy data    */
1126         if (tpmd->canvas) {
1127                 tpmd->canvas->pmd = tpmd;
1128
1129                 tpmd->canvas->ui_info[0] = '\0';
1130
1131         } else if (tpmd->brush) {
1132                 tpmd->brush->pmd = tpmd;
1133
1134                 tpmd->brush->flags = pmd->brush->flags;
1135                 tpmd->brush->collision = pmd->brush->collision;
1136
1137                 tpmd->brush->mat = pmd->brush->mat;
1138                 tpmd->brush->r = pmd->brush->r;
1139                 tpmd->brush->g = pmd->brush->g;
1140                 tpmd->brush->b = pmd->brush->b;
1141                 tpmd->brush->alpha = pmd->brush->alpha;
1142                 tpmd->brush->wetness = pmd->brush->wetness;
1143
1144                 tpmd->brush->particle_radius = pmd->brush->particle_radius;
1145                 tpmd->brush->particle_smooth = pmd->brush->particle_smooth;
1146                 tpmd->brush->paint_distance = pmd->brush->paint_distance;
1147                 tpmd->brush->psys = pmd->brush->psys;
1148
1149                 if (pmd->brush->paint_ramp)
1150                         memcpy(tpmd->brush->paint_ramp, pmd->brush->paint_ramp, sizeof(ColorBand));
1151                 if (pmd->brush->vel_ramp)
1152                         memcpy(tpmd->brush->vel_ramp, pmd->brush->vel_ramp, sizeof(ColorBand));
1153
1154                 tpmd->brush->proximity_falloff = pmd->brush->proximity_falloff;
1155                 tpmd->brush->brush_settings_context = pmd->brush->brush_settings_context;
1156                 tpmd->brush->wave_type = pmd->brush->wave_type;
1157                 tpmd->brush->ray_dir = pmd->brush->ray_dir;
1158
1159                 tpmd->brush->wave_factor = pmd->brush->wave_factor;
1160                 tpmd->brush->wave_clamp = pmd->brush->wave_clamp;
1161                 tpmd->brush->max_velocity = pmd->brush->max_velocity;
1162                 tpmd->brush->smudge_strength = pmd->brush->smudge_strength;
1163         }
1164 }
1165
1166 /* allocates surface data depending on surface type */
1167 static void dynamicPaint_allocateSurfaceType(DynamicPaintSurface *surface)
1168 {
1169         PaintSurfaceData *sData = surface->data;
1170
1171         if (surface->type == MOD_DPAINT_SURFACE_T_PAINT) {
1172                 sData->type_data = MEM_callocN(sizeof(PaintPoint)*sData->total_points, "DynamicPaintSurface Data");
1173         }
1174         else if (surface->type == MOD_DPAINT_SURFACE_T_DISPLACE) {
1175                 sData->type_data = MEM_callocN(sizeof(float)*sData->total_points, "DynamicPaintSurface DepthData");
1176         }
1177         else if (surface->type == MOD_DPAINT_SURFACE_T_WEIGHT) {
1178                 sData->type_data = MEM_callocN(sizeof(float)*sData->total_points, "DynamicPaintSurface WeightData");
1179         }
1180         else if (surface->type == MOD_DPAINT_SURFACE_T_WAVE) {
1181                 sData->type_data = MEM_callocN(sizeof(PaintWavePoint)*sData->total_points, "DynamicPaintSurface WaveData");
1182         }
1183         else return;
1184
1185         if (sData->type_data == NULL) printError(surface->canvas, "Not enough free memory!");
1186 }
1187
1188 static int surface_usesAdjDistance(DynamicPaintSurface *surface)
1189 {
1190         if (surface->type == MOD_DPAINT_SURFACE_T_PAINT && surface->effect) return 1;
1191         if (surface->type == MOD_DPAINT_SURFACE_T_WAVE) return 1;
1192         return 0;
1193 }
1194
1195 static int surface_usesAdjData(DynamicPaintSurface *surface)
1196 {
1197         if (surface_usesAdjDistance(surface)) return 1;
1198         if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX &&
1199                 surface->flags & MOD_DPAINT_ANTIALIAS) return 1;
1200
1201         return 0;
1202 }
1203
1204 /* initialize surface adjacency data */
1205 static void dynamicPaint_initAdjacencyData(DynamicPaintSurface *surface, int force_init)
1206 {
1207         PaintSurfaceData *sData = surface->data;
1208         PaintAdjData *ed;
1209         int *temp_data;
1210         int neigh_points = 0;
1211
1212         if (!surface_usesAdjData(surface) && !force_init) return;
1213
1214         if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
1215                 /* For vertex format, neighbours are connected by edges */
1216                 neigh_points = 2*surface->canvas->dm->getNumEdges(surface->canvas->dm);
1217         }
1218         else if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ)
1219                 neigh_points = sData->total_points*8;
1220
1221         if (!neigh_points) return;
1222
1223         /* allocate memory */
1224         ed = sData->adj_data = MEM_callocN(sizeof(PaintAdjData), "Surface Adj Data");
1225         if (!ed) return;
1226         ed->n_index = MEM_callocN(sizeof(int)*sData->total_points, "Surface Adj Index");
1227         ed->n_num = MEM_callocN(sizeof(int)*sData->total_points, "Surface Adj Counts");
1228         temp_data = MEM_callocN(sizeof(int)*sData->total_points, "Temp Adj Data");
1229         ed->n_target = MEM_callocN(sizeof(int)*neigh_points, "Surface Adj Targets");
1230         ed->flags = MEM_callocN(sizeof(int)*sData->total_points, "Surface Adj Flags");
1231         ed->total_targets = neigh_points;
1232
1233         /* in case of error, free allocated memory */
1234         if (!ed->n_index || !ed->n_num || !ed->n_target || !temp_data) {
1235                 dynamicPaint_freeAdjData(sData);
1236                 MEM_freeN(temp_data);
1237                 printError(surface->canvas, "Not enough free memory.");
1238                 return;
1239         }
1240
1241         if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
1242                 int i;
1243                 int n_pos;
1244
1245                 /* For vertex format, count every vertex that is connected by an edge */
1246                 int numOfEdges = surface->canvas->dm->getNumEdges(surface->canvas->dm);
1247                 int numOfFaces = surface->canvas->dm->getNumFaces(surface->canvas->dm);
1248                 struct MEdge *edge =  surface->canvas->dm->getEdgeArray(surface->canvas->dm);
1249                 struct MFace *face =  surface->canvas->dm->getFaceArray(surface->canvas->dm);
1250
1251                 /* count number of edges per vertex */
1252                 for (i=0; i<numOfEdges; i++) {
1253                         ed->n_num[edge[i].v1]++;
1254                         ed->n_num[edge[i].v2]++;
1255
1256                         temp_data[edge[i].v1]++;
1257                         temp_data[edge[i].v2]++;
1258                 }
1259
1260                 /* to locate points on "mesh edge" */
1261                 for (i=0; i<numOfFaces; i++) {
1262                         temp_data[face[i].v1]++;
1263                         temp_data[face[i].v2]++;
1264                         temp_data[face[i].v3]++;
1265                         if (face[i].v4)
1266                                 temp_data[face[i].v4]++;
1267                 }
1268
1269                 /* now check if total number of edges+faces for
1270                 *  each vertex is even, if not -> vertex is on mesh edge */
1271                 for (i=0; i<sData->total_points; i++) {
1272                         if ((temp_data[i]%2) ||
1273                                 temp_data[i] < 4)
1274                                 ed->flags[i] |= ADJ_ON_MESH_EDGE;
1275                                 
1276                         /* reset temp data */ 
1277                         temp_data[i] = 0;
1278                 }
1279
1280                 /* order n_index array */
1281                 n_pos = 0;
1282                 for (i=0; i<sData->total_points; i++) {
1283                         ed->n_index[i] = n_pos;
1284                         n_pos += ed->n_num[i];
1285                 }
1286
1287                 /* and now add neighbour data using that info */
1288                 for (i=0; i<numOfEdges; i++) {
1289                         /* first vertex */
1290                         int index = edge[i].v1;
1291                         n_pos = ed->n_index[index]+temp_data[index];
1292                         ed->n_target[n_pos] = edge[i].v2;
1293                         temp_data[index]++;
1294
1295                         /* second vertex */
1296                         index = edge[i].v2;
1297                         n_pos = ed->n_index[index]+temp_data[index];
1298                         ed->n_target[n_pos] = edge[i].v1;
1299                         temp_data[index]++;
1300                 }
1301         }
1302         else if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) {
1303                 /* for image sequences, only allocate memory.
1304                 *  bake initialization takes care of rest */
1305         }
1306
1307         MEM_freeN(temp_data);
1308 }
1309
1310 void dynamicPaint_setInitialColor(DynamicPaintSurface *surface)
1311 {
1312         PaintSurfaceData *sData = surface->data;
1313         PaintPoint* pPoint = (PaintPoint*)sData->type_data;
1314         DerivedMesh *dm = surface->canvas->dm;
1315         int i;
1316
1317         if (surface->type != MOD_DPAINT_SURFACE_T_PAINT)
1318                 return;
1319
1320         if (surface->init_color_type == MOD_DPAINT_INITIAL_NONE)
1321                 return;
1322         /* Single color */
1323         else if (surface->init_color_type == MOD_DPAINT_INITIAL_COLOR) {
1324                 /* apply color to every surface point */
1325                 #ifdef _OPENMP
1326                 #pragma omp parallel for schedule(static)
1327                 #endif
1328                 for (i=0; i<sData->total_points; i++) {
1329                         VECCOPY(pPoint[i].color, surface->init_color);
1330                         pPoint[i].alpha = surface->init_color[3];
1331                 }
1332         }
1333         /* UV mapped texture */
1334         else if (surface->init_color_type == MOD_DPAINT_INITIAL_TEXTURE) {
1335                 Tex *tex = surface->init_texture;
1336                 MTFace *tface;
1337                 MFace *mface = dm->getFaceArray(dm);
1338                 int numOfFaces = dm->getNumFaces(dm);
1339                 char uvname[40];
1340
1341                 if (!tex) return;
1342
1343                 /* get uv layer */
1344                 validate_layer_name(&dm->faceData, CD_MTFACE, surface->init_layername, uvname);
1345                 tface = CustomData_get_layer_named(&dm->faceData, CD_MTFACE, uvname);
1346                 if (!tface) return;
1347
1348                 /* for vertex surface loop through tfaces and find uv color
1349                 *  that provides highest alpha */
1350                 if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
1351                         #ifdef _OPENMP
1352                         #pragma omp parallel for schedule(static)
1353                         #endif
1354                         for (i=0; i<numOfFaces; i++) {
1355                                 int numOfVert = (mface[i].v4) ? 4 : 3;
1356                                 float uv[3] = {0.0f};
1357                                 int j;
1358                                 for (j=0; j<numOfVert; j++) {
1359                                         TexResult texres = {0};
1360                                         unsigned int *vert = (&mface[i].v1)+j;
1361
1362                                         /* remap to -1.0 to 1.0 */
1363                                         uv[0] = tface[i].uv[j][0]*2.0f - 1.0f;
1364                                         uv[1] = tface[i].uv[j][1]*2.0f - 1.0f;
1365
1366                                         multitex_ext_safe(tex, uv, &texres);
1367
1368                                         if (texres.tin > pPoint[*vert].alpha) {
1369                                                 VECCOPY(pPoint[*vert].color, &texres.tr);
1370                                                 pPoint[*vert].alpha = texres.tin;
1371                                         }
1372                                 }
1373                         }
1374                 }
1375                 else if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) {
1376                         ImgSeqFormatData *f_data = (ImgSeqFormatData*)sData->format_data;
1377                         int samples = (surface->flags & MOD_DPAINT_ANTIALIAS) ? 5 : 1;
1378
1379                         #ifdef _OPENMP
1380                         #pragma omp parallel for schedule(static)
1381                         #endif
1382                         for (i=0; i<sData->total_points; i++) {
1383                                 float uv[9] = {0.0f};
1384                                 float uv_final[3] = {0.0f};
1385                                 int j;
1386                                 TexResult texres = {0};
1387
1388                                 /* collect all uvs */
1389                                 for (j=0; j<3; j++) {
1390                                         int v=(f_data->uv_p[i].quad && j>0) ? j+1 : j;
1391                                         VECCOPY2D(&uv[j*3], tface[f_data->uv_p[i].face_index].uv[v]);
1392                                 }
1393
1394                                 /* interpolate final uv pos */
1395                                 interp_v3_v3v3v3(       uv_final, &uv[0], &uv[3], &uv[6],
1396                                         f_data->barycentricWeights[i*samples].v);
1397                                 /* remap to -1.0 to 1.0 */
1398                                 uv_final[0] = uv_final[0]*2.0f - 1.0f;
1399                                 uv_final[1] = uv_final[1]*2.0f - 1.0f;
1400                                         
1401                                 multitex_ext_safe(tex, uv_final, &texres);
1402
1403                                 /* apply color */
1404                                 VECCOPY(pPoint[i].color, &texres.tr);
1405                                 pPoint[i].alpha = texres.tin;
1406                         }
1407                 }
1408         }
1409         /* vertex color layer */
1410         else if (surface->init_color_type == MOD_DPAINT_INITIAL_VERTEXCOLOR) {
1411                 MCol *col = CustomData_get_layer_named(&dm->faceData, CD_MCOL, surface->init_layername);
1412                 if (!col) return;
1413
1414                 /* for vertex surface, just copy colors from mcol */
1415                 if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
1416                         MFace *mface = dm->getFaceArray(dm);
1417                         int numOfFaces = dm->getNumFaces(dm);
1418
1419                         #ifdef _OPENMP
1420                         #pragma omp parallel for schedule(static)
1421                         #endif
1422                         for (i=0; i<numOfFaces; i++) {
1423                                 int numOfVert = (mface[i].v4) ? 4 : 3;
1424                                 int j;
1425                                 for (j=0; j<numOfVert; j++) {
1426                                         unsigned int *vert = ((&mface[i].v1)+j);
1427
1428                                         pPoint[*vert].color[0] = 1.0f/255.f*(float)col[i*4+j].b;
1429                                         pPoint[*vert].color[1] = 1.0f/255.f*(float)col[i*4+j].g;
1430                                         pPoint[*vert].color[2] = 1.0f/255.f*(float)col[i*4+j].r;
1431                                         pPoint[*vert].alpha = 1.0f/255.f*(float)col[i*4+j].a;
1432                                 }
1433                         }
1434                 }
1435                 else if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) {
1436                         ImgSeqFormatData *f_data = (ImgSeqFormatData*)sData->format_data;
1437                         int samples = (surface->flags & MOD_DPAINT_ANTIALIAS) ? 5 : 1;
1438
1439                         #ifdef _OPENMP
1440                         #pragma omp parallel for schedule(static)
1441                         #endif
1442                         for (i=0; i<sData->total_points; i++) {
1443                                 int face_ind = f_data->uv_p[i].face_index;
1444                                 float colors[3][4] = {{0.0f,0.0f,0.0f,0.0f}};
1445                                 float final_color[4];
1446                                 int j;
1447                                 /* collect color values */
1448                                 for (j=0; j<3; j++) {
1449                                         int v=(f_data->uv_p[i].quad && j>0) ? j+1 : j;
1450                                         colors[j][0] = 1.0f/255.f*(float)col[face_ind*4+v].b;
1451                                         colors[j][1] = 1.0f/255.f*(float)col[face_ind*4+v].g;
1452                                         colors[j][2] = 1.0f/255.f*(float)col[face_ind*4+v].r;
1453                                         colors[j][3] = 1.0f/255.f*(float)col[face_ind*4+v].a;
1454                                 }
1455                                 
1456                                 /* interpolate final color */
1457                                 interp_v4_v4v4v4(       final_color, colors[0], colors[1], colors[2],
1458                                                 f_data->barycentricWeights[i*samples].v);
1459
1460                                 VECCOPY(pPoint[i].color, final_color);
1461                                 pPoint[i].alpha = final_color[3];
1462                         }
1463                 }
1464         }
1465 }
1466
1467 /* clears surface data back to zero */
1468 void dynamicPaint_clearSurface(DynamicPaintSurface *surface)
1469 {
1470         PaintSurfaceData *sData = surface->data;
1471         if (sData && sData->type_data) {
1472                 unsigned int data_size;
1473
1474                 if (surface->type == MOD_DPAINT_SURFACE_T_PAINT)
1475                         data_size = sizeof(PaintPoint);
1476                 else if (surface->type == MOD_DPAINT_SURFACE_T_WAVE)
1477                         data_size = sizeof(PaintWavePoint);
1478                 else
1479                         data_size = sizeof(float);
1480
1481                 memset(sData->type_data, 0, data_size * sData->total_points);
1482
1483                 /* set initial color */
1484                 if (surface->type == MOD_DPAINT_SURFACE_T_PAINT)
1485                         dynamicPaint_setInitialColor(surface);
1486
1487                 if (sData->bData)
1488                         sData->bData->clear = 1;
1489         }
1490 }
1491
1492 /* completely (re)initializes surface (only for point cache types)*/
1493 int dynamicPaint_resetSurface(DynamicPaintSurface *surface)
1494 {
1495         int numOfPoints = dynamicPaint_surfaceNumOfPoints(surface);
1496         /* dont touch image sequence types. they get handled only on bake */
1497         if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) return 1;
1498
1499         if (surface->data) dynamicPaint_freeSurfaceData(surface);
1500         if (numOfPoints < 1) return 0;
1501
1502         /* allocate memory */
1503         surface->data = MEM_callocN(sizeof(PaintSurfaceData), "PaintSurfaceData");
1504         if (!surface->data) return 0;
1505
1506         /* allocate data depending on surface type and format */
1507         surface->data->total_points = numOfPoints;
1508         dynamicPaint_allocateSurfaceType(surface);
1509         dynamicPaint_initAdjacencyData(surface, 0);
1510
1511         /* set initial color */
1512         if (surface->type == MOD_DPAINT_SURFACE_T_PAINT)
1513                 dynamicPaint_setInitialColor(surface);
1514
1515         return 1;
1516 }
1517
1518 /* make sure allocated surface size matches current requirements */
1519 static void dynamicPaint_checkSurfaceData(DynamicPaintSurface *surface)
1520 {
1521         if (!surface->data || ((dynamicPaint_surfaceNumOfPoints(surface) != surface->data->total_points))) {
1522                 dynamicPaint_resetSurface(surface);
1523         }
1524 }
1525
1526
1527 /***************************** Modifier processing ******************************/
1528
1529
1530 /* apply displacing vertex surface to the derived mesh */
1531 static void dynamicPaint_applySurfaceDisplace(DynamicPaintSurface *surface, DerivedMesh *result, int update_normals)
1532 {
1533         PaintSurfaceData *sData = surface->data;
1534
1535         if (!sData || surface->format != MOD_DPAINT_SURFACE_F_VERTEX) return;
1536
1537         /* displace paint */
1538         if (surface->type == MOD_DPAINT_SURFACE_T_DISPLACE) {
1539                 MVert *mvert = result->getVertArray(result);
1540                 int i;
1541                 float* value = (float*)sData->type_data;
1542
1543                 #ifdef _OPENMP
1544                 #pragma omp parallel for schedule(static)
1545                 #endif
1546                 for (i=0; i<sData->total_points; i++) {
1547                         float normal[3];
1548                         normal_short_to_float_v3(normal, mvert[i].no);
1549                         normalize_v3(normal);
1550
1551                         mvert[i].co[0] -= normal[0]*value[i];
1552                         mvert[i].co[1] -= normal[1]*value[i];
1553                         mvert[i].co[2] -= normal[2]*value[i];
1554                 }
1555         }
1556         else return;
1557
1558         if (update_normals)
1559                 CDDM_calc_normals(result);
1560 }
1561
1562 /*
1563 *       Apply canvas data to the object derived mesh
1564 */
1565 static struct DerivedMesh *dynamicPaint_Modifier_apply(DynamicPaintModifierData *pmd, Scene *scene, Object *ob, DerivedMesh *dm)
1566 {       
1567         DerivedMesh *result = CDDM_copy(dm);
1568
1569         if(pmd->canvas && !(pmd->canvas->flags & MOD_DPAINT_BAKING)) {
1570
1571                 DynamicPaintSurface *surface = pmd->canvas->surfaces.first;
1572                 pmd->canvas->flags &= ~MOD_DPAINT_PREVIEW_READY;
1573
1574                 /* loop through surfaces */
1575                 for (; surface; surface=surface->next) {
1576                         PaintSurfaceData *sData = surface->data;
1577
1578                         if (surface && surface->format != MOD_DPAINT_SURFACE_F_IMAGESEQ && sData) {
1579                                 if (!(surface->flags & (MOD_DPAINT_ACTIVE))) continue;
1580
1581                                 /* process vertex surface previews */
1582                                 if (surface->format == MOD_DPAINT_SURFACE_F_VERTEX) {
1583
1584                                         /* vertex color paint */
1585                                         if (surface->type == MOD_DPAINT_SURFACE_T_PAINT) {
1586
1587                                                 MFace *mface = result->getFaceArray(result);
1588                                                 int numOfFaces = result->getNumFaces(result);
1589                                                 int i;
1590                                                 PaintPoint* pPoint = (PaintPoint*)sData->type_data;
1591                                                 MCol *col;
1592
1593                                                 /* paint is stored on dry and wet layers, so mix final color first */
1594                                                 float *fcolor = MEM_callocN(sizeof(float)*sData->total_points*4, "Temp paint color");
1595
1596                                                 #ifdef _OPENMP
1597                                                 #pragma omp parallel for schedule(static)
1598                                                 #endif
1599                                                 for (i=0; i<sData->total_points; i++) {
1600                                                         int j=i*4;
1601
1602                                                         fcolor[j]   = pPoint[i].color[0];
1603                                                         fcolor[j+1] = pPoint[i].color[1];
1604                                                         fcolor[j+2] = pPoint[i].color[2];
1605                                                         /* mix colors */
1606                                                         if (pPoint[i].e_alpha) mixColors(&fcolor[j], pPoint[i].alpha, pPoint[i].e_color, pPoint[i].e_alpha);
1607
1608                                                         /* Use highest alpha    */
1609                                                         fcolor[j+3] = (pPoint[i].e_alpha > pPoint[i].alpha) ? pPoint[i].e_alpha : pPoint[i].alpha;
1610                                                 }
1611
1612                                                 /* viewport preview */
1613                                                 if (surface->flags & MOD_DPAINT_PREVIEW) {
1614                                                         /* Save preview results to weight layer, to be
1615                                                         *   able to share same drawing methods */
1616                                                         col = result->getFaceDataArray(result, CD_WEIGHT_MCOL);
1617                                                         if (!col) col = CustomData_add_layer(&result->faceData, CD_WEIGHT_MCOL, CD_CALLOC, NULL, numOfFaces);
1618
1619                                                         if (col) {
1620                                                                 #ifdef _OPENMP
1621                                                                 #pragma omp parallel for schedule(static)
1622                                                                 #endif
1623                                                                 for (i=0; i<numOfFaces; i++) {
1624                                                                         int j=0;
1625                                                                         float invAlpha;
1626                                                                         Material *material = give_current_material(ob, mface[i].mat_nr+1);
1627
1628                                                                         for (; j<((mface[i].v4)?4:3); j++) {
1629                                                                                 int index = (j==0)?mface[i].v1: (j==1)?mface[i].v2: (j==2)?mface[i].v3: mface[i].v4;
1630
1631                                                                                 if (surface->preview_id == MOD_DPAINT_SURFACE_PREV_PAINT) {
1632                                                                                         index *= 4;
1633                                                                                         invAlpha = 1.0f - fcolor[index+3];
1634
1635                                                                                         /* Apply material color as base vertex color for preview */
1636                                                                                         col[i*4+j].a = 255;
1637                                                                                         if (material) {
1638                                                                                                 col[i*4+j].r = (unsigned char)(material->b*255);
1639                                                                                                 col[i*4+j].g = (unsigned char)(material->g*255);
1640                                                                                                 col[i*4+j].b = (unsigned char)(material->r*255);
1641                                                                                         }
1642                                                                                         else {
1643                                                                                                 col[i*4+j].r = 165;
1644                                                                                                 col[i*4+j].g = 165;
1645                                                                                                 col[i*4+j].b = 165;
1646                                                                                         }
1647
1648                                                                                         /* mix surface color */
1649                                                                                         col[i*4+j].r = (char)(((float)col[i*4+j].r)*invAlpha + (fcolor[index+2]*255*fcolor[index+3]));
1650                                                                                         col[i*4+j].g = (char)(((float)col[i*4+j].g)*invAlpha + (fcolor[index+1]*255*fcolor[index+3]));
1651                                                                                         col[i*4+j].b = (char)(((float)col[i*4+j].b)*invAlpha + (fcolor[index]*255*fcolor[index+3]));
1652                                                                                 }
1653                                                                                 else {
1654                                                                                         float wetness = (pPoint[index].wetness<1.0f) ? pPoint[index].wetness : 1.0f;
1655                                                                                         col[i*4+j].a = 255;
1656                                                                                         col[i*4+j].r = (char)(wetness*255);
1657                                                                                         col[i*4+j].g = (char)(wetness*255);
1658                                                                                         col[i*4+j].b = (char)(wetness*255);
1659                                                                                 }
1660                                                                         }
1661                                                                 }
1662                                                                 pmd->canvas->flags |= MOD_DPAINT_PREVIEW_READY;
1663                                                         }
1664                                                 }
1665
1666
1667                                                 /* save layer data to output layer */
1668
1669                                                 /* paint layer */
1670                                                 col = CustomData_get_layer_named(&result->faceData, CD_MCOL, surface->output_name);
1671                                                 if (col) {
1672                                                         #ifdef _OPENMP
1673                                                         #pragma omp parallel for schedule(static)
1674                                                         #endif
1675                                                         for (i=0; i<numOfFaces; i++) {
1676                                                                 int j=0;
1677                                                                 for (; j<((mface[i].v4)?4:3); j++) {
1678                                                                         int index = (j==0)?mface[i].v1: (j==1)?mface[i].v2: (j==2)?mface[i].v3: mface[i].v4;
1679                                                                         index *= 4;
1680
1681                                                                         col[i*4+j].a = (char)(fcolor[index+3]*255);
1682                                                                         col[i*4+j].r = (char)(fcolor[index+2]*255);
1683                                                                         col[i*4+j].g = (char)(fcolor[index+1]*255);
1684                                                                         col[i*4+j].b = (char)(fcolor[index]*255);
1685                                                                 }
1686                                                         }
1687                                                 }
1688                                                 
1689                                                 MEM_freeN(fcolor);
1690
1691                                                 /* wet layer */
1692                                                 col = CustomData_get_layer_named(&result->faceData, CD_MCOL, surface->output_name2);
1693                                                 if (col) {
1694                                                         #ifdef _OPENMP
1695                                                         #pragma omp parallel for schedule(static)
1696                                                         #endif
1697                                                         for (i=0; i<numOfFaces; i++) {
1698                                                                 int j=0;
1699
1700                                                                 for (; j<((mface[i].v4)?4:3); j++) {
1701                                                                         int index = (j==0)?mface[i].v1: (j==1)?mface[i].v2: (j==2)?mface[i].v3: mface[i].v4;
1702
1703                                                                         float wetness = (pPoint[index].wetness<1.0f) ? pPoint[index].wetness : 1.0f;
1704                                                                         col[i*4+j].a = 255;
1705                                                                         col[i*4+j].r = (char)(wetness*255);
1706                                                                         col[i*4+j].g = (char)(wetness*255);
1707                                                                         col[i*4+j].b = (char)(wetness*255);
1708                                                                 }
1709                                                         }
1710                                                 }
1711                                         }
1712                                         /* vertex group paint */
1713                                         else if (surface->type == MOD_DPAINT_SURFACE_T_WEIGHT) {
1714                                                 int defgrp_index = defgroup_name_index(ob, surface->output_name);
1715                                                 MDeformVert *dvert = result->getVertDataArray(result, CD_MDEFORMVERT);
1716                                                 float *weight = (float*)sData->type_data;
1717                                                 /* viewport preview */
1718                                                 if (surface->flags & MOD_DPAINT_PREVIEW) {
1719                                                         /* Save preview results to weight layer, to be
1720                                                         *   able to share same drawing methods */
1721                                                         MFace *mface = result->getFaceArray(result);
1722                                                         int numOfFaces = result->getNumFaces(result);
1723                                                         int i;
1724                                                         MCol *col = result->getFaceDataArray(result, CD_WEIGHT_MCOL);
1725                                                         if (!col) col = CustomData_add_layer(&result->faceData, CD_WEIGHT_MCOL, CD_CALLOC, NULL, numOfFaces);
1726
1727                                                         if (col) {
1728                                                                 #ifdef _OPENMP
1729                                                                 #pragma omp parallel for schedule(static)
1730                                                                 #endif
1731                                                                 for (i=0; i<numOfFaces; i++) {
1732                                                                         float temp_color[3];
1733                                                                         int j=0;
1734                                                                         for (; j<((mface[i].v4)?4:3); j++) {
1735                                                                                 int index = (j==0)?mface[i].v1: (j==1)?mface[i].v2: (j==2)?mface[i].v3: mface[i].v4;
1736
1737                                                                                 col[i*4+j].a = 255;
1738
1739                                                                                 weight_to_rgb(weight[index], temp_color, temp_color+1, temp_color+2);
1740                                                                                 col[i*4+j].r = (char)(temp_color[2]*255);
1741                                                                                 col[i*4+j].g = (char)(temp_color[1]*255);
1742                                                                                 col[i*4+j].b = (char)(temp_color[0]*255);
1743                                                                         }
1744                                                                 }
1745                                                                 pmd->canvas->flags |= MOD_DPAINT_PREVIEW_READY;
1746                                                         }
1747                                                 }
1748
1749                                                 /* apply weights into a vertex group, if doesnt exists add a new layer */
1750                                                 if (defgrp_index >= 0 && !dvert && strlen(surface->output_name)>0)
1751                                                         dvert = CustomData_add_layer_named(&result->vertData, CD_MDEFORMVERT, CD_CALLOC,
1752                                                                                                                                 NULL, sData->total_points, surface->output_name);
1753                                                 if (defgrp_index >= 0 && dvert) {
1754                                                         int i;
1755                                                         for(i=0; i<sData->total_points; i++) {
1756                                                                 int j;
1757                                                                 MDeformVert *dv= &dvert[i];
1758                                                                 MDeformWeight *def_weight = NULL;
1759
1760                                                                 /* check if this vertex has a weight */
1761                                                                 for (j=0; j<dv->totweight; j++) {
1762                                                                         if (dv->dw[j].def_nr == defgrp_index) {
1763                                                                                 def_weight = &dv->dw[j];
1764                                                                                 break;
1765                                                                         }
1766                                                                 }
1767
1768                                                                 /* skip if weight value is 0 and no existing weight is found */
1769                                                                 if (!def_weight && !weight[i])
1770                                                                         continue;
1771
1772                                                                 /* if not found, add a weight for it */
1773                                                                 if (!def_weight) {
1774                                                                         MDeformWeight *newdw = MEM_callocN(sizeof(MDeformWeight)*(dv->totweight+1), 
1775                                                                                                                  "deformWeight");
1776                                                                         if(dv->dw){
1777                                                                                 memcpy(newdw, dv->dw, sizeof(MDeformWeight)*dv->totweight);
1778                                                                                 MEM_freeN(dv->dw);
1779                                                                         }
1780                                                                         dv->dw=newdw;
1781                                                                         dv->dw[dv->totweight].def_nr=defgrp_index;
1782                                                                         def_weight = &dv->dw[dv->totweight];
1783                                                                         dv->totweight++;
1784                                                                 }
1785
1786                                                                 /* set weight value */
1787                                                                 def_weight->weight = weight[i];
1788                                                         }
1789                                                 }
1790                                         }
1791                                         /* wave simulation */
1792                                         else if (surface->type == MOD_DPAINT_SURFACE_T_WAVE) {
1793                                                 MVert *mvert = result->getVertArray(result);
1794                                                 int i;
1795                                                 PaintWavePoint* wPoint = (PaintWavePoint*)sData->type_data;
1796
1797                                                 #ifdef _OPENMP
1798                                                 #pragma omp parallel for schedule(static)
1799                                                 #endif
1800                                                 for (i=0; i<sData->total_points; i++) {
1801                                                         float normal[3];
1802                                                         normal_short_to_float_v3(normal, mvert[i].no);
1803                                                         normalize_v3(normal);
1804
1805                                                         mvert[i].co[0] += normal[0]*wPoint[i].height;
1806                                                         mvert[i].co[1] += normal[1]*wPoint[i].height;
1807                                                         mvert[i].co[2] += normal[2]*wPoint[i].height;
1808                                                 }
1809                                                 CDDM_calc_normals(result);
1810                                         }
1811
1812                                         /* displace */
1813                                         dynamicPaint_applySurfaceDisplace(surface, result, 1);
1814                                 }
1815                         }
1816                 }
1817         }
1818         /* make a copy of dm to use as brush data */
1819         if (pmd->brush) {
1820                 if (pmd->brush->dm) pmd->brush->dm->release(pmd->brush->dm);
1821                 pmd->brush->dm = CDDM_copy(result);
1822         }
1823
1824         return result;
1825 }
1826
1827 /* update cache frame range */
1828 void dynamicPaint_cacheUpdateFrames(DynamicPaintSurface *surface)
1829 {
1830         if (surface->pointcache) {
1831                 surface->pointcache->startframe = surface->start_frame;
1832                 surface->pointcache->endframe = surface->end_frame;
1833         }
1834 }
1835
1836 void canvas_copyDerivedMesh(DynamicPaintCanvasSettings *canvas, DerivedMesh *dm)
1837 {
1838         if (canvas->dm) canvas->dm->release(canvas->dm);
1839         canvas->dm = CDDM_copy(dm);
1840 }
1841
1842 /*
1843 *       Updates derived mesh copy and processes dynamic paint step / caches.
1844 */
1845 static void dynamicPaint_frameUpdate(DynamicPaintModifierData *pmd, Scene *scene, Object *ob, DerivedMesh *dm)
1846 {
1847         if(pmd->canvas) {
1848                 DynamicPaintCanvasSettings *canvas = pmd->canvas;
1849                 DynamicPaintSurface *surface = canvas->surfaces.first;
1850
1851                 /* update derived mesh copy */
1852                 canvas_copyDerivedMesh(canvas, dm);
1853
1854                 /* in case image sequence baking, stop here */
1855                 if (canvas->flags & MOD_DPAINT_BAKING) return;
1856
1857                 /* loop through surfaces */
1858                 for (; surface; surface=surface->next) {
1859                         int current_frame = (int)scene->r.cfra;
1860
1861                         /* free bake data if not required anymore */
1862                         surface_freeUnusedData(surface);
1863
1864                         /* image sequences are handled by bake operator */
1865                         if (surface->format == MOD_DPAINT_SURFACE_F_IMAGESEQ) continue;
1866                         if (!(surface->flags & MOD_DPAINT_ACTIVE)) continue;
1867
1868                         /* make sure surface is valid */
1869                         dynamicPaint_checkSurfaceData(surface);
1870
1871                         /* limit frame range */
1872                         CLAMP(current_frame, surface->start_frame, surface->end_frame);
1873
1874                         if (current_frame != surface->current_frame || (int)scene->r.cfra == surface->start_frame) {
1875                                 PointCache *cache = surface->pointcache;
1876                                 PTCacheID pid;
1877                                 surface->current_frame = current_frame;
1878
1879                                 /* read point cache */
1880                                 BKE_ptcache_id_from_dynamicpaint(&pid, ob, surface);
1881                                 pid.cache->startframe = surface->start_frame;
1882                                 pid.cache->endframe = surface->end_frame;
1883                                 BKE_ptcache_id_time(&pid, scene, scene->r.cfra, NULL, NULL, NULL);
1884
1885                                 /* reset non-baked cache at first frame */
1886                                 if((int)scene->r.cfra == surface->start_frame && !(cache->flag & PTCACHE_BAKED))
1887                                 {
1888                                         cache->flag |= PTCACHE_REDO_NEEDED;
1889                                         BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
1890                                         cache->flag &= ~PTCACHE_REDO_NEEDED;
1891                                 }
1892
1893                                 /* try to read from cache */
1894                                 if(BKE_ptcache_read(&pid, (float)scene->r.cfra)) {
1895                                         BKE_ptcache_validate(cache, (int)scene->r.cfra);
1896                                 }
1897                                 /* if read failed and we're on surface range do recalculate */
1898                                 else if ((int)scene->r.cfra == current_frame
1899                                         && !(cache->flag & PTCACHE_BAKED)) {
1900                                         /* calculate surface frame */
1901                                         canvas->flags |= MOD_DPAINT_BAKING;
1902                                         dynamicPaint_calculateFrame(surface, scene, ob, current_frame);
1903                                         canvas->flags &= ~MOD_DPAINT_BAKING;
1904
1905                                         /* restore canvas derivedmesh if required */
1906                                         if (surface->type == MOD_DPAINT_SURFACE_T_DISPLACE &&
1907                                                 surface->flags & MOD_DPAINT_DISP_INCREMENTAL && surface->next)
1908                                                 canvas_copyDerivedMesh(canvas, dm);
1909
1910                                         BKE_ptcache_validate(cache, surface->current_frame);
1911                                         BKE_ptcache_write(&pid, surface->current_frame);
1912                                 }
1913                         }
1914                 }
1915         }
1916 }
1917
1918 /* Modifier call. Processes dynamic paint modifier step. */
1919 struct DerivedMesh *dynamicPaint_Modifier_do(DynamicPaintModifierData *pmd, Scene *scene, Object *ob, DerivedMesh *dm)
1920 {       
1921         /* Update canvas data for a new frame */
1922         dynamicPaint_frameUpdate(pmd, scene, ob, dm);
1923
1924         /* Return output mesh */
1925         return dynamicPaint_Modifier_apply(pmd, scene, ob, dm);
1926 }
1927
1928
1929 /***************************** Image Sequence / UV Image Surface Calls ******************************/
1930
1931 /*
1932 *       Tries to find the neighbouring pixel in given (uv space) direction.
1933 *       Result is used by effect system to move paint on the surface.
1934 *
1935 *   px,py : origin pixel x and y
1936 *       n_index : lookup direction index (use neighX,neighY to get final index)
1937 */
1938 static int dynamicPaint_findNeighbourPixel(PaintUVPoint *tempPoints, DerivedMesh *dm, char *uvname, int w, int h, int px, int py, int n_index)
1939 {
1940         /* Note: Current method only uses polygon edges to detect neighbouring pixels.
1941         *  -> It doesn't always lead to the optimum pixel but is accurate enough
1942         *  and faster/simplier than including possible face tip point links)
1943         */
1944
1945         int x,y;
1946         PaintUVPoint *tPoint = NULL;
1947         PaintUVPoint *cPoint = NULL;
1948
1949         /* shift position by given n_index */
1950         x = px + neighX[n_index];
1951         y = py + neighY[n_index];
1952
1953         if (x<0 || x>=w) return -1;
1954         if (y<0 || y>=h) return -1;
1955
1956         tPoint = &tempPoints[x+w*y];            /* UV neighbour */
1957         cPoint = &tempPoints[px+w*py];          /* Origin point */
1958
1959         /*
1960         *       Check if shifted point is on same face -> it's a correct neighbour
1961         *   (and if it isn't marked as an "edge pixel")
1962         */
1963         if ((tPoint->face_index == cPoint->face_index) && (tPoint->neighbour_pixel == -1))
1964                 return (x+w*y);
1965
1966         /*
1967         *       Even if shifted point is on another face
1968         *       -> use this point.
1969         *       
1970         *       !! Replace with "is uv faces linked" check !!
1971         *       This should work fine as long as uv island
1972         *       margin is > 1 pixel.
1973         */
1974         if ((tPoint->face_index != -1) && (tPoint->neighbour_pixel == -1)) {
1975                 return (x+w*y);
1976         }
1977
1978         /*
1979         *       If we get here, the actual neighbouring pixel
1980         *       is located on a non-linked uv face, and we have to find
1981         *       it's "real" position.
1982         *
1983         *       Simple neighbouring face finding algorithm:
1984         *       - find closest uv edge to shifted pixel and get
1985         *         the another face that shares that edge
1986         *       - find corresponding position of that new face edge
1987         *         in uv space
1988         *
1989         *       TODO: Implement something more accurate / optimized?
1990         */
1991         {
1992                 int numOfFaces = dm->getNumFaces(dm);
1993                 MFace *mface = dm->getFaceArray(dm);
1994                 MTFace *tface =  CustomData_get_layer_named(&dm->faceData, CD_MTFACE, uvname);
1995
1996                 /* Get closest edge to that subpixel on UV map  */
1997                 {
1998                         float pixel[2], dist, t_dist;
1999                         int i, uindex[2], edge1_index, edge2_index,
2000                                 e1_index, e2_index, target_face;
2001                         float closest_point[2], lambda, dir_vec[2];
2002                         int target_uv1, target_uv2, final_pixel[2], final_index;
2003
2004                         float *s_uv1, *s_uv2, *t_uv1, *t_uv2;
2005
2006                         pixel[0] = ((float)(px + neighX[n_index]) + 0.5f) / (float)w;
2007                         pixel[1] = ((float)(py + neighY[n_index]) + 0.5f) / (float)h;
2008
2009                         /* Get uv indexes for current face part */
2010                         if (cPoint->quad) {
2011                                 uindex[0] = 0; uindex[1] = 2; uindex[2] = 3;
2012                         }
2013                         else {
2014                                 uindex[0] = 0; uindex[1] = 1; uindex[2] = 2;
2015                         }
2016
2017                         /*
2018                         *       Find closest edge to that pixel
2019                         */
2020                         /* Dist to first edge   */
2021                         e1_index = cPoint->v1; e2_index = cPoint->v2; edge1_index = uindex[0]; edge2_index = uindex[1];
2022                         dist = dist_to_line_segment_v2(pixel, tface[cPoint->face_index].uv[edge1_index], tface[cPoint->face_index].uv[edge2_index]);
2023
2024                         /* Dist to second edge  */
2025                         t_dist = dist_to_line_segment_v2(pixel, tface[cPoint->face_index].uv[uindex[1]], tface[cPoint->face_index].uv[uindex[2]]);
2026                         if (t_dist < dist) {e1_index = cPoint->v2; e2_index = cPoint->v3; edge1_index = uindex[1]; edge2_index = uindex[2]; dist = t_dist;}
2027
2028                         /* Dist to third edge   */
2029                         t_dist = dist_to_line_segment_v2(pixel, tface[cPoint->face_index].uv[uindex[2]], tface[cPoint->face_index].uv[uindex[0]]);
2030                         if (t_dist < dist) {e1_index = cPoint->v3; e2_index = cPoint->v1;  edge1_index = uindex[2]; edge2_index = uindex[0]; dist = t_dist;}
2031
2032
2033                         /*
2034                         *       Now find another face that is linked to that edge
2035                         */
2036                         target_face = -1;
2037
2038                         for (i=0; i<numOfFaces; i++) {
2039                                 /*
2040                                 *       Check if both edge vertices share this face
2041                                 */
2042                                 int v4 = (mface[i].v4) ? mface[i].v4 : -1;
2043
2044                                 if ((e1_index == mface[i].v1 || e1_index == mface[i].v2 || e1_index == mface[i].v3 || e1_index == v4) &&
2045                                         (e2_index == mface[i].v1 || e2_index == mface[i].v2 || e2_index == mface[i].v3 || e2_index == v4)) {
2046                                         if (i == cPoint->face_index) continue;
2047
2048                                         target_face = i;
2049
2050                                         /*
2051                                         *       Get edge UV index
2052                                         */
2053                                         if (e1_index == mface[i].v1) target_uv1 = 0;
2054                                         else if (e1_index == mface[i].v2) target_uv1 = 1;
2055                                         else if (e1_index == mface[i].v3) target_uv1 = 2;
2056                                         else target_uv1 = 3;
2057
2058                                         if (e2_index == mface[i].v1) target_uv2 = 0;
2059                                         else if (e2_index == mface[i].v2) target_uv2 = 1;
2060                                         else if (e2_index == mface[i].v3) target_uv2 = 2;
2061                                         else target_uv2 = 3;
2062
2063                                         break;
2064                                 }
2065                         }
2066
2067                         /* If none found return -1      */
2068                         if (target_face == -1) return -1;
2069
2070                         /*
2071                         *       If target face is connected in UV space as well, just use original index
2072                         */
2073                         s_uv1 = (float *)tface[cPoint->face_index].uv[edge1_index];
2074                         s_uv2 = (float *)tface[cPoint->face_index].uv[edge2_index];
2075                         t_uv1 = (float *)tface[target_face].uv[target_uv1];
2076                         t_uv2 = (float *)tface[target_face].uv[target_uv2];
2077
2078                         //printf("connected UV : %f,%f & %f,%f - %f,%f & %f,%f\n", s_uv1[0], s_uv1[1], s_uv2[0], s_uv2[1], t_uv1[0], t_uv1[1], t_uv2[0], t_uv2[1]);
2079
2080                         if (((s_uv1[0] == t_uv1[0] && s_uv1[1] == t_uv1[1]) &&
2081                                  (s_uv2[0] == t_uv2[0] && s_uv2[1] == t_uv2[1]) ) ||
2082                                 ((s_uv2[0] == t_uv1[0] && s_uv2[1] == t_uv1[1]) &&
2083                                  (s_uv1[0] == t_uv2[0] && s_uv1[1] == t_uv2[1]) )) return ((px+neighX[n_index]) + w*(py+neighY[n_index]));
2084
2085                         /*
2086                         *       Find a point that is relatively at same edge position
2087                         *       on this other face UV
2088                         */
2089                         lambda = closest_to_line_v2(closest_point, pixel, tface[cPoint->face_index].uv[edge1_index], tface[cPoint->face_index].uv[edge2_index]);
2090                         if (lambda < 0.0f) lambda = 0.0f;
2091                         if (lambda > 1.0f) lambda = 1.0f;
2092
2093                         sub_v2_v2v2(dir_vec, tface[target_face].uv[target_uv2], tface[target_face].uv[target_uv1]);
2094
2095                         mul_v2_fl(dir_vec, lambda);
2096
2097                         copy_v2_v2(pixel, tface[target_face].uv[target_uv1]);
2098                         add_v2_v2(pixel, dir_vec);
2099                         pixel[0] = (pixel[0] * (float)w) - 0.5f;
2100                         pixel[1] = (pixel[1] * (float)h) - 0.5f;
2101
2102                         final_pixel[0] = (int)floor(pixel[0]);
2103                         final_pixel[1] = (int)floor(pixel[1]);
2104
2105                         /* If current pixel uv is outside of texture    */
2106                         if (final_pixel[0] < 0 || final_pixel[0] >= w) return -1;
2107                         if (final_pixel[1] < 0 || final_pixel[1] >= h) return -1;
2108
2109                         final_index = final_pixel[0] + w * final_pixel[1];
2110
2111                         /* If we ended up to our origin point ( mesh has smaller than pixel sized faces)        */
2112                         if (final_index == (px+w*py)) return -1;
2113                         /* If found pixel still lies on wrong face ( mesh has smaller than pixel sized faces)   */
2114                         if (tempPoints[final_index].face_index != target_face) return -1;
2115
2116                         /*
2117                         *       If final point is an "edge pixel", use it's "real" neighbour instead
2118                         */
2119                         if (tempPoints[final_index].neighbour_pixel != -1) final_index = cPoint->neighbour_pixel;
2120
2121                         return final_index;
2122                 }
2123         }
2124 }
2125
2126 /*
2127 *       Create a surface for uv image sequence format
2128 */
2129 static int dynamicPaint_createUVSurface(DynamicPaintSurface *surface)
2130 {
2131         /* Antialias jitter point relative coords       */
2132         float jitter5sample[10] =  {0.0f, 0.0f,
2133                                                         -0.2f, -0.4f,
2134                                                         0.2f, 0.4f,
2135                                                         0.4f, -0.2f,
2136                                                         -0.4f, 0.3f};
2137         int ty;
2138         int w,h;
2139         int numOfFaces;
2140         char uvname[32];
2141         int active_points = 0;
2142         int error = 0;
2143
2144         PaintSurfaceData *sData;
2145         DynamicPaintCanvasSettings *canvas = surface->canvas;
2146         DerivedMesh *dm = canvas->dm;
2147
2148         PaintUVPoint *tempPoints = NULL;
2149         Vec3f *tempWeights = NULL;
2150         MVert *mvert = NULL;
2151         MFace *mface = NULL;
2152         MTFace *tface = NULL;
2153         Bounds2D *faceBB = NULL;
2154         int *final_index;
2155         int aa_samples;
2156
2157         if (!dm) return printError(canvas, "Canvas mesh not updated.");
2158         if (surface->format != MOD_DPAINT_SURFACE_F_IMAGESEQ) return printError(canvas, "Can't bake non-\"image sequence\" formats.");
2159
2160         numOfFaces = dm->getNumFaces(dm);
2161         mvert = dm->getVertArray(dm);
2162         mface = dm->getFaceArray(dm);
2163
2164         /* get uv layer */
2165         validate_layer_name(&dm->faceData, CD_MTFACE, surface->uvlayer_name, uvname);
2166         tface = CustomData_get_layer_named(&dm->faceData, CD_MTFACE, uvname);
2167
2168         /* Check for validity   */
2169         if (!tface) return printError(canvas, "No UV data on canvas.");
2170         if (surface->image_resolution < 16 || surface->image_resolution > 8096) return printError(canvas, "Invalid resolution.");
2171
2172         w = h = surface->image_resolution;
2173
2174         /*
2175         *       Start generating the surface
2176         */
2177         printf("DynamicPaint: Preparing UV surface of %ix%i pixels and %i faces.\n", w, h, numOfFaces);
2178
2179         /* Init data struct */
2180         if (surface->data) dynamicPaint_freeSurfaceData(surface);
2181         sData = surface->data = MEM_callocN(sizeof(PaintSurfaceData), "PaintSurfaceData");
2182         if (!surface->data) return printError(canvas, "Not enough free memory.");
2183
2184         aa_samples = (surface->flags & MOD_DPAINT_ANTIALIAS) ? 5 : 1;
2185         tempPoints = (struct PaintUVPoint *) MEM_callocN(w*h*sizeof(struct PaintUVPoint), "Temp PaintUVPoint");
2186         if (!tempPoints) error=1;
2187
2188         final_index = (int *) MEM_callocN(w*h*sizeof(int), "Temp UV Final Indexes");
2189         if (!final_index) error=1;
2190
2191         if (!error) {
2192                 tempWeights = (struct Vec3f *) MEM_mallocN(w*h*aa_samples*sizeof(struct Vec3f), "Temp bWeights");
2193                 if (!tempWeights) error=1;
2194         }
2195
2196         /*
2197         *       Generate a temporary bounding box array for UV faces to optimize
2198         *       the pixel-inside-a-face search.
2199         */
2200         if (!error) {
2201                 faceBB = (struct Bounds2D *) MEM_mallocN(numOfFaces*sizeof(struct Bounds2D), "MPCanvasFaceBB");
2202                 if (!faceBB) error=1;
2203         }
2204
2205         if (!error)
2206         for (ty=0; ty<numOfFaces; ty++) {
2207                 int numOfVert = (mface[ty].v4) ? 4 : 3;
2208                 int i;
2209
2210                 VECCOPY2D(faceBB[ty].min, tface[ty].uv[0]);
2211                 VECCOPY2D(faceBB[ty].max, tface[ty].uv[0]);
2212
2213                 for (i = 1; i<numOfVert; i++) {
2214                         if (tface[ty].uv[i][0] < faceBB[ty].min[0]) faceBB[ty].min[0] = tface[ty].uv[i][0];
2215                         if (tface[ty].uv[i][1] < faceBB[ty].min[1]) faceBB[ty].min[1] = tface[ty].uv[i][1];
2216                         if (tface[ty].uv[i][0] > faceBB[ty].max[0]) faceBB[ty].max[0] = tface[ty].uv[i][0];
2217                         if (tface[ty].uv[i][1] > faceBB[ty].max[1]) faceBB[ty].max[1] = tface[ty].uv[i][1];
2218
2219                 }
2220         }
2221
2222         /*
2223         *       Loop through every pixel and check
2224         *       if pixel is uv-mapped on a canvas face.
2225         */
2226         if (!error) {
2227                 #ifdef _OPENMP
2228                 #pragma omp parallel for schedule(static)
2229                 #endif
2230                 for (ty = 0; ty < h; ty++)
2231                 {
2232                         int tx;
2233                         for (tx = 0; tx < w; tx++)
2234                         {
2235                                 int i, sample;
2236                                 int index = tx+w*ty;
2237                                 PaintUVPoint *tPoint = (&tempPoints[index]);
2238
2239                                 short isInside = 0;     /* if point is inside a uv face */
2240
2241                                 float d1[2], d2[2], d3[2], point[5][2];
2242                                 float dot00,dot01,dot02,dot11,dot12, invDenom, u,v;
2243
2244                                 /* Init per pixel settings */
2245                                 tPoint->face_index = -1;
2246                                 tPoint->neighbour_pixel = -1;
2247                                 tPoint->pixel_index = index;
2248
2249                                 /* Actual pixel center, used when collision is found    */
2250                                 point[0][0] = ((float)tx + 0.5f) / w;
2251                                 point[0][1] = ((float)ty + 0.5f) / h;
2252
2253                                 /*
2254                                 * A pixel middle sample isn't enough to find very narrow polygons
2255                                 * So using 4 samples of each corner too
2256                                 */
2257                                 point[1][0] = ((float)tx) / w;
2258                                 point[1][1] = ((float)ty) / h;
2259
2260                                 point[2][0] = ((float)tx+1) / w;
2261                                 point[2][1] = ((float)ty) / h;
2262
2263                                 point[3][0] = ((float)tx) / w;
2264                                 point[3][1] = ((float)ty+1) / h;
2265
2266                                 point[4][0] = ((float)tx+1) / w;
2267                                 point[4][1] = ((float)ty+1) / h;
2268
2269
2270                                 /* Loop through samples, starting from middle point     */
2271                                 for (sample=0; sample<5; sample++) {
2272                                         
2273                                         /* Loop through every face in the mesh  */
2274                                         for (i=0; i<numOfFaces; i++) {
2275
2276                                                 /* Check uv bb  */
2277                                                 if (faceBB[i].min[0] > (point[sample][0])) continue;
2278                                                 if (faceBB[i].min[1] > (point[sample][1])) continue;
2279                                                 if (faceBB[i].max[0] < (point[sample][0])) continue;
2280                                                 if (faceBB[i].max[1] < (point[sample][1])) continue;
2281
2282                                                 /*  Calculate point inside a triangle check
2283                                                 *       for uv0,1,2 */
2284                                                 VECSUB2D(d1,  tface[i].uv[2], tface[i].uv[0]);  // uv2 - uv0
2285                                                 VECSUB2D(d2,  tface[i].uv[1], tface[i].uv[0]);  // uv1 - uv0
2286                                                 VECSUB2D(d3,  point[sample], tface[i].uv[0]);   // point - uv0
2287
2288                                                 dot00 = d1[0]*d1[0] + d1[1]*d1[1];
2289                                                 dot01 = d1[0]*d2[0] + d1[1]*d2[1];
2290                                                 dot02 = d1[0]*d3[0] + d1[1]*d3[1];
2291                                                 dot11 = d2[0]*d2[0] + d2[1]*d2[1];
2292                                                 dot12 = d2[0]*d3[0] + d2[1]*d3[1];
2293
2294                                                 invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
2295                                                 u = (dot11 * dot02 - dot01 * dot12) * invDenom;
2296                                                 v = (dot00 * dot12 - dot01 * dot02) * invDenom;
2297
2298                                                 if ((u > 0) && (v > 0) && (u + v < 1)) {isInside=1;} /* is inside a triangle */
2299
2300                                                 /*  If collision wasn't found but the face is a quad
2301                                                 *       do another check for the second half */
2302                                                 if ((!isInside) && mface[i].v4)
2303                                                 {
2304
2305                                                         /* change d2 to test the other half     */
2306                                                         VECSUB2D(d2,  tface[i].uv[3], tface[i].uv[0]);  // uv3 - uv0
2307
2308                                                         /* test again   */
2309                                                         dot00 = d1[0]*d1[0] + d1[1]*d1[1];
2310                                                         dot01 = d1[0]*d2[0] + d1[1]*d2[1];
2311                                                         dot02 = d1[0]*d3[0] + d1[1]*d3[1];
2312                                                         dot11 = d2[0]*d2[0] + d2[1]*d2[1];
2313                                                         dot12 = d2[0]*d3[0] + d2[1]*d3[1];
2314
2315                                                         invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
2316                                                         u = (dot11 * dot02 - dot01 * dot12) * invDenom;
2317                                                         v = (dot00 * dot12 - dot01 * dot02) * invDenom;
2318
2319                                                         if ((u > 0) && (v > 0) && (u + v < 1)) {isInside=2;} /* is inside the second half of the quad */
2320
2321                                                 }
2322
2323                                                 /*
2324                                                 *       If point was inside the face
2325                                                 */
2326                                                 if (isInside != 0) {
2327
2328                                                         float uv1co[2], uv2co[2], uv3co[2], uv[2];
2329                                                         int j;
2330
2331                                                         /* Get triagnle uvs     */
2332                                                         if (isInside==1) {
2333                                                                 VECCOPY2D(uv1co, tface[i].uv[0]);
2334                                                                 VECCOPY2D(uv2co, tface[i].uv[1]);
2335                                                                 VECCOPY2D(uv3co, tface[i].uv[2]);
2336                                                         }
2337                                                         else {
2338                                                                 VECCOPY2D(uv1co, tface[i].uv[0]);
2339                                                                 VECCOPY2D(uv2co, tface[i].uv[2]);
2340                                                                 VECCOPY2D(uv3co, tface[i].uv[3]);
2341                                                         }
2342
2343                                                         /* Add b-weights per anti-aliasing sample       */
2344                                                         for (j=0; j<aa_samples; j++) {
2345                                                                 uv[0] = point[0][0] + jitter5sample[j*2] / w;
2346                                                                 uv[1] = point[0][1] + jitter5sample[j*2+1] / h;
2347
2348                                                                 barycentric_weights_v2(uv1co, uv2co, uv3co, uv, tempWeights[index*aa_samples+j].v);
2349                                                         }
2350
2351                                                         /* Set surface point face values        */
2352                                                         tPoint->face_index = i;                                                 /* face index */
2353                                                         tPoint->quad = (isInside == 2) ? 1 : 0;         /* quad or tri part*/
2354
2355                                                         /* save vertex indexes  */
2356                                                         tPoint->v1 = (isInside == 2) ? mface[i].v1 : mface[i].v1;
2357                                                         tPoint->v2 = (isInside == 2) ? mface[i].v3 : mface[i].v2;
2358                                                         tPoint->v3 = (isInside == 2) ? mface[i].v4 : mface[i].v3;
2359                                                         
2360                                                         sample = 5;     /* make sure we exit sample loop as well */
2361                                                         break;
2362                                                 }
2363                                         }
2364                                 } /* sample loop */
2365                         }
2366                 }
2367
2368                 /*
2369                 *       Now loop through every pixel that was left without index
2370                 *       and find if they have neighbouring pixels that have an index.
2371                 *       If so use that polygon as pixel surface.
2372                 *       (To avoid seams on uv island edges)
2373                 */
2374                 #ifdef _OPENMP
2375                 #pragma omp parallel for schedule(static)
2376                 #endif
2377                 for (ty = 0; ty < h; ty++)
2378                 {
2379                         int tx;
2380                         for (tx = 0; tx < w; tx++)
2381                         {
2382                                 int index = tx+w*ty;
2383                                 PaintUVPoint *tPoint = (&tempPoints[index]);
2384
2385                                 /* If point isnt't on canvas mesh       */
2386                                 if (tPoint->face_index == -1) {
2387                                         int u_min, u_max, v_min, v_max;
2388                                         int u,v, ind;
2389                                         float point[2];
2390
2391                                         /* get loop area        */
2392                                         u_min = (tx > 0) ? -1 : 0;
2393                                         u_max = (tx < (w-1)) ? 1 : 0;
2394                                         v_min = (ty > 0) ? -1 : 0;
2395                                         v_max = (ty < (h-1)) ? 1 : 0;
2396
2397                                         point[0] = ((float)tx + 0.5f) / w;
2398                                         point[1] = ((float)ty + 0.5f) / h;
2399
2400                                         /* search through defined area for neighbour    */
2401                                         for (u=u_min; u<=u_max; u++)
2402                                                 for (v=v_min; v<=v_max; v++) {
2403                                                         /* if not this pixel itself     */
2404                                                         if (u!=0 || v!=0) {
2405                                                                 ind = (tx+u)+w*(ty+v);
2406
2407                                                                 /* if neighbour has index       */
2408                                                                 if (tempPoints[ind].face_index != -1) {
2409
2410                                                                         float uv1co[2], uv2co[2], uv3co[2], uv[2];
2411                                                                         int i = tempPoints[ind].face_index, j;
2412
2413                                                                         /* Now calculate pixel data for this pixel as it was on polygon surface */
2414                                                                         if (!tempPoints[ind].quad) {
2415                                                                                 VECCOPY2D(uv1co, tface[i].uv[0]);
2416                                                                                 VECCOPY2D(uv2co, tface[i].uv[1]);
2417                                                                                 VECCOPY2D(uv3co, tface[i].uv[2]);
2418                                                                         }
2419                                                                         else {
2420                                                                                 VECCOPY2D(uv1co, tface[i].uv[0]);
2421                                                                                 VECCOPY2D(uv2co, tface[i].uv[2]);
2422                                                                                 VECCOPY2D(uv3co, tface[i].uv[3]);
2423                                                                         }
2424
2425                                                                         /* Add b-weights per anti-aliasing sample       */
2426                                                                         for (j=0; j<aa_samples; j++) {
2427
2428                                                                                 uv[0] = point[0] + jitter5sample[j*2] / w;
2429                                                                                 uv[1] = point[1] + jitter5sample[j*2+1] / h;
2430                                                                                 barycentric_weights_v2(uv1co, uv2co, uv3co, uv, tempWeights[index*aa_samples+j].v);
2431                                                                         }
2432
2433                                                                         /* Set values   */
2434                                                                         tPoint->neighbour_pixel = ind;                          // face index
2435                                                                         tPoint->quad = tempPoints[ind].quad;            // quad or tri
2436
2437                                                                         /* save vertex indexes  */
2438                                                                         tPoint->v1 = (tPoint->quad) ? mface[i].v1 : mface[i].v1;
2439                                                                         tPoint->v2 = (tPoint->quad) ? mface[i].v3 : mface[i].v2;
2440                                                                         tPoint->v3 = (tPoint->quad) ? mface[i].v4 : mface[i].v3;
2441
2442                                                                         u = u_max + 1;  /* make sure we exit outer loop as well */
2443                                                                         break;
2444                                                                 }
2445                                                 }
2446                                         }
2447                                 }
2448                         }
2449                 }
2450
2451                 /*
2452                 *       When base loop is over convert found neighbour indexes to real ones
2453                 *       Also count the final number of active surface points
2454                 */
2455                 for (ty = 0; ty < h; ty++)
2456                 {
2457                         int tx;
2458                         for (tx = 0; tx < w; tx++)
2459                         {
2460                                 int index = tx+w*ty;
2461                                 PaintUVPoint *tPoint = (&tempPoints[index]);
2462
2463                                 if (tPoint->face_index == -1 && tPoint->neighbour_pixel != -1) tPoint->face_index = tempPoints[tPoint->neighbour_pixel].face_index;
2464                                 if (tPoint->face_index != -1) active_points++;
2465                         }
2466                 }
2467
2468                 /*      If any effect enabled, create surface effect / wet layer
2469                 *       neighbour lists. Processes possibly moving data. */
2470                 if (surface_usesAdjData(surface)) {
2471
2472                         int i, cursor=0;
2473
2474                         /* Create a temporary array of final indexes (before unassigned
2475                         *  pixels have been dropped) */
2476                         for (i=0; i<w*h; i++) {
2477                                 if (tempPoints[i].face_index != -1) {
2478                                         final_index[i] = cursor;
2479                                         cursor++;
2480                                 }
2481                         }
2482                         /* allocate memory */
2483                         sData->total_points = w*h;
2484                         dynamicPaint_initAdjacencyData(surface, 0);
2485
2486                         if (sData->adj_data) {
2487                                 PaintAdjData *ed = sData->adj_data;
2488                                 unsigned int n_pos = 0;
2489                                 //#pragma omp parallel for schedule(static)
2490                                 for (ty = 0; ty < h; ty++)
2491                                 {
2492                                         int tx;
2493                                         for (tx = 0; tx < w; tx++)
2494                                         {
2495                                                 int i, index = tx+w*ty;
2496
2497                                                 if (tempPoints[index].face_index != -1) {
2498                                                         ed->n_index[final_index[index]] = n_pos;
2499                                                         ed->n_num[final_index[index]] = 0;
2500
2501                                                         for (i=0; i<8; i++) {
2502
2503                                                                 /* Try to find a neighbouring pixel in defined direction
2504                                                                 *  If not found, -1 is returned */
2505                                                                 int n_target = dynamicPaint_findNeighbourPixel(tempPoints, dm, uvname, w, h, tx, ty, i);
2506
2507                                                                 if (n_target != -1) {
2508                                                                         ed->n_target[n_pos] = final_index[n_target];
2509                                                                         ed->n_num[final_index[index]]++;
2510                                                                         n_pos++;
2511                                                                 }
2512                                                         }
2513                                                 }
2514                                         }
2515                                 }
2516                         }
2517                 }
2518
2519                 /* Create final surface data without inactive points */
2520                 {
2521                         ImgSeqFormatData *f_data = MEM_callocN(sizeof(struct ImgSeqFormatData), "ImgSeqFormatData");
2522                         if (f_data) {
2523                                 f_data->uv_p = MEM_callocN(active_points*sizeof(struct PaintUVPoint), "PaintUVPoint");
2524                                 f_data->barycentricWeights = MEM_callocN(active_points*aa_samples*sizeof(struct Vec3f), "PaintUVPoint");
2525
2526                                 if (!f_data->uv_p || !f_data->barycentricWeights) error=1;
2527                         }
2528                         else error=1;
2529
2530                         sData->total_points = active_points;
2531                         
2532                         /* in case of allocation error, free everything */
2533                         if (error) {
2534                                 if (f_data) {
2535                                         if (f_data->uv_p) MEM_freeN(f_data->uv_p);
2536                                         if (f_data->barycentricWeights) MEM_freeN(f_data->barycentricWeights);
2537                                         MEM_freeN(f_data);
2538                                 }
2539                         }
2540                         else {
2541                                 int index, cursor = 0;
2542                                 sData->total_points = active_points;
2543                                 sData->format_data = f_data;
2544
2545                                 for(index = 0; index < (w*h); index++) {
2546                                         if (tempPoints[index].face_index != -1) {
2547                                                 memcpy(&f_data->uv_p[cursor], &tempPoints[index], sizeof(PaintUVPoint));
2548                                                 memcpy(&f_data->barycentricWeights[cursor*aa_samples], &tempWeights[index*aa_samples], sizeof(Vec3f)*aa_samples);
2549                                                 cursor++;
2550                                         }
2551                                 }
2552                         }
2553                 }
2554         }
2555         if (error==1) printError(canvas, "Not enough free memory.");
2556
2557         if (faceBB) MEM_freeN(faceBB);
2558         if (tempPoints) MEM_freeN(tempPoints);
2559         if (tempWeights) MEM_freeN(tempWeights);
2560         if (final_index) MEM_freeN(final_index);
2561
2562         /* Init surface type data */
2563         if (!error) {
2564                 dynamicPaint_allocateSurfaceType(surface);
2565
2566 #if 0
2567                 /*  -----------------------------------------------------------------
2568                 *       For debug, output pixel statuses to the color map
2569                 *       -----------------------------------------------------------------*/
2570                 #ifdef _OPENMP
2571                 #pragma omp parallel for schedule(static)
2572                 #endif
2573                 for (index = 0; index < sData->total_points; index++)
2574                 {
2575                         ImgSeqFormatData *f_data = (ImgSeqFormatData*)sData->format_data;
2576                         PaintUVPoint *uvPoint = &((PaintUVPoint*)f_data->uv_p)[index];
2577                         PaintPoint *pPoint = &((PaintPoint*)sData->type_data)[index];
2578                         pPoint->alpha=1.0f;
2579
2580                         /* Every pixel that is assigned as "edge pixel" gets blue color */
2581                         if (uvPoint->neighbour_pixel != -1) pPoint->color[2] = 1.0f;
2582                         /* and every pixel that finally got an polygon gets red color   */
2583                         if (uvPoint->face_index != -1) pPoint->color[0] = 1.0f;
2584                         /* green color shows pixel face index hash      */
2585                         if (uvPoint->face_index != -1) pPoint->color[1] = (float)(uvPoint->face_index % 255)/256.0f;
2586                 }
2587
2588 #endif
2589                 dynamicPaint_setInitialColor(surface);
2590         }
2591
2592         return (error == 0);
2593 }
2594
2595 #define DPOUTPUT_PAINT 0
2596 #define DPOUTPUT_WET 1
2597 #define DPOUTPUT_DISPLACE 2
2598 #define DPOUTPUT_WAVES 3
2599
2600 /*
2601 *       Outputs an image file from uv surface data.
2602 */
2603 void dynamicPaint_outputImage(DynamicPaintSurface *surface, char* filename, short format, short type)
2604 {
2605         int index;
2606         ImBuf* mhImgB = NULL;
2607         PaintSurfaceData *sData = surface->data;
2608         ImgSeqFormatData *f_data = (ImgSeqFormatData*)sData->format_data;
2609         char output_file[250];
2610
2611         if (sData == NULL || sData->type_data == NULL) {printError(surface->canvas, "Image save failed: Invalid surface.");return;}
2612
2613         if (format == DPOUTPUT_JPEG) sprintf(output_file,"%s.jpg",filename);
2614         else if (format == DPOUTPUT_OPENEXR) sprintf(output_file,"%s.exr",filename);
2615         else sprintf(output_file,"%s.png",filename);
2616
2617         /* Validate output file path    */
2618         BLI_path_abs(output_file, G.main->name);
2619         BLI_make_existing_file(output_file);
2620
2621         /* Init image buffer    */
2622         mhImgB = IMB_allocImBuf(surface->image_resolution, surface->image_resolution, 32, IB_rectfloat);
2623         if (mhImgB == NULL) {printError(surface->canvas, "Image save failed: Not enough free memory.");return;}
2624
2625         #ifdef _OPENMP
2626         #pragma omp parallel for schedule(static)
2627         #endif
2628         for (index = 0; index < sData->total_points; index++)
2629         {
2630                 int pos=f_data->uv_p[index].pixel_index*4;      /* image buffer position */
2631
2632                 /* Set values of preferred type */
2633                 if (type == DPOUTPUT_WET) {
2634                         PaintPoint *point = &((PaintPoint*)sData->type_data)[index];
2635                         float value = (point->wetness > 1.0f) ? 1.0f : point->wetness;
2636
2637                         mhImgB->rect_float[pos]=value;
2638                         mhImgB->rect_float[pos+1]=value;
2639                         mhImgB->rect_float[pos+2]=value;
2640                         mhImgB->rect_float[pos+3]=1.0f;
2641                 }
2642                 else if (type == DPOUTPUT_PAINT) {
2643                         PaintPoint *point = &((PaintPoint*)sData->type_data)[index];
2644
2645                         mhImgB->rect_float[pos]   = point->color[0];
2646                         mhImgB->rect_float[pos+1] = point->color[1];
2647                         mhImgB->rect_float[pos+2] = point->color[2];
2648                         /* mix wet layer */
2649                         if (point->e_alpha) mixColors(&mhImgB->rect_float[pos], point->alpha, point->e_color, point->e_alpha);
2650
2651                         /* use highest alpha    */
2652                         mhImgB->rect_float[pos+3] = (point->e_alpha > point->alpha) ? point->e_alpha : point->alpha;
2653
2654                         /* Multiply color by alpha if enabled   */
2655                         if (surface->flags & MOD_DPAINT_MULALPHA) {
2656                                 mhImgB->rect_float[pos]   *= mhImgB->rect_float[pos+3];
2657                                 mhImgB->rect_float[pos+1] *= mhImgB->rect_float[pos+3];
2658                                 mhImgB->rect_float[pos+2] *= mhImgB->rect_float[pos+3];
2659                         }
2660                 }
2661                 else if (type == DPOUTPUT_DISPLACE) {
2662                         float depth = ((float*)sData->type_data)[index];
2663
2664                         if (surface->disp_type == MOD_DPAINT_DISP_DISPLACE) {
2665                                 if (surface->disp_clamp)
2666                                         depth /= surface->disp_clamp*2.0f;
2667                                 depth = (0.5f - depth);
2668                                 CLAMP(depth, 0.0f, 1.0f);
2669                         }
2670
2671                         mhImgB->rect_float[pos]=depth;
2672                         mhImgB->rect_float[pos+1]=depth;
2673                         mhImgB->rect_float[pos+2]=depth;
2674                         mhImgB->rect_float[pos+3]=1.0f;
2675                 }
2676                 else if (type == DPOUTPUT_WAVES) {
2677                         PaintWavePoint *wPoint = &((PaintWavePoint*)sData->type_data)[index];
2678                         float depth = wPoint->height/2.0f+0.5f;
2679
2680                         mhImgB->rect_float[pos]=depth;
2681                         mhImgB->rect_float[pos+1]=depth;
2682                         mhImgB->rect_float[pos+2]=depth;
2683                         mhImgB->rect_float[pos+3]=1.0f;
2684                 }
2685         }
2686
2687         /* Save image buffer    */
2688         if (format == DPOUTPUT_JPEG) {  /* JPEG */
2689                 mhImgB->ftype= JPG|95;
2690                 IMB_rect_from_float(mhImgB);
2691                 imb_savejpeg(mhImgB, output_file, IB_rectfloat);
2692         }
2693 #ifdef WITH_OPENEXR
2694         else if (format == DPOUTPUT_OPENEXR) {  /* OpenEXR 32-bit float */
2695                 mhImgB->ftype = OPENEXR | OPENEXR_COMPRESS;
2696                 IMB_rect_from_float(mhImgB);
2697                 imb_save_openexr(mhImgB, output_file, IB_rectfloat);
2698         }
2699 #endif
2700         else {  /* DPOUTPUT_PNG */
2701                 mhImgB->ftype= PNG|95;
2702                 IMB_rect_from_float(mhImgB);
2703                 imb_savepng(mhImgB, output_file, IB_rectfloat);
2704         }
2705
2706         IMB_freeImBuf(mhImgB);
2707 }
2708
2709
2710 /***************************** Material / Texture Sampling ******************************/
2711
2712 /* stores a copy of required materials to allow doing adjustments
2713 *  without interfering the render/preview */
2714 typedef struct BrushMaterials {
2715         Material *mat;
2716         Material **ob_mats;
2717         int tot;
2718 } BrushMaterials;
2719
2720 /* A modified part of shadeinput.c -> shade_input_set_uv()
2721 *  Used for sampling UV mapped texture color */
2722 static void textured_face_generate_uv(float *uv, float *normal, float *hit, float *v1, float *v2, float *v3)
2723 {
2724
2725         float detsh, t00, t10, t01, t11, xn, yn, zn;
2726         int axis1, axis2;
2727
2728         /* find most stable axis to project */
2729         xn= fabs(normal[0]);
2730         yn= fabs(normal[1]);
2731         zn= fabs(normal[2]);
2732
2733         if(zn>=xn && zn>=yn) { axis1= 0; axis2= 1; }
2734         else if(yn>=xn && yn>=zn) { axis1= 0; axis2= 2; }
2735         else { axis1= 1; axis2= 2; }
2736
2737         /* compute u,v and derivatives */
2738         t00= v3[axis1]-v1[axis1]; t01= v3[axis2]-v1[axis2];
2739         t10= v3[axis1]-v2[axis1]; t11= v3[axis2]-v2[axis2];
2740
2741         detsh= 1.0f/(t00*t11-t10*t01);
2742         t00*= detsh; t01*=detsh; 
2743         t10*=detsh; t11*=detsh;
2744
2745         uv[0] = (hit[axis1]-v3[axis1])*t11-(hit[axis2]-v3[axis2])*t10;
2746         uv[1] = (hit[axis2]-v3[axis2])*t00-(hit[axis1]-v3[axis1])*t01;
2747
2748         /* u and v are in range -1 to 0, we allow a little bit extra but not too much, screws up speedvectors */
2749         CLAMP(uv[0], -2.0f, 1.0f);
2750         CLAMP(uv[1], -2.0f, 1.0f);
2751 }
2752
2753 /* a modified part of shadeinput.c -> shade_input_set_shade_texco()
2754 *  Used for sampling UV mapped texture color */
2755 static void textured_face_get_uv(float *uv_co, float *normal, float *uv, int faceIndex, short quad, MTFace *tface)
2756 {
2757         float *uv1, *uv2, *uv3;
2758         float l;
2759
2760         l= 1.0f+uv[0]+uv[1];
2761                 
2762         uv1= tface[faceIndex].uv[0];
2763         uv2= (quad) ? tface[faceIndex].uv[2] : tface[faceIndex].uv[1];
2764         uv3= (quad) ? tface[faceIndex].uv[3] : tface[faceIndex].uv[2];
2765                                 
2766         uv_co[0]= -1.0f + 2.0f*(l*uv3[0]-uv[0]*uv1[0]-uv[1]*uv2[0]);
2767         uv_co[1]= -1.0f + 2.0f*(l*uv3[1]-uv[0]*uv1[1]-uv[1]*uv2[1]);
2768         uv_co[2]= 0.0f; /* texture.c assumes there are 3 coords */
2769 }
2770
2771 /*
2772 *       Generate an updated copy of material to use for brush sampling.
2773 *       Updates animated textures and calculates inverse matrices