DIP_UN_06

---

课程目标

掌握图像增强的一般方法。重点掌握灰度变换、直方图的生成以及各类滤波方法（包括中值滤波、高斯滤波、边界保持滤波、同态滤波等）；了解常用的锐化算子；对简单的噪声图像进行处理。

一、图像基本运算

逐点运算示例：

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt

img = cv.imread("D:\\DIP_Photo\\Lenna_RGB.tif", 0)

height = img.shape[0]
width = img.shape[1]

img01 = np.zeros((height, width, 3), np.uint8)
img02 = np.zeros((height, width, 3), np.uint8)
img03 = np.zeros((height, width, 3), np.uint8)
img04 = np.zeros((height, width, 3), np.uint8)
img05 = np.zeros((height, width, 3), np.uint8)

for i in range(height):
    for j in range(width):
        if int(img[i, j] + 50) > 255:
            gray = 255
        else:
            gray = int(img[i, j] + 50)
        img01[i, j] = gray

for i in range(height):
    for j in range(width):
        if int(img[i, j] - 50) < 0:
            gray = 0
        else:
            gray = int(img[i, j] - 50)
        img02[i, j] = gray

for i in range(height):
    for j in range(width):
        if int(img[i, j] * 2) > 255:
            gray = 255
        else:
            gray = int(img[i, j] * 2)
        img03[i, j] = gray

for i in range(height):
    for j in range(width):
        gray = int(img[i, j] * 0.8)
        img04[i, j] = gray

for i in range(height):
    for j in range(width):
        gray = 255 - img[i, j]
        img05[i, j] = gray

titles = ['Lenna', '+50', '-50', '*2', "*0.8", '255-']
images = [img, img01, img02, img03, img04, img05]


