## Archive for October, 2014

### Triangle mesh for image as surface

Sunday, October 26th, 2014

Here’s a little chuck of matlab code I use to make a surface mesh out of an image:

% load image as grayscale
% create triangle mesh grid
[V,F] = create_regular_grid(size(im,2),size(im,1),0,0);
% Scale (X,Y) to fit image
V(:,1) = V(:,1)*size(im,2);
V(:,2) = (1-V(:,2))*size(im,1);
V(:,3) = im(:);


### A very simple method for approximate blue noise sampling on a triangle mesh

Saturday, October 25th, 2014

This morning I was looking over the recent methods for generating blue-noise random samples on an arbitrary surface. Many of them involve dart throwing with complicated rejection schemes to achieve a Poisson-disk sampling (e.g. “Feature-preserving direct blue noise sampling for surface meshes” [Peyrot et al. 2014]), some involve hefty optimization procedures (e.g. “Blue Noise through Optimal Transport” [de Goes et al. 12]). A lot of the effort is spent trying to preserve features. I wasn’t particularly interested in this at the moment, and I just wanted something simple I could quickly code up and play with.

My idea stemmed from the observation that if my surface is a sphere, and I have a enough samples then the distance between samples will approximate their distance in Euclidean space. So an algorithm for Poisson-disk sampling the sphere looks like:

1. Generate n uniformly random samples on the sphere
2. For each sample i find distance dij to closest other sample j, in usual R3 sense
3. Repeat steps 1 and 2 recursively for all samples i such that i < j and dij < re.

where re is the expected radius between samples on the surface. On the unit sphere, this is easy to compute: something like re = sqrt(4/n).

This works quite well:

on the left is the original “white noise” uniform sampling, and on the right I show the iterations of this algorithm converging to a Poisson-disk-like sampling.

Now, for an arbitrary surface my original assumption is no longer true. There may be parts very close to each other in R3 but very far away in terms of the geodesic distance on the surface.

However, there exist a family of methods that lift a surface to an (usually higher dimensional) embedding Rn so that Euclidean distance in Rn approximates geodesic distance (e.g. “Biharmonic Distance” [Lipman et al. 2010]). Using this I can alter my method only slightly to work with arbitrary surfaces:

1. Lift surface to Rn using “Biharmonic embedding”
2. Generate n uniformly random samples on the original surface (the one still in R3)
3. For each sample i find distance dij to closest other sample j, after lifting to Rn
4. Repeat steps 2 and 3 recursively for all samples i such that i < j and dij < re,

but this time re is not so obviously computed. It would be easy to compute the expected geodesic distance radius, but unfortunately the biharmonic embedding I’m using does not preserve units (i.e. dij is only relative to the geodesic distance between samples i and j on the mesh). I could probably compute the expected distortion during the lift to Rn and derive re from that. But I found it much easier to just take the average closest distance between the original n samples after being lifted to Rn. If I recall correctly, this should converge to the correct expected radius as n goes to infinity. So for large n this should be a reasonable thing to do.

It’s working quite well:

again, on the left is the original “white noise” uniform sampling, and on the right I show the iterations of this algorithm converging to a Poisson-disk-like sampling.

I’ve put up matlab code for this in the mesh/random_points_on_mesh.m in my gptoolbox on git. Here’s a use example:

[V,F] = load_mesh('/usr/local/igl/libigl/examples/shared/cheburashka.off');
n = 3000;
Pwhite = random_points_on_mesh(V,F,n);
Pblue = random_points_on_mesh(V,F,n,'Color','blue');
t = tsurf([F;size(V,1)+F],[V;V(:,1)+1.1*(max(V(:,1))-min(V(:,1))) V(:,2:3)],'FaceColor','w','EdgeColor','none','FaceAlpha',0.8);
apply_ambient_occlusion(t);
hold on;
scatter3(Pblue(:,1),Pblue(:,2),Pblue(:,3));
scatter3(Pwhite(:,1)+1.1*(max(V(:,1))-min(V(:,1))),Pwhite(:,2),Pwhite(:,3));
hold off;
axis equal;


Note: The current code is rebuilding a k-nearest neighbor acceleration data structure for all samples every iteration. Surely, this could be optimized.

