In [446]:

    

# raytracing tutorial
# 16 - light fall off with distance


    

In [447]:

    

import numpy
import matplotlib.pyplot as plt
import math

from opensimplex import OpenSimplex

# plot images in this notebook
%matplotlib inline


    

In [448]:

    

# axes x to the right, y upwards. z into the screen (left hand rule)


    

In [449]:

    

# texture function, takes coordinates, returns colour

noise_generator = OpenSimplex()

def texture_noise(location):
    x = location[0]
    y = location[1]
    z = location[2]
    
    colour = numpy.array([noise_generator.noise3d(x*10,y*10,z*10) + 1.0,
                          noise_generator.noise3d(x*10,y*10,z*10) + 1.0,
                          noise_generator.noise3d(x*10,y*10,z*10) + 1.0])
    
    return colour


    

In [450]:

    

# texture based on opemsimplex noise

def texture_wavy(location):
    x = location[0]
    y = location[1]
    z = location[2]
    
    #colour = (numpy.array([math.sin(x*2), 0, 0]) + 1) / 2.0
    #colour = (numpy.array([math.sin(y*2), math.sin(y*2), 0]) + 1) / 2.0
    colour = (numpy.array([math.sin((x*3) + math.sin(y*4) + math.sin(z*4)), 0, 0]) + 1) / 2.0
    
    return colour


    

In [451]:

    

# sphere object

class Sphere():
    
    def __init__(self, centre, r, colour):
        self.name = "sphere"
        self.centre = centre
        self.radius = r
        self.colour = colour
        pass
    
    def status(self):
        print(self.name)
        print("centre = ", self.centre)
        print("radius = ", self.radius)
        print("colour = ", self.colour)
        print("")
        pass

    def intersection(self, camera_location, ray_direction_vector):
        # calculate quadratic determinant "b^2 - 4ac" for ray intersecting circle
        b = numpy.dot(2 * ray_direction_vector,(camera_location - self.centre))
        b2 = b*b
        a = numpy.dot(ray_direction_vector, ray_direction_vector)
        c = numpy.dot((self.centre - camera_location), (self.centre - camera_location)) - (self.radius * self.radius)
        delta = b2 - (4 * a * c)
        #print(delta)
        
        if (delta >= 0):
            # calculate nearest point (lowest t)
            t = (-b - math.sqrt(delta)) / (2 * a)
            intersection_point = camera_location + (t * ray_direction_vector)
            #print("sphere")
            #print("camera_location, ray_direction_vector", camera_location, ray_direction_vector)
            #print("t =", t)
            #print("intersection_point ", intersection_point)
            
            # calculate normal at surface
            normal = (intersection_point - self.centre) / numpy.linalg.norm(intersection_point - self.centre)
            
            # additional check to ensure intersection was in direction of ray, t>0
            if (t<0):
                # set delta negative as if not intersected
                delta = -1
                pass
            
            # apply texture
            ## colour = self.colour + texture_wavy(intersection_point)
            colour = self.colour + texture_noise(intersection_point)
            
            # return tuple (intersection yes/no, nearest point, normal)
            return (delta>0, intersection_point, normal, colour)
            pass
        
        # return tuple (intersection yes/no, nearest point, norm, colour)
        return (delta >= 0, 0, 0, numpy.array([0,0,0]))


    

In [452]:

    

# flat plane object

