## Posts Tagged ‘matlab’

### Offset surface of triangle mesh in matlab

Wednesday, March 22nd, 2017

Here’s a little demonstration of how to use gptoolbox and MATLAB to generate an offset surfaces of a triangle mesh. This takes a mesh in V,F and creates a mesh SV,SF of the isosurface at signed distance iso:

% Extract offset at minus 3% of bounind box diagonal length
iso = -0.03;
% Resolution grid → resolution of extracted offset surface
side = 60;
% Amount of smoothing to apply to distance field
sigma = 1.4;
bbd = norm(max(V)-min(V));
% Pad generously for positive iso values
[BC,side,r] = voxel_grid([V;max(V)+iso*1;min(V)-iso*1],side);
D = signed_distance(BC,V,F);
D = reshape(D,side([2 1 3]));
% Smooth signed distance field
D = imfilter(D,fspecial('gaussian',9,sigma),'replicate');
BC3 = reshape(BC,[side([2 1 3]) 3]);
% Use matlab's built-in marching cubes iso-surface mesher (works most of the time)
surf = isosurface(BC3(:,:,:,1),BC3(:,:,:,2),BC3(:,:,:,3),D,iso*bbd);
SV = surf.vertices;
SF = surf.faces;


Here’s a blue bunny with a positive offset surface, an orange “cage”:

Here’s a blue bunny with a negative offset surface. This is useful for hollowing out objects to 3d print:

Because the iso-surface extraction will over tesselate low curvature patches of the output surface, it would make a lot of sense to remesh/decimate this mesh.

(to create these fancy renderings:)

clf;
hold on;
t  = tsurf(F,V,'EdgeColor','none',fsoft,  'FaceVertexCData',repmat(blue,size(V,1),1),'FaceAlpha',1+(iso<0)*(0.35-1),fphong);
ts = tsurf(SF,SV,'EdgeAlpha',0.2+(iso<0)*(0-0.2),fsoft,'FaceVertexCData',repmat(orange,size(SV,1),1),fphong,'FaceAlpha',1+(iso>0)*(0.2-1));
apply_ambient_occlusion(ts);
hold off;
axis equal;
view(-20,20)
camlight;
t.SpecularStrength = 0.04;
l = light('Position',[5 -5 10],'Style','local');
set(gca,'pos',[0 0 1 1])
set(gca,'Visible','off');
set(gcf,'Color','w');
drawnow;


### Quad meshing in matlab using gptoolbox

Friday, February 17th, 2017

This morning I played around with a very simple way of quad meshing a triangle mesh. This is more or less the “ideal pipeline” for quad meshing via parameterization. It will likely only work for very simple meshes. I’m pleasantly surprised that it works at all.

