Pythonによるレイトレーシングを高速化してみる

例

プリミティブは球体のみ、シェーディングはフラットシェーディング、カメラは正投影、解像度固定。

• 最初のバージョン(ray1.py)

  import math
  from itertools import product
  from functools import partial
  
  
  def range2(rx, ry):
      return product(range(rx), range(ry))
  
  
  class Vector(object):
      def __init__(self, x, y, z):
          self.x = x
          self.y = y
          self.z = z
  
      def normalize(self):
          v = math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
          self.x /= v
          self.y /= v
          self.z /= v
  
      def __mul__(self, v):
          return self.x*v.x + self.y*v.y + self.z*v.z
  
      def __sub__(self, v):
          return Vector(self.x - v.x, self.y - v.y, self.z - v.z)
  
      def printval(self):
          print self.x, self.y, self.z
  
  
  class Color(Vector):
      def __init__(self, r, g, b):
          super(Color, self).__init__(r, g, b)
  
      def printval(self):
          print '%d %d %d' % (self.x*255, self.y*255, self.z*255)
  
  
  class Point(Vector):
      def __init__(self, x, y, z):
          super(Point, self).__init__(x, y, z)
  
  
  class Ray(object):
      def __init__(self, p, v):
          self.p = p
          self.v = v
          self.v.normalize()
  
  
  class Primitive(object):
      def __init__(self):
          pass
  
      def intersect(self, ray):
          return (None, None)
  
      def color(self):
          pass
  
  
  class Back(Primitive):
      def __init__(self):
          super(Back, self).__init__()
  
      def intersect(self, ray):
          return (self, float('inf'))
  
      def color(self):
          return Color(0.5, 0.5, 0.5)
  
  
  class Sphere(Primitive):
      def __init__(self, p, rad, col):
          self.cp = p
          self.rad = rad
          self.col = col
  
      def _pmin(self, a, b):
          t = [0, a, b]
          t.sort()
          t.append(None)
          return t[t.index(0) + 1]
  
      def intersect(self, ray):
          v = ray.p - self.cp
          n = v*ray.v
          det = n*n - v*v + self.rad*self.rad
          if det < 0.0:
              return (self, None)
          d = math.sqrt(det)
          t = self._pmin(-n - d, -n + d)
          return (self, t)
  
      def color(self):
          return self.col
  
  
  def create_ray(rx, ry):
      f = lambda r, x: 1.0 - 2.0/r*x
      v = Vector(0.0, 0.0, -1.0)
      rs = [Ray(Point(f(rx, x), f(ry, y), 2.0), v) for x, y in range2(rx, ry)]
      return rs
  
  
  def get_objs():
      return [Sphere(Point(0.0, 0.0, 0.0), 0.2, Color(1.0, 0.0, 0.0)),
              Sphere(Point(0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Back()]
  
  
  def trace(ray):
      objs = get_objs()
      func = lambda ray, obj: obj.intersect(ray)
      ts = map(partial(func, ray), objs)
      ts = filter(lambda x: x[1], ts)
      ts.sort(key=lambda x: x[1])
      return ts[0][0].color()
  
  
  def output(rx, ry, cs):
      print 'P3\n' + str(rx) + ' ' + str(ry) + '\n255'
      map(lambda c: c.printval(), cs)
  
  
  def main():
      rs = create_ray(400, 400)
      cs = map(trace, rs)
      output(400, 400, cs)
  
  main()
  

• 交差判定部分をC言語にしてみたバージョン(ray2.py)

  import math
  import myray
  from itertools import product
  from functools import partial
  
  
  def range2(rx, ry):
      return product(range(rx), range(ry))
  
  
  class Vector(object):
      def __init__(self, x, y, z):
          self.x = x
          self.y = y
          self.z = z
  
      def normalize(self):
          v = math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
          self.x /= v
          self.y /= v
          self.z /= v
  
      def __mul__(self, v):
          return self.x*v.x + self.y*v.y + self.z*v.z
  
      def __sub__(self, v):
          return Vector(self.x - v.x, self.y - v.y, self.z - v.z)
  
      def printval(self):
          print self.x, self.y, self.z
  
  
  class Color(Vector):
      def __init__(self, r, g, b):
          super(Color, self).__init__(r, g, b)
  
      def printval(self):
          print '%d %d %d' % (self.x*255, self.y*255, self.z*255)
  
  
  class Point(Vector):
      def __init__(self, x, y, z):
          super(Point, self).__init__(x, y, z)
  
  
  class Ray(object):
      def __init__(self, p, v):
          self.p = p
          self.v = v
          self.v.normalize()
  
  
  class Primitive(object):
      def __init__(self):
          pass
  
      def intersect(self, ray):
          return (None, None)
  
      def color(self):
          pass
  
  
  class Back(Primitive):
      def __init__(self):
          super(Back, self).__init__()
  
      def intersect(self, ray):
          return (self, float('inf'))
  
      def color(self):
          return Color(0.5, 0.5, 0.5)
  
  
  class Sphere(Primitive):
      def __init__(self, p, rad, col):
          self.cp = p
          self.rad = rad
          self.col = col
  
      def intersect(self, ray):
          t = myray.intersect(ray.p.x, ray.p.y, ray.p.z,
                              ray.v.x, ray.v.y, ray.v.z,
                              self.cp.x, self.cp.y, self.cp.z,
                              self.rad)
          return (self, t)
  
      def color(self):
          return self.col
  
  
  def create_ray(rx, ry):
      f = lambda r, x: 1.0 - 2.0/r*x
      v = Vector(0.0, 0.0, -1.0)
      rs = [Ray(Point(f(rx, x), f(ry, y), 2.0), v) for x, y in range2(rx, ry)]
      return rs
  
  
  def get_objs():
      return [Sphere(Point(0.0, 0.0, 0.0), 0.2, Color(1.0, 0.0, 0.0)),
              Sphere(Point(0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Back()]
  
  
  def trace(ray):
      objs = get_objs()
      func = lambda ray, obj: obj.intersect(ray)
      ts = map(partial(func, ray), objs)
      ts = filter(lambda x: x[1], ts)
      ts.sort(key=lambda x: x[1])
      return ts[0][0].color()
  
  
  def output(rx, ry, cs):
      print 'P3\n' + str(rx) + ' ' + str(ry) + '\n255'
      map(lambda c: c.printval(), cs)
  
  
  def main():
      rs = create_ray(400, 400)
      cs = map(trace, rs)
      output(400, 400, cs)
  
  main()
  

  C言語の部分(myray.c)

  #include <Python.h>
  #include <math.h>
  
  typedef struct {
      double x, y, z;
      } VECTOR;
  
  VECTOR vminus(VECTOR v1, VECTOR v2)
  {
      VECTOR v;
  
      v.x = v1.x - v2.x;
      v.y = v1.y - v2.y;
      v.z = v1.z - v2.z;
      return v;
  }
  
  double innerp(VECTOR v1, VECTOR v2)
  {
      return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
  }
  
  static PyObject *
  intersect(PyObject *self, PyObject *args)
  {
  	double rad;
  	double n;
  	double d, d1, d2;
  	double det;
  	VECTOR ray, rp, sp, v;
  
  	if(!PyArg_ParseTuple(args, "dddddddddd", &rp.x, &rp.y, &rp.z, &ray.x, &ray.y, &ray.z, &sp.x, &sp.y, &sp.z, &rad))
  		return NULL;
  	v = vminus(rp, sp);
  	n = innerp(v, ray);
  	det = n*n - innerp(v, v) + rad*rad;
  	if(det < 0.0)
  		return Py_BuildValue("");
  	d = sqrt(det);
  	d1 = -n - d;
  	d2 = -n + d;
  	if(d2 < 0.0)
  		return Py_BuildValue("");
  	if(d1 < 0.0)
  		return Py_BuildValue("d", d2);
  	return Py_BuildValue("d", d1);
  }
  
  static PyMethodDef myray_methods[] = {
  	{"intersect", intersect, METH_VARARGS, "sphere intersect."},
  	{NULL, NULL, 0, NULL}
  	};
  
  void initmyray()
  {
  	Py_InitModule("myray", myray_methods);
  }
  
  

  (setup.py)

  from distutils.core import setup, Extension
  
  module = Extension('myray', ['myray.c'])
  setup(name='myray', version='1.0', ext_modules=[module])
  

• データのインスタンス作成を少なくしたバージョン(ray3.py)

  import myray
  import math
  from itertools import product
  from functools import partial
  
  
  def range2(rx, ry):
      return product(range(rx), range(ry))
  
  
  class Vector(object):
      def __init__(self, x, y, z):
          self.x = x
          self.y = y
          self.z = z
  
      def normalize(self):
          v = math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
          self.x /= v
          self.y /= v
          self.z /= v
  
      def __mul__(self, v):
          return self.x*v.x + self.y*v.y + self.z*v.z
  
      def __sub__(self, v):
          return Vector(self.x - v.x, self.y - v.y, self.z - v.z)
  
      def printval(self):
          print self.x, self.y, self.z
  
  
  class Color(Vector):
      def __init__(self, r, g, b):
          super(Color, self).__init__(r, g, b)
  
      def printval(self):
          print '%d %d %d' % (self.x*255, self.y*255, self.z*255)
  
  
  class Point(Vector):
      def __init__(self, x, y, z):
          super(Point, self).__init__(x, y, z)
  
  
  class Ray(object):
      def __init__(self, p, v):
          self.p = p
          self.v = v
          self.v.normalize()
  
  
  class Primitive(object):
      def __init__(self):
          pass
  
      def intersect(self, ray):
          return (None, None)
  
      def color(self):
          pass
  
  
  class Back(Primitive):
      def __init__(self):
          super(Back, self).__init__()
  
      def intersect(self, ray):
          return (self, float('inf'))
  
      def color(self):
          return Color(0.5, 0.5, 0.5)
  
  
  class Sphere(Primitive):
      def __init__(self, p, rad, col):
          self.cp = p
          self.rad = rad
          self.col = col
  
      def intersect(self, ray):
          t = myray.intersect(ray.p.x, ray.p.y, ray.p.z,
                              ray.v.x, ray.v.y, ray.v.z,
                              self.cp.x, self.cp.y, self.cp.z,
                              self.rad)
          return (self, t)
  
      def color(self):
          return self.col
  
  
  class Objects(object):
      data = [Sphere(Point(0.0, 0.0, 0.0), 0.2, Color(1.0, 0.0, 0.0)),
              Sphere(Point(0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, 0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(-0.2, -0.2, -0.5), 0.2, Color(0.0, 1.0, 0.0)),
              Sphere(Point(0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.5, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.2, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.2, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, 0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, 0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(0.5, -0.4, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Sphere(Point(-0.4, -0.5, -1.0), 0.2, Color(0.0, 0.0, 1.0)),
              Back()]
  
      def __init__(self):
          pass
  
      def get(self):
          return self.data
  
  
  def create_ray(rx, ry):
      f = lambda r, x: 1.0 - 2.0/r*x
      v = Vector(0.0, 0.0, -1.0)
      rs = [Ray(Point(f(rx, x), f(ry, y), 2.0), v) for x, y in range2(rx, ry)]
      return rs
  
  
  def trace(ray):
      objs = Objects().get()
      func = lambda ray, obj: obj.intersect(ray)
      ts = map(partial(func, ray), objs)
      ts = filter(lambda x: x[1], ts)
      ts.sort(key=lambda x: x[1])
      return ts[0][0].color()
  
  
  def output(rx, ry, cs):
      print 'P3\n' + str(rx) + ' ' + str(ry) + '\n255'
      map(lambda c: c.printval(), cs)
  
  
  def main():
      rs = create_ray(400, 400)
      cs = map(trace, rs)
      output(400, 400, cs)
  
  main()
  

  C言語部分はray2.pyと同じ。

時間測定

以下のコマンドで計測した結果
time python ray*.py > tmp*.ppm

ray1.py

56.656u 0.493s 0:57.25 99.8% 0+0k 0+8io 0pf+0w

ray2.py

42.021u 0.459s 0:42.50 99.9% 0+0k 0+5io 0pf+0w

ray3.py

21.830u 0.384s 0:22.23 99.9% 0+0k 0+3io 0pf+0w