ulabrgb.pas

图形中LAB的算法

unit uLabRgb;

interface

type
  TVector3 = array[1..3] of Double;

function LabToRgb(Lab: TVector3): TVector3;
function RgbToLab(Rgb: TVector3): TVector3;

implementation

type
  TMatrix3 = array[1..3, 1..3] of Double;

var
  RgbXyz: TMatrix3 = ((1, 0, 0), (0, 1, 0), (0, 0, 1));
  XyzRgb: TMatrix3 = ((1, 0, 0), (0, 1, 0), (0, 0, 1));

  { CCIR recommended values }
  PhosphorX: TVector3 = (0.64, 0.30, 0.15);
  PhosphorY: TVector3 = (0.33, 0.60, 0.06);
  WhitePoint: TVector3 = (0.95, 1.0000, 1.09);
  Gamma: Double = 1 / 0.45;

function MultiplyMatrix3ByVector3(const M: TMatrix3;
  const V: TVector3): TVector3;
var
  I: Integer;
  J: Integer;
begin
  for I := 1 to 3 do begin
    Result[I] := 0.0;
    for J := 1 to 3 do
      Result[I] := Result[I] + M[I, J] * V[J] end end;

function MultiplyMatrix3ByMatrix3(const M1, M2: TMatrix3): TMatrix3;
var
  I: Integer;
  J: Integer;
  K: Integer;
begin
  for I := 1 to 3 do
    for J := 1 to 3 do begin
      Result[I, J] := 0.0;
      for K := 1 to 3 do
        Result[I, J] := Result[I, J] + M1[I, K] * M2[K, J] end end;

function InvertMatrix3(const M: TMatrix3): TMatrix3;
var
  I: Integer;
  J: Integer;
  D: Double;

  function Next(I: Integer): Integer;
  begin
    Result := I + 1;
    if Result > 3 then
      Result := Result - 3 end;

  function Prev(I: Integer): Integer;
  begin
    Result := I - 1;
    if Result < 1 then
      Result := Result + 3 end;

begin
  D := 0;
  for I := 1 to 3 do
    D := D + M[1, I] * (M[2, Next(I)] * M[3, Prev(I)] -
      M[2, Prev(I)] * M[3, Next(I)]);
  FillChar(Result, SizeOf(Result), 0);
  for I := 1 to 3
    do for J := 1 to 3 do
      Result[J, I] := (M[Next(I), Next(J)] * M[Prev(I), Prev(J)] -
        M[Next(I), Prev(J)] * M[Prev(I), Next(J)]) / D end;

function LabToXyz(const Lab: TVector3): TVector3;
var
  LL: Double;

  function Cube(X: Double): Double;
  begin
    if X >= (6 / 29) then
      Result := X * X * X
    else
      Result := (108 / 841) * (X - (4 / 29)) end;

begin
  LL := (Lab[1] + 16) / 116;
  Result[1] := WhitePoint[1] * Cube(LL + Lab[2] / 500);
  Result[2] := WhitePoint[2] * Cube(LL);
  Result[3] := WhitePoint[3] * Cube(LL - Lab[3] / 200) end;

function XyzToRgb(const Xyz: TVector3): TVector3;
var
  I: Integer;
begin
  Result := MultiplyMatrix3ByVector3(XyzRgb, Xyz);
  for I := 1 to 3 do
    if Result[I] <= 0.0 then
      Result[I] := 0
    else Result[I] := Exp(Ln(Result[I]) / Gamma) end;

function LabToRgb(Lab: TVector3): TVector3;
begin
  Result := XyzToRgb(LabToXyz(Lab)) end;

function RgbToXyz(const Rgb: TVector3): TVector3;
var
  I: Integer;
begin
  Result := Rgb;
  for I := 1 to 3 do
    if Result[I] <= 0.0 then
      Result[I] := 0
    else Result[I] := Exp(Ln(Result[I]) * Gamma);
  Result := MultiplyMatrix3ByVector3(RgbXyz, Result) end;

function XyzToLab(const Xyz: TVector3): TVector3;
var
  YY: Double;

  function CubeRoot(X: Double): Double;
  begin
    if X >= (216 / 24389) then
      Result := Exp(Ln(X) / 3)
    else
      Result := (841 / 108) * X + (4 / 29) end;

begin
  YY := CubeRoot(Xyz[2] / WhitePoint[2]);
  Result[1] := 116 * YY - 16;
  Result[2] := 500 * (CubeRoot(Xyz[1] / WhitePoint[1]) - YY);
  Result[3] := 200 * (YY - CubeRoot(Xyz[3] / WhitePoint[3])) end;

function RgbToLab(Rgb: TVector3): TVector3;
begin
  Result := XyzToLab(RgbToXyz(Rgb)) end;

procedure InitTransformationMatrices;
var
  I: Integer;
  J: Integer;
  C: TVector3;
  CToXyz: TMatrix3;
  XyzToC: TMatrix3;
begin
  for I := 1 to 3 do begin
    CToXyz[1, I] := PhosphorX[I];
    CToXyz[2, I] := PhosphorY[I];
    CToXyz[3, I] := 1 - PhosphorX[I] - PhosphorY[I] end;
  XyzToC := InvertMatrix3(CToXyz);
  C := MultiplyMatrix3ByVector3(XyzToC, WhitePoint);

  for I := 1 to 3 do
    for J := 1 to 3 do
      RgbXyz[I, J] := CToXyz[I, J] * C[J];

  XyzRgb := InvertMatrix3(RgbXyz)
end;

initialization
  InitTransformationMatrices

end.