class Plane():
    
    def __init__(self, norm, x, colour):
        self.name = "plane"
        self.norm = norm
        self.x = x
        self.colour = colour
        pass
    
    def status(self):
        print(self.name)
        print("norm = ", self.norm)
        print("point x =", self.x)
        print("colour = ", self.colour)
        print("")
        pass
    
    def intersection(self, camera_location, ray_direction_vector):
        intersects = False
        
        denominator = numpy.dot(ray_direction, self.norm)
        
        # check ray is not parallel to the plane,and is in fact facing ray (cos < 0)
        if (denominator < 0):
            # when does it intersect
            t = numpy.dot((self.x - camera_location), self.norm) / denominator
            
            # check direction of ray is forward, t>0 .. but with tolerance for floating point errors
            if (t > 0.001):
                intersects = True
                intersection_point = camera_location + (t * ray_direction_vector)
                # return tuple (intersection yes/no, nearest point, normal)
                #print("plane")
                #print("camera_location, ray_direction_vector", camera_location, ray_direction_vector)
                #print("t =", t)
                #print("intersection_point ", intersection_point)
                return (intersects, intersection_point, self.norm, self.colour)
                pass
            
            pass
        
        # return tuple (intersection yes/no, nearest point, norm, colour)
        return (intersects, 0, 0, numpy.array([0,0,0]))


    

In [453]:

    

# camera location
camera_location = numpy.array([0,100,-400])

# view port
view_port_location = numpy.array([-10, 0, 0])
view_port_width = 20
view_port_height = 20

# resolution (pixels per unit distance)
resolution = 40


    

In [454]:

    

# light source

# light is at above right, and a bit forward
light_location = numpy.array([10, 20, 0])

# ambient light factor
ambient = 0.25


    

In [455]:

    

# scene is a list of objects
scene_objects = []

# add spheres
for x in range(-8, 9, 4):
    for z in range(-16, 26, 8):
        scene_objects.append(Sphere(numpy.array([x,10,z]), 2, numpy.array([1.0, 0.0, 0.0])))
        pass
    pass

# add plane
#scene_objects.append(Plane(numpy.array([0.0,1.0,0.0]), numpy.array([0.0,6.0,0.0]), numpy.array([0.0,0.5,0.0])))

# get status of each object
for obj in scene_objects:
    obj.status()
    pass


    

    




sphere
centre =  [ -8  10 -16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-8 10 -8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-8 10  0]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-8 10  8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-8 10 16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-8 10 24]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ -4  10 -16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-4 10 -8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-4 10  0]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-4 10  8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-4 10 16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [-4 10 24]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [  0  10 -16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 0 10 -8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 0 10  0]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 0 10  8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 0 10 16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 0 10 24]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [  4  10 -16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 4 10 -8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 4 10  0]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 4 10  8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 4 10 16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 4 10 24]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [  8  10 -16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 8 10 -8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 8 10  0]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 8 10  8]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 8 10 16]
radius =  2
colour =  [ 1.  0.  0.]

sphere
centre =  [ 8 10 24]
radius =  2
colour =  [ 1.  0.  0.]

In [456]:

    

# create image

image = numpy.zeros([view_port_width * resolution, view_port_height * resolution, 3], dtype='float64')
print("image shape = ", image.shape)


    

    




image shape =  (800, 800, 3)

In [457]:

    

# max ray tracing depth for scene
max_depth = 3