plt.figure(figsize=(12, 6))
for i in range(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

---

二、非线性运算

1.简单非线性变换

import cv2 as cv
import numpy as np

grayImage = cv.imread("D:\\DIP_Photo\\Lenna_RGB.tif", 0)

height = grayImage.shape[0]
width = grayImage.shape[1]

NLImage = np.zeros((height, width), np.uint8)

for i in range(height):
    for j in range(width):
        gray = int(grayImage[i, j])*int(grayImage[i, j]) / 255
        NLImage[i, j] = np.uint8(gray)

cv.imshow('Input', grayImage)
cv.imshow('NLImage', NLImage)

cv.waitKey(0)
cv.destroyAllWindows()

2.对数变换

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt


def log_plot(Cons):
    x = np.arange(0, 256, 0.01)
    y = Cons * np.log(1 + x)
    plt.plot(x, y, 'r', linewidth=1)
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.title('对数变换函数')
    plt.xlim(0, 255), plt.ylim(0, 255)
    plt.show()


C = 43

grayImage = cv.imread("D:\\DIP_Photo\\Lenna_RGB.tif", 0)

LogImage = C * np.log(1.0 + grayImage)
LogImage = np.uint8(LogImage + 0.5)

cv.imshow('Input', grayImage)
cv.imshow('LogImage', LogImage)
log_plot(C)

cv.waitKey(0)
cv.destroyAllWindows()

3.伽马变换

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt


def gamma_plot(c, v):
    x = np.arange(0, 256, 0.01)
    y = c*x**v
    plt.plot(x, y, 'r', linewidth=1)
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.title('伽马变换函数')
    plt.xlim([0, 255]), plt.ylim([0, 255])
    plt.show()


def gamma(img, c, v):
    lut = np.zeros(256, dtype=np.float32)
    for i in range(256):
        lut[i] = c * i ** v
    output_img = cv.LUT(img, lut)
    output_img = np.uint8(output_img+0.5)
    return output_img


grayImage = cv.imread("D:\\DIP_Photo\\Bubbles.tif", 0)
GammaImage = gamma(grayImage, 0.00000005, 4.0)

cv.imshow('Input', grayImage)
cv.imshow('GammaImage', GammaImage)
gamma_plot(0.00000005, 4.0)

cv.waitKey(0)
cv.destroyAllWindows()

---

三、直方图

1.直方图计算

import numpy as np
import cv2
import matplotlib.pyplot as plt


def my_hist(img, bins=256, normalized=False):
    """
    create a hist of uint8 image value range[0, 255]
    param: input img: grayscale image range[0, 255]
    param: input bins: bins for the image, range[0, 255]
    return hist and bins for the image, hist -> counts for all the values
    """
    # initializ a list for values and counts
    list1 = list([x, y] for x in np.arange(bins) for y in np.arange(bins) if y == 0)
    data = img.flatten()
    for i in range(img.size):
        list1[data[i]][1] += 1
    dst = np.array(list1)
    bins, hist= dst[:, 0], dst[:, 1]
    if normalized:
        hist = hist / img.size  #nomalized hist
    return hist, bins


# 直方图
img_1st = cv2.imread('D:\\DIP_Photo\\Fig0316(1)(top_left).tif', 0)
img_2nd = cv2.imread('D:\\DIP_Photo\\Fig0316(2)(2nd_from_top).tif', 0)
img_3rd = cv2.imread('D:\\DIP_Photo\\Fig0316(3)(third_from_top).tif', 0)
img_4th = cv2.imread('D:\\DIP_Photo\\Fig0316(4)(bottom_left).tif', 0)

img_list =['img_1st', 'img_2nd', 'img_3rd', 'img_4th']

fig = plt.figure(figsize=(15, 9))
for i in range(len(img_list)):
    ax = fig.add_subplot(3, 4, i+1)
    ax.imshow(eval(img_list[i]), cmap='gray', vmin=0, vmax=255)
    ax1 = fig.add_subplot(3, 4, i+5)
    hist, bins = my_hist(eval(img_list[i]), bins=256)
    ax1.bar(bins, hist), ax1.set_title('My Hist')
    ax2 = fig.add_subplot(3, 4, i+9)
    hist, bins = np.histogram(eval(img_list[i]), bins=256, range=[0, 256])
    ax2.bar(bins[:-1], hist), ax2.set_title('Numpy Hist')  # numpy 的返回bins 比hist多一点数，需要去掉最后一个

plt.tight_layout()
plt.show()

2.matplotlib直方图

import cv2 as cv
import matplotlib.pyplot as plt

src = cv.imread("D:\\DIP_Photo\\RENZhengfei.bmp")

b, g, r = cv.split(src)

plt.figure("RZF_hist")
plt.hist(b.ravel(), bins=256, density=1, facecolor='b', edgecolor='b', alpha=0.5)
plt.hist(g.ravel(), bins=256, density=1, facecolor='g', edgecolor='g', alpha=0.5)
plt.hist(r.ravel(), bins=256, density=1, facecolor='r', edgecolor='r', alpha=0.5)

plt.xlabel("x")
plt.ylabel("y")

cv.imshow("Renzhengfei", src)
plt.show()

cv.waitKey(0)
cv.destroyAllWindows()

3.opencv直方图

import cv2
import matplotlib.pyplot as plt
import matplotlib

src = cv2.imread("D:\\DIP_Photo\\JinChen01.jpg")

img_rgb = cv2.cvtColor(src, cv2.COLOR_BGR2RGB)

histb = cv2.calcHist([src], [0], None, [256], [0, 255])
histg = cv2.calcHist([src], [1], None, [256], [0, 255])
histr = cv2.calcHist([src], [2], None, [256], [0, 255])

matplotlib.rcParams['font.sans-serif'] = ['SimHei']

plt.subplot(121)
plt.imshow(img_rgb, 'gray')
plt.axis('off')
plt.title("(a)原始图像")

plt.subplot(122)
plt.plot(histb, color='b')
plt.plot(histg, color='g')
plt.plot(histr, color='r')
plt.xlabel("x")
plt.ylabel("y")
plt.title("(b)直方图曲线")
plt.show()

4.掩膜直方图

import cv2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

img = cv2.imread("D:\\DIP_Photo\\RENZhengfei.bmp")

img_rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

mask = np.zeros(img.shape[:2], np.uint8)
mask[80:380, 140:440] = 255
masked_img = cv2.bitwise_and(img_rgb, img_rgb, mask=mask)

hist_full = cv2.calcHist([img_rgb], [0], None, [256], [0, 256])  # 通道[0]-灰度图
hist_mask = cv2.calcHist([img_rgb], [0], mask, [256], [0, 256])  # 掩膜灰度图

plt.figure(figsize=(8, 6))
matplotlib.rcParams['font.sans-serif'] = ['SimHei']

plt.subplot(221)
plt.imshow(img_rgb, 'gray')
plt.axis('off')
plt.title("(a)原始图像")

plt.subplot(222)
plt.imshow(mask, 'gray')
plt.axis('off')
plt.title("(b)掩膜")

plt.subplot(223)
plt.imshow(masked_img, 'gray')
plt.axis('off')
plt.title("(c)图像掩膜处理")

plt.subplot(224)
plt.plot(hist_full)
plt.plot(hist_mask)
plt.title("(d)直方图曲线")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

5.全局直方图均衡和及自适应直方图均衡

import cv2 as cv
import numpy as np


# 灰度图像的全局直方图均衡化和及自适应直方图均衡化
img = cv.imread("D:\\DIP_Photo\\Girl.jpg", 0)
img1 = cv.equalizeHist(img)     # 全局直方图均衡化
clahe = cv.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))      # 自适应直方图均衡化
cll = clahe.apply(img)
res = np.hstack((img, img1, cll))

