The discrete Laplacian matrix:
$L_{ij} = $ $ \begin{cases} w_{ij} = \frac{1}{2} \left(cot\; \alpha_{ij} + cot \;\beta_{ij}\right) \text{if $j$ is adjacent to $i$}\\ -\sum_{j \in \mathcal{N(i)}} w_{ij} \; \text{if $i=j$}\\ 0 \; \text{otherwise} \end{cases} $
I understand that the laplacian is the divergence of the gradient. I intuitively interpret the laplacian as a function outputting the a change required to smooth the original function in a local neighborhood, but I do not understand why the discretization of the laplacian uses cotangent weights and I do not know how to intuitively interpret them.
You can make sense of that Laplacian matrix in the context of the finite element method (FEM). Imagine you have a function $f_i$ defined on each node $i$. Then, you can interpolate $$ \bar{f}(r) = \sum_i f_i \phi_i(r) , $$ where $\phi_i(r)$ is an interpolating function. You may take tent-shaped functions, in 2D $\phi_i(r)$ looks like a camping tent, with value $1$ at $i$ and a base made up of all triangles incident to $i$ (you need an underlying triangulation, and the Delaunay is a good candidate.)
If you want to obtain the Laplacian of a function, you may want to approximate the Laplacian as $$ \nabla^2 \bar{f}(r) = \sum_i f_i \nabla^2 \phi_i(r) , $$ but these functions cannot be differentiated twice! One way to circumvent this is by noting that you usually want the Laplacian at some node: $$ \nabla^2 f_i = \int dr \phi_i(r) (\nabla^2 \bar{f}(r) ) = \int dr \phi_i(r) \sum_j f_j \nabla^2 \phi_j(r) , $$ where you are now using your interpolating functions to do the opposite task of going from a continuum function to a value defined on nodes only. (Actually, this formula lacks something called "the mass matrix", which in its simplest form means you have to divide by $1/3$ the area of your tent.) Now, the trick is to use integration by parts, and write $$ \nabla^2 f_i = \sum_j f_j \int dr \phi_i(r) \nabla^2 \phi_j(r) = - \sum_j f_j \int dr \nabla \phi_i(r) \cdot \nabla \phi_j(r) . $$ That integrand will only be non-zero in the region on which the shape functions $\phi_i$ and $\phi_j$ overlap, which is the two triangles that lie on either side of the edge connecting them. The integration is easy, since both gradients are just constant vectors. For instance, in one triangle $\phi_i$ will be a vector pointing directly towards vertex $i$ from opposite side (i.e. in the direction of the altitude), and with modulus $1/h_1$, where $h_1$ is the length of the altitude. With this, and taking into account the scalar product, you can simplify that integral and show that $$ \int dr \nabla \phi_i(r) \cdot \nabla \phi_j(r) = -\frac12 \cot\alpha_{ij} . $$
This is a nice set of notes for these issues. Also, notice the true Laplacian still lacks that "mass", see mesh Laplacians on wikipedia.