Add transform stack with matrices

Scratch.ipynb

@@ -4,11 +4,28 @@
    "cell_type": "code",
    "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/usr/local/lib/python3.7/site-packages/quaternion/calculus.py:475: UserWarning: \n",
+      "\n",
+      "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
+      "Could not import from scipy, which means that derivatives\n",
+      "and integrals will use less accurate finite-differencing\n",
+      "techniques.  You may want to install scipy.\n",
+      "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
+      "\n",
+      "  warnings.warn(warning_text)\n"
+     ]
+    }
+   ],
    "source": [
     "import meshutil\n",
     "import stl.mesh\n",
-    "import numpy"
+    "import numpy\n",
+    "import trimesh"
    ]
   },
   {
@@ -19,9 +36,9 @@
    "source": [
     "c = meshutil.cube(open_xz=False)\n",
     "mesh_fname = \"cube_test.stl\"\n",
-    "c2 = meshutil.cube_distort(0.3)\n",
+    "#c2 = meshutil.cube_distort(0.3)\n",
-    "c2.v = c2.v/2 + [0,1,0]\n",
+    "#c2.v = c2.v/2 + [0,1,0]\n",
-    "c = c.concat(c2)\n",
+    "#c = c.concat(c2)\n",
     "m = c.to_stl_mesh()\n",
     "\"\"\"\n",
     "c1 = meshutil.cube_distort(0.3)\n",
@@ -38,34 +55,52 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 9,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 4.,  0.,  0., -4.],\n",
+       "       [ 0.,  4.,  0.,  0.],\n",
+       "       [ 0.,  0.,  4.,  0.],\n",
+       "       [ 0.,  0.,  0.,  1.]])"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "#d = numpy.concatenate([c,c1,c2,c3,c4])"
+    "s = meshutil.TransformStack()\n",
+    "mtx = s.translate(-1, 0, 0).scale(4).rotate().mtx\n",
+    "#mtx = meshutil.mtx_rotate(numpy.array([0,1,0]), numpy.pi/4) @ meshutil.mtx_scale(2) @ meshutil.mtx_translate(-1, 0, 0)\n",
+    "#mtx = meshutil.mtx_translate(-1, 0, 0)\n",
+    "mtx"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import trimesh"
+    "ct = c.transform(mtx)\n",
+    "c2 = ct.concat(c)\n",
+    "c2.to_stl_mesh().save(mesh_fname)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "face_normals didn't match triangles, ignoring!\n",
+      "face_normals didn't match triangles, ignoring!\n"
-      "/home/hodapp/.local/lib/python3.6/site-packages/IPython/core/display.py:694: UserWarning: Consider using IPython.display.IFrame instead\n",
-      "  warnings.warn(\"Consider using IPython.display.IFrame instead\")\n"
      ]
     },
     {
@@ -263,7 +298,7 @@
        "Math.pow(distanceToCamera,2)+\n",
        "Math.pow(distanceToCamera,2));camera.position.set(len,len,len);controls.update();camera.lookAt(boundingSphere.center);controls.target.set(boundingSphere.center.x,boundingSphere.center.y,boundingSphere.center.z);camera.updateProjectionMatrix();}\n",
        "function centerControls(obj,camera,controls){const boundingBox=new THREE.Box3().setFromObject(obj);const boundingSphere=new THREE.Sphere();boundingBox.getBoundingSphere((target=boundingSphere));controls.update();controls.target.set(boundingSphere.center.x,boundingSphere.center.y,boundingSphere.center.z);}\n",
-       "function init(){scene=new THREE.Scene();scene.background=new THREE.Color(0xffffff);tracklight=new THREE.DirectionalLight(0xffffff,1.75);scene.add(tracklight);base64_data="Z2xURgIAAABwBQAAsAMAAEpTT057InNjZW5lIjogMCwgInNjZW5lcyI6IFt7Im5vZGVzIjogWzBdfV0sICJhc3NldCI6IHsidmVyc2lvbiI6ICIyLjAiLCAiZ2VuZXJhdG9yIjogImh0dHBzOi8vZ2l0aHViLmNvbS9taWtlZGgvdHJpbWVzaCJ9LCAiYWNjZXNzb3JzIjogW3siYnVmZmVyVmlldyI6IDAsICJjb21wb25lbnRUeXBlIjogNTEyNSwgImNvdW50IjogNjAsICJtYXgiOiBbMTRdLCAibWluIjogWzBdLCAidHlwZSI6ICJTQ0FMQVIifSwgeyJidWZmZXJWaWV3IjogMSwgImNvbXBvbmVudFR5cGUiOiA1MTI2LCAiY291bnQiOiAxNSwgInR5cGUiOiAiVkVDMyIsICJieXRlT2Zmc2V0IjogMCwgIm1heCI6IFsxLjAsIDEuNzYyMDIyNDk1MjY5Nzc1NCwgMS4wXSwgIm1pbiI6IFstMC4xNzAwOTE4NTI1NDU3MzgyMiwgMC4wLCAtMC4xNzAwOTE4NTI1NDU3MzgyMl19XSwgIm1lc2hlcyI6IFt7Im5hbWUiOiAiY3ViZV90ZXN0LnN0bCIsICJwcmltaXRpdmVzIjogW3siYXR0cmlidXRlcyI6IHsiUE9TSVRJT04iOiAxfSwgImluZGljZXMiOiAwLCAibW9kZSI6IDR9XX1dLCAiY2FtZXJhcyI6IFt7Im5hbWUiOiAiY2FtZXJhXzBHUVdIQyIsICJ0eXBlIjogInBlcnNwZWN0aXZlIiwgInBlcnNwZWN0aXZlIjogeyJhc3BlY3RSYXRpbyI6IDEuMzMzMzMzMzMzMzMzMzMzMywgInlmb3YiOiAwLjc4NTM5ODE2MzM5NzQ0ODMsICJ6bmVhciI6IDAuMDF9fV0sICJub2RlcyI6IFt7Im5hbWUiOiAid29ybGQiLCAiY2hpbGRyZW4iOiBbMV19LCB7Im5hbWUiOiAiY3ViZV90ZXN0LnN0bF8wVkZQWUxVMENMWFMiLCAibWVzaCI6IDB9XSwgImJ1ZmZlcnMiOiBbeyJieXRlTGVuZ3RoIjogNDIwfV0sICJidWZmZXJWaWV3cyI6IFt7ImJ1ZmZlciI6IDAsICJieXRlT2Zmc2V0IjogMCwgImJ5dGVMZW5ndGgiOiAyNDB9LCB7ImJ1ZmZlciI6IDAsICJieXRlT2Zmc2V0IjogMjQwLCAiYnl0ZUxlbmd0aCI6IDE4MH1dfSAgIKQBAABCSU4AAAAAAAcAAAAFAAAAAAAAAAIAAAAHAAAABQAAAAgAAAAGAAAABQAAAAcAAAAIAAAABgAAAAMAAAABAAAABgAAAAgAAAADAAAAAQAAAAIAAAAAAAAAAQAAAAMAAAACAAAAAgAAAAgAAAAHAAAAAgAAAAMAAAAIAAAAAAAAAAUAAAAGAAAAAAAAAAYAAAABAAAAAgAAAAkAAAALAAAAAgAAAAoAAAAJAAAACwAAAA0AAAAMAAAACwAAAAkAAAANAAAADAAAAA4AAAAEAAAADAAAAA0AAAAOAAAABAAAAAoAAAACAAAABAAAAA4AAAAKAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIA/AAAAAAAAgD8AAAAAAAAAAAAAgD8AAIA/AAAAAAAAgD8AAAA/AACAPwAAAAAAAAAAAACAPwAAAAAAAIA/AACAPwAAgD8AAAAAAACAPwAAgD8AAIA/QaGuPiegzj+PLC6+bU4Xvlq2uz9tThe+AAAAPwAAgD8AAAAAAAAAPwAAgD8AAAA/MDKjPvSJ4T8wMqM+jywuviegzj9Boa4+";;renderer=new THREE.WebGLRenderer({antialias:true});renderer.setPixelRatio(window.devicePixelRatio);renderer.setSize(window.innerWidth,window.innerHeight);document.body.appendChild(renderer.domElement);loader=new THREE.GLTFLoader();loader.load("data:text/plain;base64,"+base64_data,function(gltf){scene.add(gltf.scene);camera=gltf.cameras[0];controls=new THREE.TrackballControls(camera);controls.rotateSpeed=1.0;controls.zoomSpeed=1.2;controls.panSpeed=0.8;controls.noZoom=false;controls.noPan=false;controls.staticMoving=true;controls.dynamicDampingFactor=0.3;controls.keys=[65,83,68];controls.addEventListener("change",render);centerControls(scene,camera,controls);render();window.addEventListener("resize",onWindowResize,false);animate();});}\n",
+       "function init(){scene=new THREE.Scene();scene.background=new THREE.Color(0xffffff);tracklight=new THREE.DirectionalLight(0xffffff,1.75);scene.add(tracklight);base64_data="Z2xURgIAAABwBQAAgAMAAEpTT057InNjZW5lIjogMCwgInNjZW5lcyI6IFt7Im5vZGVzIjogWzBdfV0sICJhc3NldCI6IHsidmVyc2lvbiI6ICIyLjAiLCAiZ2VuZXJhdG9yIjogImh0dHBzOi8vZ2l0aHViLmNvbS9taWtlZGgvdHJpbWVzaCJ9LCAiYWNjZXNzb3JzIjogW3siYnVmZmVyVmlldyI6IDAsICJjb21wb25lbnRUeXBlIjogNTEyNSwgImNvdW50IjogNzIsICJtYXgiOiBbMTRdLCAibWluIjogWzBdLCAidHlwZSI6ICJTQ0FMQVIifSwgeyJidWZmZXJWaWV3IjogMSwgImNvbXBvbmVudFR5cGUiOiA1MTI2LCAiY291bnQiOiAxNSwgInR5cGUiOiAiVkVDMyIsICJieXRlT2Zmc2V0IjogMCwgIm1heCI6IFsxLjAsIDQuMCwgNC4wXSwgIm1pbiI6IFstNC4wLCAwLjAsIDAuMF19XSwgIm1lc2hlcyI6IFt7Im5hbWUiOiAiY3ViZV90ZXN0LnN0bCIsICJwcmltaXRpdmVzIjogW3siYXR0cmlidXRlcyI6IHsiUE9TSVRJT04iOiAxfSwgImluZGljZXMiOiAwLCAibW9kZSI6IDR9XX1dLCAiY2FtZXJhcyI6IFt7Im5hbWUiOiAiY2FtZXJhX01MOVJUSiIsICJ0eXBlIjogInBlcnNwZWN0aXZlIiwgInBlcnNwZWN0aXZlIjogeyJhc3BlY3RSYXRpbyI6IDEuMzMzMzMzMzMzMzMzMzMzMywgInlmb3YiOiAwLjc4NTM5ODE2MzM5NzQ0ODMsICJ6bmVhciI6IDAuMDF9fV0sICJub2RlcyI6IFt7Im5hbWUiOiAid29ybGQiLCAiY2hpbGRyZW4iOiBbMV19LCB7Im5hbWUiOiAiY3ViZV90ZXN0LnN0bF9NMVlPU082U0RIQVAiLCAibWVzaCI6IDB9XSwgImJ1ZmZlcnMiOiBbeyJieXRlTGVuZ3RoIjogNDY4fV0sICJidWZmZXJWaWV3cyI6IFt7ImJ1ZmZlciI6IDAsICJieXRlT2Zmc2V0IjogMCwgImJ5dGVMZW5ndGgiOiAyODh9LCB7ImJ1ZmZlciI6IDAsICJieXRlT2Zmc2V0IjogMjg4LCAiYnl0ZUxlbmd0aCI6IDE4MH1dfSAgINQBAABCSU4ABwAAAAMAAAAAAAAABwAAAAkAAAADAAAAAAAAAAQAAAABAAAAAAAAAAMAAAAEAAAAAQAAAAoAAAAIAAAAAQAAAAQAAAAKAAAACAAAAAkAAAAHAAAACAAAAAoAAAAJAAAACQAAAAQAAAADAAAACQAAAAoAAAAEAAAABwAAAAAAAAABAAAABwAAAAEAAAAIAAAAAAAAAA0AAAALAAAAAAAAAAUAAAANAAAACwAAAA4AAAAMAAAACwAAAA0AAAAOAAAADAAAAAYAAAACAAAADAAAAA4AAAAGAAAAAgAAAAUAAAAAAAAAAgAAAAYAAAAFAAAABQAAAA4AAAANAAAABQAAAAYAAAAOAAAAAAAAAAsAAAAMAAAAAAAAAAwAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIBAAAAAAAAAAAAAAIA/AAAAAAAAgEAAAAAAAAAAAAAAgEAAAIBAAAAAAAAAgD8AAAAAAAAAAAAAgD8AAIA/AACAwAAAAAAAAAAAAACAwAAAAAAAAIBAAACAwAAAgEAAAAAAAACAwAAAgEAAAIBAAACAPwAAAAAAAAAAAACAPwAAAAAAAIA/AACAPwAAgD8AAAAAAACAPwAAgD8AAIA/";;renderer=new THREE.WebGLRenderer({antialias:true});renderer.setPixelRatio(window.devicePixelRatio);renderer.setSize(window.innerWidth,window.innerHeight);document.body.appendChild(renderer.domElement);loader=new THREE.GLTFLoader();loader.load("data:text/plain;base64,"+base64_data,function(gltf){scene.add(gltf.scene);camera=gltf.cameras[0];controls=new THREE.TrackballControls(camera);controls.rotateSpeed=1.0;controls.zoomSpeed=1.2;controls.panSpeed=0.8;controls.noZoom=false;controls.noPan=false;controls.staticMoving=true;controls.dynamicDampingFactor=0.3;controls.keys=[65,83,68];controls.addEventListener("change",render);centerControls(scene,camera,controls);render();window.addEventListener("resize",onWindowResize,false);animate();});}\n",
        "function onWindowResize(){camera.aspect=window.innerWidth/window.innerHeight;camera.updateProjectionMatrix();renderer.setSize(window.innerWidth,window.innerHeight);controls.handleResize();render();}\n",
        "function animate(){requestAnimationFrame(animate);controls.update();}\n",
        "function render(){tracklight.position.copy(camera.position);renderer.render(scene,camera);}\n",
@@ -274,7 +309,7 @@
        "<IPython.core.display.HTML object>"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 11,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -308,7 +343,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.2"
   }
  },
  "nbformat": 4,

