xingma0910

学习中……

导航

Cv运动分析与对象跟踪(转)

Cv运动分析与对象跟踪

 

目录

[隐藏]

 

背景统计量的累积

Acc

将帧叠加到累积器(accumulator)中

void cvAcc( const CvArr* image, CvArr* sum, const CvArr* mask=NULL );
image
输入图像, 1- 或 3-通道, 8-比特或32-比特浮点数. (多通道的每一个通道都单独处理).
sum
同一个输入图像通道的累积,32-比特或64-比特浮点数数组.
mask
可选的运算 mask.

函数 cvAcc 将整个图像 image 或某个选择区域叠加到 sum 中:

sum(x,y)=sum(x,y)+image(x,y) if mask(x,y)!=0

SquareAcc

叠加输入图像的平方到累积器中

void cvSquareAcc( const CvArr* image, CvArr* sqsum, const CvArr* mask=NULL );
image
输入图像, 1- 或 3-通道, 8-比特或32-比特浮点数 (多通道的每一个通道都单独处理)
sqsum
同一个输入图像通道的累积,32-比特或64-比特浮点数数组.
mask
可选的运算 mask.

函数 cvSquareAcc 叠加输入图像 image 或某个选择区域的二次方,到累积器 sqsum 中

sqsum(x,y)=sqsum(x,y)+image(x,y)2 if mask(x,y)!=0

MultiplyAcc

将两幅输入图像的乘积叠加到累积器中

void cvMultiplyAcc( const CvArr* image1, const CvArr* image2, CvArr* acc, const CvArr* mask=NULL );
image1
第一个输入图像, 1- or 3-通道, 8-比特 or 32-比特 浮点数 (多通道的每一个通道都单独处理)
image2
第二个输入图像, 与第一个图像的格式一样
acc
同一个输入图像通道的累积,32-比特或64-比特浮点数数组.
mask
可选的运算 mask.

函数 cvMultiplyAcc 叠加两个输入图像的乘积到累积器 acc:

acc(x,y)=acc(x,y) + image1(x,y)•image2(x,y) if mask(x,y)!=0

RunningAvg

更新 running average 滑动平均( Hunnish: 不知道 running average 如何翻译才恰当)

void cvRunningAvg( const CvArr* image, CvArr* acc, double alpha, const CvArr* mask=NULL );
image
输入图像, 1- or 3-通道, 8-比特 or 32-比特 浮点数 (each channel of multi-channel image is processed independently).
acc
同一个输入图像通道的累积,32-比特或64-比特浮点数数组.
alpha
输入图像权重
mask
可选的运算 mask

函数 cvRunningAvg 计算输入图像 image 的加权和,以及累积器 acc 使得 acc 成为帧序列的一个 running average:

acc(x,y)=(1-α)•acc(x,y) + α•image(x,y) if mask(x,y)!=0

其中 α (alpha) 调节更新速率 (累积器以多快的速率忘掉前面的帧).

运动模板

UpdateMotionHistory

去掉影像(silhouette) 以更新运动历史图像

void cvUpdateMotionHistory( const CvArr* silhouette, CvArr* mhi,
                            double timestamp, double duration );
silhouette
影像 mask,运动发生地方具有非零象素
mhi
运动历史图像(单通道, 32-比特 浮点数),为本函数所更新
timestamp
当前时间,毫秒或其它单位
duration
运动跟踪的最大持续时间,用 timestamp 一样的时间单位

函数 cvUpdateMotionHistory 用下面方式更新运动历史图像:

mhi(x,y)=timestamp  if silhouette(x,y)!=0
         0          if silhouette(x,y)=0 and mhi(x,y)<timestamp-duration
         mhi(x,y)   otherwise

也就是,MHI(motion history image) 中在运动发生的象素点被设置为当前时间戳,而运动发生较久的象素点被清除。

CalcMotionGradient

计算运动历史图像的梯度方向

void cvCalcMotionGradient( const CvArr* mhi, CvArr* mask, CvArr* orientation,
                           double delta1, double delta2, int aperture_size=3 );