cv.imshow("Hist", res)
cv.waitKey()
cv.destroyAllWindows()


# 彩色图像全局直方图均衡化


def hisEqulColor1(img):
    # 将RGB图像转换到YCrCb空间中
    ycrcb = cv.cvtColor(img, cv.COLOR_BGR2YCR_CB)
    # 将YCrCb图像通道分离
    channels = cv.split(ycrcb)
    # 对第1个通道即亮度通道进行全局直方图均衡化并保存
    cv.equalizeHist(channels[0], channels[0])
    # 将处理后的通道和没有处理的两个通道合并，命名为ycrcb
    cv.merge(channels, ycrcb)
    # 将YCrCb图像转换回RGB图像
    cv.cvtColor(ycrcb, cv.COLOR_YCR_CB2BGR, img)
    return img


# 彩色图像进行自适应直方图均衡化，代码同上的地方不再添加注释
def hisEqulColor2(img):
    ycrcb = cv.cvtColor(img, cv.COLOR_BGR2YCR_CB)
    channels = cv.split(ycrcb)

    # 以下代码详细注释见官网：
    # https://docs.opencv.org/4.1.0/d5/daf/tutorial_py_histogram_equalization.html
    clahe = cv.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
    clahe.apply(channels[0], channels[0])

    cv.merge(channels, ycrcb)
    cv.cvtColor(ycrcb, cv.COLOR_YCR_CB2BGR, img)
    return img


img = cv.imread('D:\\DIP_Photo\\Einstein.jpg')
img = cv.resize(img, (480, 480))
img1 = img.copy()
img2 = img.copy()

res1 = hisEqulColor1(img1)
res2 = hisEqulColor2(img2)

res = np.hstack((img, res1, res2))
cv.imshow("Einstein", res)

cv.waitKey()
cv.destroyAllWindows()

6.局部直方图处理

import numpy as np
import cv2
import matplotlib.pyplot as plt


def my_hist(img, bins=256, normalized=False):
    """
    create a hist of uint8 image value range[0, 255]
    param: input img: grayscale image range[0, 255]
    param: input bins: bins for the image, range[0, 255]
    return hist and bins for the image, hist -> counts for all the values
    """
    # initializ a list for values and counts
    list1 = list([x, y] for x in np.arange(bins) for y in np.arange(bins) if y == 0)
    data = img.flatten()
    for i in range(img.size):
        list1[data[i]][1] += 1
    dst = np.array(list1)
    bins, hist = dst[:, 0], dst[:, 1]
    if normalized:
        hist = hist / img.size  # Nomalized hist
    return hist, bins


def my_calhist(img):
    """
    histogram equalization
    param: input img: uint8[0, 255] grayscale image
    return uint8[0, 255] grayscale image after histogram equalization
    """
    hist, bins = my_hist(img, bins=256, normalized=True)
    hist_cumsum = np.round(np.cumsum(hist) * 255).astype(int)

    height, width = img.shape[:2]
    img_dst = np.zeros([height, width], np.uint8)
    for h in range(height):
        for w in range(width):
            img_dst[h, w] = hist_cumsum[img[h, w]]  # dict 用[ ]
    img_dst = np.clip(img_dst, 0, 255).astype(np.uint8)
    return img_dst, hist_cumsum


