C#绘制CIE1931色度图

CIE 1931 色度图
CIE 1931 色度图是一个理想的图形，如图 2.所示。理论上的马蹄形曲线内区域包括了一切物理上能实
现的颜色。在此二维色度图中，X 轴色度坐标相当于红原色的比例，Y 轴色度坐标相当于绿原色的比例，并
且有 X＋Y＋Z＝1，由 X、Y 值可得出 Z 值。图中没有 Z 轴色度坐标，但 Z 相当于蓝原色的比例。

描绘马蹄形曲线
描绘马蹄形曲线的方法：选取 46 个光谱轨迹色度坐标。
由于马蹄形曲线并非随波长变化而均匀变化，所以选取坐标时并未按照等波长间隔选取，而是根据观
察马蹄形曲线，得出马蹄形曲线上长度大概相同的点，对照 CIE 1931 光谱轨迹色度坐标表选取。
将波长 380nm 点作为纯蓝色点。实际上该点应该为蓝色与红色的混合色：紫色点，但由于紫色是混合
色，绘制时不好控制计算机上 R、B 的比例，所以绘制马蹄形曲线的时候将该点作为纯蓝色点。设色度图中
纯蓝色点为 B 点、纯绿色点（520nm）为 G 点、纯红色点（780nm）为 R 点；则 B 点在计算机上的颜色设为
（0，0，255）、G 点为（0，255，0）、R 点为（255，0，0）。
在波长 380nm 到 520nm 之间，选取 25 个点，然后，用 B 样条曲线算法用圆滑连续曲线连接这 25 个点。
连接的时候，颜色变化规律为：B 点到 498nm（由光谱轨迹色度坐标表得到此处为 B、G 比例近似相同点）
颜色从（0，0，255）变化到（0，255，255）；498nm 到 G 点颜色从（0，255，255）变化到（0，255，0）。
同理，在波长 520nm 到 780nm 之间，选取 20 个点，然后，用 B 样条曲线算法用圆滑连续曲线连接这 20
个点。连接的时候，颜色变化规律为： G 点到 575nm（由光谱轨迹色度坐标表得到此处为 G、R 比例近似相
同点）颜色从（0，255，0）变化到（255，255，0）；575nm 到 R 点颜色从（255，255，0）变化到（255，0，
0）。
由于连接波长 780 到 380nm 两点的是直线，所以在这两点间设 511 个象素点，然后以画点的方式画从
红到蓝的直线。其变化规律相似，颜色从（255，0，0）变化到（255，0，255）；再从（255，0，255）变
化到（0，0，255）。
为了表现尽可能多的颜色， 在颜色渐变过程中采用了一种“准”线性变换的办法。即颜色在 B 点到 498nm、
498nm 到 G 点、G 点到 575nm、575nm 到 R 点以及 R 点到 B 点变化时，均为均匀变化。
这样，便绘出了马蹄形曲线。有了曲线后，以等能白点 C 点(255，255，255)（色度坐标为 X＝0．3333，
Y＝0．3333）为中心，向色度图边缘上的各点作射线，射线上的点的颜色则在白色与这条射线和马蹄形曲
线边界交点的颜色之间连续变化。此时的变化仍采用线性均匀变化。为了去除中间不连续部分，射线扫描
先由绿→蓝、绿→红→蓝，再由红→蓝、红→绿做了两次扫描，基本做到了连续。

private void XYZ2RGB(float x, float y, float z, ref int r, ref int g, ref int b)
{
    double dr, dg, db;
    dr = 0.4185 * x - 0.1587 * y - 0.0828 * z;
    dg = -0.0912 * x + 0.2524 * y + 0.0157 * z;
    db = 0.0009 * x - 0.0025 * y + 0.1786 * z;
    double max = 0;
    max = dr > dg ? dr : dg;
    max = max > db ? max : db;
    dr = dr / max * 255;
    dg = dg / max * 255;
    db = db / max * 255;
    dr = dr > 0 ? dr : 0;
        dg = dg > 0 ? dg : 0;
        db = db > 0 ? db : 0;
        if (dr > 255)
        {
            dr = 255;
        }
        if (dg > 255)
        {
            dg = 255;
        }
        if (db > 255)
        {
            db = 255;
        }
        r = (int)(dr + 0.5);
        g = (int)(dg + 0.5);
        b = (int)(db + 0.5);
    }