mhi
运动历史图像
mask
Mask 图像;用来标注运动梯度数据正确的点,为输出参数。
orientation
运动梯度的方向图像,包含从 0 到 360 角度
delta1, delta2
函数在每个象素点 (x,y) 邻域寻找 MHI 的最小值 (m(x,y)) 和最大值 (M(x,y)),并且假设梯度是正确的,当且仅当:
min(delta1,delta2) <= M(x,y)-m(x,y) <= max(delta1,delta2).
aperture_size
函数所用微分算子的开孔尺寸 CV_SCHARR, 1, 3, 5 or 7 (见 cvSobel).

函数 cvCalcMotionGradient 计算 MHI 的差分 Dx 和 Dy ,然后计算梯度方向如下式:

orientation(x,y)=arctan(Dy(x,y)/Dx(x,y))

其中都要考虑 Dx(x,y)' 和 Dy(x,y)' 的符号 (如 cvCartToPolar 类似). 然后填充 mask 以表示哪些方向是正确的(见 delta1 和delta2 的描述).

CalcGlobalOrientation

计算某些选择区域的全局运动方向

double cvCalcGlobalOrientation( const CvArr* orientation, const CvArr* mask, const CvArr* mhi,
                                double timestamp, double duration );
orientation
运动梯度方向图像,由函数 cvCalcMotionGradient 得到
mask
Mask 图像. 它可以是正确梯度 mask (由函数 cvCalcMotionGradient 得到)与区域 mask 的结合,其中区域 mask 确定哪些方向需要计算。
mhi
运动历史图象
timestamp
当前时间(单位毫秒或其它)最好在传递它到函数 cvUpdateMotionHistory 之前存储一下以便以后的重用,因为对大图像运行 cvUpdateMotionHistory 和 cvCalcMotionGradient 会花费一些时间
duration
运动跟踪的最大持续时间,用法与 cvUpdateMotionHistory 中的一致

函数 cvCalcGlobalOrientation 在选择的区域内计算整个运动方向,并且返回 0° 到 360° 之间的角度值。首先函数创建运动直方图,寻找基本方向做为直方图最大值的坐标。然后函数计算与基本方向的相对偏移量,做为所有方向向量的加权和:运行越近,权重越大。得到的角度是基本方向和偏移量的循环和。

SegmentMotion

将整个运动分割为独立的运动部分

CvSeq* cvSegmentMotion( const CvArr* mhi, CvArr* seg_mask, CvMemStorage* storage,
                        double timestamp, double seg_thresh );
mhi
运动历史图像
seg_mask
发现应当存储的 mask 的图像, 单通道, 32bits, 浮点数.
storage
包含运动连通域序列的内存存储仓
timestamp
当前时间,毫秒单位
seg_thresh
分割阈值,推荐等于或大于运动历史“每步”之间的间隔。

函数 cvSegmentMotion 寻找所有的运动分割,并且在seg_mask 用不同的单独数字(1,2,...)标识它们。它也返回一个具有 CvConnectedComp 结构的序列,其中每个结构对应一个运动部件。在这之后,每个运动部件的运动方向就可以被函数 cvCalcGlobalOrientation 利用提取的特定部件的掩模(mask)计算出来(使用 cvCmp)

对象跟踪

MeanShift

在反向投影图中发现目标中心

int cvMeanShift( const CvArr* prob_image, CvRect window,
                 CvTermCriteria criteria, CvConnectedComp* comp );
prob_image
目标直方图的反向投影(见 cvCalcBackProject).
window
初始搜索窗口
criteria
确定窗口搜索停止的准则
comp
生成的结构,包含收敛的搜索窗口坐标 (comp->rect 字段) 与窗口内部所有象素点的和 (comp->area 字段).

函数 cvMeanShift 在给定反向投影和初始搜索窗口位置的情况下,用迭代方法寻找目标中心。当搜索窗口中心的移动小于某个给定值时或者函数已经达到最大迭代次数时停止迭代。 函数返回迭代次数。

CamShift

发现目标中心,尺寸和方向

int cvCamShift( const CvArr* prob_image, CvRect window, CvTermCriteria criteria,
                CvConnectedComp* comp, CvBox2D* box=NULL );
prob_image
目标直方图的反向投影 (见 cvCalcBackProject).
window
初始搜索窗口
criteria
确定窗口搜索停止的准则
comp
生成的结构,包含收敛的搜索窗口坐标 (comp->rect 字段) 与窗口内部所有象素点的和 (comp->area 字段).
box
目标的带边界盒子。如果非 NULL, 则包含目标的尺寸和方向。

