【Python CUDA版】河北工业大学计算机图像处理实验二：图像直方图及灰度变换

一、实验目的与要求

1．掌握图像灰度直方图的概念及其计算方法，编写灰度直方图统计程序。

2．通过对图像直方图的分析，学习应用直方图法解决诸如图像二值化等具体问题。

3．熟悉直方图均衡化的计算过程及其应用。

4．掌握图像灰度变换技术，通过调整图像的对比度和亮度等参数，改善视觉效果。

二、实验相关知识

直方图是图像最基本的统计特征，是图像亮度分布的概率密度函数，反映了图像灰度值的分布情况。直方图是多种空间域处理技术的基础。直方图操作能有效地用于图像增强，如通过直方图均衡化处理，可使图像在整个灰度级范围内的分布均匀化，即在每个灰度级上都具有相同的像素点数，从而获得较好的视觉效果。另外，直方图固有的信息也可用在图像分割等其它图像处理的应用中。

灰度级变换技术可用g(x,y)=T[f(x,y)]的形式表示，其中f(x,y)为输入图像，g(x,y)为输出图像，T是对图像f进行某种处理的操作。由于(x,y)处的g值仅由f在该点处的亮度决定，T也称为一个亮度或灰度级变换函数，它与位置坐标(x,y)无关，所以通常写成如下的简化形式s=T(r)，其中r和s分别表示图像f和g在相应点(x,y)的亮度。利用灰度变换可调整一幅图像的明暗、对比度等。

三、实验内容

1、编写一个图像灰度直方图统计函数my_imhist，选择一幅图像利用my_imhist显示其直方图，将结果与MATLAB图像处理工具箱中提供的灰度直方图函数imhist的处理结果进行比较，并在同一窗口中显示出来。

2、利用以上编写的函数my_imhist或imhist，估算图像iris.tif中瞳孔的半径（以像素为单位）。

3、按照教材68页上的公式（4.1.6）

        

编程实现图像的分段线性灰度变换.

4、编写一个灰度图像的直方图均衡化函数（不可使用库函数）。（可使用的灰度等级数量不变即可），并对下图例题进行测试，给出测试结果。

四、实验代码及结果 

1、灰度直方图统计 

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


def my_imhist(image):
    count_pixel = np.zeros(256)
    size = image.shape
    for i in range(0, size[0]):
        for j in range(0, size[1]):
            count_pixel[image[i][j]] = count_pixel[image[i][j]] + 1
    return count_pixel


if __name__ == '__main__':
    plt.switch_backend('agg')
    plt.figure(dpi=500)
    # 读入原始图像
    img = cv2.imread('22_Cage2_4k.jpg')
    plt.subplot(221)
    # 通道转换
    plt.imshow(img[:, :, ::-1])
    plt.axis(False)
    plt.title("origin")
    # 读入灰度图像
    img = cv2.imread('22_Cage2_4k.jpg', 0)
    plt.subplot(222)
    plt.imshow(img, cmap="gray")
    plt.axis(False)
    plt.title("gray")
    # my_imhist 结果
    plt.subplot(223)
    plt.bar(range(256), my_imhist(img), width=1)
    plt.title("my_imhist")
    # 使用OpenCV库函数 结果
    plt.subplot(224)
    greyScale_map = cv2.calcHist([img], [0], None, [256], [0, 256])
    greyScale_map = np.array(greyScale_map).reshape(256)
    plt.bar(range(256), greyScale_map, width=1)
    plt.title("use opencv")

    plt.tight_layout()
    plt.show()

    plt.savefig("Exp2-1.png")

 

2、估算半径

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

def my_imhist(image):
    size = image.shape
    count_pixel = np.zeros(256)
    for i in range(0, size[0]):
        for j in range(0, size[1]):
            count_pixel[image[i][j]] = count_pixel[image[i][j]] + 1
    print(count_pixel)
    return count_pixel


if __name__ == "__main__":
    img = cv.imread("iris.png", 0)  # 以灰度图像读入
    # 找出最小像素值
    calHist = np.array(my_imhist(img))
    print(calHist)
    plt.switch_backend('agg')
    plt.figure(dpi=300)
    plt.bar(range(256), my_imhist(img), width=1)
    plt.savefig('Exp2-2.png')
    min_pixel = np.min(calHist)  # 找到列表中第一个非0 的数
    # 根据像素数量算出半径
    size = img.shape
    count = 0
    for i in range(size[0]):
        for j in range(size[1]):
            if img[i][j] - min_pixel < 30:
                count += 1
    print("瞳孔半径为：", np.sqrt(count / np.pi), "个像素值")

 