% Load a mesh with boundary and no topological (donut) holes
V = V*axisangle2matrix([1 0 0],-pi/2);
% Construct the LSCM matrix
[~,Q] = lscm(V,F,[],[]);
M = massmatrix(V,F);
% Use Fiedler vector as parameterization
[EV,ED] = eigs(Q,repdiag(M,2),5,'sm');
U = reshape(EV(:,end-3),[],2);
% try to find "canonical" rotation
[sU,sS,sV] = svd(U'*U);
U = U*sU;
% Lay down a grid in 2D
h = ((max(U(:,1))-min(U(:,1)))/(160-2));
[X,Y] = meshgrid(min(U(:,1))-h:h:max(U(:,1))+h,min(U(:,2))-h:h:max(U(:,2))+h);
% Find all quads that are strictly inside the parameterized mesh
[QQ,QU] = surf2patch(X,Y,0*Y);QU = QU(:,1:2);
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
QQ = QQ(all(abs(sqrD(QQ))<1e-10,2),:);
[QU,I] = remove_unreferenced(QU,QQ);
QQ = fliplr(I(QQ));
% Snap boundary vertices of quad mesh to boundary of mesh, and
% solve little Dirichlet problem to diffuse the displacements smoothly
QF = [QQ(:,[1 2 3]);QQ(:,[1 3 4])];
L = cotmatrix(QU,QF);
b = unique(outline(QQ));
[~,~,QU(b,:)] = point_mesh_squared_distance(QU(b,:),U,outline(F));
% Locate each quad mesh vertex in the parameterized mesh
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
B = barycentric_coordinates(C(:,1:2),U(F(I,1),:),U(F(I,2),:),U(F(I,3),:));
% Interpolate 3D positions at quad mesh vertex locations
QV = V(F(I,1),:).*B(:,1) + V(F(I,2),:).*B(:,2) + V(F(I,3),:).*B(:,3);


To create the image above use:

clf;
hold on;
tq = tsurf(QQ,QV,'FaceVertexCData',repmat(blue,size(QV,1),1),'EdgeColor','none',fsoft,fphong);
tt = tsurf(F,V-[0.8 0 0],'FaceVertexCData',repmat(1-0.8*(1-blue),size(V,1),1),'EdgeAlpha',0.5,fsoft,fphong);
axis equal;
camproj('persp');
view(-38,17);
drawn;
camlight;
plot_edges(QV,[QQ(:,2:3);QQ(:,[4 1])],'Color',orange);
plot_edges(QV,[QQ(:,1:2);QQ(:,3:4)],'w');hold off;
l = light('Position',[-10 0 10],'Style','infinite');
set(gca,'Visible','off');
set(gcf,'Color','w');


There’s a lot to improve about the quad mesh: the boundaries are not aligned with the meshing direction, the meshing directions are not always aligned with principle curvature, etc. Not bad for 30 mins of coding using tools that are just lying around.

### How does matlab’s hypot function work?

Saturday, January 7th, 2017

The hypot function in matlab purports to be more numerically stable at computing the hypotenuse of a (right-)triangle in 2D. The example the help hypot gives is:

a = 3*[1e300 1e-300];
b = 4*[1e300 1e-300];
c1 = sqrt(a.^2 + b.^2)
c2 = hypot(a,b)


where you see as output:

c1 =

Inf     0


and

c2 =

5e+300       5e-300


this is a compiled built-in function so you can’t just open hypot to find out what its doing. It might just be pre-dividing by the maximum absolute value. There’s probably a better reference for this, but I found it in: “Vector length and normalization difficulties” by Mike Day.

Continuing this example:

m = max(abs([a;b]))
c3 = m.*sqrt((a./m).^2 + (b./m).^2)


produces

c3 =

5e+300       5e-300


While matlab’s hypot only accepts two inputs, pre-dividing by the maximum obviously extends to any dimension.

Friday, January 6th, 2017

Like clockwork, my matlab updates have gotten out of sync with Xcode updates. It seems like fixing this SDK error always requires a different hack each year. This year I got the error:

No supported compiler or SDK was found. For options, visit
http://www.mathworks.com/support/compilers/current_release/.


To fix it, I replaced all occurrences of 10.9 with 10.11 in /Applications/MATLAB_R2017a.app/bin/maci64/mexopts/clang{++,}_maci64.xml

ld: warning: object file was built for newer OSX version (10.11) than being linked (10.9)


For now, I’m assuming that I can ignore them. We’ll see how far that gets me.

### Plot a bunch of edges in matlab with per-edge color data

Monday, December 19th, 2016

Suppose you have a list of vertices V (#V by dim) and a list of edges E (#E by 2) indexing V, then one hacky to draw all of these edges with per-vertex color data DV (#V by 1) is to use trisurf:

trisurf(E(:,[1 1 2]),V(:,1),V(:,2),V(:,3),'CData',DV,'EdgeColor','interp');


However, if you try to use this trick to draw per-edge color data DE (#E by 1), you’ll find that 'EdgeColor','flat' does not work. Instead you can explode you compact “mesh” of edges into individual edges and repeat you edge color data for each edge point:

trisurf(reshape([1:size(E,1),1:size(E,1)*2],[],3),V(E(:),1),V(E(:),2),V(E(:),3),'CData',repmat(DE,2,1),'EdgeColor','interp');


### Signed polygon area in matlab

Sunday, October 30th, 2016

Suppose you have a polygon’s 2D corners stored in P, then the signed area is given by:

signed_polyarea = @(P) 0.5*(P(1:end,1)'*P([2:end 1],2)-P(1:end,2)'*P([2:end 1],1));


### Fix dyld linker errors when installing new mosek toolbox

Monday, October 3rd, 2016

Each time I upgrade my mosek library, matlab panics and can’t find it. The old solution was to monkey with the DYLD_LIBRARY_PATH in all sorts of funny places: ~/.profile, /etc/launch.d, ~/Library/LaunchAgents/environment.plist

These all worked at some point, but as far as I can tell, no longer do. They’re also the wrong way to install libraries.

Fortunately, mosek has made life easy. Just

cd /usr/local/mosek/8/tools/platform/osx64x86/bin/
python install.py


This will actually fix all of the binaries and mex files using otool and install_name_tool to find dynamic libraries in their installed locations.

### (Slightly) Faster way to compute number of unique elements in matlab matrix

Wednesday, August 31st, 2016

The standard way to compute the number of unique entries in a matlab matrix A is:

numel(unique(A))


If our entries are positive integers, we can try to do the same thing using sparse with:

nnz(sparse(A(:),1,1,size(A,1),1))


but actually this is slower.

I don’t see a way to avoid a sort. I came up with this,

sum(diff(sort(F(:)))~=0)+1


As far as I can tell, this will work for matrices that don’t have infs and nans. It’s slightly faster than number(unique(A)). I have a feeling I’m only winning anything here because I’m avoiding overhead within unique

### Inverse of common ease curves

Wednesday, July 13th, 2016

Awkwardly I ended up needing the inverse of an ease curve (or S-curve) used in tweening animations. For trigonometric ease curves, this is easy. For example, if your ease filter is:

f = 0.5-cos(x*pi)*0.5;


then the inverse is:

x = acos(-2*(f-0.5))/pi


But if you’re using the famous cubic ease curve,

f = 3.*x.^2 - 2.*x.^3


then the relevant inverse is a complex function (involving the i = sqrt(-1)) that produces real values for f in [0,1]:

x = -1./4.*(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)-1./4./(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)+1./2-1./2.*i.*3.^(1./2).*(1./2.*(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)-1./2./(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3));


I haven’t bothered to see if this can be simplified into something tidy. It’d be great to get rid of the i.

### A simple, cheapskate super resolution implementation in pure matlab

Thursday, May 26th, 2016

As a proof of concept, I implemented a very simple image upsampling/super resolution algorithm. This method uses local self-examples: it replaces each pixel with a 2×2 block of pixels interpolated between blocks from the same image. The blocks are identified based on matching the local regions around.

If the window radius =1, then I first create a point in Wij in R9 for every pixel (i,j) so that W22 = [I11,I12,I13,I21,I22,I23,I31,I32,I33], where Iij is the color at pixel ij.

Then for every 2×2 block Dij centered at the bottom right corner of pixel (i,j) I create a point in R9, so that D22 = [0.25(I00 + I01 + I10 + I11), 0.25(I02 + I03 + I12 + I13), … , 0.25*(I44 + I45 + I54 + I55)]. So the points Dij represent downsampled blocks.

The above, treats all coordinates (neighboring pixel values) equally. I found that it’s best to weight the influence by a Gaussian so that the center pixel/block is worth more.

Then for each Wij I identify the closest K Dij using knnsearch. I replace the pixel (i,j) with an interpolation of the K blocks. This is done using Shepard interpolation with distances measured in the weighted R9 space.

Flipping around through the literature, this seems in spirit with 2000-era Freeman-type super-resolution algorithms. And the results are comparable. But by no means are these as good as state of the art heavy-learning based techniques around today.

Find the code for imupsample in gptoolbox.

Here’re some examples:

And here’s a bigger example:

There’s still a bathroom window effect happening, but with more local samples to choose from, up sampling bigger images seems to work better.