# opencv实现的局部直方图处理
img_ori = cv2.imread('D:\\DIP_Photo\\Fig0326(a)(embedded_square_noisy_512).tif', 0)
img_transform, _ = my_calhist(img_ori)

plt.figure(figsize=(15, 6))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=255), plt.title('Original')
plt.subplot(1, 3, 2), plt.imshow(img_transform, cmap='gray', vmin=0, vmax=255), plt.title(f'Global Equalize Hist')

clahe = cv2.createCLAHE(clipLimit=255, tileGridSize=(3, 3))
img_transform = clahe.apply(img_ori)

plt.subplot(1, 3, 3), plt.imshow(img_transform, cmap='gray', vmin=0, vmax=255), plt.title(f'Local Equalize Hist')
plt.tight_layout()
plt.show()

7.彩色图像的直方图规定化

import cv2
import numpy as np


img1 = cv2.imread('D:\\DIP_Photo\\Belle.jpg')
# img2 = cv2.imread('D:\\DIP_Photo\\Bird.jpg')
# img2 = cv2.imread('D:\\DIP_Photo\\Mask.bmp')
img2 = cv2.imread('D:\\DIP_Photo\\colorblocks.jpeg')


img_hsv1 = cv2.cvtColor(img1, cv2.COLOR_BGR2HSV)  # bgr转hsv
img_hsv2 = cv2.cvtColor(img2, cv2.COLOR_BGR2HSV)

color = ('h', 's', 'v')

for i, col in enumerate(color):
    # histr = cv2.calcHist([img_hsv1], [i], None, [256], [0, 256])
    hist1, bins = np.histogram(img_hsv1[:, :, i].ravel(), 256, [0, 256])
    hist2, bins = np.histogram(img_hsv2[:, :, i].ravel(), 256, [0, 256])
    cdf1 = hist1.cumsum()  # 灰度值0-255的累计值数组
    cdf2 = hist2.cumsum()
    cdf1_hist = hist1.cumsum() / cdf1.max()  # 灰度值的累计值的比率
    cdf2_hist = hist2.cumsum() / cdf2.max()

    diff_cdf = [[0 for j in range(256)] for k in range(256)]  # diff_cdf 里是每2个灰度值比率间的差值
    for j in range(256):
        for k in range(256):
            diff_cdf[j][k] = abs(cdf1_hist[j] - cdf2_hist[k])

    lut = [0 for j in range(256)]  # 映射表
    for j in range(256):
        min = diff_cdf[j][0]
        index = 0
        for k in range(256):  # 直方图规定化的映射原理
            if min > diff_cdf[j][k]:
                min = diff_cdf[j][k]
                index = k
        lut[j] = ([j, index])

    h = int(img_hsv1.shape[0])
    w = int(img_hsv1.shape[1])
    for j in range(h):  # 对原图像进行灰度值的映射
        for k in range(w):
            img_hsv1[j, k, i] = lut[img_hsv1[j, k, i]][1]

hsv_img1 = cv2.cvtColor(img_hsv1, cv2.COLOR_HSV2BGR)  # hsv转bgr

cv2.imshow('firstpic', img1)
cv2.imshow('targetpic', img2)
cv2.imshow('defpic', hsv_img1)

cv2.waitKey(0)
cv2.destroyAllWindows()

四、同态滤波

import cv2
import numpy as np


