function id_segments, segments, idx, nmatch=nmatch
s_seg=size(segments)
n_idx=n_elements(idx)
if(s_seg[1] eq 2) then begin
    if(n_idx eq 1) then return,where(segments[0,*] eq idx or segments[1,*] eq idx,nmatch)
    if(n_idx eq 2) then $
        return,where((segments[0,*] eq idx[0] and segments[1,*] eq idx[1]) or $
                     (segments[0,*] eq idx[1] and segments[1,*] eq idx[0]) ,nmatch)
    message,'not supported'
endif else if(s_seg[1] eq 3) then begin
    if(n_idx eq 1) then return,where(segments[0,*] eq idx or $
                                     segments[1,*] eq idx or $
                                     segments[2,*] eq idx,nmatch)
    if(n_idx eq 2) then $
        return,where((segments[0,*] eq idx[0] and segments[1,*] eq idx[1]) or $
                     (segments[0,*] eq idx[1] and segments[1,*] eq idx[0]) or $
                     (segments[0,*] eq idx[0] and segments[2,*] eq idx[1]) or $
                     (segments[0,*] eq idx[1] and segments[2,*] eq idx[0]) or $
                     (segments[1,*] eq idx[0] and segments[2,*] eq idx[1]) or $
                     (segments[1,*] eq idx[1] and segments[2,*] eq idx[0]) $
                     ,nmatch)
    if(n_idx eq 3) then begin
        segments1=segments
        ;for i=0,s_seg[2]-1 do segments1[*,i]=segments[sort(segments[*,i]),i]
        idx1=idx[sort(idx)]
        return,where(segments1[0,*] eq idx1[0] and $
                     segments1[1,*] eq idx1[1] and $
                     segments1[2,*] eq idx1[2],nmatch)
    endif
    message,'not supported'
endif else if(s_seg[1] eq 4) then begin
    if(n_idx eq 4) then begin
        ;for i=0,s_seg[2]-1 do segments1[*,i]=segments[sort(segments[*,i]),i]
        idx1=idx[sort(idx)]
        return,where(segments[0,*] eq idx1[0] and $
                     segments[1,*] eq idx1[1] and $
                     segments[2,*] eq idx1[2] and $
                     segments[3,*] eq idx1[3],nmatch)
    endif
endif else message,'s_seg '+string(s_seg[1])+' is not supported'

end

function concave_hull, v1_0, v2_0, v3_0, $
    dmin=dmin, dmax=dmax, tolerance=tolerance, bounds=bounds, plot=plot

if(n_elements(dmin) ne 1) then dmin=0d
if(n_elements(tolerance) ne 1) then tolerance=1d
n_p_0=n_elements(v1_0)
l_idx_0=lindgen(n_p_0)

if(n_params() eq 3) then begin ;;;; 3d
    sv12=string(v1_0)+' '+string(v2_0)+' '+string(v3_0)
    u_sv1=uniq(sv12,sort(sv12))
    v1=v1_0[u_sv1]
    v2=v2_0[u_sv1]
    v3=v3_0[u_sv1]
    l_idx=l_idx_0[u_sv1]
    n_p=n_elements(v1)
    print,'n_p_0,n_p=',n_p_0,n_p

    qhull, v1, v2, v3, tri, /del, conn=conn_arr
    qhull, v1, v2, v3, cvh
    n_del=n_elements(tri[0,*])
    for i=0L,n_del-1L do tri[*,i]=tri[sort(tri[*,i]),i]
    n_tri=n_elements(cvh[0,*])
    for i=0L,n_tri-1L do cvh[*,i]=cvh[sort(cvh[*,i]),i]

    res=lonarr(3,n_p*10L)-1L
    res[*,0:n_tri-1]=cvh
    i_max=n_tri
    n_btri=1L
    n_isotri=0L
    flag_tri=intarr(n_p*2)-1
    flag_tri[0:n_tri-1L]=0
    n_flagtri=0L
    n_black=0L
    tri_bl=res*0L-1L ;;; blacklisted faces
    n_iter=0L
    while (n_btri gt n_flagtri and i_max lt 4000) do begin
