Split a long mp3 audio file into 3 min files

July 29th, 2015

After some frustration trying to get mp3splt to work, I caved in and wrote a script to split apart a large audio file into many 3min chunks. Save this in mp3split.sh:

#!/bin/bash
big="$1"
duration_stamp=$(ffmpeg -i "$big" 2>&1 | grep Duration | sed 's/^.*Duration: *\([^ ,]*\),.*/\1/g')
title=$(ffmpeg -i "$big" 2>&1  | grep "title *:" | sed 's/^.*title *: *\(.*\)/\1/g')
# get minutes as a raw integer number (rounded up)
prefix=$(basename "$big" .mp3)
echo $duration_stamp
mins=$(echo "$duration_stamp" | sed 's/\([0-9]*\):\([0-9]*\):\([0-9]*\)\.\([0-9]*\)/\1*60+\2+\3\/60+\4\/60\/100/g' | bc -l | python -c "import math; print int(math.ceil(float(raw_input())))")
ss="0"
count="1"
total_count=$(echo "$mins/3+1" | bc)
while [ "$ss" -lt "$mins" ]
do
  zcount=$(printf "%05d" $count)
  ss_hours=$(echo "$ss/60" | bc)
  ss_mins=$(echo "$ss%60" | bc)
  ss_stamp=$(printf "%02d:%02d:00" $ss_hours $ss_mins)
  ffmpeg -i "$big" -acodec copy -t 00:03:00 -ss $ss_stamp -metadata track="$count/$total_count" -metadata title="$title $zcount" $prefix-$zcount.mp3 
  ss=$[$ss+3]
  count=$[$count+1]
done

The execute mp3split.sh my-long-file.mp3. This will output a sequence of files:

my-long-file-00001.mp3
my-long-file-00002.mp3
my-long-file-00003.mp3
my-long-file-00004.mp3
...

Each will retain the meta data from the original file except the file number will be appended to the track name and the track number will be set accordingly (i.e. this will work well for splitting enormous audiobook files into file lists that play in the correct sequence on an iphone).

Note: mp3splt really seems like the right tool for this. It supposedly has fancy features like silence detection and presumably won’t reload the file for each new split.

Remove annotations from a pdf

July 29th, 2015

I lost track of where I found this online and had to dig it out of my bash history:

perl -pi -e 's:/Annots \[[^]]+\]::g' input-and-output.pdf

Libigl tutorial entry for physics upsampling with biharmonic coordinates

July 28th, 2015

I’ve added a tutorial entry describing biharmonic coordinates and linking to an example using it to run a “physics” simulation on a coarse mesh and upsample it to a high-resolution mesh:

octopus physics biharmonic coordinates

Read more about in our upcoming SIGGRAPH paper Linear Subspace Design for Real-Time Shape Deformation.

Sanity check that constructing a parameter inline circumvents aliasing functions with Eigen

July 27th, 2015

Suppose I have a function with Eigen types as inputs and outputs:

void aliasing_function(
  const Eigen::MatrixXd & A,
  Eigen::MatrixXd & B)
{
  for(int i = 0;i<B.rows();i++)
  {
    B.row(i) = A.row(A.rows()-i-1);
  }
}

This function is dangerous despite the const on the input parameter A because if A and B reference the same place in memory, writing to B will also write onto what this function thinks is A. Usually this is undesirable and it’s called aliasing. For example consider calling with the same input as output:

aliasing_function(C,C);

Supposing that I can’t/won’t/refuse to change the implementation of aliasing_function, here’s proof that you can avoid aliasing problems by using the copy constructor of the Eigen type directly in the parameter (avoiding the explicitly named temporary variable, but still costing a copy probably). The call above is replaced with:

aliasing_function(Eigen::MatrixXd(C),C);

To have a tangible example:

int main()
{
  using namespace std;
  Eigen::MatrixXd C = (Eigen::MatrixXd(3,3)<<1,2,3,4,5,6,7,8,9).finished();
  cout<<"C: "<<endl<<C<<endl;
  aliasing_function(Eigen::MatrixXd(C),C);
  cout<<"C: "<<endl<<C<<endl;
  aliasing_function(C,C);
  cout<<"C: "<<endl<<C<<endl;
}

This outputs:

C: 
1 2 3
4 5 6
7 8 9
C: 
7 8 9
4 5 6
1 2 3
C: 
1 2 3
4 5 6
1 2 3