3、线性变换 

import math
import time
import matplotlib.pyplot as plt
import cv2 as cv
import argparse
import numpy as np
from numba import cuda

parser = argparse.ArgumentParser(description='')
parser.add_argument("--a", type=float, default=255 / 3, help="numerical value of a")
parser.add_argument("--b", type=float, default=255 / 3, help="numerical value of b")
parser.add_argument("--c", type=float, default=255 / 3, help="numerical value of c")
parser.add_argument("--d", type=float, default=255 / 3, help="numerical value of d")
args = parser.parse_args()


@cuda.jit()
def gpu_process(img):
    tx = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x
    ty = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y

    if 0 <= img[tx, ty] < a:
        img[tx, ty] = c / a * img[tx, ty]
    if a <= img[tx, ty] < b:
        img[tx, ty] = ((d - c) / (b - a)) * (img[tx, ty] - a) + c
    else:
        img[tx, ty] = ((m_g - d) / (m_f - b)) * (img[tx, ty] - b) + d


if __name__ == "__main__":
    img = cv.imread("22_Cage2_4k.jpg", 0)
    a, b, c, d = args.a, args.b, args.c, args.d
    m_f = np.max(img.ravel())
    m_g = m_f
    print(m_f)
    size = img.shape
    print(size)

    dImg = cuda.to_device(img)
    threadsperblock = (32, 32)
    blockspergrid_x = int(math.ceil(size[0] / threadsperblock[0]))
    blockspergrid_y = int(math.ceil(size[1] / threadsperblock[1]))
    blockspergrid = (blockspergrid_x, blockspergrid_y)
    cuda.synchronize()
    gpu_process[blockspergrid, threadsperblock](dImg)
    img_new = dImg.copy_to_host()
    plt.switch_backend('agg')
    plt.figure(dpi=500)
    plt.subplot(121)
    plt.title("origin")
    plt.imshow(img, cmap="gray")
    plt.subplot(122)
    plt.title("changed")
    plt.imshow(img_new, cmap="gray")
    plt.savefig('Exp2-3.png')
    plt.clf()
    plt.xlabel("f(i,j)")
    plt.ylabel("g(i,j)")
    plt.plot([0, m_f], [0, m_g], label="f=g")
    plt.plot([0, a, b, m_f], [0, c, d, m_g], label="f -> g")
    plt.legend()
    plt.savefig('Exp2-3_fig.png')

 

 4、灰度直方图均衡化

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

parser = argparse.ArgumentParser(description='')
parser.add_argument("--file", type=str)
parser.add_argument("--scale", type=int)
args = parser.parse_args()

if __name__ == "__main__":
    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 显示中文
    data = np.loadtxt(args.file).astype(int)
    shape = data.shape
    print(data)
    image = np.zeros((data.shape[0] * 10, data.shape[1] * 10))
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            for m in range(i * 10, i * 10 + 10):
                for n in range(j * 10, j * 10 + 10):
                    image[m, n] = data[i, j] * 255 / 7
    plt.subplot(121)
    plt.imshow(image, cmap="gray")
    plt.axis(False)
    plt.title("原始图像")

    data_f = data.flatten().astype(int)
    n_k = np.zeros(args.scale)
    p_r = []
    s_k = []
    new = np.zeros(shape).astype(int)
    for i in range(data_f.size):
        n_k[int(data_f[i])] += 1
    for i in range(args.scale):
        p_r.append(round(n_k[i] / (shape[0] * shape[1]), 4))
    for i in range(args.scale):
        s_k.append(round(sum(p_r[:(i + 1)]), 4))
    print("s_k: ", s_k)
    new_array = np.around(s_k * np.array([7])).astype(int)
    print("new scale: ", new_array)
    for i in range(shape[0]):
        for j in range(shape[1]):
            new[i][j] = new_array[data[i][j]]
    print("new image: \n", new)

    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            for m in range(i * 10, i * 10 + 10):
                for n in range(j * 10, j * 10 + 10):
                    image[m, n] = new[i, j] * 255 / 7
    plt.subplot(122)
    plt.imshow(image, cmap="gray")
    plt.axis(False)
    plt.title("新图像")
    plt.show()