Skip to content

Instantly share code, notes, and snippets.

@hmaarrfk

hmaarrfk/otsu_failer.ipynb

Last active November 3, 2018 14:22
Show Gist options
  • Save hmaarrfk/4afae1cfded1d78e44c9e4f58285d552 to your computer and use it in GitHub Desktop.
Save hmaarrfk/4afae1cfded1d78e44c9e4f58285d552 to your computer and use it in GitHub Desktop.
Download ZIP
Summary of the OTSU Failer
Display the source blob
Display the rendered blob
Raw
otsu_failer.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"# !conda install numpy -y\n",
"import numpy as np\n",
"\n",
"# OTSU Failer on 32 bit can be summarized as this\n",
"pixel = (3, 5)\n",
"\n",
"pixel = (19, 18)\n",
"\n",
"if pixel == (3, 5):\n",
" possible_values = (41, 81)\n",
" img = np.array([[ 53, 41, 167],\n",
" [ 30, 81, 106],\n",
" [ 63, 147, 151]], dtype=np.uint8)\n",
"elif pixel == (19, 18):\n",
" poxxible_values = (141, 172)\n",
" img = np.array([[ 78, 214, 61],\n",
" [229, 104, 141],\n",
" [ 87, 172, 224]], dtype=np.uint8)\n",
" \n",
"selem = np.array([[0, 1, 0],\n",
" [1, 1, 1],\n",
" [0, 1, 0]], dtype=np.uint8)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[104 141 172 214 229]\n"
]
}
],
"source": [
"# in 0.14.1 after passing through OTSU, this would return 81 on 64 bit, and 41 on 32 bit\n",
"# for the [1, 1] pixel.\n",
"# To understand this, lets go through the OTSU algorithm at that pixel\n",
"\n",
"img_center = img[selem==1]\n",
"img_center.sort()\n",
"\n",
"histo, index = np.histogram(img_center, bins=img_center[-1]+1, range=(-0.5,img_center[-1]+0.5), density=False)\n",
"index = index[:-1]+0.5\n",
"print(img_center)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"229.0\n",
"1\n",
"[1 1 1 1 1]\n",
"230\n"
]
}
],
"source": [
"# Just proving that I know how to use histogram (fun fact: I didn't)\n",
"print(index[-1])\n",
"print(histo[-1])\n",
"print(histo[img_center])\n",
"print(len(histo))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/srv/conda/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
" \"\"\"\n",
"/srv/conda/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide\n",
" \n"
]
}
],
"source": [
"# Following wikipedia's notation\n",
"w0 = np.cumsum(histo)\n",
"w1 = w0[-1] - w0\n",
"# \n",
"mu0 = np.cumsum(histo * index) / w0\n",
"mu0[w0 == 0] = 0\n",
"muT = np.sum(histo * index)\n",
"mu1 = (muT - w0 * mu0) / w1\n",
"mu1[w1 == 0] = 0\n",
"\n",
"# The algorithm in 0.14.1 didn't take care of the true divide case at the very end either\n",
"# This wasn't actually the issue though"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"sigma2 = w0 * w1 * (mu0 - mu1)**2"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 0. 0. 0. 0. 0. 0. 0. 0.\n",
" 28900. 28900. 28900. 28900. 28900. 28900. 28900. 28900.\n",
" 28900. 28900. 28900. 28900. 28900. 28900. 28900. 28900.\n",
" 28900. 28900. 28900. 28900. 28900. 28900. 28900. 28900.\n",
" 28900. 28900. 28900. 28900. 28900. 28900. 28900. 28900.\n",
" 28900. 28900. 28900. 28900. 28900. 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 20306.25 20306.25\n",
" 20306.25 20306.25 20306.25 20306.25 20306.25 20306.25 20306.25 20306.25\n",
" 20306.25 20306.25 20306.25 20306.25 20306.25 0. ]\n"
]
}
],
"source": [
"print(sigma2)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5 40837.5\n",
" 40837.5 40837.5 40837.5 40837.5 40837.5]\n"
]
}
],
"source": [
"print(sigma2[poxxible_values[0]:poxxible_values[1]+1])\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"# Clearly, the issue is that this is a pathological case, where depending on the\n",
"# rounding error in the particular hardware implementation and software optimization\n",
"# the algorithm will either find 41, or 81 optimal"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"i=104, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=105, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=106, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=107, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=108, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=109, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=110, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=111, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=112, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=113, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=114, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=115, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=116, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=117, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=118, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=119, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=120, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=121, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=122, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=123, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=124, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=125, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=126, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=127, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=128, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=129, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=130, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=131, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=132, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=133, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=134, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=135, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=136, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=137, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=138, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=139, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=140, mu1=104.0 mu2=756.0 pop=5, new_q1 = 1, t=-340.0, sigma_b=28900.0\n",
"i=141, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=142, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=143, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=144, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=145, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=146, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=147, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=148, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=149, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=150, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=151, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=152, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=153, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=154, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=155, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=156, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=157, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=158, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=159, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=160, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=161, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=162, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=163, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=164, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=165, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=166, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=167, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=168, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=169, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=170, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=171, mu1=245.0 mu2=615.0 pop=5, new_q1 = 2, t=-495.0, sigma_b=40837.5\n",
"i=172, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=173, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=174, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=175, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=176, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=177, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=178, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=179, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=180, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=181, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=182, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=183, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=184, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=185, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=186, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=187, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=188, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=189, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=190, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=191, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=192, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=193, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=194, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=195, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=196, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=197, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=198, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=199, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=200, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=201, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=202, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=203, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=204, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=205, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=206, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=207, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=208, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=209, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=210, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=211, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=212, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=213, mu1=417.0 mu2=443.0 pop=5, new_q1 = 3, t=-495.0, sigma_b=40837.5\n",
"i=214, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=215, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=216, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=217, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=218, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=219, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=220, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=221, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=222, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=223, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=224, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=225, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=226, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=227, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n",
"i=228, mu1=631.0 mu2=229.0 pop=5, new_q1 = 4, t=-285.0, sigma_b=20306.25\n"
]
}
],
"source": [
"# Here is the Cython implementation written in python to show you what happens\n",
"# You may choose to add the print statement in cython,\n",
"# just remember to use `with gil:`\n",
"pop = len(img_center)\n",
"mu = np.sum(histo*index)\n",
"max_i = 0\n",
"q1 = histo[0]\n",
"mu1 = 0.\n",
"max_sigma_b = 0.\n",
"\n",
"for i in range(1, len(histo)):\n",
" P = histo[i]\n",
" new_q1 = q1 + P\n",
" if new_q1 == pop:\n",
" break\n",
" if new_q1 > 0:\n",
" # Rearrange the formula so that you only have to divide once per loop\n",
" # All other operations are additions and multiplications\n",
" mu1 = mu1 + i * P\n",
" mu2 = mu - mu1\n",
" t = (pop - new_q1) * mu1 - mu2 * new_q1\n",
" sigma_b = (t * t) / (new_q1 * (pop - new_q1))\n",
" print(f\"i={i}, mu1={mu1} mu2={mu2} pop={pop}, new_q1 = {new_q1}, t={t}, sigma_b={sigma_b}\")\n",
" if sigma_b > max_sigma_b:\n",
" max_sigma_b = sigma_b\n",
" max_i = i\n",
" q1 = new_q1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment