INTERNAL: How to calculate and visualize a 3-D Voronoi diagram (using QHULL)
Topic:
This techtip is extending 3326 to three dimensions.Discussion:
Because of a bug in QHULL, this does not work exactly as it should...
For more info see CR 30832 and CR 30833.
This example code includes a workaround, and produces a reasonable 3-D Voronoi diagram given random data as input.Solution:
pro test_event, ev
widget_control, ev.id, get_uvalue=uval
widget_control, ev.handler, get_uvalue=model
obj=model->get(position=uval)
obj->getproperty, hide=hide
obj->setproperty, hide=1-hide
xobjview, refresh=ev.top
end
pro test6d,vv,vd,vn
compile_opt idl2
n=7
data=randomn(s,3,n)
qhull,data,tr,connectivity=conn,vvertices=vv,vnormals=vn,vdiagram=vd
n=n_elements(data)/3
;
; find missing connections (CR 30832)
bconn=byte(conn)*0b
i=0
while (i lt n_elements(vd)) do begin
iv1=vd[i+1]
iv2=vd[i+2]
ind=where(conn[conn[iv1]:conn[iv1+1]-1] eq iv2, count)
if count gt 0 then bconn[ind+conn[iv1]]=1b
ind=where(conn[conn[iv2]:conn[iv2+1]-1] eq iv1, count)
if count gt 0 then bconn[ind+conn[iv2]]=1b
i+=vd[i]+1
endwhile
bconn[0:conn[0]-1]=1b
ind=where(bconn eq 0, count)
help, count, ind
vextra=fltarr(4,count/2)
j=0
for i=0,count-1 do begin
if bconn[ind[i]] then continue
iv2=conn[ind[i]]
for iv1=0,n do if conn[iv1] gt ind[i] then break
iv1--
indy=where(conn[conn[iv2]:conn[iv2+1]-1] eq iv1,count)
if count ne 0 then bconn[conn[iv2]+indy]=1b
print, iv1,iv2
vextra[0:2,j]=data[*,iv2]-data[*,iv1]
vextra[3,j]=-total(vextra[0:2,j]*(data[*,iv2]+data[*,iv1]))/2
j++
endfor
vn=[[vn],[vextra]]
;
; bounding box
xmin=min([[data[0,*]],[vv[0,*]]], max=xmax)
ymin=min([[data[1,*]],[vv[1,*]]], max=ymax)
zmin=min([[data[2,*]],[vv[2,*]]], max=zmax)
eps=0.3
xmin=xmin-(xmax-xmin)*eps
xmax=xmax+(xmax-xmin)*eps
ymin=ymin-(ymax-ymin)*eps
ymax=ymax+(ymax-ymin)*eps
zmin=zmin-(zmax-zmin)*eps
zmax=zmax+(zmax-zmin)*eps
;
; define 6 sides of bounding box:
bn=[[1,0,0,-xmin],[0,1,0,-ymin],[0,0,1,-zmin],$
[1,0,0,-xmax],[0,1,0,-ymax],[0,0,1,-zmax]]
; combine all plane equations
an=[[vn],[bn]]
; find all vertices that are intersections between 3 planes
nv=n_elements(an)/4
vx=0
center=([xmax,ymax,zmax]+[xmin,ymin,zmin])/2
range=[xmax,ymax,zmax]-[xmin,ymin,zmin]
for i=2,nv-1 do begin
for j=1,i-1 do begin
for k=0,j-1 do begin
a=an[0:2,[k,j,i]]
ai=invert(a,status)
if (status eq 0) then begin
b=-an[3,[k,j,i]]
x=reform(ai ## b)
if (max(abs(x-center)/range) le 0.501) then begin
if n_elements(vx) eq 1 then vx=x else vx=[[vx],[x]]
endif
endif
endfor
endfor
endfor
; assign to respective data points
vdata=ptrarr(n)
for i=0,n-1 do vdata[i]=ptr_new(-1)
va=[[vv],[vx]]
help,va
nva=n_elements(va)/3
for i=0,nva-1 do begin
; euclidian distances^2
d=(data[0,*]-va[0,i])^2+(data[1,*]-va[1,i])^2+(data[2,*]-va[2,i])^2
ind=where(d-min(d) lt 1e-3,count) ; allow for round off errors
for j=0,count-1 do *vdata[ind[j]]=[*vdata[ind[j]], i]
endfor
for i=0,n-1 do *vdata[i]=(*vdata[i])[1:*]
sym=obj_new('idlgrsymbol',1, size=[.1,.1])
pline=obj_new('idlgrpolyline',va,symbol=sym, linestyle=6, /hide)
; visualize
model=obj_new('idlgrmodel')
model->add, obj_new('idlgraxis',0,range=[-1,2],location=[-1,-1,-1],/hide)
model->add, obj_new('idlgraxis',1,range=[-1,2],location=[-1,-1,-1],/hide)
model->add, obj_new('idlgraxis',2,range=[-1,2],location=[-1,-1,-1],/hide)
model->add, pline
color=byte(randomu(3,3,n)*256)
tc=randomu(1,2,nva)
model->add, obj_new('idlgrtext','.'+strtrim(indgen(n),2),locations=data, $
char_dimensions=[0.2,0.2])
for i=0,n-1 do begin
qhull,va[*,*vdata[i]], tr
conn=[replicate(3,1,n_elements(tr)/3),(*vdata[i])[tr]]
conn=reform(conn, n_elements(conn),/overwrite)
;
alpha=200
im=obj_new('idlgrimage',[[[color[*,i],alpha],[color[*,i],alpha]],$
[[color[*,i],alpha],[color[*,i],alpha]]], $
blend_function=[3,4])
;
model->add, obj_new('idlgrpolygon', va, polygons=conn, $
color=[255,255,255], style=2, $
texture_map=im, texture_coord=tc)
endfor
device, get_screen_size=gss
xobjview, model, tlb=tlb, xsize=gss[0]*0.7, ysize=gss[1]*0.7
base=widget_base(tlb,/row,/nonexclusive, event_pro='test_event',uvalue=model)
for i=0,n-1 do begin
widget_control, widget_button(base,value=strtrim(i,2), uvalue=i+5), $
/set_button
endfor
end