函数 cvCamShift 实现了 CAMSHIFT 目标跟踪算法([Bradski98]). 首先它调用函数 cvMeanShift 寻找目标中心,然后计算目标尺寸和方向。最后返回函数 cvMeanShift 中的迭代次数。

CvCamShiftTracker 类在 cv.hpp 中被声明,函数实现了彩色目标的跟踪。

SnakeImage

改变轮廓位置使得它的能量最小

void cvSnakeImage( const IplImage* image, CvPoint* points, int length,
                   float* alpha, float* beta, float* gamma, int coeff_usage,
                   CvSize win, CvTermCriteria criteria, int calc_gradient=1 );
image
输入图像或外部能量域
points
轮廓点 (snake).
length
轮廓点的数目
alpha
连续性能量的权 Weight[s],单个浮点数或长度为 length 的浮点数数组,每个轮廓点有一个权
beta
曲率能量的权 Weight[s],与 alpha 类似
gamma
图像能量的权 Weight[s],与 alpha 类似
coeff_usage
前面三个参数的不同使用方法:
  • CV_VALUE 表示每个 alpha, beta, gamma 都是指向为所有点所用的一个单独数值;
  • CV_ARRAY 表示每个 alpha, beta, gamma 是一个指向系数数组的指针,snake 上面各点的系数都不相同。因此,各个系数数组必须与轮廓具有同样的大小。所有数组必须与轮廓具有同样大小
win
每个点用于搜索最小值的邻域尺寸,两个 win.width 和 win.height 都必须是奇数
criteria
终止条件
calc_gradient
梯度符号。如果非零,函数为每一个图像象素计算梯度幅值,且把它当成能量场,否则考虑输入图像本身。

函数 cvSnakeImage 更新 snake 是为了最小化 snake 的整个能量,其中能量是依赖于轮廓形状的内部能量(轮廓越光滑,内部能量越小)以及依赖于能量场的外部能量之和,外部能量通常在哪些局部能量极值点中达到最小值(这些局部能量极值点与图像梯度表示的图像边缘相对应)。

参数 criteria.epsilon 用来定义必须从迭代中除掉以保证迭代正常运行的点的最少数目。

如果在迭代中去掉的点数目小于 criteria.epsilon 或者函数达到了最大的迭代次数 criteria.max_iter ,则终止函数。

光流

CalcOpticalFlowHS

计算两幅图像的光流

void cvCalcOpticalFlowHS( const CvArr* prev, const CvArr* curr, int use_previous,
                          CvArr* velx, CvArr* vely, double lambda,
                          CvTermCriteria criteria );
prev
第一幅图像, 8-比特, 单通道.
curr
第二幅图像, 8-比特, 单通道.
use_previous
使用以前的 (输入) 速度域
velx
光流的水平部分,与输入图像大小一样, 32-比特,浮点数, 单通道.
vely
光流的垂直部分,与输入图像大小一样, 32-比特, 浮点数, 单通道.
lambda
Lagrangian 乘子
criteria
速度计算的终止条件

函数 cvCalcOpticalFlowHS 为输入图像的每一个象素计算光流,使用 Horn & Schunck 算法 [Horn81].

CalcOpticalFlowLK

计算两幅图像的光流

void cvCalcOpticalFlowLK( const CvArr* prev, const CvArr* curr, CvSize win_size,
                          CvArr* velx, CvArr* vely );
prev
第一幅图像, 8-比特, 单通道.
curr
第二幅图像, 8-比特, 单通道.
win_size
用来归类象素的平均窗口尺寸 (Size of the averaging window used for grouping pixels)
velx
光流的水平部分,与输入图像大小一样, 32-比特, 浮点数, 单通道.
vely
光流的垂直部分,与 输入图像大小一样, 32-比特, 浮点数, 单通道.

函数 cvCalcOpticalFlowLK 为输入图像的每一个象素计算光流,使用 Lucas & Kanade 算法 [Lucas81].

CalcOpticalFlowBM

用块匹配方法计算两幅图像的光流

void cvCalcOpticalFlowBM( const CvArr* prev, const CvArr* curr, CvSize block_size,
                          CvSize shift_size, CvSize max_range, int use_previous,
                          CvArr* velx, CvArr* vely );