The last print out shows the effect of aliasing.

Retrieve system volume and mute status via objective-c

July 25th, 2015

Here’s a little objective-c program to grab the Mac OS X system volume and mute status (and then multiply them to get the effective volume).

#import <AudioToolbox/AudioServices.h>
#import <Foundation/NSObject.h>

@interface Volume : NSObject
{
}
 +(AudioDeviceID)defaultOutputDeviceID;
 +(float)volume;
 +(bool)mute;
 +(float)effective_volume;
@end

@implementation Volume : NSObject
+(AudioDeviceID)defaultOutputDeviceID;
{
    AudioDeviceID   outputDeviceID = kAudioObjectUnknown;
    // get output device device
    UInt32 propertySize = 0;
    OSStatus status = noErr;
    AudioObjectPropertyAddress propertyAOPA;
    propertyAOPA.mScope = kAudioObjectPropertyScopeGlobal;
    propertyAOPA.mElement = kAudioObjectPropertyElementMaster;
    propertyAOPA.mSelector = kAudioHardwarePropertyDefaultOutputDevice;
    if (!AudioHardwareServiceHasProperty(kAudioObjectSystemObject, &propertyAOPA))
    {
        NSLog(@"Cannot find default output device!");
        return outputDeviceID;
    }
    propertySize = sizeof(AudioDeviceID);
    status = AudioHardwareServiceGetPropertyData(kAudioObjectSystemObject, &propertyAOPA, 0, NULL, &propertySize, &outputDeviceID);
    if(status) 
    {
        NSLog(@"Cannot find default output device!");
    }
    return outputDeviceID;
}

+(float)volume 
{
    Float32         outputVolume;
    UInt32 propertySize = 0;
    OSStatus status = noErr;
    AudioObjectPropertyAddress propertyAOPA;
    propertyAOPA.mElement = kAudioObjectPropertyElementMaster;
    propertyAOPA.mSelector = kAudioHardwareServiceDeviceProperty_VirtualMasterVolume;
    propertyAOPA.mScope = kAudioDevicePropertyScopeOutput;
    AudioDeviceID outputDeviceID = [[self class] defaultOutputDeviceID];
    if (outputDeviceID == kAudioObjectUnknown)
    {
        NSLog(@"Unknown device");
        return 0.0;
    }
    if (!AudioHardwareServiceHasProperty(outputDeviceID, &propertyAOPA))
    {
        NSLog(@"No volume returned for device 0x%0x", outputDeviceID);
        return 0.0;
    }
    propertySize = sizeof(Float32);
    status = AudioHardwareServiceGetPropertyData(outputDeviceID, &propertyAOPA, 0, NULL, &propertySize, &outputVolume);
    if (status)
    {
        NSLog(@"No volume returned for device 0x%0x", outputDeviceID);
        return 0.0;
    }
    if (outputVolume < 0.0 || outputVolume > 1.0) return 0.0;
    return outputVolume;
}

+(bool)mute
{
    UInt32 mute;
    UInt32 propertySize = 0;
    OSStatus status = noErr;
    AudioObjectPropertyAddress propertyAOPA;
    propertyAOPA.mElement = kAudioObjectPropertyElementMaster;
    propertyAOPA.mSelector = kAudioDevicePropertyMute;
    propertyAOPA.mScope = kAudioDevicePropertyScopeOutput;
    AudioDeviceID outputDeviceID = [[self class] defaultOutputDeviceID];
    if (outputDeviceID == kAudioObjectUnknown)
    {
        NSLog(@"Unknown device");
        return 0.0;
    }
    if (!AudioHardwareServiceHasProperty(outputDeviceID, &propertyAOPA))
    {
        NSLog(@"No volume returned for device 0x%0x", outputDeviceID);
        return 0.0;
    }
    propertySize = sizeof(UInt32);
    status = AudioHardwareServiceGetPropertyData(outputDeviceID, &propertyAOPA, 0, NULL, &propertySize, &mute);
    if (status)
    {
        NSLog(@"No volume returned for device 0x%0x", outputDeviceID);
        return 0.0;
    }
    return mute;
}
+(float)effective_volume
{
  return [[self class] volume] * (![[self class] mute]);
}

@end

int main(int argc, const char * argv[])
{
  printf("%g",[Volume effective_volume]);
}

Compile with something like:

clang volume.m -o volume -framework Cocoa -framework AudioToolbox

Then when you run it will simply return a number representing the current effective volume:

0.25

3D printable stencil from binary image

July 21st, 2015

Using code lying around in gptoolbox it was easy to whip up a little function that converts a binary image into a 3d-printable stencil.

hans hass as 3d printed stencil

I’ve added a new function bwstencil.m to gptoolbox. Here’re the simple contents:

  imh = imfill(im,'holes');
  im = imresize(im,2,'nearest');
  [V,F] = bwmesh(~im,varargin{:});
  [V,F] = extrude(V,F);
  [V,F] = remesh_planar_patches(V,F);

Notice that I’m removing “islands” (regions that would be part of the stencil but not connected to the boundary). You can see that the eye of Hans Hass above has been removed in the stencil.

For this example, I’ve 3D printed the result as a proof of concept:

3d printed stencil and result

In retrospect, for this example I should have created the opposite stencil if I’m painting black on white:

hans hass as 3d printed stencil

and the physical proof:

3d printed stencil and result

Wave equation on surfaces

July 21st, 2015

goldfish wave equation

Inspired by a course Misha Kazhdan gave at SGP this year, I started playing around with the wave equation on surfaces. Below is a simple matlab demo:

The core formula here is to update the “height” H_i (extrusion along normal) at vertices for each time step i based on the last two steps:

((1 + b * dt) * M - a * dt^2 *L) H_i = (M * (( 2 - b * dt) * H_{i-1} - H_{i-2}));

b is the dampening factor, dt is the time step size, M is the mass matrix, 1/sqrt(a) is the speed of the wave, L is the cotangent Laplacian, H_i is the unknown heights, H_{i-1} and H_{i-2} are the heights at the previous two time steps. This is discretization with implicit time stepping of the wave equation:

∂²/∂t² u = a∆u - b ∂u/∂t

There’s surely a more canonical place to find this, but I liked the description in “Fast and Exact (Poisson) Solvers on Symmetric Geometries” [Kazhdan 2015].

Save in wave_demo.m

function wave_demo(V,F)
  V = bsxfun(@minus,V,mean(V));
  V = V./max(abs(V(:)));
  V(:,end+1:3) = 0;
  N = per_vertex_normals(V,F);
  L = cotmatrix(V,F);
  M = massmatrix(V,F);
  M = M / max(abs(M(:)));

  dt = 0.0075;
  % 1/sqrt(a) is the speed of the wave
  a = 1;
  H = zeros(size(V,1),1);
  % dampening
  b = 0.05;
  A = ((1+b*dt)*M - a*dt^2*L);
  % Pre-factor
  [cL,p,cS] = chol(A,'lower');
  cU = cL';
  cP = cS';
  cQ = cS;
  chol_solve = @(X) cQ * (cU \ (cL \ ( cP * X)));
  G = H;
  G_m1 = G;      

  [VV,FF,S] = loop([V+bsxfun(@times,G,N) H],F);
  tsh = trisurf(FF,VV(:,1),VV(:,2),VV(:,3),  ...
    'FaceColor','interp','FaceLighting','phong', ...
    'EdgeColor','none','CData',VV(:,4));
  tsh.SpecularStrength = 0.2;
  tsh.SpecularExponent = 100;
  tsh.ButtonDownFcn = @ondown;

  CM = parula(30);colormap(CM(end-11:end,:));
  set(gcf,'Color',[0.3 0.4 1.0]);
  set(gca,'Visible','off');
  l = light('Position',[-0.4 -0.4 1.0],'Style','infinite');
  axis equal;
  axis manual;
  expand_axis(1.2);
  view(2);
  caxis([-max(abs(H)) max(abs(H))]);
  T = get(gca,'tightinset');
  set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

  t = 0;
  draw_t = 0;
  tic;
  while true
    G_m2 = G_m1;
    G_m1 = G;
    G = chol_solve(M*((2-b*dt)*G_m1 - G_m2 + dt*0));

    t = t+dt;
    draw_t = draw_t+toc;
    if draw_t > 0.033
      VV = S*[V+bsxfun(@times,G,N) G];
      tsh.Vertices = VV(:,1:3);
      tsh.CData = VV(:,4);
      title(sprintf('t: %g, max(H): %g',t,max(abs(G))),'FontSize',15);
      drawnow;
      tic;
      draw_t = 0;
    end
  end

  function ondown(src,ev)
    H = a/10*mean(diag(M))*(M\(S'*sparse( ...
      knnsearch(VV(:,1:3),ev.IntersectionPoint,'k',2),1,[1;-1],size(VV,1),1)));
    G = G+H;
  end