print,'n_btri,n_flagtri,i_max=',n_btri,n_flagtri,i_max
flag_tri[*]=0
        gres=where(res[0,*] gt -1L,cgres)

        ;;;; computing the areas of the faces using the vector cross-product
        a1=v1[res[0,gres]]-v1[res[1,gres]]
        a2=v2[res[0,gres]]-v2[res[1,gres]]
        a3=v3[res[0,gres]]-v3[res[1,gres]]
        b1=v1[res[0,gres]]-v1[res[2,gres]]
        b2=v2[res[0,gres]]-v2[res[2,gres]]
        b3=v3[res[0,gres]]-v3[res[2,gres]]
        c1=b1-a1
        c2=b2-a2
        c3=b3-a3
        d2arr_1=a1^2+a2^2+a3^2
        d2arr_2=b1^2+b2^2+b3^2
        d2arr_3=c1^2+c2^2+c3^2
        d2arr=[d2arr_1,d2arr_2,d2arr_3]
        d2arr_min=min(d2arr,max=d2arr_max,dim=1)
        s2arr=((a2*b3-a3*b2)^2+(a3*b1-a1*b3)^2+(a1*b2-a2*b1)^2)/4d
        res_tmp=lonarr(4,cgres*2L)-1L
        i_tmp=0L
;;;        btri=where(d2arr_min lt dmin^2 or d2arr_max gt dmax^2 and flag_tri eq 0, n_btri)
        btri=where(s2arr lt dmin^2 or s2arr gt dmax^2 and flag_tri eq 0, n_btri)
        for i=0,n_btri-1 do begin
;if i gt 0 then continue
;            print,'possibly bad: ',[res[*,gres[btri[i]]],gres[btri[i]]]
            p1=res[0,gres[btri[i]]]
            p2=res[1,gres[btri[i]]]
            p3=res[2,gres[btri[i]]]
            tri_cur=[p1,p2,p3]
            tri_cur=tri_cur[sort(tri_cur)]
            sarr_cur=norm(crossp([v1[p2]-v1[p1],$
                                 v2[p2]-v2[p1],$
                                 v3[p2]-v3[p1]],$
                                [v1[p3]-v1[p2],$
                                 v2[p3]-v2[p2],$
                                 v3[p3]-v3[p2]]))
            p1_nodes_ori=conn_arr[conn_arr[p1]:conn_arr[p1+1]-1L]
            real_p1_nodes=where(p1_nodes_ori ge 0 and $
                                p1_nodes_ori ne p2 and $
                                p1_nodes_ori ne p3, creal_p1_nodes)
            if(creal_p1_nodes gt 1) then begin
                p1_nodes=p1_nodes_ori[real_p1_nodes]