prev
第一幅图像, 8-比特, 单通道.
curr
第二幅图像, 8-比特, 单通道.
block_size
比较的基本块尺寸
shift_size
块坐标的增量
max_range
块周围象素的扫描邻域的尺寸
use_previous
使用以前的 (输入) 速度域
velx
光流的水平部分,尺寸为 floor((prev->width - block_size.width)/shiftSize.width) × floor((prev->height - block_size.height)/shiftSize.height) , 32-比特,浮点数, 单通道.
vely
光流的垂直部分,与 velx 大小一样,32-比特,浮点数, 单通道.

函数 cvCalcOpticalFlowBM 为重叠块 block_size.width×block_size.height 中的每一个象素计算光流,因此其速度域小于整个图像的速度域。对每一个在图像 prev 中的块,函数试图在 curr 中某些原始块或其偏移 (velx(x0,y0),vely(x0,y0)) 块的邻域里寻找类似的块,如同在前一个函数调用中所计算的类似(如果 use_previous=1)

CalcOpticalFlowPyrLK

计算一个稀疏特征集的光流,使用金字塔中的迭代 Lucas-Kanade 方法

void cvCalcOpticalFlowPyrLK( const CvArr* prev, const CvArr* curr, CvArr* prev_pyr, CvArr* curr_pyr,
                             const CvPoint2D32f* prev_features, CvPoint2D32f* curr_features,
                             int count, CvSize win_size, int level, char* status,
                             float* track_error, CvTermCriteria criteria, int flags );
prev
在时间 t 的第一帧
curr
在时间 t + dt 的第二帧
prev_pyr
第一帧的金字塔缓存. 如果指针非 NULL , 则缓存必须有足够的空间来存储金字塔从层 1 到层 #level 的内容。尺寸 (image_width+8)*image_height/3 比特足够了
curr_pyr
与 prev_pyr 类似, 用于第二帧
prev_features
需要发现光流的点集
curr_features
包含新计算出来的位置的 点集
count
特征点的数目
win_size
每个金字塔层的搜索窗口尺寸
level
最大的金字塔层数。如果为 0 , 不使用金字塔 (即金字塔为单层), 如果为 1 , 使用两层,下面依次类推。
status
数组。如果对应特征的光流被发现,数组中的每一个元素都被设置为 1, 否则设置为 0。
error
双精度数组,包含原始图像碎片与移动点之间的差。为可选参数,可以是 NULL .
criteria
准则,指定在每个金字塔层,为某点寻找光流的迭代过程的终止条件。
flags
其它选项:
  • CV_LKFLOW_PYR_A_READY , 在调用之前,第一帧的金字塔已经准备好
  • CV_LKFLOW_PYR_B_READY , 在调用之前,第二帧的金字塔已经准备好
  • CV_LKFLOW_INITIAL_GUESSES , 在调用之前,数组 B 包含特征的初始坐标 (Hunnish: 在本节中没有出现数组 B,不知是指的哪一个)

函数 cvCalcOpticalFlowPyrLK 实现了金字塔中 Lucas-Kanade 光流计算的稀疏迭代版本 ([Bouguet00])。 它根据给出的前一帧特征点坐标计算当前视频帧上的特征点坐标。 函数寻找具有子象素精度的坐标值。

两个参数 prev_pyr 和 curr_pyr 都遵循下列规则: 如果图像指针为 0, 函数在内部为其分配缓存空间,计算金字塔,然后再处理过后释放缓存。 否则,函数计算金字塔且存储它到缓存中,除非设置标识 CV_LKFLOW_PYR_A[B]_READY 。 图像应该足够大以便能够容纳 Gaussian 金字塔数据。 调用函数以后,金字塔被计算而且相应图像的标识可以被设置,为下一次调用准备就绪 (比如:对除了第一个图像的所有图像序列,标识 CV_LKFLOW_PYR_A_READY 被设置).

预估器

CvKalman

Kalman 滤波器状态