end

Libigl wins SGP Software Award

July 20th, 2015

libigl wins SGP software award

At SGP 2015 in Graz, Austria, our C++ library libigl was awarded the SGP Software award. Libigl has been a blast to develop so far, and we’re gaining momentum, both in terms of new users and new contributors. If you’re not using libigl, yet check it out. If you are using libigl, let me know. I’d be very happy to add you or your institution to our list.

Physical simulation “skinning” with biharmonic coordinates

June 25th, 2015

This is a cute little example we didn’t have time to put into our upcoming SIGGRAPH paper Linear Subspace Design for Real-Time Shape Deformation.

octopus upsampling

In this example I first decimate the high resolution (blue) octopus to create the low-res orange octopus. I’m careful to use a decimate that keeps vertices at (or at least near) their original positions. Then I compute “biharmonic coordinates” for vertices of the blue mesh with respect to the vertices of the coarse orange mesh. Notice that the orange mesh does not need to be an enclosing cage. Rather its points just need to lie in or on the blue octopus (there’re even ways to deal with points outside the blue octopus, but right now I’m skipping that). These weights are computed in gptoolbox with:

b = knnsearch(orange_V,blue_tet_V);
W = biharmonic_coordinates(blue_tet_V,blue_tet_T,b);

Then I run a sort of co-rotational linear-elasticity simulation on the orange mesh with collisions activated for the floor. The blue mesh tracks the orange mesh by simply interpolating the orange vertex positions via the coordinates:

blue_V = W * orange_V;

Simple as that. Smooth physics “skinning” with a matrix-matrix multiplication. This is due to the affine precision of our coordinates and their squared Laplacian energy minimization. The “standard” approach to physics “skinning” would be to create the orange tet-mesh so that it envelops the blue mesh and then use piecewise-linear (i.e. non-smooth) interpolation to move the blue mesh vertices.

Update: I should make it super clear that this is not a subspace (aka model reduction) simulation. The simulation of the orange mesh is not aware of the blue mesh at all. Sometimes subspace physics would be more appropriate, though other times pushing the subspace through all parts of the simulation pipeline is impractical/impossible. Hence this prototypical result.

For posterity here’s the entire example code:

% load hi-res mesh
[uV,uF] = load_mesh('~/Dropbox/models/octopus-med.ply');
% load low-res mesh
[~,CV,CF] = qslim(uV,uF,900,'QSlimFlags','-O 0');
[CV,CF] = meshfix(CV,CF);
% Find closest points
b = knnsearch(uV,CV);
CV = uV(b,:);

% Tet-mesh fine resolution mesh and compute weights
[dV,dT,dF] = tetgen(uV,uF,'Flags','-q2');
W = biharmonic_coordinates(dV,dT,b);
W = W(1:size(uV,1),:);

[VV,VT,VF] = tetgen(CV,CF,'Flags','-q100 -Y');
off = mean(VV);
VV=bsxfun(@minus,VV,off);
UU = bsxfun(@plus,VV,[0 0 1]);
vel = zeros(size(UU));
data = [];
M = massmatrix(VV,VT);
% This normalization is very important
M = M/max(M(:));
dt = 1e-2;

