しかし,これらの関数は 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/ からダウンロードできます.他の用途への改造も比較的容易だと思うので,よければ参考にしてください.
