Fast算法原理：fastica算法步骤详解

1. Fast算法原理

我们前面已经介绍过几个特征检测器，它们的效果都很好，特别是SIFT和SURF算法，但是从实时处理的角度来看，效率还是太低了。为了解决这个问题，Edward Rosten和Tom Drummond在2006年提出了FAST算法，并在2010年对其进行了修正。

FAST (全称Features from accelerated segment test)是一种用于角点检测的算法，该算法的原理是取图像中检测点，以该点为圆心的周围邻域内像素点判断检测点是否为角点，通俗的讲就是若一个像素周围有一定数量的像素与该点像素值不同，则认为其为角点。

1. 1 FAST算法的基本流程

在图像中选取一个像素点 p，来判断它是不是关键点。$$I_p$$等于像素点 p的灰度值。

以r为半径画圆，覆盖p点周围的M个像素，通常情狂下，设置 r=3，则 M=16，如下图所示：

设置一个阈值t，如果在这 16 个像素点中存在 n 个连续像素点的灰度值都高于$$I_p + t$$，或者低于$$I_p - t$$，那么像素点 p 就被认为是一个角点。如上图中的虚线所示，n 一般取值为 12。

由于在检测特征点时是需要对图像中所有的像素点进行检测，然而图像中的绝大多数点都不是特征点，如果对每个像素点都进行上述的检测过程，那显然会浪费许多时间，因此采用一种进行非特征点判别的方法：首先对候选点的周围每个 90 度的点：1，9，5，13 进行测试(先测试 1 和 19, 如果它们符合阈值要求再测试 5 和 13)。如果 p 是角点，那么这四个点中至少有 3 个要符合阈值要求，否则直接剔除。对保留下来的点再继续进行测试(是否有 12 的点符合阈值要求)。

虽然这个检测器的效率很高，但它有以下几条缺点：

·获得的候选点比较多

·特征点的选取不是最优的，因为它的效果取决与要解决的问题和角点的分布情况。

·进行非特征点判别时大量的点被丢弃

·检测到的很多特征点都是相邻的

前 3 个问题可以通过机器学习的方法解决，最后一个问题可以使用非最大值抑制的方法解决。

1. 2 机器学习的角点检测器

选择一组训练图片(最好是跟最后应用相关的图片)

使用 FAST 算法找出每幅图像的特征点，对图像中的每一个特征点，将其周围的 16 个像素存储构成一个向量P。

每一个特征点的 16 像素点都属于下列三类中的一种

根据这些像素点的分类，特征向量 P 也被分为 3 个子集：Pd ，Ps ，Pb，

定义一个新的布尔变量$$K_p$$，如果 p 是角点就设置为 Ture，如果不是就设置为 False。

利用特征值向量p，目标值是$K_p$，训练ID3 树(决策树分类器)。

将构建好的决策树运用于其他图像的快速的检测。

1. 3 非极大值抑制

在筛选出来的候选角点中有很多是紧挨在一起的，需要通过非极大值抑制来消除这种影响。

为所有的候选角点都确定一个打分函数$$V $$ ， $$V $$的值可这样计算：先分别计算$$I_p$$与圆上16个点的像素值差值，取绝对值，再将这16个绝对值相加，就得到了$$V $$的值

最后比较毗邻候选角点的 V 值，把V值较小的候选角点pass掉。

FAST算法的思想与我们对角点的直观认识非常接近，化繁为简。FAST算法比其它角点的检测算法快，但是在噪声较高时不够稳定，这需要设置合适的阈值。

2.Fast实现

OpenCV中的FAST检测算法是用传统方法实现的，

1.实例化fast

参数：

·threshold：阈值t，有默认值10

·nonmaxSuppression：是否进行非极大值抑制，默认值True

返回：

Fast：创建的FastFeatureDetector对象

2.利用fast.detect检测关键点，没有对应的关键点描述

参数：

gray: 进行关键点检测的图像，注意是灰度图像

返回：

kp: 关键点信息，包括位置，尺度，方向信息

3.将关键点检测结果绘制在图像上，与在sift中是一样的

示例：

结果：

% function [Ahat2, shat, n_iteration Test] = nc_fastica_svd(xold,typeStr,N,A)

function [shat Ahat2] = nc_fastica_svd(xold,typeStr,N)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% non-circular complex FastICA算法，基于Newton迭代法，类似与fastICA

% ************************input***************************

% xold：混合信号，m*n,m为阵元数，n为快拍数

% typeStr： 非线性函数，'log', 'kurt', or 'sqrt'

% **************************output**************************

% Ahat：解混矩阵

% shat: 估计的源信号

% ********************************************************

% Reference

