图像处理: 可见光波长(ｗavelength)与RGB之间的转换

野性酷女 2022-01-17 09:15 399阅读 0赞

关于颜色转换，有时间看标准当然最好，不过我手头一时没有相关标准，所以参考了一些网上的资料，

380nm 到 760nm的可见光 对应的RGB大概是这个样子的，

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3Rhbm14MjE5_size_16_color_FFFFFF_t_70][]

C语言版本程序，

//指定波长转换成RGBA颜色
    std::vector<int> lambdaToColor(double lambda,double gamma = 0.8,double intensityMax = 255.0)
    {
        double r, g, b, alpha;
        if (lambda >= 380.0 && lambda < 440.0) {
            r = -1.0 * (lambda - 440.0) / (440.0 - 380.0);
            g = 0.0;
            b = 1.0;
        }else if (lambda >= 440.0 && lambda < 490.0) {
            r = 0.0;
            g = (lambda - 440.0) / (490.0 - 440.0);
            b = 1.0;
        }else if (lambda >= 490.0 && lambda < 510.0) {
            r = 0.0;
            g = 1.0;
            b = -1.0 * (lambda - 510.0) / (510.0 - 490.0);
        }else if (lambda >= 510.0 && lambda < 580.0) {
            r = (lambda - 510.0) / (580.0 - 510.0);
            g = 1.0;
            b = 0.0;
        }else if (lambda >= 580.0 && lambda < 645.0) {
            r = 1.0;
            g = -1.0 * (lambda - 645.0) / (645.0 - 580.0);
            b = 0.0;
        }else if (lambda >= 645.0 && lambda <= 780.0) {
            r = 1.0;
            g = 0.0;
            b = 0.0;
        }else {
            r = 0.0;
            g = 0.0;
            b = 0.0;
        }
    
        //在可见光谱的边缘处强度较低。
        if (lambda >= 380.0 && lambda < 420.0) {
            alpha = 0.30 + 0.70 * (lambda - 380.0) / (420.0 - 380.0);
        }else if (lambda >= 420.0 && lambda < 701.0) {
            alpha = 1.0;
        }else if (lambda >= 701.0 && lambda < 780.0) {
            alpha = 0.30 + 0.70 * (780.0 - lambda) / (780.0 - 700.0);
        }else {
            alpha = 0.0;
        }
    
        //1953年在引入NTSC电视时,计算具有荧光体的监视器的亮度公式如下
        int Y = static_cast<int>(0.212671*r + 0.715160*g + 0.072169*b);
    
        //伽马射线 gamma
        //照明强度 intensityMax
        int R = r == 0.0 ? 0 : static_cast<int>(std::round(intensityMax * std::pow(r * alpha, gamma)));
        int G = g == 0.0 ? 0 : static_cast<int>(std::round(intensityMax * std::pow(g * alpha, gamma)));
        int B = b == 0.0 ? 0 : static_cast<int>(std::round(intensityMax * std::pow(b * alpha, gamma)));
        int A = static_cast<int>(alpha);
    
        return std::vector<int>{R, G, B, A, Y};
    }

另一个转换是通过建立一个事先计算好的表来计算的，这个还没细看，不过感觉通过建立查找表，牺牲一点点空间，可以极大减少运算量，这也是高速程序的一般处理方式。

static class RgbCalculator {
    
        const int
             LEN_MIN = 380,
             LEN_MAX = 780,
             LEN_STEP = 5;
    
        static readonly double[]
            X = {
                    0.000160, 0.000662, 0.002362, 0.007242, 0.019110, 0.043400, 0.084736, 0.140638, 0.204492, 0.264737,
                    0.314679, 0.357719, 0.383734, 0.386726, 0.370702, 0.342957, 0.302273, 0.254085, 0.195618, 0.132349,
                    0.080507, 0.041072, 0.016172, 0.005132, 0.003816, 0.015444, 0.037465, 0.071358, 0.117749, 0.172953,
                    0.236491, 0.304213, 0.376772, 0.451584, 0.529826, 0.616053, 0.705224, 0.793832, 0.878655, 0.951162,
                    1.014160, 1.074300, 1.118520, 1.134300, 1.123990, 1.089100, 1.030480, 0.950740, 0.856297, 0.754930,
                    0.647467, 0.535110, 0.431567, 0.343690, 0.268329, 0.204300, 0.152568, 0.112210, 0.081261, 0.057930,
                    0.040851, 0.028623, 0.019941, 0.013842, 0.009577, 0.006605, 0.004553, 0.003145, 0.002175, 0.001506,
                    0.001045, 0.000727, 0.000508, 0.000356, 0.000251, 0.000178, 0.000126, 0.000090, 0.000065, 0.000046,
                    0.000033
                },
    
            Y = {
                    0.000017, 0.000072, 0.000253, 0.000769, 0.002004, 0.004509, 0.008756, 0.014456, 0.021391, 0.029497,
                    0.038676, 0.049602, 0.062077, 0.074704, 0.089456, 0.106256, 0.128201, 0.152761, 0.185190, 0.219940,
                    0.253589, 0.297665, 0.339133, 0.395379, 0.460777, 0.531360, 0.606741, 0.685660, 0.761757, 0.823330,
                    0.875211, 0.923810, 0.961988, 0.982200, 0.991761, 0.999110, 0.997340, 0.982380, 0.955552, 0.915175,
                    0.868934, 0.825623, 0.777405, 0.720353, 0.658341, 0.593878, 0.527963, 0.461834, 0.398057, 0.339554,
                    0.283493, 0.228254, 0.179828, 0.140211, 0.107633, 0.081187, 0.060281, 0.044096, 0.031800, 0.022602,
                    0.015905, 0.011130, 0.007749, 0.005375, 0.003718, 0.002565, 0.001768, 0.001222, 0.000846, 0.000586,
                    0.000407, 0.000284, 0.000199, 0.000140, 0.000098, 0.000070, 0.000050, 0.000036, 0.000025, 0.000018,
                    0.000013
                },
    
            Z = {
                    0.000705, 0.002928, 0.010482, 0.032344, 0.086011, 0.197120, 0.389366, 0.656760, 0.972542, 1.282500,
                    1.553480, 1.798500, 1.967280, 2.027300, 1.994800, 1.900700, 1.745370, 1.554900, 1.317560, 1.030200,
                    0.772125, 0.570060, 0.415254, 0.302356, 0.218502, 0.159249, 0.112044, 0.082248, 0.060709, 0.043050,
                    0.030451, 0.020584, 0.013676, 0.007918, 0.003988, 0.001091, 0.000000, 0.000000, 0.000000, 0.000000,
                    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                    0.000000
                };
    
        static readonly double[]
            MATRIX_SRGB_D65 = {
                 3.2404542, -1.5371385, -0.4985314,
                -0.9692660,  1.8760108,  0.0415560,
                 0.0556434, -0.2040259,  1.0572252
            };
    
        public static byte[] Calc(double len) {
            if(len < LEN_MIN || len > LEN_MAX)
                return new byte[3];
    
            len -= LEN_MIN;
            var index = (int)Math.Floor(len / LEN_STEP);
            var offset = len - LEN_STEP * index;
    
            var x = Interpolate(X, index, offset);
            var y = Interpolate(Y, index, offset);
            var z = Interpolate(Z, index, offset);
    
            var m = MATRIX_SRGB_D65;
    
            var r = m[0] * x + m[1] * y + m[2] * z;
            var g = m[3] * x + m[4] * y + m[5] * z;
            var b = m[6] * x + m[7] * y + m[8] * z;
    
            r = Clip(GammaCorrect_sRGB(r));
            g = Clip(GammaCorrect_sRGB(g));
            b = Clip(GammaCorrect_sRGB(b));
    
            return new[] { 
                (byte)(255 * r),
                (byte)(255 * g),
                (byte)(255 * b)
            };
        }
    
        static double Interpolate(double[] values, int index, double offset) {
            if(offset == 0)
                return values[index];
    
            var x0 = index * LEN_STEP;
            var x1 = x0 + LEN_STEP;
            var y0 = values[index];
            var y1 = values[1 + index];
    
            return y0 + offset * (y1 - y0) / (x1 - x0);
        }
    
        static double GammaCorrect_sRGB(double c) {
            if(c <= 0.0031308)
                return 12.92 * c;
    
            var a = 0.055;
            return (1 + a) * Math.Pow(c, 1 / 2.4) - a;
        }
    
        static double Clip(double c) {
            if(c < 0)
                return 0;
            if(c > 1)
                return 1;
            return c;
        }
    }

另外有一个python版本在这里，暂时用不到，挂个号在这里吧

