title: Improving performance of raytracing hit function
body: <p>I have a simple raytracer in python. rendering an image 200x200 takes 4 minutes, which is definitely too much for my taste. I want to improve the situation.</p>
  
  <p>Some points: I shoot multiple rays per each pixel (to provide antialiasing) for a grand total of 16 rays per pixel. 200x200x16 is a grand total of 640000 rays. Each ray must be tested for impact on multiple Sphere objects in the scene. Ray is also a rather trivial object</p>
  
  <pre><code>class Ray(object):
      def __init__(self, origin, direction):
          self.origin = numpy.array(origin)
          self.direction = numpy.array(direction)
  </code></pre>
  
  <p>Sphere is slightly more complex, and carries the logic for hit/nohit:</p>
  
  <pre><code>class Sphere(object):
      def __init__(self, center, radius, color):
          self.center = numpy.array(center)
          self.radius = numpy.array(radius)
          self.color = color
  
      @profile 
      def hit(self, ray):
          temp = ray.origin - self.center
          a = numpy.dot(ray.direction, ray.direction)
          b = 2.0 * numpy.dot(temp, ray.direction)
          c = numpy.dot(temp, temp) - self.radius * self.radius
          disc = b * b - 4.0 * a * c
  
          if (disc &lt; 0.0):
              return None
          else:
              e = math.sqrt(disc)
              denom = 2.0 * a
              t = (-b - e) / denom 
              if (t &gt; 1.0e-7):
                  normal = (temp + t * ray.direction) / self.radius
                  hit_point = ray.origin + t * ray.direction
                  return ShadeRecord.ShadeRecord(normal=normal, 
                                                 hit_point=hit_point, 
                                                 parameter=t, 
                                                 color=self.color)           
  
              t = (-b + e) / denom
  
              if (t &gt; 1.0e-7):
                  normal = (temp + t * ray.direction) / self.radius                hit_point = ray.origin + t * ray.direction
                  return ShadeRecord.ShadeRecord(normal=normal, 
                                                 hit_point=hit_point, 
                                                 parameter=t, 
                                                 color=self.color)       
  
          return None    
  </code></pre>
  
  <p>Now, I ran some profiling, and it appears that the longest processing time is in the hit() function</p>
  
  <pre><code>   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    2560000  118.831    0.000  152.701    0.000 raytrace/objects/Sphere.py:12(hit)
    1960020   42.989    0.000   42.989    0.000 {numpy.core.multiarray.array}
          1   34.566   34.566  285.829  285.829 raytrace/World.py:25(render)
    7680000   33.796    0.000   33.796    0.000 {numpy.core._dotblas.dot}
    2560000   11.124    0.000  163.825    0.000 raytrace/World.py:63(f)
     640000   10.132    0.000  189.411    0.000 raytrace/World.py:62(hit_bare_bones_object)
     640023    6.556    0.000  170.388    0.000 {map}
  </code></pre>
  
  <p>This does not surprise me, and I want to reduce this value as much as possible. I pass to line profiling, and the result is </p>
  
  <pre><code>Line #      Hits         Time  Per Hit   % Time  Line Contents
  ==============================================================
      12                                               @profile
      13                                               def hit(self, ray):
      14   2560000     27956358     10.9     19.2          temp = ray.origin - self.center
      15   2560000     17944912      7.0     12.3          a = numpy.dot(ray.direction, ray.direction)
      16   2560000     24132737      9.4     16.5          b = 2.0 * numpy.dot(temp, ray.direction)
      17   2560000     37113811     14.5     25.4          c = numpy.dot(temp, temp) - self.radius * self.radius
      18   2560000     20808930      8.1     14.3          disc = b * b - 4.0 * a * c
      19                                                   
      20   2560000     10963318      4.3      7.5          if (disc &lt; 0.0):
      21   2539908      5403624      2.1      3.7              return None
      22                                                   else:
      23     20092        75076      3.7      0.1              e = math.sqrt(disc)
      24     20092       104950      5.2      0.1              denom = 2.0 * a
      25     20092       115956      5.8      0.1              t = (-b - e) / denom
      26     20092        83382      4.2      0.1              if (t &gt; 1.0e-7):
      27     20092       525272     26.1      0.4                  normal = (temp + t * ray.direction) / self.radius
      28     20092       333879     16.6      0.2                  hit_point = ray.origin + t * ray.direction
      29     20092       299494     14.9      0.2                  return ShadeRecord.ShadeRecord(normal=normal, hit_point=hit_point, parameter=t, color=self.color)
  </code></pre>
  
  <p>So, it appears that most of the time is spent in this chunk of code:</p>
  
  <pre><code>        temp = ray.origin - self.center
          a = numpy.dot(ray.direction, ray.direction)
          b = 2.0 * numpy.dot(temp, ray.direction)
          c = numpy.dot(temp, temp) - self.radius * self.radius
          disc = b * b - 4.0 * a * c
  </code></pre>
  
  <p>Where I don't really see a lot to optimize. Do you have any idea how to make this code more performant without going C ?</p>