# Assignment 4 (Coding): Conformal Parameterization

For the coding portion of your assignment on conformal parameterization, you will implement the Spectral Conformal Parameterization (SCP) algorithm as described in the course notes.

Please implement the following routines in

1. core/geometry.js:
1. complexLaplaceMatrix
2. projects/parameterization/spectral-conformal-parameterization.js:
1. buildConformalEnergy
2. flatten
3. utils/solvers.js:
1. solveInversePowerMethod
2. residual

Notes

• The complex version of the cotan-Laplace matrix can be built in exactly the same manner as its real counterpart. The only difference now is that the cotan values of our matrix will be complex numbers with a zero imaginary component.
• The buildConformalEnergy function builds a $|V| \times |V|$ complex matrix corresponding to the conformal energy $E_C(z) = E_D(z) – A(z)$, where $E_D$ is the Dirichlet energy (given by our complex cotan-Laplace matrix) and $A$ is the total signed area of the region $z(M)$ in the complex plane (Please refer to Section 7.4 of the notes for more details). You may find it easiest to iterate over the halfedges of the mesh boundaries to construct the area matrix.
• The flatten function returns a dictionary mapping each vertex to a vector (linear-algebra/vector.js) of planar coordinates by finding the eigenvector corresponding to the smallest eigenvalue of the conformal energy matrix. You can compute this eigenvector by using solveInversePowerMethod (described below).
• Your solveInversePowerMethod function should implement Algorithm 1 in Section 7.5 of the course notes with one modification – after computing $Ay_i = y_{i-1}$, center $y_i$ around the origin by subtracting its mean. Important: Terminate your iterations when your residual is smaller that 10^-10.
• The parameterization project directory also contains a basic implementation of the Boundary First Flattening (BFF) algorithm. Feel free to play around with it in the viewer and compare the results to your SCP implementation.

Submission Instructions

Please rename your geometry.js, spectral-conformal-parameterization.js and solvers.js files to geometry.txt, spectral-conformal-parameterization.txt and solvers.txt(respectively) and put them in a single zip file called solution.zip. These files and your solution to the written exercises should be submitted together in a single email to Geometry.Collective@gmail.com with the subject line DDG17A4.

## 49 thoughts on “Assignment 4 (Coding): Conformal Parameterization”

1. yangxio says:

When I run buildConformalEnergy, I see my memory usage soaring to several gigabytes and obviously cannot get any result before Chrome halts the evaluation. I am pretty sure that it relates to complexLaplaceMatrix instead of my implementation of the complex area since when I take out the complexLaplaceMatrix I can at least get something out of my browser. I also try calling complexLaplaceMatrix in laplaceMatrix and use the test script for geometry.js but get the message that ComplexTriplet is not found. I am wondering what is going on here technically. (And my laptop is a beast with a memory of 16GB…)
In particular, would you please add a test for complexLaplaceMatrix in the test script for geometry? I understand it is unnecessary for numerics, as it is exactly the same as the real one. But apparently there seems to be a technical problem here.

1. Keenan says:

The most common error is that you are filling $O(n^2)$ entries of a matrix, e.g., you are visiting every single entry of an $n \times n$ matrix and setting its value, rather than just setting the values of the nonzeros. For instance, if you have any kind of double loop where each loop ranges over all the vertices (or edges, or faces), then you are likely doing this.

1. yangxio says:

Well, it becomes quite a delicate issue now: there is evidence supporting your diagnostics but it is definitely not convincing.

I simply adjust my laplaceMatrix code from Assignment 3, which is of order $O(nd)$ where $d$ is the (maximum) degree of vertices in the mesh, which is pretty much $O(n)$ if the mesh is regular enough. Moreover, I do similar thing with the solution to A3 posted on GitHub and still run into the same situation. If the meshes for A3 and A4 are roughly of the same magnitude, then I should have run into the problem in A3 weeks ago. In short, I am sure I am not applying quadratic operations even if mine is superlinear.

Another test I do is to comment out all the addEntry’s in my loop and run the test script. The memory activity stays calm but the code is still running slow. I agree that I haven’t fully exploited the symmetry of the laplacian, but this can only save half of the operations which is of the same order. Also, I simply run the index.html in the folder projects/parametrization. It takes a long time to load and the BFF gives me nothing but a complete darkness background or foreground. It also consumes quite some memory but not as crazy as the test script. Should the BFF run anyway even if SCF is not or ill implemented?

