<jupyter_start><jupyter_text>Assigment 4<jupyter_code>import math
import numpy as np
from numpy import matlib
import scipy.sparse as sp

import igl
import meshplot as mp

from math import sqrt
v, f = igl.read_triangle_mesh("data/irr4-cyl2.off")
tt, _ = igl.triangle_triangle_adjacency(f)
c = np.loadtxt("data/irr4-cyl2.constraints")
cf = c[:, 0].astype(np.int64)
c = c[:, 1:]<jupyter_output><empty_output><jupyter_text>1. Tangent vector fields for scalar field design ![title](img_ipynb/1.png)<jupyter_code>def align_field(V, F, TT, hard_id, hard_value):
    assert(hard_id[0] > 0)
    assert(hard_id.shape[0] == hard_value.shape[0])

    
    # Edges
    e1 = V[F[:, 1], :] - V[F[:, 0], :]
    e2 = V[F[:, 2], :] - V[F[:, 0], :]

    # Compute the local reference systems for each face, T1, T2
    T1 = e1 / np.linalg.norm(e1, axis=1)[:,None]
        
    T2 =  np.cross(T1, np.cross(T1, e2))
    T2 /= np.linalg.norm(T2, axis=1)[:,None]
    
    # Project on the local frame
    c = np.zeros(hard_id.shape[0], dtype=np.complex128)
    
    for ci in range(hard_id.shape[0]):
        f = hard_id[ci]
        v = hard_value[ci, :]
        
        c[ci] = np.dot(v, T1[f, :]) + np.dot(v, T2[f, :])* 1j
  
    # Arrays for the entries of the matrix
    data = []
    ii = []
    jj = []

    index = 0
    for f in range(F.shape[0]):
        for ei in range(3): # Loop over the edges
            
            # Look up the opposite face
            g = TT[f, ei]
            
            # If it is a boundary edge, it does not contribute to the energy
            # or avoid to count every edge twice
            if g == -1 or f > g:
                continue
                
            # Compute the complex representation of the common edge
            e  = V[F[f, (ei+1)%3], :] - V[F[f, ei], :]
            
            vef = np.array([np.dot(e, T1[f, :]), np.dot(e, T2[f, :])])
            vef /= np.linalg.norm(vef)
            ef = (vef[0] + vef[1]*1j).conjugate()
            
            veg = np.array([np.dot(e, T1[g, :]), np.dot(e, T2[g, :])])
            veg /= np.linalg.norm(veg)
            eg = (veg[0] + veg[1]*1j).conjugate()
            
            # Create the matrix Q and b
            data.append(ef);  ii.append(index); jj.append(f)
            data.append(-eg); ii.append(index); jj.append(g)

            index += 1
        
    # Lhs of the system
    A = sp.coo_matrix((data, (ii, jj)), shape=(index, F.shape[0])).asformat("csr") 
    Q = A.H @ A
    Q_prime = Q.copy().asformat("lil")
    Q_prime[:, hard_id] = 0
    Q_prime[hard_id, :] = 0
    Q_prime[hard_id, hard_id] = 1
    # Rhs of the system
    b = np.zeros(F.shape[0], dtype=np.complex128)
    b = -Q[:, hard_id] @ c
    b[hard_id] = c
    # Solve the linear system
    u = sp.linalg.spsolve(Q_prime.asformat("csr"), b)
    R = T1 * u.real[:,None] + T2 * u.imag[:,None]

    return R
def plot_mesh_field(V, F, R, constrain_faces):
    # Highlight in red the constrained faces
    col = np.ones_like(F)
    col[constrain_faces, 1:] = 0
    
    # Scaling of the representative vectors
    avg = igl.avg_edge_length(V, F)/2

    #Plot from face barycenters
    B = igl.barycenter(V, F)

    p = mp.plot(V, F, c=col)
    p.add_lines(B, B + R * avg)
    
    return p
R = align_field(v, f, tt, cf, c)
plot_mesh_field(v, f, R, cf)<jupyter_output>D:\Anaconda\envs\GM\lib\site-packages\jupyter_client\session.py:716: UserWarning: Message serialization failed with:
Out of range float values are not JSON compliant
Supporting this message is deprecated in jupyter-client 7, please make sure your message is JSON-compliant
  content = self.pack(content)<jupyter_text>2. Reconstructing a scalar field from a vector field ![title](img_ipynb/2.png)<jupyter_code>def scalar_field(V, F, u):
    A = sp.diags(np.tile(0.5 * igl.doublearea(V, F), 3))
    G = igl.grad(V, F)
    K = G.T @ A @ G
    b = -2 * G.T @ A @ u
    Lhs = K.T + K
    Lhs[0, :] = 0
    Lhs[:, 0] = 0
    Lhs[0, 0] = 1
    Rhs = -b
    Rhs[0] = 0
    s = sp.linalg.spsolve(Lhs, Rhs)
    g = G @ s
    g = np.reshape(g, f.shape, order='F')
    return s, g
def plot_scalar_field(V, F, s, g):
    # Scaling of the representative vectors
    avg = igl.avg_edge_length(V, F)/2

    #Plot from face barycenters
    B = igl.barycenter(V, F)

    p = mp.plot(V, F, c=s)
    p.add_lines(B, B + g * avg)
    
    return p
