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