# recursive ray function
def ray(ray_origin, ray_direction, depth):
    
    # start with zero colour contribution, false intersected flag and nearest point at infinity
    colour_contribution = numpy.array([0.0,0.0,0.0])    
    intersected = False
    distance_to_nearest_point = numpy.Infinity
    nearest_obj = False
    
    for obj in scene_objects:
        
        # check intersection with object
        intersection = obj.intersection(ray_origin, ray_direction)
        
        # did it intersect?
        if intersection[0]:
            intersected = True
            
            distance_to_point = numpy.linalg.norm(intersection[1] - ray_origin)
            
            if (distance_to_point < distance_to_nearest_point):
                
                # update nearest point
                distance_to_nearest_point = distance_to_point
                point_on_nearest_obj = intersection[1]
                norm = intersection[2]
                obj_colour = intersection[3]
                nearest_obj = obj
                
                # vector to light
                to_light = light_location - point_on_nearest_obj
                to_light /= numpy.linalg.norm(to_light)
                
                # illumination factor
                cos_norm_to_lightsource = numpy.dot(to_light, norm) / numpy.linalg.norm(to_light)
                
                # reflected factor (will have norm 1)
                reflected_ray = (2 * numpy.dot(norm, to_light) * norm) - to_light
                cos_ray_to_reflection = numpy.dot(-ray_direction, reflected_ray) / (numpy.linalg.norm(reflected_ray))
            
                # clip if below zero
                cos_norm_to_lightsource = numpy.clip(cos_norm_to_lightsource, 0, 1)
                cos_ray_to_reflection = numpy.clip(cos_ray_to_reflection, 0, 1)
                
                # apply power
                cos_norm_to_lightsource = math.pow(cos_norm_to_lightsource, 2)
                cos_ray_to_reflection = math.pow(cos_ray_to_reflection, 2)
                
                # candidate colour contribution
                candidate_colour_contribution = (obj_colour * cos_norm_to_lightsource * cos_ray_to_reflection)
                pass
            
            pass
        
        #next object
        pass
    
    # fire off new ray if intersected and depth is within max
    if (intersected):
        
        in_shadow = False
        for obj in scene_objects:
            # enusre object isn't current object from which shadow ray starts
            if (obj != nearest_obj):
                # check intersection with object
                if obj.intersection(point_on_nearest_obj, to_light)[0]:
                    in_shadow = True
                    pass
                pass
            pass
        
        # only add colour contribution if not in shadow
        if (not in_shadow):
            # colour contribution from nearest object
            colour_contribution = candidate_colour_contribution
            pass
        
        # light falloff factor
        distance_factor = numpy.linalg.norm(light_location - point_on_nearest_obj) - 15.0
        #print("d", distance_factor)
        falloff = (1.0 - math.tanh(distance_factor)) / 2.0
        #print("f", falloff)
        #print("")
        
        
        # add ambient light, diminish with depth
        colour_contribution += obj_colour * ambient / (depth) * falloff
        
        # recurse colour contribution from bounced ray
        if (depth < max_depth):
            # bounced ray (from source, not from light source)
            bounced_ray = (2 * numpy.dot(norm, -ray_direction) * norm) + ray_direction
            # divide by depth to reduce strength of light after many reflections
            colour_contribution += (ray(point_on_nearest_obj, bounced_ray, depth + 1)[1] / ( depth))
            pass
        
        pass
    
    return intersected, colour_contribution


    

In [458]:

    

# main loop is to consider every pixel of the viewport

for pixel_ix in range(image.shape[0]):
    print(pixel_ix, image.shape[0])
    for pixel_iy in range(image.shape[1]):
#for pixel_ix in range(20,21):
#    for pixel_iy in range(15,16):
        
        # ray direction
        current_position = view_port_location + numpy.array([pixel_ix/resolution, pixel_iy/resolution, 0])
        ray_direction = current_position - camera_location
        ray_direction /= numpy.linalg.norm(ray_direction)
        
        # pixel is set to colour contribution from (recursive) ray
        intersected, colour_contribution = ray(camera_location, ray_direction, 1)
        if (intersected):
            image[pixel_ix, pixel_iy] = colour_contribution
        else:
            image[pixel_ix, pixel_iy] = [0.3 + ray_direction[1], 0.3 + ray_direction[1], 0.6 + (ray_direction[1] * 5)]
            pass
        
        pass
    pass


    

    