Just some timing from the test scripts for my previous homework here for further diagnostics:
laplaceMatrix: 216ms in Chrome, 221ms in Edge
massMatrix: 82ms in Chrome, 137ms in Edge
MeanCurvatureFlow.integrate: 1247ms in Chrome, 101ms in Edge
ModifiedMCF.integrate: 689ms in Chrome, 628ms in Edge
I understand these time vary from time to time, but it should give a rough idea. Also some specs here: i7-7500U, 16GB RAM, windows 10 home edition

1. rohansawhney says:

It seems like you are building the laplace matrix as expected – the only difference compared to the real case should be that the entries of the Laplace matrix for this assignment are complex.

The parameterization module will take longer to load compared to the rest. This is an artifact of the linear algebra module being used – it was compiled from C++ code using emscripten. This module sets up its own heap (every time you launch index.html or test.html) and does not use Javascript’s memory management system. The increase in memory consumption you are noticing comes from allocating complex objects (that belong to the linear algebra module) on this heap. However, this memory is cleared every time you flatten the mesh.

I am able to run BFF on my computer using the latest codebase on github that does not have SCP implemented. The memory consumption is large on my machine, but nothing Chrome, Safari and Firefox can’t handle (I’m using a MacBook Pro, which is far from beefy). I recommend pulling/cloning the latest version of the repository. If this still doesn’t resolve the issue, then maybe we should meet in person (I will be in Smith from Wednesday to Friday) – the problem might be platform specific.

1. yangxio says:

Thanks a lot, this is helpful. I haven’t pulled from the repository since your “minor” update on October 2nd, and clearly it becomes quite a major issue now. Now the index.html loads instantly fast and I can pass the buildConformalEnergy test now.

However, I do have two more questions regarding the conformal energy matrix and the spectral method.

first a minor one, usually the matrix $A$ encoding the quadratic form $E$ is not halved, i.e. $E(z) = \frac{1}{2}\langle Az, z\rangle$, but here it seems that the test script is looking for $\frac{1}{2}A$ instead of $A$ or otherwise I could not pass the test?

second, according to the notes and the link, one should incorporate the geometry of the surface and hence consider the mass matrix and accordingly the generalized eigenvalue problem with the generalized (centered) inverse power method, which could not be reduced to the usual one and apply the naïve “normalize” function as in geometry.js. So I am skeptical of the hint mentioned above. Also the massMatrix we built in A3 is really of the real sparse matrix class and we cannot manipulate it with the complex counterpart without a bottom-up construction given the current library. Is there something I am missing here?

2. rohansawhney says:

Yes, the reference solution accounts for 1/2 factor, even though it does not affect the final result – just being consistent with the notes.

You should subtract the “simple” mean, not the area weighted mean. For this assignment, we require you to implement the inverse power method given by Algorithm 1 in the notes (+ modification that subtracts the mean from the solution), not the generalized power method given by Algorithm 2.

2. Leslie says:

Hi, in “boundary-first-flattening”, the function “invert2x2” in class “Solvers” is called by “let m = Solvers.invert2x2(T.timesDense(Ninv.timesDense(TT)));”.

I have been trying to call the function “solveInversePowerMethod(A)” in class Solvers by this same command:
“let x = Solvers.solveInversePowerMethod(A)”, but it doesn’t work. The error returned is “Solvers is not defined”.

Besides, I also tried to play around with the given BFF algorithm but got nothing in the viewer.

Do we need to change anything to make it work or is there anything wrong with my code?

Thanks!

1. rohansawhney says:

Try pulling the latest version of the codebase from GitHub. You should be able to play around with BFF without implementing SCP.

1. Leslie says:

Thanks! I have solved the “Solvers is not defined” problem and I can get the final results. But then I am in a kind of tricky situation.
If I set $\epsilon = 1e^{-10}$ or $1e^{-12}$ and run the test, I get the results which are pretty closed to the “solution” file but the norm of difference is larger then the pre-set threshold $1e^{-6}$.
If I loose the threshold to $1e^{-4}$, all the results are correct but I will see the error because of Timeout of 2000ms exceeded.
If I set $\epsilon < 1e^{-12}$, I will get results quite different from the solution. My initial guess is the largest 2 eigenvalues of the matrix we deal with are pretty closed, so we will get something like "local minima" if $\epsilon$ is not sufficiently small.
Do you think is there anything wrong with the whole procedure or the code?
Thanks!

