s7

import%20marimo%0A%0A__generated_with%20%3D%20%220.19.11%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20rng%20%3D%20np.random.default_rng()%0A%20%20%20%20return%20np%2C%20rng%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20matplotlib%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20plt.ion()%3B%0A%20%20%20%20return%20(plt%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20from%20skimage%20import%20io%2C%20restoration%0A%20%20%20%20import%20scipy.ndimage%20as%20nd%0A%0A%20%20%20%20return%20io%2C%20nd%2C%20restoration%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20scipy%20as%20si%0A%20%20%20%20import%20scipy.optimize%20as%20sio%0A%0A%20%20%20%20return%20si%2C%20sio%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20timeit%0A%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20from%20pathlib%20import%20Path%0A%20%20%20%20mo.pdf(src%3DPath(%22Homework7.pdf%22)%2C%20width%3D%22100%25%22%2C%20height%3D%2250vh%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Gaussian%20MLE%20and%20number%20of%20photons%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20si)%3A%0A%20%20%20%20def%20airypsf(xy%2C%20l%3D0.51%2C%20NA%3D0.9%2C%20side%3D0.7%2C%20N%3D210)%3A%0A%20%20%20%20%20%20%20%20center%20%3D%20xy%0A%20%20%20%20%20%20%20%20xs%2C%20ys%20%3D%20np.meshgrid(np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%2F2%2C%20N)%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%2F2%2C%20N))%0A%20%20%20%20%20%20%20%20r2s%20%3D%20np.square(xs%20-%20center%5B0%5D)%20%2B%20np.square(ys%20-%20center%5B1%5D)%0A%20%20%20%20%20%20%20%20v%20%3D%202%20*%20np.pi%20*%20NA%20*%20np.sqrt(r2s)%20%2F%20l%0A%20%20%20%20%20%20%20%20psf%20%3D%204%20*%20np.square(si.special.j1(v)%20%2F%20v)%0A%20%20%20%20%20%20%20%20nans%20%3D%20np.isnan(psf)%0A%20%20%20%20%20%20%20%20psf%5Bnans%5D%20%3D%201%0A%20%20%20%20%20%20%20%20scale%20%3D%20side%20%2F%20N%0A%20%20%20%20%20%20%20%20return%20psf%20%2F%20np.sum(psf)%2C%20scale%0A%0A%20%20%20%20return%20(airypsf%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20pixelate(image%2C%20inscale%2C%20shape%3D(7%2C%207))%3A%0A%20%20%20%20%20%20%20%20outscale%20%3D%20image.shape%5B0%5D%20*%20inscale%20%2F%20shape%5B0%5D%0A%20%20%20%20%20%20%20%20assert%20image.shape%5B0%5D%20%2F%20image.shape%5B1%5D%20%3D%3D%20shape%5B0%5D%20%2F%20shape%5B1%5D%2C%20%22Aspect%20Ratio%20must%20be%20preserved!%22%0A%0A%20%20%20%20%20%20%20%20out%20%3D%20np.sum(image.reshape(shape%5B0%5D%2C%20int(image.shape%5B0%5D%20%2F%20shape%5B0%5D)%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20shape%5B1%5D%2C%20int(image.shape%5B1%5D%20%2F%20shape%5B1%5D))%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20axis%3D(1%2C%203))%0A%0A%20%20%20%20%20%20%20%20return%20out%2C%20outscale%0A%0A%20%20%20%20return%20(pixelate%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20fishify(image%2C%20N%3D1%2C%20rng%3Dnp.random.default_rng())%3A%0A%20%20%20%20%20%20%20%20return%20rng.poisson(lam%3DN%20*%20(image%20%2F%20np.sum(image)))%0A%0A%20%20%20%20return%20(fishify%2C)%0A%0A%0A%40app.cell%0Adef%20_(airypsf%2C%20fishify%2C%20np%2C%20pixelate)%3A%0A%20%20%20%20def%20simage4(xc%2C%20yc%2C%20N%3D7%2C%20Np%20%3D%201000%2C%20B%3D10)%3A%0A%20%20%20%20%20%20%20%20a%2C%20iscale%20%3D%20airypsf(np.array(%5Bxc%2C%20yc%5D))%0A%20%20%20%20%20%20%20%20p%2C%20scale%20%3D%20pixelate(a%2C%20iscale%2C%20shape%3D(N%2C%20N))%0A%20%20%20%20%20%20%20%20return%20fishify(p%2C%20Np)%20%2B%20fishify(np.ones((N%2C%20N))%2CN%3D(B%20*%20(N%20**%202)))%2C%20scale%0A%0A%20%20%20%20return%20(simage4%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20gausspsf(xc%2C%20yc%2C%20A0%2C%20s%2C%20B%2C%20N%3D210%2C%20side%3D0.7)%3A%0A%20%20%20%20%20%20%20%20xs%20%3D%20np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20N)%0A%20%20%20%20%20%20%20%20ys%20%3D%20np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20N)%0A%20%20%20%20%20%20%20%20xx%2C%20yy%20%3D%20np.meshgrid(xs%2C%20ys)%0A%20%20%20%20%20%20%20%20g%20%3D%20%20A0%20*%20np.exp(-(np.square(xx%20-%20xc)%20%2B%20np.square(yy%20-%20yc))%20%2F%20(2.0%20*%20np.square(s)))%20%2B%20B%0A%20%20%20%20%20%20%20%20return%20g%2C%20(side%20%2F%20N)%0A%0A%20%20%20%20return%20(gausspsf%2C)%0A%0A%0A%40app.cell%0Adef%20_(gausspsf%2C%20np)%3A%0A%20%20%20%20def%20objgauss(params%2C%20im%2C%20side)%3A%0A%20%20%20%20%20%20%20%20xc%2C%20yc%2C%20A0%2C%20s%2C%20B%20%3D%20params%0A%20%20%20%20%20%20%20%20p%2C%20sc%20%3D%20gausspsf(xc%2C%20yc%2C%20A0%2C%20s%2C%20B%2C%20side%3Dside%2C%20N%3Dim.shape%5B0%5D)%0A%20%20%20%20%20%20%20%20%23A%2C%20_%20%3D%20pixelate(p%2C%20sc%2C%20shape%3Dim.shape)%0A%0A%20%20%20%20%20%20%20%20mlogp%20%3D%20np.sum(p%20-%20im%20*%20np.log(p))%0A%20%20%20%20%20%20%20%20%23mlogp%20%3D%20np.sum(A%20-%20im%20*%20np.log(A))%0A%20%20%20%20%20%20%20%20return%20mlogp%0A%0A%20%20%20%20return%20(objgauss%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20objgauss%2C%20sio)%3A%0A%20%20%20%20def%20MLElocalize_gauss(im%2C%20side%3D0.7%2C%20guess%3DNone)%3A%0A%20%20%20%20%20%20%20%20if%20guess%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20guess%20%3D%20np.array(%5B0%2C%200%2C%20np.max(im)%2C%200.1%2C%20np.min(im)%5D)%0A%20%20%20%20%20%20%20%20bnd%20%3D%20((-side%20%2F%202%2C%20side%20%2F%202)%2C%20(-side%20%2F%202%2C%20side%20%2F%202)%2C%20(0%2C%20np.max(im))%2C%20(1e-2%2C%20None)%2C%20(0%2C%20None))%0A%20%20%20%20%20%20%20%20res%20%3D%20sio.minimize(objgauss%2C%20guess%2C%20args%20%3D%20(im%2C%20side)%2C%20bounds%3D%20bnd)%0A%20%20%20%20%20%20%20%20return%20res%0A%0A%20%20%20%20return%20(MLElocalize_gauss%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20a.%0A%20%20%20%20See%20solution%20from%20last%20week!%20Looks%20unbiased%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(MLElocalize_gauss%2C%20np%2C%20plt%2C%20rng%2C%20simage4)%3A%0A%20%20%20%20figr%2C%20axr%20%3D%20plt.subplots()%0A%20%20%20%20xs%20%3D%20rng.uniform(-0.5%20*%200.1%2C%200.5%20*%200.1%2C%20100)%0A%20%20%20%20ys%20%3D%20rng.uniform(-0.5%20*%200.1%2C%200.5%20*%200.1%2C%20100)%0A%20%20%20%20ims%20%3D%20%5Bsimage4(xs%5Bi%5D%2C%20ys%5Bi%5D)%5B0%5D%20for%20i%20in%20range(100)%5D%0A%20%20%20%20es%20%3D%20np.array(%5BMLElocalize_gauss(image).x%5B0%5D%20for%20image%20in%20ims%5D)%0A%20%20%20%20axr.scatter(xs%2C%20es%20-%20xs)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20b.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20RMSE(xs%2C%20ts%3DNone)%3A%0A%20%20%20%20%20%20%20%20if%20ts%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20ts%20%3D%20np.zeros_like(xs)%0A%20%20%20%20%20%20%20%20return%20np.sqrt(np.mean(np.square(xs%20-%20ts)%2C%20axis%3D1))%0A%0A%20%20%20%20return%20(RMSE%2C)%0A%0A%0A%40app.cell%0Adef%20_(MLElocalize_gauss%2C%20np%2C%20simage4)%3A%0A%20%20%20%20Nps%20%3D%20np.logspace(1%2C%205%2C%2030)%0A%20%20%20%20imss%20%3D%20%5Bsimage4(0%2C%200%2C%20Np%3Di)%5B0%5D%20for%20i%20in%20Nps%5D%0A%20%20%20%20ess%20%3D%20np.array(%5BMLElocalize_gauss(image).x%20for%20image%20in%20imss%5D)%5B%3A%2C0%3A2%5D%0A%20%20%20%20return%20Nps%2C%20ess%0A%0A%0A%40app.cell%0Adef%20_(Nps%2C%20RMSE%2C%20ess%2C%20np%2C%20plt)%3A%0A%20%20%20%20fig2%2C%20ax2%20%3D%20plt.subplots()%0A%20%20%20%20ax2.loglog(Nps%2C%20RMSE(ess))%0A%20%20%20%20ax2.loglog(Nps%2C%200.51%20%2F%202%20%2F%200.9%20%2F%20np.sqrt(Nps))%20%23%20theoretical%20baseline%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20TODO%3A%20Need%20to%20RMSE%20over%20multiple%20samples%20for%20each%20Np%2C%20not%20just%20one%20like%20this%2C%20but%20trend%20is%20already%20looking%20on-point.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Radial-symmetry-based%20particle%20localization%0A%20%20%20%20TODO%3A%20Mostly%20running%20Prof.%20code.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Assessing%20Deconvolution%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20gausspsff(xc%2C%20yc%2C%20A0%2C%20s%2C%20B%2C%20N%3D210%2C%20side%3D0.7)%3A%0A%20%20%20%20%20%20%20%20xs%20%3D%20np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20N)%0A%20%20%20%20%20%20%20%20ys%20%3D%20np.linspace(-%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20side%20*%20(1%20-%201%2FN)%20%2F%202%2C%20N)%0A%20%20%20%20%20%20%20%20xx%2C%20yy%20%3D%20np.meshgrid(xs%2C%20ys)%0A%20%20%20%20%20%20%20%20g%20%3D%20%20A0%20*%20np.exp(-(np.square(xx%20-%20xc)%20%2B%20np.square(yy%20-%20yc))%20%2F%20(2.0%20*%20np.square(s)))%20%2B%20B%0A%20%20%20%20%20%20%20%20return%20g%20%2F%20np.sum(g)%2C%20(side%20%2F%20N)%0A%0A%20%20%20%20return%20(gausspsff%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20fishifyf(image%2C%20N%3D1%2C%20rng%3Dnp.random.default_rng())%3A%0A%20%20%20%20%20%20%20%20return%20rng.poisson(lam%3DN%20*%20(image%20%2F%20np.sum(image)))%0A%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gausspsff%2C%20plt)%3A%0A%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots()%0A%20%20%20%20gpsf%20%3D%20gausspsff(0%2C%200%2C%201%2C%207%2C%200%2C%2035%2C%2035)%5B0%5D%0A%20%20%20%20gpsf2%20%3D%20gausspsff(0%2C%200%2C%201%2C%205%2C%200%2C%2035%2C%2035)%5B0%5D%0A%20%20%20%20ax.imshow(gpsf)%0A%20%20%20%20return%20(gpsf%2C)%0A%0A%0A%40app.cell%0Adef%20_(io)%3A%0A%20%20%20%20mouse_glial_cells_crop_png%20%3D%20io.imread(%22mouse_glial_cells_RBurdan_crop.png%22%2C%20as_gray%3DTrue)%0A%20%20%20%20return%20(mouse_glial_cells_crop_png%2C)%0A%0A%0A%40app.cell%0Adef%20_(mouse_glial_cells_crop_png%2C%20plt)%3A%0A%20%20%20%20fig3%2C%20ax3%20%3D%20plt.subplots()%0A%20%20%20%20ax3.imshow(mouse_glial_cells_crop_png%2C%20cmap%3D%22gray%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gpsf%2C%20mouse_glial_cells_crop_png%2C%20nd)%3A%0A%20%20%20%20mouse_glial_cells_crop_conv%20%3D%20nd.convolve(mouse_glial_cells_crop_png%2C%20gpsf%2C%20mode%3D%22constant%22)%0A%20%20%20%20return%20(mouse_glial_cells_crop_conv%2C)%0A%0A%0A%40app.cell%0Adef%20_(mouse_glial_cells_crop_conv%2C%20plt)%3A%0A%20%20%20%20fig4%2C%20ax4%20%3D%20plt.subplots()%0A%20%20%20%20ax4.imshow(mouse_glial_cells_crop_conv%2C%20cmap%3D%22gray%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mouse_glial_cells_crop_conv%2C%20rng)%3A%0A%20%20%20%20mouse_glial_cells_crop_fish%20%3D%20rng.poisson(mouse_glial_cells_crop_conv)%0A%20%20%20%20return%20(mouse_glial_cells_crop_fish%2C)%0A%0A%0A%40app.cell%0Adef%20_(mouse_glial_cells_crop_fish%2C%20plt)%3A%0A%20%20%20%20fig5%2C%20ax5%20%3D%20plt.subplots()%0A%20%20%20%20ax5.imshow(mouse_glial_cells_crop_fish%2C%20cmap%3D%22gray%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gpsf%2C%20mouse_glial_cells_crop_fish%2C%20restoration)%3A%0A%20%20%20%20mouse_glial_cells_crop_restored%20%3D%20restoration.richardson_lucy(mouse_glial_cells_crop_fish%2C%20gpsf%2C%20num_iter%3D20%2C%20clip%3DFalse)%0A%20%20%20%20return%20(mouse_glial_cells_crop_restored%2C)%0A%0A%0A%40app.cell%0Adef%20_(mouse_glial_cells_crop_restored%2C%20plt)%3A%0A%20%20%20%20fig6%2C%20ax6%20%3D%20plt.subplots()%0A%20%20%20%20ax6.imshow(mouse_glial_cells_crop_restored%2C%20cmap%3D%22gray%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20gpsf%2C%0A%20%20%20%20mouse_glial_cells_crop_fish%2C%0A%20%20%20%20mouse_glial_cells_crop_png%2C%0A%20%20%20%20np%2C%0A%20%20%20%20restoration%2C%0A)%3A%0A%20%20%20%20rmss%20%3D%20%5B%5D%0A%20%20%20%20imsi%20%3D%20%5B%5D%0A%20%20%20%20for%20i%20in%20range(1%2C%20250%2C%2010)%3A%0A%20%20%20%20%20%20%20%20mouse_glial_cells_crop_restorei%20%3D%20restoration.richardson_lucy(mouse_glial_cells_crop_fish%2C%20gpsf%2C%20num_iter%3Di%2C%20clip%3DFalse)%0A%20%20%20%20%20%20%20%20mouse_glial_cells_crop_restorei%20%3D%20mouse_glial_cells_crop_restorei%5B35%3A-35%2C%2035%3A-35%5D%0A%20%20%20%20%20%20%20%20mouse_glial_cells_crop_restorei%20-%3D%20np.min(mouse_glial_cells_crop_restorei)%0A%20%20%20%20%20%20%20%20mouse_glial_cells_crop_restorei%20*%3D%20(np.max(mouse_glial_cells_crop_png%5B35%3A-35%2C%2035%3A-35%5D)%20-%20np.min(mouse_glial_cells_crop_png%5B35%3A-35%2C%2035%3A-35%5D))%20%2F%20np.max(mouse_glial_cells_crop_restorei)%0A%20%20%20%20%20%20%20%20mouse_glial_cells_crop_restorei%20%2B%3D%20np.min(mouse_glial_cells_crop_png%5B35%3A-35%2C%2035%3A-35%5D)%0A%20%20%20%20%20%20%20%20imsi.append(mouse_glial_cells_crop_restorei)%0A%20%20%20%20%20%20%20%20rms%20%3D%20np.sqrt(np.mean(np.square(mouse_glial_cells_crop_restorei%20-%20mouse_glial_cells_crop_png%5B35%3A-35%2C%2035%3A-35%5D)))%0A%20%20%20%20%20%20%20%20rmss.append(rms)%0A%20%20%20%20return%20imsi%2C%20rmss%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20rmss)%3A%0A%20%20%20%20fig7%2C%20ax7%20%3D%20plt.subplots()%0A%20%20%20%20ax7.scatter(range(1%2C%20250%2C%2010)%2C%20rmss)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(imsi%2C%20np%2C%20plt%2C%20rmss)%3A%0A%20%20%20%20fig8%2C%20ax8%20%3D%20plt.subplots()%0A%20%20%20%20ax8.imshow(imsi%5Bnp.argmin(rmss)%5D%2C%20cmap%3D%22gray%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A