#!/usr/bin/env python
    # vim:set ft=python fileencoding=utf-8 sr et ts=4 sw=4 : See help 'modeline'
    
    '''
        == A few notes about color ==
    
        Color   Wavelength(nm) Frequency(THz)
        Red     620-750        484-400
        Orange  590-620        508-484
        Yellow  570-590        526-508
        Green   495-570        606-526
        Blue    450-495        668-606
        Violet  380-450        789-668
    
        f is frequency (cycles per second)
        l (lambda) is wavelength (meters per cycle)
        e is energy (Joules)
        h (Plank's constant) = 6.6260695729 x 10^-34 Joule*seconds
                             = 6.6260695729 x 10^-34 m^2*kg/seconds
        c = 299792458 meters per second
        f = c/l
        l = c/f
        e = h*f
        e = c*h/l
    
        List of peak frequency responses for each type of 
        photoreceptor cell in the human eye:
            S cone: 437 nm
            M cone: 533 nm
            L cone: 564 nm
            rod:    550 nm in bright daylight, 498 nm when dark adapted. 
                    Rods adapt to low light conditions by becoming more sensitive.
                    Peak frequency response shifts to 498 nm.
    
    '''
    
    import sys
    import os
    import traceback
    import optparse
    import time
    import logging
    
    
    def wavelength_to_rgb(wavelength, gamma=0.8):
    
        '''This converts a given wavelength of light to an 
        approximate RGB color value. The wavelength must be given
        in nanometers in the range from 380 nm through 750 nm
        (789 THz through 400 THz).
    
        Based on code by Dan Bruton
        http://www.physics.sfasu.edu/astro/color/spectra.html
        '''
    
        wavelength = float(wavelength)
        if wavelength >= 380 and wavelength <= 440:
            attenuation = 0.3 + 0.7 * (wavelength - 380) / (440 - 380)
            R = ((-(wavelength - 440) / (440 - 380)) * attenuation) ** gamma
            G = 0.0
            B = (1.0 * attenuation) ** gamma
        elif wavelength >= 440 and wavelength <= 490:
            R = 0.0
            G = ((wavelength - 440) / (490 - 440)) ** gamma
            B = 1.0
        elif wavelength >= 490 and wavelength <= 510:
            R = 0.0
            G = 1.0
            B = (-(wavelength - 510) / (510 - 490)) ** gamma
        elif wavelength >= 510 and wavelength <= 580:
            R = ((wavelength - 510) / (580 - 510)) ** gamma
            G = 1.0
            B = 0.0
        elif wavelength >= 580 and wavelength <= 645:
            R = 1.0
            G = (-(wavelength - 645) / (645 - 580)) ** gamma
            B = 0.0
        elif wavelength >= 645 and wavelength <= 750:
            attenuation = 0.3 + 0.7 * (750 - wavelength) / (750 - 645)
            R = (1.0 * attenuation) ** gamma
            G = 0.0
            B = 0.0
        else:
            R = 0.0
            G = 0.0
            B = 0.0
        R *= 255
        G *= 255
        B *= 255
        return (int(R), int(G), int(B))
    
    
    def main(options=None, args=None):
    
    #    import ppm_dump
    #    import png_canvas
        import canvas
        if options.ppm:
            canvas = canvas.ppm_canvas(371, 278)
            canvas.is_ascii = True
        else:
            canvas = canvas.png_canvas(371, 278)
        for wl in range(380, 751):
            r, g, b = wavelength_to_rgb(wl)
            for yy in range(0, 278):
                canvas.pixel(wl - 380, yy, r, g, b)
        sys.stdout.write(str(canvas))
    
    if __name__ == '__main__':
        try:
            start_time = time.time()
            parser = optparse.OptionParser(
                formatter=optparse.TitledHelpFormatter(),
                usage=globals()['__doc__'],
                version='1'
            )
            parser.add_option(
                '-v', '--verbose', action='store_true',
                default=False, help='verbose output'
            )
            parser.add_option(
                '--png', action='store_true',
                default=True, help='Output as PNG.'
            )
            parser.add_option(
                '--ppm', action='store_true',
                default=False, help='Output as PPM ASCII (Portable Pixmap).'
            )
            (options, args) = parser.parse_args()
            #if len(args) < 1:
            #    parser.error ('missing argument')
            if options.verbose:
                print(time.asctime())
            exit_code = main(options, args)
            if exit_code is None:
                exit_code = 0
            if options.verbose:
                print(time.asctime())
                print('TOTAL TIME IN MINUTES: %f'
                      % ((time.time() - start_time) / 60.0))
            sys.exit(exit_code)
        except KeyboardInterrupt as e:  # The user pressed Ctrl-C.
            raise e
        except SystemExit as e:  # The script called sys.exit() somewhere.
            raise e
        except Exception as e:
            print('ERROR: Unexpected Exception')
            print(str(e))
            traceback.print_exc()
            os._exit(2)

【1】 [https://stackoverflow.com/questions/1472514/convert-light-frequency-to-rgb][https_stackoverflow.com_questions_1472514_convert-light-frequency-to-rgb]

【2】[http://www.noah.org/wiki/Wavelength\_to\_RGB\_in\_Python][http_www.noah.org_wiki_Wavelength_to_RGB_in_Python]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3Rhbm14MjE5_size_16_color_FFFFFF_t_70]: /images/20220117/0581dc675dc243699f99a03d20a95622.png
[https_stackoverflow.com_questions_1472514_convert-light-frequency-to-rgb]: https://stackoverflow.com/questions/1472514/convert-light-frequency-to-rgb
[http_www.noah.org_wiki_Wavelength_to_RGB_in_Python]: http://www.noah.org/wiki/Wavelength_to_RGB_in_Python