;            print,'p1_nodes=',p1_nodes
                p1dp1=sqrt((v1[p1_nodes]-v1[p1])^2+(v2[p1_nodes]-v2[p1])^2+(v3[p1_nodes]-v3[p1])^2)
                p1dp2=sqrt((v1[p1_nodes]-v1[p2])^2+(v2[p1_nodes]-v2[p2])^2+(v3[p1_nodes]-v3[p2])^2)
                p1dp3=sqrt((v1[p1_nodes]-v1[p3])^2+(v2[p1_nodes]-v2[p3])^2+(v3[p1_nodes]-v3[p3])^2)
                
                f_p1_nodes=bytarr(creal_p1_nodes)+1
                sarr_p1=dblarr(creal_p1_nodes)
                sarr_p2=dblarr(creal_p1_nodes)
                sarr_p3=dblarr(creal_p1_nodes)
                for kk=0,creal_p1_nodes-1 do begin
                    sarr_p1[kk]=norm(crossp([v1[p1_nodes[kk]]-v1[p1],$
                                        v2[p1_nodes[kk]]-v2[p1],$
                                        v3[p1_nodes[kk]]-v3[p1]],$
                                       [v1[p1_nodes[kk]]-v1[p2],$
                                        v2[p1_nodes[kk]]-v2[p2],$
                                        v3[p1_nodes[kk]]-v3[p2]]))
                    sarr_p2[kk]=norm(crossp([v1[p1_nodes[kk]]-v1[p1],$
                                        v2[p1_nodes[kk]]-v2[p1],$
                                        v3[p1_nodes[kk]]-v3[p1]],$
                                       [v1[p1_nodes[kk]]-v1[p3],$
                                        v2[p1_nodes[kk]]-v2[p3],$
                                        v3[p1_nodes[kk]]-v3[p3]]))
                    sarr_p3[kk]=norm(crossp([v1[p1_nodes[kk]]-v1[p3],$
                                        v2[p1_nodes[kk]]-v2[p3],$
                                        v3[p1_nodes[kk]]-v3[p3]],$
                                       [v1[p1_nodes[kk]]-v1[p2],$
                                        v2[p1_nodes[kk]]-v2[p2],$
                                        v3[p1_nodes[kk]]-v3[p2]]))
                    bhb=where(res eq p1_nodes[kk],cbhb1)
                    bhb=where(res_tmp eq p1_nodes[kk],cbhb2)
                    if(cbhb1 eq 3 or cbhb2 eq 3) then begin ;;;; no points are allowed to be contained in > 3 segments
                        f_p1_nodes[kk]=0
                    endif
                    tet_tmp=[p1_nodes[kk],p1,p2,p3]
                    tet_tmp=tet_tmp[sort(tet_tmp)]
                    bhb=id_segments(tri,tet_tmp,nm=cbhb)
                    if(cbhb ne 1) then f_p1_nodes[kk]=0
                    tri_tmp=[p1_nodes[kk],p1,p2]
                    tri_tmp=tri_tmp[sort(tri_tmp)]
;                    ;;;bhb=id_segments(res[*,0:i_max-1L],tri_tmp,nm=cbhb)
;                    bhb=id_segments(res[*,0:i_max-1L],p1_nodes[kk],nm=cbhb)
;                    if(cbhb gt 0) then f_p1_nodes[kk]=0
                    if(n_black gt 0) then begin ;;; checking for already used ("blacklisted") faces
                        bhb=id_segments(tri_bl[*,0:n_black-1L],tri_tmp,nm=cbhb)
                        if(cbhb gt 0) then f_p1_nodes[kk]=0
                    endif
                endfor
                sarr_p1=transpose(sarr_p1)
                sarr_p2=transpose(sarr_p2)
                sarr_p3=transpose(sarr_p3)
                g_p1_nodes=where(max([p1dp1,p1dp2,p1dp3],dim=1) lt d2arr_max[btri[i]]*tolerance and $
                                 f_p1_nodes eq 1, cg_p1_nodes)
                if(cg_p1_nodes gt 0) then begin
                print,'CG_P1_NODES=',cg_p1_nodes,g_p1_nodes
                    p1_nodes=p1_nodes[g_p1_nodes]
                    p1dp1=p1dp1[g_p1_nodes]
                    p1dp2=p1dp2[g_p1_nodes]
                    p1dp3=p1dp3[g_p1_nodes]
                    sarr_p1=sarr_p1[g_p1_nodes]
                    sarr_p2=sarr_p2[g_p1_nodes]
                    sarr_p3=sarr_p3[g_p1_nodes]
                    ;read,aaaa ;;;;;; cos between the planes should be here:
                    
                    dplane=dblarr(n_elements(p1_nodes))
                    x1vec=[v1[p1],v2[p1],v3[p1]]
                    x2vec=[v1[p2],v2[p2],v3[p2]]
                    x3vec=[v1[p3],v2[p3],v3[p3]]
