From dfbfbe19f2a46a86402de35c23a93a370d8dee21 Mon Sep 17 00:00:00 2001 From: hodapp Date: Tue, 17 Sep 2019 12:29:41 -0400 Subject: [PATCH] Add transform stack with matrices --- Scratch.ipynb | 69 ++++++++++++++++++++++++++++++++++++++------------- meshutil.py | 52 +++++++++++++++++++++++++++++++++++--- quat.py | 11 +++++++- 3 files changed, 111 insertions(+), 21 deletions(-) diff --git a/Scratch.ipynb b/Scratch.ipynb index 9a9493d..2192455 100644 --- a/Scratch.ipynb +++ b/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.v = c2.v/2 + [0,1,0]\n", - "c = c.concat(c2)\n", + "#c2 = meshutil.cube_distort(0.3)\n", + "#c2.v = c2.v/2 + [0,1,0]\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", - "/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" + "face_normals didn't match triangles, ignoring!\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 @@ "" ] }, - "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, diff --git a/meshutil.py b/meshutil.py index 1809825..b3d8d53 100644 --- a/meshutil.py +++ b/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: diff --git a/quat.py b/quat.py index d8cbb90..7f7c9bb 100644 --- a/quat.py +++ b/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], + ])