1. rohansawhney says:

You should set epsilon to 10^-10 because the reference solution was computed with that value. You probably have a missing factor of 1/2 or a sign error in your code. Make sure your implementation matches the formulas in the notes and takes the clockwise orientation of the boundary halfedges into account.

1. David Zeng says:

Since the residual is $Ax – (x^TAx)x$, a vector rather than a number, is there a specific norm we should use when checking if the residual is smaller than $10^{-10}$?

2. rohansawhney says:

Yes, use the L2 norm.

3. anobani says:

i am a little confused and lost about the co formal energy function, and i don’t know how to proceed.
so we need to solve $\dfrac{1}{2} \Delta – A$, such as $A = \dfrac{-i}{4}\sum_{e_{ij} \in E} \bar{z}_i z_j – \bar{z}_j z_i$, since we already have the complex Laplace matrix, we build A by iterating through the edges of the boundary, and we can extract the two vertices from each edge in the boundary.

the question is i can’t understand how to apply the $\bar{z}$ and $z$, should we compute them somehow ? i know this is trivial, but feel stuck translating it into code

1. anobani says:

since the conformal map preserve the angle, can we relate to the original mesh and find the difference between the angles there that reflects on the conformal map, the same as using $\bar{z}$ and $z$ ?

1. rohansawhney says:

The buildConformalEnergy function requires you to return the matrix $E_C = \dfrac{1}{2} \Delta – A$ such that $\bar{z}^T E_C z = \dfrac{1}{2}\bar{z}^T \Delta z + \dfrac{i}{4}\sum_{e_{ij} \in E} \bar{z}_i z_j – \bar{z}_j z_i$, where $z$ is a |V| x 1 complex vector containing the coordinates of your flattening. This step does not require you to do anything with $z$, you are just supposed to construct the matrix $E_C$. Next, you actually compute $z$ by solving the optimization problem $min_z \bar{z}^T E_C z$, subject to $||z|| = 1$ and $\sum_i z_i = 0$.This can be done using the inverse power method procedure described in Section 7.5 of the notes.

Hope this clears things up, but feel free to ask follow up questions.

1. rohansawhney says:

To construct the matrix $A$, start by considering the simple case of a triangle with some complex coordinates $z_1, z_2, z_3$. Then, $z$ is a 3 x 1 column vector $[z_1; z_2; z_3]$. Can you come with the 3 x 3 matrix A such that $\bar{z}^T A z = -\frac{i}{4}(\bar{z_1}z_2 – \bar{z_2}{z_1} + \bar{z_2}z_3 – \bar{z_3}{z_2} + \bar{z_3}z_1 – \bar{z_1}{z_3})$?

1. anobani says:

yes thank you, i started first with the power method and back to the conformal map,
but still i think there is something i don’t quit understand,
through the buildConformalEnergy i call flatten to create the $z$ in which i can iterate, from flatten i need to call solveInversePowerMethod to find my Eigenvectors,, what should i understand from the input of solveInversePowerMethod matrix, can i feed it the complex Laplace matrix? and from that i would have only one Eigenvectors,, should i build another matrix $A$, to find another Eigenvectors,?

2. anobani says:

Thank you, i think i am on track of reasoning over it

3. rohansawhney says:

It seems like you have the order of operations wrong. You should first construct the conformal energy matrix $E_C$, and then compute its eigenvector $z$ (corresponding to the smallest eigenvalue) with the inverse power method. This $z$ is the solution to the optimization problem $min_z \bar{z}^T E_C z$, subject to $||z|| = 1$ and $\sum_i z_i = 0$ and represents the coordinates of your conformal mapping.

2. kye says:

I’m confused about why $E_c$ is applied as $\overline{z}^T E_c z$ and not $E_c z$. From reading the code for the Laplacian, it seems that it was being applied as $L z$. Where in the notes does it specify how the matrices are applied to $\overline{z}$ and $z$?