% Mike Novey and T. Adali, "On Extending the complex FastICA algorithm

% to noncircular sources" in

% (To appear 2007/2008) IEEE Journel on Signal Processing.,

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

type = 0

if strcmp(typeStr,'log') == 1

type = 1

elseif strcmp(typeStr,'kurt') == 1

type = 2

elseif strcmp(typeStr,'sqrt') == 1

type = 3

end

tol = 1e-5

a2 = 0.1

defl = 1% components are estimated one by one in a deflationary mannerset this to 0 if you want them all estimated simultaneously

maxcounter = 50

[n,m] = size(xold)

% Whitening of s:

yyy = zeros(1,m)

[Ex, Dx] = svd(cov(xold'))

E = Ex(:,1:N)

D = Dx(1:N,1:N)

Q = mtimes(sqrt(inv(D)),E')

x = Q * xold

%Pseudo-covariance

pC = (x*transpose(x))/m

% FIXED POINT ALGORITHM

W = eye(N)

Wold = zeros(N)

k = 0

while (norm(abs(Wold'*W)-eye(N),'fro')>(N*1e-12) &&k <15*N)

k = k+1

Wold = W

for kk=1:N %Loop thru sources

yy = W(:,kk)'*x

absy =abs(yy).^2

%%Fixed point

if type == 1 %%log

g = 1./(a2 + absy)

gp = -1./(a2 + absy).^2

elseif type == 2 %Kurt

g = absy

gp = ones(size(absy))

elseif type == 3 %sqrt

g = 1./(2*sqrt(a2 + absy))

gp = -1./(4*(a2 + absy).^(3/2))

end

gRad = mean(ones(N,1)*(g.*conj(yy)).*x,2)

ggg = mean(gp.*absy + g)

B = mean(gp .* conj(yy).^2)*pC

W(:,kk) = Wold(:,kk)*(ggg) -(gRad) + (B*conj(Wold(:,kk)))

end

%Orthonormalization

[E,D] = eig(W'*W)

W = W * E * inv(sqrt(D)) * E'

end%Loop thru sources

n_iteration = k

shat = W'*x%Estimated sources

% Ahat1 = inv(Q)*W

Ahat2 = W'*Q

这个是NC-fastica，可以用。稍微注释了些

原始程序，不知道是谁写的了

1、二维点云上的FastICA

import numpy as np

import matplotlib.pyplot as plt

from sklearn.decomposition import PCA, FastICA

plt.rcParams['font.sans-serif'] = ['SimHei']

plt.rcParams['axes.unicode_minus'] = False

# 生成样本数据

rng = np.random.RandomState(42)

S = rng.standard_t(1.5, size=(20000, 2))

S[:, 0] *= 2.0

A = np.array([[1, 1], [0, 2]])  # 混合矩阵

X = np.dot(S, A.T)  # 生成观测结果

pca = PCA()

S_pca_ = pca.fit(X).transform(X)

ica = FastICA(random_state=rng)

S_ica_ = ica.fit(X).transform(X)  # 估计来源

S_ica_ /= S_ica_.std(axis=0)

# 描绘结果

def plot_samples(S, axis_list=None):

    plt.scatter(

        S[:, 0], S[:, 1], s=2, marker="o", zorder=10, color="steelblue", alpha=0.5

    )

    if axis_list is not None:

        colors = ["orange", "red"]

        for color, axis in zip(colors, axis_list):

            axis /= axis.std()

            x_axis, y_axis = axis

            plt.plot(0.1 * x_axis, 0.1 * y_axis, linewidth=2, color=color)

            plt.quiver(

                (0, 0),

                (0, 0),

                x_axis,

                y_axis,

                zorder=11,

                width=0.01,

                scale=6,

                color=color,

            )

    plt.hlines(0, -3, 3)

    plt.vlines(0, -3, 3)

    plt.xlim(-3, 3)

    plt.ylim(-3, 3)

    plt.xlabel("x")

    plt.ylabel("y")

plt.figure()

plt.subplot(2, 2, 1)

plot_samples(S / S.std())

plt.title("真实独立源")

axis_list = [pca.components_.T, ica.mixing_]

plt.subplot(2, 2, 2)

plot_samples(X / np.std(X), axis_list=axis_list)

legend = plt.legend(["PCA", "ICA"], loc="upper right")

legend.set_zorder(100)

plt.title("观察")

plt.subplot(2, 2, 3)

plot_samples(S_pca_ / np.std(S_pca_, axis=0))

plt.title("PCA恢复信号")

plt.subplot(2, 2, 4)

plot_samples(S_ica_ / np.std(S_ica_))

plt.title("ICA恢复信号")

plt.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.36)

plt.show()

---