svn merge https://svn.blender.org/svnroot/bf-blender/trunk/blender -r22130:22205
[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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, 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 #include "DNA_image_types.h"
44 #include "DNA_texture_types.h"
45
46 #include "BKE_colortools.h"
47 #include "BKE_curve.h"
48 #include "BKE_global.h"
49 #include "BKE_ipo.h"
50 #include "BKE_image.h"
51 #include "BKE_main.h"
52 #include "BKE_utildefines.h"
53
54 #include "BLI_blenlib.h"
55 #include "BLI_arithb.h"
56 #include "BLI_threads.h"
57
58 #include "IMB_imbuf.h"
59 #include "IMB_imbuf_types.h"
60
61 /* ********************************* color transforms ********************************* */
62
63 /*Transform linear RGB values to nonlinear RGB values. Rec.
64   709 is ITU-R Recommendation BT. 709 (1990) ``Basic
65   Parameter Values for the HDTV Standard for the Studio and
66   for International Programme Exchange'', formerly CCIR Rec.
67   709.*/
68 void gamma_correct_rec709(float *c, float gamma)
69 {
70         /* Rec. 709 gamma correction. */
71         float cc = 0.018f;
72         
73         if (*c < cc)
74             *c *= ((1.099f * (float)pow(cc, gamma)) - 0.099f) / cc;
75         else 
76             *c = (1.099f * (float)pow(*c, gamma)) - 0.099f;
77 }
78
79 void gamma_correct(float *c, float gamma)
80 {
81         *c = pow((*c), gamma);
82 }
83
84 float srgb_to_linearrgb(float c)
85 {
86         if (c < 0.04045f)
87                 return (c < 0.f)?0.f:c / 12.92;
88         else
89                 return pow((c + 0.055)/1.055, 2.4);
90 }
91
92 float linearrgb_to_srgb(float c)
93 {
94         if (c < 0.0031308)
95                 return (c < 0.f)?0.f:c * 12.92;
96         else
97                 return  1.055 * pow(c, 1.0/2.4) - 0.055;
98 }
99
100 /* utility function convert an RGB triplet from sRGB to linear RGB color space */
101 void color_manage_linearize(float *col_to, float *col_from)
102 {
103         col_to[0] = srgb_to_linearrgb(col_from[0]);
104         col_to[1] = srgb_to_linearrgb(col_from[1]);
105         col_to[2] = srgb_to_linearrgb(col_from[2]);
106 }
107
108 void floatbuf_to_srgb_byte(float *rectf, unsigned char *rectc, int x1, int x2, int y1, int y2, int w)
109 {
110         int x, y;
111         float *rf= rectf;
112         float srgb[3];
113         unsigned char *rc= rectc;
114         
115         for(y=y1; y<y2; y++) {
116                 for(x=x1; x<x2; x++, rf+=4, rc+=4) {
117                         srgb[0]= linearrgb_to_srgb(rf[0]);
118                         srgb[1]= linearrgb_to_srgb(rf[1]);
119                         srgb[2]= linearrgb_to_srgb(rf[2]);
120
121                         rc[0]= FTOCHAR(srgb[0]);
122                         rc[1]= FTOCHAR(srgb[1]);
123                         rc[2]= FTOCHAR(srgb[2]);
124                         rc[3]= FTOCHAR(rf[3]);
125                 }
126         }
127 }
128
129 void floatbuf_to_byte(float *rectf, unsigned char *rectc, int x1, int x2, int y1, int y2, int w)
130 {
131         int x, y;
132         float *rf= rectf;
133         unsigned char *rc= rectc;
134         
135         for(y=y1; y<y2; y++) {
136                 for(x=x1; x<x2; x++, rf+=4, rc+=4) {
137                         rc[0]= FTOCHAR(rf[0]);
138                         rc[1]= FTOCHAR(rf[1]);
139                         rc[2]= FTOCHAR(rf[2]);
140                         rc[3]= FTOCHAR(rf[3]);
141                 }
142         }
143 }
144
145
146 /* ********************************* color curve ********************* */
147
148 /* ***************** operations on full struct ************* */
149
150 CurveMapping *curvemapping_add(int tot, float minx, float miny, float maxx, float maxy)
151 {
152         CurveMapping *cumap;
153         int a;
154         float clipminx, clipminy, clipmaxx, clipmaxy;
155         
156         cumap= MEM_callocN(sizeof(CurveMapping), "new curvemap");
157         cumap->flag= CUMA_DO_CLIP;
158         if(tot==4) cumap->cur= 3;               /* rhms, hack for 'col' curve? */
159         
160         clipminx = MIN2(minx, maxx);
161         clipminy = MIN2(miny, maxy);
162         clipmaxx = MAX2(minx, maxx);
163         clipmaxy = MAX2(miny, maxy);
164         
165         BLI_init_rctf(&cumap->curr, clipminx, clipmaxx, clipminy, clipmaxy);
166         cumap->clipr= cumap->curr;
167         
168         cumap->white[0]= cumap->white[1]= cumap->white[2]= 1.0f;
169         cumap->bwmul[0]= cumap->bwmul[1]= cumap->bwmul[2]= 1.0f;
170         
171         for(a=0; a<tot; a++) {
172                 cumap->cm[a].flag= CUMA_EXTEND_EXTRAPOLATE;
173                 cumap->cm[a].totpoint= 2;
174                 cumap->cm[a].curve= MEM_callocN(2*sizeof(CurveMapPoint), "curve points");
175                 
176                 cumap->cm[a].curve[0].x= minx;
177                 cumap->cm[a].curve[0].y= miny;
178                 cumap->cm[a].curve[1].x= maxx;
179                 cumap->cm[a].curve[1].y= maxy;
180         }       
181         return cumap;
182 }
183
184 void curvemapping_free(CurveMapping *cumap)
185 {
186         int a;
187         
188         if(cumap) {
189                 for(a=0; a<CM_TOT; a++) {
190                         if(cumap->cm[a].curve) MEM_freeN(cumap->cm[a].curve);
191                         if(cumap->cm[a].table) MEM_freeN(cumap->cm[a].table);
192                 }
193                 MEM_freeN(cumap);
194         }
195 }
196
197 CurveMapping *curvemapping_copy(CurveMapping *cumap)
198 {
199         int a;
200         
201         if(cumap) {
202                 CurveMapping *cumapn= MEM_dupallocN(cumap);
203                 for(a=0; a<CM_TOT; a++) {
204                         if(cumap->cm[a].curve) 
205                                 cumapn->cm[a].curve= MEM_dupallocN(cumap->cm[a].curve);
206                         if(cumap->cm[a].table) 
207                                 cumapn->cm[a].table= MEM_dupallocN(cumap->cm[a].table);
208                 }
209                 return cumapn;
210         }
211         return NULL;
212 }
213
214 void curvemapping_set_black_white(CurveMapping *cumap, float *black, float *white)
215 {
216         int a;
217         
218         if(white)
219                 VECCOPY(cumap->white, white);
220         if(black)
221                 VECCOPY(cumap->black, black);
222         
223         for(a=0; a<3; a++) {
224                 if(cumap->white[a]==cumap->black[a])
225                         cumap->bwmul[a]= 0.0f;
226                 else
227                         cumap->bwmul[a]= 1.0f/(cumap->white[a] - cumap->black[a]);
228         }       
229 }
230
231 /* ***************** operations on single curve ************* */
232 /* ********** NOTE: requires curvemapping_changed() call after ******** */
233
234 /* removes with flag set */
235 void curvemap_remove(CurveMap *cuma, int flag)
236 {
237         CurveMapPoint *cmp= MEM_mallocN((cuma->totpoint)*sizeof(CurveMapPoint), "curve points");
238         int a, b, removed=0;
239         
240         /* well, lets keep the two outer points! */
241         cmp[0]= cuma->curve[0];
242         for(a=1, b=1; a<cuma->totpoint-1; a++) {
243                 if(!(cuma->curve[a].flag & flag)) {
244                         cmp[b]= cuma->curve[a];
245                         b++;
246                 }
247                 else removed++;
248         }
249         cmp[b]= cuma->curve[a];
250         
251         MEM_freeN(cuma->curve);
252         cuma->curve= cmp;
253         cuma->totpoint -= removed;
254 }
255
256 void curvemap_insert(CurveMap *cuma, float x, float y)
257 {
258         CurveMapPoint *cmp= MEM_callocN((cuma->totpoint+1)*sizeof(CurveMapPoint), "curve points");
259         int a, b, foundloc= 0;
260                 
261         /* insert fragments of the old one and the new point to the new curve */
262         cuma->totpoint++;
263         for(a=0, b=0; a<cuma->totpoint; a++) {
264                 if((x < cuma->curve[a].x) && !foundloc) {
265                         cmp[a].x= x;
266                         cmp[a].y= y;
267                         cmp[a].flag= CUMA_SELECT;
268                         foundloc= 1;
269                 }
270                 else {
271                         cmp[a].x= cuma->curve[b].x;
272                         cmp[a].y= cuma->curve[b].y;
273                         cmp[a].flag= cuma->curve[b].flag;
274                         cmp[a].flag &= ~CUMA_SELECT; /* make sure old points don't remain selected */
275                         cmp[a].shorty= cuma->curve[b].shorty;
276                         b++;
277                 }
278         }
279
280         /* free old curve and replace it with new one */
281         MEM_freeN(cuma->curve);
282         cuma->curve= cmp;
283 }
284
285 void curvemap_reset(CurveMap *cuma, rctf *clipr)
286 {
287         cuma->totpoint= 2;
288         
289         cuma->curve[0].x= clipr->xmin;
290         cuma->curve[0].y= clipr->ymin;
291         cuma->curve[0].flag= 0;
292         cuma->curve[1].x= clipr->xmax;
293         cuma->curve[1].y= clipr->ymax;
294         cuma->curve[1].flag= 0;
295         
296         if(cuma->table) {
297                 MEM_freeN(cuma->table);
298                 cuma->table= NULL;
299         }
300 }
301
302 /* if type==1: vector, else auto */
303 void curvemap_sethandle(CurveMap *cuma, int type)
304 {
305         int a;
306         
307         for(a=0; a<cuma->totpoint; a++) {
308                 if(cuma->curve[a].flag & CUMA_SELECT) {
309                         if(type) cuma->curve[a].flag |= CUMA_VECTOR;
310                         else cuma->curve[a].flag &= ~CUMA_VECTOR;
311                 }
312         }
313 }
314
315 /* *********************** Making the tables and display ************** */
316
317 /* reduced copy of garbled calchandleNurb() code in curve.c */
318 static void calchandle_curvemap(BezTriple *bezt, BezTriple *prev, BezTriple *next, int mode)
319 {
320         float *p1,*p2,*p3,pt[3];
321         float dx1,dy1, dx,dy, vx,vy, len,len1,len2;
322         
323         if(bezt->h1==0 && bezt->h2==0) return;
324         
325         p2= bezt->vec[1];
326         
327         if(prev==NULL) {
328                 p3= next->vec[1];
329                 pt[0]= 2*p2[0]- p3[0];
330                 pt[1]= 2*p2[1]- p3[1];
331                 p1= pt;
332         }
333         else p1= prev->vec[1];
334         
335         if(next==NULL) {
336                 p1= prev->vec[1];
337                 pt[0]= 2*p2[0]- p1[0];
338                 pt[1]= 2*p2[1]- p1[1];
339                 p3= pt;
340         }
341         else p3= next->vec[1];
342         
343         dx= p2[0]- p1[0];
344         dy= p2[1]- p1[1];
345
346         len1= (float)sqrt(dx*dx+dy*dy);
347         
348         dx1= p3[0]- p2[0];
349         dy1= p3[1]- p2[1];
350
351         len2= (float)sqrt(dx1*dx1+dy1*dy1);
352         
353         if(len1==0.0f) len1=1.0f;
354         if(len2==0.0f) len2=1.0f;
355         
356         if(bezt->h1==HD_AUTO || bezt->h2==HD_AUTO) {    /* auto */
357                 vx= dx1/len2 + dx/len1;
358                 vy= dy1/len2 + dy/len1;
359                 
360                 len= 2.5614f*(float)sqrt(vx*vx + vy*vy);
361                 if(len!=0.0f) {
362                         
363                         if(bezt->h1==HD_AUTO) {
364                                 len1/=len;
365                                 *(p2-3)= *p2-vx*len1;
366                                 *(p2-2)= *(p2+1)-vy*len1;
367                         }
368                         if(bezt->h2==HD_AUTO) {
369                                 len2/=len;
370                                 *(p2+3)= *p2+vx*len2;
371                                 *(p2+4)= *(p2+1)+vy*len2;
372                         }
373                 }
374         }
375
376         if(bezt->h1==HD_VECT) { /* vector */
377                 dx/=3.0; 
378                 dy/=3.0; 
379                 *(p2-3)= *p2-dx;
380                 *(p2-2)= *(p2+1)-dy;
381         }
382         if(bezt->h2==HD_VECT) {
383                 dx1/=3.0; 
384                 dy1/=3.0; 
385                 *(p2+3)= *p2+dx1;
386                 *(p2+4)= *(p2+1)+dy1;
387         }
388 }
389
390 /* in X, out Y. 
391    X is presumed to be outside first or last */
392 static float curvemap_calc_extend(CurveMap *cuma, float x, float *first, float *last)
393 {
394         if(x <= first[0]) {
395                 if((cuma->flag & CUMA_EXTEND_EXTRAPOLATE)==0) {
396                         /* no extrapolate */
397                         return first[1];
398                 }
399                 else {
400                         if(cuma->ext_in[0]==0.0f)
401                                 return first[1] + cuma->ext_in[1]*10000.0f;
402                         else
403                                 return first[1] + cuma->ext_in[1]*(x - first[0])/cuma->ext_in[0];
404                 }
405         }
406         else if(x >= last[0]) {
407                 if((cuma->flag & CUMA_EXTEND_EXTRAPOLATE)==0) {
408                         /* no extrapolate */
409                         return last[1];
410                 }
411                 else {
412                         if(cuma->ext_out[0]==0.0f)
413                                 return last[1] - cuma->ext_out[1]*10000.0f;
414                         else
415                                 return last[1] + cuma->ext_out[1]*(x - last[0])/cuma->ext_out[0];
416                 }
417         }
418         return 0.0f;
419 }
420
421 /* only creates a table for a single channel in CurveMapping */
422 static void curvemap_make_table(CurveMap *cuma, rctf *clipr)
423 {
424         CurveMapPoint *cmp= cuma->curve;
425         BezTriple *bezt;
426         float *fp, *allpoints, *lastpoint, curf, range;
427         int a, totpoint;
428         
429         if(cuma->curve==NULL) return;
430         
431         /* default rect also is table range */
432         cuma->mintable= clipr->xmin;
433         cuma->maxtable= clipr->xmax;
434         
435         /* hrmf... we now rely on blender ipo beziers, these are more advanced */
436         bezt= MEM_callocN(cuma->totpoint*sizeof(BezTriple), "beztarr");
437         
438         for(a=0; a<cuma->totpoint; a++) {
439                 cuma->mintable= MIN2(cuma->mintable, cmp[a].x);
440                 cuma->maxtable= MAX2(cuma->maxtable, cmp[a].x);
441                 bezt[a].vec[1][0]= cmp[a].x;
442                 bezt[a].vec[1][1]= cmp[a].y;
443                 if(cmp[a].flag & CUMA_VECTOR)
444                         bezt[a].h1= bezt[a].h2= HD_VECT;
445                 else
446                         bezt[a].h1= bezt[a].h2= HD_AUTO;
447         }
448         
449         for(a=0; a<cuma->totpoint; a++) {
450                 if(a==0)
451                         calchandle_curvemap(bezt, NULL, bezt+1, 0);
452                 else if(a==cuma->totpoint-1)
453                         calchandle_curvemap(bezt+a, bezt+a-1, NULL, 0);
454                 else
455                         calchandle_curvemap(bezt+a, bezt+a-1, bezt+a+1, 0);
456         }
457         
458         /* first and last handle need correction, instead of pointing to center of next/prev, 
459                 we let it point to the closest handle */
460         if(cuma->totpoint>2) {
461                 float hlen, nlen, vec[3];
462                 
463                 if(bezt[0].h2==HD_AUTO) {
464                         
465                         hlen= VecLenf(bezt[0].vec[1], bezt[0].vec[2]);  /* original handle length */
466                         /* clip handle point */
467                         VECCOPY(vec, bezt[1].vec[0]);
468                         if(vec[0] < bezt[0].vec[1][0])
469                                 vec[0]= bezt[0].vec[1][0];
470                         
471                         VecSubf(vec, vec, bezt[0].vec[1]);
472                         nlen= VecLength(vec);
473                         if(nlen>FLT_EPSILON) {
474                                 VecMulf(vec, hlen/nlen);
475                                 VecAddf(bezt[0].vec[2], vec, bezt[0].vec[1]);
476                                 VecSubf(bezt[0].vec[0], bezt[0].vec[1], vec);
477                         }
478                 }
479                 a= cuma->totpoint-1;
480                 if(bezt[a].h2==HD_AUTO) {
481                         
482                         hlen= VecLenf(bezt[a].vec[1], bezt[a].vec[0]);  /* original handle length */
483                         /* clip handle point */
484                         VECCOPY(vec, bezt[a-1].vec[2]);
485                         if(vec[0] > bezt[a].vec[1][0])
486                                 vec[0]= bezt[a].vec[1][0];
487                         
488                         VecSubf(vec, vec, bezt[a].vec[1]);
489                         nlen= VecLength(vec);
490                         if(nlen>FLT_EPSILON) {
491                                 VecMulf(vec, hlen/nlen);
492                                 VecAddf(bezt[a].vec[0], vec, bezt[a].vec[1]);
493                                 VecSubf(bezt[a].vec[2], bezt[a].vec[1], vec);
494                         }
495                 }
496         }       
497         /* make the bezier curve */
498         if(cuma->table)
499                 MEM_freeN(cuma->table);
500         totpoint= (cuma->totpoint-1)*CM_RESOL;
501         fp= allpoints= MEM_callocN(totpoint*2*sizeof(float), "table");
502         
503         for(a=0; a<cuma->totpoint-1; a++, fp += 2*CM_RESOL) {
504                 correct_bezpart(bezt[a].vec[1], bezt[a].vec[2], bezt[a+1].vec[0], bezt[a+1].vec[1]);
505                 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); 
506                 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);
507         }
508         
509         /* store first and last handle for extrapolation, unit length */
510         cuma->ext_in[0]= bezt[0].vec[0][0] - bezt[0].vec[1][0];
511         cuma->ext_in[1]= bezt[0].vec[0][1] - bezt[0].vec[1][1];
512         range= sqrt(cuma->ext_in[0]*cuma->ext_in[0] + cuma->ext_in[1]*cuma->ext_in[1]);
513         cuma->ext_in[0]/= range;
514         cuma->ext_in[1]/= range;
515         
516         a= cuma->totpoint-1;
517         cuma->ext_out[0]= bezt[a].vec[1][0] - bezt[a].vec[2][0];
518         cuma->ext_out[1]= bezt[a].vec[1][1] - bezt[a].vec[2][1];
519         range= sqrt(cuma->ext_out[0]*cuma->ext_out[0] + cuma->ext_out[1]*cuma->ext_out[1]);
520         cuma->ext_out[0]/= range;
521         cuma->ext_out[1]/= range;
522         
523         /* cleanup */
524         MEM_freeN(bezt);
525
526         range= CM_TABLEDIV*(cuma->maxtable - cuma->mintable);
527         cuma->range= 1.0f/range;
528         
529         /* now make a table with CM_TABLE equal x distances */
530         fp= allpoints;
531         lastpoint= allpoints + 2*(totpoint-1);
532         cmp= MEM_callocN((CM_TABLE+1)*sizeof(CurveMapPoint), "dist table");
533         
534         for(a=0; a<=CM_TABLE; a++) {
535                 curf= cuma->mintable + range*(float)a;
536                 cmp[a].x= curf;
537                 
538                 /* get the first x coordinate larger than curf */
539                 while(curf >= fp[0] && fp!=lastpoint) {
540                         fp+=2;
541                 }
542                 if(fp==allpoints || (curf >= fp[0] && fp==lastpoint))
543                         cmp[a].y= curvemap_calc_extend(cuma, curf, allpoints, lastpoint);
544                 else {
545                         float fac1= fp[0] - fp[-2];
546                         float fac2= fp[0] - curf;
547                         if(fac1 > FLT_EPSILON)
548                                 fac1= fac2/fac1;
549                         else
550                                 fac1= 0.0f;
551                         cmp[a].y= fac1*fp[-1] + (1.0f-fac1)*fp[1];
552                 }
553         }
554         
555         MEM_freeN(allpoints);
556         cuma->table= cmp;
557 }
558
559 /* call when you do images etc, needs restore too. also verifies tables */
560 /* it uses a flag to prevent premul or free to happen twice */
561 void curvemapping_premultiply(CurveMapping *cumap, int restore)
562 {
563         int a;
564         
565         if(restore) {
566                 if(cumap->flag & CUMA_PREMULLED) {
567                         for(a=0; a<3; a++) {
568                                 MEM_freeN(cumap->cm[a].table);
569                                 cumap->cm[a].table= cumap->cm[a].premultable;
570                                 cumap->cm[a].premultable= NULL;
571                         }
572                         
573                         cumap->flag &= ~CUMA_PREMULLED;
574                 }
575         }
576         else {
577                 if((cumap->flag & CUMA_PREMULLED)==0) {
578                         /* verify and copy */
579                         for(a=0; a<3; a++) {
580                                 if(cumap->cm[a].table==NULL)
581                                         curvemap_make_table(cumap->cm+a, &cumap->clipr);
582                                 cumap->cm[a].premultable= cumap->cm[a].table;
583                                 cumap->cm[a].table= MEM_mallocN((CM_TABLE+1)*sizeof(CurveMapPoint), "premul table");
584                                 memcpy(cumap->cm[a].table, cumap->cm[a].premultable, (CM_TABLE+1)*sizeof(CurveMapPoint));
585                         }
586                         
587                         if(cumap->cm[3].table==NULL)
588                                 curvemap_make_table(cumap->cm+3, &cumap->clipr);
589                 
590                         /* premul */
591                         for(a=0; a<3; a++) {
592                                 int b;
593                                 for(b=0; b<=CM_TABLE; b++) {
594                                         cumap->cm[a].table[b].y= curvemap_evaluateF(cumap->cm+3, cumap->cm[a].table[b].y);
595                                 }
596                         }
597                         
598                         cumap->flag |= CUMA_PREMULLED;
599                 }
600         }
601 }
602
603 static int sort_curvepoints(const void *a1, const void *a2)
604 {
605         const struct CurveMapPoint *x1=a1, *x2=a2;
606         
607         if( x1->x > x2->x ) return 1;
608         else if( x1->x < x2->x) return -1;
609         return 0;
610 }
611
612 /* ************************ more CurveMapping calls *************** */
613
614 /* note; only does current curvemap! */
615 void curvemapping_changed(CurveMapping *cumap, int rem_doubles)
616 {
617         CurveMap *cuma= cumap->cm+cumap->cur;
618         CurveMapPoint *cmp= cuma->curve;
619         rctf *clipr= &cumap->clipr;
620         float thresh= 0.01f*(clipr->xmax - clipr->xmin);
621         float dx= 0.0f, dy= 0.0f;
622         int a;
623         
624         /* clamp with clip */
625         if(cumap->flag & CUMA_DO_CLIP) {
626                 for(a=0; a<cuma->totpoint; a++) {
627                         if(cmp[a].flag & CUMA_SELECT) {
628                                 if(cmp[a].x < clipr->xmin)
629                                         dx= MIN2(dx, cmp[a].x - clipr->xmin);
630                                 else if(cmp[a].x > clipr->xmax)
631                                         dx= MAX2(dx, cmp[a].x - clipr->xmax);
632                                 if(cmp[a].y < clipr->ymin)
633                                         dy= MIN2(dy, cmp[a].y - clipr->ymin);
634                                 else if(cmp[a].y > clipr->ymax)
635                                         dy= MAX2(dy, cmp[a].y - clipr->ymax);
636                         }
637                 }
638                 for(a=0; a<cuma->totpoint; a++) {
639                         if(cmp[a].flag & CUMA_SELECT) {
640                                 cmp[a].x -= dx;
641                                 cmp[a].y -= dy;
642                         }
643                 }
644         }
645         
646         
647         qsort(cmp, cuma->totpoint, sizeof(CurveMapPoint), sort_curvepoints);
648         
649         /* remove doubles, threshold set on 1% of default range */
650         if(rem_doubles && cuma->totpoint>2) {
651                 for(a=0; a<cuma->totpoint-1; a++) {
652                         dx= cmp[a].x - cmp[a+1].x;
653                         dy= cmp[a].y - cmp[a+1].y;
654                         if( sqrt(dx*dx + dy*dy) < thresh ) {
655                                 if(a==0) {
656                                         cmp[a+1].flag|= 2;
657                                         if(cmp[a+1].flag & CUMA_SELECT)
658                                                 cmp[a].flag |= CUMA_SELECT;
659                                 }
660                                 else {
661                                         cmp[a].flag|= 2;
662                                         if(cmp[a].flag & CUMA_SELECT)
663                                                 cmp[a+1].flag |= CUMA_SELECT;
664                                 }
665                                 break;  /* we assume 1 deletion per edit is ok */
666                         }
667                 }
668                 if(a != cuma->totpoint-1)
669                         curvemap_remove(cuma, 2);
670         }       
671         curvemap_make_table(cuma, clipr);
672 }
673
674 /* table should be verified */
675 float curvemap_evaluateF(CurveMap *cuma, float value)
676 {
677         float fi;
678         int i;
679
680         /* index in table */
681         fi= (value-cuma->mintable)*cuma->range;
682         i= (int)fi;
683         
684         /* fi is table float index and should check against table range i.e. [0.0 CM_TABLE] */
685         if(fi<0.0f || fi>CM_TABLE)
686                 return curvemap_calc_extend(cuma, value, &cuma->table[0].x, &cuma->table[CM_TABLE].x);
687         else {
688                 if(i<0) return cuma->table[0].y;
689                 if(i>=CM_TABLE) return cuma->table[CM_TABLE].y;
690                 
691                 fi= fi-(float)i;
692                 return (1.0f-fi)*cuma->table[i].y + (fi)*cuma->table[i+1].y; 
693         }
694 }
695
696 /* works with curve 'cur' */
697 float curvemapping_evaluateF(CurveMapping *cumap, int cur, float value)
698 {
699         CurveMap *cuma= cumap->cm+cur;
700         
701         /* allocate or bail out */
702         if(cuma->table==NULL) {
703                 curvemap_make_table(cuma, &cumap->clipr);
704                 if(cuma->table==NULL)
705                         return value;
706         }
707         return curvemap_evaluateF(cuma, value);
708 }
709
710 /* vector case */
711 void curvemapping_evaluate3F(CurveMapping *cumap, float *vecout, const float *vecin)
712 {
713         vecout[0]= curvemapping_evaluateF(cumap, 0, vecin[0]);
714         vecout[1]= curvemapping_evaluateF(cumap, 1, vecin[1]);
715         vecout[2]= curvemapping_evaluateF(cumap, 2, vecin[2]);
716 }
717
718 /* RGB case, no black/white points, no premult */
719 void curvemapping_evaluateRGBF(CurveMapping *cumap, float *vecout, const float *vecin)
720 {
721         vecout[0]= curvemapping_evaluateF(cumap, 0, curvemapping_evaluateF(cumap, 3, vecin[0]));
722         vecout[1]= curvemapping_evaluateF(cumap, 1, curvemapping_evaluateF(cumap, 3, vecin[1]));
723         vecout[2]= curvemapping_evaluateF(cumap, 2, curvemapping_evaluateF(cumap, 3, vecin[2]));
724 }
725
726
727 /* RGB with black/white points and premult. tables are checked */
728 void curvemapping_evaluate_premulRGBF(CurveMapping *cumap, float *vecout, const float *vecin)
729 {
730         float fac;
731         
732         fac= (vecin[0] - cumap->black[0])*cumap->bwmul[0];
733         vecout[0]= curvemap_evaluateF(cumap->cm, fac);
734         
735         fac= (vecin[1] - cumap->black[1])*cumap->bwmul[1];
736         vecout[1]= curvemap_evaluateF(cumap->cm+1, fac);
737         
738         fac= (vecin[2] - cumap->black[2])*cumap->bwmul[2];
739         vecout[2]= curvemap_evaluateF(cumap->cm+2, fac);
740 }
741
742 void colorcorrection_do_ibuf(ImBuf *ibuf, const char *profile)
743 {
744         if (ibuf->crect == NULL)
745         {
746 #ifdef WITH_LCMS
747                 cmsHPROFILE imageProfile, proofingProfile;
748                 cmsHTRANSFORM hTransform;
749                 
750                 ibuf->crect = MEM_mallocN(ibuf->x*ibuf->y*sizeof(int), "imbuf crect");
751
752                 imageProfile  = cmsCreate_sRGBProfile();
753                 proofingProfile = cmsOpenProfileFromFile(profile, "r");
754                 
755                 cmsErrorAction(LCMS_ERROR_SHOW);
756         
757                 hTransform = cmsCreateProofingTransform(imageProfile, TYPE_RGBA_8, imageProfile, TYPE_RGBA_8, 
758                                                   proofingProfile,
759                                                   INTENT_ABSOLUTE_COLORIMETRIC,
760                                                   INTENT_ABSOLUTE_COLORIMETRIC,
761                                                   cmsFLAGS_SOFTPROOFING);
762         
763                 cmsDoTransform(hTransform, ibuf->rect, ibuf->crect, ibuf->x * ibuf->y);
764         
765                 cmsDeleteTransform(hTransform);
766                 cmsCloseProfile(imageProfile);
767                 cmsCloseProfile(proofingProfile);
768 #else
769                 ibuf->crect = ibuf->rect;
770 #endif
771         }
772 }
773
774
775 void curvemapping_do_ibuf(CurveMapping *cumap, ImBuf *ibuf)
776 {
777         int pixel;
778         
779         if(ibuf==NULL)
780                 return;
781         if(ibuf->rect_float==NULL)
782                 IMB_float_from_rect(ibuf);
783         else if(ibuf->rect==NULL)
784                 imb_addrectImBuf(ibuf);
785         
786         curvemapping_premultiply(cumap, 0);
787         
788         if(ibuf->rect_float && ibuf->rect) {
789                 float *pixf= ibuf->rect_float;
790                 float col[3];
791                 int stride= 4;
792                 char *pixc= (char *)ibuf->rect;
793                 
794                 if(ibuf->channels)
795                         stride= ibuf->channels;
796                 
797                 for(pixel= ibuf->x*ibuf->y; pixel>0; pixel--, pixf+=stride, pixc+=4) {
798                         if(stride<3) {
799                                 col[0]= curvemap_evaluateF(cumap->cm, *pixf);
800                                 pixc[1]= pixc[2]= pixc[3]= pixc[0]= FTOCHAR(col[0]);
801                         }
802                         else {
803                                 curvemapping_evaluate_premulRGBF(cumap, col, pixf);
804                                 pixc[0]= FTOCHAR(col[0]);
805                                 pixc[1]= FTOCHAR(col[1]);
806                                 pixc[2]= FTOCHAR(col[2]);
807                                 if(stride>3)
808                                         pixc[3]= FTOCHAR(pixf[3]);
809                                 else
810                                         pixc[3]= 255;
811                         }
812                 }
813         }
814         
815         curvemapping_premultiply(cumap, 1);
816 }
817
818 int curvemapping_RGBA_does_something(CurveMapping *cumap)
819 {
820         int a;
821         
822         if(cumap->black[0]!=0.0f) return 1;
823         if(cumap->black[1]!=0.0f) return 1;
824         if(cumap->black[2]!=0.0f) return 1;
825         if(cumap->white[0]!=1.0f) return 1;
826         if(cumap->white[1]!=1.0f) return 1;
827         if(cumap->white[2]!=1.0f) return 1;
828         
829         for(a=0; a<CM_TOT; a++) {
830                 if(cumap->cm[a].curve) {
831                         if(cumap->cm[a].totpoint!=2)  return 1;
832                         
833                         if(cumap->cm[a].curve[0].x != 0.0f) return 1;
834                         if(cumap->cm[a].curve[0].y != 0.0f) return 1;
835                         if(cumap->cm[a].curve[1].x != 1.0f) return 1;
836                         if(cumap->cm[a].curve[1].y != 1.0f) return 1;
837                 }
838         }
839         return 0;
840 }
841
842 void curvemapping_initialize(CurveMapping *cumap)
843 {
844         int a;
845         
846         if(cumap==NULL) return;
847         
848         for(a=0; a<CM_TOT; a++) {
849                 if(cumap->cm[a].table==NULL)
850                         curvemap_make_table(cumap->cm+a, &cumap->clipr);
851         }
852 }
853
854 void curvemapping_table_RGBA(CurveMapping *cumap, float **array, int *size)
855 {
856         int a;
857         
858         *size = CM_TABLE+1;
859         *array = MEM_callocN(sizeof(float)*(*size)*4, "CurveMapping");
860         curvemapping_initialize(cumap);
861
862         for(a=0; a<*size; a++) {
863                 if(cumap->cm[0].table)
864                         (*array)[a*4+0]= cumap->cm[0].table[a].y;
865                 if(cumap->cm[1].table)
866                         (*array)[a*4+1]= cumap->cm[1].table[a].y;
867                 if(cumap->cm[2].table)
868                         (*array)[a*4+2]= cumap->cm[2].table[a].y;
869                 if(cumap->cm[3].table)
870                         (*array)[a*4+3]= cumap->cm[3].table[a].y;
871         }
872 }
873