basic course content added, more TODO

courses/uoph410-510a_Image-Analysis/wk7/radialcenter_ImAnClass.py

@@ -0,0 +1,133 @@
+# -*- coding: utf-8 -*-
+# radialcenter_ImAnClass.py
+"""
+Author:   Raghuveer Parthasarathy
+Created on Mon Oct 31 13:33:05 2022
+Last modified on Mon Oct 31 13:33:05 2022
+
+Description
+-----------
+
+Particle localization by radial symmetry
+Python translation of MATLAB radialcenter.m
+
+** Version for Image Analysis Class**
+Same as radialcenter.py , but with sigma and meand2 outputs removed, 
+   and output positions returned relative to image center.
+
+Uses 0 indexing of positions (unlike MATLAB)
+NOTE: *Does not* optimize for image stacks (like radialcenter_stk.m);
+      just single image
+
+Copyright 2011-2022, Raghuveer Parthasarathy, The University of Oregon
+
+Calculates the center of a 2D intensity distribution.
+Method: Considers lines passing through each half-pixel point with slope
+parallel to the gradient of the intensity at that point.  Considers the
+distance of closest approach between these lines and the coordinate
+origin, and determines (analytically) the origin that minimizes the
+weighted sum of these distances-squared.
+Applies simple smoothing if size > 3x3
+
+Inputs
+  I  : 2D intensity distribution (i.e. a grayscale image)
+       Size need not be an odd number of pixels along each dimension
+
+Outputs
+  xc, yc : the center of radial symmetry, px, relative to image center
+           Note that y increases with increasing row number (i.e. "downward")
+
+To do:
+    - Test more (like MATLAB version)
+    - Faster grid creation than meshgrid? (like in MATLAB code)
+    
+see notes August 19-25, Sept. 9, Sept. 19-20 2011
+Raghuveer Parthasarathy
+The University of Oregon
+August 21, 2011 (begun)
+
+"""
+
+import numpy as np   # Will assume numpy is already imported as np !
+
+
+def radialcenter(I):
+    # The main function -- see header comments for details
+    
+    (Ny, Nx) = I.shape
+    # grid coordinates are -n:n, where Nx (or Ny) = 2*n+1
+    # grid midpoint coordinates are -n+0.5:n-0.5
+    xm, ym = np.meshgrid(np.arange(-(Nx-1)/2.0 + 0.5, (Nx-1)/2.0+0.5), 
+                         np.arange(-(Ny-1)/2.0 + 0.5, (Ny-1)/2.0+0.5))
+    # Calculate derivatives along 45-degree shifted coordinates (u and v)
+    # Note that y increases "downward" (increasing row number) -- we'll deal
+    # with this when calculating "m" below.
+    dIdu = I[0:Ny-1,1:]-I[1:,0:Nx-1]
+    dIdv = I[0:Ny-1,0:Nx-1]-I[1:,1:]
+    # Smoothing -- perhaps should be optional
+    fdu = dIdu # will overwrite if smoothing
+    fdv = dIdv
+    if np.min((Nx, Ny))>3:
+        # Only smooth if image is >3px in the smallest dimension
+        # Smooth by simple 3x3 boxcar, which I'll code directly rather than
+        #    calling a convolution.
+        # Zero-pad (expand by 1 on each side)
+        dIdu_pad = np.zeros((Ny+1,Nx+1)) # dIdu array is size Ny-1, Nx-1
+        dIdv_pad = np.zeros((Ny+1,Nx+1)) # dIdv array is size Ny-1, Nx-1
+        dIdu_pad[1:Ny, 1:Nx] = dIdu
+        dIdv_pad[1:Ny, 1:Nx] = dIdv
+        fdu = np.zeros_like(dIdu)
+        fdv = np.zeros_like(dIdv)
+        for j in range(Ny-1):
+            for k in range(Nx-1):
+                fdu[j,k] = np.mean(dIdu_pad[j:j+3,k:k+3])
+                fdv[j,k] = np.mean(dIdv_pad[j:j+3,k:k+3])
+    dImag2 = fdu*fdu + fdv*fdv # gradient magnitude, squared
+    
+    # Slope of the gradient .  Note that we need a 45 degree rotation of 
+    # the u,v components to express the slope in the x-y coordinate system.
+    # The negative sign "flips" the array to account for y increasing
+    # "downward"
+    m = -(fdv + fdu) / (fdu-fdv)
+    # Not smoothed version: m = -(dIdv + dIdu) ./ (dIdu-dIdv)
+    
+    infslope = 9e9 #replace infinite slope values with this extremely large number
+    m[np.isinf(m)] = infslope
+    
+    # Shorthand "b", which also happens to be the
+    # y intercept of the line of slope m that goes through each grid midpoint
+    b = ym - m*xm
+    
+    # Weighting: weight by square of gradient magnitude and inverse 
+    # distance to gradient intensity centroid.
+    sdI2 = np.sum(dImag2)
+    xcentroid = np.sum(dImag2*xm)/sdI2
+    ycentroid = np.sum(dImag2*ym)/sdI2
+    w  = dImag2/np.sqrt((xm-xcentroid)*(xm-xcentroid) + 
+                        (ym-ycentroid)*(ym-ycentroid))
+    
+    # if the intensity is completely flat, m will be NaN (0/0)
+    # give these points zero weight (and set m, b = 0 to avoid 0*NaN=NaN)
+    w[np.isnan(m)]=0
+    b[np.isnan(m)]=0
+    m[np.isnan(m)]=0
+    
+    # least-squares minimization to determine the translated coordinate
+    # system origin (xc, yc) such that lines y = mx+b have
+    # the minimal total distance^2 to the origin:
+    # Unilke the MATLAB version, where I have a separate function
+    #   for this (lsradialcenterfit), I'll just write the calculation here:
+    # Note m, b, w are defined on a grid;  w are the weights for each point
+    wm2p1 = w/(m*m+1)
+    sw  = np.sum(wm2p1)
+    smmw = np.sum(m*m*wm2p1)
+    smw  = np.sum(m*wm2p1)
+    smbw = np.sum(m*b*wm2p1)
+    sbw  = np.sum(b*wm2p1)
+    det = smw*smw - smmw*sw
+    xc = (smbw*sw - smw*sbw)/det # relative to image center
+    yc = (smbw*smw - smmw*sbw)/det # relative to image center
+    
+    
+    return xc, yc
+

