rgb2lab2(r, g, b)
{
        local rgb, ll, aa, bb, m, xyz, xyz_norm, xw, yw, zw;
        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 /= 255;
                g /= 255;
                b /= 255;
        }

        rgb = ravel(r[..], g[..], b[..]);

        /* 1. Convert RGB to XYZ */
        rgb = (rgb > 0.04045) ?  ((rgb + 0.055) / 1.055)^2.4 : rgb / 12.92;

        m = {{0.4124, 0.3576, 0.1805},
             {0.2126, 0.7152, 0.0722},
             {0.0193, 0.1192, 0.9505}};

        m = m'; //

        xyz = rgb *^ m;

        /* Convert XYZ to XYZ' */
        (xw, yw, zw) = (95.047, 100.000, 108.883);  // D65 white point
        xyz_norm = xyz / transpose({xw, yw, zw});

          /* 3. Calculate L*, a*, and b* */
        ll = 116 * rgb2lab2_cie(xyz_norm[.., 2]) - 16;
        aa = 500 * (rgb2lab2_cie(xyz_norm[.., 1]) - rgb2lab2_cie(xyz_norm[.., 2]));
        bb = 200 * (rgb2lab2_cie(xyz_norm[.., 2]) - rgb2lab2_cie(xyz_norm[.., 3]));

        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 (nr > 1)
                {
                        ll = ravel(ll, nr);
                        aa = ravel(aa, nr);
                        bb = ravel(bb, nr);
                }
        }

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



rgb2lab2_f(t)
{
        local f;

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

        return(f);
}


rgb2lab2_cie(t)
{
        local epsilon, kappa, f;

        epsilon = 216.0 / 24389.0;
        kappa = 24389.0 / 27.0;

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

        return(f);
}