PCA sine curve example using OpenCV

#include "opencv.hpp"

#include <iostream>

#include <time.h>

#define _USE_MATH_DEFINES

#include <math.h>

int main(int argc, char* argv[])

{

cv::Mat mSineWave;

cv::Mat mSineWaveTest;

cv::Mat mShow;

cv::Mat mShowReconstrction;

mSineWave.create(1000, 512, CV_32FC1);

mSineWaveTest.create(100, 512, CV_32FC1);

mShow.create(200, 512, CV_8UC1);

mShowReconstrction.create(200, 512, CV_8UC1);

float freq_mean = 5.f;

float mag_mean = 60.f;

float bias_mean = 100.f;

srand(time(0));

cv::RNG rnd;

for (int i = 0; i < 1000; ++i)

{

float freq = freq_mean + rnd.gaussian(2.f);

float mag = mag_mean + rnd.gaussian(10.f);

float bias = bias_mean + rnd.gaussian(5.f);

float phase = rnd.gaussian(3.f);

for (float x = 0; x < 512.f; x+=1.f)

{

float _x = freq*x*M_PI / 180.f + phase;

float y = mag*sin(_x) + bias;

mSineWave.at<float>(i, (int)x) = y;

}

}

for (int i = 0; i < 100; ++i)

{

float freq = freq_mean + rnd.gaussian(2.2f);

float mag = mag_mean + rnd.gaussian(11.f);

float bias = bias_mean + rnd.gaussian(5.5f);

float phase = rnd.gaussian(3.f);

for (float x = 0; x < 512.f; x += 1.f)

{

float _x = freq*x*M_PI / 180.f + phase;

float y = mag*sin(_x) + bias;

mSineWaveTest.at<float>(i, (int)x) = y;

}

}

cv::PCA pca(mSineWave, cv::Mat(), 0, 0.9999);

pca.mean;// row vector

pca.eigenvalues; // col vector 200 x 1

pca.eigenvectors; // 200 x 512 matrix

std::cout << " # of eigen : " << pca.eigenvalues.rows<< std::endl;

char szT[256] = {0,};

for (int i = 0; i < pca.eigenvalues.rows; ++i)

{

float st_eigen = -3.f*sqrt(abs(pca.eigenvalues.at<float>(i, 0))); // variance -> std deviation

float end_eigen = -st_eigen;

float tick_eigen = (end_eigen - st_eigen) / 50.f;

for (float ev = st_eigen; ev <= end_eigen; ev += tick_eigen)

{

cv::Mat mEV = pca.eigenvectors.row(i);

cv::Mat mEigenSine = pca.mean + ev*mEV;

mShow = 0;

for (float x = 0.f; x < 512; x += 1.f)

{

float y = MAX(0, MIN(199, mEigenSine.at<float>(0, x)));

mShow.at<uchar>((int)y, (int)x) = 255;

}

sprintf(szT, "#  : %d , V : %f ", i, ev);

cv::putText(mShow, szT, cv::Point(5, 13), 2, 0.5, cv::Scalar::all(255));

cv::imshow("sine wave", mShow);

cv::waitKey();

}

}

 

int row = 0;

while (1)

{

mShow = 0;

mShowReconstrction = 0;

cv::Mat mCoef = pca.project( mSineWaveTest.row(row));

cv::Mat mRec = pca.backProject(mCoef);

std::cout << "coefficient" << std::endl;

std::cout << mCoef << std::endl;

for (float x = 0.f; x < 512; x += 1.f)

{

float y = MAX(0, MIN(199, mSineWaveTest.at<float>(row, x)));

//float y = pca.mean.at<float>(0, x); // mean check

mShow.at<uchar>((int)y, (int)x) = 255;

float yr = MAX( 0, MIN(199,mRec.at<float>(0, x)));

mShowReconstrction.at<uchar>((int)yr, (int)x) = 255;

}

printf(" diff = %g\n",  cv::norm(mSineWaveTest.row(row), mRec, cv::NORM_L2));

cv::imshow("sine wave", mShow);

cv::imshow("reconstruction sine wave", mShowReconstrction);

cv::waitKey();

row++;

row %= mSineWaveTest.rows;

}

return 0;

}