Archive for July, 2015

Split a long mp3 audio file into 3 min files

Wednesday, 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

Wednesday, 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

Tuesday, 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

Monday, 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

Saturday, 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

Tuesday, 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

Tuesday, 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

Monday, 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.