color.cc

一个用MATLAB语言编写的摄像机标定工具箱,内容丰富

// 
// Color.cc
//
// $Id: color.cc,v 1.1.1.1 2001/02/28 00:28:35 cstolte Exp $
//

#include <sgl/color.h>
#include <sgl/rgb.h>
#include <sgl/matrix.h>

const Primary SMPTE_Primary(
	Chroma(.635, .340),	/* R */
	Chroma(.305, .595),	/* G */
	Chroma(.155, .070),	/* B */
	Chroma(.3127, .3291)	/* W */
);

//
// Chroma
//

Chroma::Chroma(const XYZ &XYZ) {
  double t = XYZ.X() + XYZ.Y() + XYZ.Z();

  if (zero(t))
      X = Y = 0;
  else {
      X = XYZ.X() / t;
      Y = XYZ.Y() / t;
  }
}

Chroma 
Chroma::operator+(const Chroma &b) const {
  return Chroma(X+b.X, Y+b.Y);
}

Chroma &
Chroma::operator+=(const Chroma &b) {
  X += b.X;
  Y += b.Y;
  return *this;
}

int 
Chroma::operator==(const Chroma &b) const {
  return (X == b.X && Y == b.Y);
}

//
// Primary
//

Primary::Primary(const Chroma &r, const Chroma &g, 
			const Chroma &b, const Chroma &w)
  : white(w) {
    p[0] = r;
    p[1] = g;
    p[2] = b;
    init();
}

Primary::Primary(const Chroma c[4])
  : white(c[3]) {
    p[0] = c[0];
    p[1] = c[1];
    p[2] = c[2];
    init();
}

XYZ 
Primary::toXYZ(const RGBPrimary &c) const {
#ifdef __GNUG__
    Vec<3, double> cvec;
    cvec[0] = c.R();
    cvec[1] = c.G();
    cvec[2] = c.B();
    Vec<3, double> ret = cvec * RGB2XYZ;
    return XYZ(ret[0], ret[1], ret[2]);
#else
    return XYZ(c.R()*RGB2XYZ[0][0] + c.G()*RGB2XYZ[0][1] +
	       c.B()*RGB2XYZ[0][2],
	       c.R()*RGB2XYZ[1][0] + c.G()*RGB2XYZ[1][1] +
	       c.B()*RGB2XYZ[1][2],
	       c.R()*RGB2XYZ[2][0] + c.G()*RGB2XYZ[2][1] +
	       c.B()*RGB2XYZ[2][2]);
#endif // __GNUG__
}

RGBPrimary 
Primary::RGB(const XYZ &c) const {
#ifdef __GNUG__
    Vec<3, double> cvec;
    cvec[0] = c.X();
    cvec[1] = c.Y();
    cvec[2] = c.Z();
    Vec<3, double> ret = cvec * XYZ2RGB;
    return RGBPrimary(ret[0], ret[1], ret[2], this);
#else
    return RGBPrimary(c.X()*XYZ2RGB[0][0] + c.Y()*XYZ2RGB[0][1] +
		      c.Z()*XYZ2RGB[0][2],
		      c.X()*XYZ2RGB[1][0] + c.Y()*XYZ2RGB[1][1] + 
		      c.Z()*XYZ2RGB[1][2],
		      c.X()*XYZ2RGB[2][0] + c.Y()*XYZ2RGB[2][1] + 
		      c.Z()*XYZ2RGB[2][2],
		      this);
#endif // __GNUG__
}

void
Primary::init()
{
  int i, j, k, idx[3], tmp;
  double inv[3][3], w[3], s[3], fac;

  /*
   * Build RGB->XYZ and inverse transformation matrices.
   *
   * RGB->XYZ is R/G/B Chromaticities in matrix form, scaled
   * by sr, sg, sb as to ensure that (1,1,1) transforms to
   * normalized XYZ of white.   White XYZ is normalized to have
   * Y value of 1.0.
   *
   * XYZ->RGB is the inverse of the above.
   */

  if (white.y() < 0.0 || zero(white.y())) {
      /* punt */
      return;
  }

  for (i = 0; i < 3; i++) {
      idx[i] = i;
      inv[0][i] = p[i].x();
      inv[1][i] = p[i].y();
      inv[2][i] = p[i].z(); 
  }

  w[0] = white.x() / white.y();
  w[1] = 1.0;
  w[2] = white.z() / white.y();

  for (i = 0; i < 2; i++)  {
      for (j = i+1; j < 3; j++) {
	  if (fabs(inv[idx[j]][i]) > fabs(inv[idx[i]][i])) {
	      tmp = idx[i];
	      idx[i] = idx[j];
	      idx[j] = tmp;
	  }
      }
    
      assert(!zero(inv[idx[i]][i])); /* singular */

      for (j = i+1; j < 3; j++) {
	  fac = inv[idx[j]][i] / inv[idx[i]][i];
	  w[idx[j]] -= w[idx[i]]*fac;
	  for (k = i+1; k < 3; k++)
	      inv[idx[j]][k] -= inv[idx[i]][k] * fac;
      }
  }
  assert(!zero(inv[idx[2]][2])); /* singular */

  s[idx[2]] = w[idx[2]] / inv[idx[2]][2];
  s[idx[1]] = (w[idx[1]] - (inv[idx[1]][2] * s[idx[2]]))/inv[idx[1]][1];
  s[idx[0]] = (w[idx[0]] - (inv[idx[0]][1] * s[idx[1]]) -
	       (inv[idx[0]][2] * s[idx[2]])) / inv[idx[0]][0];

  for (i = 0; i < 3; i++) {
      RGB2XYZ[0][i] = p[i].x() * s[i];
      RGB2XYZ[1][i] = p[i].y() * s[i];
      RGB2XYZ[2][i] = p[i].z() * s[i];
  }

#ifdef __GNUG__
  XYZ2RGB = RGB2XYZ.inverse();
#else
  i = InverseMatrix3x3(RGB2XYZ, XYZ2RGB);
  assert(i);
#endif // __GNUG__
}

