#include "kalman.h"
#include <stdio.h>


/* tester de printer toutes les valeurs des vecteurs...*/
/* tester de changer les matrices du noises */
/* replace state by cvkalman->state_post ??? */


CvRandState rng;
const double T = 0.1;
kalman::kalman(int x,int xv,int y,int yv)
{					
    cvkalman = cvCreateKalman( 4, 4, 0 );
    state = cvCreateMat( 4, 1, CV_32FC1 );
    process_noise = cvCreateMat( 4, 1, CV_32FC1 );
    measurement = cvCreateMat( 4, 1, CV_32FC1 );
    int code = -1;
    
    /* create matrix data */
     const float A[] = { 
			1, T, 0, 0,
			0, 1, 0, 0,
			0, 0, 1, T,
			0, 0, 0, 1
		};
					
     const float H[] = { 
 			1, 0, 0, 0,
 			0, 0, 0, 0,
			0, 0, 1, 0,
			0, 0, 0, 0
		};
		     
     const float P[] = {
 			pow(320,2), pow(320,2)/T, 0, 0,
			pow(320,2)/T, pow(320,2)/pow(T,2), 0, 0,
			0, 0, pow(240,2), pow(240,2)/T,
			0, 0, pow(240,2)/T, pow(240,2)/pow(T,2)
 			};

     const float Q[] = {
			pow(T,3)/3, pow(T,2)/2, 0, 0,
			pow(T,2)/2, T, 0, 0,
			0, 0, pow(T,3)/3, pow(T,2)/2,
			0, 0, pow(T,2)/2, T
			};
			
     const float R[] = {
			1, 0, 0, 0,
			0, 0, 0, 0,
			0, 0, 1, 0,
			0, 0, 0, 0
			};
			
    
    cvRandInit( &rng, 0, 1, -1, CV_RAND_UNI );

    cvZero( measurement );
    
    cvRandSetRange( &rng, 0, 0.1, 0 );
    rng.disttype = CV_RAND_NORMAL;

    cvRand( &rng, state );

    memcpy( cvkalman->transition_matrix->data.fl, A, sizeof(A));
    memcpy( cvkalman->measurement_matrix->data.fl, H, sizeof(H));
    memcpy( cvkalman->process_noise_cov->data.fl, Q, sizeof(Q));
    memcpy( cvkalman->error_cov_post->data.fl, P, sizeof(P));
    memcpy( cvkalman->measurement_noise_cov->data.fl, R, sizeof(R));
    //cvSetIdentity( cvkalman->process_noise_cov, cvRealScalar(1e-5) );    
    //cvSetIdentity( cvkalman->error_cov_post, cvRealScalar(1));
	//cvSetIdentity( cvkalman->measurement_noise_cov, cvRealScalar(1e-1) );

    /* choose initial state */

    state->data.fl[0]=x;
    state->data.fl[1]=xv;
    state->data.fl[2]=y;
    state->data.fl[3]=yv;
    cvkalman->state_post->data.fl[0]=x;
    cvkalman->state_post->data.fl[1]=xv;
    cvkalman->state_post->data.fl[2]=y;
    cvkalman->state_post->data.fl[3]=yv;

	cvRandSetRange( &rng, 0, sqrt(cvkalman->process_noise_cov->data.fl[0]), 0 );
    cvRand( &rng, process_noise );


    }

     
CvPoint2D32f kalman::get_predict(float x, float y){
    

    /* update state with current position */
    state->data.fl[0]=x;
    state->data.fl[2]=y;

    
    /* predict point position */
    /* x'k=A?xk+B?uk
       P'k=A?Pk-1*AT + Q */
    cvRandSetRange( &rng, 0, sqrt(cvkalman->measurement_noise_cov->data.fl[0]), 0 );
    cvRand( &rng, measurement );
    
     /* xk=A?xk-1+B?uk+wk */
    cvMatMulAdd( cvkalman->transition_matrix, state, process_noise, cvkalman->state_post );
    
    /* zk=H?xk+vk */
    cvMatMulAdd( cvkalman->measurement_matrix, cvkalman->state_post, measurement, measurement );
    
    /* adjust Kalman filter state */
    /* Kk=P'k?HT?(H?P'k?HT+R)-1
       xk=x'k+Kk?(zk-H?x'k)
       Pk=(I-Kk?H)?P'k */
    cvKalmanCorrect( cvkalman, measurement );
    float measured_value_x = measurement->data.fl[0];
    float measured_value_y = measurement->data.fl[2];

    
	const CvMat* prediction = cvKalmanPredict( cvkalman, 0 );
    float predict_value_x = prediction->data.fl[0];
    float predict_value_y = prediction->data.fl[2];

    return(cvPoint2D32f(predict_value_x,predict_value_y));
}

void kalman::init_kalman(int x,int xv,int y,int yv)
{
	state->data.fl[0]=x;
    state->data.fl[1]=xv;
    state->data.fl[2]=y;
    state->data.fl[3]=yv;
    cvkalman->state_post->data.fl[0]=x;
    cvkalman->state_post->data.fl[1]=xv;
    cvkalman->state_post->data.fl[2]=y;
    cvkalman->state_post->data.fl[3]=yv;
}