t = tsurf(VF,UU, ...
  'FaceColor',[0.8 0.5 0.2], ...
  'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
%t.Visible = 'off';
l = light('Style','local','Position',[5 0 20]);
l2 = light('Style','infinite','Position',[-1 -1 2]);
hold on;
up = @(UU) bsxfun(@plus,[1 0 0],W*UU(1:numel(b),:));
u = tsurf(uF,up(UU),'EdgeColor','none', ...
  'FaceVertexCData',repmat([0.3 0.4 1.0],size(uV,1),1), ...
  fphong,'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
shadows = false;
if shadows
  e = -1e-4;
  QV = [-0.8,-0.5 e; -0.8  0.5 e; 1.5  0.5 e;1.5,-0.5 e];
  QQ = [1 2 3 4];
  trisurf(QQ,QV(:,1),QV(:,2),QV(:,3),'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);
  st = add_shadow(t,l,'Ground',[0 0 -1 0]);
  su = add_shadow(u,l,'Ground',[0 0 -1 0]);
  st.FaceColor = [0.5 0.5 0.5];
  su.FaceColor = [0.5 0.5 0.5];
  set(gca,'Visible','off');
  set(gcf,'Color','w');
end
hold off;
axis equal;
axis([-0.8 1.8 -0.5 0.5 0 1.5]);
axis manual;
camproj('persp');
T = get(gca,'tightinset');
set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

while true
  dur = 0;
  UU0 = UU;
  [UU,data] = arap(VV,VT,[],[], ...
    'Energy','elements', ...
    'V0',UU0, ...
    'Data',data, ...
    'Dynamic',M*repmat([0 0 -9.8],size(VV,1),1), ...
    'MaxIter',100, ...
    'TimeStep',dt, ...
    'Vm1',UU-vel);
  vel = UU-UU0;
  C = UU(:,3)<0;
  UU(C,3) = -UU(C,3);
  vel(C,3) = -vel(C,3);
  dur = dur+dt;
  t.Vertices = UU;
  u.Vertices = up(UU);

  if shadows
    light_pos = [l.Position strcmp(l.Style,'local')];
    ground = [0 0 -1 0];
    d = ground * light_pos';
    shadow_mat = d*eye(4) - light_pos'*ground;
    H = [t.Vertices ones(size(t.Vertices,1),1)]*shadow_mat';
    st.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
    H = [u.Vertices ones(size(u.Vertices,1),1)]*shadow_mat';
    su.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
  end

  drawnow;
end

Hausdorff distance between two triangle meshes

June 23rd, 2015

hausdorff distance

Unless I’m missing something the Hausdorff distance between two triangle meshes will be the maximum over the maximum minimum distance from the vertices of mesh A to the surface of mesh B and the maximum minimum distance from the vertices of mesh B to the surface of mesh A. That is, if p on A and q on B determine the Hausdorff distance between A and B, then it follows that p and/or q is a vertex of A or B respectively. This means you can compute Hausdorff distances easily if you already have a function to compute distances of points (vertices) to a triangle mesh.

So, I’ve added Hausdorff distance calculation between two meshes to my gptoolbox as the function hausdorff. Above I’m drawing the pair of points that determine the Hausdorff distance between the Cheburashka and the knight. I’m also drawing the intersection curves of these two meshes. Here’s how I generate this image:

% compute Hausdorff distance and "determining" pair
[d,pair] = hausdorff(VA,FA,VB,FB);
% Mesh intersection, reindex as a single mesh
[SV,SF,~,J,IM] = selfintersect([VA;VB],[FA;size(VA,1)+FB]);
SF = IM(SF);
% Extract edges along intersection
allE = [SF(:,[2 3]); SF(:,[3 1]); SF(:,[1 2])];
sortallE = sort(allE,2);
[uE,~,IC] = unique(sortallE,'rows');
uE2F = sparse(IC(:),repmat(1:size(SF,1),1,3)',1);
BE = uE(any(uE2F(:,J<=size(FA,1)),2)&any(uE2F(:,J>size(FA,1)),2),:);

% Draw meshes, intersection lines and Hausdorff distance "determiner"
t = tsurf(SF,SV,'CData',1*(J<=size(FA,1)),'EdgeColor','none','FaceAlpha',0.6);
colormap([0.3 0.4 1.0;0.75 0.8 1.0]);axis equal;view(2);
t.SpecularStrength = 0;
t.DiffuseStrength = 0.7;
t.AmbientStrength = 0.3;
hold on;
plot_edges(SV,BE,'LineWidth',2,'Color',[0.7 0.65 0.55]);
p = tsurf([1 2 2],pair,'EdgeColor',[0.9 0.5 0.1],'LineWidth',6);
hold off;

drawnow; 
l = light('Position',[1 -1.5 1.8],'Style','infinite');
ht = add_shadow(t,l);
hp = add_shadow(p,l,'Ground',[0 0 -1 min(SV(:,3))]);
ht.FaceColor = [0.9 0.9 0.9];
hp.FaceColor = ht.FaceColor;
hp.EdgeColor = hp.FaceColor;
hp.LineWidth = p.LineWidth;
apply_ambient_occlusion(t);
set(gcf,'Color','w');
set(gca,'Visible','off');
camproj('persp');
view(-27,6);