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