MATLABでQhullを使う

Qhull （http://www.qhull.org/） は計算幾何に関する有用なライブラリーで，MATLABの幾つかの関数，delaunayn や voronoin は Qhull を利用しています．

しかし，これらの関数は Qhull が持つすべての機能を実装しているわけではなく，例えばドロネー分割における法線や隣接単体の情報を出力するオプションはMATLABでの実装ではサポートされていません．

本来なら Qhull を呼び出すmexファイルを書くのが正攻法ですが，面倒なのでコンパイル済みのQhullの実行ファイルを呼び出す関数を作ってみました．

僕が欲しかったのはドロネー単体分割時の隣接単体の情報だったので，以下の関数はドロネー単体の情報に加えて単体間の隣接行列を返す仕様になっています．

function [tes,adj]=mydelaunayn(X)
%MYDELAUNAYN N-D Delaunay triangulation with adjacency matrix
%   [tes,adj]=mydelaunayn(X)
%   tes: sparse matrix which represents tessellation
%      tes(i,j)==1 if j-th simplex has vertex X(i,:)
%   adj: adjacency matrix among Delaunay simplices
%      adj(i,j)==1 if i-th and j-th simplices has common hyper plane
%
%   NOTE:
%     This routine requires excutable of Qhull.
%     http://www.qhull.org/


%% Compose input
infile=tempname;
fid=fopen(infile,'w');
d=size(X,2);
N=size(X,1);
fprintf(fid,'%d\n',d);
fprintf(fid,'%d\n',N);
fprintf(fid,strcat('%e',repmat('\t%e',1,d-1),'\n'),X');
fclose(fid);

%% Invoke Qhull
outfile=tempname;
opts={'Qt','Qbb','Qc'};
if d>3
    opts{end+1}='Qx';
end
opts=strcat(opts,{' '});
system(sprintf('qhull d s i %s Fn < %s > %s',cat(2,opts{:}),infile,outfile));

%% Read Results
fid=fopen(outfile,'r');

% compose tessellation matrix
Ntes=fscanf(fid,'%d',1);
tesmat=fscanf(fid,'%d',[d+1,Ntes])'+1;
idx=1:Ntes;
tes=sparse(reshape(tesmat',[],1),reshape(idx(ones(d+1,1),:),[],1),1,N,Ntes);

% compose adjancy matrix
Nadj=fscanf(fid,'%d',1);
adjout=fscanf(fid,'%d',[d+2,Nadj])';
assert(all(adjout(:,1)==d+1));
assert(Nadj==Ntes);
adjmat=[reshape(idx(ones(d+1,1),:),[],1),reshape(adjout(:,2:end)',[],1)];
adjmat(adjmat(:,2)<0,:)=[];
adj=sparse(adjmat(:,1),adjmat(:,2)+1,1,Ntes,Ntes);

%% Cleanup
delete(infile);
delete(outfile);
end

上の関数の実行にはパスに Qhull の実行ファイルがあることが必要で，これは http://www.qhull.org/ からダウンロードできます．他の用途への改造も比較的容易だと思うので，よければ参考にしてください．