header cleanup and typo's
[blender.git] / source / blender / nodes / composite / node_composite_util.c
1 /*
2  * $Id$
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version. 
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19  *
20  * The Original Code is Copyright (C) 2006 Blender Foundation.
21  * All rights reserved.
22  *
23  * The Original Code is: all of this file.
24  *
25  * Contributor(s): none yet.
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29
30 /** \file blender/nodes/composite/node_composite_util.c
31  *  \ingroup nodes
32  */
33
34
35 #include "node_composite_util.h"
36
37 CompBuf *alloc_compbuf(int sizex, int sizey, int type, int alloc)
38 {
39         CompBuf *cbuf= MEM_callocN(sizeof(CompBuf), "compbuf");
40         
41         cbuf->x= sizex;
42         cbuf->y= sizey;
43         cbuf->xrad= sizex/2;
44         cbuf->yrad= sizey/2;
45         
46         cbuf->type= type;
47         if(alloc) {
48                 if(cbuf->type==CB_RGBA)
49                         cbuf->rect= MEM_mapallocN(4*sizeof(float)*sizex*sizey, "compbuf RGBA rect");
50                 else if(cbuf->type==CB_VEC3)
51                         cbuf->rect= MEM_mapallocN(3*sizeof(float)*sizex*sizey, "compbuf Vector3 rect");
52                 else if(cbuf->type==CB_VEC2)
53                         cbuf->rect= MEM_mapallocN(2*sizeof(float)*sizex*sizey, "compbuf Vector2 rect");
54                 else
55                         cbuf->rect= MEM_mapallocN(sizeof(float)*sizex*sizey, "compbuf Fac rect");
56                 cbuf->malloc= 1;
57         }
58         cbuf->disprect.xmin= 0;
59         cbuf->disprect.ymin= 0;
60         cbuf->disprect.xmax= sizex;
61         cbuf->disprect.ymax= sizey;
62         
63         return cbuf;
64 }
65
66 CompBuf *dupalloc_compbuf(CompBuf *cbuf)
67 {
68         CompBuf *dupbuf= alloc_compbuf(cbuf->x, cbuf->y, cbuf->type, 1);
69         if(dupbuf) {
70                 memcpy(dupbuf->rect, cbuf->rect, cbuf->type*sizeof(float)*cbuf->x*cbuf->y);
71         
72                 dupbuf->xof= cbuf->xof;
73                 dupbuf->yof= cbuf->yof;
74         }       
75         return dupbuf;
76 }
77
78 /* instead of reference counting, we create a list */
79 CompBuf *pass_on_compbuf(CompBuf *cbuf)
80 {
81         CompBuf *dupbuf= (cbuf)? alloc_compbuf(cbuf->x, cbuf->y, cbuf->type, 0): NULL;
82         CompBuf *lastbuf;
83         
84         if(dupbuf) {
85                 dupbuf->rect= cbuf->rect;
86                 dupbuf->xof= cbuf->xof;
87                 dupbuf->yof= cbuf->yof;
88                 dupbuf->malloc= 0;
89                 
90                 /* get last buffer in list, and append dupbuf */
91                 for(lastbuf= cbuf; lastbuf; lastbuf= lastbuf->next)
92                         if(lastbuf->next==NULL)
93                                 break;
94                 lastbuf->next= dupbuf;
95                 dupbuf->prev= lastbuf;
96         }       
97         return dupbuf;
98 }
99
100
101 void free_compbuf(CompBuf *cbuf)
102 {
103         /* check referencing, then remove from list and set malloc tag */
104         if(cbuf->prev || cbuf->next) {
105                 if(cbuf->prev)
106                         cbuf->prev->next= cbuf->next;
107                 if(cbuf->next)
108                         cbuf->next->prev= cbuf->prev;
109                 if(cbuf->malloc) {
110                         if(cbuf->prev)
111                                 cbuf->prev->malloc= 1;
112                         else
113                                 cbuf->next->malloc= 1;
114                         cbuf->malloc= 0;
115                 }
116         }
117         
118         if(cbuf->malloc && cbuf->rect)
119                 MEM_freeN(cbuf->rect);
120
121         MEM_freeN(cbuf);
122 }
123
124 void print_compbuf(char *str, CompBuf *cbuf)
125 {
126         printf("Compbuf %s %d %d %p\n", str, cbuf->x, cbuf->y, (void *)cbuf->rect);
127         
128 }
129
130 void compbuf_set_node(CompBuf *cbuf, bNode *node)
131 {
132         if (cbuf) cbuf->node = node;
133 }
134
135 /* used for disabling node  (similar code in node_draw.c for disable line and node_edit for untangling nodes) */
136 void node_compo_pass_on(bNode *node, bNodeStack **nsin, bNodeStack **nsout)
137 {
138         CompBuf *valbuf= NULL, *colbuf= NULL, *vecbuf= NULL;
139         bNodeSocket *sock;
140         int a;
141         
142         /* connect the first value buffer in with first value out */
143         /* connect the first RGBA buffer in with first RGBA out */
144         
145         /* test the inputs */
146         for(a=0, sock= node->inputs.first; sock; sock= sock->next, a++) {
147                 if(nsin[a]->data) {
148                         CompBuf *cbuf= nsin[a]->data;
149                         if(cbuf->type==1 && valbuf==NULL) valbuf= cbuf;
150                         if(cbuf->type==3 && vecbuf==NULL) vecbuf= cbuf;
151                         if(cbuf->type==4 && colbuf==NULL) colbuf= cbuf;
152                 }
153         }
154         
155         /* outputs */
156         if(valbuf || colbuf || vecbuf) {
157                 for(a=0, sock= node->outputs.first; sock; sock= sock->next, a++) {
158                         if(nsout[a]->hasoutput) {
159                                 if(sock->type==SOCK_FLOAT && valbuf) {
160                                         nsout[a]->data= pass_on_compbuf(valbuf);
161                                         valbuf= NULL;
162                                 }
163                                 if(sock->type==SOCK_VECTOR && vecbuf) {
164                                         nsout[a]->data= pass_on_compbuf(vecbuf);
165                                         vecbuf= NULL;
166                                 }
167                                 if(sock->type==SOCK_RGBA && colbuf) {
168                                         nsout[a]->data= pass_on_compbuf(colbuf);
169                                         colbuf= NULL;
170                                 }
171                         }
172                 }
173         }
174 }
175
176
177 CompBuf *get_cropped_compbuf(rcti *drect, float *rectf, int rectx, int recty, int type)
178 {
179         CompBuf *cbuf;
180         rcti disprect= *drect;
181         float *outfp;
182         int dx, y;
183         
184         if(disprect.xmax>rectx) disprect.xmax= rectx;
185         if(disprect.ymax>recty) disprect.ymax= recty;
186         if(disprect.xmin>= disprect.xmax) return NULL;
187         if(disprect.ymin>= disprect.ymax) return NULL;
188         
189         cbuf= alloc_compbuf(disprect.xmax-disprect.xmin, disprect.ymax-disprect.ymin, type, 1);
190         outfp= cbuf->rect;
191         rectf += type*(disprect.ymin*rectx + disprect.xmin);
192         dx= type*cbuf->x;
193         for(y=cbuf->y; y>0; y--, outfp+=dx, rectf+=type*rectx)
194                 memcpy(outfp, rectf, sizeof(float)*dx);
195         
196         return cbuf;
197 }
198
199 CompBuf *scalefast_compbuf(CompBuf *inbuf, int newx, int newy)
200 {
201         CompBuf *outbuf; 
202         float *rectf, *newrectf, *rf;
203         int x, y, c, pixsize= inbuf->type;
204         int ofsx, ofsy, stepx, stepy;
205         
206         if(inbuf->x==newx && inbuf->y==newy)
207                 return dupalloc_compbuf(inbuf);
208         
209         outbuf= alloc_compbuf(newx, newy, inbuf->type, 1);
210         newrectf= outbuf->rect;
211         
212         stepx = (65536.0 * (inbuf->x - 1.0) / (newx - 1.0)) + 0.5;
213         stepy = (65536.0 * (inbuf->y - 1.0) / (newy - 1.0)) + 0.5;
214         ofsy = 32768;
215         
216         for (y = newy; y > 0 ; y--){
217                 rectf = inbuf->rect;
218                 rectf += pixsize * (ofsy >> 16) * inbuf->x;
219
220                 ofsy += stepy;
221                 ofsx = 32768;
222                 
223                 for (x = newx ; x>0 ; x--) {
224                         
225                         rf= rectf + pixsize*(ofsx >> 16);
226                         for(c=0; c<pixsize; c++)
227                                 newrectf[c] = rf[c];
228                         
229                         newrectf+= pixsize;
230                         
231                         ofsx += stepx;
232                 }
233         }
234         
235         return outbuf;
236 }
237
238 void typecheck_compbuf_color(float *out, float *in, int outtype, int intype)
239 {
240         if(intype == outtype) {
241                 memcpy(out, in, sizeof(float)*outtype);
242         }
243         else if(outtype==CB_VAL) {
244                 if(intype==CB_VEC2) {
245                         *out= 0.5f*(in[0]+in[1]);
246                 }
247                 else if(intype==CB_VEC3) {
248                         *out= 0.333333f*(in[0]+in[1]+in[2]);
249                 }
250                 else if(intype==CB_RGBA) {
251                         *out= in[0]*0.35f + in[1]*0.45f + in[2]*0.2f;
252                 }
253         }
254         else if(outtype==CB_VEC2) {
255                 if(intype==CB_VAL) {
256                         out[0]= in[0];
257                         out[1]= in[0];
258                 }
259                 else if(intype==CB_VEC3) {
260                         out[0]= in[0];
261                         out[1]= in[1];
262                 }
263                 else if(intype==CB_RGBA) {
264                         out[0]= in[0];
265                         out[1]= in[1];
266                 }
267         }
268         else if(outtype==CB_VEC3) {
269                 if(intype==CB_VAL) {
270                         out[0]= in[0];
271                         out[1]= in[0];
272                         out[2]= in[0];
273                 }
274                 else if(intype==CB_VEC2) {
275                         out[0]= in[0];
276                         out[1]= in[1];
277                         out[2]= 0.0f;
278                 }
279                 else if(intype==CB_RGBA) {
280                         out[0]= in[0];
281                         out[1]= in[1];
282                         out[2]= in[2];
283                 }
284         }
285         else if(outtype==CB_RGBA) {
286                 if(intype==CB_VAL) {
287                         out[0]= in[0];
288                         out[1]= in[0];
289                         out[2]= in[0];
290                         out[3]= 1.0f;
291                 }
292                 else if(intype==CB_VEC2) {
293                         out[0]= in[0];
294                         out[1]= in[1];
295                         out[2]= 0.0f;
296                         out[3]= 1.0f;
297                 }
298                 else if(intype==CB_VEC3) {
299                         out[0]= in[0];
300                         out[1]= in[1];
301                         out[2]= in[2];
302                         out[3]= 1.0f;
303                 }
304         }
305 }
306
307 CompBuf *typecheck_compbuf(CompBuf *inbuf, int type)
308 {
309         if(inbuf && inbuf->type!=type) {
310                 CompBuf *outbuf;
311                 float *inrf, *outrf;
312                 int x;
313
314                 outbuf= alloc_compbuf(inbuf->x, inbuf->y, type, 1); 
315
316                 /* warning note: xof and yof are applied in pixelprocessor, but should be copied otherwise? */
317                 outbuf->xof= inbuf->xof;
318                 outbuf->yof= inbuf->yof;
319
320                 if(inbuf->rect_procedural) {
321                         outbuf->rect_procedural= inbuf->rect_procedural;
322                         VECCOPY(outbuf->procedural_size, inbuf->procedural_size);
323                         VECCOPY(outbuf->procedural_offset, inbuf->procedural_offset);
324                         outbuf->procedural_type= inbuf->procedural_type;
325                         outbuf->node= inbuf->node;
326                         return outbuf;
327                 }
328
329                 inrf= inbuf->rect;
330                 outrf= outbuf->rect;
331                 x= inbuf->x*inbuf->y;
332                 
333                 if(type==CB_VAL) {
334                         if(inbuf->type==CB_VEC2) {
335                                 for(; x>0; x--, outrf+= 1, inrf+= 2)
336                                         *outrf= 0.5f*(inrf[0]+inrf[1]);
337                         }
338                         else if(inbuf->type==CB_VEC3) {
339                                 for(; x>0; x--, outrf+= 1, inrf+= 3)
340                                         *outrf= 0.333333f*(inrf[0]+inrf[1]+inrf[2]);
341                         }
342                         else if(inbuf->type==CB_RGBA) {
343                                 for(; x>0; x--, outrf+= 1, inrf+= 4)
344                                         *outrf= inrf[0]*0.35f + inrf[1]*0.45f + inrf[2]*0.2f;
345                         }
346                 }
347                 else if(type==CB_VEC2) {
348                         if(inbuf->type==CB_VAL) {
349                                 for(; x>0; x--, outrf+= 2, inrf+= 1) {
350                                         outrf[0]= inrf[0];
351                                         outrf[1]= inrf[0];
352                                 }
353                         }
354                         else if(inbuf->type==CB_VEC3) {
355                                 for(; x>0; x--, outrf+= 2, inrf+= 3) {
356                                         outrf[0]= inrf[0];
357                                         outrf[1]= inrf[1];
358                                 }
359                         }
360                         else if(inbuf->type==CB_RGBA) {
361                                 for(; x>0; x--, outrf+= 2, inrf+= 4) {
362                                         outrf[0]= inrf[0];
363                                         outrf[1]= inrf[1];
364                                 }
365                         }
366                 }
367                 else if(type==CB_VEC3) {
368                         if(inbuf->type==CB_VAL) {
369                                 for(; x>0; x--, outrf+= 3, inrf+= 1) {
370                                         outrf[0]= inrf[0];
371                                         outrf[1]= inrf[0];
372                                         outrf[2]= inrf[0];
373                                 }
374                         }
375                         else if(inbuf->type==CB_VEC2) {
376                                 for(; x>0; x--, outrf+= 3, inrf+= 2) {
377                                         outrf[0]= inrf[0];
378                                         outrf[1]= inrf[1];
379                                         outrf[2]= 0.0f;
380                                 }
381                         }
382                         else if(inbuf->type==CB_RGBA) {
383                                 for(; x>0; x--, outrf+= 3, inrf+= 4) {
384                                         outrf[0]= inrf[0];
385                                         outrf[1]= inrf[1];
386                                         outrf[2]= inrf[2];
387                                 }
388                         }
389                 }
390                 else if(type==CB_RGBA) {
391                         if(inbuf->type==CB_VAL) {
392                                 for(; x>0; x--, outrf+= 4, inrf+= 1) {
393                                         outrf[0]= inrf[0];
394                                         outrf[1]= inrf[0];
395                                         outrf[2]= inrf[0];
396                                         outrf[3]= 1.0f;
397                                 }
398                         }
399                         else if(inbuf->type==CB_VEC2) {
400                                 for(; x>0; x--, outrf+= 4, inrf+= 2) {
401                                         outrf[0]= inrf[0];
402                                         outrf[1]= inrf[1];
403                                         outrf[2]= 0.0f;
404                                         outrf[3]= 1.0f;
405                                 }
406                         }
407                         else if(inbuf->type==CB_VEC3) {
408                                 for(; x>0; x--, outrf+= 4, inrf+= 3) {
409                                         outrf[0]= inrf[0];
410                                         outrf[1]= inrf[1];
411                                         outrf[2]= inrf[2];
412                                         outrf[3]= 1.0f;
413                                 }
414                         }
415                 }
416                 
417                 return outbuf;
418         }
419         return inbuf;
420 }
421
422 static float *compbuf_get_pixel(CompBuf *cbuf, float *defcol, float *use, int x, int y, int xrad, int yrad)
423 {
424         if(cbuf) {
425                 if(cbuf->rect_procedural) {
426                         cbuf->rect_procedural(cbuf, use, (float)x/(float)xrad, (float)y/(float)yrad);
427                         return use;
428                 }
429                 else {
430                         static float col[4]= {0.0f, 0.0f, 0.0f, 0.0f};
431                         
432                         /* map coords */
433                         x-= cbuf->xof;
434                         y-= cbuf->yof;
435                         
436                         if(y<-cbuf->yrad || y>= -cbuf->yrad+cbuf->y) return col;
437                         if(x<-cbuf->xrad || x>= -cbuf->xrad+cbuf->x) return col;
438                         
439                         return cbuf->rect + cbuf->type*( (cbuf->yrad+y)*cbuf->x + (cbuf->xrad+x) );
440                 }
441         }
442         else return defcol;
443 }
444
445 /* **************************************************** */
446
447 /* Pixel-to-Pixel operation, 1 Image in, 1 out */
448 void composit1_pixel_processor(bNode *node, CompBuf *out, CompBuf *src_buf, float *src_col,
449                                                                           void (*func)(bNode *, float *, float *), 
450                                                                           int src_type)
451 {
452         CompBuf *src_use;
453         float *outfp=out->rect, *srcfp;
454         float color[4]; /* local color if compbuf is procedural */
455         int xrad, yrad, x, y;
456         
457         src_use= typecheck_compbuf(src_buf, src_type);
458         
459         xrad= out->xrad;
460         yrad= out->yrad;
461         
462         for(y= -yrad; y<-yrad+out->y; y++) {
463                 for(x= -xrad; x<-xrad+out->x; x++, outfp+=out->type) {
464                         srcfp= compbuf_get_pixel(src_use, src_col, color, x, y, xrad, yrad);
465                         func(node, outfp, srcfp);
466                 }
467         }
468         
469         if(src_use!=src_buf)
470                 free_compbuf(src_use);
471 }
472
473 /* Pixel-to-Pixel operation, 2 Images in, 1 out */
474 void composit2_pixel_processor(bNode *node, CompBuf *out, CompBuf *src_buf, float *src_col,
475                                                                           CompBuf *fac_buf, float *fac, void (*func)(bNode *, float *, float *, float *), 
476                                                                           int src_type, int fac_type)
477 {
478         CompBuf *src_use, *fac_use;
479         float *outfp=out->rect, *srcfp, *facfp;
480         float color[4]; /* local color if compbuf is procedural */
481         int xrad, yrad, x, y;
482         
483         src_use= typecheck_compbuf(src_buf, src_type);
484         fac_use= typecheck_compbuf(fac_buf, fac_type);
485
486         xrad= out->xrad;
487         yrad= out->yrad;
488         
489         for(y= -yrad; y<-yrad+out->y; y++) {
490                 for(x= -xrad; x<-xrad+out->x; x++, outfp+=out->type) {
491                         srcfp= compbuf_get_pixel(src_use, src_col, color, x, y, xrad, yrad);
492                         facfp= compbuf_get_pixel(fac_use, fac, color, x, y, xrad, yrad);
493                         
494                         func(node, outfp, srcfp, facfp);
495                 }
496         }
497         if(src_use!=src_buf)
498                 free_compbuf(src_use);
499         if(fac_use!=fac_buf)
500                 free_compbuf(fac_use);
501 }
502
503 /* Pixel-to-Pixel operation, 3 Images in, 1 out */
504 void composit3_pixel_processor(bNode *node, CompBuf *out, CompBuf *src1_buf, float *src1_col, CompBuf *src2_buf, float *src2_col, 
505                                                                           CompBuf *fac_buf, float *fac, void (*func)(bNode *, float *, float *, float *, float *), 
506                                                                           int src1_type, int src2_type, int fac_type)
507 {
508         CompBuf *src1_use, *src2_use, *fac_use;
509         float *outfp=out->rect, *src1fp, *src2fp, *facfp;
510         float color[4]; /* local color if compbuf is procedural */
511         int xrad, yrad, x, y;
512         
513         src1_use= typecheck_compbuf(src1_buf, src1_type);
514         src2_use= typecheck_compbuf(src2_buf, src2_type);
515         fac_use= typecheck_compbuf(fac_buf, fac_type);
516         
517         xrad= out->xrad;
518         yrad= out->yrad;
519         
520         for(y= -yrad; y<-yrad+out->y; y++) {
521                 for(x= -xrad; x<-xrad+out->x; x++, outfp+=out->type) {
522                         src1fp= compbuf_get_pixel(src1_use, src1_col, color, x, y, xrad, yrad);
523                         src2fp= compbuf_get_pixel(src2_use, src2_col, color, x, y, xrad, yrad);
524                         facfp= compbuf_get_pixel(fac_use, fac, color, x, y, xrad, yrad);
525                         
526                         func(node, outfp, src1fp, src2fp, facfp);
527                 }
528         }
529         
530         if(src1_use!=src1_buf)
531                 free_compbuf(src1_use);
532         if(src2_use!=src2_buf)
533                 free_compbuf(src2_use);
534         if(fac_use!=fac_buf)
535                 free_compbuf(fac_use);
536 }
537
538 /* Pixel-to-Pixel operation, 4 Images in, 1 out */
539 void composit4_pixel_processor(bNode *node, CompBuf *out, CompBuf *src1_buf, float *src1_col, CompBuf *fac1_buf, float *fac1, 
540                                                                           CompBuf *src2_buf, float *src2_col, CompBuf *fac2_buf, float *fac2, 
541                                                                           void (*func)(bNode *, float *, float *, float *, float *, float *), 
542                                                                           int src1_type, int fac1_type, int src2_type, int fac2_type)
543 {
544         CompBuf *src1_use, *src2_use, *fac1_use, *fac2_use;
545         float *outfp=out->rect, *src1fp, *src2fp, *fac1fp, *fac2fp;
546         float color[4]; /* local color if compbuf is procedural */
547         int xrad, yrad, x, y;
548         
549         src1_use= typecheck_compbuf(src1_buf, src1_type);
550         src2_use= typecheck_compbuf(src2_buf, src2_type);
551         fac1_use= typecheck_compbuf(fac1_buf, fac1_type);
552         fac2_use= typecheck_compbuf(fac2_buf, fac2_type);
553         
554         xrad= out->xrad;
555         yrad= out->yrad;
556         
557         for(y= -yrad; y<-yrad+out->y; y++) {
558                 for(x= -xrad; x<-xrad+out->x; x++, outfp+=out->type) {
559                         src1fp= compbuf_get_pixel(src1_use, src1_col, color, x, y, xrad, yrad);
560                         src2fp= compbuf_get_pixel(src2_use, src2_col, color, x, y, xrad, yrad);
561                         fac1fp= compbuf_get_pixel(fac1_use, fac1, color, x, y, xrad, yrad);
562                         fac2fp= compbuf_get_pixel(fac2_use, fac2, color, x, y, xrad, yrad);
563                         
564                         func(node, outfp, src1fp, fac1fp, src2fp, fac2fp);
565                 }
566         }
567         
568         if(src1_use!=src1_buf)
569                 free_compbuf(src1_use);
570         if(src2_use!=src2_buf)
571                 free_compbuf(src2_use);
572         if(fac1_use!=fac1_buf)
573                 free_compbuf(fac1_use);
574         if(fac2_use!=fac2_buf)
575                 free_compbuf(fac2_use);
576 }
577
578
579 CompBuf *valbuf_from_rgbabuf(CompBuf *cbuf, int channel)
580 {
581         CompBuf *valbuf= alloc_compbuf(cbuf->x, cbuf->y, CB_VAL, 1);
582         float *valf, *rectf;
583         int tot;
584         
585         /* warning note: xof and yof are applied in pixelprocessor, but should be copied otherwise? */
586         valbuf->xof= cbuf->xof;
587         valbuf->yof= cbuf->yof;
588         
589         valf= valbuf->rect;
590
591         /* defaults to returning alpha channel */
592         if ((channel < CHAN_R) || (channel > CHAN_A)) channel = CHAN_A;
593
594         rectf= cbuf->rect + channel;
595         
596         for(tot= cbuf->x*cbuf->y; tot>0; tot--, valf++, rectf+=4)
597                 *valf= *rectf;
598         
599         return valbuf;
600 }
601
602 static CompBuf *generate_procedural_preview(CompBuf *cbuf, int newx, int newy)
603 {
604         CompBuf *outbuf;
605         float *outfp;
606         int xrad, yrad, x, y;
607         
608         outbuf= alloc_compbuf(newx, newy, CB_RGBA, 1);
609
610         outfp= outbuf->rect;
611         xrad= outbuf->xrad;
612         yrad= outbuf->yrad;
613         
614         for(y= -yrad; y<-yrad+outbuf->y; y++)
615                 for(x= -xrad; x<-xrad+outbuf->x; x++, outfp+=outbuf->type)
616                         cbuf->rect_procedural(cbuf, outfp, (float)x/(float)xrad, (float)y/(float)yrad);
617
618         return outbuf;
619 }
620
621 void generate_preview(void *data, bNode *node, CompBuf *stackbuf)
622 {
623         RenderData *rd= data;
624         bNodePreview *preview= node->preview;
625         int xsize, ysize;
626         int color_manage= rd->color_mgt_flag & R_COLOR_MANAGEMENT;
627         unsigned char *rect;
628         
629         if(preview && stackbuf) {
630                 CompBuf *cbuf, *stackbuf_use;
631                 
632                 if(stackbuf->rect==NULL && stackbuf->rect_procedural==NULL) return;
633                 
634                 stackbuf_use= typecheck_compbuf(stackbuf, CB_RGBA);
635
636                 if(stackbuf->x > stackbuf->y) {
637                         xsize= 140;
638                         ysize= (140*stackbuf->y)/stackbuf->x;
639                 }
640                 else {
641                         ysize= 140;
642                         xsize= (140*stackbuf->x)/stackbuf->y;
643                 }
644                 
645                 if(stackbuf_use->rect_procedural)
646                         cbuf= generate_procedural_preview(stackbuf_use, xsize, ysize);
647                 else
648                         cbuf= scalefast_compbuf(stackbuf_use, xsize, ysize);
649
650                 /* convert to byte for preview */
651                 rect= MEM_callocN(sizeof(unsigned char)*4*xsize*ysize, "bNodePreview.rect");
652
653                 if(color_manage)
654                         floatbuf_to_srgb_byte(cbuf->rect, rect, 0, xsize, 0, ysize, xsize);
655                 else
656                         floatbuf_to_byte(cbuf->rect, rect, 0, xsize, 0, ysize, xsize);
657                 
658                 free_compbuf(cbuf);
659                 if(stackbuf_use!=stackbuf)
660                         free_compbuf(stackbuf_use);
661
662                 BLI_lock_thread(LOCK_PREVIEW);
663
664                 if(preview->rect)
665                         MEM_freeN(preview->rect);
666                 preview->xsize= xsize;
667                 preview->ysize= ysize;
668                 preview->rect= rect;
669
670                 BLI_unlock_thread(LOCK_PREVIEW);
671         }
672 }
673
674 void do_rgba_to_yuva(bNode *UNUSED(node), float *out, float *in)
675 {
676         rgb_to_yuv(in[0],in[1],in[2], &out[0], &out[1], &out[2]);
677         out[3]=in[3];
678 }
679
680 void do_rgba_to_hsva(bNode *UNUSED(node), float *out, float *in)
681 {
682         rgb_to_hsv(in[0],in[1],in[2], &out[0], &out[1], &out[2]);
683         out[3]=in[3];
684 }
685
686 void do_rgba_to_ycca(bNode *UNUSED(node), float *out, float *in)
687 {
688         rgb_to_ycc(in[0],in[1],in[2], &out[0], &out[1], &out[2], BLI_YCC_ITU_BT601);
689         out[3]=in[3];
690 }
691
692 void do_yuva_to_rgba(bNode *UNUSED(node), float *out, float *in)
693 {
694         yuv_to_rgb(in[0],in[1],in[2], &out[0], &out[1], &out[2]);
695         out[3]=in[3];
696 }
697
698 void do_hsva_to_rgba(bNode *UNUSED(node), float *out, float *in)
699 {
700         hsv_to_rgb(in[0],in[1],in[2], &out[0], &out[1], &out[2]);
701         out[3]=in[3];
702 }
703
704 void do_ycca_to_rgba(bNode *UNUSED(node), float *out, float *in)
705 {
706         ycc_to_rgb(in[0],in[1],in[2], &out[0], &out[1], &out[2], BLI_YCC_ITU_BT601);
707         out[3]=in[3];
708 }
709
710 void do_copy_rgba(bNode *UNUSED(node), float *out, float *in)
711 {
712         QUATCOPY(out, in);
713 }
714
715 void do_copy_rgb(bNode *UNUSED(node), float *out, float *in)
716 {
717         VECCOPY(out, in);
718         out[3]= 1.0f;
719 }
720
721 void do_copy_value(bNode *UNUSED(node), float *out, float *in)
722 {
723         out[0]= in[0];
724 }
725
726 void do_copy_a_rgba(bNode *UNUSED(node), float *out, float *in, float *fac)
727 {
728         VECCOPY(out, in);
729         out[3]= *fac;
730 }
731
732 /* only accepts RGBA buffers */
733 void gamma_correct_compbuf(CompBuf *img, int inversed)
734 {
735         float *drect;
736         int x;
737
738         if(img->type!=CB_RGBA) return;
739
740         drect= img->rect;
741         if(inversed) {
742                 for(x=img->x*img->y; x>0; x--, drect+=4) {
743                         if(drect[0]>0.0f) drect[0]= sqrt(drect[0]); else drect[0]= 0.0f;
744                         if(drect[1]>0.0f) drect[1]= sqrt(drect[1]); else drect[1]= 0.0f;
745                         if(drect[2]>0.0f) drect[2]= sqrt(drect[2]); else drect[2]= 0.0f;
746                 }
747         }
748         else {
749                 for(x=img->x*img->y; x>0; x--, drect+=4) {
750                         if(drect[0]>0.0f) drect[0]*= drect[0]; else drect[0]= 0.0f;
751                         if(drect[1]>0.0f) drect[1]*= drect[1]; else drect[1]= 0.0f;
752                         if(drect[2]>0.0f) drect[2]*= drect[2]; else drect[2]= 0.0f;
753                 }
754         }
755 }
756
757 void premul_compbuf(CompBuf *img, int inversed)
758 {
759         float *drect;
760         int x;
761
762         if(img->type!=CB_RGBA) return;
763
764         drect= img->rect;
765         if(inversed) {
766                 for(x=img->x*img->y; x>0; x--, drect+=4) {
767                         if(fabs(drect[3]) < 1e-5f) {
768                                 drect[0]= 0.0f;
769                                 drect[1]= 0.0f;
770                                 drect[2]= 0.0f;
771                         }
772                         else {
773                                 drect[0] /= drect[3];
774                                 drect[1] /= drect[3];
775                                 drect[2] /= drect[3];
776                         }
777                 }
778         }
779         else {
780                 for(x=img->x*img->y; x>0; x--, drect+=4) {
781                         drect[0] *= drect[3];
782                         drect[1] *= drect[3];
783                         drect[2] *= drect[3];
784                 }
785         }
786 }
787
788
789
790 /*
791  *  2D Fast Hartley Transform, used for convolution
792  */
793
794 typedef float fREAL;
795
796 // returns next highest power of 2 of x, as well it's log2 in L2
797 static unsigned int nextPow2(unsigned int x, unsigned int* L2)
798 {
799         unsigned int pw, x_notpow2 = x & (x-1);
800         *L2 = 0;
801         while (x>>=1) ++(*L2);
802         pw = 1 << (*L2);
803         if (x_notpow2) { (*L2)++;  pw<<=1; }
804         return pw;
805 }
806
807 //------------------------------------------------------------------------------
808
809 // from FXT library by Joerg Arndt, faster in order bitreversal
810 // use: r = revbin_upd(r, h) where h = N>>1
811 static unsigned int revbin_upd(unsigned int r, unsigned int h)
812 {
813         while (!((r^=h)&h)) h >>= 1;
814         return r;
815 }
816 //------------------------------------------------------------------------------
817 static void FHT(fREAL* data, unsigned int M, unsigned int inverse)
818 {
819         double tt, fc, dc, fs, ds, a = M_PI;
820         fREAL t1, t2;
821         int n2, bd, bl, istep, k, len = 1 << M, n = 1;
822
823         int i, j = 0;
824         unsigned int Nh = len >> 1;
825         for (i=1;i<(len-1);++i) {
826                 j = revbin_upd(j, Nh);
827                 if (j>i) {
828                         t1 = data[i];
829                         data[i] = data[j];
830                         data[j] = t1;
831                 }
832         }
833
834         do {
835                 fREAL* data_n = &data[n];
836
837                 istep = n << 1;
838                 for (k=0; k<len; k+=istep) {
839                         t1 = data_n[k];
840                         data_n[k] = data[k] - t1;
841                         data[k] += t1;
842                 }
843
844                 n2 = n >> 1;
845                 if (n>2) {
846                         fc = dc = cos(a);
847                         fs = ds = sqrt(1.0 - fc*fc); //sin(a);
848                         bd = n-2;
849                         for (bl=1; bl<n2; bl++) {
850                                 fREAL* data_nbd = &data_n[bd];
851                                 fREAL* data_bd = &data[bd];
852                                 for (k=bl; k<len; k+=istep) {
853                                         t1 = fc*data_n[k] + fs*data_nbd[k];
854                                         t2 = fs*data_n[k] - fc*data_nbd[k];
855                                         data_n[k] = data[k] - t1;
856                                         data_nbd[k] = data_bd[k] - t2;
857                                         data[k] += t1;
858                                         data_bd[k] += t2;
859                                 }
860                                 tt = fc*dc - fs*ds;
861                                 fs = fs*dc + fc*ds;
862                                 fc = tt;
863                                 bd -= 2;
864                         }
865                 }
866
867                 if (n>1) {
868                         for (k=n2; k<len; k+=istep) {
869                                 t1 = data_n[k];
870                                 data_n[k] = data[k] - t1;
871                                 data[k] += t1;
872                         }
873                 }
874
875                 n = istep;
876                 a *= 0.5;
877         } while (n<len);
878
879         if (inverse) {
880                 fREAL sc = (fREAL)1 / (fREAL)len;
881                 for (k=0; k<len; ++k)
882                         data[k] *= sc;
883         }
884 }
885 //------------------------------------------------------------------------------
886 /* 2D Fast Hartley Transform, Mx/My -> log2 of width/height,
887         nzp -> the row where zero pad data starts,
888         inverse -> see above */
889 static void FHT2D(fREAL *data, unsigned int Mx, unsigned int My,
890                 unsigned int nzp, unsigned int inverse)
891 {
892         unsigned int i, j, Nx, Ny, maxy;
893         fREAL t;
894
895         Nx = 1 << Mx;
896         Ny = 1 << My;
897
898         // rows (forward transform skips 0 pad data)
899         maxy = inverse ? Ny : nzp;
900         for (j=0; j<maxy; ++j)
901                 FHT(&data[Nx*j], Mx, inverse);
902
903         // transpose data
904         if (Nx==Ny) {  // square
905                 for (j=0; j<Ny; ++j)
906                         for (i=j+1; i<Nx; ++i) {
907                                 unsigned int op = i + (j << Mx), np = j + (i << My);
908                                 t=data[op], data[op]=data[np], data[np]=t;
909                         }
910         }
911         else {  // rectangular
912                 unsigned int k, Nym = Ny-1, stm = 1 << (Mx + My);
913                 for (i=0; stm>0; i++) {
914                         #define pred(k) (((k & Nym) << Mx) + (k >> My))
915                         for (j=pred(i); j>i; j=pred(j));
916                         if (j < i) continue;
917                         for (k=i, j=pred(i); j!=i; k=j, j=pred(j), stm--)
918                                 { t=data[j], data[j]=data[k], data[k]=t; }
919                         #undef pred
920                         stm--;
921                 }
922         }
923         // swap Mx/My & Nx/Ny
924         i = Nx, Nx = Ny, Ny = i;
925         i = Mx, Mx = My, My = i;
926
927         // now columns == transposed rows
928         for (j=0; j<Ny; ++j)
929                 FHT(&data[Nx*j], Mx, inverse);
930
931         // finalize
932         for (j=0; j<=(Ny >> 1); j++) {
933                 unsigned int jm = (Ny - j) & (Ny-1);
934                 unsigned int ji = j << Mx;
935                 unsigned int jmi = jm << Mx;
936                 for (i=0; i<=(Nx >> 1); i++) {
937                         unsigned int im = (Nx - i) & (Nx-1);
938                         fREAL A = data[ji + i];
939                         fREAL B = data[jmi + i];
940                         fREAL C = data[ji + im];
941                         fREAL D = data[jmi + im];
942                         fREAL E = (fREAL)0.5*((A + D) - (B + C));
943                         data[ji + i] = A - E;
944                         data[jmi + i] = B + E;
945                         data[ji + im] = C + E;
946                         data[jmi + im] = D - E;
947                 }
948         }
949
950 }
951
952 //------------------------------------------------------------------------------
953
954 /* 2D convolution calc, d1 *= d2, M/N - > log2 of width/height */
955 static void fht_convolve(fREAL* d1, fREAL* d2, unsigned int M, unsigned int N)
956 {
957         fREAL a, b;
958         unsigned int i, j, k, L, mj, mL;
959         unsigned int m = 1 << M, n = 1 << N;
960         unsigned int m2 = 1 << (M-1), n2 = 1 << (N-1);
961         unsigned int mn2 = m << (N-1);
962
963         d1[0] *= d2[0];
964         d1[mn2] *= d2[mn2];
965         d1[m2] *= d2[m2];
966         d1[m2 + mn2] *= d2[m2 + mn2];
967         for (i=1; i<m2; i++) {
968                 k = m - i;
969                 a = d1[i]*d2[i] - d1[k]*d2[k];
970                 b = d1[k]*d2[i] + d1[i]*d2[k];
971                 d1[i] = (b + a)*(fREAL)0.5;
972                 d1[k] = (b - a)*(fREAL)0.5;
973                 a = d1[i + mn2]*d2[i + mn2] - d1[k + mn2]*d2[k + mn2];
974                 b = d1[k + mn2]*d2[i + mn2] + d1[i + mn2]*d2[k + mn2];
975                 d1[i + mn2] = (b + a)*(fREAL)0.5;
976                 d1[k + mn2] = (b - a)*(fREAL)0.5;
977         }
978         for (j=1; j<n2; j++) {
979                 L = n - j;
980                 mj = j << M;
981                 mL = L << M;
982                 a = d1[mj]*d2[mj] - d1[mL]*d2[mL];
983                 b = d1[mL]*d2[mj] + d1[mj]*d2[mL];
984                 d1[mj] = (b + a)*(fREAL)0.5;
985                 d1[mL] = (b - a)*(fREAL)0.5;
986                 a = d1[m2 + mj]*d2[m2 + mj] - d1[m2 + mL]*d2[m2 + mL];
987                 b = d1[m2 + mL]*d2[m2 + mj] + d1[m2 + mj]*d2[m2 + mL];
988                 d1[m2 + mj] = (b + a)*(fREAL)0.5;
989                 d1[m2 + mL] = (b - a)*(fREAL)0.5;
990         }
991         for (i=1; i<m2; i++) {
992                 k = m - i;
993                 for (j=1; j<n2; j++) {
994                         L = n - j;
995                         mj = j << M;
996                         mL = L << M;
997                         a = d1[i + mj]*d2[i + mj] - d1[k + mL]*d2[k + mL];
998                         b = d1[k + mL]*d2[i + mj] + d1[i + mj]*d2[k + mL];
999                         d1[i + mj] = (b + a)*(fREAL)0.5;
1000                         d1[k + mL] = (b - a)*(fREAL)0.5;
1001                         a = d1[i + mL]*d2[i + mL] - d1[k + mj]*d2[k + mj];
1002                         b = d1[k + mj]*d2[i + mL] + d1[i + mL]*d2[k + mj];
1003                         d1[i + mL] = (b + a)*(fREAL)0.5;
1004                         d1[k + mj] = (b - a)*(fREAL)0.5;
1005                 }
1006         }
1007 }
1008
1009 //------------------------------------------------------------------------------
1010
1011 void convolve(CompBuf* dst, CompBuf* in1, CompBuf* in2)
1012 {
1013         fREAL *data1, *data2, *fp;
1014         unsigned int w2, h2, hw, hh, log2_w, log2_h;
1015         fRGB wt, *colp;
1016         int x, y, ch;
1017         int xbl, ybl, nxb, nyb, xbsz, ybsz;
1018         int in2done = 0;
1019
1020         CompBuf* rdst = alloc_compbuf(in1->x, in1->y, in1->type, 1);
1021
1022         // convolution result width & height
1023         w2 = 2*in2->x - 1;
1024         h2 = 2*in2->y - 1;
1025         // FFT pow2 required size & log2
1026         w2 = nextPow2(w2, &log2_w);
1027         h2 = nextPow2(h2, &log2_h);
1028
1029         // alloc space
1030         data1 = (fREAL*)MEM_callocN(3*w2*h2*sizeof(fREAL), "convolve_fast FHT data1");
1031         data2 = (fREAL*)MEM_callocN(w2*h2*sizeof(fREAL), "convolve_fast FHT data2");
1032
1033         // normalize convolutor
1034         wt[0] = wt[1] = wt[2] = 0.f;
1035         for (y=0; y<in2->y; y++) {
1036                 colp = (fRGB*)&in2->rect[y*in2->x*in2->type];
1037                 for (x=0; x<in2->x; x++)
1038                         fRGB_add(wt, colp[x]);
1039         }
1040         if (wt[0] != 0.f) wt[0] = 1.f/wt[0];
1041         if (wt[1] != 0.f) wt[1] = 1.f/wt[1];
1042         if (wt[2] != 0.f) wt[2] = 1.f/wt[2];
1043         for (y=0; y<in2->y; y++) {
1044                 colp = (fRGB*)&in2->rect[y*in2->x*in2->type];
1045                 for (x=0; x<in2->x; x++)
1046                         fRGB_colormult(colp[x], wt);
1047         }
1048
1049         // copy image data, unpacking interleaved RGBA into separate channels
1050         // only need to calc data1 once
1051
1052         // block add-overlap
1053         hw = in2->x >> 1;
1054         hh = in2->y >> 1;
1055         xbsz = (w2 + 1) - in2->x;
1056         ybsz = (h2 + 1) - in2->y;
1057         nxb = in1->x / xbsz;
1058         if (in1->x % xbsz) nxb++;
1059         nyb = in1->y / ybsz;
1060         if (in1->y % ybsz) nyb++;
1061         for (ybl=0; ybl<nyb; ybl++) {
1062                 for (xbl=0; xbl<nxb; xbl++) {
1063
1064                         // each channel one by one
1065                         for (ch=0; ch<3; ch++) {
1066                                 fREAL* data1ch = &data1[ch*w2*h2];
1067
1068                                 // only need to calc fht data from in2 once, can re-use for every block
1069                                 if (!in2done) {
1070                                         // in2, channel ch -> data1
1071                                         for (y=0; y<in2->y; y++) {
1072                                                 fp = &data1ch[y*w2];
1073                                                 colp = (fRGB*)&in2->rect[y*in2->x*in2->type];
1074                                                 for (x=0; x<in2->x; x++)
1075                                                         fp[x] = colp[x][ch];
1076                                         }
1077                                 }
1078
1079                                 // in1, channel ch -> data2
1080                                 memset(data2, 0, w2*h2*sizeof(fREAL));
1081                                 for (y=0; y<ybsz; y++) {
1082                                         int yy = ybl*ybsz + y;
1083                                         if (yy >= in1->y) continue;
1084                                         fp = &data2[y*w2];
1085                                         colp = (fRGB*)&in1->rect[yy*in1->x*in1->type];
1086                                         for (x=0; x<xbsz; x++) {
1087                                                 int xx = xbl*xbsz + x;
1088                                                 if (xx >= in1->x) continue;
1089                                                 fp[x] = colp[xx][ch];
1090                                         }
1091                                 }
1092
1093                                 // forward FHT
1094                                 // zero pad data start is different for each == height+1
1095                                 if (!in2done) FHT2D(data1ch, log2_w, log2_h, in2->y+1, 0);
1096                                 FHT2D(data2, log2_w, log2_h, in2->y+1, 0);
1097
1098                                 // FHT2D transposed data, row/col now swapped
1099                                 // convolve & inverse FHT
1100                                 fht_convolve(data2, data1ch, log2_h, log2_w);
1101                                 FHT2D(data2, log2_h, log2_w, 0, 1);
1102                                 // data again transposed, so in order again
1103
1104                                 // overlap-add result
1105                                 for (y=0; y<(int)h2; y++) {
1106                                         const int yy = ybl*ybsz + y - hh;
1107                                         if ((yy < 0) || (yy >= in1->y)) continue;
1108                                         fp = &data2[y*w2];
1109                                         colp = (fRGB*)&rdst->rect[yy*in1->x*in1->type];
1110                                         for (x=0; x<(int)w2; x++) {
1111                                                 const int xx = xbl*xbsz + x - hw;
1112                                                 if ((xx < 0) || (xx >= in1->x)) continue;
1113                                                 colp[xx][ch] += fp[x];
1114                                         }
1115                                 }
1116
1117                         }
1118                         in2done = 1;
1119                 }
1120         }
1121
1122         MEM_freeN(data2);
1123         MEM_freeN(data1);
1124         memcpy(dst->rect, rdst->rect, sizeof(float)*dst->x*dst->y*dst->type);
1125         free_compbuf(rdst);
1126 }
1127
1128
1129 /*
1130  *
1131  * Utility functions qd_* should probably be intergrated better with other functions here.
1132  *
1133  */
1134 // sets fcol to pixelcolor at (x, y)
1135 void qd_getPixel(CompBuf* src, int x, int y, float* col)
1136 {
1137         if(src->rect_procedural) {
1138                 float bc[4];
1139                 src->rect_procedural(src, bc, (float)x/(float)src->xrad, (float)y/(float)src->yrad);
1140
1141                 switch(src->type){
1142                         /* these fallthrough to get all the channels */
1143                         case CB_RGBA: col[3]=bc[3]; 
1144                         case CB_VEC3: col[2]=bc[2];
1145                         case CB_VEC2: col[1]=bc[1];
1146                         case CB_VAL: col[0]=bc[0];
1147                 }
1148         }
1149         else if ((x >= 0) && (x < src->x) && (y >= 0) && (y < src->y)) {
1150                 float* bc = &src->rect[(x + y*src->x)*src->type];
1151                 switch(src->type){
1152                         /* these fallthrough to get all the channels */
1153                         case CB_RGBA: col[3]=bc[3]; 
1154                         case CB_VEC3: col[2]=bc[2];
1155                         case CB_VEC2: col[1]=bc[1];
1156                         case CB_VAL: col[0]=bc[0];
1157                 }
1158         }
1159         else {
1160                 switch(src->type){
1161                         /* these fallthrough to get all the channels */
1162                         case CB_RGBA: col[3]=0.0; 
1163                         case CB_VEC3: col[2]=0.0; 
1164                         case CB_VEC2: col[1]=0.0; 
1165                         case CB_VAL: col[0]=0.0; 
1166                 }
1167         }
1168 }
1169
1170 // sets pixel (x, y) to color col
1171 void qd_setPixel(CompBuf* src, int x, int y, float* col)
1172 {
1173         if ((x >= 0) && (x < src->x) && (y >= 0) && (y < src->y)) {
1174                 float* bc = &src->rect[(x + y*src->x)*src->type];
1175                 switch(src->type){
1176                         /* these fallthrough to get all the channels */
1177                         case CB_RGBA: bc[3]=col[3]; 
1178                         case CB_VEC3: bc[2]=col[2];
1179                         case CB_VEC2: bc[1]=col[1];
1180                         case CB_VAL: bc[0]=col[0];
1181                 }
1182         }
1183 }
1184
1185 // adds fcol to pixelcolor (x, y)
1186 void qd_addPixel(CompBuf* src, int x, int y, float* col)
1187 {
1188         if ((x >= 0) && (x < src->x) && (y >= 0) && (y < src->y)) {
1189                 float* bc = &src->rect[(x + y*src->x)*src->type];
1190                 bc[0] += col[0], bc[1] += col[1], bc[2] += col[2];
1191         }
1192 }
1193
1194 // multiplies pixel by factor value f
1195 void qd_multPixel(CompBuf* src, int x, int y, float f)
1196 {
1197         if ((x >= 0) && (x < src->x) && (y >= 0) && (y < src->y)) {
1198                 float* bc = &src->rect[(x + y*src->x)*src->type];
1199                 bc[0] *= f, bc[1] *= f, bc[2] *= f;
1200         }
1201 }
1202
1203 // bilinear interpolation with wraparound
1204 void qd_getPixelLerpWrap(CompBuf* src, float u, float v, float* col)
1205 {
1206         const float ufl = floor(u), vfl = floor(v);
1207         const int nx = (int)ufl % src->x, ny = (int)vfl % src->y;
1208         const int x1 = (nx < 0) ? (nx + src->x) : nx;
1209         const int y1 = (ny < 0) ? (ny + src->y) : ny;
1210         const int x2 = (x1 + 1) % src->x, y2 = (y1 + 1) % src->y;
1211         const float* c00 = &src->rect[(x1 + y1*src->x)*src->type];
1212         const float* c10 = &src->rect[(x2 + y1*src->x)*src->type];
1213         const float* c01 = &src->rect[(x1 + y2*src->x)*src->type];
1214         const float* c11 = &src->rect[(x2 + y2*src->x)*src->type];
1215         const float uf = u - ufl, vf = v - vfl;
1216         const float w00=(1.f-uf)*(1.f-vf), w10=uf*(1.f-vf), w01=(1.f-uf)*vf, w11=uf*vf;
1217         col[0] = w00*c00[0] + w10*c10[0] + w01*c01[0] + w11*c11[0];
1218         if (src->type != CB_VAL) {
1219                 col[1] = w00*c00[1] + w10*c10[1] + w01*c01[1] + w11*c11[1];
1220                 col[2] = w00*c00[2] + w10*c10[2] + w01*c01[2] + w11*c11[2];
1221                 col[3] = w00*c00[3] + w10*c10[3] + w01*c01[3] + w11*c11[3];
1222         }
1223 }
1224
1225 // as above, without wrap around
1226 void qd_getPixelLerp(CompBuf* src, float u, float v, float* col)
1227 {
1228         const float ufl = floor(u), vfl = floor(v);
1229         const int x1 = (int)ufl, y1 = (int)vfl;
1230         const int x2 = (int)ceil(u), y2 = (int)ceil(v);
1231         if ((x2 >= 0) && (y2 >= 0) && (x1 < src->x) && (y1 < src->y)) {
1232                 const float B[4] = {0,0,0,0};
1233                 const int ox1 = (x1 < 0), oy1 = (y1 < 0), ox2 = (x2 >= src->x), oy2 = (y2 >= src->y);
1234                 const float* c00 = (ox1 || oy1) ? B : &src->rect[(x1 + y1*src->x)*src->type];
1235                 const float* c10 = (ox2 || oy1) ? B : &src->rect[(x2 + y1*src->x)*src->type];
1236                 const float* c01 = (ox1 || oy2) ? B : &src->rect[(x1 + y2*src->x)*src->type];
1237                 const float* c11 = (ox2 || oy2) ? B : &src->rect[(x2 + y2*src->x)*src->type];
1238                 const float uf = u - ufl, vf = v - vfl;
1239                 const float w00=(1.f-uf)*(1.f-vf), w10=uf*(1.f-vf), w01=(1.f-uf)*vf, w11=uf*vf;
1240                 col[0] = w00*c00[0] + w10*c10[0] + w01*c01[0] + w11*c11[0];
1241                 if (src->type != CB_VAL) {
1242                         col[1] = w00*c00[1] + w10*c10[1] + w01*c01[1] + w11*c11[1];
1243                         col[2] = w00*c00[2] + w10*c10[2] + w01*c01[2] + w11*c11[2];
1244                         col[3] = w00*c00[3] + w10*c10[3] + w01*c01[3] + w11*c11[3];
1245                 }
1246         }
1247         else col[0] = col[1] = col[2] = col[3] = 0.f;
1248 }
1249
1250 // as above, sampling only one channel
1251 void qd_getPixelLerpChan(CompBuf* src, float u, float v, int chan, float* out)
1252 {
1253         const float ufl = floor(u), vfl = floor(v);
1254         const int x1 = (int)ufl, y1 = (int)vfl;
1255         const int x2 = (int)ceil(u), y2 = (int)ceil(v);
1256         if (chan >= src->type) chan = 0;
1257         if ((x2 >= 0) && (y2 >= 0) && (x1 < src->x) && (y1 < src->y)) {
1258                 const float B[4] = {0,0,0,0};
1259                 const int ox1 = (x1 < 0), oy1 = (y1 < 0), ox2 = (x2 >= src->x), oy2 = (y2 >= src->y);
1260                 const float* c00 = (ox1 || oy1) ? B : &src->rect[(x1 + y1*src->x)*src->type + chan];
1261                 const float* c10 = (ox2 || oy1) ? B : &src->rect[(x2 + y1*src->x)*src->type + chan];
1262                 const float* c01 = (ox1 || oy2) ? B : &src->rect[(x1 + y2*src->x)*src->type + chan];
1263                 const float* c11 = (ox2 || oy2) ? B : &src->rect[(x2 + y2*src->x)*src->type + chan];
1264                 const float uf = u - ufl, vf = v - vfl;
1265                 const float w00=(1.f-uf)*(1.f-vf), w10=uf*(1.f-vf), w01=(1.f-uf)*vf, w11=uf*vf;
1266                 out[0] = w00*c00[0] + w10*c10[0] + w01*c01[0] + w11*c11[0];
1267         }
1268         else *out = 0.f;
1269 }
1270
1271
1272 CompBuf* qd_downScaledCopy(CompBuf* src, int scale)
1273 {
1274         CompBuf* fbuf;
1275         if (scale <= 1)
1276                 fbuf = dupalloc_compbuf(src);
1277         else {
1278                 int nw = src->x/scale, nh = src->y/scale;
1279                 if ((2*(src->x % scale)) > scale) nw++;
1280                 if ((2*(src->y % scale)) > scale) nh++;
1281                 fbuf = alloc_compbuf(nw, nh, src->type, 1);
1282                 {
1283                         int x, y, xx, yy, sx, sy, mx, my;
1284                         float colsum[4] = {0.0f, 0.0f, 0.0f, 0.0f};
1285                         float fscale = 1.f/(float)(scale*scale);
1286                         for (y=0; y<nh; y++) {
1287                                 fRGB* fcolp = (fRGB*)&fbuf->rect[y*fbuf->x*fbuf->type];
1288                                 yy = y*scale;
1289                                 my = yy + scale;
1290                                 if (my > src->y) my = src->y;
1291                                 for (x=0; x<nw; x++) {
1292                                         xx = x*scale;
1293                                         mx = xx + scale;
1294                                         if (mx > src->x) mx = src->x;
1295                                         colsum[0] = colsum[1] = colsum[2] = 0.f;
1296                                         for (sy=yy; sy<my; sy++) {
1297                                                 fRGB* scolp = (fRGB*)&src->rect[sy*src->x*src->type];
1298                                                 for (sx=xx; sx<mx; sx++)
1299                                                         fRGB_add(colsum, scolp[sx]);
1300                                         }
1301                                         fRGB_mult(colsum, fscale);
1302                                         fRGB_copy(fcolp[x], colsum);
1303                                 }
1304                         }
1305                 }
1306         }
1307         return fbuf;
1308 }
1309
1310 // fast g.blur, per channel
1311 // xy var. bits 1 & 2 ca be used to blur in x or y direction separately
1312 void IIR_gauss(CompBuf* src, float sigma, int chan, int xy)
1313 {
1314         double q, q2, sc, cf[4], tsM[9], tsu[3], tsv[3];
1315         double *X, *Y, *W;
1316         int i, x, y, sz;
1317
1318         // <0.5 not valid, though can have a possibly useful sort of sharpening effect
1319         if (sigma < 0.5) return;
1320         
1321         if ((xy < 1) || (xy > 3)) xy = 3;
1322         
1323         // XXX The YVV macro defined below explicitly expects sources of at least 3x3 pixels,
1324         //     so just skiping blur along faulty direction if src's def is below that limit!
1325         if (src->x < 3) xy &= ~(int) 1;
1326         if (src->y < 3) xy &= ~(int) 2;
1327         if (xy < 1) return;
1328
1329         // see "Recursive Gabor Filtering" by Young/VanVliet
1330         // all factors here in double.prec. Required, because for single.prec it seems to blow up if sigma > ~200
1331         if (sigma >= 3.556)
1332                 q = 0.9804*(sigma - 3.556) + 2.5091;
1333         else // sigma >= 0.5
1334                 q = (0.0561*sigma + 0.5784)*sigma - 0.2568;
1335         q2 = q*q;
1336         sc = (1.1668 + q)*(3.203729649  + (2.21566 + q)*q);
1337         // no gabor filtering here, so no complex multiplies, just the regular coefs.
1338         // all negated here, so as not to have to recalc Triggs/Sdika matrix
1339         cf[1] = q*(5.788961737 + (6.76492 + 3.0*q)*q)/ sc;
1340         cf[2] = -q2*(3.38246 + 3.0*q)/sc;
1341         // 0 & 3 unchanged
1342         cf[3] = q2*q/sc;
1343         cf[0] = 1.0 - cf[1] - cf[2] - cf[3];
1344
1345         // Triggs/Sdika border corrections,
1346         // it seems to work, not entirely sure if it is actually totally correct,
1347         // Besides J.M.Geusebroek's anigauss.c (see http://www.science.uva.nl/~mark),
1348         // found one other implementation by Cristoph Lampert,
1349         // but neither seem to be quite the same, result seems to be ok so far anyway.
1350         // Extra scale factor here to not have to do it in filter,
1351         // though maybe this had something to with the precision errors
1352         sc = cf[0]/((1.0 + cf[1] - cf[2] + cf[3])*(1.0 - cf[1] - cf[2] - cf[3])*(1.0 + cf[2] + (cf[1] - cf[3])*cf[3]));
1353         tsM[0] = sc*(-cf[3]*cf[1] + 1.0 - cf[3]*cf[3] - cf[2]);
1354         tsM[1] = sc*((cf[3] + cf[1])*(cf[2] + cf[3]*cf[1]));
1355         tsM[2] = sc*(cf[3]*(cf[1] + cf[3]*cf[2]));
1356         tsM[3] = sc*(cf[1] + cf[3]*cf[2]);
1357         tsM[4] = sc*(-(cf[2] - 1.0)*(cf[2] + cf[3]*cf[1]));
1358         tsM[5] = sc*(-(cf[3]*cf[1] + cf[3]*cf[3] + cf[2] - 1.0)*cf[3]);
1359         tsM[6] = sc*(cf[3]*cf[1] + cf[2] + cf[1]*cf[1] - cf[2]*cf[2]);
1360         tsM[7] = sc*(cf[1]*cf[2] + cf[3]*cf[2]*cf[2] - cf[1]*cf[3]*cf[3] - cf[3]*cf[3]*cf[3] - cf[3]*cf[2] + cf[3]);
1361         tsM[8] = sc*(cf[3]*(cf[1] + cf[3]*cf[2]));
1362
1363 #define YVV(L)\
1364 {\
1365         W[0] = cf[0]*X[0] + cf[1]*X[0] + cf[2]*X[0] + cf[3]*X[0];\
1366         W[1] = cf[0]*X[1] + cf[1]*W[0] + cf[2]*X[0] + cf[3]*X[0];\
1367         W[2] = cf[0]*X[2] + cf[1]*W[1] + cf[2]*W[0] + cf[3]*X[0];\
1368         for (i=3; i<L; i++)\
1369                 W[i] = cf[0]*X[i] + cf[1]*W[i-1] + cf[2]*W[i-2] + cf[3]*W[i-3];\
1370         tsu[0] = W[L-1] - X[L-1];\
1371         tsu[1] = W[L-2] - X[L-1];\
1372         tsu[2] = W[L-3] - X[L-1];\
1373         tsv[0] = tsM[0]*tsu[0] + tsM[1]*tsu[1] + tsM[2]*tsu[2] + X[L-1];\
1374         tsv[1] = tsM[3]*tsu[0] + tsM[4]*tsu[1] + tsM[5]*tsu[2] + X[L-1];\
1375         tsv[2] = tsM[6]*tsu[0] + tsM[7]*tsu[1] + tsM[8]*tsu[2] + X[L-1];\
1376         Y[L-1] = cf[0]*W[L-1] + cf[1]*tsv[0] + cf[2]*tsv[1] + cf[3]*tsv[2];\
1377         Y[L-2] = cf[0]*W[L-2] + cf[1]*Y[L-1] + cf[2]*tsv[0] + cf[3]*tsv[1];\
1378         Y[L-3] = cf[0]*W[L-3] + cf[1]*Y[L-2] + cf[2]*Y[L-1] + cf[3]*tsv[0];\
1379         for (i=L-4; i>=0; i--)\
1380                 Y[i] = cf[0]*W[i] + cf[1]*Y[i+1] + cf[2]*Y[i+2] + cf[3]*Y[i+3];\
1381 }
1382
1383         // intermediate buffers
1384         sz = MAX2(src->x, src->y);
1385         X = MEM_callocN(sz*sizeof(double), "IIR_gauss X buf");
1386         Y = MEM_callocN(sz*sizeof(double), "IIR_gauss Y buf");
1387         W = MEM_callocN(sz*sizeof(double), "IIR_gauss W buf");
1388         if (xy & 1) {   // H
1389                 for (y=0; y<src->y; ++y) {
1390                         const int yx = y*src->x;
1391                         for (x=0; x<src->x; ++x)
1392                                 X[x] = src->rect[(x + yx)*src->type + chan];
1393                         YVV(src->x);
1394                         for (x=0; x<src->x; ++x)
1395                                 src->rect[(x + yx)*src->type + chan] = Y[x];
1396                 }
1397         }
1398         if (xy & 2) {   // V
1399                 for (x=0; x<src->x; ++x) {
1400                         for (y=0; y<src->y; ++y)
1401                                 X[y] = src->rect[(x + y*src->x)*src->type + chan];
1402                         YVV(src->y);
1403                         for (y=0; y<src->y; ++y)
1404                                 src->rect[(x + y*src->x)*src->type + chan] = Y[y];
1405                 }
1406         }
1407
1408         MEM_freeN(X);
1409         MEM_freeN(W);
1410         MEM_freeN(Y);
1411 #undef YVV
1412 }
1413