diff --git a/.gitignore b/.gitignore
index 3001f4b..1e3ed9c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -71,4 +71,5 @@ build
!asset/armadillo.obj
-venv
\ No newline at end of file
+venv
+task09/output.png
\ No newline at end of file
diff --git a/README.md b/README.md
index 61ce5db..d02a69e 100644
--- a/README.md
+++ b/README.md
@@ -55,8 +55,8 @@ Topics:
|(8)
June 3| **Character animation**
Linear blend skinning | [task07](task07) | [[21]](http://nobuyuki-umetani.com/acg2024s/jacobian.pdf) |
|(9)
June 10| Guest lecture by Dr. Rex West | | |
|(10)
June 17| **Character animation2**
Inverse kinematic | [task08](task08) | [[20]](http://nobuyuki-umetani.com/acg2024s/character_deformation.pdf) |
-|(11)
June 24| **Optimization** | | |
-|(12)
July 12| **Laplacian mesh deformation**
| task09 | |
+|(11)
June 24| **Optimization**
Newton-Raphson method, gradient descent | | [[22]](http://nobuyuki-umetani.com/acg2024s/optimization.pdf) |
+|(12)
July 1| **Laplacian mesh deformation**
Sparse linear system | [task09](task09) | [[23]](http://nobuyuki-umetani.com/acg2024s/mesh_laplacian.pdf) [[24]](http://nobuyuki-umetani.com/acg2024s/sparse_linear_system.pdf) |
|(13)
July 8| **Grid-based Fluid**
Poisson equation, Stam fluid | | |
@@ -85,8 +85,8 @@ Look at the following document.
| [task05](task05) | **Fragment shader practice**
Ray marching method, CSG modeling, implicit modeling |
|
| [task06](task06) | **Monte Carlo integration1**
Ambient occlusion, importance sampling, BVH |
|
| [task07](task07) | **Monte Carlo integration2**
Multiple importance sampling |
|
-| task08 | **Skeletal Character Animation**
Linear blend skinning, articulated rigid body |
|
-| task09 | TBD | |
+| [task08](task08) | **Skeletal Character Animation**
Linear blend skinning, articulated rigid body |
|
+| [task09](task09) | **Laplacian Mesh Deformation**
Quadratic programming, sparse linear system |
|
### Policy
@@ -136,9 +136,10 @@ Look at the following document.
- [[19]Ray Triangle Collision](http://nobuyuki-umetani.com/acg2024s/ray_triangle_collision.pdf)
- [[20]Character deformation](http://nobuyuki-umetani.com/acg2024s/character_deformation.pdf)
-- [[21]Jacobian&Hessian](http://nobuyuki-umetani.com/acg2024s/jacobian.pdf)
-
-
+- [[21] Jacobian&Hessian](http://nobuyuki-umetani.com/acg2024s/jacobian.pdf)
+- [[22] Optimization](http://nobuyuki-umetani.com/acg2024s/optimization.pdf)
+- [[23] Mesh laplacian](http://nobuyuki-umetani.com/acg2024s/mesh_laplacian.pdf)
+- [[24] Sparse linear system](http://nobuyuki-umetani.com/acg2024s/sparse_linear_system.pdf)
diff --git a/task09/README.md b/task09/README.md
new file mode 100644
index 0000000..b45bcd9
--- /dev/null
+++ b/task09/README.md
@@ -0,0 +1,110 @@
+# Task09: Laplacian Mesh Deformation (Quadratic Programming, Sparse Matrix)
+
+
+
+**Deadline: July 4th (Thu) at 15:00pm**
+
+----
+
+## Before Doing Assignment
+
+### Install Python (if necessary)
+We use Python for this assignment.
+This assigment only supports Python ver. 3.
+
+To check if Python 3.x is installed, launch a command prompt and type `python3 --version` and see the version.
+
+For MacOS and Ubuntu you have Python installed by default.
+For Windows, you may need to install the Python by yourself.
+[This document](https://docs.python.org/3/using/windows.html) show how to install Python 3.x on Windows.
+
+
+### Virtual environment
+
+We want to install dependency ***locally*** for this assignment.
+
+```bash
+cd acg-
+python3 -m venv venv # make a virtual environment named "venv"
+```
+
+Then, start the virtual environment.
+For Mac or Linux, type
+
+```bash
+source venv/bin/activate # start virtual environment
+```
+
+For Windows, type
+
+```bash
+venv\Scripts\activate.bat # start virtual environment
+```
+
+In the command prompt, you will see `(venv)` at the beginning of each line.
+There will be `venv` folder under `acg-`.
+
+### Install dependency
+
+In this assignment we use many external library. We use `pip` to install these.
+
+```bash
+pip3 install numpy
+pip3 install moderngl
+pip3 install moderngl_window
+pip3 install scipy
+```
+
+Alternatively, you can install above dependency at once by
+
+```bash
+cd acg-/task09
+pip3 install -r requirements.txt
+```
+
+type `pip3 list` and then confirm you have libraries such as `moderngl`, `numpy`, `pillow`, `pyglet`, `pyrr`, `scipy` etc.
+
+### Make branch
+
+Follow [this document](../doc/submit.md) to submit the assignment, In a nutshell, before doing the assignment,
+- make sure you synchronized the `main ` branch of your local repository to that of remote repository.
+- make sure you created branch `task09` from `main` branch.
+- make sure you are currently in the `task09` branch (use `git branch -a` command).
+
+Now you are ready to go!
+
+---
+
+## Problem 1 (Python execution practice)
+
+Run the code with `python3 main.py`
+
+The program will output `out.png`. Rename it to `out_default.png` then it replaces the image below.
+
+
+
+
+## Problem 2 (smooth deformation with Laplacian)
+
+Write a few lines of code around line #134 to implement smooth mesh deformation using Laplacian.
+
+The program will output `output.png`. Rename it to `out_laplacian.png` then it replaces the image below.
+
+
+
+
+## Problem 3 (even smoother deformation with Bi-Laplacian)
+
+Write a few lines of code around line #134 to implement smooth mesh deformation using BiLaplacian.
+
+The program will output `out.png`. Rename it to `out_bilaplacian.png` then it replaces the image below.
+
+
+
+
+## After Doing the Assignment
+
+After modify the code, push the code and submit a pull request. Make sure your pull request only contains the files you edited. Good luck!
+
+BTW, You can exit the virtual environment by typing `deactivate` in the command prompt.
+
diff --git a/task09/main.py b/task09/main.py
new file mode 100644
index 0000000..63702b8
--- /dev/null
+++ b/task09/main.py
@@ -0,0 +1,210 @@
+import os
+import math
+
+import numpy
+#
+import pyrr
+import numpy as np
+import moderngl
+import moderngl_window as mglw
+from PIL import Image, ImageOps
+from scipy import sparse
+from scipy.sparse.linalg import spsolve
+#
+import util_for_task09
+
+
+def matrix_graph_laplacian(tri2vtx, num_vtx):
+ '''
+ function to make graph laplacian matrix
+ :param tri2vtx: index of vertex for triangles
+ :param num_vtx: number of vertex
+ :return: sparse matrix using scipy.sparse
+ '''
+ tri2ones = np.ones((tri2vtx.shape[0],))
+ W = np.empty(0, np.float32)
+ I = np.empty(0, np.uint32)
+ J = np.empty(0, np.uint32)
+ #
+ for (i, j) in ((0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)):
+ tri2vi = tri2vtx[:, i]
+ tri2vj = tri2vtx[:, j]
+ if i == j:
+ coeff = 1.
+ else:
+ coeff = -0.5
+ W = np.append(W, tri2ones * coeff)
+ I = np.append(I, tri2vi)
+ J = np.append(J, tri2vj)
+ return sparse.csr_matrix((W, (I, J)), shape=(num_vtx, num_vtx))
+
+
+class HelloWorld(mglw.WindowConfig):
+ '''
+ Window to show the mesh deformation
+ '''
+ gl_version = (3, 3)
+ title = "task09: laplacian mesh deformation"
+ window_size = (500, 500)
+ aspect_ratio = float(window_size[0]) / float(window_size[1])
+ resizable = False
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ self.is_screenshot_taken = False
+
+ # initialize mesh
+ asset_path = os.path.normpath(os.path.join(__file__, '../../asset/bunny.obj'))
+ self.tri2vtx, self.vtx2xyz_ini = util_for_task09.load_triangle_mesh_from_wavefront_obj_file(asset_path)
+ util_for_task09.normalize_vtx2xyz(self.vtx2xyz_ini)
+ self.vtx2xyz_def = self.vtx2xyz_ini.copy()
+
+ # make a list for fixed vertices
+ self.fixbase2vtx = (self.vtx2xyz_ini[:, 2] < -0.48).nonzero()[0].tolist()
+ self.fixear2vtx = (self.vtx2xyz_ini[:, 2] > +0.48).nonzero()[0].tolist()
+ self.fixback2vtx = (np.logical_and(
+ self.vtx2xyz_ini[:, 0] > +0.15,
+ self.vtx2xyz_ini[:, 2] > +0.115)).nonzero()[0].tolist()
+
+ # make laplacian, bi-laplacian matrices
+ num_vtx = self.vtx2xyz_ini.shape[0]
+ self.matrix_laplace = matrix_graph_laplacian(self.tri2vtx, num_vtx)
+ self.matrix_bilaplace = self.matrix_laplace * self.matrix_laplace
+
+ # make matrix for fixed points
+ coeff_penalty = 100.0
+ fix2vtx = numpy.concatenate([self.fixback2vtx, self.fixear2vtx, self.fixbase2vtx])
+ self.matrix_fix = sparse.csr_matrix((
+ np.ones(fix2vtx.shape[0], np.float32) * coeff_penalty, (fix2vtx, fix2vtx)),
+ shape=(num_vtx, num_vtx))
+
+ # -----------------------------
+ # below: visualization related
+
+ self.prog = self.ctx.program(
+ vertex_shader='''
+ #version 330
+ uniform mat4 matrix;
+ in vec3 in_vert;
+ in vec3 in_nrm;
+ out vec3 out_nrm;
+ void main() {
+ out_nrm = normalize((matrix * vec4(in_nrm, 0.0)).xyz);
+ gl_Position = matrix * vec4(in_vert, 1.0);
+ }
+ ''',
+ fragment_shader='''
+ #version 330
+ uniform bool is_shading; // variable of the program
+ uniform vec3 color;
+ out vec4 f_color;
+ in vec3 out_nrm;
+ void main() {
+ float ratio = abs(out_nrm.z);
+ if( !is_shading ){ ratio = 1.0; }
+ f_color = vec4(color*ratio, 1.0);
+ }
+ ''',
+ )
+
+ # initialize visualization for mesh
+ ebo = self.ctx.buffer(self.tri2vtx) # send triangle index data to GPU (element buffer object)
+ self.vbo_vtx2xyz_def = self.ctx.buffer(self.vtx2xyz_def) # send initial vertex coordinates data to GPU
+ vtx2nrm = util_for_task09.vtx2nrm(self.tri2vtx, self.vtx2xyz_def)
+ self.vbo_vtx2nrm_def = self.ctx.buffer(vtx2nrm) # send initial vertex coordinates data to GPU
+ self.vao_ini = self.ctx.vertex_array(
+ self.prog, [
+ (self.vbo_vtx2xyz_def, '3f', 'in_vert'),
+ (self.vbo_vtx2nrm_def, '3f', 'in_nrm')],
+ ebo, 4, mode=moderngl.TRIANGLES) # tell gpu about the mesh information and how to draw it
+
+ # initialize visualization for points
+ tri2vtx_octa, vtx2xyz_octa = util_for_task09.octahedron()
+ self.vao_octa = self.ctx.vertex_array(
+ self.prog, [
+ (self.ctx.buffer(vtx2xyz_octa), '3f', 'in_vert')],
+ self.ctx.buffer(tri2vtx_octa), 4, mode=moderngl.TRIANGLES)
+
+ def render(self, time, frame_time):
+ self.vtx2xyz_def[:] = self.vtx2xyz_ini[:]
+ # set fixed deformation
+ self.vtx2xyz_def[self.fixear2vtx, 0] += -0.5 * math.sin(time)
+ self.vtx2xyz_def[self.fixback2vtx, 2] += -0.3 * math.sin(2. * time)
+
+ # write a few line code below to implement laplacian mesh deformation
+ # Problem 2: Laplacian deformation, which minimizes (x-x_def)D(x-x_def) + (x-x_ini)L(x-x_ini) w.r.t x,
+ # Problem 3: Bi-Laplacian deformation, which minimizes (x-x_def)D(x-x_def) + (x-x_ini)L^2(x-x_ini) w.r.t x,
+ # where D is the diagonal matrix with entry 1 if the vertex is fixed, 0 if the vertex is free, a.k.a `self.matrix_fix`
+ # L is the graph Laplacian matrix a.k.a `self.matrix_laplace`
+ # you may use `spsolve` to solve the liner system
+ # spsolve: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve
+
+
+ # do not edit beyond here
+ # above: deformation
+ # ---------------------
+ # below: visualization
+
+ self.ctx.clear(1.0, 1.0, 1.0) # initizlize frame buffer with white
+ self.ctx.enable(moderngl.DEPTH_TEST)
+
+ self.vbo_vtx2xyz_def.write(self.vtx2xyz_def) # send vertex information to GPU
+ vtx2nrm = util_for_task09.vtx2nrm(self.tri2vtx, self.vtx2xyz_def)
+ self.vbo_vtx2nrm_def.write(vtx2nrm) # send normal information to GPU
+
+ # make view transformation
+ view_rot_z = pyrr.Matrix44.from_z_rotation(-np.pi * 0.2)
+ view_rot_x = pyrr.Matrix44.from_x_rotation(np.pi * 0.4)
+ view_z_flip = pyrr.Matrix44.from_scale((1.3, 1.3, -1.3, 1.))
+ view_transf = view_z_flip * view_rot_x * view_rot_z
+
+ # render triangle mesh
+ self.prog['matrix'].value = tuple(view_transf.flatten())
+ self.prog['color'].value = (0.8, 0.8, 0.8)
+ self.prog['is_shading'].value = 1
+ self.vao_ini.render()
+
+ r = 0.02
+ self.prog['is_shading'].value = 0
+
+ # render fixed bases
+ for i_vtx in self.fixbase2vtx:
+ pos = self.vtx2xyz_def[i_vtx].copy()
+ model_transf = pyrr.Matrix44.from_translation(pos) * pyrr.Matrix44.from_scale((r, r, r, 1.))
+ self.prog['matrix'].value = tuple((view_transf * model_transf).flatten())
+ self.prog['color'].value = (0.0, 0.0, 1.0)
+ self.vao_octa.render()
+
+ # render fixed ear
+ for i_vtx in self.fixear2vtx:
+ pos = self.vtx2xyz_def[i_vtx].copy()
+ model_transf = pyrr.Matrix44.from_translation(pos) * pyrr.Matrix44.from_scale((r, r, r, 1.))
+ self.prog['matrix'].value = tuple((view_transf * model_transf).flatten())
+ self.prog['color'].value = (0.0, 1.0, 0.0)
+ self.vao_octa.render()
+
+ # render fixed back
+ for i_vtx in self.fixback2vtx:
+ pos = self.vtx2xyz_def[i_vtx].copy()
+ model_transf = pyrr.Matrix44.from_translation(pos) * pyrr.Matrix44.from_scale((r, r, r, 1.))
+ self.prog['matrix'].value = tuple((view_transf * model_transf).flatten())
+ self.prog['color'].value = (1.0, 0.0, 0.0)
+ self.vao_octa.render()
+
+ # take a screenshot
+ if not self.is_screenshot_taken and time > 2.2:
+ self.is_screenshot_taken = True
+ rgb = np.frombuffer(self.ctx.fbo.read(), dtype=np.uint8)
+ rgb = rgb.reshape(self.ctx.fbo.size[0], self.ctx.fbo.size[1], 3)
+ if rgb.shape[0] == 1000:
+ rgb = rgb[::2, ::2, :].copy('C')
+ rgb = Image.fromarray(rgb)
+ ImageOps.flip(rgb).save("output.png")
+
+
+def main():
+ HelloWorld.run()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/task09/out_bilaplacian.png b/task09/out_bilaplacian.png
new file mode 100644
index 0000000..6680934
Binary files /dev/null and b/task09/out_bilaplacian.png differ
diff --git a/task09/out_default.png b/task09/out_default.png
new file mode 100644
index 0000000..2bbe647
Binary files /dev/null and b/task09/out_default.png differ
diff --git a/task09/out_laplacian.png b/task09/out_laplacian.png
new file mode 100644
index 0000000..c0eca84
Binary files /dev/null and b/task09/out_laplacian.png differ
diff --git a/task09/preview.png b/task09/preview.png
new file mode 100644
index 0000000..81fa83a
Binary files /dev/null and b/task09/preview.png differ
diff --git a/task09/requirements.txt b/task09/requirements.txt
new file mode 100644
index 0000000..aaa64f3
--- /dev/null
+++ b/task09/requirements.txt
@@ -0,0 +1,9 @@
+glcontext==2.5.0
+moderngl==5.10.0
+moderngl-window==2.4.6
+multipledispatch==1.0.0
+numpy==1.26.4
+pillow==10.3.0
+pyglet==2.0.15
+pyrr==0.10.3
+scipy==1.14.0
diff --git a/task09/util_for_task09.py b/task09/util_for_task09.py
new file mode 100644
index 0000000..866b01e
--- /dev/null
+++ b/task09/util_for_task09.py
@@ -0,0 +1,76 @@
+import numpy
+import numpy.typing
+import math
+
+
+def load_triangle_mesh_from_wavefront_obj_file(file_path: str):
+ try:
+ vtx2xyz = []
+ tri2vtx = []
+ with open(file_path) as f:
+ for line in f:
+ if line[0] == "v":
+ xyz = list(map(float, line[2:].strip().split()))
+ vtx2xyz.append(xyz)
+ elif line[0] == "f":
+ vtx = list(map(int, line[2:].strip().split()))
+ tri2vtx.append(vtx)
+ tri2vtx = numpy.array(tri2vtx, dtype=numpy.uint32) - 1
+ vtx2xyz = numpy.array(vtx2xyz, dtype=numpy.float32)
+ return (tri2vtx, vtx2xyz)
+
+ except FileNotFoundError:
+ print(f"{file_path} not found.")
+ except:
+ print("An error occurred while loading the shape.")
+
+
+def normalize_vtx2xyz(vtx2xyz):
+ """
+ fit the points inside unit cube [0,1]^3
+ :param tri2center:
+ :return:
+ """
+ vmin = numpy.min(vtx2xyz, axis=0)
+ vmax = numpy.max(vtx2xyz, axis=0)
+ vtx2xyz -= (vmin+vmax)*0.5
+ vtx2xyz *= 1.0/(vmax - vmin).max()
+
+
+def vtx2nrm(tri2vtx, vtx2xyz) -> numpy.typing.NDArray:
+ v0 = tri2vtx[:, 0]
+ v1 = tri2vtx[:, 1]
+ v2 = tri2vtx[:, 2]
+ u = vtx2xyz[v1] - vtx2xyz[v0]
+ v = vtx2xyz[v2] - vtx2xyz[v0]
+ c = numpy.cross(u, v)
+ c = c / numpy.linalg.norm(c, axis=1)[:, numpy.newaxis]
+ vtx2nrm = numpy.zeros_like(vtx2xyz)
+ vtx2nrm[v0] += c
+ vtx2nrm[v1] += c
+ vtx2nrm[v2] += c
+ vtx2nrm = vtx2nrm / numpy.linalg.norm(vtx2nrm, axis=1)[:, numpy.newaxis]
+ return vtx2nrm
+
+
+def octahedron() -> (numpy.typing.NDArray, numpy.typing.NDArray):
+ tri2vtx = [
+ [0, 1, 2],
+ [0, 2, 3],
+ [0, 3, 4],
+ [0, 4, 1],
+ [5, 4, 3],
+ [5, 3, 2],
+ [5, 2, 1],
+ [5, 1, 4]]
+ vtx2xyz = [
+ [0., 0., -1.],
+ [1., 0., 0.],
+ [0., 1., 0.],
+ [-1., 0., 0.],
+ [0., -1., 0.],
+ [0., 0., 1.]]
+ tri2vtx = numpy.array(tri2vtx, dtype=numpy.uint32)
+ vtx2xyz = numpy.array(vtx2xyz, dtype=numpy.float32)
+
+ return (tri2vtx, vtx2xyz)