corner detection notebook

import%20marimo%0A%0A__generated_with%20%3D%20%220.23.2%22%0Aapp%20%3D%20marimo.App(width%3D%22full%22)%0A%0A%0A%40app.cell%0Adef%20imports()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20scipy.ndimage%20import%20gaussian_filter%2C%20map_coordinates%0A%0A%20%20%20%20return%20gaussian_filter%2C%20map_coordinates%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell%0Adef%20intro(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Camera%20calibration%20part%203%20of%20%3F%3A%20But%20where%20are%20the%20features%3F%0A%0A%0A%20%20%20%20%23%23%20Pseudo%20philosophy%0A%0A%20%20%20%20We're%20still%20on%20our%20quest%20to%20turn%20an%20artifact%20of%20the%20photo-electric%20effect%20into%20something%20we%20can%20comfortably%20call%20%22a%20measurement%22.%20Philosophically%20these%20are%20two%20things%20which%20are%20connected%2C%20but%20are%20not%20the%20same.%20The%20former%20is%20just%20a%20matter-of-fact%20phenomenon%20which%20happens%20when%20light%20falls%2C%20in%20a%20meticulously%20engineered%20and%20precogitated%20fashion%2C%20on%20a%20sensor%20and%20is%20converted%20to%20electrical%20charge.%20It%20could%20equally%20have%20been%20a%20chemically%20activated%20photographic%20plate%2C%20but%20we%20are%20dealing%20with%20a%20digital%20image%20here.%20This%20artifact%20is%20a%20recording%20of%20the%20scene%20from%20which%20the%20light%20came.%20But%20it%20is%20not%20something%20we%20can%20directly%20observe%20to%20conclude%20anything%20about%20distances%20and%20angles%20in%20that%20scene.%20It's%20also%20not%20something%20we%20can%20immediately%20use%20to%20conclude%20anything%20about%20material%20properties%20of%20objects%20in%20the%20scene.%20In%20fact%2C%20any%20conclusion%20we%20reach%20about%20the%20world%20would%20require%20we%20use%20some%20modelling%20to%20connect%20the%20%22real%20world%22%20to%20the%20recording%20we%20made.%20It's%20not%20unlike%20the%20shadows%20on%20the%20side%20of%20the%20cave%20from%20%5BPlato's%20Allegory%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FAllegory_of_the_cave).%0A%0A%20%20%20%20To%20conclude%20our%20modelling%20we%20need%20to%20get%20data%20from%20the%20image.%20In%20a%20previous%20post%2C%20we've%20already%20refered%20to%20such%20data%2C%20namely%3A%20the%20coordinates%20of%20point%20on%20the%20image%20of%20a%20checkerboard%20grid.%20But%20we've%20not%20said%20anything%20about%20HOW%20we%20get%20them%20from%20the%20image.%20Surely%2C%20there%20has%20to%20be%20some%20procedure%20from%20which%20we%20identify%20a%20corner%20between%20checkerboard%20squares%20and%20extract%20the%20coordinates%20on%20the%20sensor.%20Clearly%2C%20the%20precision%20with%20which%20we%20can%20get%20these%20coordinates%20is%20key%20for%20the%20precision%20with%20which%20we%20will%20calibrate%20our%20camera.%0A%0A%0A%20%20%20%20%23%23%20Precision%20in%20sampling%0A%0A%20%20%20%20Our%20accuracy%20depends%20on%20multiple%20system%20parameters.%20For%20example%2C%20the%20pixel%20size%20of%20the%20CCD%20we%20sample%20the%20image.%20Different%20properties%20of%20the%20CCD%20determine%20many%20other%20things%2C%20but%20for%20this%20simplified%20discussion%20we%20will%20limit%20ourselves%20to%20the%20how%20much%20information%20we%20get%20in%20SPACE%20and%20only%20due%20to%20the%20pixel%20size.%20If%20we%20have%20a%20pixel%20that's%2010%20microns%20across%2C%20then%20we%20will%20only%20be%20able%20to%20get%20a%20data%20point%20on%20a%2010X10%20micron%20grid.%20But%20is%20that%20the%20precision%20with%20which%20we%20measure%20the%20world%3F%20no!%20To%20get%20a%20handle%20on%20that%20we%20need%20to%20bring%20in%20the%20optics%20and%20see%20how%20the%20world%20projects%20onto%20the%20sensor.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20plate_scale_controls(mo)%3A%0A%20%20%20%20k1_sl%20%20%20%3D%20mo.ui.slider(-0.50%2C%200.50%2C%20step%3D0.05%2C%20value%3D-0.15%2C%20label%3D%22k%E2%82%81%20%20(dimensionless)%22)%0A%20%20%20%20pf_sl%20%20%20%3D%20mo.ui.slider(0.5%2C%20%20%205.0%2C%20%20step%3D0.25%2C%20value%3D2.0%2C%20%20%20label%3D%22p%2Ff%20%20(pixel%20angular%20size%2C%20%20mrad%2Fpixel)%20%20%5B2%20mrad%20%E2%89%88%2010%20%C2%B5m%20pixel%2C%205%20mm%20EFL%5D%22)%0A%20%20%20%20tilt_sl%20%3D%20mo.ui.slider(0%2C%20%20%20%20%2060%2C%20%20%20step%3D5%2C%20%20%20%20value%3D0%2C%20%20%20%20%20label%3D%22board%20tilt%20%20(%C2%B0)%22)%0A%20%20%20%20return%20k1_sl%2C%20pf_sl%2C%20tilt_sl%0A%0A%0A%40app.cell%0Adef%20checkerboard_functions(np)%3A%0A%20%20%20%20BOARD_ROWS%20%3D%206%0A%20%20%20%20BOARD_COLS%20%3D%209%0A%20%20%20%20SQUARE_SIZE%20%3D%2025.0%0A%20%20%20%20PATTERN_SIZE%20%3D%20(BOARD_COLS%20-%201%2C%20BOARD_ROWS%20-%201)%20%20%23%20inner%20corners%3A%20(8%2C%205)%0A%0A%20%20%20%20def%20generate_checkerboard(rows%3DBOARD_ROWS%2C%20cols%3DBOARD_COLS%2C%20square_px%3D60)%3A%0A%20%20%20%20%20%20%20%20h%20%3D%20rows%20*%20square_px%0A%20%20%20%20%20%20%20%20w%20%3D%20cols%20*%20square_px%0A%20%20%20%20%20%20%20%20img%20%3D%20np.zeros((h%2C%20w)%2C%20dtype%3Dnp.uint8)%0A%20%20%20%20%20%20%20%20for%20r%20in%20range(rows)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20c%20in%20range(cols)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20(r%20%2B%20c)%20%25%202%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20img%5Br%20*%20square_px%3A(r%20%2B%201)%20*%20square_px%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%20*%20square_px%3A(c%20%2B%201)%20*%20square_px%5D%20%3D%20255%0A%20%20%20%20%20%20%20%20return%20img%0A%0A%20%20%20%20def%20get_checkerboard_3d_points(rows%3DBOARD_ROWS%2C%20cols%3DBOARD_COLS%2C%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%20square_size%3DSQUARE_SIZE)%3A%0A%20%20%20%20%20%20%20%20pattern%20%3D%20(cols%20-%201%2C%20rows%20-%201)%0A%20%20%20%20%20%20%20%20objp%20%3D%20np.zeros((pattern%5B0%5D%20*%20pattern%5B1%5D%2C%203)%2C%20dtype%3Dnp.float64)%0A%20%20%20%20%20%20%20%20objp%5B%3A%2C%20%3A2%5D%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20np.mgrid%5B0%3Apattern%5B0%5D%2C%200%3Apattern%5B1%5D%5D.T.reshape(-1%2C%202)%20*%20square_size%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20objp%0A%0A%20%20%20%20ideal_board%20%3D%20generate_checkerboard()%0A%20%20%20%20obj_points_3d%20%3D%20get_checkerboard_3d_points()%0A%20%20%20%20return%20PATTERN_SIZE%2C%20ideal_board%2C%20obj_points_3d%0A%0A%0A%40app.cell%0Adef%20plate_scale_plot(k1_sl%2C%20mo%2C%20np%2C%20pf_sl%2C%20plt%2C%20tilt_sl)%3A%0A%20%20%20%20_expl%20%3D%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20What%20is%20the%20position%20uncertainty%20for%20integer-pixel%20coordinates%3F%0A%0A%20%20%20%20The%20pixel%20pitch%20sets%20how%20finely%20the%20sensor%20samples%20the%20image.%20What%20if%20we%20work%20with%20this%20and%20only%20report%20a%20corner's%20position%20at%20integer%20pixel%20coordinate%3F%20The%20true%20location%20could%20be%20anywhere%20within%20that%20pixel%20%E2%80%94%20an%20uncertainty%0A%20%20%20%20of%20%C2%B10.5%20px%20in%20each%20direction.%20How%20bad%20is%20that%3F%20%0A%0A%20%20%20%20The%20angular%20footprint%20of%20a%20pixel%20that's%20on%20the%20optical%20axis%20in%20world%20space%20is%20given%20by%3A%0A%0A%20%20%20%20%24%24s_0%20%3D%20%5Cfrac%7Bp%20%5Ccdot%20%5Ctext%7BWD%7D%7D%7Bf%7D%20%5Cquad%20%5Ctext%7B%5Bmm%2Fpixel%5D%7D%24%24%0A%0A%20%20%20%20This%20gives%20a%20pixel-grid%20resolution%20limit%20of%20**%C2%B1%24s_0%2F2%24**%20mm%20at%20the%20optical%20axis.%0A%0A%20%20%20%20**Two%20effects%20make%20this%20worse%20away%20from%20the%20centre%3A**%0A%0A%20%20%20%201.%20Lens%20distortion%20-%20the%20local%20mapping%20stretches%20due%20to%20distortion.%20For%20simplicity%20let's%20select%20a%20distortion%20model%20with%20just%20one%20coefficient%20%24k_1%24%20with%20normalised%20radius%20%24r_n%24%3A%0A%20%20%20%20%20%20%20%24s_%7B%5Ctext%7Beff%7D%7D(r)%20%3D%20s_0%20%2F%20(1%20%2B%203k_1%20r_n%5E2)%24.%0A%20%20%20%20%20%20%20Barrel%20distortion%20(%24k_1%20%3C%200%24)%20means%20edge%20pixels%20cover%20more%20object%20area.%0A%0A%20%20%20%202.%20Board%20tilt%20%E2%80%94%20when%20the%20calibration%20target%20is%20tilted%20by%20angle%20%24%5Ctheta%24%2C%0A%20%20%20%20%20%20%20different%20parts%20of%20the%20board%20sit%20at%20different%20depths.%20A%20ray%20through%20pixel%0A%20%20%20%20%20%20%20%24(u%2C%20v)%24%20has%20normalised%20vertical%20coordinate%20%24y_n%20%3D%20(v%20-%20c_y)%20%2F%20f_x%24.%0A%20%20%20%20%20%20%20Intersecting%20this%20ray%20with%20a%20plane%20tilted%20%24%5Ctheta%24%20around%20the%20horizontal%0A%20%20%20%20%20%20%20axis%20at%20distance%20WD%20gives%20a%20local%20depth%3A%20%24Z_%7B%5Ctext%7Blocal%7D%7D%20%3D%20%5Cfrac%7B%5Ctext%7BWD%7D%20%5Ccdot%20%5Ccos%5Ctheta%7D%7B%5Csin%5Ctheta%20%5Ccdot%20y_n%20%2B%20%5Ccos%5Ctheta%7D%24.%20The%20plate%20scale%20at%20each%20pixel%20scales%20proportionally%20to%20this%20depth%20(farther%20%3D%20more%20mm%20per%20pixel)%2C%20so%3A%20%24s_%7B%5Ctext%7Btilt%7D%7D(y_n)%20%3D%20s_0%20%5Ccdot%20%5Cfrac%7B%5Ccos%5Ctheta%7D%7B%5Csin%5Ctheta%20%5Ccdot%20y_n%20%2B%20%5Ccos%5Ctheta%7D%24.%20At%20%24%5Ctheta%20%3D%200%24%20this%20reduces%20to%20%24s_0%24%20everywhere.%20At%20%24%5Ctheta%20%3D%2030%C2%B0%24%20the%20top%20and%20bottom%20of%20the%20sensor%20can%20differ%20by%20tens%20of%20percent%20%E2%80%94%20visible%20as%20vertical%20asymmetry%20in%20the%20heatmap.%0A%0A%20%20%20%20We%20draw%20a%20heat%20map%20showing%20how%20this%20scale%20changes%20across%20the%20FOV.%20This%20gives%20us%20the%20position%20uncertainty.%20%0A%0A%20%20%20%20Note%20that%20areas%20of%20the%20map%20that%20are%20white%20(blank)%20indicate%20where%20the%20%24k_1%20r%5E2%24%20distortion%20mode%20breaks%20down.%20The%20radial%20distortion%20expression%20%241%20%2B%203k_1%20r_n%5E2%24%20passes%20through%20zero%20at%20%24r_n%20%3D%201%2F%5Csqrt%7B3%7Ck_1%7C%7D%24%2C%20beyond%20which%20the%20polynomial%20predicts%20the%20mapping%20folds%0A%20%20%20%20back%20on%20itself%20%E2%80%94%20two%20different%20world%20points%20projecting%20to%20the%20same%20pixel.%0A%20%20%20%20Because%20this%20is%20not%20physical%2C%20we%20mask%20these%20regions.%20%0A%0A%20%20%20%20What%20we%20see%20by%20playing%20around%20with%20the%20sliders%20is%20that%20we%20typically%20get%20a%20position%20uncertainty%20of%20hundreds%20of%20microns.%20Is%20this%20good%3F%20bad%3F%20Depends%20on%20what%20you%20want%20to%20do%2C%20but%20if%20you%20want%20for%20example%20to%20use%20machine%20vision%20for%20assembly%20inspection%20with%20tight%20tollerances%2C%20that's%20not%20good%20enough.%0A%0A%20%20%20%20%22%22%22)%0A%0A%0A%20%20%20%20_k1%20%20%20%20%20%20%20%3D%20float(k1_sl.value)%0A%20%20%20%20_pf_mrad%20%20%3D%20float(pf_sl.value)%0A%20%20%20%20_tilt_deg%20%3D%20float(tilt_sl.value)%0A%20%20%20%20_tilt_rad%20%3D%20np.deg2rad(_tilt_deg)%0A%0A%20%20%20%20_WD%20%20%3D%20200.0%0A%20%20%20%20_pf%20%20%3D%20_pf_mrad%20*%201e-3%0A%20%20%20%20_fx%20%20%3D%201.0%20%2F%20_pf%0A%20%20%20%20_s0%20%20%3D%20_WD%20*%20_pf%0A%20%20%20%20_W%2C%20_H%20%20%20%3D%20640%2C%20480%0A%20%20%20%20_cx%2C%20_cy%20%3D%20_W%20%2F%202.0%2C%20_H%20%2F%202.0%0A%20%20%20%20_u_g%2C%20_v_g%20%3D%20np.meshgrid(np.arange(_W%2C%20dtype%3Dfloat)%2C%20np.arange(_H%2C%20dtype%3Dfloat))%0A%20%20%20%20_xn%20%20%3D%20(_u_g%20-%20_cx)%20%2F%20_fx%0A%20%20%20%20_yn%20%20%3D%20(_v_g%20-%20_cy)%20%2F%20_fx%0A%20%20%20%20_rn2%20%3D%20_xn%20**%202%20%2B%20_yn%20**%202%0A%0A%20%20%20%20%23%20Distortion%20Jacobian%20and%20tilt%20denominator%20%E2%80%94%20no%20clipping%2C%20NaN%20for%20invalid%0A%20%20%20%20_jac%20%20%20%20%20%20%20%20%3D%201.0%20%2B%203.0%20*%20_k1%20*%20_rn2%0A%20%20%20%20_cos_t%20%20%20%20%20%20%3D%20np.cos(_tilt_rad)%0A%20%20%20%20_sin_t%20%20%20%20%20%20%3D%20np.sin(_tilt_rad)%0A%20%20%20%20_tilt_denom%20%3D%20_sin_t%20*%20_yn%20%2B%20_cos_t%0A%0A%20%20%20%20%23%20Valid%20where%20both%20denominators%20are%20positive%20(model%20is%20monotonic%20%26%20in%20front%20of%20camera)%0A%20%20%20%20_valid%20%3D%20(_jac%20%3E%200.01)%20%26%20(_tilt_denom%20%3E%200.01)%0A%20%20%20%20_VMAX%20%3D%200.8%0A%20%20%20%20_uncert_raw%20%3D%20np.where(_valid%2C%200.5%20*%20_s0%20*%20_cos_t%20%2F%20(_jac%20*%20_tilt_denom)%2C%20np.nan)%0A%20%20%20%20_uncert%20%3D%20np.where(_uncert_raw%20%3C%3D%20_VMAX%2C%20_uncert_raw%2C%20np.nan)%0A%0A%20%20%20%20_r_corner%20%3D%20float(np.sqrt(_cx%20**%202%20%2B%20_cy%20**%202))%0A%0A%20%20%20%20%23%20Scalar%20summaries%0A%20%20%20%20_rn2_c%20%20%20%20%3D%20(_cx%20**%202%20%2B%20_cy%20**%202)%20%2F%20_fx%20**%202%0A%20%20%20%20_jac_c%20%20%20%20%3D%201.0%20%2B%203.0%20*%20_k1%20*%20_rn2_c%0A%20%20%20%20_u_ctr%20%20%20%20%3D%200.5%20*%20_s0%0A%20%20%20%20_u_cor%20%20%20%20%3D%200.5%20*%20_s0%20%2F%20_jac_c%20if%20_jac_c%20%3E%200.01%20else%20float('nan')%0A%20%20%20%20_pct%20%20%20%20%20%20%3D%20(1.0%20%2F%20_jac_c%20-%201.0)%20*%20100.0%20if%20_jac_c%20%3E%200.01%20else%20float('nan')%0A%20%20%20%20_fov_deg%20%20%3D%20float(np.degrees(np.arctan(_r_corner%20%2F%20_fx)))%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Figure%20%E2%94%80%E2%94%80%0A%20%20%20%20_fig%2C%20(_ax1%2C%20_ax2)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205.5))%0A%0A%20%20%20%20%23%20Left%3A%202-D%20heatmap%0A%20%20%20%20_im%20%20%20%3D%20_ax1.imshow(_uncert%2C%20cmap%3D%22YlOrRd%22%2C%20origin%3D%22upper%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20extent%3D%5B-_cx%2C%20_cx%2C%20_cy%2C%20-_cy%5D%2C%20vmin%3D0%2C%20vmax%3D_VMAX)%0A%20%20%20%20plt.colorbar(_im%2C%20ax%3D_ax1%2C%20label%3D%22Pixel-grid%20resolution%20limit%20(mm)%22)%0A%20%20%20%20_ax1.plot(0%2C%200%2C%20%22b%2B%22%2C%20markersize%3D14%2C%20markeredgewidth%3D2.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22centre%20%20%C2%B1%7B_u_ctr%3A.3f%7D%20mm%22)%0A%20%20%20%20_ax1.text(-_cx%20%2B%204%2C%20_cy%20-%206%2C%20f%22corner%20%20%C2%B1%7B_u_cor%3A.3f%7D%20mm%20(%7B_pct%3A%2B.1f%7D%25)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20color%3D%22navy%22%2C%20fontsize%3D9%2C%20fontweight%3D%22bold%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20bbox%3Ddict(facecolor%3D%22white%22%2C%20alpha%3D0.75%2C%20boxstyle%3D%22round%2Cpad%3D0.2%22))%0A%20%20%20%20_ax1.legend(fontsize%3D9%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_ax1.set_xlabel(%22u%20%E2%88%92%20centre%20(pixels)%22)%0A%20%20%20%20_ax1.set_ylabel(%22v%20%E2%88%92%20centre%20(pixels)%22)%0A%20%20%20%20_ax1.set_title(%22Pixel-grid%20resolution%20limit%5Cnin%20object%20space%20(mm)%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%2B%20(f%22%20%20%5Btilt%20%3D%20%7B_tilt_deg%3A.0f%7D%C2%B0%5D%22%20if%20_tilt_deg%20%3E%200%20else%20%22%22))%0A%20%20%20%20_ax1.axvline(0%2C%20color%3D%22steelblue%22%2C%20lw%3D1%2C%20ls%3D%22--%22%2C%20alpha%3D0.6%2C%20label%3D%22cut%20%E2%86%92%22)%0A%0A%20%20%20%20%23%20Right%3A%20vertical%20cross-section%20through%20the%20centre%20column%20(u%20%3D%20cx)%0A%20%20%20%20_cut_v%20%3D%20_uncert%5B%3A%2C%20int(_cx)%5D%0A%20%20%20%20_v_centered%20%3D%20np.arange(_H%2C%20dtype%3Dfloat)%20-%20_cy%20%20%23%200%20at%20centre%0A%20%20%20%20_ax2.fill_between(_v_centered%2C%200%2C%20_cut_v%2C%20alpha%3D0.3%2C%20color%3D%22crimson%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20where%3D~np.isnan(_cut_v))%0A%20%20%20%20_ax2.plot(_v_centered%2C%20_cut_v%2C%20color%3D%22crimson%22%2C%20lw%3D2.5)%0A%20%20%20%20_ax2.axhline(_u_ctr%2C%20color%3D%22steelblue%22%2C%20lw%3D1%2C%20ls%3D%22%3A%22%2C%20label%3Df%22on-axis%20s%E2%82%80%2F2%20%3D%20%7B_u_ctr%3A.3f%7D%20mm%22)%0A%20%20%20%20_ax2.axvline(0%2C%20color%3D%22gray%22%2C%20lw%3D0.5%2C%20alpha%3D0.4)%0A%20%20%20%20_ax2.set_xlabel(%22v%20%E2%88%92%20centre%20(pixels)%22)%0A%20%20%20%20_ax2.set_ylabel(%22Pixel-grid%20resolution%20limit%20(mm)%22)%0A%20%20%20%20_ax2.set_title(%22Vertical%20cut%20through%20centre%20column%22)%0A%20%20%20%20_ax2.legend(fontsize%3D9)%0A%20%20%20%20_ax2.set_xlim(-_cy%2C%20_cy)%0A%20%20%20%20_ax2.set_ylim(0%2C%20_VMAX)%0A%20%20%20%20_ax2.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Distortion%20description%20%E2%94%80%E2%94%80%0A%20%20%20%20if%20abs(_k1)%20%3C%200.01%3A%0A%20%20%20%20%20%20%20%20_dist_desc%20%3D%20%22no%20distortion%20%E2%80%94%20plate%20scale%20is%20flat%20across%20the%20field%22%0A%20%20%20%20elif%20_k1%20%3C%200%3A%0A%20%20%20%20%20%20%20%20_dist_desc%20%3D%20(f%22barrel%20distortion%20(k%E2%82%81%20%3D%20%7B_k1%3A.2f%7D)%3A%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22edge%20pixels%20cover%20%7Babs(_pct)%3A.1f%7D%25%20**more**%20object%20area%20than%20the%20centre%22)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20_dist_desc%20%3D%20(f%22pincushion%20distortion%20(k%E2%82%81%20%3D%20%7B_k1%3A.2f%7D)%3A%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22edge%20pixels%20cover%20%7Babs(_pct)%3A.1f%7D%25%20**less**%20object%20area%20than%20the%20centre%22)%0A%0A%0A%0A%20%20%20%20mo.vstack(%5B_expl%2C%20mo.hstack(%5Bpf_sl%2C%20k1_sl%2C%20tilt_sl%5D%2C%20justify%3D%22start%22)%2C%20_fig%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20camera_math(map_coordinates%2C%20np)%3A%0A%20%20%20%20%22%22%22Pure-numpy%20%2F%20scipy%20camera%20geometry%20primitives%20(no%20OpenCV).%22%22%22%0A%0A%20%20%20%20def%20rodrigues_to_matrix(rvec)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Rotation%20vector%20%E2%86%92%203%C3%973%20rotation%20matrix%20(Rodrigues%20formula).%22%22%22%0A%20%20%20%20%20%20%20%20r%20%3D%20np.asarray(rvec%2C%20dtype%3Dnp.float64).ravel()%0A%20%20%20%20%20%20%20%20theta%20%3D%20np.linalg.norm(r)%0A%20%20%20%20%20%20%20%20if%20theta%20%3C%201e-12%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20np.eye(3)%0A%20%20%20%20%20%20%20%20k%20%3D%20r%20%2F%20theta%0A%20%20%20%20%20%20%20%20K%20%3D%20np.array(%5B%5B0%2C%20-k%5B2%5D%2C%20k%5B1%5D%5D%2C%20%5Bk%5B2%5D%2C%200%2C%20-k%5B0%5D%5D%2C%20%5B-k%5B1%5D%2C%20k%5B0%5D%2C%200%5D%5D)%0A%20%20%20%20%20%20%20%20return%20np.eye(3)%20%2B%20np.sin(theta)%20*%20K%20%2B%20(1%20-%20np.cos(theta))%20*%20(K%20%40%20K)%0A%0A%20%20%20%20def%20project_points(obj_pts%2C%20rvec%2C%20tvec%2C%20K%2C%20dist)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Project%203D%20world%20points%20to%20distorted%20pixel%20coordinates.%0A%0A%20%20%20%20%20%20%20%20dist%20%3D%20%5Bk1%2C%20k2%2C%20p1%2C%20p2%2C%20k3%5D%0A%20%20%20%20%20%20%20%20Returns%20shape%20(N%2C%202).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20R%20%3D%20rodrigues_to_matrix(rvec)%0A%20%20%20%20%20%20%20%20t%20%3D%20np.asarray(tvec%2C%20dtype%3Dnp.float64).ravel()%0A%20%20%20%20%20%20%20%20pts_cam%20%3D%20(R%20%40%20obj_pts.T).T%20%2B%20t%0A%20%20%20%20%20%20%20%20xn%20%3D%20pts_cam%5B%3A%2C%200%5D%20%2F%20pts_cam%5B%3A%2C%202%5D%0A%20%20%20%20%20%20%20%20yn%20%3D%20pts_cam%5B%3A%2C%201%5D%20%2F%20pts_cam%5B%3A%2C%202%5D%0A%20%20%20%20%20%20%20%20r2%20%3D%20xn%20**%202%20%2B%20yn%20**%202%0A%20%20%20%20%20%20%20%20d%20%3D%20np.asarray(dist%2C%20dtype%3Dnp.float64).ravel()%0A%20%20%20%20%20%20%20%20k1%2C%20k2%20%3D%20d%5B0%5D%2C%20d%5B1%5D%0A%20%20%20%20%20%20%20%20p1%2C%20p2%20%3D%20(d%5B2%5D%2C%20d%5B3%5D)%20if%20len(d)%20%3E%203%20else%20(0.0%2C%200.0)%0A%20%20%20%20%20%20%20%20k3%20%3D%20d%5B4%5D%20if%20len(d)%20%3E%204%20else%200.0%0A%20%20%20%20%20%20%20%20radial%20%3D%201%20%2B%20k1%20*%20r2%20%2B%20k2%20*%20r2%20**%202%20%2B%20k3%20*%20r2%20**%203%0A%20%20%20%20%20%20%20%20xd%20%3D%20xn%20*%20radial%20%2B%202%20*%20p1%20*%20xn%20*%20yn%20%2B%20p2%20*%20(r2%20%2B%202%20*%20xn%20**%202)%0A%20%20%20%20%20%20%20%20yd%20%3D%20yn%20*%20radial%20%2B%20p1%20*%20(r2%20%2B%202%20*%20yn%20**%202)%20%2B%202%20*%20p2%20*%20xn%20*%20yn%0A%20%20%20%20%20%20%20%20u%20%3D%20K%5B0%2C%200%5D%20*%20xd%20%2B%20K%5B0%2C%202%5D%0A%20%20%20%20%20%20%20%20v%20%3D%20K%5B1%2C%201%5D%20*%20yd%20%2B%20K%5B1%2C%202%5D%0A%20%20%20%20%20%20%20%20return%20np.stack(%5Bu%2C%20v%5D%2C%20axis%3D1)%0A%0A%20%20%20%20def%20undistort_points_iterative(pts_dist%2C%20K%2C%20dist%2C%20max_iter%3D20)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Map%20distorted%20pixel%20coords%20to%20undistorted%20pixel%20coords.%0A%0A%20%20%20%20%20%20%20%20Uses%20iterative%20Newton%20refinement%20on%20the%20distortion%20model.%0A%20%20%20%20%20%20%20%20Returns%20shape%20(N%2C%202).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20d%20%3D%20np.asarray(dist%2C%20dtype%3Dnp.float64).ravel()%0A%20%20%20%20%20%20%20%20k1%2C%20k2%20%3D%20d%5B0%5D%2C%20d%5B1%5D%0A%20%20%20%20%20%20%20%20p1%2C%20p2%20%3D%20(d%5B2%5D%2C%20d%5B3%5D)%20if%20len(d)%20%3E%203%20else%20(0.0%2C%200.0)%0A%20%20%20%20%20%20%20%20k3%20%3D%20d%5B4%5D%20if%20len(d)%20%3E%204%20else%200.0%0A%20%20%20%20%20%20%20%20fx%2C%20fy%2C%20cx%2C%20cy%20%3D%20K%5B0%2C%200%5D%2C%20K%5B1%2C%201%5D%2C%20K%5B0%2C%202%5D%2C%20K%5B1%2C%202%5D%0A%20%20%20%20%20%20%20%20pts%20%3D%20np.asarray(pts_dist%2C%20dtype%3Dnp.float64).reshape(-1%2C%202)%0A%20%20%20%20%20%20%20%20xn%20%3D%20(pts%5B%3A%2C%200%5D%20-%20cx)%20%2F%20fx%0A%20%20%20%20%20%20%20%20yn%20%3D%20(pts%5B%3A%2C%201%5D%20-%20cy)%20%2F%20fy%0A%20%20%20%20%20%20%20%20%23%20Initial%20guess%3A%20normalized%20distorted%20coords%0A%20%20%20%20%20%20%20%20x0%2C%20y0%20%3D%20xn.copy()%2C%20yn.copy()%0A%20%20%20%20%20%20%20%20for%20_%20in%20range(max_iter)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20r2%20%3D%20x0%20**%202%20%2B%20y0%20**%202%0A%20%20%20%20%20%20%20%20%20%20%20%20radial%20%3D%201%20%2B%20k1%20*%20r2%20%2B%20k2%20*%20r2%20**%202%20%2B%20k3%20*%20r2%20**%203%0A%20%20%20%20%20%20%20%20%20%20%20%20dx%20%3D%202%20*%20p1%20*%20x0%20*%20y0%20%2B%20p2%20*%20(r2%20%2B%202%20*%20x0%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20dy%20%3D%20p1%20*%20(r2%20%2B%202%20*%20y0%20**%202)%20%2B%202%20*%20p2%20*%20x0%20*%20y0%0A%20%20%20%20%20%20%20%20%20%20%20%20x0%20%3D%20(xn%20-%20dx)%20%2F%20radial%0A%20%20%20%20%20%20%20%20%20%20%20%20y0%20%3D%20(yn%20-%20dy)%20%2F%20radial%0A%20%20%20%20%20%20%20%20return%20np.stack(%5Bfx%20*%20x0%20%2B%20cx%2C%20fy%20*%20y0%20%2B%20cy%5D%2C%20axis%3D1)%0A%0A%20%20%20%20def%20undistort_image(img%2C%20K%2C%20dist)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Undistort%20a%20grayscale%20image.%0A%0A%20%20%20%20%20%20%20%20For%20each%20output%20pixel%2C%20maps%20back%20to%20the%20distorted%20source%20via%0A%20%20%20%20%20%20%20%20undistort_points_iterative%2C%20then%20samples%20with%20bilinear%20interp.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20h%2C%20w%20%3D%20img.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20v_coords%2C%20u_coords%20%3D%20np.mgrid%5B0%3Ah%2C%200%3Aw%5D%0A%20%20%20%20%20%20%20%20pts_dist%20%3D%20np.stack(%5Bu_coords.ravel()%2C%20v_coords.ravel()%5D%2C%20axis%3D1).astype(%0A%20%20%20%20%20%20%20%20%20%20%20%20np.float64%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20pts_src%20%3D%20undistort_points_iterative(pts_dist%2C%20K%2C%20dist)%0A%20%20%20%20%20%20%20%20map_u%20%3D%20pts_src%5B%3A%2C%200%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20map_v%20%3D%20pts_src%5B%3A%2C%201%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20out%20%3D%20map_coordinates(%0A%20%20%20%20%20%20%20%20%20%20%20%20img.astype(np.float64)%2C%20%5Bmap_v%2C%20map_u%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20order%3D1%2C%20mode%3D%22constant%22%2C%20cval%3D128.0%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20np.clip(out%2C%200%2C%20255).astype(np.uint8)%0A%0A%20%20%20%20def%20find_homography_dlt(src%2C%20dst)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Estimate%20a%203%C3%973%20homography%20H%20such%20that%20dst%20%E2%89%88%20H%20%40%20src%20(homogeneous).%0A%0A%20%20%20%20%20%20%20%20Uses%20the%20direct%20linear%20transform%20(DLT)%20%E2%80%94%20no%20RANSAC.%0A%20%20%20%20%20%20%20%20src%2C%20dst%3A%20(N%2C%202)%20arrays.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20N%20%3D%20len(src)%0A%20%20%20%20%20%20%20%20A%20%3D%20np.zeros((2%20*%20N%2C%209))%0A%20%20%20%20%20%20%20%20for%20i%2C%20((x%2C%20y)%2C%20(u%2C%20v))%20in%20enumerate(zip(src%2C%20dst))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20A%5B2%20*%20i%5D%20%3D%20%5B-x%2C%20-y%2C%20-1%2C%200%2C%200%2C%200%2C%20u%20*%20x%2C%20u%20*%20y%2C%20u%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20A%5B2%20*%20i%20%2B%201%5D%20%3D%20%5B0%2C%200%2C%200%2C%20-x%2C%20-y%2C%20-1%2C%20v%20*%20x%2C%20v%20*%20y%2C%20v%5D%0A%20%20%20%20%20%20%20%20_%2C%20_%2C%20Vt%20%3D%20np.linalg.svd(A)%0A%20%20%20%20%20%20%20%20H%20%3D%20Vt%5B-1%5D.reshape(3%2C%203)%0A%20%20%20%20%20%20%20%20return%20H%20%2F%20H%5B2%2C%202%5D%0A%0A%20%20%20%20def%20warp_perspective(img%2C%20H%2C%20out_size%2C%20border_val%3D128)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Apply%20a%20perspective%20warp%20H%20to%20a%20grayscale%20image.%0A%0A%20%20%20%20%20%20%20%20Uses%20backward%20mapping%3A%20for%20each%20output%20pixel%2C%20find%20its%20source%20in%20img.%0A%20%20%20%20%20%20%20%20out_size%3A%20(width%2C%20height).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20w%2C%20h%20%3D%20out_size%0A%20%20%20%20%20%20%20%20v_coords%2C%20u_coords%20%3D%20np.mgrid%5B0%3Ah%2C%200%3Aw%5D%0A%20%20%20%20%20%20%20%20pts_dst%20%3D%20np.stack(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20u_coords.ravel()%2C%20v_coords.ravel()%2C%20np.ones(h%20*%20w)%0A%20%20%20%20%20%20%20%20%5D%2C%20axis%3D0)%20%20%23%203%20%C3%97%20N%0A%20%20%20%20%20%20%20%20H_inv%20%3D%20np.linalg.inv(H)%0A%20%20%20%20%20%20%20%20pts_src%20%3D%20H_inv%20%40%20pts_dst%0A%20%20%20%20%20%20%20%20pts_src%20%2F%3D%20pts_src%5B2%3A3%2C%20%3A%5D%0A%20%20%20%20%20%20%20%20map_u%20%3D%20pts_src%5B0%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20map_v%20%3D%20pts_src%5B1%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20valid%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20(map_u%20%3E%3D%200)%20%26%20(map_u%20%3C%20img.shape%5B1%5D%20-%201)%20%26%0A%20%20%20%20%20%20%20%20%20%20%20%20(map_v%20%3E%3D%200)%20%26%20(map_v%20%3C%20img.shape%5B0%5D%20-%201)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20out%20%3D%20np.full((h%2C%20w)%2C%20border_val%2C%20dtype%3Dnp.float64)%0A%20%20%20%20%20%20%20%20out%5Bvalid%5D%20%3D%20map_coordinates(%0A%20%20%20%20%20%20%20%20%20%20%20%20img.astype(np.float64)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bmap_v%5Bvalid%5D%2C%20map_u%5Bvalid%5D%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20order%3D1%2C%20mode%3D%22constant%22%2C%20cval%3Dfloat(border_val)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20np.clip(out%2C%200%2C%20255).astype(np.uint8)%0A%0A%20%20%20%20def%20render_simulated_image(%0A%20%20%20%20%20%20%20%20board_img%2C%20obj_pts_3d%2C%20rvec%2C%20tvec%2C%20K%2C%20dist_coeffs%2C%20img_size%3D(640%2C%20480)%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Render%20a%20checkerboard%20as%20seen%20through%20a%20pinhole%2Bdistortion%20camera.%0A%0A%20%20%20%20%20%20%20%20Steps%3A%0A%20%20%20%20%20%20%20%201.%20Compute%20homography%20from%20board%20plane%20to%20undistorted%20image%20(no%20distortion)%0A%20%20%20%20%20%20%20%202.%20Warp%20board%20%E2%86%92%20undistorted%20view%0A%20%20%20%20%20%20%20%203.%20Backward-map%20distortion%3A%20for%20each%20distorted%20pixel%2C%20find%20undistorted%20source%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20w%2C%20h%20%3D%20img_size%0A%20%20%20%20%20%20%20%20board_h%2C%20board_w%20%3D%20board_img.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20src_pts%20%3D%20np.array(%5B%5B0%2C%200%5D%2C%20%5Bboard_w%2C%200%5D%2C%20%5Bboard_w%2C%20board_h%5D%2C%20%5B0%2C%20board_h%5D%5D%2C%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%20dtype%3Dnp.float64)%0A%20%20%20%20%20%20%20%20rows%2C%20cols%20%3D%206%2C%209%0A%20%20%20%20%20%20%20%20sq%20%3D%2025.0%0A%20%20%20%20%20%20%20%20dst_3d%20%3D%20np.array(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B0%2C%200%2C%200%5D%2C%20%5B(cols%20-%201)%20*%20sq%20%2B%20sq%2C%200%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B(cols%20-%201)%20*%20sq%20%2B%20sq%2C%20(rows%20-%201)%20*%20sq%20%2B%20sq%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B0%2C%20(rows%20-%201)%20*%20sq%20%2B%20sq%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%5D%2C%20dtype%3Dnp.float64)%0A%20%20%20%20%20%20%20%20%23%20Project%20board%20corners%20with%20NO%20distortion%20%E2%86%92%20undistorted%20image%20plane%0A%20%20%20%20%20%20%20%20no_dist%20%3D%20np.zeros(5)%0A%20%20%20%20%20%20%20%20proj_corners%20%3D%20project_points(dst_3d%2C%20rvec%2C%20tvec%2C%20K%2C%20no_dist)%0A%20%20%20%20%20%20%20%20H%20%3D%20find_homography_dlt(src_pts%2C%20proj_corners)%0A%20%20%20%20%20%20%20%20undistorted%20%3D%20warp_perspective(board_img%2C%20H%2C%20(w%2C%20h))%0A%20%20%20%20%20%20%20%20%23%20Backward-map%20distortion%0A%20%20%20%20%20%20%20%20v_coords%2C%20u_coords%20%3D%20np.mgrid%5B0%3Ah%2C%200%3Aw%5D%0A%20%20%20%20%20%20%20%20pts_dist%20%3D%20np.stack(%5Bu_coords.ravel()%2C%20v_coords.ravel()%5D%2C%20axis%3D1).astype(%0A%20%20%20%20%20%20%20%20%20%20%20%20np.float64%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20pts_undist%20%3D%20undistort_points_iterative(pts_dist%2C%20K%2C%20dist_coeffs)%0A%20%20%20%20%20%20%20%20%23%20Round-trip%20check%3A%20re-apply%20distortion%20and%20verify%20consistency%0A%20%20%20%20%20%20%20%20xn%20%3D%20(pts_undist%5B%3A%2C%200%5D%20-%20K%5B0%2C%202%5D)%20%2F%20K%5B0%2C%200%5D%0A%20%20%20%20%20%20%20%20yn%20%3D%20(pts_undist%5B%3A%2C%201%5D%20-%20K%5B1%2C%202%5D)%20%2F%20K%5B1%2C%201%5D%0A%20%20%20%20%20%20%20%20r2%20%3D%20xn%20**%202%20%2B%20yn%20**%202%0A%20%20%20%20%20%20%20%20d%20%3D%20np.asarray(dist_coeffs).ravel()%0A%20%20%20%20%20%20%20%20k1%2C%20k2%20%3D%20d%5B0%5D%2C%20d%5B1%5D%0A%20%20%20%20%20%20%20%20k3%20%3D%20d%5B4%5D%20if%20len(d)%20%3E%204%20else%200.0%0A%20%20%20%20%20%20%20%20radial%20%3D%201%20%2B%20k1%20*%20r2%20%2B%20k2%20*%20r2%20**%202%20%2B%20k3%20*%20r2%20**%203%0A%20%20%20%20%20%20%20%20u_pred%20%3D%20K%5B0%2C%200%5D%20*%20xn%20*%20radial%20%2B%20K%5B0%2C%202%5D%0A%20%20%20%20%20%20%20%20v_pred%20%3D%20K%5B1%2C%201%5D%20*%20yn%20*%20radial%20%2B%20K%5B1%2C%202%5D%0A%20%20%20%20%20%20%20%20invalid%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.abs(u_pred%20-%20pts_dist%5B%3A%2C%200%5D)%20%3E%201.0)%20%7C%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.abs(v_pred%20-%20pts_dist%5B%3A%2C%201%5D)%20%3E%201.0)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20map_u%20%3D%20pts_undist%5B%3A%2C%200%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20map_v%20%3D%20pts_undist%5B%3A%2C%201%5D.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20out%20%3D%20np.full((h%2C%20w)%2C%20128.0)%0A%20%20%20%20%20%20%20%20valid_mask%20%3D%20~invalid.reshape(h%2C%20w)%0A%20%20%20%20%20%20%20%20out%5Bvalid_mask%5D%20%3D%20map_coordinates(%0A%20%20%20%20%20%20%20%20%20%20%20%20undistorted.astype(np.float64)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bmap_v%5Bvalid_mask%5D%2C%20map_u%5Bvalid_mask%5D%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20order%3D1%2C%20mode%3D%22constant%22%2C%20cval%3D128.0%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20np.clip(out%2C%200%2C%20255).astype(np.uint8)%0A%0A%20%20%20%20return%20project_points%2C%20render_simulated_image%0A%0A%0A%40app.cell%0Adef%20corner_detection_markdown(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Doing%20better%20than%20the%20grid%3A%20Sub-Pixel%20refinement%0A%0A%20%20%20%20If%20we%20want%20to%20do%20better%20than%20the%20half-pixel%20provided%20to%20us%20by%20the%20pixel%20grid%2C%20we%20need%20a%20new%2C%20sub-pixel%2C%20method.%20There%20are%20a%20few%20ways%20to%20do%20this%2C%20but%20the%20idea%20is%20to%20look%20at%20the%20non%20local%20behavior%3A%20gradients%2C%20or%20correlation%20with%20a%20kernel%20with%20appropriate%20shape.%20This%20non%20local%20search%20brings%20in%20information%20that's%20not%20just%20from%20a%20region%20that%20has%20the%20size%20of%20the%20pixel%20(or%20angular-size%20of%20the%20sensor)%20and%20can%20therefore%20locate%20us%20much%20more%20accurately.%20It's%20like%20fitting%20a%20gaussian%20to%20a%20histogram%2C%20if%20that%20analogy%20helps%20you.%20The%20histogram%20has%20bins%20(maybe%20just%20a%20few)%2C%20but%20by%20a%20fit%20to%20a%20continuous%20(non-local)%20gaussian%20we%20can%20locate%20the%20mean%20with%20much%20higher%20precision%20than%20the%20bin-bin%20distance%20which%20is%20analogous%20to%20our%20pixel%20size.%0A%0A%20%20%20%20We%20compare%20two%20approaches%3A%0A%0A%20%20%20%20%7C%20Method%20%7C%20Origin%20%7C%20Principle%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20**%60cornerSubPix%60**%20%7C%20OpenCV%20built-in%20(%60cv2.cornerSubPix%60)%20%7C%20Gradient%20perpendicularity%20condition%20%E2%80%94%20general-purpose%20%7C%0A%20%20%20%20%7C%20**Saddle%20fitting**%20%7C%20Custom%20implementation%20%7C%20Quadratic%20surface%20model%20%E2%80%94%20checkerboard-specific%20%7C%0A%0A%20%20%20%20We%20implement%20both%20from%20scratch%20in%20numpy%20so%20we%20can%20study%20exactly%20what%20they%20do%0A%20%20%20%20and%20compare%20them%20on%20data%20where%20we%20know%20the%20ground%20truth.%0A%0A%20%20%20%20---%0A%0A%20%20%20%20%23%23%23%20Method%201%20%E2%80%94%20%60cornerSubPix%60%20(OpenCV's%20approach)%0A%0A%20%20%20%20%60cv2.cornerSubPix%60%20is%20the%20standard%20sub-pixel%20corner%20refiner%20in%20OpenCV%2C%20used%20in%0A%20%20%20%20virtually%20every%20checkerboard%20calibration%20pipeline.%20It%20works%20for%20*any*%20corner%20type%2C%0A%20%20%20%20not%20just%20checkerboards.%0A%0A%20%20%20%20**Setup%3A**%20we%20start%20with%20an%20integer-pixel%20seed%20%24%5Cmathbf%7Bp%7D_0%24%20(our%20initial%20guess%20for%0A%20%20%20%20the%20corner%20location)%20and%20define%20a%20search%20window%20of%20size%20%24(2w%2B1)%20%5Ctimes%20(2w%2B1)%24%20pixels%0A%20%20%20%20centred%20on%20it.%20The%20points%20%24%5Cmathbf%7Bq%7D_i%24%20are%20simply%20**every%20pixel%20position%20inside%0A%20%20%20%20this%20window**%20%E2%80%94%20a%20regular%20grid%20of%20%24(2w%2B1)%5E2%24%20positions.%20At%20each%20%24%5Cmathbf%7Bq%7D_i%24%20we%0A%20%20%20%20compute%20the%20image%20gradient%20%24%5Cnabla%20I(%5Cmathbf%7Bq%7D_i)%24%20using%20a%20Sobel%20filter.%0A%0A%20%20%20%20**The%20gradient%20perpendicularity%20condition%3A**%20at%20the%20true%20corner%20%24%5Cmathbf%7Bp%7D%24%2C%20the%0A%20%20%20%20gradient%20at%20any%20nearby%20point%20%24%5Cmathbf%7Bq%7D%24%20should%20be%20perpendicular%20to%20the%20vector%0A%20%20%20%20from%20%24%5Cmathbf%7Bp%7D%24%20to%20%24%5Cmathbf%7Bq%7D%24.%20Why%3F%20Because%20edges%20run%20*through*%20a%20corner%2C%20and%0A%20%20%20%20the%20gradient%20is%20*normal*%20to%20the%20edge%20%E2%80%94%20so%20if%20%24%5Cmathbf%7Bq%7D%24%20sits%20on%20an%20edge%20passing%0A%20%20%20%20through%20%24%5Cmathbf%7Bp%7D%24%2C%20the%20gradient%20at%20%24%5Cmathbf%7Bq%7D%24%20must%20be%20perpendicular%20to%20the%0A%20%20%20%20edge%20direction%2C%20which%20is%20along%20%24(%5Cmathbf%7Bq%7D%20-%20%5Cmathbf%7Bp%7D)%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20(%5Cmathbf%7Bq%7D%20-%20%5Cmathbf%7Bp%7D)%20%5Ccdot%20%5Cnabla%20I(%5Cmathbf%7Bq%7D)%20%3D%200%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Writing%20%24%5Cmathbf%7BG%7D_i%20%3D%20%5Cnabla%20I(%5Cmathbf%7Bq%7D_i)%24%20and%20summing%20over%20all%20points%20in%20the%0A%20%20%20%20window%20gives%20a%20linear%20system%20for%20%24%5Cmathbf%7Bp%7D%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7Bp%7D%20%3D%20%5Cleft(%5Csum_i%20%5Cmathbf%7BG%7D_i%20%5Cmathbf%7BG%7D_i%5ET%5Cright)%5E%7B-1%7D%0A%20%20%20%20%5Csum_i%20%5Cmathbf%7BG%7D_i%20%5Cmathbf%7BG%7D_i%5ET%20%5C%2C%20%5Cmathbf%7Bq%7D_i%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**Iteration%3A**%20the%20solution%20%24%5Cmathbf%7Bp%7D%24%20from%20this%20system%20becomes%20the%20new%20window%0A%20%20%20%20centre.%20We%20re-centre%20the%20window%20on%20the%20updated%20%24%5Cmathbf%7Bp%7D%24%2C%20recompute%20the%20gradients%2C%0A%20%20%20%20solve%20again%2C%20and%20repeat%20until%20the%20position%20shift%20is%20smaller%20than%20a%20threshold%20%24%5Cvarepsilon%24%0A%20%20%20%20or%20we%20hit%20a%20maximum%20iteration%20count.%0A%0A%20%20%20%20**Where%20does%20the%20initial%20guess%20%24%5Cmathbf%7Bp%7D_0%24%20come%20from%3F**%20In%20practice%2C%20a%20coarse%0A%20%20%20%20detector%20like%20OpenCV's%20%60findChessboardCorners%60%20finds%20corners%20at%20integer-pixel%0A%20%20%20%20accuracy%20%E2%80%94%20those%20are%20the%20seeds.%20In%20this%20notebook%20we%20use%20projected%20ground-truth%0A%20%20%20%20positions%20(rounded%20to%20integers)%20since%20our%20images%20are%20synthetic%20and%20we%20know%20the%0A%20%20%20%203D%20geometry%20exactly.%0A%0A%20%20%20%20**What%20if%20there's%20no%20corner%20in%20the%20window%3F**%20The%20matrix%0A%20%20%20%20%24%5Csum_i%20%5Cmathbf%7BG%7D_i%20%5Cmathbf%7BG%7D_i%5ET%24%20encodes%20the%20gradient%20structure%20inside%20the%0A%20%20%20%20window.%20At%20a%20true%20corner%2C%20gradients%20point%20in%20at%20least%20two%20independent%20directions%2C%0A%20%20%20%20so%20this%20matrix%20is%20well-conditioned%20(both%20eigenvalues%20are%20large).%20On%20a%20flat%20region%0A%20%20%20%20the%20matrix%20is%20nearly%20zero%3B%20along%20a%20straight%20edge%20it%20has%20only%20one%20large%20eigenvalue.%0A%20%20%20%20In%20both%20cases%20the%20system%20is%20ill-conditioned%20%E2%80%94%20the%20determinant%20is%20small%20and%20the%0A%20%20%20%20solver%20either%20fails%20or%20produces%20a%20large%2C%20unreliable%20shift.%20A%20robust%20implementation%0A%20%20%20%20checks%20for%20this%20and%20keeps%20the%20original%20seed%20when%20the%20matrix%20is%20degenerate.%0A%0A%20%20%20%20**Key%20parameters%3A**%0A%20%20%20%20-%20%60win_size%60%20(%24w%24)%3A%20half-size%20of%20the%20search%20window%20(e.g.%20%24w%3D5%24%20%E2%86%92%20%2411%20%5Ctimes%2011%24%20window)%0A%20%20%20%20-%20%60max_iter%60%2C%20%60eps%60%3A%20stopping%20criteria%20for%20the%20iteration%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20gradient_structure_demo(%0A%20%20%20%20K_true%2C%0A%20%20%20%20dist_true%2C%0A%20%20%20%20gaussian_filter%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obj_points_3d%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20poses%2C%0A%20%20%20%20project_points%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20from%20scipy.signal%20import%20convolve2d%20as%20_conv2d%0A%0A%20%20%20%20_img%20%3D%20simulated_images%5B0%5D%0A%20%20%20%20_rvec%2C%20_tvec%20%3D%20poses%5B0%5D%0A%20%20%20%20_img_f%20%3D%20gaussian_filter(_img.astype(np.float64)%2C%20sigma%3D0.5)%0A%20%20%20%20_kx%20%3D%20np.array(%5B%5B-1%2C%200%2C%201%5D%2C%20%5B-2%2C%200%2C%202%5D%2C%20%5B-1%2C%200%2C%201%5D%5D%2C%20dtype%3Dnp.float64)%0A%20%20%20%20_Gx%20%3D%20_conv2d(_img_f%2C%20_kx%2C%20mode%3D%22same%22)%0A%20%20%20%20_Gy%20%3D%20_conv2d(_img_f%2C%20_kx.T%2C%20mode%3D%22same%22)%0A%0A%20%20%20%20%23%20Find%20a%20corner%2C%20an%20edge%20midpoint%2C%20and%20a%20flat%20region%0A%20%20%20%20_gt%20%3D%20project_points(obj_points_3d%2C%20_rvec%2C%20_tvec%2C%20K_true%2C%20dist_true)%0A%20%20%20%20_h%2C%20_w%20%3D%20_img.shape%5B%3A2%5D%0A%20%20%20%20_cx%2C%20_cy%20%3D%20_w%20%2F%202.0%2C%20_h%20%2F%202.0%0A%20%20%20%20%23%20Pick%20a%20corner%20near%20the%20centre%0A%20%20%20%20_dists%20%3D%20np.sqrt((_gt%5B%3A%2C%200%5D%20-%20_cx)%20**%202%20%2B%20(_gt%5B%3A%2C%201%5D%20-%20_cy)%20**%202)%0A%20%20%20%20_ci%20%3D%20int(np.argmin(_dists))%0A%20%20%20%20_corner_pt%20%3D%20_gt%5B_ci%5D%0A%20%20%20%20%23%20Edge%20midpoint%3A%20halfway%20between%20this%20corner%20and%20the%20next%20one%20in%20the%20same%20row%0A%20%20%20%20_edge_pt%20%3D%20(_gt%5B_ci%5D%20%2B%20_gt%5B_ci%20%2B%201%5D)%20%2F%202.0%0A%20%20%20%20%23%20Flat%3A%20offset%20half%20a%20square%20from%20the%20corner%20(into%20a%20white%2Fblack%20square)%0A%20%20%20%20_flat_pt%20%3D%20_corner_pt%20%2B%20np.array(%5B15.0%2C%2015.0%5D)%0A%0A%20%20%20%20_win%20%3D%206%0A%20%20%20%20_locations%20%3D%20%5B%0A%20%20%20%20%20%20%20%20(%22Corner%22%2C%20_corner_pt)%2C%0A%20%20%20%20%20%20%20%20(%22Edge%20(no%20corner)%22%2C%20_edge_pt)%2C%0A%20%20%20%20%20%20%20%20(%22Flat%20region%22%2C%20_flat_pt)%2C%0A%20%20%20%20%5D%0A%0A%20%20%20%20_fig_gs%2C%20_axes_gs%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(15%2C%205))%0A%0A%20%20%20%20for%20_ax%2C%20(_label%2C%20_pt)%20in%20zip(_axes_gs%2C%20_locations)%3A%0A%20%20%20%20%20%20%20%20_ix%2C%20_iy%20%3D%20int(round(_pt%5B0%5D))%2C%20int(round(_pt%5B1%5D))%0A%20%20%20%20%20%20%20%20_x0%20%3D%20max(0%2C%20_ix%20-%20_win)%0A%20%20%20%20%20%20%20%20_x1%20%3D%20min(_w%2C%20_ix%20%2B%20_win)%0A%20%20%20%20%20%20%20%20_y0%20%3D%20max(0%2C%20_iy%20-%20_win)%0A%20%20%20%20%20%20%20%20_y1%20%3D%20min(_h%2C%20_iy%20%2B%20_win)%0A%0A%20%20%20%20%20%20%20%20_ax.imshow(_img%5B_y0%3A_y1%2C%20_x0%3A_x1%5D%2C%20cmap%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20extent%3D%5B_x0%2C%20_x1%2C%20_y1%2C%20_y0%5D%2C%20interpolation%3D%22nearest%22)%0A%0A%20%20%20%20%20%20%20%20%23%20Draw%20gradient%20arrows%0A%20%20%20%20%20%20%20%20_step%20%3D%201%0A%20%20%20%20%20%20%20%20for%20_y%20in%20range(_y0%2C%20_y1%2C%20_step)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_x%20in%20range(_x0%2C%20_x1%2C%20_step)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_gx%20%3D%20_Gx%5B_y%2C%20_x%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_gy%20%3D%20_Gy%5B_y%2C%20_x%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_mag%20%3D%20np.sqrt(_gx%20**%202%20%2B%20_gy%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20_mag%20%3E%205%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_sc%20%3D%201.2%20%2F%20max(_mag%2C%201e-6)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_ax.arrow(_x%2C%20_y%2C%20_gx%20*%20_sc%2C%20_gy%20*%20_sc%2C%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%20head_width%3D0.25%2C%20head_length%3D0.15%2C%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%20fc%3D%22cyan%22%2C%20ec%3D%22cyan%22%2C%20alpha%3D0.8)%0A%0A%20%20%20%20%20%20%20%20%23%20Compute%20%CE%A3G%E1%B5%A2G%E1%B5%A2%E1%B5%80%20and%20its%20eigenvalues%0A%20%20%20%20%20%20%20%20_gx_win%20%3D%20_Gx%5B_y0%3A_y1%2C%20_x0%3A_x1%5D.ravel()%0A%20%20%20%20%20%20%20%20_gy_win%20%3D%20_Gy%5B_y0%3A_y1%2C%20_x0%3A_x1%5D.ravel()%0A%20%20%20%20%20%20%20%20_A%20%3D%20np.array(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bnp.sum(_gx_win%20**%202)%2C%20np.sum(_gx_win%20*%20_gy_win)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bnp.sum(_gx_win%20*%20_gy_win)%2C%20np.sum(_gy_win%20**%202)%5D%2C%0A%20%20%20%20%20%20%20%20%5D)%0A%20%20%20%20%20%20%20%20_eigvals%20%3D%20np.sort(np.linalg.eigvalsh(_A))%5B%3A%3A-1%5D%0A%20%20%20%20%20%20%20%20_det%20%3D%20_eigvals%5B0%5D%20*%20_eigvals%5B1%5D%0A%20%20%20%20%20%20%20%20_ratio%20%3D%20_eigvals%5B0%5D%20%2F%20max(_eigvals%5B1%5D%2C%201e-6)%0A%0A%20%20%20%20%20%20%20%20if%20_eigvals%5B1%5D%20%3E%201e3%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_verdict%20%3D%20%22CORNER%22%0A%20%20%20%20%20%20%20%20%20%20%20%20_color%20%3D%20%22lime%22%0A%20%20%20%20%20%20%20%20elif%20_eigvals%5B0%5D%20%3E%201e3%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_verdict%20%3D%20%22EDGE%22%0A%20%20%20%20%20%20%20%20%20%20%20%20_color%20%3D%20%22orange%22%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_verdict%20%3D%20%22FLAT%22%0A%20%20%20%20%20%20%20%20%20%20%20%20_color%20%3D%20%22red%22%0A%0A%20%20%20%20%20%20%20%20_ax.set_title(f%22%7B_label%7D%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%CE%BB%E2%82%81%3D%7B_eigvals%5B0%5D%3A.0f%7D%20%20%CE%BB%E2%82%82%3D%7B_eigvals%5B1%5D%3A.0f%7D%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%E2%86%92%20%7B_verdict%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D10%2C%20color%3D_color%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_pt%5B0%5D%2C%20_pt%5B1%5D%2C%20%22%2B%22%2C%20color%3D_color%2C%20markersize%3D14%2C%20markeredgewidth%3D2)%0A%20%20%20%20%20%20%20%20_ax.axis(%22off%22)%0A%0A%20%20%20%20plt.suptitle(%22Gradient%20structure%20at%20three%20locations%3A%20why%20the%20detector%20needs%20a%20corner%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D12)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20_expl_gs%20%3D%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Why%20does%20%60cornerSubPix%60%20need%20a%20corner%3F%0A%0A%20%20%20%20The%20matrix%20%24%5Cmathbf%7BA%7D%20%3D%20%5Csum_i%20%5Cmathbf%7BG%7D_i%20%5Cmathbf%7BG%7D_i%5ET%24%20encodes%20the%20gradient%0A%20%20%20%20structure%20inside%20the%20search%20window.%20Its%20two%20eigenvalues%20%24%5Clambda_1%20%5Cgeq%20%5Clambda_2%24%0A%20%20%20%20tell%20us%20what%20kind%20of%20feature%20is%20present%3A%0A%0A%20%20%20%20%7C%20%24%5Clambda_1%24%20%7C%20%24%5Clambda_2%24%20%7C%20Structure%20%7C%20Can%20we%20localise%3F%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20Large%20%7C%20Large%20%7C%20**Corner**%20%E2%80%94%20gradients%20in%20two%20directions%20%7C%20Yes%20%E2%80%94%20system%20is%20well-conditioned%20%7C%0A%20%20%20%20%7C%20Large%20%7C%20Small%20%7C%20**Edge**%20%E2%80%94%20gradients%20in%20one%20direction%20only%20%7C%20Only%20along%20the%20edge%20normal%2C%20not%20along%20it%20%7C%0A%20%20%20%20%7C%20Small%20%7C%20Small%20%7C%20**Flat**%20%E2%80%94%20no%20gradient%20signal%20%7C%20No%20%E2%80%94%20matrix%20is%20nearly%20singular%20%7C%0A%0A%20%20%20%20The%20three%20patches%20below%20show%20this%20on%20our%20checkerboard%3A%20the%20corner%20has%20strong%20gradients%0A%20%20%20%20in%20two%20directions%20(both%20eigenvalues%20large)%2C%20the%20edge%20has%20gradients%20only%20perpendicular%0A%20%20%20%20to%20itself%20(one%20eigenvalue%20dominates)%2C%20and%20the%20flat%20square%20interior%20has%20almost%20no%20gradient%0A%20%20%20%20at%20all.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20mo.vstack(%5B_expl_gs%2C%20_fig_gs%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20ground_truth_params(np)%3A%0A%20%20%20%20IMG_W%2C%20IMG_H%20%3D%20640%2C%20480%0A%20%20%20%20FX%2C%20FY%20%3D%20520.0%2C%20520.0%0A%20%20%20%20CX%2C%20CY%20%3D%20320.0%2C%20240.0%0A%20%20%20%20K1_TRUE%2C%20K2_TRUE%20%3D%20-0.15%2C%200.05%0A%0A%20%20%20%20K_true%20%3D%20np.array(%5B%5BFX%2C%200%2C%20CX%5D%2C%20%5B0%2C%20FY%2C%20CY%5D%2C%20%5B0%2C%200%2C%201%5D%5D%2C%20dtype%3Dnp.float64)%0A%20%20%20%20dist_true%20%3D%20np.array(%5BK1_TRUE%2C%20K2_TRUE%2C%200.0%2C%200.0%2C%200.0%5D%2C%20dtype%3Dnp.float64)%0A%0A%20%20%20%20pose_configs%20%3D%20%5B%0A%20%20%20%20%20%20%20%20((0%2C%200%2C%200)%2C%20(0%2C%200%2C%20220))%2C%0A%20%20%20%20%20%20%20%20((12%2C%200%2C%200)%2C%20(5%2C%20-5%2C%20230))%2C%0A%20%20%20%20%20%20%20%20((-8%2C%2012%2C%200)%2C%20(-10%2C%208%2C%20210))%2C%0A%20%20%20%20%20%20%20%20((0%2C%20-15%2C%208)%2C%20(8%2C%200%2C%20235))%2C%0A%20%20%20%20%20%20%20%20((15%2C%208%2C%20-4)%2C%20(-5%2C%2010%2C%20205))%2C%0A%20%20%20%20%20%20%20%20((-12%2C%20-8%2C%204)%2C%20(3%2C%20-8%2C%20220))%2C%0A%20%20%20%20%20%20%20%20((8%2C%2015%2C%208)%2C%20(-8%2C%205%2C%20225))%2C%0A%20%20%20%20%20%20%20%20((-4%2C%20-12%2C%20-8)%2C%20(10%2C%20-3%2C%20215))%2C%0A%20%20%20%20%5D%0A%0A%20%20%20%20def%20make_pose(rot_deg%2C%20trans_mm)%3A%0A%20%20%20%20%20%20%20%20rot_rad%20%3D%20np.array(rot_deg%2C%20dtype%3Dnp.float64)%20*%20np.pi%20%2F%20180.0%0A%20%20%20%20%20%20%20%20rvec%20%3D%20rot_rad.reshape(3%2C%201)%0A%20%20%20%20%20%20%20%20tx%20%3D%20trans_mm%5B0%5D%20-%20112.5%0A%20%20%20%20%20%20%20%20ty%20%3D%20trans_mm%5B1%5D%20-%2075.0%0A%20%20%20%20%20%20%20%20tz%20%3D%20trans_mm%5B2%5D%0A%20%20%20%20%20%20%20%20tvec%20%3D%20np.array(%5Btx%2C%20ty%2C%20tz%5D%2C%20dtype%3Dnp.float64).reshape(3%2C%201)%0A%20%20%20%20%20%20%20%20return%20rvec%2C%20tvec%0A%0A%20%20%20%20poses%20%3D%20%5Bmake_pose(r%2C%20t)%20for%20r%2C%20t%20in%20pose_configs%5D%0A%20%20%20%20return%20IMG_H%2C%20IMG_W%2C%20K_true%2C%20dist_true%2C%20poses%0A%0A%0A%40app.cell%0Adef%20saddle_theory_markdown(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Method%202%20%E2%80%94%20Saddle-Point%20Fitting%20(checkerboard-specific)%0A%0A%20%20%20%20At%20a%20checkerboard%20corner%20the%20intensity%20surface%20forms%20a%20**saddle%20point**%20%E2%80%94%20if%20you%0A%20%20%20%20walk%20in%20one%20direction%20you%20go%20from%20black%20to%20white%2C%20and%20in%20the%20perpendicular%20direction%0A%20%20%20%20from%20white%20to%20black.%20This%20is%20a%20very%20specific%20geometric%20structure%20that%20we%20can%20model%0A%20%20%20%20directly.%0A%0A%20%20%20%20%23%23%23%23%20Step%201%3A%20Fit%20a%20quadratic%20surface%0A%0A%20%20%20%20Extract%20a%20%24(2w%2B1)%20%5Ctimes%20(2w%2B1)%24%20image%20patch%20centred%20on%20the%20integer-pixel%20seed.%0A%20%20%20%20In%20local%20coordinates%20%24(x%2C%20y)%24%20relative%20to%20the%20patch%20centre%2C%20model%20the%20intensity%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20I(x%2C%20y)%20%3D%20ax%5E2%20%2B%20bxy%20%2B%20cy%5E2%20%2B%20dx%20%2B%20ey%20%2B%20f%0A%20%20%20%20%24%24%0A%0A%20%20%20%20This%20has%206%20unknowns.%20We%20have%20%24(2w%2B1)%5E2%24%20pixel%20intensities%20(e.g.%20121%20for%20%24w%3D5%24)%2C%0A%20%20%20%20so%20the%20system%20is%20heavily%20overdetermined.%20We%20solve%20it%20by%20ordinary%20least-squares%3A%0A%20%20%20%20stack%20all%20pixel%20coordinates%20into%20a%20design%20matrix%20%24%5Cmathbf%7BA%7D%24%20where%20each%20row%20is%0A%20%20%20%20%24%5Bx_i%5E2%2C%5C%3B%20x_i%20y_i%2C%5C%3B%20y_i%5E2%2C%5C%3B%20x_i%2C%5C%3B%20y_i%2C%5C%3B%201%5D%24%2C%20and%20solve%0A%20%20%20%20%24%5Cmathbf%7BA%7D%20%5Cboldsymbol%7B%5Cbeta%7D%20%3D%20%5Cmathbf%7Bz%7D%24%20for%20the%206%20coefficients%0A%20%20%20%20%24%5Cboldsymbol%7B%5Cbeta%7D%20%3D%20%5Ba%2C%20b%2C%20c%2C%20d%2C%20e%2C%20f%5D%5ET%24%20via%20the%20pseudo-inverse.%0A%0A%20%20%20%20%23%23%23%23%20Step%202%3A%20Find%20the%20critical%20point%0A%0A%20%20%20%20The%20saddle%20(or%20extremum)%20is%20where%20%24%5Cnabla%20I%20%3D%200%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cfrac%7B%5Cpartial%20I%7D%7B%5Cpartial%20x%7D%20%3D%202ax%20%2B%20by%20%2B%20d%20%3D%200%20%5Cqquad%0A%20%20%20%20%5Cfrac%7B%5Cpartial%20I%7D%7B%5Cpartial%20y%7D%20%3D%20bx%20%2B%202cy%20%2B%20e%20%3D%200%0A%20%20%20%20%24%24%0A%0A%20%20%20%20This%20is%20a%20%242%20%5Ctimes%202%24%20linear%20system%20(the%20Hessian%20matrix%20times%20the%20position%0A%20%20%20%20equals%20the%20negative%20of%20the%20linear%20terms)%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cunderbrace%7B%5Cbegin%7Bbmatrix%7D%202a%20%26%20b%20%5C%5C%20b%20%26%202c%20%5Cend%7Bbmatrix%7D%7D_%7B%5Cmathbf%7BH%7D%7D%0A%20%20%20%20%5Cbegin%7Bbmatrix%7D%20x_0%20%5C%5C%20y_0%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%3D%20-%5Cbegin%7Bbmatrix%7D%20d%20%5C%5C%20e%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20solution%20%24(x_0%2C%20y_0)%24%20is%20the%20sub-pixel%20offset%20from%20the%20patch%20centre.%0A%0A%20%20%20%20%23%23%23%23%20Step%203%3A%20Verify%20it's%20actually%20a%20saddle%0A%0A%20%20%20%20A%20critical%20point%20can%20be%20a%20minimum%2C%20maximum%2C%20or%20saddle.%20For%20a%20checkerboard%20corner%0A%20%20%20%20we%20specifically%20need%20a%20saddle%2C%20which%20is%20characterised%20by%20the%20Hessian%20having%20one%0A%20%20%20%20positive%20and%20one%20negative%20eigenvalue.%20Equivalently%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cdet(%5Cmathbf%7BH%7D)%20%3D%204ac%20-%20b%5E2%20%3C%200%20%5Cquad%20%5Ctext%7B(saddle%20condition)%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20If%20this%20fails%20(e.g.%20we're%20on%20a%20flat%20region%20or%20a%20single%20edge)%2C%20the%20fit%20is%20rejected%0A%20%20%20%20and%20we%20keep%20the%20original%20seed%20position.%0A%0A%20%20%20%20%23%23%23%23%20Why%20is%20this%20checkerboard-specific%3F%0A%0A%20%20%20%20The%20quadratic%20saddle%20model%20*assumes*%20the%20intensity%20surface%20looks%20like%20a%20hyperbolic%0A%20%20%20%20paraboloid%20%E2%80%94%20exactly%20what%20a%20checkerboard%20corner%20produces.%20On%20other%20corner%20types%0A%20%20%20%20(e.g.%20a%20T-junction%2C%20or%20a%20corner%20of%20a%20polygon)%20the%20surface%20might%20not%20be%20a%20clean%0A%20%20%20%20saddle%2C%20and%20this%20method%20would%20either%20reject%20the%20point%20or%20give%20a%20biased%20result.%0A%20%20%20%20%60cornerSubPix%60%20(Method%201)%20makes%20no%20such%20assumption%20%E2%80%94%20it%20just%20looks%20for%20where%0A%20%20%20%20gradients%20converge%2C%20which%20works%20for%20any%20corner%20geometry.%0A%0A%20%20%20%20---%0A%0A%20%20%20%20%23%23%23%20Why%20compare%20the%20two%20methods%3F%0A%0A%20%20%20%20-%20**%60cornerSubPix%60**%20is%20general%3A%20it%20works%20for%20any%20corner%20type%20but%20has%20no%20knowledge%0A%20%20%20%20%20%20that%20the%20pattern%20is%20a%20checkerboard.%0A%20%20%20%20-%20**Saddle%20fitting**%20is%20purpose-built%3A%20it%20directly%20models%20the%20physics%20of%20the%0A%20%20%20%20%20%20intensity%20surface%20at%20a%20checkerboard%20corner.%0A%0A%20%20%20%20On%20synthetic%20data%20with%20known%20ground%20truth%20we%20can%20measure%20whether%20the%20specialised%0A%20%20%20%20model%20offers%20any%20accuracy%20advantage%20%E2%80%94%20and%20whether%20that%20advantage%20grows%20under%0A%20%20%20%20challenging%20conditions%20such%20as%20oblique%20viewing%20angles%20and%20image%20noise.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20corner_subpix_impl(gaussian_filter%2C%20np)%3A%0A%20%20%20%20%22%22%22%0A%20%20%20%20Numpy%20implementation%20of%20cornerSubPix%20(gradient%20perpendicularity%20condition).%0A%0A%20%20%20%20Equivalent%20to%20cv2.cornerSubPix%20with%20a%20Sobel-gradient%20window.%0A%20%20%20%20%22%22%22%0A%0A%20%20%20%20def%20corner_subpix(img%2C%20corners%2C%20win_size%3D5%2C%20max_iter%3D30%2C%20eps%3D1e-3)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Refine%20corner%20positions%20using%20the%20gradient%20perpendicularity%20condition.%0A%0A%20%20%20%20%20%20%20%20For%20each%20seed%20corner%2C%20iteratively%20solves%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20p%20%3D%20(%CE%A3G%E1%B5%A2G%E1%B5%A2%E1%B5%80)%E2%81%BB%C2%B9%20%C2%B7%20%CE%A3G%E1%B5%A2G%E1%B5%A2%E1%B5%80q%E1%B5%A2%0A%0A%20%20%20%20%20%20%20%20where%20G%E1%B5%A2%20%3D%20%E2%88%87I(q%E1%B5%A2)%20for%20points%20q%E1%B5%A2%20in%20the%20(2*win_size%2B1)%C2%B2%20search%20window.%0A%0A%20%20%20%20%20%20%20%20Parameters%0A%20%20%20%20%20%20%20%20----------%0A%20%20%20%20%20%20%20%20img%20%3A%20ndarray%2C%20shape%20(H%2C%20W)%0A%20%20%20%20%20%20%20%20corners%20%3A%20ndarray%2C%20shape%20(N%2C%202)%20%20float%20pixel%20positions%0A%20%20%20%20%20%20%20%20win_size%20%3A%20int%20%20%20half-window%0A%20%20%20%20%20%20%20%20max_iter%2C%20eps%20%3A%20convergence%20criteria%0A%0A%20%20%20%20%20%20%20%20Returns%0A%20%20%20%20%20%20%20%20-------%0A%20%20%20%20%20%20%20%20refined%20%3A%20ndarray%2C%20shape%20(N%2C%202)%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20img_f%20%3D%20gaussian_filter(img.astype(np.float64)%2C%20sigma%3D0.5)%0A%20%20%20%20%20%20%20%20h%2C%20w%20%3D%20img_f.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20n%20%3D%20corners.shape%5B0%5D%0A%20%20%20%20%20%20%20%20refined%20%3D%20corners.astype(np.float64).copy()%0A%0A%20%20%20%20%20%20%20%20%23%20Sobel%20kernels%0A%20%20%20%20%20%20%20%20kx%20%3D%20np.array(%5B%5B-1%2C%200%2C%201%5D%2C%20%5B-2%2C%200%2C%202%5D%2C%20%5B-1%2C%200%2C%201%5D%5D%2C%20dtype%3Dnp.float64)%0A%20%20%20%20%20%20%20%20ky%20%3D%20kx.T%0A%0A%20%20%20%20%20%20%20%20from%20scipy.signal%20import%20convolve2d%0A%20%20%20%20%20%20%20%20Gx%20%3D%20convolve2d(img_f%2C%20kx%2C%20mode%3D%22same%22)%0A%20%20%20%20%20%20%20%20Gy%20%3D%20convolve2d(img_f%2C%20ky%2C%20mode%3D%22same%22)%0A%0A%20%20%20%20%20%20%20%20%23%20Build%20local%20coordinate%20offsets%20for%20the%20window%0A%20%20%20%20%20%20%20%20r%20%3D%20win_size%0A%20%20%20%20%20%20%20%20ys_off%2C%20xs_off%20%3D%20np.mgrid%5B-r%3Ar%20%2B%201%2C%20-r%3Ar%20%2B%201%5D%20%20%23%20(2r%2B1%2C%202r%2B1)%0A%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20px%2C%20py%20%3D%20refined%5Bi%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_%20in%20range(max_iter)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ix%2C%20iy%20%3D%20int(round(px))%2C%20int(round(py))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20ix%20-%20r%20%3C%200%20or%20ix%20%2B%20r%20%3E%3D%20w%20or%20iy%20-%20r%20%3C%200%20or%20iy%20%2B%20r%20%3E%3D%20h%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Sample%20gradient%20at%20each%20window%20point%20(use%20integer%20pixel%20centre%20of%20win)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20gx%20%3D%20Gx%5Biy%20-%20r%3Aiy%20%2B%20r%20%2B%201%2C%20ix%20-%20r%3Aix%20%2B%20r%20%2B%201%5D.ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20gy%20%3D%20Gy%5Biy%20-%20r%3Aiy%20%2B%20r%20%2B%201%2C%20ix%20-%20r%3Aix%20%2B%20r%20%2B%201%5D.ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20World%20coords%20of%20window%20points%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20qx%20%3D%20(ix%20%2B%20xs_off).ravel().astype(np.float64)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20qy%20%3D%20(iy%20%2B%20ys_off).ravel().astype(np.float64)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Build%20system%3A%20A%20%3D%20%CE%A3G%E1%B5%A2G%E1%B5%A2%E1%B5%80%2C%20b%20%3D%20%CE%A3G%E1%B5%A2G%E1%B5%A2%E1%B5%80q%E1%B5%A2%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20A00%20%3D%20np.sum(gx%20*%20gx)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20A01%20%3D%20np.sum(gx%20*%20gy)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20A11%20%3D%20np.sum(gy%20*%20gy)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20b0%20%3D%20np.sum(gx%20*%20gx%20*%20qx%20%2B%20gx%20*%20gy%20*%20qy)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20b1%20%3D%20np.sum(gy%20*%20gx%20*%20qx%20%2B%20gy%20*%20gy%20*%20qy)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20A%20%3D%20np.array(%5B%5BA00%2C%20A01%5D%2C%20%5BA01%2C%20A11%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20b_vec%20%3D%20np.array(%5Bb0%2C%20b1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20det%20%3D%20A00%20*%20A11%20-%20A01%20**%202%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20abs(det)%20%3C%201e-10%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20new_p%20%3D%20np.linalg.solve(A%2C%20b_vec)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20np.linalg.norm(new_p%20-%20np.array(%5Bpx%2C%20py%5D))%20%3C%20eps%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20px%2C%20py%20%3D%20new_p%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20px%2C%20py%20%3D%20new_p%0A%20%20%20%20%20%20%20%20%20%20%20%20refined%5Bi%5D%20%3D%20%5Bpx%2C%20py%5D%0A%20%20%20%20%20%20%20%20return%20refined%0A%0A%20%20%20%20return%20(corner_subpix%2C)%0A%0A%0A%40app.cell%0Adef%20saddle_function(np)%3A%0A%20%20%20%20def%20saddle_refine_corners(img%2C%20corners%2C%20win_size%3D5)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Refine%20corners%20by%20fitting%20a%20quadratic%20saddle%20surface.%0A%0A%20%20%20%20%20%20%20%20I(x%2Cy)%20%3D%20ax%C2%B2%2Bbxy%2Bcy%C2%B2%2Bdx%2Bey%2Bf%20%E2%80%94%20saddle%20point%20where%20%E2%88%87I%3D0.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20img_f%20%3D%20img.astype(np.float64)%0A%20%20%20%20%20%20%20%20h%2C%20w%20%3D%20img_f.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20refined%20%3D%20corners.astype(np.float64).copy()%0A%20%20%20%20%20%20%20%20win%20%3D%20win_size%0A%20%20%20%20%20%20%20%20ys%2C%20xs%20%3D%20np.mgrid%5B-win%3Awin%20%2B%201%2C%20-win%3Awin%20%2B%201%5D%0A%20%20%20%20%20%20%20%20xs_flat%20%3D%20xs.ravel().astype(np.float64)%0A%20%20%20%20%20%20%20%20ys_flat%20%3D%20ys.ravel().astype(np.float64)%0A%20%20%20%20%20%20%20%20A%20%3D%20np.column_stack(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20xs_flat%20**%202%2C%20xs_flat%20*%20ys_flat%2C%20ys_flat%20**%202%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20xs_flat%2C%20ys_flat%2C%20np.ones_like(xs_flat)%2C%0A%20%20%20%20%20%20%20%20%5D)%0A%20%20%20%20%20%20%20%20AtA_inv_At%20%3D%20np.linalg.pinv(A)%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(len(corners))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20cx_f%2C%20cy_f%20%3D%20corners%5Bi%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20ix%2C%20iy%20%3D%20int(round(cx_f))%2C%20int(round(cy_f))%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20ix%20-%20win%20%3C%200%20or%20ix%20%2B%20win%20%3E%3D%20w%20or%20iy%20-%20win%20%3C%200%20or%20iy%20%2B%20win%20%3E%3D%20h%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20%20%20%20%20patch%20%3D%20img_f%5Biy%20-%20win%3Aiy%20%2B%20win%20%2B%201%2C%20ix%20-%20win%3Aix%20%2B%20win%20%2B%201%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20coeffs%20%3D%20AtA_inv_At%20%40%20patch.ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20a%2C%20b%2C%20c%2C%20d%2C%20e%2C%20_f%20%3D%20coeffs%0A%20%20%20%20%20%20%20%20%20%20%20%20if%204%20*%20a%20*%20c%20-%20b%20*%20b%20%3E%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20%20%20%20%20H_mat%20%3D%20np.array(%5B%5B2%20*%20a%2C%20b%5D%2C%20%5Bb%2C%202%20*%20c%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20delta%20%3D%20np.linalg.solve(H_mat%2C%20np.array(%5B-d%2C%20-e%5D))%0A%20%20%20%20%20%20%20%20%20%20%20%20except%20np.linalg.LinAlgError%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20%20%20%20%20dx%2C%20dy%20%3D%20delta%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20abs(dx)%20%3C%20win%20and%20abs(dy)%20%3C%20win%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20refined%5Bi%5D%20%3D%20%5Bix%20%2B%20dx%2C%20iy%20%2B%20dy%5D%0A%20%20%20%20%20%20%20%20return%20refined%0A%0A%20%20%20%20return%20(saddle_refine_corners%2C)%0A%0A%0A%40app.cell%0Adef%20render_all_views(%0A%20%20%20%20IMG_H%2C%0A%20%20%20%20IMG_W%2C%0A%20%20%20%20K_true%2C%0A%20%20%20%20dist_true%2C%0A%20%20%20%20ideal_board%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20obj_points_3d%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20poses%2C%0A%20%20%20%20render_simulated_image%2C%0A)%3A%0A%20%20%20%20simulated_images%20%3D%20%5B%5D%0A%20%20%20%20for%20_rvec%2C%20_tvec%20in%20poses%3A%0A%20%20%20%20%20%20%20%20_img%20%3D%20render_simulated_image(%0A%20%20%20%20%20%20%20%20%20%20%20%20ideal_board%2C%20obj_points_3d%2C%20_rvec%2C%20_tvec%2C%20K_true%2C%20dist_true%2C%20(IMG_W%2C%20IMG_H)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20simulated_images.append(_img)%0A%0A%20%20%20%20_n_cols%20%3D%204%0A%20%20%20%20_n_rows%20%3D%20(len(simulated_images)%20%2B%20_n_cols%20-%201)%20%2F%2F%20_n_cols%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(_n_rows%2C%20_n_cols%2C%20figsize%3D(16%2C%204%20*%20_n_rows))%0A%20%20%20%20_flat%20%3D%20_axes.flatten()%0A%20%20%20%20for%20_i%2C%20_img%20in%20enumerate(simulated_images)%3A%0A%20%20%20%20%20%20%20%20_flat%5B_i%5D.imshow(_img%2C%20cmap%3D%22gray%22)%0A%20%20%20%20%20%20%20%20_flat%5B_i%5D.set_title(f%22View%20%7B_i%20%2B%201%7D%22)%0A%20%20%20%20%20%20%20%20_flat%5B_i%5D.axis(%22off%22)%0A%20%20%20%20for%20_i%20in%20range(len(simulated_images)%2C%20len(_flat))%3A%0A%20%20%20%20%20%20%20%20_flat%5B_i%5D.axis(%22off%22)%0A%20%20%20%20plt.suptitle(%22Simulated%20Camera%20Views%20with%20Radial%20Distortion%22%2C%20fontsize%3D14)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20mo.vstack(%5Bmo.md('''%23%23%20Testing%20corner%20detection%20on%20synthetic%20images%0A%0A%20%20%20%20We%20will%20now%20generate%20a%20set%20of%20checkeboard%20images%20at%20different%20pose%20angles%2C%20and%20with%20distortion%2C%20and%20run%20our%20two%20methods%20to%20detect%20corners.%20%0A%20%20%20%20''')%2C_fig%5D)%0A%20%20%20%20return%20(simulated_images%2C)%0A%0A%0A%40app.cell%0Adef%20run_detection(%0A%20%20%20%20K_true%2C%0A%20%20%20%20PATTERN_SIZE%2C%0A%20%20%20%20corner_subpix%2C%0A%20%20%20%20dist_true%2C%0A%20%20%20%20gaussian_filter%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obj_points_3d%2C%0A%20%20%20%20poses%2C%0A%20%20%20%20project_points%2C%0A%20%20%20%20saddle_refine_corners%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20%22%22%22%0A%20%20%20%20Seed%20corners%20by%20projecting%20known%203D%20points%2C%20then%20apply%20both%20refinement%20methods.%0A%0A%20%20%20%20This%20replaces%20cv2.findChessboardCorners%20for%20synthetic%20images.%0A%20%20%20%20%22%22%22%0A%20%20%20%20all_corners_pixel%20%3D%20%5B%5D%0A%20%20%20%20all_corners_subpix%20%3D%20%5B%5D%0A%20%20%20%20all_corners_saddle%20%3D%20%5B%5D%0A%20%20%20%20detection_results%20%3D%20%5B%5D%0A%0A%20%20%20%20for%20_i%2C%20(_img%2C%20(_rvec%2C%20_tvec))%20in%20enumerate(zip(simulated_images%2C%20poses))%3A%0A%20%20%20%20%20%20%20%20_img_blur%20%3D%20gaussian_filter(_img.astype(np.float64)%2C%20sigma%3D0.5)%0A%20%20%20%20%20%20%20%20_img_blur%20%3D%20np.clip(_img_blur%2C%200%2C%20255).astype(np.uint8)%0A%0A%20%20%20%20%20%20%20%20%23%20Project%203D%20corners%20to%20get%20integer-accurate%20seed%20positions%0A%20%20%20%20%20%20%20%20_proj%20%3D%20project_points(obj_points_3d%2C%20_rvec%2C%20_tvec%2C%20K_true%2C%20dist_true)%0A%20%20%20%20%20%20%20%20%23%20Filter%20corners%20that%20fall%20inside%20the%20image%0A%20%20%20%20%20%20%20%20h%2C%20w%20%3D%20_img.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20_valid%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20(_proj%5B%3A%2C%200%5D%20%3E%3D%205)%20%26%20(_proj%5B%3A%2C%200%5D%20%3C%20w%20-%205)%20%26%0A%20%20%20%20%20%20%20%20%20%20%20%20(_proj%5B%3A%2C%201%5D%20%3E%3D%205)%20%26%20(_proj%5B%3A%2C%201%5D%20%3C%20h%20-%205)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20if%20_valid.sum()%20%3C%20PATTERN_SIZE%5B0%5D%20*%20PATTERN_SIZE%5B1%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20all_corners_pixel.append(None)%0A%20%20%20%20%20%20%20%20%20%20%20%20all_corners_subpix.append(None)%0A%20%20%20%20%20%20%20%20%20%20%20%20all_corners_saddle.append(None)%0A%20%20%20%20%20%20%20%20%20%20%20%20detection_results.append(False)%0A%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%0A%20%20%20%20%20%20%20%20%23%20Use%20projected%20positions%20as%20integer-pixel%20seeds%0A%20%20%20%20%20%20%20%20_corners_pixel%20%3D%20np.round(_proj%5B_valid%5D).astype(np.float64)%0A%0A%20%20%20%20%20%20%20%20%23%20Sub-pixel%20refinement%0A%20%20%20%20%20%20%20%20_corners_subpix%20%3D%20corner_subpix(_img_blur%2C%20_corners_pixel%2C%20win_size%3D5)%0A%20%20%20%20%20%20%20%20_corners_saddle%20%3D%20saddle_refine_corners(_img_blur%2C%20_corners_pixel%2C%20win_size%3D5)%0A%0A%20%20%20%20%20%20%20%20all_corners_pixel.append(_corners_pixel)%0A%20%20%20%20%20%20%20%20all_corners_subpix.append(_corners_subpix)%0A%20%20%20%20%20%20%20%20all_corners_saddle.append(_corners_saddle)%0A%20%20%20%20%20%20%20%20detection_results.append(True)%0A%0A%20%20%20%20_n_detected%20%3D%20sum(detection_results)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20all_corners_pixel%2C%0A%20%20%20%20%20%20%20%20all_corners_saddle%2C%0A%20%20%20%20%20%20%20%20all_corners_subpix%2C%0A%20%20%20%20%20%20%20%20detection_results%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20challenging_view(%0A%20%20%20%20K_true%2C%0A%20%20%20%20corner_subpix%2C%0A%20%20%20%20dist_true%2C%0A%20%20%20%20gaussian_filter%2C%0A%20%20%20%20ideal_board%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obj_points_3d%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20project_points%2C%0A%20%20%20%20render_simulated_image%2C%0A%20%20%20%20saddle_refine_corners%2C%0A)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Corner%20Detection%20on%20a%20Challenging%20Oblique%20View%0A%0A%20%20%20%20The%20previous%20views%20use%20modest%20tilts%20(%E2%89%A415%C2%B0).%20Here%20we%20push%20to%20**~45%C2%B0%20combined%20tilt**%20%E2%80%94%0A%20%20%20%20severe%20perspective%20foreshortening%20compresses%20the%20far%20rows%20of%20corners%20and%20increases%0A%20%20%20%20sub-pixel%20shift%20magnitudes%2C%20making%20both%20the%20problem%20and%20the%20refinement%20more%20visible.%0A%20%20%20%20%22%22%22)%0A%0A%20%20%20%20%23%20Challenging%20pose%3A%20heavy%20tilt%0A%20%20%20%20_rot_deg%20%3D%20np.array(%5B42.0%2C%2022.0%2C%205.0%5D)%20*%20np.pi%20%2F%20180.0%0A%20%20%20%20_rvec_c%20%3D%20_rot_deg.reshape(3%2C%201)%0A%20%20%20%20_tvec_c%20%3D%20np.array(%5B-30.0%2C%20-20.0%2C%20290.0%5D).reshape(3%2C%201)%0A%0A%20%20%20%20_img_c%20%3D%20render_simulated_image(%0A%20%20%20%20%20%20%20%20ideal_board%2C%20obj_points_3d%2C%20_rvec_c%2C%20_tvec_c%2C%20K_true%2C%20dist_true%2C%20(640%2C%20480)%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Add%20mild%20noise%20to%20make%20it%20more%20realistic%0A%20%20%20%20_rng%20%3D%20np.random.RandomState(7)%0A%20%20%20%20_img_noisy%20%3D%20np.clip(%0A%20%20%20%20%20%20%20%20_img_c.astype(np.float64)%20%2B%20_rng.randn(*_img_c.shape)%20*%203.0%2C%200%2C%20255%0A%20%20%20%20).astype(np.uint8)%0A%0A%20%20%20%20_img_blur%20%3D%20gaussian_filter(_img_noisy.astype(np.float64)%2C%20sigma%3D0.5)%0A%20%20%20%20_img_blur_u8%20%3D%20np.clip(_img_blur%2C%200%2C%20255).astype(np.uint8)%0A%0A%20%20%20%20%23%20Seed%20corners%20from%20projected%20ground%20truth%0A%20%20%20%20_proj%20%3D%20project_points(obj_points_3d%2C%20_rvec_c%2C%20_tvec_c%2C%20K_true%2C%20dist_true)%0A%20%20%20%20_h%2C%20_w%20%3D%20_img_c.shape%5B%3A2%5D%0A%20%20%20%20_valid%20%3D%20(%0A%20%20%20%20%20%20%20%20(_proj%5B%3A%2C%200%5D%20%3E%3D%208)%20%26%20(_proj%5B%3A%2C%200%5D%20%3C%20_w%20-%208)%20%26%0A%20%20%20%20%20%20%20%20(_proj%5B%3A%2C%201%5D%20%3E%3D%208)%20%26%20(_proj%5B%3A%2C%201%5D%20%3C%20_h%20-%208)%0A%20%20%20%20)%0A%20%20%20%20_seeds%20%3D%20np.round(_proj%5B_valid%5D).astype(np.float64)%0A%20%20%20%20_proj_valid%20%3D%20_proj%5B_valid%5D%0A%0A%20%20%20%20_csp%20%3D%20corner_subpix(_img_blur_u8%2C%20_seeds%2C%20win_size%3D5)%0A%20%20%20%20_csd%20%3D%20saddle_refine_corners(_img_blur_u8%2C%20_seeds%2C%20win_size%3D5)%0A%0A%20%20%20%20%23%20Shift%20magnitudes%0A%20%20%20%20_shifts_sp%20%3D%20np.linalg.norm(_csp%20-%20_seeds%2C%20axis%3D1)%0A%20%20%20%20_shifts_sd%20%3D%20np.linalg.norm(_csd%20-%20_seeds%2C%20axis%3D1)%0A%0A%20%20%20%20%23%20---%20Figure%20---%0A%20%20%20%20_fig_c%2C%20_axes_c%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(17%2C%205))%0A%0A%20%20%20%20%23%20Panel%201%3A%20full%20image%20with%20all%20detected%20corners%0A%20%20%20%20_axes_c%5B0%5D.imshow(_img_noisy%2C%20cmap%3D%22gray%22%2C%20vmin%3D0%2C%20vmax%3D255)%0A%20%20%20%20_axes_c%5B0%5D.scatter(_seeds%5B%3A%2C%200%5D%2C%20_seeds%5B%3A%2C%201%5D%2C%20s%3D12%2C%20c%3D%22red%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Pixel%20seed%22%2C%20zorder%3D4)%0A%20%20%20%20_axes_c%5B0%5D.scatter(_csp%5B%3A%2C%200%5D%2C%20_csp%5B%3A%2C%201%5D%2C%20s%3D8%2C%20c%3D%22lime%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22cornerSubPix%22%2C%20zorder%3D5)%0A%20%20%20%20_axes_c%5B0%5D.legend(fontsize%3D8%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_axes_c%5B0%5D.set_title(%22Oblique%20view%20(~45%C2%B0%20tilt)%5CnAll%20detected%20corners%22)%0A%20%20%20%20_axes_c%5B0%5D.axis(%22off%22)%0A%0A%20%20%20%20%23%20Panel%202%3A%20quiver%20of%20sub-pixel%20shifts%20coloured%20by%20magnitude%0A%20%20%20%20_ax2%20%3D%20_axes_c%5B1%5D%0A%20%20%20%20_ax2.imshow(_img_noisy%2C%20cmap%3D%22gray%22%2C%20vmin%3D0%2C%20vmax%3D255)%0A%20%20%20%20_sc%20%3D%20_ax2.scatter(_seeds%5B%3A%2C%200%5D%2C%20_seeds%5B%3A%2C%201%5D%2C%20c%3D_shifts_sp%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22plasma%22%2C%20s%3D30%2C%20vmin%3D0%2C%20vmax%3D_shifts_sp.max()%2C%20zorder%3D5)%0A%20%20%20%20_scale%20%3D%208.0%0A%20%20%20%20for%20_i%20in%20range(len(_seeds))%3A%0A%20%20%20%20%20%20%20%20_dx%20%3D%20_csp%5B_i%2C%200%5D%20-%20_seeds%5B_i%2C%200%5D%0A%20%20%20%20%20%20%20%20_dy%20%3D%20_csp%5B_i%2C%201%5D%20-%20_seeds%5B_i%2C%201%5D%0A%20%20%20%20%20%20%20%20_ax2.arrow(_seeds%5B_i%2C%200%5D%2C%20_seeds%5B_i%2C%201%5D%2C%20_dx%20*%20_scale%2C%20_dy%20*%20_scale%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20head_width%3D1.5%2C%20head_length%3D1.0%2C%20fc%3D%22yellow%22%2C%20ec%3D%22yellow%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.7%2C%20length_includes_head%3DTrue)%0A%20%20%20%20plt.colorbar(_sc%2C%20ax%3D_ax2%2C%20label%3D%22Shift%20magnitude%20(px)%22%2C%20shrink%3D0.8)%0A%20%20%20%20_ax2.set_title(f%22Sub-pixel%20shift%20vectors%20(%C3%97%7B_scale%3A.0f%7D)%5Cncolour%20%3D%20magnitude%22)%0A%20%20%20%20_ax2.axis(%22off%22)%0A%0A%20%20%20%20%23%20Panel%203%3A%20zoomed%20corner%20%E2%80%94%20pick%20the%20one%20with%20largest%20shift%0A%20%20%20%20_worst%20%3D%20int(np.argmax(_shifts_sp))%0A%20%20%20%20_px_pt%20%3D%20_seeds%5B_worst%5D%0A%20%20%20%20_sp_pt%20%3D%20_csp%5B_worst%5D%0A%20%20%20%20_sd_pt%20%3D%20_csd%5B_worst%5D%0A%20%20%20%20_zoom%20%3D%2018%0A%20%20%20%20_x_lo%20%3D%20max(0%2C%20int(_px_pt%5B0%5D)%20-%20_zoom)%0A%20%20%20%20_x_hi%20%3D%20min(_img_c.shape%5B1%5D%2C%20int(_px_pt%5B0%5D)%20%2B%20_zoom)%0A%20%20%20%20_y_lo%20%3D%20max(0%2C%20int(_px_pt%5B1%5D)%20-%20_zoom)%0A%20%20%20%20_y_hi%20%3D%20min(_img_c.shape%5B0%5D%2C%20int(_px_pt%5B1%5D)%20%2B%20_zoom)%0A%20%20%20%20_axes_c%5B2%5D.imshow(%0A%20%20%20%20%20%20%20%20_img_noisy%5B_y_lo%3A_y_hi%2C%20_x_lo%3A_x_hi%5D%2C%20cmap%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20extent%3D%5B_x_lo%2C%20_x_hi%2C%20_y_hi%2C%20_y_lo%5D%2C%20interpolation%3D%22nearest%22%2C%0A%20%20%20%20)%0A%20%20%20%20_axes_c%5B2%5D.plot(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D%2C%20%22rs%22%2C%20markersize%3D13%2C%20label%3D%22Pixel%20seed%22%2C%20zorder%3D6)%0A%20%20%20%20_axes_c%5B2%5D.plot(_sp_pt%5B0%5D%2C%20_sp_pt%5B1%5D%2C%20%22g%5E%22%2C%20markersize%3D13%2C%20label%3D%22cornerSubPix%22%2C%20zorder%3D6)%0A%20%20%20%20_axes_c%5B2%5D.plot(_sd_pt%5B0%5D%2C%20_sd_pt%5B1%5D%2C%20%22bo%22%2C%20markersize%3D11%2C%20label%3D%22Saddle%20fit%22%2C%20zorder%3D6)%0A%20%20%20%20for%20(_pt_a%2C%20_pt_b%2C%20_col)%20in%20%5B%0A%20%20%20%20%20%20%20%20(_px_pt%2C%20_sp_pt%2C%20%22lime%22)%2C%20(_px_pt%2C%20_sd_pt%2C%20%22dodgerblue%22)%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_axes_c%5B2%5D.annotate(%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%2C%20xy%3D(_pt_b%5B0%5D%2C%20_pt_b%5B1%5D)%2C%20xytext%3D(_pt_a%5B0%5D%2C%20_pt_a%5B1%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D_col%2C%20lw%3D2)%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20_axes_c%5B2%5D.text(%0A%20%20%20%20%20%20%20%20_x_lo%20%2B%201%2C%20_y_lo%20%2B%202%2C%0A%20%20%20%20%20%20%20%20f%22SubPix%20shift%3A%20%7B_shifts_sp%5B_worst%5D%3A.3f%7D%20px%5CnSaddle%20shift%3A%20%7B_shifts_sd%5B_worst%5D%3A.3f%7D%20px%22%2C%0A%20%20%20%20%20%20%20%20color%3D%22white%22%2C%20fontsize%3D8%2C%20fontweight%3D%22bold%22%2C%0A%20%20%20%20%20%20%20%20bbox%3Ddict(boxstyle%3D%22round%22%2C%20facecolor%3D%22black%22%2C%20alpha%3D0.7)%2C%0A%20%20%20%20)%0A%20%20%20%20_axes_c%5B2%5D.legend(fontsize%3D8%2C%20loc%3D%22lower%20right%22)%0A%20%20%20%20_axes_c%5B2%5D.set_title(%22Zoomed%3A%20corner%20with%20largest%20shift%5Cn(far%20end%20of%20oblique%20board)%22)%0A%0A%20%20%20%20plt.suptitle(%0A%20%20%20%20%20%20%20%20f%22Oblique%20view%20%E2%80%94%20%7B_valid.sum()%7D%20corners%20detected%2C%20%22%0A%20%20%20%20%20%20%20%20f%22mean%20shift%20%7B_shifts_sp.mean()%3A.3f%7D%20px%20(SubPix)%20%2F%20%7B_shifts_sd.mean()%3A.3f%7D%20px%20(Saddle)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D12%2C%0A%20%20%20%20)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20mo.vstack(%5Bmo.md('''%23%23%20Detection%20examples%0A%0A%20%20%20%20We%20select%20a%20particularly%20nasty%20pose%20and%20run%20the%20detection%20algorithms.%20We%20can%20see%20they%20all%20align%20quite%20well.%20There's%20no%20drama%20here.%20%0A%20%20%20%20''')%2C_fig_c%5D)%20%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20subpix_visualization(%0A%20%20%20%20all_corners_pixel%2C%0A%20%20%20%20all_corners_subpix%2C%0A%20%20%20%20detection_results%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20_view_idx%20%3D%20next(_i%20for%20_i%2C%20_d%20in%20enumerate(detection_results)%20if%20_d)%0A%20%20%20%20_img_show%20%3D%20simulated_images%5B_view_idx%5D%0A%20%20%20%20_cpx%20%3D%20all_corners_pixel%5B_view_idx%5D%0A%20%20%20%20_csp%20%3D%20all_corners_subpix%5B_view_idx%5D%0A%0A%20%20%20%20_fig_corners%2C%20_axes_c%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205))%0A%0A%20%20%20%20_img_rgb%20%3D%20np.stack(%5B_img_show%5D%20*%203%2C%20axis%3D-1)%0A%20%20%20%20_axes_c%5B0%5D.imshow(_img_rgb)%0A%20%20%20%20_axes_c%5B0%5D.plot(_cpx%5B%3A%2C%200%5D%2C%20_cpx%5B%3A%2C%201%5D%2C%20%22r%2B%22%2C%20markersize%3D8%2C%20label%3D%22Pixel%20corners%22)%0A%20%20%20%20_axes_c%5B0%5D.plot(_csp%5B%3A%2C%200%5D%2C%20_csp%5B%3A%2C%201%5D%2C%20%22gx%22%2C%20markersize%3D6%2C%20label%3D%22Sub-pixel%20corners%22)%0A%20%20%20%20_axes_c%5B0%5D.legend()%0A%20%20%20%20_axes_c%5B0%5D.set_title(f%22View%20%7B_view_idx%20%2B%201%7D%3A%20Detected%20Corners%22)%0A%20%20%20%20_axes_c%5B0%5D.axis(%22off%22)%0A%0A%20%20%20%20_corner_idx%20%3D%20len(_cpx)%20%2F%2F%202%0A%20%20%20%20_px_pt%20%3D%20_cpx%5B_corner_idx%5D%0A%20%20%20%20_sp_pt%20%3D%20_csp%5B_corner_idx%5D%0A%20%20%20%20_dx%20%3D%20_sp_pt%5B0%5D%20-%20_px_pt%5B0%5D%0A%20%20%20%20_dy%20%3D%20_sp_pt%5B1%5D%20-%20_px_pt%5B1%5D%0A%20%20%20%20_zoom_r%20%3D%2015%0A%20%20%20%20_x_lo%20%3D%20max(0%2C%20int(_px_pt%5B0%5D)%20-%20_zoom_r)%0A%20%20%20%20_x_hi%20%3D%20min(_img_show.shape%5B1%5D%2C%20int(_px_pt%5B0%5D)%20%2B%20_zoom_r)%0A%20%20%20%20_y_lo%20%3D%20max(0%2C%20int(_px_pt%5B1%5D)%20-%20_zoom_r)%0A%20%20%20%20_y_hi%20%3D%20min(_img_show.shape%5B0%5D%2C%20int(_px_pt%5B1%5D)%20%2B%20_zoom_r)%0A%20%20%20%20_axes_c%5B1%5D.imshow(%0A%20%20%20%20%20%20%20%20_img_show%5B_y_lo%3A_y_hi%2C%20_x_lo%3A_x_hi%5D%2C%20cmap%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20extent%3D%5B_x_lo%2C%20_x_hi%2C%20_y_hi%2C%20_y_lo%5D%2C%20interpolation%3D%22nearest%22%2C%0A%20%20%20%20)%0A%20%20%20%20_axes_c%5B1%5D.plot(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D%2C%20%22rs%22%2C%20markersize%3D12%2C%20label%3D%22Pixel%22)%0A%20%20%20%20_axes_c%5B1%5D.plot(_sp_pt%5B0%5D%2C%20_sp_pt%5B1%5D%2C%20%22g%5E%22%2C%20markersize%3D12%2C%20label%3D%22Sub-pixel%22)%0A%20%20%20%20_axes_c%5B1%5D.annotate(%0A%20%20%20%20%20%20%20%20%22%22%2C%20xy%3D(_sp_pt%5B0%5D%2C%20_sp_pt%5B1%5D)%2C%20xytext%3D(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D)%2C%0A%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D%22yellow%22%2C%20lw%3D2)%2C%0A%20%20%20%20)%0A%20%20%20%20_axes_c%5B1%5D.text(%0A%20%20%20%20%20%20%20%20_px_pt%5B0%5D%20%2B%201%2C%20_px_pt%5B1%5D%20-%202%2C%0A%20%20%20%20%20%20%20%20f%22dx%3D%7B_dx%3A.3f%7D%5Cndy%3D%7B_dy%3A.3f%7D%22%2C%20color%3D%22yellow%22%2C%20fontsize%3D9%2C%0A%20%20%20%20%20%20%20%20fontweight%3D%22bold%22%2C%20bbox%3Ddict(boxstyle%3D%22round%22%2C%20facecolor%3D%22black%22%2C%20alpha%3D0.7)%2C%0A%20%20%20%20)%0A%20%20%20%20_axes_c%5B1%5D.legend(fontsize%3D8)%0A%20%20%20%20_axes_c%5B1%5D.set_title(%22Zoomed%3A%20Sub-pixel%20Shift%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20mo.md(%22%23%23%20Pixel%20vs%20Sub-Pixel%20Corner%20Localization%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20three_way_comparison(%0A%20%20%20%20all_corners_pixel%2C%0A%20%20%20%20all_corners_saddle%2C%0A%20%20%20%20all_corners_subpix%2C%0A%20%20%20%20detection_results%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20_view_idx%20%3D%20next(_i%20for%20_i%2C%20_d%20in%20enumerate(detection_results)%20if%20_d)%0A%20%20%20%20_img_show%20%3D%20simulated_images%5B_view_idx%5D%0A%20%20%20%20_cpx%20%3D%20all_corners_pixel%5B_view_idx%5D%0A%20%20%20%20_csp%20%3D%20all_corners_subpix%5B_view_idx%5D%0A%20%20%20%20_csd%20%3D%20all_corners_saddle%5B_view_idx%5D%0A%0A%20%20%20%20_corner_idx%20%3D%20len(_cpx)%20%2F%2F%202%0A%20%20%20%20_px_pt%20%3D%20_cpx%5B_corner_idx%5D%0A%20%20%20%20_sp_pt%20%3D%20_csp%5B_corner_idx%5D%0A%20%20%20%20_sd_pt%20%3D%20_csd%5B_corner_idx%5D%0A%0A%20%20%20%20_zoom_r%20%3D%2015%0A%20%20%20%20_cx_i%20%3D%20int(_px_pt%5B0%5D)%0A%20%20%20%20_cy_i%20%3D%20int(_px_pt%5B1%5D)%0A%20%20%20%20_x_lo%20%3D%20max(0%2C%20_cx_i%20-%20_zoom_r)%0A%20%20%20%20_x_hi%20%3D%20min(_img_show.shape%5B1%5D%2C%20_cx_i%20%2B%20_zoom_r)%0A%20%20%20%20_y_lo%20%3D%20max(0%2C%20_cy_i%20-%20_zoom_r)%0A%20%20%20%20_y_hi%20%3D%20min(_img_show.shape%5B0%5D%2C%20_cy_i%20%2B%20_zoom_r)%0A%0A%20%20%20%20_fig_3w%2C%20_ax_3w%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(7%2C%207))%0A%20%20%20%20_ax_3w.imshow(%0A%20%20%20%20%20%20%20%20_img_show%5B_y_lo%3A_y_hi%2C%20_x_lo%3A_x_hi%5D%2C%20cmap%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20extent%3D%5B_x_lo%2C%20_x_hi%2C%20_y_hi%2C%20_y_lo%5D%2C%20interpolation%3D%22nearest%22%2C%0A%20%20%20%20)%0A%20%20%20%20_ax_3w.plot(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D%2C%20%22rs%22%2C%20markersize%3D14%2C%20label%3D%22Pixel%22%2C%20zorder%3D5)%0A%20%20%20%20_ax_3w.plot(_sp_pt%5B0%5D%2C%20_sp_pt%5B1%5D%2C%20%22g%5E%22%2C%20markersize%3D14%2C%20label%3D%22cornerSubPix%22%2C%20zorder%3D5)%0A%20%20%20%20_ax_3w.plot(_sd_pt%5B0%5D%2C%20_sd_pt%5B1%5D%2C%20%22bo%22%2C%20markersize%3D14%2C%20label%3D%22Saddle%20fit%22%2C%20zorder%3D5)%0A%20%20%20%20_ax_3w.annotate(%0A%20%20%20%20%20%20%20%20%22%22%2C%20xy%3D(_sp_pt%5B0%5D%2C%20_sp_pt%5B1%5D)%2C%20xytext%3D(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D)%2C%0A%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D%22lime%22%2C%20lw%3D2)%2C%0A%20%20%20%20)%0A%20%20%20%20_dx_sp%20%3D%20_sp_pt%5B0%5D%20-%20_px_pt%5B0%5D%0A%20%20%20%20_dy_sp%20%3D%20_sp_pt%5B1%5D%20-%20_px_pt%5B1%5D%0A%20%20%20%20_ax_3w.text(%0A%20%20%20%20%20%20%20%20_px_pt%5B0%5D%20-%205%2C%20_px_pt%5B1%5D%20-%203%2C%0A%20%20%20%20%20%20%20%20f%22SubPix%3A%20(%7B_dx_sp%3A%2B.3f%7D%2C%20%7B_dy_sp%3A%2B.3f%7D)%22%2C%0A%20%20%20%20%20%20%20%20color%3D%22lime%22%2C%20fontsize%3D8%2C%20fontweight%3D%22bold%22%2C%0A%20%20%20%20%20%20%20%20bbox%3Ddict(boxstyle%3D%22round%22%2C%20facecolor%3D%22black%22%2C%20alpha%3D0.7)%2C%0A%20%20%20%20)%0A%20%20%20%20_ax_3w.annotate(%0A%20%20%20%20%20%20%20%20%22%22%2C%20xy%3D(_sd_pt%5B0%5D%2C%20_sd_pt%5B1%5D)%2C%20xytext%3D(_px_pt%5B0%5D%2C%20_px_pt%5B1%5D)%2C%0A%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D%22dodgerblue%22%2C%20lw%3D2)%2C%0A%20%20%20%20)%0A%20%20%20%20_dx_sd%20%3D%20_sd_pt%5B0%5D%20-%20_px_pt%5B0%5D%0A%20%20%20%20_dy_sd%20%3D%20_sd_pt%5B1%5D%20-%20_px_pt%5B1%5D%0A%20%20%20%20_ax_3w.text(%0A%20%20%20%20%20%20%20%20_px_pt%5B0%5D%20-%205%2C%20_px_pt%5B1%5D%20%2B%204%2C%0A%20%20%20%20%20%20%20%20f%22Saddle%3A%20(%7B_dx_sd%3A%2B.3f%7D%2C%20%7B_dy_sd%3A%2B.3f%7D)%22%2C%0A%20%20%20%20%20%20%20%20color%3D%22dodgerblue%22%2C%20fontsize%3D8%2C%20fontweight%3D%22bold%22%2C%0A%20%20%20%20%20%20%20%20bbox%3Ddict(boxstyle%3D%22round%22%2C%20facecolor%3D%22black%22%2C%20alpha%3D0.7)%2C%0A%20%20%20%20)%0A%20%20%20%20_dist%20%3D%20np.sqrt((_sp_pt%5B0%5D%20-%20_sd_pt%5B0%5D)%20**%202%20%2B%20(_sp_pt%5B1%5D%20-%20_sd_pt%5B1%5D)%20**%202)%0A%20%20%20%20_ax_3w.set_title(%0A%20%20%20%20%20%20%20%20f%22Three-way%20corner%20comparison%5CnSubPix%20vs%20Saddle%20distance%3A%20%7B_dist%3A.4f%7D%20px%22%0A%20%20%20%20)%0A%20%20%20%20_ax_3w.legend(fontsize%3D9%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20gradient_visualization(%0A%20%20%20%20all_corners_subpix%2C%0A%20%20%20%20detection_results%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20from%20scipy.ndimage%20import%20sobel%20as%20_sobel%0A%0A%20%20%20%20_view_idx_g%20%3D%20next(_i%20for%20_i%2C%20_d%20in%20enumerate(detection_results)%20if%20_d)%0A%20%20%20%20_img_g%20%3D%20simulated_images%5B_view_idx_g%5D.astype(np.float64)%0A%20%20%20%20_corners_g%20%3D%20all_corners_subpix%5B_view_idx_g%5D%0A%0A%20%20%20%20_gx%20%3D%20_sobel(_img_g%2C%20axis%3D1)%0A%20%20%20%20_gy%20%3D%20_sobel(_img_g%2C%20axis%3D0)%0A%0A%20%20%20%20_corner_idx_g%20%3D%20len(_corners_g)%20%2F%2F%202%0A%20%20%20%20_cp%20%3D%20_corners_g%5B_corner_idx_g%5D%0A%20%20%20%20_zoom_g%20%3D%2012%0A%20%20%20%20_x_lo_g%20%3D%20max(0%2C%20int(_cp%5B0%5D)%20-%20_zoom_g)%0A%20%20%20%20_x_hi_g%20%3D%20min(_img_g.shape%5B1%5D%2C%20int(_cp%5B0%5D)%20%2B%20_zoom_g)%0A%20%20%20%20_y_lo_g%20%3D%20max(0%2C%20int(_cp%5B1%5D)%20-%20_zoom_g)%0A%20%20%20%20_y_hi_g%20%3D%20min(_img_g.shape%5B0%5D%2C%20int(_cp%5B1%5D)%20%2B%20_zoom_g)%0A%0A%20%20%20%20_fig_grad%2C%20_ax_grad%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(7%2C%207))%0A%20%20%20%20_ax_grad.imshow(%0A%20%20%20%20%20%20%20%20_img_g%5B_y_lo_g%3A_y_hi_g%2C%20_x_lo_g%3A_x_hi_g%5D%2C%20cmap%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20extent%3D%5B_x_lo_g%2C%20_x_hi_g%2C%20_y_hi_g%2C%20_y_lo_g%5D%2C%20interpolation%3D%22nearest%22%2C%0A%20%20%20%20)%0A%20%20%20%20_step%20%3D%202%0A%20%20%20%20for%20_y%20in%20range(_y_lo_g%2C%20_y_hi_g%2C%20_step)%3A%0A%20%20%20%20%20%20%20%20for%20_x%20in%20range(_x_lo_g%2C%20_x_hi_g%2C%20_step)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_gx_val%20%3D%20_gx%5B_y%2C%20_x%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_gy_val%20%3D%20_gy%5B_y%2C%20_x%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_mag%20%3D%20np.sqrt(_gx_val%20**%202%20%2B%20_gy_val%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_mag%20%3E%2010%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_scale%20%3D%201.5%20%2F%20max(_mag%2C%201e-6)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_ax_grad.arrow(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_x%2C%20_y%2C%20_gx_val%20*%20_scale%2C%20_gy_val%20*%20_scale%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20head_width%3D0.3%2C%20head_length%3D0.2%2C%20fc%3D%22cyan%22%2C%20ec%3D%22cyan%22%2C%20alpha%3D0.8%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20_ax_grad.plot(_cp%5B0%5D%2C%20_cp%5B1%5D%2C%20%22r*%22%2C%20markersize%3D15%2C%20label%3D%22Sub-pixel%20corner%22)%0A%20%20%20%20_ax_grad.legend()%0A%20%20%20%20_ax_grad.set_title(%0A%20%20%20%20%20%20%20%20%22Gradient%20Vectors%20Near%20Corner%5Cn(tangent%20to%20edges%2C%20perpendicular%20to%20corner%20direction)%22%0A%20%20%20%20)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20corner_difference_histogram(%0A%20%20%20%20all_corners_saddle%2C%0A%20%20%20%20all_corners_subpix%2C%0A%20%20%20%20detection_results%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A)%3A%0A%20%20%20%20_dists%20%3D%20%5B%5D%0A%20%20%20%20_dx_all%2C%20_dy_all%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20for%20_det%2C%20_sp%2C%20_sd%20in%20zip(detection_results%2C%20all_corners_subpix%2C%20all_corners_saddle)%3A%0A%20%20%20%20%20%20%20%20if%20_det%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_d%20%3D%20np.sqrt(np.sum((_sp%20-%20_sd)%20**%202%2C%20axis%3D1))%0A%20%20%20%20%20%20%20%20%20%20%20%20_dists.extend(_d)%0A%20%20%20%20%20%20%20%20%20%20%20%20_dx_all.extend(_sd%5B%3A%2C%200%5D%20-%20_sp%5B%3A%2C%200%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_dy_all.extend(_sd%5B%3A%2C%201%5D%20-%20_sp%5B%3A%2C%201%5D)%0A%0A%20%20%20%20_dists%20%3D%20np.array(_dists)%0A%20%20%20%20_dx_all%20%3D%20np.array(_dx_all)%0A%20%20%20%20_dy_all%20%3D%20np.array(_dy_all)%0A%0A%20%20%20%20_fig_hist%2C%20(_ax_h%2C%20_ax_sc)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205))%0A%20%20%20%20_p99%20%3D%20np.percentile(_dists%2C%2099)%0A%20%20%20%20_n_clip%20%3D%20int(np.sum(_dists%20%3E%20_p99))%0A%20%20%20%20_ax_h.hist(_dists%2C%20bins%3D40%2C%20edgecolor%3D%22black%22%2C%20alpha%3D0.7%2C%20color%3D%22steelblue%22)%0A%20%20%20%20_ax_h.axvline(_dists.mean()%2C%20color%3D%22red%22%2C%20linestyle%3D%22--%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Mean%3A%20%7B_dists.mean()%3A.4f%7D%20px%22)%0A%20%20%20%20_ax_h.axvline(np.median(_dists)%2C%20color%3D%22orange%22%2C%20linestyle%3D%22--%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Median%3A%20%7Bnp.median(_dists)%3A.4f%7D%20px%22)%0A%20%20%20%20_ax_h.set_xlabel(%22Distance%20between%20cornerSubPix%20and%20Saddle%20(px)%22)%0A%20%20%20%20_ax_h.set_ylabel(%22Count%22)%0A%20%20%20%20_title%20%3D%20%22Refinement%20Differences%22%0A%20%20%20%20if%20_n_clip%20%3E%200%3A%0A%20%20%20%20%20%20%20%20_title%20%2B%3D%20f%22%20(%7B_n_clip%7D%20pts%20beyond%20p99%20clipped)%22%0A%20%20%20%20_ax_h.set_title(_title)%0A%20%20%20%20_ax_h.set_xlim(0%2C%20_p99)%0A%20%20%20%20_ax_h.legend()%0A%0A%20%20%20%20_ax_sc.scatter(_dx_all%2C%20_dy_all%2C%20s%3D3%2C%20alpha%3D0.5%2C%20color%3D%22steelblue%22)%0A%20%20%20%20_ax_sc.set_xlabel(%22dx%3A%20Saddle%20%E2%88%92%20SubPix%20(px)%22)%0A%20%20%20%20_ax_sc.set_ylabel(%22dy%3A%20Saddle%20%E2%88%92%20SubPix%20(px)%22)%0A%20%20%20%20_ax_sc.set_title(%22Per-Corner%20Offset%3A%20Saddle%20vs%20cornerSubPix%22)%0A%20%20%20%20_ax_sc.set_aspect(%22equal%22)%0A%20%20%20%20_ax_sc.axhline(0%2C%20color%3D%22gray%22%2C%20linewidth%3D0.5)%0A%20%20%20%20_ax_sc.axvline(0%2C%20color%3D%22gray%22%2C%20linewidth%3D0.5)%0A%20%20%20%20_lim%20%3D%20max(np.percentile(np.abs(_dx_all)%2C%2099)%2C%20np.percentile(np.abs(_dy_all)%2C%2099))%0A%20%20%20%20_ax_sc.set_xlim(-_lim%2C%20_lim)%0A%20%20%20%20_ax_sc.set_ylim(-_lim%2C%20_lim)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20noise_controls(mo)%3A%0A%20%20%20%20sigma_sl%20%3D%20mo.ui.slider(0%2C%2025%2C%20step%3D1%2C%20value%3D0%2C%20label%3D%22Image%20noise%20%CF%83%20(grey%20levels)%22)%0A%20%20%20%20return%20(sigma_sl%2C)%0A%0A%0A%40app.cell%0Adef%20noise_sensitivity(%0A%20%20%20%20K_true%2C%0A%20%20%20%20corner_subpix%2C%0A%20%20%20%20dist_true%2C%0A%20%20%20%20gaussian_filter%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obj_points_3d%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20poses%2C%0A%20%20%20%20project_points%2C%0A%20%20%20%20saddle_refine_corners%2C%0A%20%20%20%20sigma_sl%2C%0A%20%20%20%20simulated_images%2C%0A)%3A%0A%20%20%20%20_sigma%20%3D%20float(sigma_sl.value)%0A%20%20%20%20_margin%20%3D%2015%0A%20%20%20%20_n_trials%20%3D%205%0A%0A%20%20%20%20%23%20Interior%20corners%20only%3A%20exclude%20first%2Flast%20row%20and%20column%20of%20the%208%C3%975%20grid%0A%20%20%20%20_n_cols_grid%20%3D%208%0A%20%20%20%20_n_rows_grid%20%3D%205%0A%20%20%20%20_interior%20%3D%20np.zeros(len(obj_points_3d)%2C%20dtype%3Dbool)%0A%20%20%20%20for%20_row%20in%20range(1%2C%20_n_rows_grid%20-%201)%3A%0A%20%20%20%20%20%20%20%20for%20_col%20in%20range(1%2C%20_n_cols_grid%20-%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_interior%5B_row%20*%20_n_cols_grid%20%2B%20_col%5D%20%3D%20True%0A%0A%20%20%20%20_errs_sp%2C%20_errs_sd%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20_total_corners%20%3D%200%0A%0A%20%20%20%20for%20_vi%2C%20(_img%2C%20(_rvec%2C%20_tvec))%20in%20enumerate(zip(simulated_images%2C%20poses))%3A%0A%20%20%20%20%20%20%20%20_gt%20%3D%20project_points(obj_points_3d%2C%20_rvec%2C%20_tvec%2C%20K_true%2C%20dist_true)%0A%20%20%20%20%20%20%20%20_h%2C%20_w%20%3D%20_img.shape%5B%3A2%5D%0A%20%20%20%20%20%20%20%20_mask%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20(_gt%5B%3A%2C%200%5D%20%3E%3D%20_margin)%20%26%20(_gt%5B%3A%2C%200%5D%20%3C%20_w%20-%20_margin)%20%26%0A%20%20%20%20%20%20%20%20%20%20%20%20(_gt%5B%3A%2C%201%5D%20%3E%3D%20_margin)%20%26%20(_gt%5B%3A%2C%201%5D%20%3C%20_h%20-%20_margin)%20%26%0A%20%20%20%20%20%20%20%20%20%20%20%20_interior%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_gt_valid%20%3D%20_gt%5B_mask%5D%0A%20%20%20%20%20%20%20%20if%20len(_gt_valid)%20%3C%205%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20_total_corners%20%2B%3D%20len(_gt_valid)%0A%0A%20%20%20%20%20%20%20%20for%20_t%20in%20range(_n_trials)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_rng%20%3D%20np.random.RandomState(_vi%20*%2010000%20%2B%20_t%20*%201000%20%2B%20int(_sigma%20*%20100))%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_sigma%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_noisy%20%3D%20np.clip(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_img.astype(np.float64)%20%2B%20_rng.randn(*_img.shape)%20*%20_sigma%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%200%2C%20255%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20).astype(np.uint8)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_noisy%20%3D%20_img.copy()%0A%20%20%20%20%20%20%20%20%20%20%20%20_blur%20%3D%20gaussian_filter(_noisy.astype(np.float64)%2C%20sigma%3D0.5)%0A%20%20%20%20%20%20%20%20%20%20%20%20_blur_u8%20%3D%20np.clip(_blur%2C%200%2C%20255).astype(np.uint8)%0A%20%20%20%20%20%20%20%20%20%20%20%20_seeds%20%3D%20np.round(_gt_valid).astype(np.float64)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ref_sp%20%3D%20corner_subpix(_blur_u8%2C%20_seeds%2C%20win_size%3D5)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ref_sd%20%3D%20saddle_refine_corners(_blur_u8%2C%20_seeds%2C%20win_size%3D5)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20_err_sp%20%3D%20np.sqrt((_ref_sp%5B%3A%2C%200%5D%20-%20_gt_valid%5B%3A%2C%200%5D)%20**%202%20%2B%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(_ref_sp%5B%3A%2C%201%5D%20-%20_gt_valid%5B%3A%2C%201%5D)%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20_err_sd%20%3D%20np.sqrt((_ref_sd%5B%3A%2C%200%5D%20-%20_gt_valid%5B%3A%2C%200%5D)%20**%202%20%2B%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(_ref_sd%5B%3A%2C%201%5D%20-%20_gt_valid%5B%3A%2C%201%5D)%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20_errs_sp.extend(_err_sp)%0A%20%20%20%20%20%20%20%20%20%20%20%20_errs_sd.extend(_err_sd)%0A%0A%20%20%20%20_errs_sp%20%3D%20np.array(_errs_sp)%0A%20%20%20%20_errs_sd%20%3D%20np.array(_errs_sd)%0A%20%20%20%20_rms_sp%20%3D%20float(np.sqrt(np.mean(_errs_sp%20**%202)))%20if%20len(_errs_sp)%20else%200.0%0A%20%20%20%20_rms_sd%20%3D%20float(np.sqrt(np.mean(_errs_sd%20**%202)))%20if%20len(_errs_sd)%20else%200.0%0A%20%20%20%20_n_views%20%3D%20len(simulated_images)%0A%0A%20%20%20%20_fig_ns%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%205))%0A%20%20%20%20_xmax%20%3D%20max(0.15%2C%20float(max(_errs_sp.max()%2C%20_errs_sd.max()))%20*%201.05)%20if%20len(_errs_sp)%20else%200.5%0A%20%20%20%20_bins%20%3D%20np.linspace(0%2C%20_xmax%2C%2045)%0A%20%20%20%20_ax.hist(_errs_sp%2C%20bins%3D_bins%2C%20color%3D%22green%22%2C%20alpha%3D0.55%2C%20edgecolor%3D%22white%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%20label%3Df%22cornerSubPix%20%20(RMS%20%3D%20%7B_rms_sp%3A.4f%7D%20px)%22)%0A%20%20%20%20_ax.hist(_errs_sd%2C%20bins%3D_bins%2C%20color%3D%22steelblue%22%2C%20alpha%3D0.55%2C%20edgecolor%3D%22white%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%20label%3Df%22Saddle%20fit%20%20(RMS%20%3D%20%7B_rms_sd%3A.4f%7D%20px)%22)%0A%20%20%20%20_ax.set_xlabel(%22Corner%20position%20error%20(pixels)%22)%0A%20%20%20%20_ax.set_ylabel(%22Density%22)%0A%20%20%20%20_ax.set_title(%0A%20%20%20%20%20%20%20%20f%22Sub-pixel%20refinement%20error%20at%20%CF%83%20%3D%20%7B_sigma%3A.0f%7D%20%20%22%0A%20%20%20%20%20%20%20%20f%22(%7B_total_corners%7D%20corners%20%C3%97%20%7B_n_views%7D%20views%20%C3%97%20%7B_n_trials%7D%20trials)%22%0A%20%20%20%20)%0A%20%20%20%20_ax.set_xlim(0%2C%20_xmax)%0A%20%20%20%20_ax.legend(fontsize%3D10)%0A%20%20%20%20_ax.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20_expl_ns%20%3D%20mo.md(f%22%22%22%0A%20%20%20%20%23%23%23%20Noise%20Sensitivity%20of%20Sub-Pixel%20Corner%20Detection%0A%0A%20%20%20%20What's%20a%20bit%20more%20interesting%20is%20to%20compare%20the%20methods%20across%20all%20poses%20and%20see%20how%20we%20err%20with%20respect%20to%20the%20ground%20truth.%20And%20even%20more%20so%2C%20how%20noise%20affects%20this.%20Drag%20the%20**%CF%83%20slider**%20to%20add%20Gaussian%20noise%20and%20compare%20the%20two%20refinement%20methods%0A%20%20%20%20across%20**all%20%7B_n_views%7D%20board%20poses**%20(various%20tilts%20and%20rotations)%3A%0A%0A%20%20%20%20-%20**cornerSubPix**%20(green)%20%E2%80%94%20gradient%20perpendicularity%2C%20general-purpose%0A%20%20%20%20-%20**Saddle%20fit**%20(blue)%20%E2%80%94%20quadratic%20surface%20model%2C%20checkerboard-specific%0A%0A%20%20%20%20%7C%20Method%20%7C%20RMS%20error%20at%20%CF%83%20%3D%20%7B_sigma%3A.0f%7D%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20cornerSubPix%20%7C%20**%7B_rms_sp%3A.4f%7D%20px**%20%7C%0A%20%20%20%20%7C%20Saddle%20fit%20%7C%20**%7B_rms_sd%3A.4f%7D%20px**%20%7C%0A%0A%20%20%20%20We%20can%20see%20that%20in%20this%20case%20the%20saddle%20fit%20performed%20a%20bit%20better%2C%20probably%20due%20to%20having%20useful%20prior%20information.%20%0A%20%20%20%20%22%22%22)%0A%20%20%20%20mo.vstack(%5B_expl_ns%2C%20sigma_sl%2C%20_fig_ns%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20summary(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Summary%0A%0A%20%20%20%20Both%20methods%20converge%20to%20essentially%20the%20same%20location%20on%20clean%20images%20%E2%80%94%20the%20advantage%20of%0A%20%20%20%20saddle%20fitting%20shows%20up%20mainly%20under%20noise%2C%20where%20the%20checker%20geometry%20gives%20it%20a%20slight%0A%20%20%20%20edge.%20On%20real%20cameras%20with%20lens%20blur%20and%20sensor%20noise%20the%20difference%20widens.%0A%0A%20%20%20%20**Next**%3A%20next%20time%20I%20think%20(like%20I%20did%20multiple%20times%20before)%20that%20we%20are%20ready%20to%20explor%20Zhang's%20Calibration%20method%20and%20get%20this%20saga%20over%20with.%0A%20%20%20%20%22%22%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