diff -urN rt_v12/ax.py rt_v13/ax.py --- rt_v12/ax.py 2018-03-21 21:15:26.000000000 +0900 +++ rt_v13/ax.py 2018-03-30 16:27:00.000000000 +0900 @@ -1,6 +1,7 @@ #!/usr/bin/env python import ut +import mt import v import line import vecs @@ -9,13 +10,15 @@ def new(p, vs): # vs: vecs + p = mt.to_arr(p) + e = ut.Empty() e.typ = 'ax' e.p = p e.vs = vs - l2g = lambda p: v.op2( '+', vs.tr('l2g', p), e.p ) - g2l = lambda p: vs.tr( 'g2l', v.op2('-', p, e.p) ) + l2g = lambda p: vs.tr( 'l2g', p ) + e.p + g2l = lambda p: vs.tr( 'g2l', p - e.p ) e.tr_p = lambda d, p: { 'l2g': l2g, 'g2l': g2l }.get(d, l2g)(p) e.tr_v = lambda d, v: vs.tr(d, v) e.tr_line = lambda d, l: line.new( e.tr_p(d, l.p), e.tr_v(d, l.v) ) @@ -46,7 +49,7 @@ vs_ = vs.rev() if not vs_: return [] - p_ = vs_.tr( 'l2g', v.nega(p) ) + p_ = vs_.tr( 'l2g', -p ) return new(p_, vs_) e.rev = rev diff -urN rt_v12/cross.py rt_v13/cross.py --- rt_v12/cross.py 2018-03-26 18:27:58.000000000 +0900 +++ rt_v13/cross.py 2018-04-01 10:49:17.000000000 +0900 @@ -1,6 +1,7 @@ #!/usr/bin/env python import math +import numpy as np import mt import v import line @@ -12,12 +13,11 @@ def plane_z0(l, d): # z=0 (x-y) - (vx, vy, vz) = l.v - (px, py, pz) = l.p - + vz = l.v[2] if vz == 0: return {} + pz = l.p[2] t = -float(pz) / vz if t <= 0: return {} @@ -65,9 +65,11 @@ # r = 1, c = zero # l = p + v t # v^2 t^2 + 2 v p t + p^2 = 1^2 - a = sum( v.op2('*', l.v, l.v) ) - b = 2 * sum( v.op2('*', l.v, l.p) ) - c = sum( v.op2('*', l.p, l.p) ) - 1 + + a = np.dot( l.v, l.v ) + b = 2 * np.dot( l.v, l.p ) + c = np.dot( l.p, l.p ) - 1 + ts = mt.quadratic_formula(a, b, c) ts = sorted( filter( lambda v: v > 0, ts ) ) if not ts: diff -urN rt_v12/cylx.py rt_v13/cylx.py --- rt_v12/cylx.py 2018-03-16 20:38:43.000000000 +0900 +++ rt_v13/cylx.py 2018-03-30 14:51:29.000000000 +0900 @@ -2,6 +2,7 @@ import math import ut +import mt import v import line @@ -14,13 +15,13 @@ (x, y, z) = p rad = x/e.r r = z - return [ r*math.cos(r), r*math.sin(r), y ] + return mt.arr( [ r*math.cos(r), r*math.sin(r), y ] ) def g2l(p): (x, y, z) = p rad = math.atan2(y, x) r = v.len1( [x,y,0] ) - return [ e.r*rad, z, r ] + return mt.arr( [ e.r*rad, z, r ] ) dval = lambda d: d if not rev else { 'l2g': 'g2l', 'g2l': 'l2g' }.get(d, 'g2l') e.tr_p = lambda d, p: { 'l2g': l2g, 'g2l': g2l }.get( dval(d), l2g)(p) diff -urN rt_v12/dat.py rt_v13/dat.py --- rt_v12/dat.py 2018-03-28 22:31:01.000000000 +0900 +++ rt_v13/dat.py 2018-03-30 16:35:54.000000000 +0900 @@ -1,6 +1,7 @@ #!/usr/bin/env python import ut +import mt import v import line import ax @@ -161,8 +162,8 @@ data = data_dic.get( ut.arg_s('data_name'), [] ) lights = [ - { 'p': v.op1( '*', [1,-1,-1], 30 ), 'e': 2.5 }, - { 'p': v.op1( '*', [-1,1,1], 30 ), 'e': 2.5 }, + { 'p': mt.arr( [1,-1,-1] ) * 30, 'e': 2.5 }, + { 'p': mt.arr( [-1,1,1] ) * 30, 'e': 2.5 }, ] # EOF diff -urN rt_v12/fcx.py rt_v13/fcx.py --- rt_v12/fcx.py 2018-03-16 20:38:43.000000000 +0900 +++ rt_v13/fcx.py 2018-03-30 14:51:18.000000000 +0900 @@ -1,6 +1,7 @@ #!/usr/bin/env python import ut +import mt import v import line @@ -13,11 +14,11 @@ def l2g(p): (x, y, z) = p - return [ x*z/dz, y*z/dz, z ] if dz != 0 else [] + return mt.arr( [ x*z/dz, y*z/dz, z ] if dz != 0 else [] ) def g2l(p): (x, y, z) = p - return [ x*dz/z, y*dz/z, z ] if z != 0 else [] + return mt.arr( [ x*dz/z, y*dz/z, z ] if z != 0 else [] ) dval = lambda d: d if not rev else { 'l2g': 'g2l', 'g2l': 'l2g' }.get(d, 'g2l') e.tr_p = lambda d, p: { 'l2g': l2g, 'g2l': g2l }.get(d, l2g)(p) diff -urN rt_v12/line.py rt_v13/line.py --- rt_v12/line.py 2018-03-16 20:38:43.000000000 +0900 +++ rt_v13/line.py 2018-04-01 11:47:50.000000000 +0900 @@ -1,14 +1,18 @@ #!/usr/bin/env python import ut +import mt import v def new(p, v_): + p = mt.to_arr(p) + v_ = mt.to_arr(v_) + e = ut.Empty() e.typ = 'line' e.p = p e.v = v_ - e.on_line_p = lambda t: v.op2( '+', p, v.op1('*', e.v, t) ) + e.on_line_p = lambda t: p + e.v * t e.get_p_t = lambda : ( p, e.on_line_p(1) ) def tr_by_func(f): @@ -22,21 +26,21 @@ # (e.p, e.v) is plane # - e.is_plane_side = lambda p: v.dot_product( e.v, v.sub(p, e.p) ) >= 0 + e.is_plane_side = lambda p: np.dot( e.v, p - e.p ) >= 0 def cross_plane_line(other_line): l = other_line - base = v.dot_product(l.v, e.v) + base = np.dot(l.v, e.v) if base == 0: return [] - t = float( v.dot_product( v.sub(e.p, l.p), e.v ) ) / base + t = float( np.dot( e.p - l.p, e.v ) ) / base return l.on_line_p(t) e.cross_plane_line = cross_plane_line return e -new_p2_no_unit = lambda p, t: new( p, v.sub(t, p) ) -new_p2 = lambda p, t: new( p, v.sub_unit(t, p) ) +new_p2_no_unit = lambda p, t: new( p, t - p ) +new_p2 = lambda p, t: new( p, v.unit( t - p ) ) x1 = new( [0,0,0], v.x1 ) y1 = new( [0,0,0], v.y1 ) diff -urN rt_v12/lstx.py rt_v13/lstx.py --- rt_v12/lstx.py 2018-03-28 16:34:54.000000000 +0900 +++ rt_v13/lstx.py 2018-03-31 17:25:01.000000000 +0900 @@ -18,6 +18,8 @@ dval = lambda d, rev: d if not rev else { 'l2g': 'g2l', 'g2l': 'l2g' }.get(d, 'g2l') def tr_p(lst, d, p): + if len(lst) == 1: + return lst[0].tr_p(d, p) if d == 'g2l': lst = reversed(lst) return reduce( lambda p, x: x.tr_p(d, p), lst, p ) diff -urN rt_v12/mt.py rt_v13/mt.py --- rt_v12/mt.py 2018-03-19 14:49:47.000000000 +0900 +++ rt_v13/mt.py 2018-03-31 16:53:31.000000000 +0900 @@ -1,35 +1,21 @@ #!/usr/bin/env python import math +import numpy as np + +arr = lambda lst: np.array( lst, dtype=np.float32 ) +to_arr = lambda lst: arr(lst) if type(lst) != np.ndarray else lst +T = lambda v: v.reshape(v.size, 1) linear_conv = lambda s, sa, sb, da, db: da + (s - sa) * (db - da) / float(sb - sa) if sa != sb else (da + db) * 0.5 deg_to_rad = lambda deg: deg * math.pi / 180 rad_to_deg = lambda rad: rad * 180 / math.pi -def v3_det(v3): - (x, y, z) = v3 - (xx, xy, xz) = x - (yx, yy, yz) = y - (zx, zy, zz) = z - return xx*yy*zz + xy*yz*zx + xz*yx*zy - xx*yz*zy - xz*yy*zx - xy*yx*zz - def v3_rev(v3): - d = v3_det(v3) - if d == 0: + if np.linalg.det(v3) == 0: return [] - d_ = 1.0 / d - - (x, y, z) = v3 - (xx, xy, xz) = x - (yx, yy, yz) = y - (zx, zy, zz) = z - - tx = [ yy*zz - zy*yz, zy*xz - xy*zz, xy*yz - yy*xz ] - ty = [ zx*yz - yx*zz, xx*zz - zx*xz, yx*xz - xx*yz ] - tz = [ yx*zy - zx*yy, zx*xy - xx*zy, xx*yy - yx*xy ] - - return map( lambda v: map( lambda e: e * d_, v ), [tx, ty, tz] ) + return np.linalg.inv(v3) ### diff -urN rt_v12/rt.py rt_v13/rt.py --- rt_v12/rt.py 2018-03-23 11:23:36.000000000 +0900 +++ rt_v13/rt.py 2018-04-01 11:51:37.000000000 +0900 @@ -1,7 +1,10 @@ #!/usr/bin/env python +import sys import math +import numpy as np import ut +import mt import v import line import lstx @@ -18,7 +21,7 @@ xyz2g = m.get('xyz2g') nv_line = lstx.tr( xyz2g, 'g2l', crs.get('l') ) img_nv = v.z1 - if v.dot_product( nv_line.v, img_nv ) < 0: + if np.dot( nv_line.v, img_nv ) < 0: return [] fn = m.get( 'fn_r' if rev else 'fn', '' ) if not fn: @@ -29,7 +32,7 @@ cols = map( f, d.get('maps') ) cols = filter( lambda col: col != [], cols ) - col = v.lst_add(cols) if cols else def_col + col = np.sum(cols, axis=0) if cols else def_col return map( lambda v: min(v, 255), col ) def get_col(data, lights, l_g, video, nest_rate=1.0, ref_len=None): @@ -56,23 +59,24 @@ crst = crs.get('t') nv = v.unit( crs.get('nv') ) eyev = v.unit(l_g.v) - eyev_nega = v.nega(eyev) - ang_nv_eyev = v.dot_product( nv, eyev ) + eyev_nega = -eyev + ang_nv_eyev = np.dot( nv, eyev ) ang_nv_eyev_nega = -ang_nv_eyev rev = ang_nv_eyev_nega < 0 col = map_col( d, crs, rev, video ) + col = mt.arr(col) def diffusivity(): - diff = rtd.get('diff', 0.9) + diff = rtd.get('diff', 0.9) if not diff: return 0 def f(ltd): - ltv = ltd.get( 'v', [0,0,-1] ) + ltv = ltd.get( 'v', mt.arr( [0,0,-1] ) ) if 'p' in ltd: - ltv = v.sub_unit( crsp, ltd.get('p') ) - ang_nv_lt = v.dot_product( nv, ltv ) + ltv = v.unit( crsp - ltd.get('p') ) + ang_nv_lt = np.dot( nv, ltv ) ang_nv_lt_nega = -ang_nv_lt if ang_nv_eyev_nega * ang_nv_lt_nega < 0: return 0 @@ -81,7 +85,7 @@ return sum( map( f, lights ) ) * diff diff = diffusivity() - col = v.op1( '*', col, base + diff ) + col = col * (base + diff) def reflect_col(): reflect = rtd.get('reflect', 0) @@ -89,29 +93,29 @@ return [] if ang_nv_eyev_nega < 0 and not rtd.get('reflect_backside', False): return [] - ref_v = v.op2( '+', v.op1( '*', nv, 2 * ang_nv_eyev_nega ), eyev ) + ref_v = nv * ( 2 * ang_nv_eyev_nega ) + eyev ref_l = line.new( crsp, ref_v ) col = get_col( data, lights, ref_l, video, nest_rate * reflect, crst ) - return v.op1( '*', col, reflect ) + return col * reflect if col != [] else [] col_ = reflect_col() - if col_: - col = v.op2( '+', col, col_ ) + if col_ != []: + col = col + col_ def reflact_v(): density = rtd.get('density', 1) if density == 1 or abs(ang_nv_eyev) == 1 or ang_nv_eyev == 0: return eyev - cos_v = v.op1( '*', nv, ang_nv_eyev ) - sin_v = v.sub( eyev, cos_v ) + cos_v = nv * ang_nv_eyev + sin_v = eyev - cos_v if ang_nv_eyev < 0: density = 1.0 / density - r_sin_v = v.op1( '*', sin_v, density ) + r_sin_v = sin_v * density sin_len2 = v.len2( r_sin_v ) sin_len2 = min( sin_len2, 1 ) cos_len = math.sqrt( 1 - sin_len2 ) - r_cos_v = v.op1( '*', cos_v, cos_len / v.len1( cos_v ) ) - return v.unit( v.op2( '+', r_sin_v, r_cos_v ) ) + r_cos_v = cos_v * ( cos_len / v.len1( cos_v ) ) + return v.unit( r_sin_v + r_cos_v ) def reflact_col(): reflact = rtd.get('reflact', 0) @@ -120,23 +124,25 @@ ref_v = reflact_v() ref_l = line.new( crsp, ref_v ) col = get_col( data, lights, ref_l, video, nest_rate * reflact, crst ) - return v.op1( '*', col, reflact ) + return col * reflact if col != [] else [] col_ = reflact_col() - if col_: - col = v.op2( '+', col, col_ ) + if col_ != []: + col = col + col_ - return map( lambda v: min( int(v), 255 ), col ) + return col def draw(data, lights, eye2g, wh2eye, sc_sz, video): (w, h, d_) = sc_sz wh2g = lstx.opt( wh2eye + [ eye2g ] ) cnt = ut.new_cnt_wh(w, h, 'wh') - while cnt.up(): + show = '-quiet' not in sys.argv + while cnt.up(show): (ix, iy) = cnt.cur() - p = lstx.tr( wh2g, 'l2g', [ix,iy,-d_] ) + p = lstx.tr( wh2g, 'l2g', mt.arr( [ix,iy,-d_] ) ) l_g = line.new_p2(eye2g.p, p) col = get_col(data, lights, l_g, video) - if col: + col = map( lambda v: min( int(v), 255 ), col ) + if col != []: video.set(ix, iy, col) # EOF diff -urN rt_v12/ut.py rt_v13/ut.py --- rt_v12/ut.py 2018-03-23 11:00:56.000000000 +0900 +++ rt_v13/ut.py 2018-03-31 12:12:15.000000000 +0900 @@ -79,7 +79,21 @@ e.show_sec = None e.show_interval = 1.0 - def up(): + def up_no_show(): + if e.st_sec == None: + e.st_sec = now_sec() + else: + e.i += 1 + run = e.i < e.n + if not run: + now = now_sec() + e.total_sec = now - e.st_sec + show_stat(now) + return run + + def up(show=True): + if not show: + return up_no_show() now = now_sec() if e.st_sec == None: e.st_sec = now @@ -89,7 +103,7 @@ run = e.i < e.n if e.total_sec != None and run: e.total_sec = e.total_sec * n / e.i - show_stat() + show_stat(now) return run e.up = up @@ -105,8 +119,9 @@ stat_str = lambda : '{}{}/{}({}%) : {}'.format( id_str, e.i, n, 1000*e.i/n*0.1, fin_str() ) - def show_stat(): - now = now_sec() + def show_stat(now=None): + if not now: + now = now_sec() if e.show_sec == None or e.i >= n or now - e.show_sec >= e.show_interval: e.show_sec = now print stat_str() diff -urN rt_v12/v.py rt_v13/v.py --- rt_v12/v.py 2018-03-16 20:38:43.000000000 +0900 +++ rt_v13/v.py 2018-04-01 11:55:34.000000000 +0900 @@ -1,64 +1,24 @@ #!/usr/bin/env python import math +import numpy as np +import mt -zero = [0,0,0] -one = [1,1,1] -x1 = [1,0,0] -y1 = [0,1,0] -z1 = [0,0,1] -x_nega = [-1,0,0] -y_nega = [0,-1,0] -z_nega = [0,0,-1] -all = lambda e: [e,e,e] - -def op1(op, v, p): - ops = { - '+': lambda e: e + p, - '-': lambda e: e - p, - '*': lambda e: e * p, - '/': lambda e: e / float(p), - } - f = ops.get( op, lambda e: None ) - return map(f, v) - -def op2(op, va, vb): - ops = { - '+': lambda (a, b): a + b, - '-': lambda (a, b): a - b, - '*': lambda (a, b): a * b, - '/': lambda (a, b): a / float(b), - } - f = ops.get( op, lambda (a, b): None ) - return map( f, zip(va, vb) ) - -nega = lambda v: map( lambda e: -e, v ) - -lst_add = lambda vlst: map( sum, zip(*vlst) ) - -dot_product = lambda va, vb: sum( op2('*', va, vb) ) - -len2 = lambda v: dot_product(v, v) - -len1 = lambda v: math.sqrt( len2(v) ) - -unit = lambda v: ( lambda l: op1('/', v, l) if l != 0 else z1 )( len1(v) ) - -sub = lambda a, b: op2('-', a, b) - -sub_unit = lambda a, b: unit( sub(a, b) ) -sub_len = lambda a, b: len1( sub(a, b) ) - -def dot_product_len(v_targ, v_base): - base = dot_product(v_base, v_base) - return float( dot_product(v_targ, v_base) ) / base if base != 0 else 0 - -def cross_product(va, vb): - jks = map( lambda i: ( (i+1)%3, (i+2)%3 ), range(3) ) - return map( lambda (j, k): va[j] * vb[k] - va[k] * vb[j], jks ) +zero = mt.arr( [0,0,0] ) +one = mt.arr( [1,1,1] ) +x1 = mt.arr( [1,0,0] ) +y1 = mt.arr( [0,1,0] ) +z1 = mt.arr( [0,0,1] ) +all = lambda e: mt.arr( [e,e,e] ) -cross_product_unit = lambda va, vb: unit( cross_product(va, vb) ) +len2 = lambda v: np.dot(v,v) -p3_to_nv = lambda o, a, b: cross_product_unit( sub(a, o), sub(b, o) ) +len1 = lambda v: np.linalg.norm(v) + +unit = lambda v: ( lambda l: v / l if l != 0 else z1 )( np.linalg.norm( mt.to_arr(v) ) ) # ! + +cross_product_unit = lambda va, vb: unit( np.cross(va, vb) ) + +p3_to_nv = lambda o, a, b: cross_product_unit( a-o, b-o ) # EOF diff -urN rt_v12/vecs.py rt_v13/vecs.py --- rt_v12/vecs.py 2018-03-28 22:25:34.000000000 +0900 +++ rt_v13/vecs.py 2018-03-31 17:54:33.000000000 +0900 @@ -1,6 +1,7 @@ #!/usr/bin/env python import math +import numpy as np import ut import mt import v @@ -8,10 +9,12 @@ def new(v3): # v3: [ v.x1, v.y1, v.z1 ] + v3 = mt.to_arr(v3) + e = ut.Empty() e.typ = 'vecs' e.v3 = v3 - l2g = lambda v_: v.lst_add( map( lambda (a, e): v.op1('*', a, e), zip(v3, v_) ) ) + l2g = lambda v_: np.dot( v3.T, v_ ) g2l = lambda v: e.rev().tr('l2g', v) e.tr = lambda d, v: { 'l2g': l2g, 'g2l': g2l }.get(d, l2g)(v) @@ -22,7 +25,7 @@ return e.rev_vs e.rev = rev - e.zoom = lambda zm3: new( map( lambda (v_, zm): v.op1('*', v_, zm), zip(v3, zm3) ) ) + e.zoom = lambda zm3: new( v3 * mt.T( mt.to_arr(zm3) ) ) e.zoom_all = lambda zm: e.zoom( [zm,zm,zm] ) e.zoom_x = lambda zm: e.zoom( [zm,1,1] ) e.zoom_y = lambda zm: e.zoom( [1,zm,1] ) @@ -44,8 +47,8 @@ v3 = [ [],[],[] ] v3[fix_i] = fix - if fix == ref: - ref = ( lambda (x, y, z): (-y, x, z) )( ref ) if ref != v.z1 else v.x1 + if all( fix == ref ): + ref = ( lambda (x, y, z): (-y, x, z) )( ref ) if all( ref != v.z1 ) else v.x1 (a, b) = (fix, ref) if (fix_i + 1) % 3 == ref_i else (ref, fix) v3[i] = c = v.cross_product_unit(a, b) diff -urN rt_v12/way.py rt_v13/way.py --- rt_v12/way.py 2018-03-16 20:38:43.000000000 +0900 +++ rt_v13/way.py 2018-03-30 15:56:13.000000000 +0900 @@ -56,7 +56,7 @@ def liss(ps, pe, t3): lsts = map( lambda (s, e, t): [ (s, t*0.5), (e, t*0.5) ], zip(ps, pe, t3) ) - return lambda t, itp3=True: map( lambda lst: get(lst, t, itp3), lsts ) + return lambda t, itp3=True: mt.arr( map( lambda lst: get(lst, t, itp3), lsts ) ) def liss_eye(): # sys.argv @@ -75,8 +75,9 @@ def liss_get(key): (p, r, t) = p_r_t(key) - ps = v.op1('-', p, r) - pe = v.op1('+', p, r) + p = mt.arr(p) + ps = p - r + pe = p + r t3 = (t, t*4/3, t*5/3) return liss(ps, pe, t3) diff -urN rt_v12/wf.py rt_v13/wf.py --- rt_v12/wf.py 2018-03-25 12:01:36.000000000 +0900 +++ rt_v13/wf.py 2018-03-30 16:29:31.000000000 +0900 @@ -1,5 +1,6 @@ #!/usr/bin/env python +import mt import v import line import ax @@ -13,16 +14,16 @@ if kind == 'square': n = 4 - ps = map( lambda i: ax.rot_z(i*90).tr('l2g', [1,1,0]), range(n) ) + ps = map( lambda i: ax.rot_z(i*90).tr('l2g', mt.arr( [1,1,0] ) ), range(n) ) return rv_ps(ps) if kind == 'triangle': n = 3 - ps = [ [0,1,0],[-1,0,0],[1,0,0] ] + ps = mat.arr( [ [0,1,0],[-1,0,0],[1,0,0] ] ) return rv_ps(ps) if kind == 'poly_n': n = d.get('n', 3) deg = 360.0 / n - ps = map( lambda i: ax.rot_z( (i+0.5)*deg ).tr('l2g', [0,-1,0]), range(n) ) + ps = map( lambda i: ax.rot_z( (i+0.5)*deg ).tr('l2g', mt.arr( [0,-1,0] ) ), range(n) ) return rv_ps(ps) if kind == 'circle': d = { 'kind': 'poly_n', 'n': 16 } @@ -61,7 +62,7 @@ def mk_pyramid(ps): n = len(ps) idxs = len_to_idxs(n) + map( lambda i: (i,n), range(n) ) - return ( ps + [ [0,0,1] ], idxs ) + return ( ps + [ mt.arr( [0,0,1] ) ], idxs ) if kind =='poly_n_pyramid': d = { 'kind': 'poly_n', 'n': d.get('n') } @@ -91,10 +92,10 @@ def eye_clip(eye2g, wh2eye, sc_sz, pa, pb): (w, h, d) = sc_sz n = 4 - p4 = [ [0,0,-d], [w,0,-d], [w,h,-d], [0,h,-d] ] + p4 = mt.arr( [ [0,0,-d], [w,0,-d], [w,h,-d], [0,h,-d] ] ) p4 = map( lambda p: lstx.tr(wh2eye, 'l2g', p), p4 ) for i in range(n): - nv = v.p3_to_nv( [0,0,0], p4[i], p4[(i+1)%n] ) + nv = v.p3_to_nv( mt.arr( [0,0,0] ), p4[i], p4[(i+1)%n] ) r = nv_clip(nv, pa, pb) if not r: return [] @@ -102,7 +103,7 @@ return [pa, pb] def draw(data, lights, eye2g, wh2eye, sc_sz, video): - col = v.all(255) # white + col = [255,255,255] # white for d in data: ps = [] idxs = []