meshutil.py

@@ -18,7 +18,7 @@ rtb = numpy.array([1,1,1])
 
 class FaceVertexMesh(object):
     def __init__(self, v, f):
-        # v & f should both be of shape (N,3)ll
+        # v & f should both be of shape (N,3)
         self.v = v
         self.f = f
     def concat(self, other_mesh):
@@ -27,6 +27,12 @@ class FaceVertexMesh(object):
         f2 = numpy.concatenate([self.f, other_mesh.f + self.v.shape[0]])
         m2 = FaceVertexMesh(v2, f2)
         return m2
+    def transform(self, mtx):
+        vh = numpy.hstack([self.v, numpy.ones((self.v.shape[0], 1), dtype=self.v.dtype)])
+        v2 = vh.dot(mtx.T)[:,0:3]
+        # TODO: why transpose?
+        # TODO: fix homogenous (even if I don't use it)
+        return FaceVertexMesh(v2, self.f)
     def to_stl_mesh(self):
         data = numpy.zeros(self.f.shape[0], dtype=stl.mesh.Mesh.dtype)
         v = data["vectors"]
@@ -34,11 +40,51 @@ class FaceVertexMesh(object):
             v[i] = [self.v[iv0], self.v[iv1], self.v[iv2]]
         return stl.mesh.Mesh(data)
 