1. rohansawhney says:

$E_C(z)$ is a quadratic form given by $E_{C}(z) = \frac{1}{2}\langle \langle \Delta z, z \rangle \rangle – A(z)$ and we want to find a $z$ that minimizes $E_C(z)$. Another way to look at this is that we want to minimize an energy (which is a number) $\bar{z}^T E_C z$ where $E_C$ is a matrix given by $\frac{1}{2}\Delta – A$.

4. Yifanh says:

My “solveInversePowerMethod” could only converge to residual ~ 1e-4, then no matter how many iterations I run, the residual does not go down.

At each iteration, I solve $Ay_i=y_{i-1}$ using cholesky decomposition, subtract $y_i$ by its mean (y.sum() is then close to 0, something like 1e-17), then divide it by its 2-norm (y.norm(2) is 1).
The residual of first 10 iterations look like this: (dy is norm(y_new – y_old) )

iter 1 dy: 67553560841.77606 residual: 0.006827703640345173
iter 2 dy: 1.526547749542435 residual: 0.00010232609496920101
iter 3 dy: 0.03727808862025678 residual: 0.0000988644103631873
iter 4 dy: 0.02497516894581189 residual: 0.00009578056223184056
iter 5 dy: 0.024695055750158984 residual: 0.00009938449503269249
iter 6 dy: 0.024678821492846344 residual: 0.00009935638375974477
iter 7 dy: 0.024677444145042995 residual: 0.00009910923673666039
iter 8 dy: 0.02467731750973144 residual: 0.00009724250659028833
iter 9 dy: 0.024677305698380773 residual: 0.00009637130943827171
iter 10 dy: 0.024677304593815855 residual: 0.000098766343763732

then residual is always ~1e-4 for the rest 1000 iterations.

If I do not subtract y by its mean, the algorithm converge within ~10 iterations, but the result won’t pass the test.

Any clue on what could I do wrong?

1. Leslie says:

I had THE same problem previously. I think we need to call the “residual” function in each iteration. Maybe you can check whether the norm of $y$ changes or not after each time you call the “residual” function.
Actually, the norm of $y$ should not change (will be 1), but I used something like “scaleBy” in “residual” part and changed the norm, which made the algorithm not converge.
That’s my experience and you can check the code based on this. But I am not sure it will work or not.

1. Yifanh says:

Oh…you are right…thank you so much for sharing your experience.
I fixed my problem by getting rid of scaleBy().
This is so beyond my imagination..

For those who are new to javascript, here is a list of behaviors you need to be aware:
1. Numbers are assigned by value, but objects are assigned by reference.
let a = 1;
let b = a;
b = 2;
// now, a is still 1. this is sane.
let a = ComplexDenseMatrix.ones(1, 1);
let b = a;
b.scaleBy(new Complex(2,0));
console.log(‘a ‘, a.norm(2));
// now, you will get 2 ……

2. function arguments: number variables are passed by value, objects are passed by reference.
let a=1;
myfunc(a);
// now, a is still 1

let a = ComplexDenseMatrix.ones(1, 1);
myfunc(a);
// now, a may be different…

2. rohansawhney says:

This seemed to be a common problem in office hours. Keep in mind that your solution vector $y$ is complex when computing $\lambda$.

5. yyyy says:

When trying to run my flatten, I get this error:

RangeError: Source is too large
at Uint8Array.set (native)
at _emscripten_memcpy_big

that throws from multiple different locations in different runs, (ComplexDenseMatrix.transpose, ComplexTriplet.addEntry). I’m not sure exactly what this error means, am I running out of memory?

1. rohansawhney says:

This might be happening because you are allocating a lot of complex objects on the heap (see my comment above in response to yangxio’s questions). Maybe try to cut down on how many complex objects you allocate…

6. plav says:

I still need to check why I get into the infinite loop…

When I limit the loop for 100 iterations, my flattening converges to the solution in the test. I printed all the vectors from my flattening and from the solution. They have almost the same value, but the resulting format of the console.log is different:

mine -0.22430933772835607 -0.9222229495889633 0
solution Vector {x: -0.2243093978040462, y: -0.9222230128593861, z: 0}

