Archive for February, 2016

Simple “parallel for” using C++11 std::thread

Thursday, February 25th, 2016

It seems like Apple’s built-in clang will never support OpenMP. I really used OpenMP only for its #omp parallel for. Now that C++11 has support for multi-threading, I tried to write a very basic parallel for without OpenMP.

#include <vector>
#include <thread>
#include <iostream>

int main(int argc, char *argv[])
{
  const size_t nloop = 11;

  // Serial version
  {
    // Pre loop
    std::cout<<"serial:"<<std::endl;
    // loop over all items
    for(int i = 0;i<nloop;i++)
    {
      // inner loop
      {
        const int j = i*i;
        std::cout<<j<<std::endl;
      }
    }
    // Post loop
    std::cout<<std::endl;
  }

  // Parallel version
  // number of threads
  const size_t nthreads = std::thread::hardware_concurrency();
  {
    // Pre loop
    std::cout<<"parallel ("<<nthreads<<" threads):"<<std::endl;
    std::vector<std::thread> threads(nthreads);
    std::mutex critical;
    for(int t = 0;t<nthreads;t++)
    {
      threads[t] = std::thread(std::bind(
        [&](const int bi, const int ei, const int t)
        {
          // loop over all items
          for(int i = bi;i<ei;i++)
          {
            // inner loop
            {
              const int j = i*i;
              // (optional) make output critical
              std::lock_guard<std::mutex> lock(critical);
              std::cout<<j<<std::endl;
            }
          }
        },t*nloop/nthreads,(t+1)==nthreads?nloop:(t+1)*nloop/nthreads,t));
    }
    std::for_each(threads.begin(),threads.end(),[](std::thread& x){x.join();});
    // Post loop
    std::cout<<std::endl;
  }
}

It’s not as simple as slapping down #omp parallel for but it’s really just a few lines above and below the for loop. It can even determine the number of cores available and handle simple atomic operations.

Repeat each row of matlab matrix

Tuesday, February 16th, 2016

I swear I’ve posted this before, but here’s how to turn:

X = [11 12 13; ...
     21 22 23];

Into:

Y = [
  11 12 13; ...
  11 12 13; ...
  11 12 13; ...
  11 12 13; ...
  21 22 23; ...
  21 22 23; ...
  21 22 23; ...
  21 22 23];

That is, repeat each row, but maintain the order

Y = reshape(permute(repmat(X,[1 1 num_reps]),[3 1 2]),[],size(X,2));

Solid angle at mesh vertices

Wednesday, February 3rd, 2016

The solid angle (aka generalized winding number) at a point lying on a surface can be defined to be the signed area of the surface of the mesh projected onto the unit sphere centered at that point.

For triangle meshes, we can always theoretically compute this as sum of the signed solid angles of each individual triangle. Since the point lies on the mesh, we should be careful about possible numerical pitfalls. In particular, we should ignore any (potentially erroneous) contribution from faces incident on the point. For example, if summing solid angles for the barycenter of some face f we should omit the contribution of f since we know its projection ought to be exactly zero: it projects to a curve on the sphere. Theoretically we only need to ignore contributions from faces intersecting the evaluation point, but one could play it safe and ignore each face for which the evaluation point lies on its plane.

If the mesh is a closed, embedded manifold then the solid angle is a local computation involving the dihedral angles of the incident faces.

For points in the middle of a face f (i.e., not at an edge or vertex), the solid angle is as the face’s plane cuts the sphere in half. Anything projecting on the outside half must double-back on itself to cancel out since the surface is closed and embedded. For points along the middle of the edge between two faces f and g (i.e., not at a vertex) the solid angle is proportional to dihedral angle. The incident faces cut a wedge from the sphere with surface area (assuming Φ runs from 0 to . Finally, for the solid angle at a vertex, we can utilize the spherical excess formula, stating that the solid angle :

      N
Ω = ( ∑ Φi ) - (N-2)π
     i=1

Or using gptoolbox

[AD,C] = adjacency_dihedral_angle_matrix(V,F);
AF = AD~=0;
[AFI,~,AFV] = find(AF);
[~,~,CV] = find(C);
[ADI,~,ADV] = find(AD);
D = sparse( F(sub2ind(size(F),[ADI;ADI],mod([CV;CV+1],3)+1)), 1, 1*repmat(ADV,2,1),size(V,1),1)/2;
N = sparse(F,1,1,size(V,1),1);
SD = D - (N-2)*pi;

Pseudocoloring the solid angle per-vertex has an interesting effect: one can interpret the solid angle as a sort of cheap-o local ambient occlusion. Indeed it’s revealing how much of the immediate surface is occluding any incoming light.

knight solid angle ambient occlusion

Of course it only works at edges and vertices (where there could be any curvature).

Boolean operations using generalized winding numbers

Tuesday, February 2nd, 2016

booleans using generalized winding number

I posted a little tech report about how to use the generalized winding number for constructive-solid geometry (CSG) style boolean operations (union, intersection, difference etc.) on nasty meshes.

Code implementing this has existed for a little while now in the mesh_boolean_winding_number function of gptoolbox.

Code and presentation slides for Nested Cages

Monday, February 1st, 2016

nested cages teaser

We’ve released the source code and presentation slides for our paper Nested Cages, presented last November at SIGGRAPH Asia 2015.