typedef struct CvKalman
{
    int MP;                     /* 测量向量维数 */
    int DP;                     /* 状态向量维数 */
    int CP;                     /* 控制向量维数 */

    /* 向后兼容字段 */
#if 1
    float* PosterState;         /* =state_pre->data.fl */
    float* PriorState;          /* =state_post->data.fl */
    float* DynamMatr;           /* =transition_matrix->data.fl */
    float* MeasurementMatr;     /* =measurement_matrix->data.fl */
    float* MNCovariance;        /* =measurement_noise_cov->data.fl */
    float* PNCovariance;        /* =process_noise_cov->data.fl */
    float* KalmGainMatr;        /* =gain->data.fl */
    float* PriorErrorCovariance;/* =error_cov_pre->data.fl */
    float* PosterErrorCovariance;/* =error_cov_post->data.fl */
    float* Temp1;               /* temp1->data.fl */
    float* Temp2;               /* temp2->data.fl */
#endif

    CvMat* state_pre;           /* 预测状态 (x'(k)): 
                                    x(k)=A*x(k-1)+B*u(k) */
    CvMat* state_post;          /* 矫正状态 (x(k)):
                                    x(k)=x'(k)+K(k)*(z(k)-H*x'(k)) */
    CvMat* transition_matrix;   /* 状态传递矩阵 state transition matrix (A) */
    CvMat* control_matrix;      /* 控制矩阵 control matrix (B)
                                   (如果没有控制,则不使用它)*/
    CvMat* measurement_matrix;  /* 测量矩阵 measurement matrix (H) */
    CvMat* process_noise_cov;   /* 过程噪声协方差矩阵
                                        process noise covariance matrix (Q) */
    CvMat* measurement_noise_cov; /* 测量噪声协方差矩阵
                                          measurement noise covariance matrix (R) */
    CvMat* error_cov_pre;       /* 先验误差计协方差矩阵
                                        priori error estimate covariance matrix (P'(k)):
                                     P'(k)=A*P(k-1)*At + Q)*/
    CvMat* gain;                /* Kalman 增益矩阵 gain matrix (K(k)):
                                    K(k)=P'(k)*Ht*inv(H*P'(k)*Ht+R)*/
    CvMat* error_cov_post;      /* 后验错误估计协方差矩阵
                                        posteriori error estimate covariance matrix (P(k)):
                                     P(k)=(I-K(k)*H)*P'(k) */
    CvMat* temp1;               /* 临时矩阵 temporary matrices */
    CvMat* temp2;
    CvMat* temp3;
    CvMat* temp4;
    CvMat* temp5;
}
CvKalman;

结构 CvKalman 用来保存 Kalman 滤波器状态。它由函数 cvCreateKalman 创建,由函数f cvKalmanPredict 和 cvKalmanCorrect 更新,由 cvReleaseKalman 释放. 通常该结构是为标准 Kalman 所使用的 (符号和公式都借自非常优秀的 Kalman 教程 [Welch95]):

系统运动方程:x_k=A\cdot x_{k-1}+ B\cdot u_k+w_k
系统观测方程:z_k=H\cdot x_k + v_k

其中:

xk(xk − 1) - 系统在时刻 k (k-1) 的状态向量 (state of the system at the moment k (k-1))
zk - 在时刻 k 的系统状态测量向量 (measurement of the system state at the moment k)
uk - 应用于时刻 k 的外部控制 (external control applied at the moment k)
wkvk 分别为正态分布的运动和测量噪声
p(w) ~ N(0,Q)
p(v) ~ N(0,R),
即,
Q - 运动噪声的相关矩阵,常量或变量
R - 测量噪声的相关矩阵,常量或变量

对标准 Kalman 滤波器,所有矩阵: A, B, H, Q 和 R 都是通过 cvCreateKalman 在分配结构 CvKalman 时初始化一次。但是,同样的结构和函数,通过在当前系统状态邻域中线性化扩展 Kalman 滤波器方程,可以用来模拟扩展 Kalman 滤波器,在这种情况下, A, B, H (也许还有 Q 和 R) 在每一步中都被更新。

CreateKalman

分配 Kalman 滤波器结构

CvKalman* cvCreateKalman( int dynam_params, int measure_params, int control_params=0 );
dynam_params
状态向量维数
measure_params
测量向量维数
control_params
控制向量维数

函数 cvCreateKalman 分配 CvKalman 以及它的所有矩阵和初始参数

ReleaseKalman

释放 Kalman 滤波器结构

void cvReleaseKalman( CvKalman** kalman );
kalman
指向 Kalman 滤波器结构的双指针

函数 cvReleaseKalman 释放结构 CvKalman 和里面所有矩阵

KalmanPredict

估计后来的模型状态

const CvMat* cvKalmanPredict( CvKalman* kalman, const CvMat* control=NULL );
#define cvKalmanUpdateByTime cvKalmanPredict
kalman
Kalman 滤波器状态
control
控制向量 (uk), 如果没有外部控制 (control_params=0) 应该为 NULL

函数 cvKalmanPredict 根据当前状态估计后来的随机模型状态,并存储于 kalman->state_pre:

x'_k=A \cdot x_{k-1}+ B \cdot u_k
P'_k=A \cdot P_{k-1} \cdot A^T + Q,

其中

x'k 是预测状态 (kalman->state_pre),
xk − 1 是前一步的矫正状态 (kalman->state_post),应该在开始的某个地方初始化,即缺省为零向量,
uk 是外部控制(control 参数),
P'k 是先验误差相关矩阵 (kalman->error_cov_pre)
Pk − 1 是前一步的后验误差相关矩阵(kalman->error_cov_post),应该在开始的某个地方初始化,即缺省为单位矩阵.

函数返回估计得到的状态值

KalmanCorrect

调节模型状态

const CvMat* cvKalmanCorrect( CvKalman* kalman, const CvMat* measurement );
#define cvKalmanUpdateByMeasurement cvKalmanCorrect
kalman
被更新的 Kalman 结构的指针
measurement
指向测量向量的指针,向量形式为 CvMat

函数 cvKalmanCorrect 在给定的模型状态的测量基础上,调节随机模型状态:

K_k=P'_k\cdot H^T \cdot (H \cdot P'_k \cdot H^T+R)^{-1}
x_k=x'_k+K_k \cdot (z_k-H \cdot x'_k)
P_k=(I-K_k \cdot H)\cdot P'_k

其中

zk - 给定测量(mesurement parameter)
Kk - Kalman "增益" 矩阵

函数存储调节状态到 kalman->state_post 中并且输出时返回它。

例子. 使用 Kalman 滤波器跟踪一个旋转的点

#include "cv.h"
#include "highgui.h"
#include <math.h>

int main(int argc, char** argv)
{
    /* A matrix data */
    const float A[] = { 1, 1, 0, 1 };

    IplImage* img = cvCreateImage( cvSize(500,500), 8, 3 );
    CvKalman* kalman = cvCreateKalman( 2, 1, 0 );
    /* state is (phi, delta_phi) - angle and angle increment */
    CvMat* state = cvCreateMat( 2, 1, CV_32FC1 );
    CvMat* process_noise = cvCreateMat( 2, 1, CV_32FC1 );
    /* only phi (angle) is measured */
    CvMat* measurement = cvCreateMat( 1, 1, CV_32FC1 );
    CvRandState rng;
    int code = -1;

    cvRandInit( &rng, 0, 1, -1, CV_RAND_UNI );

    cvZero( measurement );
    cvNamedWindow( "Kalman", 1 );

    for(;;)
    {
        cvRandSetRange( &rng, 0, 0.1, 0 );
        rng.disttype = CV_RAND_NORMAL;

        cvRand( &rng, state );

        memcpy( kalman->transition_matrix->data.fl, A, sizeof(A));
        cvSetIdentity( kalman->measurement_matrix, cvRealScalar(1) );//初始化带尺度的单位矩阵 
        cvSetIdentity( kalman->process_noise_cov, cvRealScalar(1e-5) );
        cvSetIdentity( kalman->measurement_noise_cov, cvRealScalar(1e-1) );
        cvSetIdentity( kalman->error_cov_post, cvRealScalar(1));
        /* choose random initial state */
        cvRand( &rng, kalman->state_post );

        rng.disttype = CV_RAND_NORMAL;

        for(;;)
        {
            #define calc_point(angle)                                      \
                cvPoint( cvRound(img->width/2 + img->width/3*cos(angle)),  \
                         cvRound(img->height/2 - img->width/3*sin(angle)))

            float state_angle = state->data.fl[0];
            CvPoint state_pt = calc_point(state_angle);

            /* predict point position */
            const CvMat* prediction = cvKalmanPredict( kalman, 0 );
            float predict_angle = prediction->data.fl[0];
            CvPoint predict_pt = calc_point(predict_angle);
            float measurement_angle;
            CvPoint measurement_pt;

            cvRandSetRange( &rng, 0, sqrt(kalman->measurement_noise_cov->data.fl[0]), 0 );
            cvRand( &rng, measurement );

            /* generate measurement */
            cvMatMulAdd( kalman->measurement_matrix, state, measurement, measurement );

            measurement_angle = measurement->data.fl[0];
            measurement_pt = calc_point(measurement_angle);

            /* plot points */
            #define draw_cross( center, color, d )                                 \
                cvLine( img, cvPoint( center.x - d, center.y - d ),                \
                             cvPoint( center.x + d, center.y + d ), color, 1, 0 ); \
                cvLine( img, cvPoint( center.x + d, center.y - d ),                \
                             cvPoint( center.x - d, center.y + d ), color, 1, 0 )

            cvZero( img );
            draw_cross( state_pt, CV_RGB(255,255,255), 3 );
            draw_cross( measurement_pt, CV_RGB(255,0,0), 3 );
            draw_cross( predict_pt, CV_RGB(0,255,0), 3 );
            cvLine( img, state_pt, predict_pt, CV_RGB(255,255,0), 3, 0 );

            /* adjust Kalman filter state */
            cvKalmanCorrect( kalman, measurement );

            cvRandSetRange( &rng, 0, sqrt(kalman->process_noise_cov->data.fl[0]), 0 );
            cvRand( &rng, process_noise );
            cvMatMulAdd( kalman->transition_matrix, state, process_noise, state );

            cvShowImage( "Kalman", img );
            code = cvWaitKey( 100 );

            if( code > 0 ) /* break current simulation by pressing a key */
                break;
        }
        if( code == 27 ) /* exit by ESCAPE */
            break;
    }

    return 0;
}

CvConDensation

ConDensaation 状态

typedef struct CvConDensation
{
    int MP;     // 测量向量的维数: Dimension of measurement vector
    int DP;     // 状态向量的维数: Dimension of state vector
    float* DynamMatr;       // 线性动态系统矩阵:Matrix of the linear Dynamics system
    float* State;           // 状态向量: Vector of State
    int SamplesNum;         // 粒子数: Number of the Samples
    float** flSamples;      // 粒子向量数组: array of the Sample Vectors
    float** flNewSamples;   // 粒子向量临时数组: temporary array of the Sample Vectors
    float* flConfidence;    // 每个粒子的置信度(译者注:也就是粒子的权值):Confidence for each Sample
    float* flCumulative;    // 权值的累计: Cumulative confidence
    float* Temp;            // 临时向量:Temporary vector
    float* RandomSample;    // 用来更新粒子集的随机向量: RandomVector to update sample set
    CvRandState* RandS;     // 产生随机向量的结构数组: Array of structures to generate random vectors
} CvConDensation;

结构 CvConDensation中条件概率密度传播(译者注:粒子滤波的一种特例)(Con-Dens-Ation: 单词 CONditional DENSity propagATION 的缩写)跟踪器的状态。该算法描述可参考 http://www.dai.ed.ac.uk/CVonline/LOCAL_COPIES/ISARD1/condensation.html

CreateConDensation

分配 ConDensation 滤波器结构

CvConDensation* cvCreateConDensation( int dynam_params, int measure_params, int sample_count );
dynam_params
状态向量的维数
measure_params
测量向量的维数
sample_count
粒子数

函数 cvCreateConDensation 创建结构 CvConDensation 并且返回结构指针。

ReleaseConDensation

释放 ConDensation 滤波器结构

void cvReleaseConDensation( CvConDensation** condens );
condens
要释放结构的双指针

函数 cvReleaseConDensation 释放结构 CvConDensation (见cvConDensation) 并且清空所有事先被开辟的内存空间。

ConDensInitSampleSet

初始化 ConDensation 算法中的粒子集

void cvConDensInitSampleSet( CvConDensation* condens, CvMat* lower_bound, CvMat* upper_bound );
condens
需要初始化的结构指针
lower_bound
每一维的下界向量
upper_bound
每一维的上界向量

函数 cvConDensInitSampleSet 在指定区间内填充结构 CvConDensation 中的样例数组。

ConDensUpdateByTime

估计下个模型状态

void cvConDensUpdateByTime( CvConDensation* condens );
condens
要更新的结构指针

函数 cvConDensUpdateByTime 从当前状态估计下一个随机模型状态。

 

posted on 2013-08-19 11:15  xingma0910  阅读(3842)  评论(0编辑  收藏  举报