# 同态滤波
def homomorphic_filter(src, d0=1, r1=2, rh=2, c=4, h=2.0, l=0.5):
    # 图像灰度化处理
    gray = src.copy()
    if len(src.shape) > 2:  # 维度>2
        gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)

    gray = np.float64(gray)  # 图像格式处理
    rows, cols = gray.shape  # 设置数据维度n

    gray_fft = np.fft.fft2(gray)  # 傅里叶变换
    gray_fftshift = np.fft.fftshift(gray_fft)  # 将零频点移到频谱的中间，就是中间化处理

    # 生成一个和gray_fftshift一样的全零数据结构
    # dst_fftshift = np.zeros_like(gray_fftshift)

    # arange函数用于创建等差数组，分解f(x,y)=i(x,y)r(x,y)
    M, N = np.meshgrid(np.arange(-cols // 2, cols // 2), np.arange(-rows // 2, rows // 2))  # 注意，//就是除法

    # 使用频率增强函数处理原函数（也就是处理原图像dst_fftshift）
    D = np.sqrt(M ** 2 + N ** 2)  # **2是平方
    Z = (rh - r1) * (1 - np.exp(-c * (D ** 2 / d0 ** 2))) + r1
    dst_fftshift = Z * gray_fftshift
    dst_fftshift = (h - l) * dst_fftshift + l

    # 傅里叶反变换（之前是正变换，现在该反变换变回去了）
    dst_ifftshift = np.fft.ifftshift(dst_fftshift)
    dst_ifft = np.fft.ifft2(dst_ifftshift)

    # 选取元素的实部
    dst = np.real(dst_ifft)

    # dst中，比0小的都会变成0，比0大的都变成255
    # uint8是专门用于存储各种图像的（包括RGB，灰度图像等），范围是从0–255
    dst = np.uint8(np.clip(dst, 0, 255))
    return dst


img = cv2.imread("D:\\DIP_Photo\\Diningroom.jpg", 0)
img_new = homomorphic_filter(img)  # 将图片执行同态滤波器

result = np.hstack((img, img_new))  # 输入和输出合并在一起输出
cv2.imshow('Dining room', result)

cv2.waitKey()
cv2.destroyAllWindows()

五、各类滤波比较

import cv2 as cv
import numpy as np

img = cv.imread("D:\\DIP_Photo\\JinChen02.jpg", 1)
kernel = np.ones([5, 5], np.float32)/25  # 用户自定义模糊,除以25是防止数值溢出

AvgBlurImg = cv.blur(img, (1, 15))  # 均值模糊（去随机噪声有很好的去噪效果,（1, 15）是垂直方向模糊）
MedBlurImg = cv.medianBlur(img, 5)  # 中值模糊（对椒盐噪声有很好的去燥效果，模板大小为5*5）
CustomBlurImg = cv.filter2D(img, -1, kernel)  # 自定义滤波器，模板大小为5*5）
BiBlurImg = cv.bilateralFilter(img, 30, 50, 5)  # 双边滤波器
ShiftImg = cv.pyrMeanShiftFiltering(img, 10, 50)  # 均值迁移滤波器,注意只能是彩色图像

result01 = np.hstack((img, AvgBlurImg, MedBlurImg))
result02 = np.hstack((CustomBlurImg, BiBlurImg, ShiftImg))

cv.imshow('BlurImg01', result01)
cv.imshow('BlurImg02', result02)

cv.waitKey(0)
cv.destroyAllWindows()

六、平滑滤波器

1. 盐噪声和中值滤波

import cv2
import numpy as np


def salt(SlatImg, n):
    for k in range(n):
        i = int(np.random.random() * SlatImg.shape[1])
        j = int(np.random.random() * SlatImg.shape[0])
        if SlatImg.ndim == 2:
            SlatImg[j, i] = 255
        elif SlatImg.ndim == 3:
            SlatImg[j, i, 0] = 255
            SlatImg[j, i, 1] = 255
            SlatImg[j, i, 2] = 255
    return SlatImg


img = cv2.imread("D:\\DIP_Photo\\Lenna_RGB.tif", 0)
cv2.imshow("Img", img)

result = salt(img, 500)
median = cv2.medianBlur(result, 3)

cv2.imshow("Salt", result)
cv2.imshow("Median", median)

cv2.waitKey(0)
cv2.destroyAllWindows()

2. 高斯噪声与高斯模糊

import cv2 as cv
import numpy as np


def clamp(pv):
    if pv > 255:
        return 255
    if pv < 0:
        return 0
    else:
        return pv


def gaussian_noise(image):  # 加高斯噪声
    h, w, c = image.shape
    for row in range(h):
        for col in range(w):
            s = np.random.normal(0, 20, 3)
            b = image[row, col, 0]  # blue
            g = image[row, col, 1]  # green
            r = image[row, col, 2]  # red
            image[row, col, 0] = clamp(b + s[0])
            image[row, col, 1] = clamp(g + s[1])
            image[row, col, 2] = clamp(r + s[2])
    cv.imshow("Noise image", image)


src = cv.imread("D:\\DIP_Photo\\Lenna_RGB.tif")
cv.imshow('input_image', src)

gaussian_noise(src)
dst = cv.GaussianBlur(src, (5, 5), 0)  # 高斯模糊

cv.imshow("Gaussian_Blur", dst)

cv.waitKey(0)
cv.destroyAllWindows()

七、锐化滤波器

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt

img = cv.imread("D:\\DIP_Photo\\Lenna_RGB.tif")  # 读取图像
lenna_img = cv.cvtColor(img, cv.COLOR_BGR2RGB)

grayImage = cv.cvtColor(lenna_img, cv.COLOR_BGR2GRAY)  # 灰度化处理图像
gaussianBlur = cv.GaussianBlur(grayImage, (3, 3), 0)  # 高斯滤波
ret, binary = cv.threshold(gaussianBlur, 127, 255, cv.THRESH_BINARY)  # 阈值处理

# Roberts算子
kernelx = np.array([[-1, 0], [0, 1]], dtype=int)
kernely = np.array([[0, -1], [1, 0]], dtype=int)
x = cv.filter2D(binary, cv.CV_16S, kernelx)
y = cv.filter2D(binary, cv.CV_16S, kernely)
absX = cv.convertScaleAbs(x)
absY = cv.convertScaleAbs(y)
Roberts = cv.addWeighted(absX, 0.5, absY, 0.5, 0)

# Prewitt算子
kernelx = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
kernely = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
x = cv.filter2D(binary, cv.CV_16S, kernelx)
y = cv.filter2D(binary, cv.CV_16S, kernely)
absX = cv.convertScaleAbs(x)
absY = cv.convertScaleAbs(y)
Prewitt = cv.addWeighted(absX, 0.5, absY, 0.5, 0)

# Sobel算子
x = cv.Sobel(binary, cv.CV_16S, 1, 0)
y = cv.Sobel(binary, cv.CV_16S, 0, 1)
absX = cv.convertScaleAbs(x)
absY = cv.convertScaleAbs(y)
Sobel = cv.addWeighted(absX, 0.5, absY, 0.5, 0)

# 拉普拉斯算法
dst = cv.Laplacian(binary, cv.CV_16S, ksize=3)
Laplacian = cv.convertScaleAbs(dst)

# 效果图
titles = ['Source Image', 'Binary Image', 'Roberts Image',
          'Prewitt Image', 'Sobel Image', 'Laplacian Image']
images = [lenna_img, binary, Roberts, Prewitt, Sobel, Laplacian]

for i in np.arange(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

八、一个综合例子

import numpy as np
import cv2
import matplotlib.pyplot as plt


def normalize(img):
    return (img - img.min()) / (img.max() - img.min())


def kernel_seperate(kernel):
    """
    get kernel seperate vector if separable
    param: kernel: input kerel of the filter
    return two separable vector c, r
    """
    rank = np.linalg.matrix_rank(kernel)

    if rank == 1:

        # -----------------Numpy------------------
        c = np.min(kernel, axis=1)
        r = np.min(kernel, axis=0)

    # ------------------Loop------------------
    #     r = w[0, :]
    #     c = []
    #     for i in range(w.shape[0]):
    #         temp = int((w[i, :] / r)[0])
    #         c.append(temp)

    else:
        print('Not separable')

    #     c = np.array(c).reshape(-1, 1)
    #     r = np.array(r).reshape(-1, 1)

    return c, r


def separate_kernel_conv2D(img, kernel):
    """
    separable kernel convolution 2D, if 2D kernel not separable, then canot use this fuction. in the fuction. first need to
    seperate the 2D kernel into two 1D vectors.
    param: img: input image you want to implement 2d convolution with 2d kernel
    param: kernel: 2D separable kernel, such as Gaussian kernel, Box Kernel
    return image after 2D convolution
    """
    m, n = kernel.shape[:2]
    pad_h = int((m - 1) / 2)
    pad_w = int((n - 1) / 2)
    w_1, w_2 = kernel_seperate(kernel)
    convovle = np.vectorize(np.convolve, signature='(x),(y)->(k)')
    r_2 = convovle(img, w_2)
    r_r = np.rot90(r_2)
    r_1 = convovle(r_r, w_1)
    r_1 = r_1.T
    r_1 = r_1[:, ::-1]
    img_dst = img.copy().astype(np.float)
    img_dst = r_1[pad_h:-pad_h, pad_w:-pad_w]
    return img_dst


def gamma_transform(img, c, gamma):
    """
    gamma transform 2d grayscale image, convert uint image to float
    param: input img: input grayscale image
    param: input c: scale of the transform
    param: input gamma: gamma value of the transoform
    """
    img = img.astype(float)  # 先要把图像转换成为float，不然结果点不太相同
    epsilon = 1e-5  # 非常小的值以防出现除0的情况
    img_dst = np.zeros(img.shape[:2], dtype=np.float)
    img_dst = c * np.power(img + epsilon, gamma)
    img_dst = np.uint8(normalize(img_dst) * 255)

    return img_dst


# 1 拉普拉斯突出细节
# 2 平滑后的梯度图像来掩蔽拉普拉斯图像
# 3 灰度变换增大灰度级的动态范围
# 图1
img_ori = cv2.imread("D:\\DIP_Photo\\Fig0343(a)(skeleton_orig).tif", 0)

# 图2，拉普拉斯变换
kernel_laplacian_d = np.array([
                    [-1, -1, -1],
                    [-1, 8, -1],
                    [-1, -1, -1]
                                ],)
img_laplacian = cv2.filter2D(img_ori, ddepth=-1, kernel=kernel_laplacian_d)
img_laplacian = np.uint8(normalize(img_laplacian) * 255)

# 图3，原图+拉普拉斯
img_ori_laplacian = img_ori + img_laplacian
img_ori_laplacian = normalize(img_ori_laplacian) * 255

# 图4，原图Sobel变换
sobel_x = np.zeros([3, 3], np.int)
sobel_x[0, :] = np.array([-1, -2, -1])
sobel_x[2, :] = np.array([1, 2, 1])
sobel_y = np.zeros([3, 3], np.int)
sobel_y[:, 0] = np.array([-1, -2, -1])
sobel_y[:, 2] = np.array([1, 2, 1])

gx = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_x)
gy = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_y)

img_sobel = abs(gx) + abs(gy)
img_sobel = np.uint8(normalize(img_sobel) * 255)

#  图5， 使用5x5的盒式滤波器平滑Sobel
kernel_box = np.ones([5, 5])
kernel_box = kernel_box / kernel_box.sum()
sobel_box = separate_kernel_conv2D(img_sobel, kernel=kernel_box)
sobel_box = normalize(sobel_box)
# sobel_box = np.uint8(normalize(sobel_box) * 255)

# 图6，图2与图5相乘的模板图像
mask = img_laplacian * sobel_box
img_mask = np.uint8(normalize(mask) * 255)

# 图7，原图与图6相加
img_passi = img_ori + img_mask * 0.3
img_passi = np.uint(normalize(img_passi) * 255)

# 图8 对图7做幂律变换
img_gamma = gamma_transform(img_passi, 1, gamma=0.5)

plt.figure(figsize=(12, 9))
plt.subplot(2, 4, 1), plt.imshow(img_ori,   'gray', vmax=255), plt.title("Original A"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 2), plt.imshow(img_laplacian, 'gray', vmax=255), plt.title("Laplacian B"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 3), plt.imshow(img_ori_laplacian, 'gray', vmax=255), plt.title("Original + Laplacian C"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 4), plt.imshow(img_sobel, 'gray', vmax=255), plt.title("Sobel D"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 5), plt.imshow(sobel_box, 'gray', vmax=1), plt.title("Sobel Box filter E"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 6), plt.imshow(img_mask, 'gray', vmax=255), plt.title("Sobel mask F"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 7), plt.imshow(img_passi, 'gray', vmax=255), plt.title("Passivation G"), plt.xticks([]), plt.yticks([])
plt.subplot(2, 4, 8), plt.imshow(img_gamma, 'gray', vmax=255), plt.title("Gamma Transform H"), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

---

代码示例仅供学生学习参考，如有引用未标注请提醒。