6236708c9e13754abb95ef4340de92c1a4f77743
[blender.git] / source / blender / compositor / operations / COM_LensGhostOperation.cpp
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor: 
19  *              Jeroen Bakker 
20  *              Monique Dewanchand
21  */
22
23 #include "COM_LensGhostOperation.h"
24 #include "BLI_math.h"
25 #include "BLI_utildefines.h"
26
27 #define MAX_STEP 256
28 class Ray {
29 public:
30         float position[3];
31         float direction[3];
32         float uv[2];
33         double wavelength;
34         float intensity;
35         bool valid;
36         void copyFrom(Ray* other) {
37                 copy_v3_v3(position, other->position);
38                 copy_v3_v3(direction, other->direction);
39                 copy_v2_v2(uv, other->uv);
40                 wavelength = other->wavelength;
41                 intensity = other->intensity;
42                 this->valid = other->valid;
43         }
44 };
45
46 class Intersection {
47 public:
48         float position[3];
49         float normal[3];
50         double theta;
51         bool hit;
52         bool inverted;
53 };
54
55 class LensInterface {
56 public:
57         float position[3];
58         float radius;
59         float nominalRadius;
60         double refraction1;
61         double refraction2;
62         double refraction3;
63         float thicknessCoathing;
64         virtual bool isFlat() = 0;
65         virtual void intersect(Intersection* result, Ray* ray) = 0;
66 };
67
68 class FlatInterface: public LensInterface {
69 public:
70         bool isFlat() {return true;}
71         FlatInterface(float positionX, float positionY, float positionZ, float radius) {
72                 this->position[0] = positionX;
73                 this->position[1] = positionY;
74                 this->position[2] = positionZ;
75                 this->radius = radius;
76                 this->nominalRadius = radius;
77                 this->refraction1 = 1.0f;
78                 this->refraction2 = 1.0f;
79                 this->refraction3 = 1.0f;
80                 this->thicknessCoathing = 0.0f;
81
82         }
83         void intersect(Intersection* result, Ray* ray) {
84                 const float dz = this->position[2]-ray->position[2];
85                 result->position[0] = ray->position[0] + ray->direction[0]*(dz)/ray->direction[2];
86                 result->position[1] = ray->position[1] + ray->direction[1]*(dz)/ray->direction[2];
87                 result->position[2] = ray->position[2] + ray->direction[2]*(dz)/ray->direction[2];
88                 result->normal[0] = 0.0f;
89                 result->normal[1] = 0.0f;
90                 result->normal[2] = ray->direction[2]>0?-1.0f:1.0f;
91                 result->theta = 0.0f;
92 //              result->hit = this->nominalRadius>maxf(fabs(result->position[0]), fabs(result->position[1]));
93                 result->hit = true;
94                 result->inverted = false;
95         }
96 };
97
98 class SphereInterface: public LensInterface {
99 public:
100         SphereInterface(float positionX, float positionY, float positionZ, float radius, float nominalRadius, float n0, float n2, float coatingPhase) {
101                 this->position[0] = positionX;
102                 this->position[1] = positionY;
103                 this->position[2] = positionZ;
104                 this->radius = radius;
105                 this->nominalRadius = nominalRadius;
106                 this->refraction1 = n0;
107                 this->refraction3 = n2;
108                 this->refraction2 = maxf(sqrtf(n0*n2), 1.38);
109
110                 this->thicknessCoathing = coatingPhase/4/this->refraction2;
111         }
112         bool isFlat() {return false;}
113         void intersect(Intersection* result, Ray* ray) {
114                 float delta[3] ={ray->position[0] - this->position[0],
115                         ray->position[1] - this->position[1],
116                         ray->position[2] - this->position[2]};
117                 float b = dot_v3v3(delta, ray->direction);
118                 float c = dot_v3v3(delta, delta) - this->radius*this->radius;
119                 float b2c = b*b-c;
120                 if (b2c < 0) {
121                         result->hit = false;
122                 } else {
123                         float sgn = (this->radius*ray->direction[2])>0?1.0f:-1.0f;
124                         float t = sqrtf(b2c)*sgn-b;
125                         result->position[0] = ray->direction[0]*t+ray->position[0];
126                         result->position[1] = ray->direction[1]*t+ray->position[1];
127                         result->position[2] = ray->direction[2]*t+ray->position[2];
128
129                         float p[3] = {
130                                 result->position[0] - this->position[0],
131                                 result->position[1] - this->position[1],
132                                 result->position[2] - this->position[2]
133                         };
134                         normalize_v3(p);
135
136                         if (dot_v3v3(p, ray->direction)> 0) {
137                                 result->normal[0] = -p[0];
138                                 result->normal[1] = -p[1];
139                                 result->normal[2] = -p[2];
140                         } else {
141                                 result->normal[0] = p[0];
142                                 result->normal[1] = p[1];
143                                 result->normal[2] = p[2];
144                         }
145
146                         float inverse[3] ={
147                                 -ray->direction[0],
148                                 -ray->direction[1],
149                                 -ray->direction[2]};
150
151                         result->theta = acosf(dot_v3v3(inverse, result->normal));
152                         result->hit = this->nominalRadius>sqrt(result->position[0]*result->position[0]+result->position[1]*result->position[1]);
153 //                      result->hit = this->nominalRadius>maxf(fabs(result->position[0]), fabs(result->position[1]));
154 //                      result->hit = true;
155                         result->inverted = t < 0;
156                 }
157         }
158 };
159 class RayResult {
160 public:
161         float x;
162         float y;
163         float intensity[3];
164         float u;
165         float v;
166         float screenX;
167         float screenY;
168         bool valid;
169         bool hasIntensity;
170 };
171 class Bounce{
172 public:
173         LensInterface *interface1;
174         LensInterface *interface2;
175         RayResult *raster;
176         int length; // number of interfaces to travel
177         int rasterLength;
178         Bounce(LensInterface *interface1, LensInterface *interface2, int length, int rasterStep) {
179                 this->interface1 = interface1;
180                 this->interface2 = interface2;
181                 this->length = length;
182                 this->rasterLength = rasterStep;
183                 this->raster = new RayResult[rasterLength*rasterLength];
184                 for (int i = 0 ; i < rasterLength*rasterLength ; i++) {
185                         RayResult * res = &this->raster[i];
186                         res->intensity[0] = 0.0f;
187                         res->intensity[1] = 0.0f;
188                         res->intensity[2] = 0.0f;
189                         res->x = 0.0f;
190                         res->y = 0.0f;
191                         res->u = 0.0f;
192                         res->v = 0.0f;
193                         res->valid = false;
194                 }
195         }
196         ~Bounce() {
197                 delete raster;
198
199         }
200
201         RayResult* getRayResult(int x, int y) {
202                 return &(raster[x+y*rasterLength]);
203         }
204 };
205 class LensSystem {
206 public:
207         vector<LensInterface*> interfaces;
208         vector<Bounce*> bounces;
209         int bokehIndex;
210         int lensIndex;
211
212         ~LensSystem() {
213                 for (int index = 0 ; index <bounces.size();index++) {delete bounces[index];}
214                 for (int index = 0 ; index <interfaces.size();index++) {delete interfaces[index];}
215         }
216
217         void updateBounces(int step) {
218                 for (int i = 0; i < interfaces.size()-1 ; i ++) {
219                         if (!interfaces[i]->isFlat()) {
220                                 for (int j = i+1; j < interfaces.size()-1 ; j ++) {
221                                         if (!interfaces[j]->isFlat()) {
222                                                 int length = interfaces.size()+2*(j-i);
223                                                 Bounce* bounce = new Bounce(interfaces[j], interfaces[i], length, step);
224                                                 bounces.push_back(bounce);
225                                         }
226                                 }
227                         }
228                 }
229
230         }
231
232         void addInterface(LensInterface* pinterface) {
233                 this->interfaces.push_back(pinterface);
234                 this->lensIndex = this->interfaces.size()-1;
235         }
236
237         static int refraction(float *refract, float *n, float *view, double index)
238         {
239
240                 return 1;
241
242 //                float dot, fac;
243
244 //                VECCOPY(refract, view);
245
246 //                dot= view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
247
248 //                if(dot>0.0f) {
249 //                        index = 1.0f/index;
250 //                        fac= 1.0f - (1.0f - dot*dot)*index*index;
251 //                        if(fac<= 0.0f) return 0;
252 //                        fac= -dot*index + sqrt(fac);
253 //                }
254 //                else {
255 //                        fac= 1.0f - (1.0f - dot*dot)*index*index;
256 //                        if(fac<= 0.0f) return 0;
257 //                        fac= -dot*index - sqrt(fac);
258 //                }
259
260 //                refract[0]= index*view[0] + fac*n[0];
261 //                refract[1]= index*view[1] + fac*n[1];
262 //                refract[2]= index*view[2] + fac*n[2];
263
264 //                normalize_v3(refract);
265 //                return 1;
266                 //---
267 //              const double cosI = dot_v3v3(n, view);
268 //              const double sinT2 = index * index * (1.0 - cosI * cosI);
269 //              if (sinT2 >= 1.0f)
270 //              {
271 //                      return 0;
272 //              }
273 //              refract[0] = index*view[0] - (index + sqrt(1.0-sinT2))*n[0];
274 //              refract[1] = index*view[1] - (index + sqrt(1.0-sinT2))*n[1];
275 //              refract[2] = index*view[2] - (index + sqrt(1.0-sinT2))*n[2];
276 //              normalize_v3(refract);
277 //              return 1;
278                 //---
279
280 //                double ni = -dot_v3v3(view, n);
281 //                double test = 1.0f - index*index*(1.0f-ni*ni);
282 //                if (test < 0) {
283 //                        return 0;
284 //                } else {
285 //                        double mul = index*ni + sqrt(test);
286 //                        refract[0] = index * view[0] - mul*n[0];
287 //                        refract[1] = index * view[1] - mul*n[1];
288 //                        refract[2] = index * view[2] - mul*n[2];
289 //                        normalize_v3(refract);
290 //                        return 1;
291 //                }
292         }
293
294         /* orn = original face normal */
295         static void reflection(float *ref, float *n, float *view)
296         {
297                 float f1;
298
299                 f1= -2.0f*dot_v3v3(n, view);
300
301                 ref[0]= (view[0]+f1*n[0]);
302                 ref[1]= (view[1]+f1*n[1]);
303                 ref[2]= (view[2]+f1*n[2]);
304                 normalize_v3(ref);
305         }
306
307         static float fresnelAR(float theta0, float lambda, float d1, float n0, float n1, float n2) {
308                 // refractionangles in coating and the 2nd medium
309                 float theta1 = asin(sin(theta0)*n0/n1);
310                 float theta2 = asin(sin(theta0)*n0/n2);
311
312                 float rs01 = -sin(theta0-theta1)/sin(theta0+theta1);
313                 float rp01 = tan( theta0-theta1)/tan(theta0+theta1);
314                 float ts01 = 2 * sin ( theta1 ) * cos ( theta0 ) / sin ( theta0+theta1 ) ;
315                 float tp01 = ts01*cos(theta0-theta1);
316                 // amplitude for inner reflection
317                 float rs12 = -sin ( theta1-theta2 ) / sin ( theta1+theta2 ) ;
318                 float rp12 = +tan ( theta1-theta2 ) / tan ( theta1+theta2 ) ;
319                 // after passing through first surface twice :
320                 // 2 transmissions and 1 reflection
321                 float ris = ts01 * ts01 * rs12 ;
322                 float rip = tp01 * tp01 * rp12 ;
323                 // phase difference between outer and inner reflections
324                 float dy = d1 * n1 ;
325                 float dx = tan ( theta1 ) * dy ;
326                 float delay = sqrt ( dx * dx+dy * dy ) ;
327                 float relPhase = 4 * M_PI / lambda * ( delay-dx * sin ( theta0 ) ) ;
328                 // Add up sines of different phase and amplitude
329                 float out_s2 = rs01 * rs01 + ris * ris + 2 * rs01 * ris * cos ( relPhase ) ;
330                 float out_p2 = rp01 * rp01 + rip * rip + 2 * rp01 * rip * cos ( relPhase ) ;
331                 return ( out_s2+out_p2 ) / 2 ;
332         }
333
334         void detectHit(Ray* result, Ray* inputRay, Bounce *bounce) {
335                 int phase = 0;
336                 int delta = 1;
337                 int t = 1;
338                 int k;
339                 result->copyFrom(inputRay);
340                 result->valid = false;
341                 LensInterface* next = bounce->interface1;
342                 LensInterface* f = NULL;
343                 Intersection intersection;
344                 for (k = 0 ; k < bounce->length-1;k++, t+=delta) {
345                         f = this->interfaces[t];
346                         bool breflect = next == f;
347                         if (breflect) {
348                                 delta = -delta;
349                                 if (phase == 0) {
350                                         next = bounce->interface2;
351                                 } else {
352                                         next = NULL;
353                                 }
354                                 phase ++;
355                         }
356
357                         f->intersect(&intersection, result);
358                         if (!intersection.hit) {
359                                 break;
360                         }
361                         if (f->isFlat()) {
362                                 if (t == this->bokehIndex) {
363                                         result->uv[0] = intersection.position[0]/f->nominalRadius;
364                                         result->uv[1] = intersection.position[1]/f->nominalRadius;
365                                 }
366                         }
367
368                         float p[3] = {
369                                 intersection.position[0]-result->position[0],
370                                 intersection.position[1]-result->position[1],
371                                 intersection.position[2]-result->position[2]
372                         };
373
374                         float nfac = sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
375
376                         if (intersection.inverted) {
377                                 nfac *= -1;
378                         }
379
380                         result->direction[0] = p[0]/nfac;
381                         result->direction[1] = p[1]/nfac;
382                         result->direction[2] = p[2]/nfac;
383                         result->position[0] = intersection.position[0];
384                         result->position[1] = intersection.position[1];
385                         result->position[2] = intersection.position[2];
386
387                         if (!f->isFlat()) {
388                                 // do refraction and reflection
389                                 double n0 = result->direction[2]<0?f->refraction1:f->refraction3;
390                                 double n1 = f->refraction2;
391                                 double n2 = result->direction[2]<0?f->refraction3:f->refraction1;
392                                 if (!breflect) {
393                                         float view[3] ={
394                                                 result->direction[0],
395                                                 result->direction[1],
396                                                 result->direction[2]
397                                         };
398                                         int ref = this->refraction(result->direction, intersection.normal, view, n0/n1);
399                                         if (ref == 0) {
400                                                 break;
401                                         }
402                                 } else {
403                                         this->reflection(result->direction, intersection.normal, result->direction);
404                                         float  fresnelMultiplyer = fresnelAR(intersection.theta, result->wavelength, f->thicknessCoathing, n0, n1, n2);
405                                         if (isnan(fresnelMultiplyer)) {
406                                                 fresnelMultiplyer = 0.0f;
407                                         }
408                                         result->intensity *= fresnelMultiplyer;
409                                 }
410                         }
411
412                 }
413                 if (k < bounce->length-1) {
414                         result->intensity = 0;
415                 } else {
416                         result->valid = true;
417                 }
418         }
419 };
420
421 typedef struct LensFace {
422         RayResult* v1;
423         RayResult* v2;
424         RayResult* v3;
425 } LensFace;
426
427 LensGhostProjectionOperation::LensGhostProjectionOperation(): NodeOperation() {
428         this->addInputSocket(COM_DT_COLOR);
429         this->addInputSocket(COM_DT_COLOR, COM_SC_NO_RESIZE);
430         this->addOutputSocket(COM_DT_COLOR);
431         this->lampObject = NULL;
432         this->cameraObject = NULL;
433         this->system = NULL;
434         this->quality = COM_QUALITY_HIGH;
435         this->setComplex(false);
436 }
437
438 LensGhostOperation::LensGhostOperation(): LensGhostProjectionOperation() {
439         this->setComplex(true);
440
441 }
442
443 void LensGhostProjectionOperation::initExecution() {
444         if (this->cameraObject != NULL && this->lampObject != NULL) {
445                 if (lampObject == NULL || cameraObject == NULL) {
446                         visualLampPosition[0] = 0;
447                         visualLampPosition[1] = 0;
448                         visualLampPosition[2] = 0;
449                 } else {
450                         /* too simple, better to return the distance on the view axis only
451                          * return len_v3v3(ob->obmat[3], cam->dof_ob->obmat[3]); */
452                         float matt[4][4], imat[4][4], obmat[4][4];
453
454                         copy_m4_m4(obmat, cameraObject->obmat);
455                         normalize_m4(obmat);
456                         invert_m4_m4(imat, obmat);
457                         mult_m4_m4m4(matt, imat, lampObject->obmat);
458
459                         visualLampPosition[0] = (float)(matt[3][0]);
460                         visualLampPosition[1] = (float)(matt[3][1]);
461                         visualLampPosition[2] = (float)fabs(matt[3][2]);
462                 }
463         }
464         this->lamp = (Lamp*)lampObject->data;
465
466         this->step = this->quality==COM_QUALITY_LOW?64:this->quality==COM_QUALITY_MEDIUM?128:256;
467         this->bokehReader = this->getInputSocketReader(1);
468
469 #define MM *0.001f
470 #define CM *0.01f
471 #define NM *0.000000001
472 #define RED 650 NM
473 #define GREEN 510 NM
474 #define BLUE 475 NM
475 #define AIR 1.000293f
476 #define GLASS 1.5200f
477 #define TEST 0.000002f
478         // determine interfaces
479         LensSystem *system = new LensSystem();
480         system->addInterface(new FlatInterface(0.0f,0.0f, 6.5 CM, 30 MM)); //ENTRANCE
481         system->addInterface(new SphereInterface(0.0f,0.0f, -3 CM, 8 CM, 3 CM, AIR, GLASS , 0.f));
482         system->addInterface(new SphereInterface(0.0f,0.0f, -4 CM, 8 CM, 3 CM, GLASS, AIR, GREEN));
483         system->addInterface(new FlatInterface(0.0f,0.0f, 3.0 CM, 15 MM)); // BOKEH
484         system->addInterface(new SphereInterface(0.0f,0.0f, 6 CM, 3 CM, 2 CM, AIR, GLASS, 0.0f));
485         system->addInterface(new SphereInterface(0.0f,0.0f, 5.5 CM, 3 CM, 2 CM, GLASS, AIR, 0.f));
486         system->addInterface(new FlatInterface(0.0f,0.0f,0 CM, 30 MM)); // SENSOR
487         system->bokehIndex =3;
488
489         // determine interfaces
490 //      LensSystem *system = new LensSystem();
491 //      system->addInterface(new FlatInterface(0.0f,0.0f, 6.5 CM, 30 MM)); //ENTRANCE
492 //      system->addInterface(new SphereInterface(0.0f,0.0f, 14 CM, 8 CM, 6 CM, AIR, GLASS , 0.0f));
493 //      system->addInterface(new SphereInterface(0.0f,0.0f, 12 CM, 8 CM, 6 CM, GLASS, AIR, GREEN));
494 //      system->addInterface(new FlatInterface(0.0f,0.0f, 3.0 CM, 30 MM)); // BOKEH
495 //      system->addInterface(new SphereInterface(0.0f,0.0f, 1 CM, 3 CM, 2 CM, AIR, GLASS, GREEN));
496 //      system->addInterface(new SphereInterface(0.0f,0.0f, -2 CM, 3 CM, 2 CM, GLASS, AIR, RED));
497 //      system->addInterface(new FlatInterface(0.0f,0.0f,0 CM, 20 MM)); // SENSOR
498 //      system->bokehIndex = 3;
499 #undef CM
500 #undef MM
501         // determine bounces
502         system->updateBounces(step);
503         this->system = system;
504 }
505
506 void LensGhostOperation::initExecution() {
507         LensGhostProjectionOperation::initExecution();
508         LensSystem *system = (LensSystem*)this->system;
509         LensInterface *interface1 = system->interfaces[0];
510
511         // for every herz
512         float HERZ[3]={650 NM,510 NM,475 NM}; /// @todo use 7 for high quality?
513         for (int iw = 0 ; iw < 3 ; iw ++) {
514                 float wavelength = HERZ[iw];
515         // for every bounce
516                 for (int ib = 0 ; ib < system->bounces.size() ; ib++) {
517                         Bounce* bounce = system->bounces[ib];
518         // based on quality setting the number of iteration will be different (128^2, 64^2, 32^2)
519                         for (int xi = 0 ; xi < step ; xi ++) {
520                                 float x = -interface1->radius+xi*(interface1->radius*2/step);
521                                 for (int yi = 0 ; yi < step ; yi ++) {
522                                         float y = -interface1->radius+yi*(interface1->radius*2/step);
523                                         Ray r;
524                                         Ray result;
525                                         r.wavelength = wavelength;
526                                         r.intensity = this->lamp->energy;
527                                         r.uv[0] = 0.0f;
528                                         r.uv[1] = 0.0f;
529                                         r.position[0] = visualLampPosition[0];
530                                         r.position[1] = visualLampPosition[1];
531                                         r.position[2] = visualLampPosition[2];
532                                         r.direction[0] = interface1->position[0]+x - r.position[0];
533                                         r.direction[1] = interface1->position[1]+y - r.position[1];
534                                         r.direction[2] = interface1->position[2] - r.position[2];
535                                         normalize_v3(r.direction);
536                                         system->detectHit(&result, &r, bounce);
537                                         RayResult *res = bounce->getRayResult(xi, yi);
538                                         if (iw == 0) {
539                                                 res->x = result.position[0];
540                                                 res->y = result.position[1];
541                                                 res->u = result.uv[0];
542                                                 res->v = result.uv[1];
543                                         }
544                                         res->intensity[iw] = result.intensity;
545                                         if (result.valid) {
546                                                 res->valid = true;
547                                         }
548                                 }
549                         }
550                 }
551         }
552 #undef NM
553         const int width = this->getWidth();
554         const int height = this->getHeight();
555         const float width2 = width/2.0f;
556         const float height2 = height/2.0f;
557         float *data = new float[width*height*4];
558         for (int i = 0 ; i < width*height ; i ++) {
559                 data[i*4+0] = 0.0f;
560                 data[i*4+1] = 0.0f;
561                 data[i*4+2] = 0.0f;
562                 data[i*4+3] = 1.0f;
563         }
564         /// @todo every bounce creates own image. these images are added together at the end
565 //      LensSystem *system = (LensSystem*)this->system;
566         LensInterface * lens = system->interfaces[system->lensIndex];
567         for (int i = 0 ; i < system->bounces.size() ; i ++) {
568                 Bounce* bounce = system->bounces[i];
569                 for (int r = 0 ; r < bounce->rasterLength*bounce->rasterLength ; r ++) {
570                         RayResult *result = &bounce->raster[r];
571 //                      if (result->valid) {
572                                 float ru= result->x/lens->nominalRadius*width2+width2;
573                                 float rv  = result->y/lens->nominalRadius*height2+height2;
574                                 result->screenX = ru;
575                                 result->screenY = rv;
576                                 result->hasIntensity = result->intensity[0]>0.0f &&result->intensity[1]>0.0f&& result->intensity[2]>0.0f;
577 //                      }
578                 }
579         }
580 }
581
582 void* LensGhostOperation::initializeTileData(rcti *rect, MemoryBuffer **memoryBuffers) {
583         vector<LensFace*>* result = new vector<LensFace*>();
584         LensSystem *system = (LensSystem*)this->system;
585         const float minx = rect->xmin;
586         const float miny = rect->ymin;
587         const float maxx = rect->xmax;
588         const float maxy = rect->ymax;
589         for (int i = 0 ; i < system->bounces.size() ; i ++) {
590                 Bounce* bounce = system->bounces[i];
591                 int faceX, faceY;
592                 for (faceX = 0 ; faceX < bounce->rasterLength-1 ; faceX++) {
593                         for (faceY = 0 ; faceY < bounce->rasterLength-1 ; faceY++) {
594                                 RayResult* vertex1 = bounce->getRayResult(faceX, faceY);
595                                 RayResult* vertex2 = bounce->getRayResult(faceX+1, faceY);
596                                 RayResult* vertex3 = bounce->getRayResult(faceX+1, faceY+1);
597                                 RayResult* vertex4 = bounce->getRayResult(faceX, faceY+1);
598                                 // early hit test
599                                 if (!((vertex1->screenX < minx && vertex2->screenX < minx && vertex3->screenX < minx && vertex4->screenX < minx) ||
600                                         (vertex1->screenX > maxx && vertex2->screenX > maxx && vertex3->screenX > maxx && vertex4->screenX > maxx) ||
601                                         (vertex1->screenY < miny && vertex2->screenY < miny && vertex3->screenY < miny && vertex4->screenY < miny) ||
602                                         (vertex1->screenY > maxy && vertex2->screenY > maxy && vertex3->screenY > maxy && vertex4->screenY > maxy))) {
603                                         int number = vertex1->hasIntensity +vertex2->hasIntensity +vertex3->hasIntensity +vertex4->hasIntensity;
604                                         if (number == 4) {
605                                                 LensFace* face = new LensFace();
606                                                 face->v1 = vertex1;
607                                                 face->v2 = vertex2;
608                                                 face->v3 = vertex3;
609                                                 result->push_back(face);
610                                                 face = new LensFace();
611                                                 face->v1 = vertex3;
612                                                 face->v2 = vertex4;
613                                                 face->v3 = vertex1;
614                                                 result->push_back(face);
615                                         } else if (number == 3) {
616                                                 LensFace *face = new LensFace();
617                                                 if (!vertex1->hasIntensity) {
618                                                         face->v1 = vertex2;
619                                                         face->v2 = vertex3;
620                                                         face->v3 = vertex4;
621                                                 } else if (!vertex2->hasIntensity) {
622                                                         face->v1 = vertex1;
623                                                         face->v2 = vertex3;
624                                                         face->v3 = vertex4;
625                                                 } else if (!vertex3->hasIntensity) {
626                                                         face->v1 = vertex1;
627                                                         face->v2 = vertex2;
628                                                         face->v3 = vertex4;
629                                                 } else {
630                                                         face->v1 = vertex1;
631                                                         face->v2 = vertex2;
632                                                         face->v3 = vertex3;
633                                                 }
634                                                 result->push_back(face);
635                                         }
636                                 }
637                         }
638                 }
639         }
640
641         return result;
642 }
643
644 void LensGhostOperation::deinitializeTileData(rcti *rect, MemoryBuffer **memoryBuffers, void *data) {
645         if (data) {
646                 vector<LensFace*>* faces = (vector<LensFace*>*)data;
647                 while (faces->size() != 0) {
648                         LensFace *face = faces->back();
649                         faces->pop_back();
650                         delete face;
651                 }
652                 delete faces;
653         }
654 }
655
656
657 void LensGhostProjectionOperation::executePixel(float* color, float x, float y, PixelSampler sampler, MemoryBuffer *inputBuffers[]) {
658         float bokeh[4];
659         LensSystem *system = (LensSystem*)this->system;
660         LensInterface *interface1 = system->interfaces[0];
661         color[0] = 0.0f;
662         color[1] = 0.0f;
663         color[2] = 0.0f;
664         color[3] = 0.0f;
665         const float width = this->getWidth();
666         const float height = this->getHeight();
667         const float size = min(height, width);
668         const float width2 = width/2;
669         const float height2 = height/2;
670         const float size2 = size/2;
671
672 #define NM *0.000000001
673         float HERZ[3]={650 NM,510 NM,475 NM}; /// @todo use 7 for high quality?
674         float rx = ((x-width2)/size2) * interface1->radius;
675         float ry = ((y-height2)/size2) * interface1->radius;
676
677         for (int iw = 0 ; iw < 3 ; iw ++) {
678                 float intensity = 0.0f;
679                 float wavelength = HERZ[iw];
680                 float colorcomponent = 0.0f;
681                 if (iw ==0 ) colorcomponent = lamp->r;
682                 if (iw ==1 ) colorcomponent = lamp->g;
683                 if (iw ==2 ) colorcomponent = lamp->b;
684
685
686         // for every bounce
687                 for (int ib = 0 ; ib < system->bounces.size() ; ib++) {
688                         Bounce* bounce = system->bounces[ib];
689         // based on quality setting the number of iteration will be different (128^2, 64^2, 32^2)
690
691                         Ray r;
692                         Ray result;
693                         r.wavelength = wavelength;
694                         r.intensity = this->lamp->energy;
695                         r.uv[0] = 0.0f;
696                         r.uv[1] = 0.0f;
697                         r.position[0] = visualLampPosition[0];
698                         r.position[1] = visualLampPosition[1];
699                         r.position[2] = visualLampPosition[2];
700                         r.direction[0] = interface1->position[0]+rx - r.position[0];
701                         r.direction[1] = interface1->position[1]+ry - r.position[1];
702                         r.direction[2] = interface1->position[2] - r.position[2];
703                         normalize_v3(r.direction);
704                         system->detectHit(&result, &r, bounce);
705                         if (result.valid) {
706                                 float u = ((result.uv[0]+1.0f)/2)*bokehReader->getWidth();
707                                 float v = ((result.uv[1]+1.0f)/2)*bokehReader->getHeight();
708
709                                 bokehReader->read(bokeh, u, v, sampler, inputBuffers);
710
711                                 intensity += result.intensity *bokeh[iw];
712                         }
713                 }
714                 intensity = maxf(0.0f, intensity);
715                 color[iw] = intensity*colorcomponent;
716         }
717         color[3] = 1.0f;
718 #undef NM
719
720 }
721
722
723
724 void LensGhostOperation::executePixel(float* color, int x, int y, MemoryBuffer *inputBuffers[], void* data) {
725         vector<LensFace*>* faces = (vector<LensFace*>*)data;
726 #if 0 /* UNUSED */
727         const float bokehWidth = bokehReader->getWidth();
728         const float bokehHeight = bokehReader->getHeight();
729         float bokeh[4];
730 #endif
731         color[0] = 0.0f;
732         color[1] = 0.0f;
733         color[2] = 0.0f;
734         color[3] = 1.0f;
735
736         unsigned int index;
737         for (index = 0 ; index < faces->size() ; index ++) {
738                 LensFace * face = faces->operator [](index);
739                 RayResult* vertex1 = face->v1;
740                 RayResult* vertex2 = face->v2;
741                 RayResult* vertex3 = face->v3;
742                 if (!((vertex1->screenX < x && vertex2->screenX < x && vertex3->screenX < x) ||
743                         (vertex1->screenX > x && vertex2->screenX > x && vertex3->screenX > x) ||
744                         (vertex1->screenY < y && vertex2->screenY < y && vertex3->screenY < y) ||
745                         (vertex1->screenY > y && vertex2->screenY > y && vertex3->screenY > y))) {
746
747                         const float v1[2] = {vertex1->screenX, vertex1->screenY};
748                         const float v2[2] = {vertex2->screenX, vertex2->screenY};
749                         const float v3[2] = {vertex3->screenX, vertex3->screenY};
750                         const float co[2] = {x, y};
751                         float weights[3];
752
753                         barycentric_weights_v2(v1, v2, v3, co, weights);
754                         if (weights[0]>=0.0f && weights[0]<=1.0f &&
755                                         weights[1]>=0.0f && weights[1]<=1.0f &&
756                                         weights[2]>=0.0f && weights[2]<=1.0f) {
757 //                            const float u = (vertex1->u*weights[0]+vertex2->u*weights[1]+vertex3->u*weights[2]);
758 //                            const float v = (vertex1->v*weights[0]+vertex2->v*weights[1]+vertex3->v*weights[2]);
759 //                            const float tu = ((u+1.0f)/2.0f)*bokehWidth;
760 //                            const float tv = ((v+1.0f)/2.0f)*bokehHeight;
761 //                            bokehReader->read(bokeh, tu, tv, inputBuffers);
762
763 //                                                      color[0] = max(color[0], bokeh[0]*(vertex1->intensity[0]*weights[0]+vertex2->intensity[0]*weights[1]+vertex3->intensity[0]*weights[2]));
764 //                            color[1] = max(color[1], bokeh[1]*(vertex1->intensity[1]*weights[0]+vertex2->intensity[1]*weights[1]+vertex3->intensity[1]*weights[2]));
765 //                            color[2] = max(color[2], bokeh[2]*(vertex1->intensity[2]*weights[0]+vertex2->intensity[2]*weights[1]+vertex3->intensity[2]*weights[2]));
766                                                         color[0] = max(color[0], (vertex1->intensity[0]*weights[0]+vertex2->intensity[0]*weights[1]+vertex3->intensity[0]*weights[2]));
767                             color[1] = max(color[1], (vertex1->intensity[1]*weights[0]+vertex2->intensity[1]*weights[1]+vertex3->intensity[1]*weights[2]));
768                             color[2] = max(color[2], (vertex1->intensity[2]*weights[0]+vertex2->intensity[2]*weights[1]+vertex3->intensity[2]*weights[2]));
769                         }
770                 }
771         }
772 }
773
774
775 void LensGhostProjectionOperation::deinitExecution() {
776         if (this->system) delete (LensSystem*)this->system;
777         this->system = NULL;
778         this->bokehReader = NULL;
779 }
780
781 bool LensGhostProjectionOperation::determineDependingAreaOfInterest(rcti *input, ReadBufferOperation *readOperation, rcti *output) {
782         rcti bokehInput;
783
784         NodeOperation *operation = this->getInputOperation(1);
785         bokehInput.xmax = operation->getWidth();
786         bokehInput.xmin = 0;
787         bokehInput.ymax = operation->getHeight();
788         bokehInput.ymin = 0;
789         if (operation->determineDependingAreaOfInterest(&bokehInput, readOperation, output) ) {
790                 return true;
791         }
792
793         return NodeOperation::determineDependingAreaOfInterest(input, readOperation, output);
794 }