However, I get the following message from the test:
“Error: Timeout of 2000ms exceeded. For async tests and hooks, ensure “done()” is called; if returning a Promise, ensure it resolves.”

1. minkyuk says:

Out of curiosity, are you centering your solution around the origin in the iteration? I fixed a similar issue by simply centering my solution vector around origin every iteration

1. plav says:

I am
1) subtracting it by a vector with the mean value
2) I am scaling the new y by 1/norm.

2. rohansawhney says:

My implementation of solveInversePowerMethod converges in about 7 iterations. Something is probably not quite right, print your residual values to see how they behave.

1. plav says:

I a stopped using the while loop, because I am getting into a n infinite loop.
I am forcing 100 iterations with a for loop. It kind of converges in the 3rd iteration:
solvers.js:45 residual 2.560910196146597
solvers.js:45 residual 0.006821166870907288
solvers.js:45 residual 0.00009702186648588244
solvers.js:45 residual 0.00009590016341776931
solvers.js:45 residual 0.00009590630261367495
solvers.js:45 residual 0.00009590743607528293
solvers.js:45 residual 0.00009590754789132559
solvers.js:45 residual 0.00009590755847338477
solvers.js:45 residual 0.00009590755946688405
solvers.js:45 residual 0.00009590755955992604
solvers.js:45 residual 0.0000959075595686597
solvers.js:45 residual 0.00009590755956946545
solvers.js:45 residual 0.00009590755956953372
solvers.js:45 residual 0.00009590755956953889
solvers.js:45 residual 0.0000959075595695423
solvers.js:45 residual 0.00009590755956953802
solvers.js:45 residual 0.00009590755956953809
solvers.js:45 residual 0.0000959075595695447
solvers.js:45 residual 0.00009590755956952814
solvers.js:45 residual 0.00009590755956953502
solvers.js:45 residual 0.00009590755956954513
solvers.js:45 residual 0.00009590755956953528
solvers.js:45 residual 0.00009590755956953042
solvers.js:45 residual 0.00009590755956954738
solvers.js:45 residual 0.00009590755956954283
solvers.js:45 residual 0.00009590755956953515
solvers.js:45 residual 0.0000959075595695483
solvers.js:45 residual 0.00009590755956953749
solvers.js:45 residual 0.00009590755956953687
solvers.js:45 residual 0.00009590755956955362
solvers.js:45 residual 0.0000959075595695343
solvers.js:45 residual 0.00009590755956954393
solvers.js:45 residual 0.00009590755956954005
solvers.js:45 residual 0.00009590755956953478
solvers.js:45 residual 0.00009590755956954202
solvers.js:45 residual 0.00009590755956953628
solvers.js:45 residual 0.00009590755956954302
solvers.js:45 residual 0.00009590755956953753
solvers.js:45 residual 0.00009590755956954419
solvers.js:45 residual 0.00009590755956953382
solvers.js:45 residual 0.00009590755956955037
solvers.js:45 residual 0.00009590755956953764
solvers.js:45 residual 0.00009590755956954192
solvers.js:45 residual 0.0000959075595695407
solvers.js:45 residual 0.00009590755956954168
solvers.js:45 residual 0.00009590755956954009
solvers.js:45 residual 0.00009590755956954328
solvers.js:45 residual 0.00009590755956953284
solvers.js:45 residual 0.00009590755956953989
solvers.js:45 residual 0.00009590755956954354
solvers.js:45 residual 0.00009590755956953389
solvers.js:45 residual 0.00009590755956954551
solvers.js:45 residual 0.00009590755956954095
solvers.js:45 residual 0.00009590755956955083
solvers.js:45 residual 0.00009590755956953806
solvers.js:45 residual 0.00009590755956953737
solvers.js:45 residual 0.00009590755956953305
solvers.js:45 residual 0.00009590755956953562
solvers.js:45 residual 0.0000959075595695452
solvers.js:45 residual 0.00009590755956953837
solvers.js:45 residual 0.00009590755956952662
solvers.js:45 residual 0.00009590755956953647
solvers.js:45 residual 0.00009590755956953727
solvers.js:45 residual 0.00009590755956953402
solvers.js:45 residual 0.0000959075595695379
solvers.js:45 residual 0.00009590755956953765
solvers.js:45 residual 0.00009590755956954824
solvers.js:45 residual 0.00009590755956954229
solvers.js:45 residual 0.00009590755956954371
solvers.js:45 residual 0.00009590755956953989
solvers.js:45 residual 0.00009590755956953321
solvers.js:45 residual 0.00009590755956954081
solvers.js:45 residual 0.00009590755956954302
solvers.js:45 residual 0.00009590755956953887
solvers.js:45 residual 0.00009590755956955168
solvers.js:45 residual 0.00009590755956953622
solvers.js:45 residual 0.00009590755956953258
solvers.js:45 residual 0.00009590755956953111
solvers.js:45 residual 0.00009590755956953924
solvers.js:45 residual 0.0000959075595695386
solvers.js:45 residual 0.00009590755956954988
solvers.js:45 residual 0.00009590755956954867
solvers.js:45 residual 0.00009590755956954951
solvers.js:45 residual 0.00009590755956954486
solvers.js:45 residual 0.0000959075595695398
solvers.js:45 residual 0.00009590755956954957
solvers.js:45 residual 0.00009590755956953872
solvers.js:45 residual 0.00009590755956953447
solvers.js:45 residual 0.00009590755956954455
solvers.js:45 residual 0.00009590755956954107
solvers.js:45 residual 0.00009590755956953473
solvers.js:45 residual 0.00009590755956954551
solvers.js:45 residual 0.0000959075595695383
solvers.js:45 residual 0.00009590755956954637
solvers.js:45 residual 0.00009590755956954558
solvers.js:45 residual 0.00009590755956954026
solvers.js:45 residual 0.00009590755956955163
solvers.js:45 residual 0.00009590755956953965
solvers.js:45 residual 0.0000959075595695458
solvers.js:45 residual 0.00009590755956954463

