Find Smallest Circle Fit to Image Blob

// calculate the smallest circle to fit boundary defined by wy, wx
Function calc_SmallestCircle(wy,wx)
    wave wy, wx
    
    variable np, npl, npr
    variable px1, py1, px2, py2, px3, py3
    variable vm1, vm2
    
    // restrict to under 1000 points as needed to avoid long times
    np = numpnts(wx)
    npl = log(np)
    if (npl > 3)
        npl = trunc(npl)
        npr = trunc(np/10^npl)
        npl = (npr*10^(npl - 2))/5
        duplicate/O wx wdx
        duplicate/O wy wdy
        Resample/DOWN=(npl) wdx
        Resample/DOWN=(npl) wdy
        np = numpnts(wdx)
    else
        duplicate/O wx, wdx
        duplicate/O wy, wdy
        np = numpnts(wx)
    endif
 
    // calculate the distances  
    make/N=(np)/FREE/D dmaxarray
    dmaxarray = 0
    
    MultiThread dmaxarray = CalcMaxDistance(wdx, wdy, np, p)
    
    // get three points 
    make/N=3/D/FREE wpx, wpy
    npr = 0.03*np
    WaveStats/Q dmaxarray
    vm1 = v_maxloc
    wpx[0] = wdx[vm1]
    wpy[0] = wdy[vm1]
    npl = (vm1 - npr) < 0 ? 0 : vm1 - npr
    dmaxarray[npl,vm1 + npr] = 0
    WaveStats/Q/R=(inf,vm1) dmaxarray
    vm2 = v_maxloc
    wpx[1] = wdx[vm2]
    wpy[1] = wdy[vm2]
    npl = (vm2 + npr) > (np - 1) ? np - 1 : vm2 + npr
    dmaxarray[vm2 - npr, npl] = 0   
    if ((vm2 - vm1) > (np - vm2)) 
        WaveStats/Q/R=(vm1,vm2) dmaxarray
    else
        WaveStats/Q/R=(vm2,np-1) dmaxarray  
    endif
    wpx[2] = wdx[v_maxloc]
    wpy[2] = wdy[v_maxloc]
    dmaxarray[v_maxloc] = 0
    
    GetCircumCenter(wpx,wpy)
    return 0
end
 
ThreadSafe Function CalcMaxDistance(wx,wy,np,pp)
    wave wx, wy
    variable np, pp
    
    variable rmax
    
    // create 2-D wave  
    make/N=(np)/FREE dd = 0
    // calculate distances and return maximum
    dd = sqrt((wx[p] - wx[pp])^2 + (wy[p] - wy[pp])^2)
    rmax = wavemax(dd)
        
    return rmax
end
 
Function GetCircumCenter(wpx,wpy)
    wave wpx, wpy
    
    variable dlen, cx, cy, rc
    
    dlen = 2*(wpy[0]*wpx[2] + wpy[1]*wpx[0] - wpy[1]*wpx[2] - wpy[0]*wpx[1] - wpy[2]*wpx[0] + wpy[2]*wpx[1])
    
    cx = wpy[1]*wpx[0]^2 - wpy[2]*wpx[0]^2 - (wpy[1]^2)*wpy[0] + (wpy[2]^2)*wpy[0] 
    cx += (wpx[1]^2)*wpy[2] + (wpy[0]^2)*wpy[1] + wpx[2]^2*wpy[0] - (wpy[2]^2)*wpy[1] 
    cx += - (wpx[2]^2)*wpy[1] - (wpx[1]^2)*wpy[0] + (wpy[1]^2)*wpy[2] - (wpy[0]^2)*wpy[2]
    cx /= dlen
    
    cy = (wpx[0]^2)*wpx[2] + (wpy[0]^2)*wpx[2] + (wpx[1]^2)*wpx[0] - (wpx[1]^2)*wpx[2]
    cy += + (wpy[1]^2)*wpx[0] - (wpy[1]^2)*wpx[2] - (wpx[0]^2)*wpx[1] - (wpy[0]^2)*wpx[1]
    cy += - (wpx[2]^2)*wpx[0] + (wpx[2]^2)*wpx[1] - (wpy[2]^2)*wpx[0] + (wpy[2]^2)*wpx[1]
    cy /= dlen
    
    rc = sqrt((wpx[0] - cx)^2 + (wpy[0] - cy)^2)
 
    print "Values (cx, cy, rc) are: ", cx, cy, rc
    
    return 0
end