;                    nvec = crossp((x2vec-x1vec),(x3vec-x1vec))
;                    nvec = nvec/norm(nvec,/double)
;                    for k=0,n_elements(p1_nodes)-1 do begin
;                        x0vec=[v1[p1_nodes[k]],v2[p1_nodes[k]],v3[p1_nodes[k]]]
;                        dplane[k]=abs(total(nvec*(x0vec-x1vec)))
;;                        print,'p1_node,dplane:',p1_nodes[k],dplane[k]
;                    endfor
;                    mindpl=min(dplane,cmdpl) ;;; min angle
                    mindpl=min(sarr_p1+sarr_p2+sarr_p3-sarr_cur,cmdpl)
                    res_tmp[*,i_tmp]=[p1,p1_nodes[cmdpl],p2,gres[btri[i]]]
                    res_tmp[*,i_tmp+1L]=[p1,p1_nodes[cmdpl],p3,i_max]
                    res_tmp[*,i_tmp+2L]=[p1_nodes[cmdpl],p2,p3,i_max+1L]

                    res_tmp[0:2,i_tmp]=res_tmp[sort(res_tmp[0:2,i_tmp]),i_tmp]
                    res_tmp[0:2,i_tmp+1L]=res_tmp[sort(res_tmp[0:2,i_tmp+1L]),i_tmp+1L]
                    res_tmp[0:2,i_tmp+2L]=res_tmp[sort(res_tmp[0:2,i_tmp+2L]),i_tmp+2L]
;                    print, 'new point:',p1_nodes[cmdpl]
;                    print, 'new point:',v1[p1_nodes[cmdpl]],v2[p1_nodes[cmdpl]],v3[p1_nodes[cmdpl]]
;                    print, 'new1:',res_tmp[0:2,i_tmp]
;                    print, 'new2:',res_tmp[0:2,i_tmp+1L]
;                    print, 'new3:',res_tmp[0:2,i_tmp+2L]

                    tri_bl[0:2,n_black:n_black+2L]=res_tmp[0:2,i_tmp:i_tmp+2L] ;; blacklisting
                    i_tmp=i_tmp+3L
                    n_black=n_black+3L
                    ;conn_arr[conn_arr[p1]+real_p1_nodes[g_p1_nodes[cmdpl]]]=-1L
                    flag_tri[btri[i]]=0
                    i_max=i_max+2L
                endif else begin
                    ;;;; other possible segments are too long
                    flag_tri[btri[i]]=1
                    tri_bl[0:2,n_black]=tri_cur
                    n_black=n_black+1L
                endelse
            endif else begin
                ;;;; external point
                flag_tri[btri[i]]=1
                tri_bl[0:2,n_black]=tri_cur
                n_black=n_black+1L
            endelse
        endfor
        nw=where(res_tmp[3,*] gt -1L,cnw)
        if(cnw gt 0) then begin
            res[*,res_tmp[3,nw]]=res_tmp[0:2,nw]
;            print,'added to the res:'
;            print,res_tmp[0:3,nw]
        endif
        bflag=where(flag_tri eq 1,n_flagtri)
        flag_tri[0:i_max-1]=0
