Fix for histogram Luma mode not working when waveform is in RGB mode.
[blender.git] / source / blender / blenkernel / intern / colortools.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) 2005 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/BL DUAL LICENSE BLOCK *****
28  */
29
30 #include <string.h>
31 #include <math.h>
32 #include <stdlib.h>
33 #include <float.h>
34
35 #ifdef WITH_LCMS
36 #include <lcms.h>
37 #endif
38
39 #include "MEM_guardedalloc.h"
40
41 #include "DNA_color_types.h"
42 #include "DNA_curve_types.h"
43
44 #include "BKE_colortools.h"
45 #include "BKE_curve.h"
46 #include "BKE_ipo.h"
47 #include "BKE_utildefines.h"
48
49 #include "BLI_blenlib.h"
50 #include "BLI_math.h"
51
52 #include "IMB_imbuf.h"
53 #include "IMB_imbuf_types.h"
54
55
56 void floatbuf_to_srgb_byte(float *rectf, unsigned char *rectc, int x1, int x2, int y1, int y2, int w)
57 {
58         int x, y;
59         float *rf= rectf;
60         float srgb[3];
61         unsigned char *rc= rectc;
62         
63         for(y=y1; y<y2; y++) {
64                 for(x=x1; x<x2; x++, rf+=4, rc+=4) {
65                         srgb[0]= linearrgb_to_srgb(rf[0]);
66                         srgb[1]= linearrgb_to_srgb(rf[1]);
67                         srgb[2]= linearrgb_to_srgb(rf[2]);
68
69                         rc[0]= FTOCHAR(srgb[0]);
70                         rc[1]= FTOCHAR(srgb[1]);
71                         rc[2]= FTOCHAR(srgb[2]);
72                         rc[3]= FTOCHAR(rf[3]);
73                 }
74         }
75 }
76
77 void floatbuf_to_byte(float *rectf, unsigned char *rectc, int x1, int x2, int y1, int y2, int w)
78 {
79         int x, y;
80         float *rf= rectf;
81         unsigned char *rc= rectc;
82         
83         for(y=y1; y<y2; y++) {
84                 for(x=x1; x<x2; x++, rf+=4, rc+=4) {
85                         rc[0]= FTOCHAR(rf[0]);
86                         rc[1]= FTOCHAR(rf[1]);
87                         rc[2]= FTOCHAR(rf[2]);
88                         rc[3]= FTOCHAR(rf[3]);
89                 }
90         }
91 }
92
93
94 /* ********************************* color curve ********************* */
95
96 /* ***************** operations on full struct ************* */
97
98 CurveMapping *curvemapping_add(int tot, float minx, float miny, float maxx, float maxy)
99 {
100         CurveMapping *cumap;
101         int a;
102         float clipminx, clipminy, clipmaxx, clipmaxy;
103         
104         cumap= MEM_callocN(sizeof(CurveMapping), "new curvemap");
105         cumap->flag= CUMA_DO_CLIP;
106         if(tot==4) cumap->cur= 3;               /* rhms, hack for 'col' curve? */
107         
108         clipminx = MIN2(minx, maxx);
109         clipminy = MIN2(miny, maxy);
110         clipmaxx = MAX2(minx, maxx);
111         clipmaxy = MAX2(miny, maxy);
112         
113         BLI_init_rctf(&cumap->curr, clipminx, clipmaxx, clipminy, clipmaxy);
114         cumap->clipr= cumap->curr;
115         
116         cumap->white[0]= cumap->white[1]= cumap->white[2]= 1.0f;
117         cumap->bwmul[0]= cumap->bwmul[1]= cumap->bwmul[2]= 1.0f;
118         
119         for(a=0; a<tot; a++) {
120                 cumap->cm[a].flag= CUMA_EXTEND_EXTRAPOLATE;
121                 cumap->cm[a].totpoint= 2;
122                 cumap->cm[a].curve= MEM_callocN(2*sizeof(CurveMapPoint), "curve points");
123                 
124                 cumap->cm[a].curve[0].x= minx;
125                 cumap->cm[a].curve[0].y= miny;
126                 cumap->cm[a].curve[1].x= maxx;
127                 cumap->cm[a].curve[1].y= maxy;
128         }       
129         return cumap;
130 }
131
132 void curvemapping_free(CurveMapping *cumap)
133 {
134         int a;
135         
136         if(cumap) {
137                 for(a=0; a<CM_TOT; a++) {
138                         if(cumap->cm[a].curve) MEM_freeN(cumap->cm[a].curve);
139                         if(cumap->cm[a].table) MEM_freeN(cumap->cm[a].table);
140                         if(cumap->cm[a].premultable) MEM_freeN(cumap->cm[a].premultable);
141                 }
142                 MEM_freeN(cumap);
143         }
144 }
145
146 CurveMapping *curvemapping_copy(CurveMapping *cumap)
147 {
148         int a;
149         
150         if(cumap) {
151                 CurveMapping *cumapn= MEM_dupallocN(cumap);
152                 for(a=0; a<CM_TOT; a++) {
153                         if(cumap->cm[a].curve) 
154                                 cumapn->cm[a].curve= MEM_dupallocN(cumap->cm[a].curve);
155                         if(cumap->cm[a].table) 
156                                 cumapn->cm[a].table= MEM_dupallocN(cumap->cm[a].table);
157                         if(cumap->cm[a].premultable) 
158                                 cumapn->cm[a].premultable= MEM_dupallocN(cumap->cm[a].premultable);
159                 }
160                 return cumapn;
161         }
162         return NULL;
163 }
164
165 void curvemapping_set_black_white(CurveMapping *cumap, float *black, float *white)
166 {
167         int a;
168         
169         if(white)
170                 VECCOPY(cumap->white, white);
171         if(black)
172                 VECCOPY(cumap->black, black);
173         
174         for(a=0; a<3; a++) {
175                 if(cumap->white[a]==cumap->black[a])
176                         cumap->bwmul[a]= 0.0f;
177                 else
178                         cumap->bwmul[a]= 1.0f/(cumap->white[a] - cumap->black[a]);
179         }       
180 }
181
182 /* ***************** operations on single curve ************* */
183 /* ********** NOTE: requires curvemapping_changed() call after ******** */
184
185 /* removes with flag set */
186 void curvemap_remove(CurveMap *cuma, int flag)
187 {
188         CurveMapPoint *cmp= MEM_mallocN((cuma->totpoint)*sizeof(CurveMapPoint), "curve points");
189         int a, b, removed=0;
190         
191         /* well, lets keep the two outer points! */
192         cmp[0]= cuma->curve[0];
193         for(a=1, b=1; a<cuma->totpoint-1; a++) {
194                 if(!(cuma->curve[a].flag & flag)) {
195                         cmp[b]= cuma->curve[a];
196                         b++;
197                 }
198                 else removed++;
199         }
200         cmp[b]= cuma->curve[a];
201         
202         MEM_freeN(cuma->curve);
203         cuma->curve= cmp;
204         cuma->totpoint -= removed;
205 }
206
207 void curvemap_insert(CurveMap *cuma, float x, float y)
208 {
209         CurveMapPoint *cmp= MEM_callocN((cuma->totpoint+1)*sizeof(CurveMapPoint), "curve points");
210         int a, b, foundloc= 0;
211                 
212         /* insert fragments of the old one and the new point to the new curve */
213         cuma->totpoint++;
214         for(a=0, b=0; a<cuma->totpoint; a++) {
215                 if((x < cuma->curve[a].x) && !foundloc) {
216                         cmp[a].x= x;
217                         cmp[a].y= y;
218                         cmp[a].flag= CUMA_SELECT;
219                         foundloc= 1;
220                 }
221                 else {
222                         cmp[a].x= cuma->curve[b].x;
223                         cmp[a].y= cuma->curve[b].y;
224                         cmp[a].flag= cuma->curve[b].flag;
225                         cmp[a].flag &= ~CUMA_SELECT; /* make sure old points don't remain selected */
226                         cmp[a].shorty= cuma->curve[b].shorty;
227                         b++;
228                 }
229         }
230
231         /* free old curve and replace it with new one */
232         MEM_freeN(cuma->curve);
233         cuma->curve= cmp;
234 }
235
236 void curvemap_reset(CurveMap *cuma, rctf *clipr, int preset)
237 {
238         if(cuma->curve)
239                 MEM_freeN(cuma->curve);
240
241         switch(preset) {
242                 case CURVE_PRESET_LINE: cuma->totpoint= 2; break;
243                 case CURVE_PRESET_SHARP: cuma->totpoint= 3; break;
244                 case CURVE_PRESET_SMOOTH: cuma->totpoint= 4; break;
245                 case CURVE_PRESET_MAX: cuma->totpoint= 2; break;
246                 case CURVE_PRESET_MID9: cuma->totpoint= 9;
247         }
248
249         cuma->curve= MEM_callocN(cuma->totpoint*sizeof(CurveMapPoint), "curve points");
250
251         switch(preset) {
252                 case CURVE_PRESET_LINE:
253                         cuma->curve[0].x= clipr->xmin;
254                         cuma->curve[0].y= clipr->ymin;
255                         cuma->curve[0].flag= 0;
256                         cuma->curve[1].x= clipr->xmax;
257                         cuma->curve[1].y= clipr->ymax;
258                         cuma->curve[1].flag= 0;
259                         break;
260                 case CURVE_PRESET_SHARP:
261                         cuma->curve[0].x= 0;
262                         cuma->curve[0].y= 1;
263                         cuma->curve[1].x= 0.33;
264                         cuma->curve[1].y= 0.33;
265                         cuma->curve[2].x= 1;
266                         cuma->curve[2].y= 0;
267                         break;
268                 case CURVE_PRESET_SMOOTH:
269                         cuma->curve[0].x= 0;
270                         cuma->curve[0].y= 1;
271                         cuma->curve[1].x= 0.25;
272                         cuma->curve[1].y= 0.92;
273                         cuma->curve[2].x= 0.75;
274                         cuma->curve[2].y= 0.08;
275                         cuma->curve[3].x= 1;
276                         cuma->curve[3].y= 0;
277                         break;
278                 case CURVE_PRESET_MAX:
279                         cuma->curve[0].x= 0;
280                         cuma->curve[0].y= 1;
281                         cuma->curve[1].x= 1;
282                         cuma->curve[1].y= 1;
283                         break;
284                 case CURVE_PRESET_MID9:
285                         {
286                                 int i;
287                                 for (i=0; i < cuma->totpoint; i++)
288                                 {
289                                         cuma->curve[i].x= i / ((float)cuma->totpoint-1);
290                                         cuma->curve[i].y= 0.5;
291                                 }
292                         }
293         }
294         
295         if(cuma->table) {
296                 MEM_freeN(cuma->table);
297                 cuma->table= NULL;
298         }
299 }
300
301 /* if type==1: vector, else auto */
302 void curvemap_sethandle(CurveMap *cuma, int type)
303 {
304         int a;
305         
306         for(a=0; a<cuma->totpoint; a++) {
307                 if(cuma->curve[a].flag & CUMA_SELECT) {
308                         if(type) cuma->curve[a].flag |= CUMA_VECTOR;
309                         else cuma->curve[a].flag &= ~CUMA_VECTOR;
310                 }
311         }
312 }
313
314 /* *********************** Making the tables and display ************** */
315
316 /* reduced copy of garbled calchandleNurb() code in curve.c */
317 static void calchandle_curvemap(BezTriple *bezt, BezTriple *prev, BezTriple *next, int mode)
318 {
319         float *p1,*p2,*p3,pt[3];
320         float dx1,dy1, dx,dy, vx,vy, len,len1,len2;
321         
322         if(bezt->h1==0 && bezt->h2==0) return;
323         
324         p2= bezt->vec[1];
325         
326         if(prev==NULL) {
327                 p3= next->vec[1];
328                 pt[0]= 2*p2[0]- p3[0];
329                 pt[1]= 2*p2[1]- p3[1];
330                 p1= pt;
331         }
332         else p1= prev->vec[1];
333         
334         if(next==NULL) {
335                 p1= prev->vec[1];
336                 pt[0]= 2*p2[0]- p1[0];
337                 pt[1]= 2*p2[1]- p1[1];
338                 p3= pt;
339         }
340         else p3= next->vec[1];
341         
342         dx= p2[0]- p1[0];
343         dy= p2[1]- p1[1];
344
345         len1= (float)sqrt(dx*dx+dy*dy);
346         
347         dx1= p3[0]- p2[0];
348         dy1= p3[1]- p2[1];
349
350         len2= (float)sqrt(dx1*dx1+dy1*dy1);
351         
352         if(len1==0.0f) len1=1.0f;
353         if(len2==0.0f) len2=1.0f;
354         
355         if(bezt->h1==HD_AUTO || bezt->h2==HD_AUTO) {    /* auto */
356                 vx= dx1/len2 + dx/len1;
357                 vy= dy1/len2 + dy/len1;
358                 
359                 len= 2.5614f*(float)sqrt(vx*vx + vy*vy);
360                 if(len!=0.0f) {
361                         
362                         if(bezt->h1==HD_AUTO) {
363                                 len1/=len;
364                                 *(p2-3)= *p2-vx*len1;
365                                 *(p2-2)= *(p2+1)-vy*len1;
366                         }
367                         if(bezt->h2==HD_AUTO) {
368                                 len2/=len;
369                                 *(p2+3)= *p2+vx*len2;
370                                 *(p2+4)= *(p2+1)+vy*len2;
371                         }
372                 }
373         }
374
375         if(bezt->h1==HD_VECT) { /* vector */
376                 dx/=3.0; 
377                 dy/=3.0; 
378                 *(p2-3)= *p2-dx;
379                 *(p2-2)= *(p2+1)-dy;
380         }
381         if(bezt->h2==HD_VECT) {
382                 dx1/=3.0; 
383                 dy1/=3.0; 
384                 *(p2+3)= *p2+dx1;
385                 *(p2+4)= *(p2+1)+dy1;
386         }
387 }
388
389 /* in X, out Y. 
390    X is presumed to be outside first or last */
391 static float curvemap_calc_extend(CurveMap *cuma, float x, float *first, float *last)
392 {
393         if(x <= first[0]) {
394                 if((cuma->flag & CUMA_EXTEND_EXTRAPOLATE)==0) {
395                         /* no extrapolate */
396                         return first[1];
397                 }
398                 else {
399                         if(cuma->ext_in[0]==0.0f)
400                                 return first[1] + cuma->ext_in[1]*10000.0f;
401                         else
402                                 return first[1] + cuma->ext_in[1]*(x - first[0])/cuma->ext_in[0];
403                 }
404         }
405         else if(x >= last[0]) {
406                 if((cuma->flag & CUMA_EXTEND_EXTRAPOLATE)==0) {
407                         /* no extrapolate */
408                         return last[1];
409                 }
410                 else {
411                         if(cuma->ext_out[0]==0.0f)
412                                 return last[1] - cuma->ext_out[1]*10000.0f;
413                         else
414                                 return last[1] + cuma->ext_out[1]*(x - last[0])/cuma->ext_out[0];
415                 }
416         }
417         return 0.0f;
418 }
419
420 /* only creates a table for a single channel in CurveMapping */
421 static void curvemap_make_table(CurveMap *cuma, rctf *clipr)
422 {
423         CurveMapPoint *cmp= cuma->curve;
424         BezTriple *bezt;
425         float *fp, *allpoints, *lastpoint, curf, range;
426         int a, totpoint;
427         
428         if(cuma->curve==NULL) return;
429         
430         /* default rect also is table range */
431         cuma->mintable= clipr->xmin;
432         cuma->maxtable= clipr->xmax;
433         
434         /* hrmf... we now rely on blender ipo beziers, these are more advanced */
435         bezt= MEM_callocN(cuma->totpoint*sizeof(BezTriple), "beztarr");
436         
437         for(a=0; a<cuma->totpoint; a++) {
438                 cuma->mintable= MIN2(cuma->mintable, cmp[a].x);
439                 cuma->maxtable= MAX2(cuma->maxtable, cmp[a].x);
440                 bezt[a].vec[1][0]= cmp[a].x;
441                 bezt[a].vec[1][1]= cmp[a].y;
442                 if(cmp[a].flag & CUMA_VECTOR)
443                         bezt[a].h1= bezt[a].h2= HD_VECT;
444                 else
445                         bezt[a].h1= bezt[a].h2= HD_AUTO;
446         }
447         
448         for(a=0; a<cuma->totpoint; a++) {
449                 if(a==0)
450                         calchandle_curvemap(bezt, NULL, bezt+1, 0);
451                 else if(a==cuma->totpoint-1)
452                         calchandle_curvemap(bezt+a, bezt+a-1, NULL, 0);
453                 else
454                         calchandle_curvemap(bezt+a, bezt+a-1, bezt+a+1, 0);
455         }
456         
457         /* first and last handle need correction, instead of pointing to center of next/prev, 
458                 we let it point to the closest handle */
459         if(cuma->totpoint>2) {
460                 float hlen, nlen, vec[3];
461                 
462                 if(bezt[0].h2==HD_AUTO) {
463                         
464                         hlen= len_v3v3(bezt[0].vec[1], bezt[0].vec[2]); /* original handle length */
465                         /* clip handle point */
466                         VECCOPY(vec, bezt[1].vec[0]);
467                         if(vec[0] < bezt[0].vec[1][0])
468                                 vec[0]= bezt[0].vec[1][0];
469                         
470                         sub_v3_v3(vec, bezt[0].vec[1]);
471                         nlen= len_v3(vec);
472                         if(nlen>FLT_EPSILON) {
473                                 mul_v3_fl(vec, hlen/nlen);
474                                 add_v3_v3v3(bezt[0].vec[2], vec, bezt[0].vec[1]);
475                                 sub_v3_v3v3(bezt[0].vec[0], bezt[0].vec[1], vec);
476                         }
477                 }
478                 a= cuma->totpoint-1;
479                 if(bezt[a].h2==HD_AUTO) {
480                         
481                         hlen= len_v3v3(bezt[a].vec[1], bezt[a].vec[0]); /* original handle length */
482                         /* clip handle point */
483                         VECCOPY(vec, bezt[a-1].vec[2]);
484                         if(vec[0] > bezt[a].vec[1][0])
485                                 vec[0]= bezt[a].vec[1][0];
486                         
487                         sub_v3_v3(vec, bezt[a].vec[1]);
488                         nlen= len_v3(vec);
489                         if(nlen>FLT_EPSILON) {
490                                 mul_v3_fl(vec, hlen/nlen);
491                                 add_v3_v3v3(bezt[a].vec[0], vec, bezt[a].vec[1]);
492                                 sub_v3_v3v3(bezt[a].vec[2], bezt[a].vec[1], vec);
493                         }
494                 }
495         }       
496         /* make the bezier curve */
497         if(cuma->table)
498                 MEM_freeN(cuma->table);
499         totpoint= (cuma->totpoint-1)*CM_RESOL;
500         fp= allpoints= MEM_callocN(totpoint*2*sizeof(float), "table");
501         
502         for(a=0; a<cuma->totpoint-1; a++, fp += 2*CM_RESOL) {
503                 correct_bezpart(bezt[a].vec[1], bezt[a].vec[2], bezt[a+1].vec[0], bezt[a+1].vec[1]);
504                 forward_diff_bezier(bezt[a].vec[1][0], bezt[a].vec[2][0], bezt[a+1].vec[0][0], bezt[a+1].vec[1][0], fp, CM_RESOL-1, 2*sizeof(float));   
505                 forward_diff_bezier(bezt[a].vec[1][1], bezt[a].vec[2][1], bezt[a+1].vec[0][1], bezt[a+1].vec[1][1], fp+1, CM_RESOL-1, 2*sizeof(float));
506         }
507         
508         /* store first and last handle for extrapolation, unit length */
509         cuma->ext_in[0]= bezt[0].vec[0][0] - bezt[0].vec[1][0];
510         cuma->ext_in[1]= bezt[0].vec[0][1] - bezt[0].vec[1][1];
511         range= sqrt(cuma->ext_in[0]*cuma->ext_in[0] + cuma->ext_in[1]*cuma->ext_in[1]);
512         cuma->ext_in[0]/= range;
513         cuma->ext_in[1]/= range;
514         
515         a= cuma->totpoint-1;
516         cuma->ext_out[0]= bezt[a].vec[1][0] - bezt[a].vec[2][0];
517         cuma->ext_out[1]= bezt[a].vec[1][1] - bezt[a].vec[2][1];
518         range= sqrt(cuma->ext_out[0]*cuma->ext_out[0] + cuma->ext_out[1]*cuma->ext_out[1]);
519         cuma->ext_out[0]/= range;
520         cuma->ext_out[1]/= range;
521         
522         /* cleanup */
523         MEM_freeN(bezt);
524
525         range= CM_TABLEDIV*(cuma->maxtable - cuma->mintable);
526         cuma->range= 1.0f/range;
527         
528         /* now make a table with CM_TABLE equal x distances */
529         fp= allpoints;
530         lastpoint= allpoints + 2*(totpoint-1);
531         cmp= MEM_callocN((CM_TABLE+1)*sizeof(CurveMapPoint), "dist table");
532         
533         for(a=0; a<=CM_TABLE; a++) {
534                 curf= cuma->mintable + range*(float)a;
535                 cmp[a].x= curf;
536                 
537                 /* get the first x coordinate larger than curf */
538                 while(curf >= fp[0] && fp!=lastpoint) {
539                         fp+=2;
540                 }
541                 if(fp==allpoints || (curf >= fp[0] && fp==lastpoint))
542                         cmp[a].y= curvemap_calc_extend(cuma, curf, allpoints, lastpoint);
543                 else {
544                         float fac1= fp[0] - fp[-2];
545                         float fac2= fp[0] - curf;
546                         if(fac1 > FLT_EPSILON)
547                                 fac1= fac2/fac1;
548                         else
549                                 fac1= 0.0f;
550                         cmp[a].y= fac1*fp[-1] + (1.0f-fac1)*fp[1];
551                 }
552         }
553         
554         MEM_freeN(allpoints);
555         cuma->table= cmp;
556 }
557
558 /* call when you do images etc, needs restore too. also verifies tables */
559 /* it uses a flag to prevent premul or free to happen twice */
560 void curvemapping_premultiply(CurveMapping *cumap, int restore)
561 {
562         int a;
563         
564         if(restore) {
565                 if(cumap->flag & CUMA_PREMULLED) {
566                         for(a=0; a<3; a++) {
567                                 MEM_freeN(cumap->cm[a].table);
568                                 cumap->cm[a].table= cumap->cm[a].premultable;
569                                 cumap->cm[a].premultable= NULL;
570                         }
571                         
572                         cumap->flag &= ~CUMA_PREMULLED;
573                 }
574         }
575         else {
576                 if((cumap->flag & CUMA_PREMULLED)==0) {
577                         /* verify and copy */
578                         for(a=0; a<3; a++) {
579                                 if(cumap->cm[a].table==NULL)
580                                         curvemap_make_table(cumap->cm+a, &cumap->clipr);
581                                 cumap->cm[a].premultable= cumap->cm[a].table;
582                                 cumap->cm[a].table= MEM_mallocN((CM_TABLE+1)*sizeof(CurveMapPoint), "premul table");
583                                 memcpy(cumap->cm[a].table, cumap->cm[a].premultable, (CM_TABLE+1)*sizeof(CurveMapPoint));
584                         }
585                         
586                         if(cumap->cm[3].table==NULL)
587                                 curvemap_make_table(cumap->cm+3, &cumap->clipr);
588                 
589                         /* premul */
590                         for(a=0; a<3; a++) {
591                                 int b;
592                                 for(b=0; b<=CM_TABLE; b++) {
593                                         cumap->cm[a].table[b].y= curvemap_evaluateF(cumap->cm+3, cumap->cm[a].table[b].y);
594                                 }
595                         }
596                         
597                         cumap->flag |= CUMA_PREMULLED;
598                 }
599         }
600 }
601
602 static int sort_curvepoints(const void *a1, const void *a2)
603 {
604         const struct CurveMapPoint *x1=a1, *x2=a2;
605         
606         if( x1->x > x2->x ) return 1;
607         else if( x1->x < x2->x) return -1;
608         return 0;
609 }
610
611 /* ************************ more CurveMapping calls *************** */
612
613 /* note; only does current curvemap! */
614 void curvemapping_changed(CurveMapping *cumap, int rem_doubles)
615 {
616         CurveMap *cuma= cumap->cm+cumap->cur;
617         CurveMapPoint *cmp= cuma->curve;
618         rctf *clipr= &cumap->clipr;
619         float thresh= 0.01f*(clipr->xmax - clipr->xmin);
620         float dx= 0.0f, dy= 0.0f;
621         int a;
622         
623         /* clamp with clip */
624         if(cumap->flag & CUMA_DO_CLIP) {
625                 for(a=0; a<cuma->totpoint; a++) {
626                         if(cmp[a].flag & CUMA_SELECT) {
627                                 if(cmp[a].x < clipr->xmin)
628                                         dx= MIN2(dx, cmp[a].x - clipr->xmin);
629                                 else if(cmp[a].x > clipr->xmax)
630                                         dx= MAX2(dx, cmp[a].x - clipr->xmax);
631                                 if(cmp[a].y < clipr->ymin)
632                                         dy= MIN2(dy, cmp[a].y - clipr->ymin);
633                                 else if(cmp[a].y > clipr->ymax)
634                                         dy= MAX2(dy, cmp[a].y - clipr->ymax);
635                         }
636                 }
637                 for(a=0; a<cuma->totpoint; a++) {
638                         if(cmp[a].flag & CUMA_SELECT) {
639                                 cmp[a].x -= dx;
640                                 cmp[a].y -= dy;
641                         }
642                 }
643         }
644         
645         
646         qsort(cmp, cuma->totpoint, sizeof(CurveMapPoint), sort_curvepoints);
647         
648         /* remove doubles, threshold set on 1% of default range */
649         if(rem_doubles && cuma->totpoint>2) {
650                 for(a=0; a<cuma->totpoint-1; a++) {
651                         dx= cmp[a].x - cmp[a+1].x;
652                         dy= cmp[a].y - cmp[a+1].y;
653                         if( sqrt(dx*dx + dy*dy) < thresh ) {
654                                 if(a==0) {
655                                         cmp[a+1].flag|= 2;
656                                         if(cmp[a+1].flag & CUMA_SELECT)
657                                                 cmp[a].flag |= CUMA_SELECT;
658                                 }
659                                 else {
660                                         cmp[a].flag|= 2;
661                                         if(cmp[a].flag & CUMA_SELECT)
662                                                 cmp[a+1].flag |= CUMA_SELECT;
663                                 }
664                                 break;  /* we assume 1 deletion per edit is ok */
665                         }
666                 }
667                 if(a != cuma->totpoint-1)
668                         curvemap_remove(cuma, 2);
669         }       
670         curvemap_make_table(cuma, clipr);
671 }
672
673 /* table should be verified */
674 float curvemap_evaluateF(CurveMap *cuma, float value)
675 {
676         float fi;
677         int i;
678
679         /* index in table */
680         fi= (value-cuma->mintable)*cuma->range;
681         i= (int)fi;
682         
683         /* fi is table float index and should check against table range i.e. [0.0 CM_TABLE] */
684         if(fi<0.0f || fi>CM_TABLE)
685                 return curvemap_calc_extend(cuma, value, &cuma->table[0].x, &cuma->table[CM_TABLE].x);
686         else {
687                 if(i<0) return cuma->table[0].y;
688                 if(i>=CM_TABLE) return cuma->table[CM_TABLE].y;
689                 
690                 fi= fi-(float)i;
691                 return (1.0f-fi)*cuma->table[i].y + (fi)*cuma->table[i+1].y; 
692         }
693 }
694
695 /* works with curve 'cur' */
696 float curvemapping_evaluateF(CurveMapping *cumap, int cur, float value)
697 {
698         CurveMap *cuma= cumap->cm+cur;
699         
700         /* allocate or bail out */
701         if(cuma->table==NULL) {
702                 curvemap_make_table(cuma, &cumap->clipr);
703                 if(cuma->table==NULL)
704                         return value;
705         }
706         return curvemap_evaluateF(cuma, value);
707 }
708
709 /* vector case */
710 void curvemapping_evaluate3F(CurveMapping *cumap, float *vecout, const float *vecin)
711 {
712         vecout[0]= curvemapping_evaluateF(cumap, 0, vecin[0]);
713         vecout[1]= curvemapping_evaluateF(cumap, 1, vecin[1]);
714         vecout[2]= curvemapping_evaluateF(cumap, 2, vecin[2]);
715 }
716
717 /* RGB case, no black/white points, no premult */
718 void curvemapping_evaluateRGBF(CurveMapping *cumap, float *vecout, const float *vecin)
719 {
720         vecout[0]= curvemapping_evaluateF(cumap, 0, curvemapping_evaluateF(cumap, 3, vecin[0]));
721         vecout[1]= curvemapping_evaluateF(cumap, 1, curvemapping_evaluateF(cumap, 3, vecin[1]));
722         vecout[2]= curvemapping_evaluateF(cumap, 2, curvemapping_evaluateF(cumap, 3, vecin[2]));
723 }
724
725
726 /* RGB with black/white points and premult. tables are checked */
727 void curvemapping_evaluate_premulRGBF(CurveMapping *cumap, float *vecout, const float *vecin)
728 {
729         float fac;
730         
731         fac= (vecin[0] - cumap->black[0])*cumap->bwmul[0];
732         vecout[0]= curvemap_evaluateF(cumap->cm, fac);
733         
734         fac= (vecin[1] - cumap->black[1])*cumap->bwmul[1];
735         vecout[1]= curvemap_evaluateF(cumap->cm+1, fac);
736         
737         fac= (vecin[2] - cumap->black[2])*cumap->bwmul[2];
738         vecout[2]= curvemap_evaluateF(cumap->cm+2, fac);
739 }
740
741 void colorcorrection_do_ibuf(ImBuf *ibuf, const char *profile)
742 {
743         if (ibuf->crect == NULL)
744         {
745 #ifdef WITH_LCMS
746                 cmsHPROFILE imageProfile, proofingProfile;
747                 cmsHTRANSFORM hTransform;
748                 
749                 ibuf->crect = MEM_mallocN(ibuf->x*ibuf->y*sizeof(int), "imbuf crect");
750
751                 imageProfile  = cmsCreate_sRGBProfile();
752                 proofingProfile = cmsOpenProfileFromFile(profile, "r");
753                 
754                 cmsErrorAction(LCMS_ERROR_SHOW);
755         
756                 hTransform = cmsCreateProofingTransform(imageProfile, TYPE_RGBA_8, imageProfile, TYPE_RGBA_8, 
757                                                                                           proofingProfile,
758                                                                                           INTENT_ABSOLUTE_COLORIMETRIC,
759                                                                                           INTENT_ABSOLUTE_COLORIMETRIC,
760                                                                                           cmsFLAGS_SOFTPROOFING);
761         
762                 cmsDoTransform(hTransform, ibuf->rect, ibuf->crect, ibuf->x * ibuf->y);
763         
764                 cmsDeleteTransform(hTransform);
765                 cmsCloseProfile(imageProfile);
766                 cmsCloseProfile(proofingProfile);
767 #else
768                 ibuf->crect = ibuf->rect;
769 #endif
770         }
771 }
772
773 /* only used for image editor curves */
774 void curvemapping_do_ibuf(CurveMapping *cumap, ImBuf *ibuf)
775 {
776         ImBuf *tmpbuf;
777         int pixel;
778         float *pix_in;
779         float col[3];
780         int stride= 4;
781         float *pix_out;
782         
783         if(ibuf==NULL)
784                 return;
785         if(ibuf->rect_float==NULL)
786                 IMB_float_from_rect(ibuf);
787         else if(ibuf->rect==NULL)
788                 imb_addrectImBuf(ibuf);
789         
790         if (!ibuf->rect || !ibuf->rect_float)
791                 return;
792         
793         /* work on a temp buffer, so can color manage afterwards.
794          * No worse off memory wise than comp nodes */
795         tmpbuf = IMB_dupImBuf(ibuf);
796         
797         curvemapping_premultiply(cumap, 0);
798         
799         pix_in= ibuf->rect_float;
800         pix_out= tmpbuf->rect_float;
801
802         if(ibuf->channels)
803                 stride= ibuf->channels;
804         
805         for(pixel= ibuf->x*ibuf->y; pixel>0; pixel--, pix_in+=stride, pix_out+=4) {
806                 if(stride<3) {
807                         col[0]= curvemap_evaluateF(cumap->cm, *pix_in);
808                         
809                         pix_out[1]= pix_out[2]= pix_out[3]= pix_out[0]= col[0];
810                 }
811                 else {
812                         curvemapping_evaluate_premulRGBF(cumap, col, pix_in);
813                         pix_out[0]= col[0];
814                         pix_out[1]= col[1];
815                         pix_out[2]= col[2];
816                         if(stride>3)
817                                 pix_out[3]= pix_in[3];
818                         else
819                                 pix_out[3]= 1.f;
820                 }
821         }
822         
823         IMB_rect_from_float(tmpbuf);
824         SWAP(unsigned int *, tmpbuf->rect, ibuf->rect);
825         IMB_freeImBuf(tmpbuf);
826         
827         curvemapping_premultiply(cumap, 1);
828 }
829
830 int curvemapping_RGBA_does_something(CurveMapping *cumap)
831 {
832         int a;
833         
834         if(cumap->black[0]!=0.0f) return 1;
835         if(cumap->black[1]!=0.0f) return 1;
836         if(cumap->black[2]!=0.0f) return 1;
837         if(cumap->white[0]!=1.0f) return 1;
838         if(cumap->white[1]!=1.0f) return 1;
839         if(cumap->white[2]!=1.0f) return 1;
840         
841         for(a=0; a<CM_TOT; a++) {
842                 if(cumap->cm[a].curve) {
843                         if(cumap->cm[a].totpoint!=2)  return 1;
844                         
845                         if(cumap->cm[a].curve[0].x != 0.0f) return 1;
846                         if(cumap->cm[a].curve[0].y != 0.0f) return 1;
847                         if(cumap->cm[a].curve[1].x != 1.0f) return 1;
848                         if(cumap->cm[a].curve[1].y != 1.0f) return 1;
849                 }
850         }
851         return 0;
852 }
853
854 void curvemapping_initialize(CurveMapping *cumap)
855 {
856         int a;
857         
858         if(cumap==NULL) return;
859         
860         for(a=0; a<CM_TOT; a++) {
861                 if(cumap->cm[a].table==NULL)
862                         curvemap_make_table(cumap->cm+a, &cumap->clipr);
863         }
864 }
865
866 void curvemapping_table_RGBA(CurveMapping *cumap, float **array, int *size)
867 {
868         int a;
869         
870         *size = CM_TABLE+1;
871         *array = MEM_callocN(sizeof(float)*(*size)*4, "CurveMapping");
872         curvemapping_initialize(cumap);
873
874         for(a=0; a<*size; a++) {
875                 if(cumap->cm[0].table)
876                         (*array)[a*4+0]= cumap->cm[0].table[a].y;
877                 if(cumap->cm[1].table)
878                         (*array)[a*4+1]= cumap->cm[1].table[a].y;
879                 if(cumap->cm[2].table)
880                         (*array)[a*4+2]= cumap->cm[2].table[a].y;
881                 if(cumap->cm[3].table)
882                         (*array)[a*4+3]= cumap->cm[3].table[a].y;
883         }
884 }
885
886 /* ***************** Histogram **************** */
887
888 #define INV_255         (1.f/255.f)
889
890 DO_INLINE int get_bin_float(float f)
891 {
892         int bin= (int)(f*255);
893
894         /* note: clamp integer instead of float to avoid problems with NaN */
895         CLAMP(bin, 0, 255);
896         
897         //return (int) (((f + 0.25) / 1.5) * 255);
898         
899         return bin;
900 }
901
902 DO_INLINE void save_sample_line(Scopes *scopes, const int idx, const float fx, float *rgb, float *ycc)
903 {
904         float yuv[3];
905
906         /* vectorscope*/
907         rgb_to_yuv(rgb[0], rgb[1], rgb[2], &yuv[0], &yuv[1], &yuv[2]);
908         scopes->vecscope[idx + 0] = yuv[1];
909         scopes->vecscope[idx + 1] = yuv[2];
910
911         /* waveform */
912         switch (scopes->wavefrm_mode) {
913                 case SCOPES_WAVEFRM_RGB:
914                         scopes->waveform_1[idx + 0] = fx;
915                         scopes->waveform_1[idx + 1] = rgb[0];
916                         scopes->waveform_2[idx + 0] = fx;
917                         scopes->waveform_2[idx + 1] = rgb[1];
918                         scopes->waveform_3[idx + 0] = fx;
919                         scopes->waveform_3[idx + 1] = rgb[2];
920                         break;
921                 case SCOPES_WAVEFRM_LUMA:
922                         scopes->waveform_1[idx + 0] = fx;
923                         scopes->waveform_1[idx + 1] = ycc[0];
924                         break;
925                 case SCOPES_WAVEFRM_YCC_JPEG:
926                 case SCOPES_WAVEFRM_YCC_709:
927                 case SCOPES_WAVEFRM_YCC_601:
928                         scopes->waveform_1[idx + 0] = fx;
929                         scopes->waveform_1[idx + 1] = ycc[0];
930                         scopes->waveform_2[idx + 0] = fx;
931                         scopes->waveform_2[idx + 1] = ycc[1];
932                         scopes->waveform_3[idx + 0] = fx;
933                         scopes->waveform_3[idx + 1] = ycc[2];
934                         break;
935         }
936 }
937
938 void scopes_update(Scopes *scopes, ImBuf *ibuf, int use_color_management)
939 {
940         int x, y, c, n, nl;
941         double div, divl;
942         float *rf=NULL;
943         unsigned char *rc=NULL;
944         unsigned int *bin_r, *bin_g, *bin_b, *bin_lum;
945         int savedlines, saveline;
946         float rgb[3], ycc[3], luma;
947         int ycc_mode=-1;
948
949         if (scopes->ok == 1 ) return;
950
951         if (scopes->hist.ymax == 0.f) scopes->hist.ymax = 1.f;
952
953         /* hmmmm */
954         if (!(ELEM(ibuf->channels, 3, 4))) return;
955         scopes->hist.channels = 3;
956         scopes->hist.x_resolution = 256;
957
958         switch (scopes->wavefrm_mode) {
959                 case SCOPES_WAVEFRM_RGB:
960                         ycc_mode = -1;
961                         break;
962                 case SCOPES_WAVEFRM_LUMA:
963                 case SCOPES_WAVEFRM_YCC_JPEG:
964                         ycc_mode = BLI_YCC_JFIF_0_255;
965                         break;
966                 case SCOPES_WAVEFRM_YCC_601:
967                         ycc_mode = BLI_YCC_ITU_BT601;
968                         break;
969                 case SCOPES_WAVEFRM_YCC_709:
970                         ycc_mode = BLI_YCC_ITU_BT709;
971                         break;
972         }
973
974         /* temp table to count pix value for histo */
975         bin_r = MEM_callocN(256 * sizeof(unsigned int), "temp historgram bins");
976         bin_g = MEM_callocN(256 * sizeof(unsigned int), "temp historgram bins");
977         bin_b = MEM_callocN(256 * sizeof(unsigned int), "temp historgram bins");
978         bin_lum = MEM_callocN(256 * sizeof(unsigned int), "temp historgram bins");
979
980         /* convert to number of lines with logarithmic scale */
981         scopes->sample_lines = (scopes->accuracy*0.01) * (scopes->accuracy*0.01) * ibuf->y;
982         
983         if (scopes->sample_full)
984                 scopes->sample_lines = ibuf->y;
985
986         /* scan the image */
987         savedlines=0;
988         for (c=0; c<3; c++) {
989                 scopes->minmax[c][0]=25500.0f;
990                 scopes->minmax[c][1]=-25500.0f;
991         }
992         
993         scopes->waveform_tot = ibuf->x*scopes->sample_lines;
994         
995         if (scopes->waveform_1)
996                 MEM_freeN(scopes->waveform_1);
997         if (scopes->waveform_2)
998                 MEM_freeN(scopes->waveform_2);
999         if (scopes->waveform_3)
1000                 MEM_freeN(scopes->waveform_3);
1001         if (scopes->vecscope)
1002                 MEM_freeN(scopes->vecscope);
1003         
1004         scopes->waveform_1= MEM_callocN(scopes->waveform_tot * 2 * sizeof(float), "waveform point channel 1");
1005         scopes->waveform_2= MEM_callocN(scopes->waveform_tot * 2 * sizeof(float), "waveform point channel 2");
1006         scopes->waveform_3= MEM_callocN(scopes->waveform_tot * 2 * sizeof(float), "waveform point channel 3");
1007         scopes->vecscope= MEM_callocN(scopes->waveform_tot * 2 * sizeof(float), "vectorscope point channel");
1008         
1009         if (ibuf->rect_float)
1010                 rf = ibuf->rect_float;
1011         else if (ibuf->rect)
1012                 rc = (unsigned char *)ibuf->rect;
1013
1014         for (y = 0; y < ibuf->y; y++) {
1015                 if (savedlines<scopes->sample_lines && y>=((savedlines)*ibuf->y)/(scopes->sample_lines+1)) {
1016                         saveline=1;
1017                 } else saveline=0;
1018                 for (x = 0; x < ibuf->x; x++) {
1019
1020                         if (ibuf->rect_float) {
1021                                 if (use_color_management)
1022                                         linearrgb_to_srgb_v3_v3(rgb, rf);
1023                                 else
1024                                         copy_v3_v3(rgb, rf);
1025                         }
1026                         else if (ibuf->rect) {
1027                                 for (c=0; c<3; c++)
1028                                         rgb[c] = rc[c] * INV_255;
1029                         }
1030
1031                         /* we still need luma for histogram */
1032                         luma = 0.299*rgb[0] + 0.587*rgb[1] + 0.114 * rgb[2];
1033
1034                         /* check for min max */
1035                         if(ycc_mode == -1 ) {
1036                                 for (c=0; c<3; c++) {
1037                                         if (rgb[c] < scopes->minmax[c][0]) scopes->minmax[c][0] = rgb[c];
1038                                         if (rgb[c] > scopes->minmax[c][1]) scopes->minmax[c][1] = rgb[c];
1039                                 }
1040                         }
1041                         else {
1042                                 rgb_to_ycc(rgb[0],rgb[1],rgb[2],&ycc[0],&ycc[1],&ycc[2], ycc_mode);
1043                                 for (c=0; c<3; c++) {
1044                                         ycc[c] *=INV_255;
1045                                         if (ycc[c] < scopes->minmax[c][0]) scopes->minmax[c][0] = ycc[c];
1046                                         if (ycc[c] > scopes->minmax[c][1]) scopes->minmax[c][1] = ycc[c];
1047                                 }
1048                         }
1049                         /* increment count for histo*/
1050                         bin_r[ get_bin_float(rgb[0]) ] += 1;
1051                         bin_g[ get_bin_float(rgb[1]) ] += 1;
1052                         bin_b[ get_bin_float(rgb[2]) ] += 1;
1053                         bin_lum[ get_bin_float(luma) ] += 1;
1054
1055                         /* save sample if needed */
1056                         if(saveline) {
1057                                 const float fx = (float)x / (float)ibuf->x;
1058                                 const int idx = 2*(ibuf->x*savedlines+x);
1059                                 save_sample_line(scopes, idx, fx, rgb, ycc);
1060                         }
1061
1062                         rf+= ibuf->channels;
1063                         rc+= ibuf->channels;
1064                 }
1065                 if (saveline)
1066                         savedlines +=1;
1067         }
1068
1069         /* convert hist data to float (proportional to max count) */
1070         n=0;
1071         nl=0;
1072         for (x=0; x<256; x++) {
1073                 if (bin_r[x] > n)
1074                         n = bin_r[x];
1075                 if (bin_g[x] > n)
1076                         n = bin_g[x];
1077                 if (bin_b[x] > n)
1078                         n = bin_b[x];
1079                 if (bin_lum[x] > nl)
1080                         nl = bin_lum[x];
1081         }
1082         div = 1.f/(double)n;
1083         divl = 1.f/(double)nl;
1084         for (x=0; x<256; x++) {
1085                 scopes->hist.data_r[x] = bin_r[x] * div;
1086                 scopes->hist.data_g[x] = bin_g[x] * div;
1087                 scopes->hist.data_b[x] = bin_b[x] * div;
1088                 scopes->hist.data_luma[x] = bin_lum[x] * divl;
1089         }
1090         MEM_freeN(bin_r);
1091         MEM_freeN(bin_g);
1092         MEM_freeN(bin_b);
1093         MEM_freeN(bin_lum);
1094
1095         scopes->ok = 1;
1096 }
1097
1098 void scopes_free(Scopes *scopes)
1099 {
1100         if (scopes->waveform_1) {
1101                 MEM_freeN(scopes->waveform_1);
1102                 scopes->waveform_1 = NULL;
1103         }
1104         if (scopes->waveform_2) {
1105                 MEM_freeN(scopes->waveform_2);
1106                 scopes->waveform_2 = NULL;
1107         }
1108         if (scopes->waveform_3) {
1109                 MEM_freeN(scopes->waveform_3);
1110                 scopes->waveform_3 = NULL;
1111         }
1112         if (scopes->vecscope) {
1113                 MEM_freeN(scopes->vecscope);
1114                 scopes->vecscope = NULL;
1115         }
1116 }
1117
1118 void scopes_new(Scopes *scopes)
1119 {
1120         scopes->accuracy=30.0;
1121         scopes->hist.mode=HISTO_MODE_RGB;
1122         scopes->wavefrm_alpha=0.3;
1123         scopes->vecscope_alpha=0.3;
1124         scopes->wavefrm_height= 100;
1125         scopes->vecscope_height= 100;
1126         scopes->hist.height= 100;
1127         scopes->ok= 0;
1128         scopes->waveform_1 = NULL;
1129         scopes->waveform_2 = NULL;
1130         scopes->waveform_3 = NULL;
1131         scopes->vecscope = NULL;
1132 }