Archive for February, 2010

The trials and tribulations of building a simple mesh viewer

Friday, February 26th, 2010

I am beginning a new project and to start I needed to install a peer’s CGAL-based mesh viewer. The mesh viewer has many dependencies which in turn have even more dependencies. I will try to recap my struggle and eventual success, here. My warned that this is not an installation guide and is probably missing many things I did and maybe even advertising bad solutions.

64-bit on Mac OS X 10.5

Short story: Don’t bother.

Long story: My peer has his whole setup all running on 64 bits. This requires building all the dependencies running on 64 bits, which for Mac OS X 10.4 (Leopard) this means compiling everything from source and being really careful about flags. It also means, as far as I know, not using macports.
I tried to mimic his compilations and installs but without knowing the exact flags I got almost to the point of having everything working but then my 32-bit installs from macports started really getting in the way. I ended up uninstalling all my macports in vain only to realize that my python(s) were not 64-bit or at least couldn’t find 32-bit framework libraries. I recompiled python and pyqt but only found new errors. In the end I gave up on 64-bit (which I’d spent 3 days on) and switched (back) to 32-bit (which only took 1 day).

Dependencies for CGAL using macports

CGAL itself has a macport but it’s poorly maintained and will often not compile correctly with your current setup. It seems that a lot of people are having trouble getting this to even build correctly much less work correctly. Instead use macports to install the dependencies then build cgal from source by hand (not so bad).

CGAL depends on or supports the following libraries that you don’t already have on your mac:

  • boost
  • Gmp
  • Mpfr
  • zlib
  • BLAS*
  • LAPACK*
  • taucs*

Qt3, Qt4 and libqglviewer are only necessary for building the examples and demos which require them and not necessary for the install. Really the only dependency is boost. The rest add more features: exact arithmetic, linear algebra solvers, etc. Those with a *, I’ll write more about later on.

To install most of the above with macports just issue:


sudo port install boost gmp mpfr zlib

Now, you should be ready to build CGAL.

Installing CGAL

cd to the CGAL directory and issue

cmake -i .

You’ll see a prompt for advanced options, I was paranoid so I typed Yes:


Would you like to see advanced options? [No]:Yes

You need to change all library and include directories to point to the macports install directory: /opt/local/lib/WHATEVER or /opt/local/include.

For example:


Variable Name: Boost_INCLUDE_DIR
Description: Path to a file.
Current Value: /usr/local/include
New Value (Enter to keep current value): /opt/local/include

Do this for boost, gmp, mpfr etc. Watch out for the CMAKE_INSTALL_PREFIX prompt. If you want CGAL to hang out with your macports then change this to /opt/local like so:


Variable Name: CMAKE_INSTALL_PREFIX
Description: Install path prefix, prepended onto install directories.
Current Value: /usr/local
New Value (Enter to keep current value): /opt/local

In the end, my hackish solution to a later problem prevents all this /opt/local business from mattering but it seems like good practice to have related software in the same place. Unless you feel like macports should be the only one touching /opt/local/ stuff, but it's your computer so why not.

Then you can make, sudo make install to finish with CGAL.

Python 2.6 and PyQt4

Simply issue:


sudo port install python26 py26-pyqt4 py26-opengl
sudo python_select python26

These should install fine and work on a simple example. Getting CGAL to recognize your python (and subsequently PyQt) is another problem.

taucs with LU

I was trying to install and unsupported version of the math software taucs. This version had LU decomposition necessary for a feature of the mesh viewer. I ran into trouble forcing the taucs build to see my blas and lapack and to make a proper 32 bit file.

After running configure if I tried to run make I would see an error:


build/darwin9.0/makefile:14: config/darwin9.0.mk: No such file or directory
make: *** No rule to make target `config/darwin9.0.mk'.  Stop.

This is because the taucs only cam with a premade make file for darwin, not darwin9.0 as configure has recognized my OS.

For the most part you can just copy the darwin make file:


cp config/darwin.mk config/darwin9.0.mk

But make a few changes, namely edit the following lines:


CFLAGS    = -arch i386 -O3 -faltivec

LIBBLAS   = -framework Accelerate
LIBLAPACK = -framework Accelerate

One more thing. The LIBF77 = -Lexternal/lib/darwin -lf2c line is pointing to the right place but the file there is wrong. If you ar -x external/lib/darwin/libf2c.a you'll find out it's full of x86_64 .o files, which will lead to a confused architecture build. My solution was to download the taucs_full from the CGAL download page and copy the libf2c.a on top of this file. Maybe the same is true for the lbmetis.a, I don't remember. Just check that after you
make your libtaus.a unarchives to i386 mach-o files.

After running make you should move your taucs to the /opt/local area. I just did


sudo mkdir /opt/local/taucs_with_lu
sudo cp -r taucs_with_lu/* /opt/local/taucs_with_lu/.

Hacks

At this point things seemed to work but our meshviewer uses swig to combine cpp code and python. swig was easy enough to install:


sudo port install swig swig-python

But CGAL was not playing nicely. When I built our CGAL meshviewer, CGAL's cmake finders were not locating the libraries I had installed with macports. It kept trying to look in /usr/local/. My hack was to move /usr/local to /usr/local-off and simple link it to /opt/local. Everyone I told this too agreed it was ugly.


sudo mv /usr/local /usr/local-off
sudo ln -s /opt/local /usr/local

This is basically saying, I agree to use macports for everything or I must be very careful.

At this point CGAL would play along and build the toy version of our meshviewer. But python would not display it. I think it was because swig would not read the headers from the right place (/opt/local/WHATEVER) instead it was looking in


/Library/Frameworks/Python.Framework/WHATEVER. To handle this I also hacked:
cd /Library/Frameworks/Python.framework/Versions
ln -s /opt/local/Library/Frameworks/Python.framework/Versions/Current Current
cd /Library/Frameworks/Python.framework/
mv Headers Headers-off
ln -s /opt/local/Library/Frameworks/Python.framework/Headers Headers

Again, everyone I told this too said it was a poor man's hack.

CGAL and a deprecated header

Finally, I had everything in order to build the full version of my meshviewer. But I got funny errors about a certain boost header no existing, namely /opt/local/include/boost/property_map.hpp.

Upon inspection I noticed that indeed this file does not exist. On my peer's boost install (he did not use macports) he had this file but opening it we saw that it had been long deprecated by boost and simple pointed to the real header in property_map/property_map.hpp one directory lower. I copied my peer's deprecated property_map.hpp and put it in /opt/local/include/boost/ and CGAL found the file correctly.

Everything is up and running currently with the exception that my LAPACK is complaining more than my peer's about instantiating primitives correctly. When I solve this (hopefully not with another hack) I will post the results.

Hope this helps somebody. Please let me know if you find non-hack solutions to any of these.

Update:
I think I figured out Blas and Lapack...

In my cmake CONFIG.cmake file I have


    SET(CMAKE_CXX_FLAGS "-arch i386")
    SET(CMAKE_C_FLAGS   "-arch i386")
    SET(BLAS_LIBRARIES "/System/Library/Frameworks/Accelerate.framework")
    SET(BLAS_DEFINITIONS "-DBLAS_USE_F2C")
    INCLUDE_DIRECTORIES("/opt/local/include")
    LINK_DIRECTORIES("/opt/local/lib" "lib")
    INCLUDE_DIRECTORIES(/usr/X11/include)
    LINK_DIRECTORIES(/usr/X11/lib)
    SET(TAUCS_INCLUDE_DIR "/opt/local/taucs_with_lu/src/" "/opt/local/taucs_with
    SET(TAUCS_LIBRARIES "/opt/local/taucs_with_lu/lib/darwin9.0/libtaucs.a")
    SET(METIS_LIBRARIES "/opt/local/lib/libmetis.a")

In some files it worked to just above the line where you include the Acceleration framework


#include 

put instead:


#define __LP64__
#include 

For other directories I added this line to the CMakeLists.txt file:


ADD_DEFINITIONS("-D__LP64__")

If I find an elegant way to do this I will repost and update...

Two thousand results on google for “How to turn off google buzz” only three days after release

Monday, February 22nd, 2010

2010 results for how to turn off google buzz three days after release

Only three days after its release (one day after google’s “skip this ad” on gmail), I searched for “How to turn of google buzz” (in quotes) and found about 2010 pages containing that quote.

Interleave rows of two n by m matrices, using matlab

Tuesday, February 16th, 2010

I have to 3 by 2 matrices:


     /1 2\
A = | 3 4 |
     \5 6/

     /9 8\
B = | 7 6 |
     \5 4/

and I want to interleave leave their rows into a 6 by 2 matrix, C:



     /1 2\
    | 9 8 |
C = | 3 4 |
    | 7 6 |
    | 5 6 |
     \5 4/

To achieve this in matlab I use the following:


reshape([A(:) B(:)]',size(A,1)+size(B,1), [])

Perhaps there is a better way?

3D plot from vertex and face list with vertex index labels, using matlab

Monday, February 15th, 2010

I’ve had to rediscover how this is done in matlab too many times, so here it is to be remembered. If you have a triangle mesh with 3d vertex positions V (number of vertices by 3 ) and face list F (number of faces by 3) you can plot and label using:


figure;trisurf( F,V(:,1), V(:,2),V(:,3));view(2);
text(V(:,1),V(:,2),V(:,3),num2str((1:size(V,1))'));

If your V’s are actual 2d (number of vertices by 2) and you just want to plot your two dimensional domain on a plane the use:


figure;trisurf( F,V(:,1), V(:,2),zeros(size(V,1),1));view(2);
text(V(:,1),V(:,2),zeros(size(V,1),1),num2str((1:size(V,1))'));

Areas of quadrilaterals within a triangle determined by its circumcenter (or perpendicular bisectors), using matlab, part 2

Saturday, February 13th, 2010

Yesterday, I posted some matlab code to determine the areas of quadrilaterals within a triangle determined by its circumcenter. Now, I’m updating to talk about handling obtuse triangles, where the circumcenter is outside of the triangle, i.e. the areas of the quads may be negative. So the formulas from yesterday will work for right angle triangles and acute triangles but might fail for obtuse triangles. We then need a set of areas for obtuse triangles. The areas for a right triangle and a barely right triangle should be at the limit the same, i.e. no jump. To handle this, I noticed that for a right angle triangle the circumcenter is the same as the midpoint along the edge opposite the right angle.
right triangle mid point along long edge and circumcenter
Connect the midpoints making three “quads”, on at each point (two are collapsed into triangles). Notice that all the inner triangles are congruent and have the same are thus the areas of the green, blue and yellow “quads” are area/2, area/4, area/4 corresponding to point A,B, and C, where area of course is the area of the entire triangle.
Now, do the same midpoint connections for an obtuse triangle:
obtuse triangle mid point along long edge and circumcenter
The circumcenter is no longer lined up with the midpoint along the longest side, but you do have the same congruent inner triangles and corresponding “quads”. Again the areas of the green, blue and yellow “quads” are area/2, area/4, area/4 corresponding to point A,B, and C, where area of course is the area of the entire triangle.
If you have already done the first step in my last post, to handle the right and acute triangles, then to zap out the obtuse triangles with the “corrected” quad areas use:


quads(cosines(:,1)<0,:) = [areas(cosines(:,1)<0,:)*0.5, ...
  areas(cosines(:,1)<0,:)*0.25, areas(cosines(:,1)<0,:)*0.25];
quads(cosines(:,2)<0,:) = [areas(cosines(:,2)<0,:)*0.25, ...
  areas(cosines(:,2)<0,:)*0.5, areas(cosines(:,2)<0,:)*0.25];
quads(cosines(:,3)<0,:) = [areas(cosines(:,3)<0,:)*0.25, ...
  areas(cosines(:,3)<0,:)*0.25, areas(cosines(:,3)<0,:)*0.5];

Areas of quadrilaterals within a triangle determined by its circumcenter (or perpendicular bisectors), using matlab

Thursday, February 11th, 2010

First, I need to locate the circumcenter of the triangle. The circumcenter is the center of the circle circumscribing the three points of the triangle. Thus, it is equidistance from all three points.
circumscribed triangle circumcenter
I do not need the actual world coordinates of the circumcenter, rather it is okay for me to just know the relative distance to each corner. The law of cosines tells me that the relative distance of the circumcenter from the line opposite point A (called a in the picture) is cos A (that is cosine of the angle at A). Likewise for B and C. Knowing these I can convert to barycentric coordinates by multiplying by the opposite edge lengths.


cosines = [ ...
  (edge_norms(:,3).^2+edge_norms(:,2).^2-edge_norms(:,1).^2)./(2*edge_norms(:,2).*edge_norms(:,3)), ...
  (edge_norms(:,1).^2+edge_norms(:,3).^2-edge_norms(:,2).^2)./(2*edge_norms(:,1).*edge_norms(:,3)), ...
  (edge_norms(:,1).^2+edge_norms(:,2).^2-edge_norms(:,3).^2)./(2*edge_norms(:,1).*edge_norms(:,2))];
barycentric = cosines.*edge_norms;
normalized_barycentric = barycentric./[sum(barycentric')' sum(barycentric')' sum(barycentric')'];

A quality of the barycentric trilinear coordinates of a point in a triangle is that they specify the relative area of inner triangles form by connecting the corners to that point. In our case, the triangles made by connecting the circumcenter to each corner. Here, the aquamarine, red, and purple triangles corresponding to points A, B, and C.
circumscribed triangle partial triangles
Now I find the area of the triangle and then multiple against the normalized barycentric coordinates of the circumcenter to get the areas of each of these inner triangles:


areas = 0.25*sqrt( ...
  (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3)).* ...
  (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
  (-edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
  (edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)));
partial_triangle_areas = normalized_barycentric.*[areas areas areas];

Finally the quadrilateral at each point is built from half of the neighboring inner triangles: the green, blue and yellow quads corresponding to A, B, and C.
circumscribed triangle quadrilaterals
Note, that the fact that these inner triangles are cut perfectly in half follows immediately from the fact the the circumcenter is equidistance from each corner. So the last step is just to add up the halves:


quads = [ (partial_triangle_areas(:,2)+ partial_triangle_areas(:,3))*0.5 ...
  (partial_triangle_areas(:,1)+ partial_triangle_areas(:,3))*0.5 ...
  (partial_triangle_areas(:,1)+ partial_triangle_areas(:,2))*0.5];

Note: I’ve suffered trying to use mathematic or maple to reduce these steps to a simple algebra formula (I must be mistyping something). As far as I can tell the above steps check out (of course, it only makes sense for cases where the circumcenter is within or on the triangle: i.e. the triangle is not obtuse). I have checked out my equations on maple with the following functions:


mycos:=(a,b,c)->(b^2+c^2-a^2)/(2*b*c);
barycentric:=(a,b,c)->a*mycos(a,b,c);
normalized_barycentric:=(a,b,c)->barycentric(a,b,c)/(barycentric(a,b,c)+barycentric(b,c,a)+barycentric(c,a,b));
area:=(a,b,c)->sqrt((a+b-c)*(a-b+c)*(-a+b+c)*(a+b+c))/4;                                            
inner_triangle:=(a,b,c)->normalized_barycentric(a,b,c)*area(a,b,c);
quad:=(a,b,c)->(inner_triangle(b,c,a)+inner_triangle(c,b,a))/2;

Which then you can run:


simplify(quad(a,b,c));

To see a beautiful long algebra formula for the quadrilateral at point A.
Which in matlab corresponds to:


quads = ...
  [-sqrt(-(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))).* ...
    (2.*edge_norms(:,2).^2.*edge_norms(:,3).^2  + edge_norms(:,1).^2.* ...
    edge_norms(:,2).^2  - edge_norms(:,2).^4  + edge_norms(:,1).^2.* ...
    edge_norms(:,3).^2  - edge_norms(:,3).^4)./ ...
    (8.*(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))), ...
  -sqrt(-(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))).* ...
    (2.*edge_norms(:,3).^2.*edge_norms(:,1).^2  + edge_norms(:,2).^2.* ...
    edge_norms(:,3).^2  - edge_norms(:,3).^4  + edge_norms(:,2).^2.* ...
    edge_norms(:,1).^2  - edge_norms(:,1).^4)./ ...
    (8.*(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))), ...
  -sqrt(-(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2))).* ...
    (2.*edge_norms(:,1).^2.*edge_norms(:,2).^2  + edge_norms(:,3).^2.* ...
    edge_norms(:,1).^2  - edge_norms(:,1).^4  + edge_norms(:,3).^2.* ...
    edge_norms(:,2).^2  - edge_norms(:,2).^4)./ ...
    (8.*(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2)))];

Note: It seems also that the wikipedia page for the circumscribed circle has a direct formula for the barycentric coordinates of the circumcenter. Though, wolfram’s mathworld has an ostensibly different formula. Upon a quick review Wolfram’s is at least easily verifiable to be the same as mine (recall these only need to be relative to each other).

How to turn off google buzz

Thursday, February 11th, 2010

Sure this is all over the web but why not make it extra clear. Scroll to the bottom of your gmail home page and in small print you fill find “turn off buzz”:
turn off google buzz in gmail screenshot

Average ratio of incenters to circumcenters, using matlab

Monday, February 8th, 2010

The incenter radius of a triangle is given by the formula:


      2*A      with A = area
r = -------
    (a+b+c)  


the circumcenter radius is given by


     abc
R = -----
     (4A)

We want the ratio of the incenter to circumcenter so that’s


r         8A2
_ = ---------------
R   ((a+b+c)*(abc))

So then if you have a triangle mesh with a variable edge_norms which is an n by 3 matrix of the lengths of edges in your mesh, a way to get the incenter radii is, first applying Heron’s formula to find the areas of each triangle:


semi_perimeters = (edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3))*0.5;
areas = sqrt( ...
  semi_perimeters.* ...
  (semi_perimeters-edge_norms(:,1)).* ...
  (semi_perimeters-edge_norms(:,2)).* ...
  (semi_perimeters-edge_norms(:,3)));
ratios = 8*areas.*areas./(sum(edge_norms')'.*prod(edge_norms')')


Then just take the mean:


mean(ratios)

Note: I coded up the calculations of the incenters and circumcenters to check this out so I figure I might as well post that, too:


r=2*areas./(edge_norms(:,1)+edge_norms(:,2)+edge_norms(:,3))
R=(edge_norms(:,1).*edge_norms(:,2).*edge_norms(:,3))./(4*areas)

Find maximum angle in a triangle mesh, using matlab

Monday, February 8th, 2010

If a,b,c are the lengths of the side of a triangle then the maximum angle of the triangle is the angle opposite the largest side, call that one c. A formula for this angle is then


     / a2 + b2 - c2 \
cos-1| ------------ |
     \    2*a*b    /

So then if you have a variable edge_norms which is an n by 3 matrix of the lengths of edges in your mesh, a way to find the maximum angle of each face is to issue:


acos(((sum(edge_norms'.^2)-2*max(edge_norms').^2)) ...
  ./(2*prod(edge_norms')./max(edge_norms')))

I like to make mine in degrees by adding


(...)./(2*pi).*360

And then you can find the maximum angle over all by issuing

max(...)

Surely this can be optimized, but I think matlab does a fair job on its own at that already.

Expected number of each coin denomination in change

Sunday, February 7th, 2010

Nickels are rare in change. This is verifiable even in just a experimental and anecdotal way. I look in my change bin at home and there are many quarters, many dimes and pennies and very few nickels. It also makes sense when thinking just cursorily at the likelihoods of each coin denomination popping up in your change at the super market. I was thinking about these likelihoods—more formally the expected amount of each coin after a transaction—a little (too) much and thought I might as well post my findings. So if E[X] is the expected number of times for event X to happen, then I’ll write E[pennies] to mean the expected number of pennies in my change after a transaction. I’ll assume that the total change is a fair, random number between 0 and 99 cents. I’ll also assume that the cashier will give you the minimal amount of coins to make the change.


E[quarters] = 0*25/100 + 1*25/100 + 2*25/100 + 3*25/100 = 1.5
E[dimes] = 0*10/25 + 1*10/25 + 2*5/25 = 0.8
E[nickels] = 0*5/10 + 1*5/10 = 0.5
E[pennies] = 0*1/5 + 1*1/5 + 2*1/5 + 3*1/5 + 4*1/5 = 2

So here’s the expected change for a random purchase at the store:
expected change

Note: You can even double check the math by noting that the expected total amount in change should be $0.50. And the total amount from the above coins is:


1.5*$0.25 + 0.8*$0.10 + 0.5*$0.05 + 2*$0.01 = $0.50