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