def plot_error_field(V, F, R, g):
    # Scaling of the representative vectors
    avg = igl.avg_edge_length(V, F)/2

    #Plot from face barycenters
    B = igl.barycenter(V, F)

    p = mp.plot(V, F, c=np.linalg.norm(R - g, ord=2, axis=1))

    print("Poisson Error: ", sum(np.linalg.norm(R - g, ord=2, axis=1)))
    
    return p
s, g = scalar_field(v, f, R.T.flatten())
plot_scalar_field(v, f, s, g)
plot_error_field(v, f, R, g)<jupyter_output><empty_output><jupyter_text>3. Harmonic and LSCM Parameterizations<jupyter_code>v_camel, f_camel = igl.read_triangle_mesh("data/camel_head.off")
def harmonic(V, F):
    # Find the open boundary
    bnd = igl.boundary_loop(F)

    # Map the boundary to a circle, preserving edge proportions
    bnd_uv = igl.map_vertices_to_circle(V, bnd)

    # Harmonic parametrization for the internal vertices
    uv = igl.harmonic_weights(V, F, bnd, bnd_uv, 1)
    
    g = np.reshape(igl.grad(V, F) @ uv[:, 1], F.shape, order='F')
    return g, uv

def LSCM(V, F):
    # Fix two points on the boundary
    b = np.array([2, 1])

    bnd = igl.boundary_loop(F)
    b[0] = bnd[0]
    b[1] = bnd[int(bnd.size / 2)]

    bc = np.array([[0.0, 0.0], [1.0, 0.0]])

    # LSCM parametrization
    _, uv = igl.lscm(V, F, b, bc)

    return uv
def plot_parameterizations(V, F, uv_harmonic, uv_LSCM, g_harmonic):
    v_p = np.hstack([uv_harmonic, np.zeros((uv_harmonic.shape[0], 1))])
    p = mp.subplot(V, F, s=[1, 4, 0], uv=uv_harmonic)
    mp.subplot(v_p, F, shading={"wireframe": True}, s=[1, 4, 1], data=p)

    avg = igl.avg_edge_length(V, F) / 2
    B = igl.barycenter(V, F)
    pp = mp.plot(V, F, uv=uv_harmonic, return_plot=True)
    pp.add_lines(B, B + g_harmonic * avg)
    
    mp.subplot(V, F, s=[1, 4, 2], uv=uv_LSCM, data=p)
    mp.subplot(uv_LSCM, F, shading={"wireframe": True}, s=[1, 4, 3], data=p)


g_harmonic, uv_harmonic = harmonic(v_camel, f_camel)
uv_LSCM = LSCM(v_camel, f_camel)
plot_parameterizations(v_camel, f_camel, uv_harmonic, uv_LSCM, g_harmonic)<jupyter_output><empty_output><jupyter_text>4. Editing a parameterization with vector fields  Editing the parameterization<jupyter_code>g_irrcyl, uv_irrcyl = harmonic(v, f)
p = mp.subplot(v, f, uv=uv_irrcyl, s=[1,2,0])
uv_edited = uv_irrcyl.copy()
uv_edited[:, 1] = s
mp.subplot(uv_edited, f, s=[1,2,1], shading={"wireframe": True}, data=p)

p = mp.plot(v, f, uv=uv_edited)
avg = igl.avg_edge_length(v, f) / 2
B = igl.barycenter(v, f)
p.add_lines(B, B + g * avg, shading={"line_color": "red"})
p.add_lines(B, B + np.reshape(igl.grad(v, f) @ uv_irrcyl[:, 0], f.shape, order='F') * avg, shading={"line_color": "blue"})
# p.add_lines(B, B + np.reshape(igl.grad(v, f) @ uv_irrcyl[:, 1], f.shape, order='F') * avg, shading={"line_color": "green"})<jupyter_output><empty_output><jupyter_text>Detecting problems with the parameterization<jupyter_code>v_p_edited = np.hstack([uv_edited, np.zeros((uv_edited.shape[0],1))])
normal = igl.per_face_normals(v, f, np.array([0.0, 0.0, 0.0]))
normal_edited = igl.per_face_normals(v_p_edited, f, np.array([0.0, 0.0, 0.0]))

# normal = np.zeros(f.shape)
# normal_edited = np.zeros(f.shape)
# for fi in range(f.shape[0]):
#     e1  = v[f[fi, 1], :] - v[f[fi, 0], :]
#     e2 = v[f[fi, 2], :] - v[f[fi, 0], :]
#     normal[fi, :] = np.cross(e1, e2)
#     # normal[fi, :] /= np.linalg.norm(normal[fi, :])


#     e3  = v_p_edited[f[fi, 1], :] - v_p_edited[f[fi, 0], :]
#     e4 = v_p_edited[f[fi, 2], :] - v_p_edited[f[fi, 0], :]
#     normal_edited[fi, :] = np.cross(e3, e4)
#     # normal_edited[fi, :] /= np.linalg.norm(normal_edited[fi, :])

col = np.ones_like(f)
for i in range(f.shape[0]):
    if np.dot(normal[i, :], normal_edited[i, :]) < 0:
        col[i, 1:] = 0

p = mp.plot(v_p_edited, f, shading={"wireframe": True}, c=col)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16875153298994358, 'fraction_numerical': 0.02146185921020358, 'mean_word_length': 2.903781713738631, 'pattern_counts': {'":': 7, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '218643', 'license': 'MIT', 'n_tokens_mistral': 3628, 'n_tokens_neox': 3432, 'n_words': 774, 'path': 'Assignment_4/Assigment4.ipynb', 'repo_name': 'AmazingJinqi/gp', 'size': '8154'}