C语言库函数：https://blog.csdn.net/weixin_44793491/article/details/107644666
原文见：https://blog.csdn.net/wxplol/article/details/81239896
C++ 实现：https://blog.csdn.net/cxf7394373/article/details/45155053

步骤1 计算区域像素积分和和积分平方和
步骤2 计算标准差，标准差的计算方法为: std = sqrt((sqdiff - diff*diff / area) / (area - 1))
步骤3 sauvola二值化算法的阈值为T = mean*(1 + k*((std / 128) - 1)).

使用Cython重写代码，优化了运行时间。
jupyter notebook :%load_ext Cython
编译命令：python setup.py build_ext --inplace

%%cython -a --cplus --compile-args=/openmp
import cv2
cimport cython
cimport numpy as np
import numpy as np
from cython.parallel import parallel, prange
from libc.math cimport pow
from libc.math cimport sqrt
from libc.stdlib cimport exit


@cython.boundscheck(False)
@cython.wraparound(False)
def Sauvola(np.ndarray[np.uint8_t, ndim=3] src, float k=0.1, int kernerl=31):
    '''
    计算图像的积分和平方积分
    :param ndarray[np.uint8,ndim=3] src  输入图像，为三通道BGR图像
    :param k:float:                      修正参数,一般0<k<1
    :param kernerl:int                   窗口大小
    :return:ndarray[np.uint8,ndim=2] img：二值化输出图像，为单通道图像；
    '''
    cdef Py_ssize_t height = src.shape[0]
    cdef Py_ssize_t width  = src.shape[1]
    cdef Py_ssize_t rows = height
    cdef Py_ssize_t cols = width
    cdef Py_ssize_t h = 0, w = 0
    cdef Py_ssize_t row=0, col=0
    cdef int sum_  = 0
    cdef double sqrt_sum  = 0
    cdef int whalf = 31 >> 1
    
    integ_sum  = np.zeros((height, width), dtype=np.uint32)
    integ_sqrt = integ_sum.copy()
    img = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
    cdef unsigned char[:,::1] img_view = img
    cdef unsigned long[:,::1] sum_view = integ_sum
    cdef unsigned long[:,::1] sqrt_view = integ_sqrt
    with nogil:
        for h in range(height):
            sum_ = 0
            sqrt_sum = 0
            for w in range(width):
                sum_ += img_view[h][w]
                sqrt_sum += sqrt(img_view[h][w])

                if h == 0:
                    sum_view[h][w] = sum_
                    sqrt_view[h][w] = <int>sqrt_sum
                else:
                    sum_view[h][w] = sum_ + sum_view[h - 1][w]
                    sqrt_view[h][w] = <int>(sqrt_sum + sqrt_view[h - 1][w])
    #==============计算积分图和积分平方和图完成========================

    diff      = np.zeros((height, width), np.float32)
    sqrt_diff = np.zeros((height, width), np.float32)
    mean      = np.zeros((height, width), np.float32)
    threshold = np.zeros((height, width), np.float32)
    std       = np.zeros((height, width), np.float32)
    cdef float[:,::1] diff_view = diff
    cdef float[:,::1] sdiff_view = sqrt_diff
    cdef float[:,::1] mean_view = mean
    cdef float[:,::1] thresh_view = threshold
    cdef float[:,::1] std_view = std
    cdef int xmin  = 0
    cdef int ymin  = 0
    cdef int xmax  = 0
    cdef int ymax  = 0
    cdef int area  = 0
    
    with nogil:
        for row in range(rows):
            for col in range(cols):
                xmin = max(0, row - whalf)
                ymin = max(0, col - whalf)
                xmax = min(rows - 1, row + whalf)
                ymax = min(cols - 1, col + whalf)
                area = (xmax - xmin + 1) * (ymax - ymin + 1)
                if area <= 0:
                    exit(1)

                if xmin == 0 and ymin == 0:
                    diff_view[row, col] = sum_view[xmax, ymax]
                    sdiff_view[row, col] = sqrt_view[xmax, ymax]
                elif xmin > 0 and ymin == 0:
                    diff_view[row, col] = sum_view[xmax, ymax] - sum_view[xmin - 1, ymax]
                    sdiff_view[row, col] = sqrt_view[xmax, ymax] - sqrt_view[xmin - 1, ymax]
                elif xmin == 0 and ymin > 0:
                    diff_view[row, col] = sum_view[xmax, ymax] - sum_view[xmax, ymax - 1]
                    sdiff_view[row, col] = sqrt_view[xmax, ymax] - sqrt_view[xmax, ymax - 1]
                else:
                    diff_view[row, col] = (sum_view[xmax, ymax] + sum_view[xmin - 1, ymin - 1]) - (sum_view[xmax, ymin - 1] + sum_view[xmin - 1, ymax])
                    sdiff_view[row, col] = (sqrt_view[xmax, ymax] + sqrt_view[xmin - 1, ymin - 1]) - (sqrt_view[xmax, ymin - 1] + sqrt_view[xmin - 1, ymax])


                mean_view[row, col] = diff_view[row, col] / area
                std_view[row, col] = sqrt((sdiff_view[row, col] - sqrt(diff_view[row, col]) / area) / (area - 1))
                thresh_view[row, col] = mean_view[row, col] * (1 + k * ((std_view[row, col] / 128) - 1))
                if img_view[row, col] < thresh_view[row, col]:
                    img_view[row, col] = 0
                else:
                    img_view[row, col] = 255
    return img

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np


ext_modules = [Extension('Sauvola', 
                         sources=['Sauvola.pyx'],
                         extra_compile_args=['/openmp'],
                         language='c')]

setup(
    name = 'Sauvola',
    cmdclass = {
    
    'build_ext': build_ext},
    ext_modules = ext_modules,
    include_dirs=[np.get_include()]
)

import time
import numpy as np
import cv2
from Sauvola import Sauvola



if __name__ == '__main__':
    t0 = time.time()
    # 请输入RGB三通道图，不支持sRGB
    input_fn = "E:/Desktop/sauvola_test.jpg"
    img = cv2.imread(input_fn)
    print("[img]",type(img),img.dtype, img.shape)
    dstimg = Sauvola(img)
    cv2.imwrite("E:/Desktop/sauvola_test_result.png", dstimg)
    t1 = time.time()
    print("Total time is:",t1-t0) 
"""
[img] <class 'numpy.ndarray'> uint8 (256, 1518, 3)
Total time is: 0.023972749710083008
"""