courses/uoph410-510a_Image-Analysis/wk7/s7.py

@@ -0,0 +1,362 @@
+import marimo
+
+__generated_with = "0.19.11"
+app = marimo.App(width="medium")
+
+
+@app.cell
+def _():
+    import marimo as mo
+
+    return (mo,)
+
+
+@app.cell
+def _():
+    import numpy as np
+    rng = np.random.default_rng()
+    return np, rng
+
+
+@app.cell
+def _():
+    import matplotlib
+    import matplotlib.pyplot as plt
+    plt.ion();
+    return (plt,)
+
+
+@app.cell
+def _():
+    from skimage import io, restoration
+    import scipy.ndimage as nd
+
+    return io, nd, restoration
+
+
+@app.cell
+def _():
+    import scipy as si
+    import scipy.optimize as sio
+
+    return si, sio
+
+
+@app.cell
+def _():
+    import timeit
+
+    return
+
+
+@app.cell
+def _(mo):
+    from pathlib import Path
+    mo.pdf(src=Path("Homework7.pdf"), width="100%", height="50vh")
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Gaussian MLE and number of photons
+    """)
+    return
+
+
+@app.cell
+def _(np, si):
+    def airypsf(xy, l=0.51, NA=0.9, side=0.7, N=210):
+        center = xy
+        xs, ys = np.meshgrid(np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N)/2, N), 
+                             np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N)/2, N))
+        r2s = np.square(xs - center[0]) + np.square(ys - center[1])
+        v = 2 * np.pi * NA * np.sqrt(r2s) / l
+        psf = 4 * np.square(si.special.j1(v) / v)
+        nans = np.isnan(psf)
+        psf[nans] = 1
+        scale = side / N
+        return psf / np.sum(psf), scale
+
+    return (airypsf,)
+
+
+@app.cell
+def _(np):
+    def pixelate(image, inscale, shape=(7, 7)):
+        outscale = image.shape[0] * inscale / shape[0]
+        assert image.shape[0] / image.shape[1] == shape[0] / shape[1], "Aspect Ratio must be preserved!"
+
+        out = np.sum(image.reshape(shape[0], int(image.shape[0] / shape[0]), 
+                                   shape[1], int(image.shape[1] / shape[1])), 
+                     axis=(1, 3))
+
+        return out, outscale
+
+    return (pixelate,)
+
+
+@app.cell
+def _(np):
+    def fishify(image, N=1, rng=np.random.default_rng()):
+        return rng.poisson(lam=N * (image / np.sum(image)))
+
+    return (fishify,)
+
+
+@app.cell
+def _(airypsf, fishify, np, pixelate):
+    def simage4(xc, yc, N=7, Np = 1000, B=10):
+        a, iscale = airypsf(np.array([xc, yc]))
+        p, scale = pixelate(a, iscale, shape=(N, N))
+        return fishify(p, Np) + fishify(np.ones((N, N)),N=(B * (N ** 2))), scale
+
+    return (simage4,)
+
+
+@app.cell
+def _(np):
+    def gausspsf(xc, yc, A0, s, B, N=210, side=0.7):
+        xs = np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N) / 2, N)
+        ys = np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N) / 2, N)
+        xx, yy = np.meshgrid(xs, ys)
+        g =  A0 * np.exp(-(np.square(xx - xc) + np.square(yy - yc)) / (2.0 * np.square(s))) + B
+        return g, (side / N)
+
+    return (gausspsf,)
+
+
+@app.cell
+def _(gausspsf, np):
+    def objgauss(params, im, side):
+        xc, yc, A0, s, B = params
+        p, sc = gausspsf(xc, yc, A0, s, B, side=side, N=im.shape[0])
+        #A, _ = pixelate(p, sc, shape=im.shape)
+
+        mlogp = np.sum(p - im * np.log(p))
+        #mlogp = np.sum(A - im * np.log(A))
+        return mlogp
+
+    return (objgauss,)
+
+
+@app.cell
+def _(np, objgauss, sio):
+    def MLElocalize_gauss(im, side=0.7, guess=None):
+        if guess is None:
+            guess = np.array([0, 0, np.max(im), 0.1, np.min(im)])
+        bnd = ((-side / 2, side / 2), (-side / 2, side / 2), (0, np.max(im)), (1e-2, None), (0, None))
+        res = sio.minimize(objgauss, guess, args = (im, side), bounds= bnd)
+        return res
+
+    return (MLElocalize_gauss,)
+
+
+@app.cell
+def _():
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ### a.
+    See solution from last week! Looks unbiased:
+    """)
+    return
+
+
+@app.cell
+def _(MLElocalize_gauss, np, plt, rng, simage4):
+    figr, axr = plt.subplots()
+    xs = rng.uniform(-0.5 * 0.1, 0.5 * 0.1, 100)
+    ys = rng.uniform(-0.5 * 0.1, 0.5 * 0.1, 100)
+    ims = [simage4(xs[i], ys[i])[0] for i in range(100)]
+    es = np.array([MLElocalize_gauss(image).x[0] for image in ims])
+    axr.scatter(xs, es - xs)
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ### b.
+    """)
+    return
+
+
+@app.cell
+def _(np):
+    def RMSE(xs, ts=None):
+        if ts is None:
+            ts = np.zeros_like(xs)
+        return np.sqrt(np.mean(np.square(xs - ts), axis=1))
+
+    return (RMSE,)
+
+
+@app.cell
+def _(MLElocalize_gauss, np, simage4):
+    Nps = np.logspace(1, 5, 30)
+    imss = [simage4(0, 0, Np=i)[0] for i in Nps]
+    ess = np.array([MLElocalize_gauss(image).x for image in imss])[:,0:2]
+    return Nps, ess
+
+
+@app.cell
+def _(Nps, RMSE, ess, np, plt):
+    fig2, ax2 = plt.subplots()
+    ax2.loglog(Nps, RMSE(ess))
+    ax2.loglog(Nps, 0.51 / 2 / 0.9 / np.sqrt(Nps)) # theoretical baseline
+    return
+
+
+@app.cell
+def _(mo):
+    mo.md(r"""
+    TODO: Need to RMSE over multiple samples for each Np, not just one like this, but trend is already looking on-point.
+    """)
+    return
+
+
+@app.cell
+def _(mo):
+    mo.md(r"""
+    ## Radial-symmetry-based particle localization
+    TODO: Mostly running Prof. code.
+    """)
+    return
+
+
+@app.cell
+def _(mo):
+    mo.md(r"""
+    ## Assessing Deconvolution
+    """)
+    return
+
+
+@app.cell
+def _(np):
+    def gausspsff(xc, yc, A0, s, B, N=210, side=0.7):
+        xs = np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N) / 2, N)
+        ys = np.linspace(- side * (1 - 1/N) / 2, side * (1 - 1/N) / 2, N)
+        xx, yy = np.meshgrid(xs, ys)
+        g =  A0 * np.exp(-(np.square(xx - xc) + np.square(yy - yc)) / (2.0 * np.square(s))) + B
+        return g / np.sum(g), (side / N)
+
+    return (gausspsff,)
+
+
+@app.cell
+def _(np):
+    def fishifyf(image, N=1, rng=np.random.default_rng()):
+        return rng.poisson(lam=N * (image / np.sum(image)))
+
+    return
+
+
+@app.cell
+def _(gausspsff, plt):
+    fig, ax = plt.subplots()
+    gpsf = gausspsff(0, 0, 1, 7, 0, 35, 35)[0]
+    gpsf2 = gausspsff(0, 0, 1, 5, 0, 35, 35)[0]
+    ax.imshow(gpsf)
+    return (gpsf,)
+
+
+@app.cell
+def _(io):
+    mouse_glial_cells_crop_png = io.imread("mouse_glial_cells_RBurdan_crop.png", as_gray=True)
+    return (mouse_glial_cells_crop_png,)
+
+
+@app.cell
+def _(mouse_glial_cells_crop_png, plt):
+    fig3, ax3 = plt.subplots()
+    ax3.imshow(mouse_glial_cells_crop_png, cmap="gray")
+    return
+
+
+@app.cell
+def _(gpsf, mouse_glial_cells_crop_png, nd):
+    mouse_glial_cells_crop_conv = nd.convolve(mouse_glial_cells_crop_png, gpsf, mode="constant")
+    return (mouse_glial_cells_crop_conv,)
+
+
+@app.cell
+def _(mouse_glial_cells_crop_conv, plt):
+    fig4, ax4 = plt.subplots()
+    ax4.imshow(mouse_glial_cells_crop_conv, cmap="gray")
+    return
+
+
+@app.cell
+def _(mouse_glial_cells_crop_conv, rng):
+    mouse_glial_cells_crop_fish = rng.poisson(mouse_glial_cells_crop_conv)
+    return (mouse_glial_cells_crop_fish,)
+
+
+@app.cell
+def _(mouse_glial_cells_crop_fish, plt):
+    fig5, ax5 = plt.subplots()
+    ax5.imshow(mouse_glial_cells_crop_fish, cmap="gray")
+    return
+
+
+@app.cell
+def _(gpsf, mouse_glial_cells_crop_fish, restoration):
+    mouse_glial_cells_crop_restored = restoration.richardson_lucy(mouse_glial_cells_crop_fish, gpsf, num_iter=20, clip=False)
+    return (mouse_glial_cells_crop_restored,)
+
+
+@app.cell
+def _(mouse_glial_cells_crop_restored, plt):
+    fig6, ax6 = plt.subplots()
+    ax6.imshow(mouse_glial_cells_crop_restored, cmap="gray")
+    return
+
+
+@app.cell
+def _(
+    gpsf,
+    mouse_glial_cells_crop_fish,
+    mouse_glial_cells_crop_png,
+    np,
+    restoration,
+):
+    rmss = []
+    imsi = []
+    for i in range(1, 250, 10):
+        mouse_glial_cells_crop_restorei = restoration.richardson_lucy(mouse_glial_cells_crop_fish, gpsf, num_iter=i, clip=False)
+        mouse_glial_cells_crop_restorei = mouse_glial_cells_crop_restorei[35:-35, 35:-35]
+        mouse_glial_cells_crop_restorei -= np.min(mouse_glial_cells_crop_restorei)
+        mouse_glial_cells_crop_restorei *= (np.max(mouse_glial_cells_crop_png[35:-35, 35:-35]) - np.min(mouse_glial_cells_crop_png[35:-35, 35:-35])) / np.max(mouse_glial_cells_crop_restorei)
+        mouse_glial_cells_crop_restorei += np.min(mouse_glial_cells_crop_png[35:-35, 35:-35])
+        imsi.append(mouse_glial_cells_crop_restorei)
+        rms = np.sqrt(np.mean(np.square(mouse_glial_cells_crop_restorei - mouse_glial_cells_crop_png[35:-35, 35:-35])))
+        rmss.append(rms)
+    return imsi, rmss
+
+
+@app.cell
+def _(plt, rmss):
+    fig7, ax7 = plt.subplots()
+    ax7.scatter(range(1, 250, 10), rmss)
+    return
+
+
+@app.cell
+def _(imsi, np, plt, rmss):
+    fig8, ax8 = plt.subplots()
+    ax8.imshow(imsi[np.argmin(rmss)], cmap="gray")
+    return
+
+
+@app.cell
+def _():
+    return
+
+
+if __name__ == "__main__":
+    app.run()