;;    if(n_iter gt -1) then break
        n_iter=n_iter+1L
    endwhile

    res=res[*,where(res[0,*] ge 0)]
    n_hull=n_elements(res[0,*])

    res=l_idx[res]
    bounds=res[uniq(res,sort(res))]
    
    if(keyword_set(plot)) then begin
         common managed, wids, widNames, widOuterModal
                                        
         oOrb1 = OBJ_NEW('orb', COLOR=[255, 128, 0]) 
         oOrb1->Scale, .04, .04, .04
         oSymbol1 = OBJ_NEW('IDLgrSymbol', oOrb1) 

         oOrb2 = OBJ_NEW('orb', COLOR=[0, 0, 255]) 
         oOrb2->Scale, .08, .08, .08
         oSymbol2 = OBJ_NEW('IDLgrSymbol', oOrb2) 

         oOrb3 = OBJ_NEW('orb', COLOR=[0, 255, 0]) 
         oOrb3->Scale, .07, .07, .07
         oSymbol3 = OBJ_NEW('IDLgrSymbol', oOrb3) 

         xplot3d,v1_0,v2_0,v3_0,symbol=oSymbol1,linest=6

         wid = KEYWORD_SET(wids) ? LookupManagedWidget('xplot3d') : 0l
         if (wid ne 0l) then begin
             tlb = wid
             widget_control, wid, get_uvalue=pState
             oWindow = (*pState).oWindow
             oScene = (*pState).oScene
             oView2 = (*pState).oView2
             (*pState).oXaxis->GetProperty, RANGE=xaxrange, XCOORD_CONV=xs
             (*pState).oYaxis->GetProperty, RANGE=yaxrange, YCOORD_CONV=ys
             (*pState).oZaxis->GetProperty, RANGE=zaxrange, ZCOORD_CONV=zs
             oModel = OBJ_NEW('IDLgrModel')
             for i=0,n_elements(res[0,*])-1 do begin
                 oPoly=obj_new('IDLgrPolygon',shading=1,$
                     v1_0[res[*,i]],v2_0[res[*,i]],v3_0[res[*,i]],col=[255,0,0],$
                     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs,alpha_channel=0.7)
                 oModel->add,oPoly
             endfor
             oView2->Add, oModel

         endif
                                                
         xplot3d,v1_0[bounds],v2_0[bounds],v3_0[bounds],symbol=oSymbol3,linest=6,/overplot

         bcvh=cvh[uniq(cvh,sort(cvh))]
         xplot3d,v1[bcvh],v2[bcvh],v3[bcvh],symbol=oSymbol2,linest=6,/overplot

         for i=0,n_elements(res[0,*])-1 do $
             xplot3d,[v1_0[res[*,i]],v1_0[res[0,i]]],$
                     [v2_0[res[*,i]],v2_0[res[0,i]]],$
                     [v3_0[res[*,i]],v3_0[res[0,i]]],/overplot,col=[255,255,0],thick=3

         
    endif

    return,res
endif else begin ;;;; 2d
    ;;; stupid check for non-unique points using strings -- weird!!!!
    
    sv12=string(v1_0)+' '+string(v2_0)
    u_sv1=uniq(sv12,sort(sv12))
    v1=v1_0[u_sv1]
    v2=v2_0[u_sv1]
    l_idx=l_idx_0[u_sv1]
    n_p=n_elements(v1)
;    print,'n_p_0,n_p=',n_p_0,n_p

    qhull, v1, v2, tri, /del, conn=conn_arr
    qhull, v1, v2, cvh

    n_seg=n_elements(cvh[0,*])
    for i=0,n_seg-1 do cvh[*,i]=cvh[sort(cvh[*,i]),i]

    res=lonarr(2,n_p)-1L
    res[*,0:n_seg-1]=cvh
    i_max=n_seg
    n_bseg=1L
    n_isoseg=0L
    flag_seg=intarr(n_p)-1
    flag_seg[0:n_seg-1L]=0
    n_flagseg=0L
    while (n_bseg gt n_flagseg) do begin
        gres=where(res[0,*] gt -1L,cgres)
        d2arr=(v1[res[0,gres]]-v1[res[1,gres]])^2+$
              (v2[res[0,gres]]-v2[res[1,gres]])^2
        res_tmp=lonarr(3,cgres*2L)-1L
        i_tmp=0L
        bseg=where(d2arr lt dmin^2 or d2arr gt dmax^2 and flag_seg eq 0, n_bseg)
        for i=0,n_bseg-1 do begin
            p1=res[0,gres[bseg[i]]]
            p2=res[1,gres[bseg[i]]]
            p1_nodes_ori=conn_arr[conn_arr[p1]:conn_arr[p1+1]-1L]
            real_p1_nodes=where(p1_nodes_ori ge 0 and p1_nodes_ori ne p2, creal_p1_nodes)
            p1_nodes=p1_nodes_ori[real_p1_nodes]
            if(n_elements(p1_nodes) gt 1) then begin
                p1da=sqrt((v1[p1_nodes]-v1[p1])^2+(v2[p1_nodes]-v2[p1])^2)
                p1db=sqrt(d2arr[bseg[i]])
                p1dc=sqrt((v1[p1_nodes]-v1[p2])^2+(v2[p1_nodes]-v2[p2])^2)
                f_p1_nodes=bytarr(creal_p1_nodes)+1
                for kk=0,creal_p1_nodes-1 do begin
                    bhb=where(res eq p1_nodes[kk],cbhb)
                    if(cbhb eq 2) then begin ;;;; no points are allowed to be contained in > 2 segments
                        f_p1_nodes[kk]=0