1. rohansawhney says:

Its not converged to the required precision – check if you have the same bug Yifanh had (see his comments above).

1. plav says:

Solved. I was subtracting the mean twice….
Thanks.

2. plav says:

Actually, it only passes because I am using a for loop. The residual fails. Anyway….

3. rohansawhney says:

Subtracting the mean a second time shouldn’t cause any problem since the mean should be zero.

4. plav says:

The fact is that my solveInversePowerMethod works now…but only when I do not use the while loop with the residual. I will double check the algebra now.

5. plav says:

I am only using non-destructive methods to do the algebra (timesDense, minus, timesComplex, transpose, get…). That is why I think my problem is distinct than Yifanh ‘s.

6. anobani says:

i have the same problem, I am using non-destructive methods, however stuck in infinite loop.

7. rohansawhney says:

I don’t know if you are already doing this, but keep in mind that your solution vector is complex, you have to perform an additional operation before/after transposing it while computing $\lambda$.

8. anobani says:

shouldn’t we evaluate $result =Ax – (x^tAx)x$, then return $result.norm()$ for comparison ?

9. rohansawhney says:

Since you working with a complex numbers, $x$ in your expression above should be conjugate transposed and not just transposed.

7. intrepidowl says:

What are the dependencies between each part here?

It seems like for buildConformalEnergy, we need the conformally mapped geometry z(M) to get the $z_i,z_j$. Is this what flatten is supposed to do? But flatten depends on other things, etc.

1. kye says:

I believe buildConformalEnergy does not depend on flatten. It builds a matrix that “waits for” z.

2. rohansawhney says:

You are supposed to build the conformal energy matrix and use the inverse power method on it to compute its smallest eigenvector $z$. Both buildConformalEnergy and solveInversePowerMethod should be called from flatten (in that order). The comments above clarify this.

8. maliang says:

I’m getting the “Solvers is not defined” error even after redownloading the GitHub repository. Trying to play with BFF just gives me a black screen. Is there a way to fix this?

1. Leslie says:

I solved same problem by checking the “residual(A, x)” function. Maybe this function does not return the correct answer and you can check it step by step.

2. rohansawhney says:

Comment out any code you’ve written yourself in the repo you downloaded and open the viewer. You should be able to interact with the BFF implementation then.

3. This error comes whenever your code breaks.. maybe this issue, (xTAx) in residual block is supposed to be a scalar! Also, while multiplication it gets returned as DenseMatrix of 1×1. Also, just check if any brackets aren’t getting missed.