+class TransformStack(object):
+    def __init__(self, mtx=None):
+        if mtx is None:
+            self.mtx = numpy.identity(4)
+        else:
+            self.mtx = mtx
+    def scale(self, *a, **kw):
+        mtx = mtx_scale(*a, **kw)
+        return TransformStack(mtx @ self.mtx)
+    def translate(self, *a, **kw):
+        mtx = mtx_translate(*a, **kw)
+        return TransformStack(mtx @ self.mtx)
+    def rotate(self, *a, **kw):
+        mtx = mtx_rotate(*a, **kw)
+        return TransformStack(mtx @ self.mtx)
+
+def mtx_scale(sx, sy=None, sz=None):
+    if sy is None:
+        sy = sx
+    if sz is None:
+        sz = sx
+    return numpy.array([
+        [sx, 0,  0,  0],
+        [0, sy,  0,  0],
+        [0,  0, sz,  0],
+        [0,  0,  0,  1],
+    ])
+
+def mtx_translate(x, y, z):
+    return numpy.array([
+        [1, 0, 0, x],
+        [0, 1, 0, y],
+        [0, 0, 1, z],
+        [0, 0, 0, 1],
+    ])
+
+def mtx_rotate(axis, angle):
+    q = quat.rotation_quaternion(axis, angle)
+    return quat.quat2mat(q)
+
 def cube(open_xz=False):
     verts = numpy.array([
         lbf, rbf, ltf, rtf,
         lbb, rbb, ltb, rtb,
-    ], dtype=numpy.float32)
+    ], dtype=numpy.float64)
     if open_xz:
         faces = numpy.zeros((8,3), dtype=int)
     else:
@@ -69,7 +115,7 @@ def cube_distort(angle):
     verts = numpy.array([
         lbf, rbf, ltf2, rtf2,
         lbb, rbb, ltb2, rtb2,
-    ], dtype=numpy.float32)
+    ], dtype=numpy.float64)
     if True: #open_xz:
         faces = numpy.zeros((8,3), dtype=int)
     else:

quat.py

@@ -15,10 +15,19 @@ def rotation_quaternion(axis, angle):
     """
     qc = numpy.cos(angle / 2)
     qs = numpy.sin(angle / 2)
-    qv = qs * axis
+    qv = qs * numpy.array(axis)
     return numpy.quaternion(qc, qv[0], qv[1], qv[2])
 
 def vec2quat(vs):
     qs = numpy.zeros(vs.shape[0], dtype=numpy.quaternion)
     quaternion.as_float_array(qs)[:,1:4] = vs
     return qs
+
+def quat2mat(q):
+    s = 1
+    return numpy.array([
+        [1-2*s*(q.y**2+q.z**2), 2*s*(q.x*q.y-q.z*q.w), 2*s*(q.x*q.z+q.y*q.w), 0],
+        [2*s*(q.x*q.y+q.z*q.w), 1-2*s*(q.x**2+q.z**2), 2*s*(q.y*q.z-q.x*q.w), 0],
+        [2*s*(q.x*q.z-q.y*q.w), 2*s*(q.y*q.z+q.x*q.w), 1-2*s*(q.x**2+q.y**2), 0],
+        [0, 0, 0, 1],
+    ])