0 800
1 800
2 800
3 800
4 800
5 800
6 800
7 800
8 800
9 800
10 800
11 800
12 800
13 800
14 800
15 800
16 800
17 800
18 800
19 800
20 800
21 800
22 800
23 800
24 800
25 800
26 800
27 800
28 800
29 800
30 800
31 800
32 800
33 800
34 800
35 800
36 800
37 800
38 800
39 800
40 800
41 800
42 800
43 800
44 800
45 800
46 800
47 800
48 800
49 800
50 800
51 800
52 800
53 800
54 800
55 800
56 800
57 800
58 800
59 800
60 800
61 800
62 800
63 800
64 800
65 800
66 800
67 800
68 800
69 800
70 800
71 800
72 800
73 800
74 800
75 800
76 800
77 800
78 800
79 800
80 800
81 800
82 800
83 800
84 800
85 800
86 800
87 800
88 800
89 800
90 800
91 800
92 800
93 800
94 800
95 800
96 800
97 800
98 800
99 800
100 800
101 800
102 800
103 800
104 800
105 800
106 800
107 800
108 800
109 800
110 800
111 800
112 800
113 800
114 800
115 800
116 800
117 800
118 800
119 800
120 800
121 800
122 800
123 800
124 800
125 800
126 800
127 800
128 800
129 800
130 800
131 800
132 800
133 800
134 800
135 800
136 800
137 800
138 800
139 800
140 800
141 800
142 800
143 800
144 800
145 800
146 800
147 800
148 800
149 800
150 800
151 800
152 800
153 800
154 800
155 800
156 800
157 800
158 800
159 800
160 800
161 800
162 800
163 800
164 800
165 800
166 800
167 800
168 800
169 800
170 800
171 800
172 800
173 800
174 800
175 800
176 800
177 800
178 800
179 800
180 800
181 800
182 800
183 800
184 800
185 800
186 800
187 800
188 800
189 800
190 800
191 800
192 800
193 800
194 800
195 800
196 800
197 800
198 800
199 800
200 800
201 800
202 800
203 800
204 800
205 800
206 800
207 800
208 800
209 800
210 800
211 800
212 800
213 800
214 800
215 800
216 800
217 800
218 800
219 800
220 800
221 800
222 800
223 800
224 800
225 800
226 800
227 800
228 800
229 800
230 800
231 800
232 800
233 800
234 800
235 800
236 800
237 800
238 800
239 800
240 800
241 800
242 800
243 800
244 800
245 800
246 800
247 800
248 800
249 800
250 800
251 800
252 800
253 800
254 800
255 800
256 800
257 800
258 800
259 800
260 800
261 800
262 800
263 800
264 800
265 800
266 800
267 800
268 800
269 800
270 800
271 800
272 800
273 800
274 800
275 800
276 800
277 800
278 800
279 800
280 800
281 800
282 800
283 800
284 800
285 800
286 800
287 800
288 800
289 800
290 800
291 800
292 800
293 800
294 800
295 800
296 800
297 800
298 800
299 800
300 800
301 800
302 800
303 800
304 800
305 800
306 800
307 800
308 800
309 800
310 800
311 800
312 800
313 800
314 800
315 800
316 800
317 800
318 800
319 800
320 800
321 800
322 800
323 800
324 800
325 800
326 800
327 800
328 800
329 800
330 800
331 800
332 800
333 800
334 800
335 800
336 800
337 800
338 800
339 800
340 800
341 800
342 800
343 800
344 800
345 800
346 800
347 800
348 800
349 800
350 800
351 800
352 800
353 800
354 800
355 800
356 800
357 800
358 800
359 800
360 800
361 800
362 800
363 800
364 800
365 800
366 800
367 800
368 800
369 800
370 800
371 800
372 800
373 800
374 800
375 800
376 800
377 800
378 800
379 800
380 800
381 800
382 800
383 800
384 800
385 800
386 800
387 800
388 800
389 800
390 800
391 800
392 800
393 800
394 800
395 800
396 800
397 800
398 800
399 800
400 800
401 800
402 800
403 800
404 800
405 800
406 800
407 800
408 800
409 800
410 800
411 800
412 800
413 800
414 800
415 800
416 800
417 800
418 800
419 800
420 800
421 800
422 800
423 800
424 800
425 800
426 800
427 800
428 800
429 800
430 800
431 800
432 800
433 800
434 800
435 800
436 800
437 800
438 800
439 800
440 800
441 800
442 800
443 800
444 800
445 800
446 800
447 800
448 800
449 800
450 800
451 800
452 800
453 800
454 800
455 800
456 800
457 800
458 800
459 800
460 800
461 800
462 800
463 800
464 800
465 800
466 800
467 800
468 800
469 800
470 800
471 800
472 800
473 800
474 800
475 800
476 800
477 800
478 800
479 800
480 800
481 800
482 800
483 800
484 800
485 800
486 800
487 800
488 800
489 800
490 800
491 800
492 800
493 800
494 800
495 800
496 800
497 800
498 800
499 800
500 800
501 800
502 800
503 800
504 800
505 800
506 800
507 800
508 800
509 800
510 800
511 800
512 800
513 800
514 800
515 800
516 800
517 800
518 800
519 800
520 800
521 800
522 800
523 800
524 800
525 800
526 800
527 800
528 800
529 800
530 800
531 800
532 800
533 800
534 800
535 800
536 800
537 800
538 800
539 800
540 800
541 800
542 800
543 800
544 800
545 800
546 800
547 800
548 800
549 800
550 800
551 800
552 800
553 800
554 800
555 800
556 800
557 800
558 800
559 800
560 800
561 800
562 800
563 800
564 800
565 800
566 800
567 800
568 800
569 800
570 800
571 800
572 800
573 800
574 800
575 800
576 800
577 800
578 800
579 800
580 800
581 800
582 800
583 800
584 800
585 800
586 800
587 800
588 800
589 800
590 800
591 800
592 800
593 800
594 800
595 800
596 800
597 800
598 800
599 800
600 800
601 800
602 800
603 800
604 800
605 800
606 800
607 800
608 800
609 800
610 800
611 800
612 800
613 800
614 800
615 800
616 800
617 800
618 800
619 800
620 800
621 800
622 800
623 800
624 800
625 800
626 800
627 800
628 800
629 800
630 800
631 800
632 800
633 800
634 800
635 800
636 800
637 800
638 800
639 800
640 800
641 800
642 800
643 800
644 800
645 800
646 800
647 800
648 800
649 800
650 800
651 800
652 800
653 800
654 800
655 800
656 800
657 800
658 800
659 800
660 800
661 800
662 800
663 800
664 800
665 800
666 800
667 800
668 800
669 800
670 800
671 800
672 800
673 800
674 800
675 800
676 800
677 800
678 800
679 800
680 800
681 800
682 800
683 800
684 800
685 800
686 800
687 800
688 800
689 800
690 800
691 800
692 800
693 800
694 800
695 800
696 800
697 800
698 800
699 800
700 800
701 800
702 800
703 800
704 800
705 800
706 800
707 800
708 800
709 800
710 800
711 800
712 800
713 800
714 800
715 800
716 800
717 800
718 800
719 800
720 800
721 800
722 800
723 800
724 800
725 800
726 800
727 800
728 800
729 800
730 800
731 800
732 800
733 800
734 800
735 800
736 800
737 800
738 800
739 800
740 800
741 800
742 800
743 800
744 800
745 800
746 800
747 800
748 800
749 800
750 800
751 800
752 800
753 800
754 800
755 800
756 800
757 800
758 800
759 800
760 800
761 800
762 800
763 800
764 800
765 800
766 800
767 800
768 800
769 800
770 800
771 800
772 800
773 800
774 800
775 800
776 800
777 800
778 800
779 800
780 800
781 800
782 800
783 800
784 800
785 800
786 800
787 800
788 800
789 800
790 800
791 800
792 800
793 800
794 800
795 800
796 800
797 800
798 800
799 800

In [459]:

    

# apply squashing function to image
# first shift data into range [0,1] asymptotically
# then remap to colour RGB range [0,255] dtype=uint8

# squash with tanh()
image = numpy.tanh(image)

# remap to RGB range
image_rgb = numpy.array(image*255, dtype='uint8')


    

In [460]:

    

# transpose array so origin is bottom left, by swapping dimensions 0 and 1, but leave dimension 3

image_rgb2 = numpy.transpose(image_rgb, (1, 0, 2))
plt.imshow(image_rgb2, origin='lower')


    

    
Out[460]:





<matplotlib.image.AxesImage at 0x119cb3ba8>






    

In [461]:

    

plt.imsave('test.png', image_rgb2, origin='lower')


    

In [ ]: