サインカーブのフィッティング

Tweet

偏光の解析

物体表面を反射した光を観測すると、偏光板を０°から１７５°まで５°おきに３６枚の画像を取得すると、サインカーブを描く。最小二乗法で画素値からサインカーブを当てはめるプログラムをここに記す。

const double PI2 = 6.28318530717958647692528676655901;
const double EPS = 1.0e-15;

double sin2x[36] =
{0.000000000000000,
0.173648177666930,
0.342020143325669,
0.500000000000000,
0.642787609686539,
0.766044443118978,
0.866025403784439,
0.939692620785908,
0.984807753012208,
1.000000000000000,
0.984807753012208,
0.939692620785908,
0.866025403784439,
0.766044443118978,
0.642787609686539,
0.500000000000000,
0.342020143325669,
0.173648177666931,
0.000000000000000,
-0.173648177666930,
-0.342020143325669,
-0.500000000000000,
-0.642787609686539,
-0.766044443118978,
-0.866025403784438,
-0.939692620785908,
-0.984807753012208,
-1.000000000000000,
-0.984807753012208,
-0.939692620785908,
-0.866025403784439,
-0.766044443118978,
-0.642787609686540,
-0.500000000000000,
-0.342020143325669,
-0.173648177666930};

double cos2x[36] =
{1.000000000000000,
0.984807753012208,
0.939692620785908,
0.866025403784439,
0.766044443118978,
0.642787609686539,
0.500000000000000,
0.342020143325669,
0.173648177666930,
0.000000000000000,
-0.173648177666930,
-0.342020143325669,
-0.500000000000000,
-0.642787609686539,
-0.766044443118978,
-0.866025403784439,
-0.939692620785908,
-0.984807753012208,
-1.000000000000000,
-0.984807753012208,
-0.939692620785908,
-0.866025403784439,
-0.766044443118978,
-0.642787609686539,
-0.500000000000000,
-0.342020143325669,
-0.173648177666930,
-0.000000000000000,
0.173648177666930,
0.342020143325669,
0.500000000000000,
0.642787609686539,
0.766044443118978,
0.866025403784438,
0.939692620785908,
0.984807753012208};

/*******************************************************************/
// y = asin2x + bcos2x + c
// S = Σ(y - asin2x - bcos2x - c)^2 = min を求めて
// y = asin2x + bcos2x + c = Asin(2x - 2B) + C = Asin2(x - B) + C
// A = sqrt(a^2 + b^2)
// sin(-2B) = b / sqrt(a^2 + b^2)
// cos(-2B) = a / sqrt(a^2 + b^2)
// sum1 = Σ
// sumy = Σy
// sums = Σsin2x
// sumc = Σcos2x
// sumys = Σysin2x
// sumyc = Σycos2x
// sumss = Σ(sin2x)^2
// sumsc = Σsin2x cos2x
// sumcc = Σ(cos2x)^2
/*
  numeratora
    = sums * sumc * sumyc
    + sumy * sumc * sumsc
    + sum1 * sumys * sumcc
    - sum1 * sumyc * sumsc
    - sumc * sumc * sumys
    - sumy * sums * sumcc;
  numeratorb
    = sum1 * sumyc * sumss
    + sums * sumc * sumys
    + sumy * sums * sumsc
    - sums * sums * sumyc
    - sumy * sumc * sumss
    - sum1 * sumys * sumsc;
  numeratorc
    = sums * sumsc * sumyc
    + sumc * sumsc * sumys
    + sumy * sumss * sumcc
    - sumc * sumss * sumyc
    - sumy * sumsc * sumsc
    - sums * sumcc * sumys;
  denominator
    = 2.0 * sums * sumc * sumsc
    + sum1 * sumss * sumcc
    - sums * sums * sumcc
    - sumc * sumc * sumss
    - sum1 * sumsc * sumsc;
    */
  /*
  a = numeratora / denominator;
  b = numeratorb / denominator;
  C = numeratorc / denominator;
  */
// 36: 0, 5, 10, 15, ... 175
// 18: 0, 10, 20, 30, ... 170
// sum1 = 36 18
// sumy = ? ?
// sums = 0 0
// sumc = 0 0
// sumys = ? ?
// sumyc = ? ?
// sumss = 18 9
// sumsc = 0 0
// sumcc = 18 9
void CalcSinusoid(double sinusoidy[36], double *resulta, double *resultb, double *resultc)
{
  double sum1;
  double sumy;
  double sumys;
  double sumyc;
  double sumss;
  double sumcc;
  double r;
  int i;
  double y;
  double sx;
  double cx;
  double a, b, A, B, C;
  double asinalpha, acosalpha;
  double minus2B;

  sum1 = 36.0;
  sumy = 0.0;
  sumys = 0.0;
  sumyc = 0.0;
  sumss = 18.0;
  sumcc = 18.0;
  r = 5.0;

  // 各種の和を求める
  for(i = 0; i < 36; i++)
  {
    y = sinusoidy[i];
    sx = sin2x[i];
    cx = cos2x[i];
    sumy += y;
    sumys += y * sx;
    sumyc += y * cx;
  }

  a = sumys / sumss;
  b = sumyc / sumcc;
  C = sumy / sum1;

  if(C == 0.0)
  {
    A = B = 0.0;
    C = EPS;
  }
  else if(a == 0.0 && b == 0.0)
  {
    A = B = 0.0;
  }
  else
  {
    A = sqrt(a * a + b * b);
    asinalpha = asin(b / A);
    acosalpha = acos(a / A);
    if(asinalpha < 0.0) minus2B = PI2 - acosalpha;
    else minus2B = acosalpha;
    B = -0.5 * minus2B;
  }

  // A, B, C の値を返す
  *resulta = A;
  *resultb = B;
  *resultc = C;
}

もどる