//
// XYZ
//

XYZ::XYZ(const RGBPrimary &p) { 
    *this = p.toXYZ(); 
}

int 
XYZ::operator==(const XYZ &b) const {
  return (x == b.x && y == b.y && z == b.z);
}

XYZ 
XYZ::operator+(const XYZ &b) const {
  return XYZ(x + b.x, y + b.y, z + b.z);
}

XYZ &
XYZ::operator+=(const XYZ &b) {
    x += b.x;
    y += b.y;
    z += b.z;
    return *this;
}

XYZ
XYZ::operator*(double a) const {
    return XYZ(x * a, y * a, z * a);
}

XYZ &
XYZ::operator*=(double a) {
    x *= a;
    y *= a;
    z *= a;
    return *this;
}

XYZ
operator*(double a, const XYZ &xyz) {
    return XYZ(xyz.x * a, xyz.y * a, xyz.z * a);
}

//
// Lab
//

Lab::Lab(const XYZ &xyz, const XYZ &white) {
    double xf = xyz.X() / white.X();
    double yf = xyz.Y() / white.Y();
    double zf = xyz.Z() / white.Z();

    if (xf > 0.008856)
	xf = pow(xf, 1./3.);
    else
	xf = 7.787 * xf + 16.0/116.0;

    if (yf > 0.008856) {
	yf = pow(yf, 1./3.);
	L = 116.0 * yf - 16.0;
    }
    else {
	  L = 903.3 * yf;
	  yf = 7.787 * yf + 16./116.;
      }

    if (zf > 0.008856)
	zf = pow(zf, 1./3.);
    else
	zf = 7.787 * zf + 16./116.;

    A = 500. * (xf - yf);
    B = 200. * (yf - zf);
}


Lab::Lab(double l, double a, double b, const XYZ &w)
  : L(l), A(a), B(b), white(w) {
}


Lab::Lab(double Lab[3], const XYZ &w)
  : L(Lab[0]), A(Lab[1]), B(Lab[2]), white(w) {
}


Lab::Lab(const Lab &, const XYZ &) {
  // ???
}

Lab 
Lab::operator+(const Lab &c) const {
    return Lab(L + c.L, A + c.A, B + c.B, white + c.white);
}

Lab &
Lab::operator+=(const Lab &c) {
    L += c.L;
    A += c.A;
    B += c.B;
    white += c.white;
    return *this;
}

Lab 
Lab::operator*(double) const {
    // ???
    assert(1 == 0);
    return Lab(*this, XYZ());
}

Lab &
Lab::operator*=(double) {
    // ???
    assert(1 == 0);
    return *this;
}

Lab 
Lab::operator/(double) const {
    // ???
    assert(1 == 0);
    return Lab(*this, XYZ());
}

Lab &
Lab::operator/=(double) {
    // ???
    assert(1 == 0);
    return *this;
}

int 
Lab::operator==(const Lab &) const {
    // ???
    assert(1 == 0);
    return 0;
}

//
// Luv
//

Luv::Luv(const XYZ &xyz, const XYZ &white) {
    double d = xyz.X() + 15*xyz.Y() + 3*xyz.Z();
    double dn = white.X() + 15*white.Y() + 3*white.Z();
    
    if (zero(white.Y()) || zero(d) || zero(dn)) 
      {
	  L = U = V = 0.0;
	  return;
      }
    
    double rat = xyz.Y() / white.Y();
    
    if (rat <= 0.008856) 
	L = 903.3 * rat;
    else 
	L = 116. * pow(rat, 1/3) - 16.;
    
    double up = 4 * xyz.X() / d;
    double upn = 4 * white.X() / dn;
    double vp = 9 * xyz.Y() / d;
    double vpn = 9 * white.Y() / dn;
    
    U = 13 * L *(up - upn);
    V = 13 * L *(vp - vpn);
}


Luv::Luv(double l, double u, double v, const XYZ &w)
: L(l), U(u), V(v), white(w) {
}


Luv::Luv(double Luv[3], const XYZ &w)
: L(Luv[0]), U(Luv[1]), V(Luv[2]), white(w) {
}


Luv::Luv(const Luv &, const XYZ &) {
    // ???
    assert(1 == 0);
}

int 
Luv::operator==(const Luv &) const {
    // ???
    // what if l.white != this->white??
    assert(1 == 0);
    return 0;
}

Luv 
Luv::operator+(const Luv &l) const {
    return Luv(L + l.L, U + l.U, V + l.V, white + l.white);
}

Luv &
Luv::operator+=(const Luv &l) {
    L += l.L;
    U += l.U;
    V += l.V;
    white += l.white;
    return *this;
}

Luv 
Luv::operator*(double) const {
    assert(1 == 0);
    return Luv(*this, XYZ());
}

Luv &
Luv::operator*=(double) {
    assert(1 == 0);
    return *this;
}