private void DrawCIE()
{
    List<float> mySmallx = new List<float>();
    List<float> mySmally = new List<float>();
    List<float> mySmallz = new List<float>();
    using (FileStream fs = new FileStream("轮廓坐标.txt", FileMode.Open, FileAccess.Read, FileShare.ReadWrite))
    {
        lock (fs)
        {
            StreamReader sr = new StreamReader(fs, Encoding.UTF8);
            while (!sr.EndOfStream)
            {
                string str = sr.ReadLine();
                string[] temp = str.Split(new char[] { }, StringSplitOptions.RemoveEmptyEntries);
                if (temp.Length < 3) continue;
                mySmallx.Add(float.Parse(temp[1]));
                mySmally.Add(float.Parse(temp[2]));
            }
        }
    }

    //1:选择3个顶点
    PointF p1 = new PointF(); // 代表最下面的点
    PointF p2 = new PointF(); // 代表最上面的点
    PointF p3 = new PointF(); // 代表最右边的点
    //代表最下面的点
    p1.X = 0.173101f;
    p1.Y = 0.004774f;

    //代表最上面的点
    p2.X = 0.082053f;
    p2.Y = 0.83409f;

    //代表最右边的点
    p3.X = 0.73469f;
    p3.Y = 0.26531f;

    float k3 = 0, b3 = 0;
    k3 = (p3.Y - p1.Y) / (p3.X - p1.X);
    b3 = p3.Y - p3.X * k3;

    //左边数组
    List<PointF> leftArray = new List<PointF>();
    List<PointF> rightArray = new List<PointF>();
    PointF tempPonit = new PointF();
    for (int i = 0; i < (int)mySmallx.Count(); i++)
    {
        tempPonit.X = mySmallx[i];
        tempPonit.Y = mySmally[i];
        if (i < 161)
        {
            leftArray.Add(tempPonit);
            if (i == 160)
            {
                rightArray.Add(tempPonit);
            }
        }
        else if (i >= 161 && i < 339)
        {
            rightArray.Add(tempPonit);
        }
    }

    List<PointF> downArray = new List<PointF>();
    for (float y = 0.001f; y < 0.265; y += 0.001f)
    {
        if (y > p1.Y && y < p3.Y)
        {
            tempPonit.X = (y - b3) / k3;
            tempPonit.Y = y;
            downArray.Add(tempPonit);
        }
    }
    Bitmap bitmap = new Bitmap(2000, 2000, System.Drawing.Imaging.PixelFormat.Format24bppRgb);
    Rectangle dimension = new Rectangle(0, 0, bitmap.Width, bitmap.Height);
    BitmapData picData = bitmap.LockBits(dimension, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
    
    float tempx = 0, tempy = 0;
    float miny = 0, maxy = 0;
    float minx = 0, maxx = 0;
    float middlex = 0;
    double minvalue1 = 999999;
    double minvalue2 = 999999;
    int maxindex = 0, minindex = 0;
    for (int i = 0; i < picData.Height; i++)
    {
        tempy = (float)i / (float)picData.Height;
        if (tempy < 0.004774 || tempy > 0.83409)
        {
            continue;
        }
        for (int j = 0; j < picData.Width; j++)
        {
            tempx = (float)j / (float)picData.Width;
            if (tempx < 0.003636 || tempx > 0.73469)
            {
                continue;
            }
            //1、判断左边的数组
            maxindex = 0;
            minindex = 0;
            minvalue2 = 999999;
            minvalue1 = 999999;
            for (int k = 0; k < (int)leftArray.Count(); k++)
            {
                double value1 = tempy - leftArray[k].Y;
                double value2 = leftArray[k].Y - tempy;
                if (value1 < minvalue1 && value1 > 0)
                {
                    minvalue1 = value1;
                    minindex = k;
                }
                if (value2 < minvalue2 && value2 > 0)
                {
                    minvalue2 = value2;
                    maxindex = k;
                }
            }
            miny = leftArray[minindex].Y;
            minx = leftArray[minindex].X;
            maxy = leftArray[maxindex].Y;
            maxx = leftArray[maxindex].X;
            middlex = minx - (tempy - miny) / (maxy - miny) * (minx - maxx);
            if (tempx < middlex || maxindex == minindex)
            {
                continue;
            }

            //2、判断右边的数组
            maxindex = 0;
            minindex = 0;
            minvalue2 = 999999;
            minvalue1 = 999999;
            for (int k = 0; k < (int)rightArray.Count(); k++)
            {
                double value1 = tempy - rightArray[k].Y;
                double value2 = rightArray[k].Y - tempy;

                if (value1 < minvalue1 && value1 > 0)
                {
                    minvalue1 = value1;
                    minindex = k;
                }
                if (value2 < minvalue2 && value2 > 0)
                {
                    minvalue2 = value2;
                    maxindex = k;
                }
            }
            miny = rightArray[minindex].Y;
            minx = rightArray[minindex].X;
            maxy = rightArray[maxindex].Y;
            maxx = rightArray[maxindex].X;
            middlex = minx - (tempy - miny) / (maxy - miny) * (minx - maxx);
            if (tempx > middlex || maxindex==minindex)
            {
                continue;
            }

            //3、判断下边的数组
            maxindex = 0;
            minindex = 0;
            minvalue2 = 999999;
            minvalue1 = 999999;
            for (int k = 0; k < (int)downArray.Count(); k++)
            {
                double value1 = tempy - downArray[k].Y;
                double value2 = downArray[k].Y - tempy;
                if (value1 < minvalue1 && value1 > 0)
                {
                    minvalue1 = value1;
                    minindex = k;
                }
                if (value2 < minvalue2 && value2 > 0)
                {
                    minvalue2 = value2;
                    maxindex = k;
                }
            }
            miny = downArray[minindex].Y;
            minx = downArray[minindex].X;
            maxy = downArray[maxindex].Y;
            maxx = downArray[maxindex].X;
            middlex = minx - (tempy - miny) / (maxy - miny) * (minx - maxx);
            if (tempx > middlex || maxindex == minindex)
            {
                continue;
            }
            float x = 0, y = 0, z = 0;
            x = tempx;
            y = tempy;
            z = 1 - x - y;
            int R = 0, G = 0, B = 0;
            XYZ2RGB(x, y, z, ref R, ref G, ref B);
            unsafe
            {
                byte* target = (byte*)picData.Scan0.ToPointer();
                target[(picData.Width - i) * picData.Stride + 3 * j] = (byte)B;
                target[(picData.Width - i) * picData.Stride + 3 * j + 1] = (byte)G;
                target[(picData.Width - i) * picData.Stride + 3 * j + 2] = (byte)R;
            }
        }
    }
    bitmap.UnlockBits(picData);
    bitmap = bitmap.Clone(new Rectangle(0, 200, 1600, 1800), bitmap.PixelFormat);
    bitmap.Save("Cie1931.bmp", ImageFormat.Bmp);
}

结果：

项目源码下载地址