### The Freedom Tower is a coarse triangle mesh

Sunday, October 19th, 2014

Let’s apply Loop subdivision to make it smoother.

### Extruding a Bezier curve into a triangle mesh in maya

Friday, October 17th, 2014

Today I struggled to convince Maya to let me extrude a Bezier Curve into a solid shape (sweep a closed curve and finish with planar end caps). I could used the Surface > Extrude tool to extrude the curve and then select the boundary edges and use the Surface > Planar tool to close the endcaps, but this just creates a group of 3 surfaces which are not topologically connected.

My end goal today was to create something to send to the 3D printer. So in this case I eventually wanted a triangle mesh. Here’re the steps I took to convert a bezier curve to a polygonal mesh:

1. draw bezier curves
2. find the Polygons > Plane tool
3. draw a plane behind the curves
4. Select one of the curves and the plane
5. Choose Edit Mesh > Project curve onto mesh
6. Select the new projected curve and the plane
7. Choose Edit Mesh > Split Mesh with projected curve
8. Right click and hold and drag to select “Face” selection mode
9. Mouse over the plane until just the filled curve appears (there seem to be many overlapping faces.
10. Choose Edit > Invert Selection
11. Then choose Edit > Delete
12. Select just the original Bezier curve and delete it.
13. Repeat steps 2-12 for the other curves (why can’t we do all curves at once?)
14. Select both filled curves and choose the Polygons > Extrude tool
15. Pull down on the widget’s arrow to extrude.
16. Select both extruded surfaces and choose Mesh > Cleanup...
17. Make sure 4-sided faces is unchecked
18. Make sure Faces with zero geometry area is checked with very small Area tolerance (e.g. 0.00001)
19. Hit Cleanup
20. The choose Mesh > Triangulate
21. The surface is now triangulated. Select everything.
22. File > Export Selection and save as an obj

Wow. So 21 steps. Not particularly easy for a task I thought would be dead simple. I must be missing some faster way to do this.

### Set default axis and figure colors

Wednesday, October 15th, 2014

Here are two commands to permanently change the axes and figure background colors.

set(groot,'DefaultAxesColor',[0.94 0.94 0.94]);
set(groot,'DefaultFigureColor',[0.75 0.75 0.75]);


I find this new combination a bit less abrasive to the eyes.

### Continuous Slider in Matlab

Wednesday, October 15th, 2014

Here’s the typical way to create a slider in matlab:

sld = uicontrol('Style', 'slider','UserData',10,'Callback',@(src,data) disp(src.UserData * src.Value));


The callback here is boring but shows how to use 'UserData' to pass information to the callback: here, it’s just multiplying the slider’s current value by 10 and displaying the result.

However, the callback only triggers on the mouse up event. This is rather unsatisfying. To fix this, add this generic listener to the slider, which bootstraps the current callback:

sld.addlistener('Value','PostSet',@(src,data) data.AffectedObject.Callback(data.AffectedObject,struct('Source',data.AffectedObject,'EventName','Action')));


Now, the callback is firing on the mouse drag events. It’s generic in the since that you shouldn’t need to change this line if you used a different callback. It will even continue to work if you change the callback during interaction.

I’m aping the data parameter in the uicontrol callback with a simple struct. I’m not really sure what this data parameter is good for anyway.

### Ambient occlusion + anti-aliasing in MATLAB

Wednesday, October 15th, 2014

I’m enjoying the new anti-aliased graphics in matlab 2014b. Here’s the xyzrgb dragon rendered with soft lighting, ambient occlusion and a simple colormap:

Here’s the code to reproduce it:

[V,F] = load_mesh('/Users/ajx/Documents/AutoDOF/Code/skinning/skinning/xyzrgb_dragon/xyzrgb_dragon.obj');
AO = ambient_occlusion(V,F,V,per_vertex_normals(V,F),1000);
colormap(parula(9));
t = tsurf(F,V,'EdgeColor','none');
C = squeeze(ind2rgb(floor(matrixnormalize(t.FaceVertexCData)*size(colormap,1))+1,colormap));
t.FaceVertexCData = bsxfun(@times,C,1-AO);
t.SpecularStrength = 0.1;
t.DiffuseStrength = 0.1;
t.AmbientStrength = 0.8;
l = light('Position',[1 1 100],'Style','infinite');
l2 = light('Position',[1 -100 1],'Style','infinite');
set(gca,'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[],'Color',[0.94 0.94 0.94]);
set(gcf,'Color','w');
axis equal
camproj('persp');
t.Vertices = V*axisangle2matrix([0 0 1],pi)*axisangle2matrix([1 0 0],pi/2);
view(-43,6);
axis tight;
drawnow;


The ambient_occlusion call takes quite some time on my little laptop. But I think the result looks nice.

### Read-only proxy class for primitives in C++

Tuesday, October 14th, 2014

One of my least favorite experiences coding in C++ is writing boilerplate setters and getters. Usually I won’t write setters unless I need side effects. For that, read-and-write case I’ll just use public member functions. The most common case I run into is when I have members I’d like to expose as read-only. That is, I want to be able to get the value of a private or protected member. So I’ll end up with something like this:

class MyFancyClass
{
private:
int m_var1;
double m_var2;
...
double m_varn;
public:
int get_var1(){return m_var1}
double get_var2(){return m_var2}
...
double get_varn(){return m_varn}
} a;


This is annoying. Anytime I add a new member m_varnplus1 I’ll need to add a corresponding MyFancyClass::get_varnplus1(). It’s also clunky from the caller’s point of view. I can’t just use something like x = a.m_var1, instead I need to use the caller x = a.get_var1(). Often I find myself getting to lazy to do all this so I just sloppily use:

class MyFancyClass
{
public:
int m_var1;
double m_var2;
...
double m_varn;
} a;


Boom. Fuck it. Everything’s gonna be public.

I’m toying with the idea of widely employing this proxy class for primitives. I’d save the following in RO.h:

#ifndef RO_H
#define RO_H
// A proxy class for creating "read only" members: member
//
template<typename Container, typename Primitive>
class RO
{
friend Container;
public:
inline operator Primitive() const             { return x; }
template<typename Other> inline bool  operator==(const Other& y) const { return x == y; }
template<typename Other> inline bool  operator!=(const Other& y) const { return x != y; }
template<typename Other> inline bool  operator< (const Other& y) const { return x < y; }
template<typename Other> inline bool  operator> (const Other& y) const { return x > y; }
template<typename Other> inline bool  operator<=(const Other& y) const { return x <= y; }
template<typename Other> inline bool  operator>=(const Other& y) const { return x >= y; }
template<typename Other> inline Other operator+ (const Other& y) const { return x + y; }
template<typename Other> inline Other operator- (const Other& y) const { return x - y; }
template<typename Other> inline Other operator* (const Other& y) const { return x * y; }
template<typename Other> inline Other operator/ (const Other& y) const { return x / y; }
template<typename Other> inline Other operator<<(const Other& y) const { return x <<y; }
template<typename Other> inline Other operator>>(const Other& y) const { return x >> y; }
template<typename Other> inline Other operator^ (const Other& y) const { return x ^ y; }
template<typename Other> inline Other operator| (const Other& y) const { return x | y; }
template<typename Other> inline Other operator& (const Other& y) const { return x & y; }
template<typename Other> inline Other operator% (const Other& y) const { return x % y; }
template<typename Other> inline Other operator&&(const Other& y) const { return x && y; }
template<typename Other> inline Other operator||(const Other& y) const { return x || y; }
template<typename Other> inline Other operator~() const                { return ~x; }
template<typename Other> inline Other operator!() const                { return !x; }
protected:
inline RO()                                  :x(){ }
template<typename Other> inline RO(const Other& y)                    :x(y) { }
template<typename Other> inline Primitive& operator= (const Other& y) { return x = y; }
template<typename Other> inline Primitive& operator+=(const Other& y) { return x += y; }
template<typename Other> inline Primitive& operator-=(const Other& y) { return x -= y; }
template<typename Other> inline Primitive& operator*=(const Other& y) { return x *= y; }
template<typename Other> inline Primitive& operator/=(const Other& y) { return x /= y; }
template<typename Other> inline Primitive& operator&=(const Other& y) { return x &= y; }
template<typename Other> inline Primitive& operator|=(const Other& y) { return x |= y; }
template<typename Other> inline Primitive& operator%=(const Other& y) { return x %= y; }
template<typename Other> inline Primitive& operator++()               { return ++x; }
inline Primitive  operator++(int)            { return x++; }
template<typename Other> inline Primitive& operator--()               { return --x; }
inline Primitive  operator--(int)            { return x--; }
inline Primitive* ptr()                      { return &x; }
inline Primitive& ref()                      { return x; }
inline const Primitive* const_ptr() const    { return &x; }
inline const Primitive& const_ref() const    { return x; }
Primitive x;
};

#endif


Then I could use it with something like:

class MyFancyClass
{
typedef MyFancyClass Self;
public:
RO<Self,int> m_var1;
RO<Self,double> m_var2;
...
RO<Self,double> m_varn;
void example_of_write_access()
{
m_var1 = 10;
m_var2 = 100;
m_varn = m_var2++;
}
} a;


Outside this class I’m only allowed read operations int x = a.m_var1 or if(a.m_var1) etc. But I can’t write to these members: can’t write a.m_var2 = 3 or a.m_varn += 7.

I can use my members as usual inside my class with full read and write access: e.g. m_var1 = 10 or m_var2++.

These read only members must be declared public to be readable outside the container class. Declaring them private would be pointless. However, declaring them protected will make them read-only to friends of the container. Thus the public/protected/private class access modifiers effectively determine the read-access, and write-access is always private.

Lingering issue: I’m not sure how you’d create a publicly read-only member who’s write-access was protected (shared with friends and derived classes). Somehow I’d have to convince RO to be friends also with the inherited class, which isn’t known during the definition of the base container class. Seems since friendship is neither inherited nor transitive, this is not going to be easy.

### Burying the rainbow

Monday, October 13th, 2014

I come to bury Caesar, not to praise him.

MATLAB’s newest version has finally tossed the jet default colormap for parula.

The visualization community has long been warning against the use of “rainbow” colormaps like jet. Today I looked around for some of the papers articulating why.

The well-cited, “Rainbow Color Map (Still) Considered Harmful” by David Borland and Russell M. Taylor II presents synthetic counterexamples, highlighting why the “Rainbow” colormap should be avoided.

The first counterexample cites the book Information Visualization: Perception for Design by C. Ware. From the text in [Borland and Taylor], it seems this book contains details of a perceptual study in which subjects correctly order grayscale paint chips consistently but often fail to order “rainbow” paint chips (presumably by hue).

Borland and Taylor go on to write,

Some even put them in alphabetical order. To put them in the so-called correct order, most people must remember Roy G. Biv (red, orange, yellow, green, blue, indigo, violet), or some other mnemonic representation of the order of colors in the rainbow.

I was a bit surprised. I’m quite accustomed to matlab’s jet and would have never chosen the order red, green, yellow, blue. I could admit that perhaps I’m perceptually faster at ordering grayscale colors that hue-based spectra. So I went digging in the Ware book to find the original study to find more details.

In Chapter 4, Application 3 “Application 3: Color Sequences for Data Maps” Ware cites an earlier paper of his called “Color Sequences for Univariate Maps: Theory, Experiments, and Principles”. This paper contains three small scale studies providing evidence against chromatic sequences: one quantitative study and two qualitative survey-style studies. The studies and results are aimed specifically at a pseudocoloring’s ability to represent a height field or surface.

The other cited study is “Face-based Luminance Matching for Perceptual Colormap Generation” by [Kindlmann et al.]. This paper seems to contain a user study regarding their proposed system for luminance calibration and doesn’t seem directly related to Ware’s claim that “Experimental studies have confirmed that grayscale maps are much better for form perception.”

Their more papers cited in “How NOT to Lie with Visualization” by [Rogowitz and Treinish] backing the finding that the “[rainbow hue colormap] produces several well-documented artifacts”. Some of these artifacts are convincing without a perceptual study: e.g. focus is drawn to bright, yellow regions of this MRI while subtleties are drown out in the large green span in the spectrum:

Question: Is there no large scale study quantitatively confirming the perceptual limitations of the rainbow colormap?

Maybe now that Matlab’s changed its default everybody (including probably me) will just stop using it out of laziness and we won’t need to worry about verifying these claims empirically.

I continued through “Rainbow Color Map (Still) Considered Harmful” by David Borland and Russell M. Taylor II to the second counterexample. Here an scalar field varying frequency is mapped to a grayscale colormap and a rainbow” one:

This example immediately puzzled me. Why is so much lost in the green? I would not expect to see this in matlab using either its jet or hsv rainbow colormaps.

Indeed, this counterexample seems particular the choice of “rainbow” colormap.

As a testament to the simplicity of the grayscale colormap, we can directly recover the underlying scalar field by taking the intensity value of the grayscale image in Fig 2 left. Using this, we can recover the colormap used for the “rainbow” visualization from their Fig 2 right. Then, to verify I reproduce the input in matlab above.

Now, I can easily try other colormaps like jet and hsv:

Neither suffer from the washout of the original “rainbow” visualization. Though there’s plenty of room to argue about whether too much attention is drawn from the vertical stripes to the extrema values on the bottom.

Looking at just this scalar field makes it difficult to judge the colormaps being used. Here I’ve taken the domain of scalar values (roughly 0 to 250) linearly and applied colormap:

The top right shows the recovered “rainbow” colormap. Notice the wide range of green values and lack of purple values.

Figure 3 in the same paper shows a “rainbow” colormap with a wider gamut (left):

The recovered colormap seems is well approximated by a truncated versioned of this wider colormap (right).

If we apply the full “rainbow” from Fig 3 to the original data we get a bit more reasonable visualization (left):

whereas the truncated colormap produces a washed out visualization (right) similar to the original one in Figure 2.

It seems rather obvious that using a colormap with such a wide span covered by essentially the same green will be troublesome. Here’re the spans in the recovered colormap, the full “rainbow”, and the truncated rainbow which are within 30.72 RGB8 units from the central “green” value:

The jet and hsv colormaps don’t contain this same green. hsv seems to suffer from a similar problem, but jet to a much lesser degree (more cyan-ish colors). This finding seems to agree with the [Kindlmann et al.] idea that there are “good” rainbows and “bad” rainbows. Though it doesn’t quite count as evidence that all rainbows are bad, (even if we all know that they must be).

To come full circle, here’s that dataset visualized with MATLAB’s new parula colormap:

I’m sold on it. I agree with Robert Kosara’s explanation why rainbows are popular:

Given the issues, why are the rainbow colormap and its variants so popular? I think the answer is quite simple: it’s attractive. Using a single hue to show the data would be reasonably effective, but much less interesting to look at.

parula seems to offer a decent trade-off. I don’t think I’ll mind staring it from now on. Perhaps, I’ll still silently long to see those other hues of the rainbow.

### MATLAB2014b features anti-aliasing

Monday, October 13th, 2014

Finally. I’m pretty happy about the results:

[V,F] = load_mesh('/usr/local/igl/libigl/examples/shared/cheburashka.off');
AO = ambient_occlusion(V,F,V,per_vertex_normals(V,F),1000);
t = tsurf(F,V,fphong,'EdgeColor','none');
C = squeeze(ind2rgb(floor(matrixnormalize(t.FaceVertexCData)*size(colormap,1))+1,colormap));
t.FaceVertexCData = bsxfun(@times,C,1-AO);
t.SpecularStrength = 0.1;
t.DiffuseStrength = 0.1;
t.AmbientStrength = 0.8;
l = light('Position',[1 1 100],'Style','infinite');
l2 = light('Position',[1 -100 1],'Style','infinite');
set(gca,'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[],'Color',[0.94 0.94 0.94]);
set(gcf,'Color','w');


And to spin the camera around:

axis equal
axis vis3d;
camproj('persp');
for f = 1:numel(T)
t = T(f);
view(-cos(t*2*pi)*45,sin(t*2*pi)*45+45);
drawnow;
frame = getframe(gcf);
[SIf,cm] = rgb2ind(frame.cdata,256);
if f == 1
imwrite(SIf,cm,filename,'Loop',Inf,'Delay',0);
else
imwrite(SIf,cm, filename,'WriteMode','append','Delay',0);
end
end


With the awesome but now obsolete myaa.m hacked anti-aliasing, creating this gif would have taken many minutes. This runs in real time.