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. 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. 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. 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. 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. 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. 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. 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. 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. 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. 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}$?

  3. 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. 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. 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. 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. 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. 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. 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. $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. 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. 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. 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…

  5. 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. 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. 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. 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

    2. 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. 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. 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.

          2. 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.

          3. 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$.

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

  7. 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. 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. 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. 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.

    2. 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.

Leave a Reply