彩色图像直方图均衡——HSI模型

最新推荐文章于 2024-05-06 16:49:26 发布

csdn_newbie_uniq

最新推荐文章于 2024-05-06 16:49:26 发布

阅读量364

点赞数

文章标签：算法

本文链接：https://blog.csdn.net/csdn_newbie_uniq/article/details/130772903

版权

前言

参考书籍：《数字图像处理（第三版）》冈萨雷斯
参考博客：参考博客
限定范围：处理 24 位色深 BMP 格式文件

RGB转HSI公式

RGB 值先除以 $255$ ，转到 $[0, 1]$ 区间。
$\begin{aligned} \theta&=\arccos\frac{\displaystyle\frac 12[(R-G)+(R-B)]}{[(R-G)^2+(R-B)(G-B)]^{1/2}}\\ H&=\begin{cases} \theta,&B\leq G\\ 2\pi-\theta,&B>G \end{cases}\\ S&=1-\frac{3}{R+G+B}[\min(R,G,B)]\\ I&=\frac 13(R+G+B) \end{aligned}$

HSI转RGB公式

当 $0\leq H< \displaystyle\frac{2\pi}3$ 时：
$\begin{aligned} B&=I(1-S)\\ R&=I[1+\frac{S\cos H}{cos(\displaystyle\frac{\pi}3-H)}]\\ G&=3I-(R+B) \end{aligned}$
当 $\displaystyle\frac{2\pi}3 \leq H< \displaystyle\frac{4\pi}3$ 时：
$\begin{aligned} H&=H-\displaystyle\frac{2\pi}3\\ R&=I(1-S)\\ G&=I[1+\frac{S\cos H}{cos(\displaystyle\frac{\pi}3-H)}]\\ B&=3I-(R+G) \end{aligned}$
当 $\displaystyle\frac{4\pi}3 \leq H< 2\pi$ 时：
$\begin{aligned} H&=H-\displaystyle\frac{4\pi}3\\ G&=I(1-S)\\ B&=I[1+\frac{S\cos H}{cos(\displaystyle\frac{\pi}3-H)}]\\ R&=3I-(G+B) \end{aligned}$
转换后的 RGB 值在 $[0, 1]$ 区间，乘以 $255$ 转到原区间。

代码

#include<iostream>
#include<windows.h>
#include<cmath>
#include<vector>
#define pi acos(-1)
using namespace std;

vector<vector<double>> model;