;                        print,'BBB10000',bhb
                    endif
                    bhb=id_segments(res[*,0:i_max-1L],[p1_nodes[kk],p1],nm=cbhb)
                    if(cbhb gt 0) then f_p1_nodes[kk]=0
                    bhb=id_segments(res[*,0:i_max-1L],[p1_nodes[kk],p2],nm=cbhb)
                    if(cbhb gt 0) then f_p1_nodes[kk]=0
                endfor
                g_p1_nodes=where(p1da lt p1db*tolerance and $
                                 p1dc lt p1db*tolerance and $
                                 p1da gt dmin and $
                                 p1dc gt dmin and $
                                 f_p1_nodes eq 1, cg_p1_nodes)
                if(cg_p1_nodes gt 0) then begin
                    p1_nodes=p1_nodes[g_p1_nodes]
                    p1da=p1da[g_p1_nodes]
                    p1db=p1db[g_p1_nodes]
                    p1dc=p1dc[g_p1_nodes]
                    cosab=(p1da^2+d2arr[bseg[i]]-p1dc^2)/(2d*p1da*p1db)
                    mcos=max(cosab,cmcos) ;;; min angle
                    res_tmp[*,i_tmp]=[p1,p1_nodes[cmcos],gres[bseg[i]]]
                    res_tmp[*,i_tmp+1L]=[p1_nodes[cmcos],p2,i_max]
                    res_tmp[0:1,i_tmp]=res_tmp[sort(res_tmp[0:1,i_tmp]),i_tmp]
                    res_tmp[0:1,i_tmp+1L]=res_tmp[sort(res_tmp[0:1,i_tmp+1L]),i_tmp+1L]
                    i_tmp=i_tmp+2L
                    conn_arr[conn_arr[p1]+real_p1_nodes[g_p1_nodes[cmcos]]]=-1L
                    flag_seg[bseg[i]]=0
                    i_max=i_max+1L
                endif else begin
                    ;;;; other possible segments are too long
                    flag_seg[bseg[i]]=1
                endelse
            endif else begin
                ;;;; external point
                flag_seg[bseg[i]]=1
            endelse
;print,res[*,sort(res[0,0:i_max-1L])]
        endfor
        nw=where(res_tmp[2,*] gt -1L,cnw)
        if(cnw gt 0) then $
            res[*,res_tmp[2,nw]]=res_tmp[0:1,nw]
        bflag=where(flag_seg eq 1,n_flagseg)
        flag_seg[0:i_max-1]=0
;        if(keyword_set(plot)) then begin
;            plot,v1,v2,psym=7,syms=0.5,xs=3,ys=3
;            oplot,v1[res[0,0:i_max-1]],v2[res[0,0:i_max-1]],psym=4,syms=2.0,col=254
;            oplot,v1[res[1,0:i_max-1]],v2[res[1,0:i_max-1]],psym=4,syms=2.0,col=254
;            for i=0,i_max-1 do $
;                plots,[v1[res[0,i]],v1[res[1,i]]],$
;                      [v2[res[0,i]],v2[res[1,i]]],col=128
;        endif
;        print,'n_bseg,n_isoseg,i_max,cgres=',n_bseg,n_isoseg,i_max,cgres
;aaa=''
;read,aaa
    endwhile

    res=res[*,where(res[0,*] ge 0)]
    n_hull=n_elements(res[0,*])
    if(keyword_set(plot)) then begin
        plot,v1,v2,psym=7,syms=0.5,xs=3,ys=3,/iso
        oplot,v1[res[0,*]],v2[res[0,*]],psym=4,syms=2.0,col=254
        oplot,v1[res[1,*]],v2[res[1,*]],psym=4,syms=2.0,col=254
        for i=0,n_hull-1 do $
            plots,[v1[res[0,i]],v1[res[1,i]]],[v2[res[0,i]],v2[res[1,i]]],col=128
    endif

    res=l_idx[res]
    bounds=res[uniq(res,sort(res))]
    
    return,res
endelse

end