rgb2lab(r, g, b)
{
        local ll, aa, bb, x, y, z, xn, yn, zn;
        local nr = -1;

        if (argc == 1)
        {
                if (length(r) == 3)
                {
                        nr = 3;
                        r = transpose(r);
                }

                if (numcols(r) == 3)
                {
                        g = col(r, 2);
                        b = col(r, 3);
                        r = col(r, 1);
                }
                else        
                {
                        /* assume image */
                        (r, g, b) = getrgb(r);

                        nr = numrows(r);
                }
        }
        else if (argc != 3)
        {
                error(sprintf("%s - 1 or 3 input arguments required", __FUNC__));
        }

        if (max(r) > 1 || max(g) > 1 || max(b) > 1)         
        {
                r = r / 255;
                g = g / 255;
                b = b / 255;
        }

        r = rgb2lab_gamma_correct(r);
        g = rgb2lab_gamma_correct(g);
        b = rgb2lab_gamma_correct(b);

        /* sRGB to XYZ */
        x = r * 0.4124564 + g * 0.3575761 + b * 0.1804375;

//        y = r * 0.2126    + g * 0.7152    + b * 0.0722;
        y = r * 0.2126729 + g * 0.7151521 + b * 0.0721750;

        z = r * 0.0193339 + g * 0.1191920 + b * 0.9503041;


//        x = r * 0.412410846488538800 + g * 0.35758456785295190 + b * 0.18045380393360833;
//        y = r * 0.212649342720652830 + g * 0.71516913570590380 + b * 0.07218152157344333;
//        z = r * 0.019331758429150258 + g * 0.11919485595098397 + b * 0.95039003405033730;

//        x = r * 0.4124564 + g * 0.3575761 + b * 0.1804375;
//        y = r * 0.2126729 + g * 0.7151522 + b * 0.0721750;
//        z = r * 0.0193339 + g * 0.1191920 + b * 0.9503041;

        /* reference white point (D65) */
        xn = 0.95047;
        yn = 1.0;
        zn = 1.08883;

        /* Convert XYZ to XYZ' */
        x = x / xn;
        y = y / yn;
        z = z / zn;
        
          /* 3. Calculate L*, a*, and b* */
        ll = 116 * rgb2lab_f(y) - 16;
        aa = 500 * (rgb2lab_f(x) - rgb2lab_f(y));
        bb = 200 * (rgb2lab_f(y) - rgb2lab_f(z));

        if (outargc <= 1)
        {
                /* cast to series */
                ll = {ll};
                aa = {aa};
                bb = {bb};
        }

        /* label */
        if (nr == 3)
        {
                setvunits(ll, "LAB");
        }
        else
        {
                if (isarray(ll)) setvunits(ll, "L✱");
                if (isarray(aa)) setvunits(aa, "a✱");
                if (isarray(bb)) setvunits(bb, "b✱");

        }

        if (outargc > 1)
        {
                return(ll, aa, bb);
        }
        else if (nr == 3)
        {
                return({ll, aa, bb});
        }
        else
        {
                return(ravel(ll, aa, bb));
        }
}


rgb2lab_gamma_correct(rgb)
{
        rgb = (rgb > 0.04045) ?  ((rgb + 0.055) / 1.055)^2.4 : rgb / 12.92;

        return(rgb);
}


rgb2lab_f0(t)
{
        local f;

        f = (t > 0.008856) ? t ^ (1/3) : 7.787 * t + 16/116;

        return(f);
}


rgb2lab_f(t)
{
        local f;
        static epsilon = (6/29)^3;
        static kappa   = (1/3) * ((29/6) ^ 2);

        f = (t > epsilon) ? t ^ (1/3) : kappa * t + 4/29;

        return(f);
}


rgb2lab_cie(t)
{
        local epsilon, kappa, f;
        static epsilon = 216.0 / 24389.0;
        static kappa   = 24389.0 / 27.0;

        f = (t > epsilon) ? t ^ (1/3) : (kappa * t + 16) / 116;

        return(f);
}