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 $
        segments1=segments
        for i=0,n_seg-1 do segments1[*,i]=segments[sort(segments[*,i]),i]
        idx1=idx[sort(idx)]
        return,where(segments1[0,*] eq idx[0] and $
                     segments1[1,*] eq idx[1] and $
                     segments1[2,*] eq idx[2],nmatch)
    message,'not supported'
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, bounds=bounds, plot=plot

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

if(n_params() eq 3) then begin ;;;; 3d
    return,!values.f_nan ;;;; not implemented yet

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
        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]
;            print,'p1,p2:',p1,p2
            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*1d and p1dc lt p1db*1d and f_p1_nodes eq 1, cg_p1_nodes)
;;;;                g_p1_nodes=where(p1da lt p1db, 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[*,gres[bseg[i]]]=[p1,p1_nodes[cmcos]]
                    res[*,i_max]=[p1_nodes[cmcos],p2]
                    res[*,gres[bseg[i]]]=res[sort(res[*,gres[bseg[i]]]),gres[bseg[i]]]
                    res[*,i_max]=res[sort(res[*,i_max]),i_max]
                    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
                    flag_seg[bseg[i]]=1
;                    print,'AAAAAAAAAAAAAAA'
                endelse
            endif else begin
                ;;;; external point
;                print,'BBBBBBBBBBBBBB'
                flag_seg[bseg[i]]=1
            endelse
;print,res[*,sort(res[0,0:i_max-1L])]
        endfor
        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