int main()
{
    //fp为原图像文件指针，fp1为处理后图像文件指针，该文件可不存在
    FILE* fp = fopen("cry.bmp", "rb");
    FILE* fp1 = fopen("test.bmp", "wb");
    if (!fp)
    {
        cout << "原图像文件不存在，请保证原图像与代码在同一目录下" << endl;
        return -1;
    }
    BITMAPFILEHEADER* fileheader = new BITMAPFILEHEADER;
    BITMAPINFOHEADER* infoheader = new BITMAPINFOHEADER;
    fread(fileheader, 14, 1, fp);
    fread(infoheader, 40, 1, fp);
    int bits = infoheader->biBitCount;
    if (bits != 24)
    {
        cout << "该文件不符合要求" << endl;
        return -1;
    }
    int width = infoheader->biWidth, height = infoheader->biHeight;
    int lineByte = (width * 3 + 3) / 4 * 4;
    fileheader->bfSize = 14 + 40 + lineByte * height;
    infoheader->biSizeImage = lineByte * height;
    fwrite(fileheader, 14, 1, fp1);
    fwrite(infoheader, 40, 1, fp1);
    model.resize(height, vector<double>(3 * width));
    unsigned char* origin = new unsigned char[lineByte * height];
    memset(origin, 0, lineByte * height);
    double counter[256] = { 0 };
    double R, G, B;
    double H, S, I;
    int x = 0, y = 0;
    double eps = 1e-6;
    double Smin = 10, Smax = -1, Imin = 10, Imax = -1;
    double Rmin = 10, Rmax = -1, Gmin = 10, Gmax = -1, Bmin = 10, Bmax = -1;
    double low, top, bottom, theta;
    fread(origin, 1, lineByte * height, fp);
    fclose(fp);

    //rgb转hsi
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            model[i][3 * j] = *(origin + i * lineByte + 3 * j + 2);
            model[i][3 * j + 1] = *(origin + i * lineByte + 3 * j + 1);
            model[i][3 * j + 2] = *(origin + i * lineByte + 3 * j);
            R = model[i][3 * j] / 255;
            G = model[i][3 * j + 1] / 255;
            B = model[i][3 * j + 2] / 255;
            model[i][3 * j + 2] = (R + G + B) / 3;
            if (R == G && G == B && B == R) model[i][3 * j + 1] = 0;
            else
            {
                low = min(min(R, G), B);
                model[i][3 * j + 1] = 1 - 3 * low / (R + G + B);
            }
            if (model[i][3 * j + 1] == 0) model[i][3 * j] = 0;
            else {
                top = 0.5 * (2 * R - G - B);
                bottom = sqrt((R - G) * (R - G) + (R - B) * (G - B));
                theta = acos(top / (bottom +eps));
                if (B <= G) model[i][3 * j] = theta;
                else model[i][3 * j] = 2 * pi - theta;
            }
            Smin = min(Smin, model[i][3 * j + 1]);
            Imin = min(Imin, model[i][3 * j + 2]);
        }

    //对s、i进行归一化
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            model[i][3 * j + 1] -= Smin;
            model[i][3 * j + 2] -= Imin;
            Smax = max(Smax, model[i][3 * j + 1]);
            Imax = max(Imax, model[i][3 * j + 2]);
        }
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            if (Smax) model[i][3 * j + 1] /= Smax;
            if (Imax) model[i][3 * j + 2] /= Imax;
            counter[(int)(model[i][3 * j + 2] * 255 + 0.5)]++;
        }

    for (int i = 1;i < 256;i++)
        counter[i] += counter[i - 1];

    //hsi转rgb
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            model[i][3 * j + 2] = counter[(int)(model[i][3 * j + 2] * 255 + 0.5)] / width / height;
            H = model[i][3 * j];
            S = model[i][3 * j + 1];
            I = model[i][3 * j + 2];

            if (H < 2 * pi / 3)
            {
                B = I * (1 - S);
                R = I * (1 + S * cos(H) / cos(pi / 3 - H));
                G = 3 * I - R - B;
            }
            else if (H < 4 * pi / 3)
            {
                H -= 2 * pi / 3;
                R = I * (1 - S);
                G = I * (1 + S * cos(H) / cos(pi / 3 - H));
                B = 3 * I - R - G;
            }
            else
            {
                H -= 4 * pi / 3;
                G = I * (1 - S);
                B = I * (1 + S * cos(H) / cos(pi / 3 - H));
                R = 3 * I - G - B;
            }
            model[i][3 * j] = R;
            model[i][3 * j + 1] = G;
            model[i][3 * j + 2] = B;
            Rmin = min(Rmin, R);
            Gmin = min(Gmin, G);
            Bmin = min(Bmin, B);
        }

    //对r、g、b进行归一化
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            model[i][3 * j] -= Rmin;
            model[i][3 * j + 1] -= Gmin;
            model[i][3 * j + 2] -= Bmin;
            Rmax = max(Rmax, model[i][3 * j]);
            Gmax = max(Gmax, model[i][3 * j + 1]);
            Bmax = max(Bmax, model[i][3 * j + 2]);
        }
    for (int i = 0;i < height;i++)
        for (int j = 0;j < width;j++)
        {
            if (Rmax) model[i][3 * j] /= Rmax;
            if (Gmax) model[i][3 * j + 1] /= Gmax;
            if (Bmax) model[i][3 * j + 2] /= Bmax;
            *(origin + i * lineByte + 3 * j) = (unsigned char)(model[i][3 * j + 2] * 255 + 0.5);
            *(origin + i * lineByte + 3 * j + 1) = (unsigned char)(model[i][3 * j + 1] * 255 + 0.5);
            *(origin + i * lineByte + 3 * j + 2) = (unsigned char)(model[i][3 * j] * 255 + 0.5);
        }
    fwrite(origin, lineByte * height, 1, fp1);
    fclose(fp1);
    return 0;
}

代码说明

model 数组中当存的是 RGB 值时按照 R、G、B 的顺序存储，存的是 HSI 值时按照 H、S、I 的顺序存储，注意 BMP 文件的 RGB 值是按照 B、G、R 的顺序存储，最后存到 origin 数组中需调换顺序。
eps 是一个极小的量，不会影响最终结果，但是它是必要的，它能防止在计算 $\theta$ 时出现 $0/0$ 的状况，而且通常 $0/0$ 不会报错，但输出图像时会出现很多无法察觉的噪声。
这是原图像：
在这里插入图片描述
这是不加 eps 时的图像（不同编译器不同，有的编译器输出结果变化特别大，此为 g++，MSVC比较逆天）：

放大到阿尼亚的额头处：

加了 eps 后恢复正常：

注意归一化过程中做除法前判断分母是否为 $0$ ，为 $0$ 的情况有原图所有像素 RGB 均相等（即实际上为灰度图时）、所有像素亮度值相等和其他的一些可能原因，总之需注意。