njustyxy

搞KLPP时用到归一化这个概念，转一个帖子：

归一化化定义：我是这样认为的，归一化化就是要把你需要处理的数据经过处理后（通过某种算法）限制在你需要的一定范围内。首先归一化是为了后面数据处理的方便，其次是保正程序运行时收敛加快。
在 matlab 里面，用于归一化的方法共有三种 :
（ 1 ） premnmx 、 postmnmx 、 tramnmx
（ 2 ） prestd 、 poststd 、 trastd
（ 3 ）是用 matlab 语言自己编程。
premnmx 指的是归一到 [ － 1 1],prestd 归一到单位方差和零均值。（ 3 ）关于自己编程一般是归一到 [0.1  0.9] 。具体用法见下面实例。
为什么要用归一化呢？首先先说一个概念，叫做奇异样本数据，所谓奇异样本数据数据指的是相对于其他输入样本特别大或特别小的样本矢量。
下面举例：
m=[0.11 0.15 0.32 0.45 30;
      0.13 0.24 0.27 0.25 45];
其中的第五列数据相对于其他 4 列数据就可以成为奇异样本数据（下面所说的网络均值 bp ）。奇异样本数据存在所引起的网络训练时间增加，并可能引起网络无法收敛，所以对于训练样本存在奇异样本数据的数据集在训练之前，最好先进形归一化，若不存在奇异样本数据，则不需要事先归一化。
具体举例：
close all
clear
echo on
clc
%BP 建模
% 原始数据归一化
m_data=[1047.92 1047.83 0.39 0.39 1.0 3500 5075;
    1047.83 1047.68 0.39 0.40  1.0 3452 4912;
    1047.68 1047.52  0.40  0.41 1.0  3404 4749;
    1047.52  1047.27  0.41  0.42 1.0  3356 4586;
    1047.27  1047.41 0.42 0.43  1.0  3308  4423;
    1046.73  1046.74 1.70 1.80 0.75  2733  2465;
    1046.74  1046.82  1.80  1.78 0.75  2419 2185;
    1046.82 1046.73  1.78  1.75  0.75 2105  1905;
    1046.73  1046.48 1.75 1.85 0.70 1791  1625;
    1046.48  1046.03  1.85  1.82  0.70 1477 1345;
    1046.03 1045.33 1.82 1.68  0.70  1163  1065;
    1045.33  1044.95  1.68  1.71 0.70  849  785;
    1044.95  1045.21 1.71  1.72  0.70  533  508;
    1045.21 1045.64  1.72  1.70 0.70 567  526;
    1045.64 1045.44 1.70  1.69  0.70  601  544;
    1045.44 1045.78  1.69  1.69 0.70  635  562;
    1045.78 1046.20  1.69  1.52 0.75  667  580];
% 定义网络输入 p 和期望输出 t
pause
clc
p1=m_data(:,1:5);
t1=m_data(:,6:7);
p=p1';t=t1';
[pn,minp,maxp,tn,mint,maxt]=premnmx(p,t)
% 设置网络隐单元的神经元数 (5~30 验证后 5 个最好）
n=5;
% 建立相应的 BP 网络
pause
clc
net=newff(minmax(pn),[n,2],{'tansig','purelin'},'traingdm');
inputWeights=net.IW{1,1};
inputbias=net.b{1};
layerWeights=net.IW{1,1};
layerbias=net.b{2};
pause
clc
% 训练网络
net.trainParam.show=50;
net.trainParam.lr=0.05;
net.trainParam.mc=0.9;
net.trainParam.epochs=200000;
net.trainParam.goal=1e-3;
pause
clc
% 调用 TRAINGDM 算法训练 BP 网络
net=train(net,pn,tn);
% 对 BP 网络进行仿真
A=sim(net,pn);
E=A-tn;
M=sse(E)
N=mse(E)
pause
clc
p2=[1046.20 1046.05 1.52 1.538 0.75;
    1046.05 1046.85 1.538 1.510 0.75;
    1046.85 1046.60 1.510 1.408 0.75;
    1046.60 1046.77 1.408 1.403 0.75;
    1046.77 1047.18 1.403 1.319 0.75];
p2=p2';
p2n=tramnmx(p2,minp,maxp);
a2n=sim(net,p2n);
a2=postmnmx(a2n,mint,maxt)
echo off
pause
clc
程序说明：所用样本数据（见 m_data ）包括输入和输出数据，都先进行归一化，还有一个问题就是你要进行预测的样本数据 ( 见本例 p2) 在进行仿真前，必须要用 tramnmx 函数 进行事先归一化处理，然后才能用于预测，最后的仿真结果要用 postmnmx 进行反归一，这时的输出数据才是您所需要的预测结果。
个人认为： tansig 、 purelin 、 logsig 是网络结构的传递函数，本身和归一化没什么直接关系，归一化只是一种数据预处理方法。

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/paincupid/archive/2009/06/16/4273970.aspx

需要说明的事并不是任何问题都必须事先把原始数据进行规范化,也就是数据规范化这一步并不是必须要做的,要具体问题具体看待,测试表明有时候规范化后的预测准确率比没有规范化的预测准确率低很多.就最大最小值法而言,当你用这种方式将原始数据规范化后,事实上意味着你承认了一个假设就是测试数据集的每一模式的所有特征分量的最大值(最小值)不会大于(小于)训练数据集的每一模式的所有特征分量的最大值(最小值),但这条假设显然过于强,实际情况并不一定会这样.使用平均数方差法也会有同样类似的问题.故数据规范化这一步并不是必须要做的,要具体问题具体看待

由于采集的各数据单位不一致，因而须对数据进行[-1，1]归一化处理，归一化方法主要有如下几种，供大家参考：（by james）
1、线性函数转换，表达式如下：
y=(x-MinValue)/(MaxValue-MinValue)
说明：x、y分别为转换前、后的值，MaxValue、MinValue分别为样本的最大值和最小值。
2、对数函数转换，表达式如下：
y=log10(x)
说明：以10为底的对数函数转换。
3、反余切函数转换，表达式如下：
y=atan(x)*2/PI
归一化是为了加快训练网络的收敛性，可以不进行归一化处理

归一化的具体作用是归纳统一样本的统计分布性。归一化在0-1之间是统计的概率分布，归一化在-1--+1之间是统计的坐标分布。归一化有同一、统一和合一的意思。无论是为了建模还是为了计算，首先基本度量单位要同一，神经网络是以样本在事件中的统计分别几率来进行训练（概率计算）和预测的，归一化是同一在0-1之间的统计概率分布；

当所有样本的输入信号都为正值时，与第一隐含层神经元相连的权值只能同时增加或减小，从而导致学习速度很慢。为了避免出现这种情况，加快网络学习速度，可以对输入信号进行归一化，使得所有样本的输入信号其均值接近于0或与其均方差相比很小。

归一化是因为sigmoid函数的取值是0到1之间的，网络最后一个节点的输出也是如此，所以经常要对样本的输出归一化处理。所以这样做分类的问题时用[0.9 0.1 0.1]就要比用[1 0 0]要好。

但是归一化处理并不总是合适的，根据输出值的分布情况，标准化等其它统计变换方法有时可能更好。
关于用premnmx语句进行归一化：
premnmx语句的语法格式是：[Pn,minp,maxp,Tn,mint,maxt]=premnmx(P,T)
其中P，T分别为原始输入和输出数据，minp和maxp分别为P中的最小值和最大值。mint和maxt分别为T的最小值和最大值。
premnmx函数用于将网络的输入数据或输出数据进行归一化，归一化后的数据将分布在[-1,1]区间内。
我们在训练网络时如果所用的是经过归一化的样本数据，那么以后使用网络时所用的新数据也应该和样本数据接受相同的预处理，这就要用到tramnmx。
下面介绍tramnmx函数：
[Pn]=tramnmx(P,minp,maxp)
其中P和Pn分别为变换前、后的输入数据，maxp和minp分别为premnmx函数找到的最大值和最小值。
（by terry2008）

matlab中的归一化处理有三种方法
1. premnmx、postmnmx、tramnmx
2. restd、poststd、trastd
3. 自己编程
具体用那种方法就和你的具体问题有关了
（by happy）
pm=max(abs(p(i,:))); p(i,:)=p(i,:)/pm;
和
for i=1:27
p(i,:)=(p(i,:)-min(p(i,:)))/(max(p(i,:))-min(p(i,:)));
end 可以归一到0 1 之间
0.1+(x-min)/(max-min)*(0.9-0.1)其中max和min分别表示样本最大值和最小值。
这个可以归一到0.1-0.9

=================================by  ratbaby
补充一个吧， 归一还可以用 mapminmax。
这个函数可以把矩阵的每一行归一到[-1 1].
[y1,PS] = mapminmax(x1). 其中x1 是需要归一的矩阵 y1是结果
当需要对另外一组数据做归一时，比如SVM 中的 training data用以上方法归一，而test data就可以用下面的方法做相同的归一了
y2 = mapminmax('apply',x2,PS)
当需要把归一的数据还原时，可以用以下命令
x1_again = mapminmax('reverse',y1,PS)

PS:归一化是指将数据正态化，也就是将原始数据转化为符合平均值为0，标准差为1的正态分布的数据。

将原始数据减去平均值，然后除以标准差就可以了。
另外一种方法，是将求出的参数，如能量谱，除以方差就可以了。

 也可参考《应用多元统计分析方法》P27，《基于KPCA法的定风量空调系统传感器故障诊断》

Abstract : The goals of this lecture is to use the level set framework in order to do curve evolution. The mean curvature motion is the basic tool, and it can be extended into edge-based (geodesic active contours) and region-based (Chan-Vese) snakes.

Setting up Matlab.

• First download the Matlab toolbox toolbox_fast_marching.zip. Unzip it into your working directory. You should have a directory toolbox_fast_marching/ in your path.
• The first thing to do is to install this toolbox in your path.

  path(path, 'toolbox_fast_marching/');
  path(path, 'toolbox_fast_marching/data/');
  path(path, 'toolbox_fast_marching/toolbox/');

• Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

  cd toolbox_fast_marching
  compile_mex;
  cd ..

Managing level set function.

• In order to perform curve evolution, we will deal with a distance function stored in a 2D image D. The curve will be embedded in the level set D=0. This curve will be evolved by modifying the image D. A curve evolution ODE can be replaced by an PDE on D. This allows to deal with topological changes when a curve split or two curves merge.

  n = 200; % size of the image
  % load a distance function
  D0 = compute_levelset_shape('circlerect2', n);
  % type 'help compute_levelset_shape' to see other
  % basic curve you can load.
  
  % display the curve
  clf; hold on;
  imagesc(D); axis image; axis off; axis([1 n 1 n]);
  [c,h] = contour(D,[0 0], 'r');
  set(h, 'LineWidth', 2);
  hold off;
  colormap gray(256);
  
  % do the union of two curves
  options.center = [0.15 0.15]*n;
  options.radius = 0.1*n;
  D1 = compute_levelset_shape('circle', n,options);
  imagesc(min(D0,D1)<0);

• During the curve evolution, the image D might become far from being a distance function. In order to stabilize the algorithm, one needs to re-compute this distance function.

  % here we simulate a modification of the distance function
  [Y,X] = meshgrid(1:n,1:n);
  D = (D0.^3) .* (X+n/3);
  D1 = perform_redistancing(D);
  % display both the original and the new,
  % redistanced, curve (should be very close)
  ...

Mean Curvature Motion.

• In order to compute differential quantities (tangent, normal, curvature, etc) on the curve, you can compute derivatives of the image D.

  % the gradient
  g0 = divgrad(D);
  % display the gradient (as arrow field with 'quiver', ...)
  ...
  % the normalized gradient
  d = max(eps, sqrt(sum(g0.^2,3)) );
  g = g0 ./ repmat( d, [1 1 2] );
  % display
  ...
  % the curvature
  K = d .* divgrad( g );
  % display
  ...

• The mean curvature motion of the level sets of some image u is driven be the following equation.

Implement this evolution explicitly in time using finite differences.

Tmax = 1000; % maximum time of evolution
dt = 0.4; % time step (should be small)
niter = round(Tmax/dt); % number of iterations
D = D0; % initialization
for i=1:niter
    % compute the right hand size of the PDE
    ...
    % update the distance field
    D = ...;
    % redistance the function from time to time
    if mod(i,30)==0
        D = perform_redistancing(D);
    end
    % display from time to time
    if mod(i,30)==1
        % display here
        ...
    end
end

 

Curve evolution under the mean curvature motion (the background is the distance function D).

Edge-based Segmentation with Geodesic Active Contour (snakes + level set).

• Given a background image M to segment, one needs to compute an edge-stopping function E. It should be small in area of high gradient, and high in area of large gradient.

  % load an image
  name = 'brain';
  M = rescale( sum( load_image(name, n), 3) );
  % display it
  ...
  % compute a smoothed gradient
  sigma = 4; % blurring size
  G = divgrad( perform_blurring(M,sigma) );
  % compute the norm of the gradient
  d = ...
  % compute the edge-stopping function
  E = 1./max(d, 1e-3);
  % rescale it so that it is in realistic ranges
  E = rescale(E,0.3,1);

• The geodesic active contour evolution is a mean curvature motion modulated by the edge-stopping function:
              
  Implement this evolution explicitly in time using finite differences [warning: the g in the formula is the edge stopping function E].

  % initialize the level set
  D0 = compute_levelset_shape('square', n);
  % precompute the gradient of E [which is g in the formula]
  ...
  % do the iterations
  ...

Region-based Segmentation with Chan-Vese (Mumord-Shah + level sets).

• The geodesic active contour uses an edge-based energy. It has lots of local minima and is very sensitive to initialization. In order to circumvent these drawbacks, one can use a region based energy like the Mumford-Shah functional. Re-casted into the level set framework, it reads:

The corresponding gradient descent is the Chan-Vese active contour method:

Implement this evolution explicitly in time using finite differences, when c1 and c2 are known in advanced.

% initialize with a complex distance function
D0 = compute_levelset_shape('small-disks', n);
% set parameters
lambda = 0.8;
c1 = ...; % black
c2 = ...; % gray
...

Segmentation with Chan-Vese active contour without edges.

没办法，常常国外网站不能访问，转一下。

Abstract : Learn the basic features of Matlab, and learn how to load and visualize signals and images.

Basic Matlab commands.

• Create a variable and an array

  % this is a comment
  a = 1; a = 2+1i; % real and complex numbers
  b = [1 2 3 4]; % row vector
  c = [1; 2; 3; 4]; % column vector
  d = 1:2:7; % here one has d=[1 3 5 7]
  A = eye(4); B = ones(4); C = rand(4); % identity, 1 and random matrices
  c = b'; % transpose

• Modification of vectors and matrices

  A(2,2) = B(1,1) + b(1); % to access an entry in a vector or matrix
  b(1:3) = 0; % to access a set of indices in a matrix
  b(end-2:end) = 1: % to access the last entries
  b = b(end:-1:1); % reverse a vector
  b = sort(b); % sort values
  b = b .* (b>2); % set to zeros (threshold) the values below 2
  b(3) = []; % suppress the 3rd entry of a vector
  B = [b; b]; % create a matrix of size 2x4
  c = B(:,2); % to access 2nd column

• Advanced instructions

  a = cos(b); a = sqrt(b); % usual function
  help perform_wavelet_transform; % print the help
  a = abs(b); a = real(b); a = imag(b); a = angle(b); % norm, real part and imaginary part of a complex
  disp('Hello'); % display a text
  disp( sprintf('Value of x=%.2f', x) ); % print a values with 2 digits
  A(A==Inf) = 3; % replace Inf values by 3
  A(:); % flatten a matrix into a column vector
  max(A(:)); % max of a matrix
  M = M .* (abs(M)>T); % threshold to 0 values below T.

• Display

  plot( 1:10, (1:10).^2 ); % display a 1D function
  title('My title'); % title
  xlabel('variable x'); ylabel('variable y'); % axis
  subplot(2, 2, 1); % divide the screen in 2x2 and select 1st quadrant

• Programmation

  for i=1:4 % repeat the loop for i=1, i=2, i=3 et i=4
     disp(i); % make here something
  end
  i = 4;
  while i>0 % while syntax
     disp(i); % do smth
     i = i-1;
  end

Load and visualize signals and images

• Load and plot a signal: (function load_signal.m should be in the toolbox of each course)

  f = load_signal('Piece-Regular', n); % signal of size n
  plot(f);

• Load and display an image (download function load_image.m should be in the toolboxes)

  I = load_image('lena');
  imagesc(I); axis image; colormap gray(256);

为了搞manifold，要用到matlab编译，网上找的

软件下载地址: https://sourceforge.net/projects/gnumex 软件安装和使用说明: http://gnumex.sourceforge.net/   软件功能: 在MATLAB上调用MinGW或Cygwin编译C和Fortran的mex代码。 这就意味着，我们可以在MATLAB上编译Linux的C程序代码和调用Linux的库编译出MATLAB能用的Mex程序。What a nice tool! 安装步骤: 1. Sourceforge上下载软件，目前最新版2.01。解压缩到任意一个固定的目录。如C:/MATLAB/gnumex 2. 运行MATLAB, 把上面的目录加入到MATLAB的path中。 3. MATLAB上输入命令gnumex运行安装设置。设置很简单，主要是一下几个方面 - MinGW的root目录 (如果使用的是MinGW的gcc，这里必须设置正确) - Cygwin的root目录 (同上，使用时设置，否则留空) - f95 和 gfortran的目录 (有的话就设置，没有的话留空。这个是编译fortran程序选用的编译器) 注意: gfortran是gcc的一个组件，这个一般都有，就在MinGW的root目录/bin下面。最好选上，有备无患。 - linking环境 这里选用mingw或cygwin根据自己的需要设置。注意 -mno-cygwin这一项是指在cygwin中link的时候使用windows的lib. 据说这个Cygwin gcc的option争议很大，功能不健全，在官方的maillist上要求取缔此项功能的言论很多，虽然目前仍然在更新和支持中。所以，建议还是使用MinGW比较好。 - 语言方便，选用自己需要的C/C++或fortran语言。 - Generate 里选择 mex dll. 这个一般是默认的。 当然mex也可以编译成 exe文件。根据需要选择。建议使用mex dll。这个可以再matlab上像一般函数一样直接调用程序，很方便。 关于最后的两个路径，建议默认。最后一个mexopts.bat的路径是MATLAB默认的mex程序设置路径，我们覆盖这个文件将使我们的设置成为mex的默认设置。这样就可以直接调用mex编译，如 mex    你的命令 而不用加参。加参使用的方法如: mex   -f    你的mex设置文件.bat   你的命令 如果你不想改变默认的mex设置，也可以把这个文件保存在其他方便的位置。需要的时候使用 -f 参数调用。 设置完成后，按make option file确认即可。 提示: 1- 如果想重新设置MATLAB的mex. 可以使用mex -setup 设置。这是matlab的默认命令和gnumex无关。 2- 如果要查看当前mex的设置情况，可以使用mex -v 来查看。 3- 你可以在gnumex的File菜单下，选择save config保存当前的设置界面信息到gnumexcfg.mat中。下次修改的时候可以使用菜单上的Load config调用。           测试运行: 先检查一下mex的设置是否生效 mex -v 我们看到       MATLAB                 = D:\MATLAB~1 ->    COMPILER               = gcc ->    Compiler flags:          COMPFLAGS           = -c -DMATLAB_MEX_FILE          OPTIMFLAGS          = -O3          DEBUGFLAGS          = -g 看到上面的gcc了吗，这说明我们的设置已经生效了，现在mex用的是gcc为编译器。 做个c程序的例子 复制MATLAB根目录下 extern\examples\mex\yprime.c 文件到自己的测试目录下。例如gnumex的examples目录下。 输入命令 mex yprime.c 就完成了编译工作，生成 yprime.mexw32 文件。 我们可以在MATLAB上调用这个程序运行看看。 yprime(1,1:4) ans =     2.0000    8.9685    4.0000   -1.0947 提醒: MATLAB对命令的优先级是，当前目录下的程序最优先。       使用技巧: 1. 如果你使用的是cygwin, 注意一定要把cygwin里的cygwin1.dll文件(在root的bin里面) 加入到windows的path里（而不是matlab的path里）。因为cygwin的程序需要调用这个动态链接库才能使用。这也是cygwin的不方便之处，而且会影响程序的运行速度。最重要的是，目前该软件对cygwin的版本支持很老(只支持cygwin的gcc 3.2)，现在的4.x都不支持。所以建议使用MinGW吧。那个没版本限制，我已经成功的和最新的gcc4.4.1兼容了。   2. 如果测试程序删不掉，如yprime.mexw32文件无法删除，说明matlab运行了这个程序而没有卸载。我们在MATLAB上输入clear yprime。然后就可以正常删除了。 3. gcc有很多的options，和注意事项。如果大家想更好学习gcc。建议参考官方的文档http://gcc.gnu.org/onlinedocs/，也可以找一份中文的参考文档慢慢研究。 gnumex的gcc使用的默认的options是 -O3 -mcpu=pentium -malign-double -fno-exceptions 需要什么option, 大家可以手动修改mexopts.bat文件。(高级用户使用，不熟悉gcc的人慎用) 4. 调用LAPACK/BLAS的方法。建议直接调用MATLAB里提供的lib静态库。位置在extern/lib/win32/microsoft/   分别是 libmwlapack.lib 和 libmwblas.lib. 使用的时候可以直接输入全地址调用，也可以用 -lmwlapack -lmwblas 作为MATLAB的option使用(放在最后)。 如果想用自己优化过的lapack和blas库，使用的时候可以直接输入全地址。也可以把他们改名，前面加lib，然后 复制到上面提到的位置下，通过-l调用。例如: 你有optlapack.lib文件，改名为liboptlapack.lib 放到extern/lib/win32/microsoft/ 下。然后mex命令中使用的时候加入 -loptlapack 就可以了。 5. 如果你要link一个fortran的obj(o)或lib(a)文件到一个c程序中，那么请注意在编译fortran的时候使用-fno-underscoring。 6. 通过autotool使用makefile调用MATLAB编译mex文件的方法和工具可以参考 http://gnumex.sourceforge.net/autotools/。这是一个非常有用的话题，等于把Windows的MATLAB移植到MinGW下使用。可以通过Makefile大批量编译文件。值得关注和研究。   kaien   26/07/2009 根据大家的回复，由于gnumex的版本到2.01就不更新了，所以不支持MATLAB2008以后的版本，因此我替作者给程序做了修正，算是版本v2.02吧。测试已经兼容了MATLAB2010。 更新文件下载: http://kaienfr.ys168.com/   替换2.01版的同名文件即可使用。 另外，对于matlab2010，使用gnumex编译前，必须用文本编辑器打开matlab目录下的bin/mex.pl文件，把两个$IMPLICIT_LIBS删掉，这样就不会出现gcc: getValidInputLinkLibraries: No such file or directory 的错误了。 kaien 01/09/2010

看到老外的一个教学，正好在学。转一下

This course is an introduction to the computational theory of manifolds. Manifold models arise in various area of mathematics, image processing, data mining or computer science. Surfaces of arbitrary dimension can be used to model non-linear datasets that one encounters in modern data processing. Numerical methods allow to exploit this geometric non-linear prior in order to extract relevant information from the data. These methods include in particular local differential computations (related to the Laplacian operator and its variants) and global distance methods (related to geodesic computations). In this course, you will learn how to perform differential and geodesic computations on images, volumes, surfaces and high dimensional graphs.

The course includes a set of Matlab experiments. These experiments give an overview of various tasks in computer vision, image processing, learning theory and mesh processing. This includes computation of shortest paths, Voronoi segmentations, geodesic Delaunay triangulations, surface flattening, dimensionality reduction and mesh processing.

One should copy/paste the provided code into a file named e.g. experiments.m, and launch the script directly from Matlab command line > experiments;. Some of the scripts contain "holes" that you should try to fill on your own.

图像处理中不适定问题（ill posed problem）或称为反问题（inverse Problem）的研究从20世纪末成为国际上的热点问题，成为现代数学家、计算机视觉和图像处理学者广为关注的研究领域。数学和物理上的反问题的研究由来已久，法国数学家阿达马早在19世纪就提出了不适定问题的概念:称一个数学物理定解问题的解存在、唯一并且稳定的则称该问题是适定的（Well Posed）.如果不满足适定性概念中的上述判据中的一条或几条，称该问题是不适定的。典型的图像处理不适定问题包括：图像去噪（Image De-nosing），图像恢复（Image Restorsion），图像放大（Image Zooming），图像修补（Image Inpainting），图像去马赛克（image Demosaicing），图像超分辨(Image super-resolution )等。
迄今为止，人们已经提出许多方法来解决图像处理中的不适定性。但是如何进一步刻画图像的边缘、纹理和角形等图像中重要视觉几何结构，提高该类方法在噪声抑制基础上有效保持结构和纹理能力是有待深入研究的问题。

1 不适定图像处理问题的国内外研究现状评述

由于图像处理中的反问题往往是不适定的。解决不适定性的有效途径是在图像处理中引入关于图像的先验信息。因此图像的先验模型对于图像反问题和其它计算机视觉还是图像处理问题至关重要。对于图像的先验模型的研究，研究者们从多个角度进行研究，其代表主要有“统计方法”和“正则化几何建模方法”，“稀疏表示方法”三种主流方法，而最近兴起的图像形态分量分析（MCA）方法吸引了大批国内外研究者的广泛关注。
1.1 正则化几何模型日新月异
关于自然图像建模的“正则化几何方法”是最近几年热点讨论的主题。其中一类方法是利用偏微分方程理论建立图像处理模型，目前的发展趋势是从有选择性非线性扩散的角度设计各类低阶、高阶或者低阶与高阶综合的偏微分方程, 或者从实扩散向复扩散推广, 从空域向空频域相结合以及不同奇异性结构的综合处理[1]。
另一类方法是基于能量泛函最优的变分方法。1992年，Rudin-Osher-Fatemi 提出图像 能被分解为一个属于有界变差空间的分量 和一个属于 的分量 的全变差模型 [2]。根据国际上及本人的研究表明：ROF模型模型较好地刻画了图像中视觉重要边缘结构，但不能描述纹理信息。2001年Meyer提出了振荡模式分解理论[2]：他认为振荡分量可以表示为某个向量函数的散度形式，而振荡分量可以属于3个可能的函数空间。首先引入有界变差(bounded variational , BV) 空间的一个近似对偶空间来表征图像的振荡分量；Meyer进一步指出John-Nirenberg的有界均值振荡空间和齐性Besov空间都是振荡分量比较合适的函数空间，由此导出了将图像分解的（BV,G）模型,（BV,F）模型和（BV,E）模型。Meyer从理论上基本解决了振荡分量的理论框架，成为纹理等振荡模式分解的奠基性工作，但是原始模型比较难计算。后来的学者大都在Meyer工作的基础上展开工作。Vese-Osher提出将振荡分量建模为的向量场的散度来逼近（BV,G）模型[3]，实质上是将G空间 近似为负Soblev空间 [4]。L.Lieu和L.Vese进一步推广到分数阶负Soblev空间 [5]。Aujol,Chamboll等人定义了G-空间中的一个子空间，并根据Chamboll早期提出的ROF模型的投影算法的基础上，提出图像的振荡分量是在该子空间上的投影分量，由此提出了著名的BV空间半范 + G空间范数 + L2 范数约束优化的A2BC模型及子空间投影算法 [6-7]。J.B.Garnet，T.M.Le,Y.Meyer, L.A.Vese提出更一般的齐性Besov空间 来刻画振荡分量 [8]。最近，J.Aujol, A.Chamboll分别对TV范数、G范数、F范数、E范数，L 2 范数对图像的卡通图像、纹理分量、高斯噪声进行数理统计和相关性分析，提出了分别运用TV范数、G范数和E范数分别来约束图像的卡通分量、纹理分量 和噪声分量 的三分量图像分解模型[9]。2007年，G.Gilboa 和S. Osher受提出了非局部化G-空间的概念，并概括性的初步提出了非局部ROF模型、非局部Meyer模型、非局部ROF+L1模型[10]，从理论上提供了图像先验模型研究的新思路。但综合目前研究来看，变分方法的主要不足是对于纹理和噪声的刻画还不够精细。

1.2 稀疏表示方兴未艾

图像的稀疏表示问题最早源于“有效编码假说”。Attneave最先提出：视觉感知的目标就是产生一个外部输入信号的有效表示。在神经生物学领域Barlow基于信息论提出了“有效编码假设”，认为初级视皮层神经细胞的主要功能就是去除输入刺激的统计相关性[11]。“有效编码假设”被提出以后，很多研究人员根据它的思想提出了不同的理论。主要思路分为两大类。直接方法是机理测试方法，即从生物机理上，在自然图像刺激条件下检测神经细胞的响应特性。著名的工作如：2001年在《Nature》上发表的研究结果表明，在冗余性测度和自然刺激条件下一组视网膜神经节对外界刺激独立编码[12]；2000年在《Science》上发表了类似的成果[13]：通过记录短尾猿 V1 区神经细胞在开放的自然场景和模拟自然场景条件下的神经细胞响应，验证了视皮层（V1 区）神经细胞用稀疏编码有效表示自然场景，稀疏编码用最小冗余度传递信息。另外一个替代的方法是模型仿真方法，即利用自然图像的统计特性，建立模型模拟早期视觉处理系统的处理机制。例如Olshausen 和 Field[14] 提出了稀疏编码模型，稀疏编码理论表明，通过寻找自然图像的稀疏编码表示，该神经网络可以学习得到类似于简单细胞感受野的结构。Bell提出了基于信息最大化的无监督算法，通过度量“因子”的联合信息熵并且使之最大化，扩展了独立成分分析（ICA）方法，成功地构建有效编码模型并得到了与上面类似的结果[15]。Hyvarinen 更进一步，应用一个两层的稀疏编码模型构造出类似于复杂细胞响应特性的基函数，而且基函数集合形成一个有规律的拓扑结构[16]。这部分表明有效编码假设也可适用于视觉系统高级区域神经细胞的处理过程。
目前关于图像稀疏表示系统的研究大体上沿着两条主线展开。其中一条是沿着多尺度几何分析理论。研究者认为图像的非平稳性和非高斯性，很难用线性算法进行处理，而应该建立合适的能够处理边缘到纹理各层面几何结构的图像模型；二维图像中的性状奇异性边缘和3-D 图像中丝状物(filaments) 和管状物(tubes)几何特征不能被各向同性的“方块基”（如小波基）表示，而最优或者 “最稀疏”的函数表示方法应该由各向异性的“锲形基”表征。因此以Ridgelet、Curvelet、Bandlet, Contourlet变换为代表的多尺度几何分析[16-22]理论成为图像稀疏表示的有效途径。图2.1.1(a)给出了二维可分离小波在不同分辨率下逼近曲线的过程，随着分辨率升高，尺度变细，最终表现为使用众多的“点”来逼近曲线。
 

与小波相比，contourlet不仅具有小波的多分辨率特性和时频局部化特性，还具有很好的方向性和各向异性，即在尺度j时，小波基的支撑域边长近似为，而Contourlet的在该尺度下的基函数支撑域的纵横比可以任意选择。图2.1.1(b)为用Contourlet基函数的支撑域来逼近曲线的过程，由于它的基函数的支撑域表现为“长方形”，因而是一种更为有效稀疏的表示法。与二维可分离小波基函数的方向支撑域的各向同性不同，Contourlet基的“长方形”支撑域表现出来的是各向异性(anisotropy)的特点。
上述稀疏表示方法都是采用“单一基”，另外一条图像稀疏表示的途径是：基函数被称之为原子库的过完备的冗余系统取代。Mallat和Zhang于1993年首先提出了信号在过完备库(over-complete dictionary)上分解的思想[23]. 通过信号在过完备库上的分解,用来表示信号的基可自适应地根据信号本身的特点灵活选取以得到信号非常稀疏的表示. 后来人们提出了诸如基追踪算法、匹配追踪算法（MP）、正交匹配追踪算法（OMP）、混合匹配追踪算法（HMP）及许多变种。涉及的原子包括多尺度Gabor函数，各向异性的精细原子，小波和正弦函数的级联[24-15]等，并通过训练方法获得结构和纹理分量稀疏表示字典[26-28] 。 
目前图像稀疏表示的研究也引起国内众多研究者的关注。中科院杨谦、汪云九等人，中科院计算所史忠植研究员，西安电子科技大学的焦李成教授、华南理工大学谢胜利教授，西南交通大学尹忠科教授等，南京理工大学韦志辉教授，肖亮博士等纷纷展开了稀疏表示的相关问题的研究。 目前图像稀疏表示的研究成为近3年国内众多研究者关注的热点问题，根据<<中国期刊全文数据库>>的检索来看，在2004年之前几乎没有相关报道，而从2004年1月至2008年2月，中国期刊发表的图像稀疏表示与多尺度几和分析应用方面的论文达到187篇,其中关于Ridgelet 56篇，关于Contourlet 63篇，关于 Curvelet 34篇，关于过完备稀疏表示34篇。西安电子科技大学的焦李成教授、华南理工大学谢胜利教授，西安交通大学尹忠科教授、国防科技大学王正明、教授及课题组成员等纷纷展开了基于稀疏表示的相关应用问题的研究[29-33]。本文作者在基于多尺度几何分析的图像增强、去噪、融合、边缘检测、感知压缩和数字水印等展开了相关应用研究，研究结果表明，基于稀疏表示的形态分量分解理论能够很好的捕获图像的几何特征，在图像建模和处理方面具有先天优势。但是综观国内的这些研究还与国外原创性成果具有很大差距。特别在稀疏表示字典的构造、高效稀疏分解算法、稀疏性重建等层面均有大量工作可做。

1.3 形态分量分析暂露头角

MCA方法是国际著名学者J.-L. Starck, M. Elad, D.L. Donoho 在2004年提出的一种将图像分解为 “几何结构”、“纹理”、“噪声”的形态分量分解方法[34]。该方法与混叠信号盲分离在本质上近乎相同，和独立分量分析（ICA）具有紧密联系。在MCA提出之前，图像分解的研究如火如荼。主要包括“基于稀疏表示的图像分解”和“基于变分方法的图像分解”。MCA方法较好的结合了变分方法和稀疏表示方法两类图像分解的优点，为不适定图像处理问题提供了良好的处理机制。
首先从关于图像形态分量分解的变分方法来看，国际上研究的研究朝着对图像结构和纹理等形态成分刻画更精细、计算更简单的方向发展。图像分解的（BV,G）模型,（BV,F）模型和（BV,E）模型在本质上就是一种形态分量分析方法。
与基于变分方法的图像分解处理思路不同，J.L.Stack,M.Elad 和 D.L.Donoho的MCA框架中，一个重要的假设是图像的几何结构和纹理分量在某个特定的基库或过完备子字典下是类内稀疏的，而各形态分量稀疏表示的基库或过完备子字典之间具有不相干性。通过关于结构分量和纹理分量的分类稀疏表示稀疏的强稀疏性（l0 范数或l1 范数度量）达到图像形态分量的有效分离。由于目前所涉及的稀疏表示系统有三类：正交系统（如DCT,DWT）；冗余系统（如Curvelet, Contoulet）；过完备字典（如AR-Gauss混合字典）。随着稀疏表示理论的发展，通过不同的分类稀疏表示字典、稀疏性度量和正则化方法，可以导出不同的图像形态分量分析算法[35]。之后学者们对MCA中形态成分稀疏性和多样性[36]、自适应投影阈值选取[37]等问题作了研究，并推广到多通道情形[38] 。
1.4 统计模型经久不衰
关于自然图像“统计建模方法”的研究由来已久。早期的研究工作，很大程度上受David Field在20世纪80年代中期的一个重要发现所推动：自然图像的滤波器响应总是呈现出较大的峰度的统计性质[39]。经典小波分析之所以在信号和图像处理中取得重大成功，一个比较重要的因素是对小波变换域统计模型的研究取得重大进展，主要包括小波域的MRF模型，小波域隐马尔科夫模型和分层隐马尔科夫模型等。随着多尺度几何分析的兴起，人们对于Ridgelet、Curvelet、Bandlet, Contourlet变换域的统计模型相当关注。事实上，稀疏表示系数的直方图的峰度要远远大于3，呈现明显的非高斯性，这表明非高斯性蕴含图像的本质特性。
例如，对Cameraman图像的Contourlet系数进行分析。观察上面的分布可以发现，Contourlet系数呈现明显的重尾分布。考察直方图的峰度(Kurtosis)

经计算，峰度值远远大于典型的高斯分布Kurtosis值（大约为3）。
       许多单变量先验模型比如广义高斯模型、学生t-distribution模型已经被成功地用于自然图像的小波系数的这种非高斯统计特性。最近，学者们开始关注自然图像滤波器响应的联合统计行为。Zhu S.C较为全面的论述了自然图像视觉模式的四种主流的统计研究方法，并从信号的稀疏表示出发论述了包括多个Markov随机场模型及其变种[40]。焦李成等比较研究了多尺度变换域包括隐马尔科夫树（HMT）、背景隐马尔科夫模型（CHMM）等在内的10种统计模型[41]。

[1] A.Buades, B.Coll, J.M.Morel, A review of image denoising algorithms, with a new one. Multiscale Modeling and Simulation, 2005,4(2) 490-530.
[2] L. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D, 1992,60:259-268.
[3] Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, University.Lecture Series, Vol. 22, Amer. Math. Soc., 2001.
[4] L. A. Vese, and S. J. Osher, Image Denoising and Decomposition with Total Variation Minimization and Oscillatory Functions. Journal of Mathematical Imaging and Vision, 2004, 20(1):7-18.
[5] S. Osher, A. Sol´e, and L. Vese, Image Decomposition and Restoration Using Total Variation Minimization and the H−1 Norm. Multiscale Modeling and Simulation, 2003, 1(3): 349-370.
[6] L. Lieu and L. Vese, Image Restoration and Decomposition via Bounded Variation and Negative Hilbert-Sobolev Spaces, UCLA CAM Report, 05-33, May 2005.
[7] J.F. Aujol, G. Aubert, L. Blanc-F´eraud, and A. Chambolle, Image decomposition into a bounded variation component and an oscillating component, Journal of Mathematical Imaging and Vision, 2005, 22(1): 71-88.
[8] J.B.Garnett, M.L.Triet, Y.Meyer, L.Vese. Image Decomposition using bounded variation generalized Homogeneous Besov spaces. 2005:UCLA CAM Report 05-57.
[9] J.F Aujol and A. Chambolle. Dual norms and image decomposition models. International Journal ofComputerVision,2005, 63(1):85-104.
[10] G.Giboa, S.Osher, Non-local linear image reconstruction and supervised segmentation. SIAM Multiscale Modeling and simulation, 2007, 6(2):595-630.
[11] H.B.Barlow, Possible principles underlying the transformation of sensory messages. Sensory Communication. Edited by W A Rosenbluth ( Cambridge, MA: MIT Press) 1961, 217-234.
[12] S. Nirenberg, S. M Carcieri, A .L Jacobs, P. E Latham. Retinal ganglion cell sact largely as independent encoders. Nature, 2001, 411: 698-701.
[13] William E. Vinje, Jack L. Gallant. Sparse Coding and Decorrelation in Primary Visual Cortex During Natural Vision. Science 18 February 2000:col. 287. no. 5456, pp. 1272-1276.
[14] Olshausen B. A, Field D. J. Sparse coding of sensory inputs. Current Opinion in Neurobiology. 2004, 14: 481-487.
[15] Bell A J and Sejnowski T J. The ‘independent components’ of natural scenes are edge filters. Vision Research. 1997,37: 3327-3338.
[16] Hyvarinen A, Hoyer P. O, A two-layer sparse coding model learns simple and complex cell receptive fields and topography from natural images. Vision Research. 2001, 41(18): 2413-2423.
[17] Candes E J. Ridgelet:theory and application. Ph.D dissertation, Stanford Univ.,1998.
[18] J.-L. Starck, E. J. Candès, and D. L. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Processing, vol. 11, pp. 670–684,June 2002.
[19] Erwan Le Pennec and Stéphane Mallat, Sparse Geometric Image Representations With Bandelets. IEEE Trans. Image Processing, 2005, 14(4):423-438.
[20] Do.M.N,Vertterli.M. Framing pyramids. IEEE Transactions on Signal Processing, 2003,14(9):2329-2342.
[21] Do.M.N,Vertterli.M. The contourlet transform: an efficient directional multiresilution image representation. IEEE Transactions on Image Processing,2005,14(12):2091-2106.
[22] 焦李成，谭山. 图像的多尺度几何分析：回顾和展望. 电子学报，2003;31(12A):1975-1981.
[23] S. Mallat and Zhang, Matching pursuit with time-frequency dictionaries. IEEE Transactions on Signal Processing. 1993, 41(12): 3397–3415.
[24] M. F.Rosa V. Pierre Low-rate and flexible image coding with redundant representation. IEEE Transactions on Image Processing,2006,15(3):726-739.
[25] Xu Peng,Yao Dezhong Two dictionaries matching pursuit for sparse decomposition of signals
Signal Processing, 2006，86(11) : 3472-3480.
[26] M. Elad, A.Michal. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image Processing, 2006, 15(12) : 3736-3745
[27] A.Michal， M.Elad,; B.Alfred. K-SVD: An algorithm for designing over-complete dictionaries for sparse representation. IEEE Transactions on Signal Processing, 2006, (54)11: 4311-4322
[28] G.Monaci, P. Jost, P.Vandergheynst, etal. Learning Multimodal Dictionaries. IEEE Transactions on Image Processing, 2007,16(9): 2273-2283.
[29] 何昭水，谢胜利，傅予力. 信号的稀疏性分析.自然科学进展，2006，16（9）：1167-1173.
[30] 谢胜利，何昭水，傅予力.基于稀疏元分析的欠定混叠自适应盲分离方法.中国科学E,2007,37(8) : 1086~1098.
[31] 尹忠科等.利用FFT实现基于MP的信号稀疏分解,电子与信息学报.2006，28(4)：614-618.
[32] 汪雄良,冉承其,王正明. 基于紧致字典的基追踪方法在SAR图像超分辨中的应用,电子学报，2006，34(6):997-1000.
[33] 杜小勇,胡卫东,郁文贤. 基于稀疏分量分析的逆合成孔径雷达成像技术.电子学报，2006，34(3):491-495.
[34] J. L. Starck, M. Elad, and D.L. Donoho. Redundant multiscale transforms and their application for morphological component analysis. Advances in Imaging and Electron Physics, 2004,132
[35] J.L. Starck, M. Elad, and D.L. Donoho.Image decomposition via the combination of sparse representation and a variational approach. IEEE Transaction on Image Processing, 2005,14(10):1570-1582.
[36] J.Bobin1, Y.Moudden1, J.Fadili2 and J-L.Starck. Morphological diversity and sparsity in blind source separation. IEEE Transactions on Image Processing, 2007, 16(1):2662 – 2674.
[37] J. Bobin, J.-L. Starck, J. Fadili, Y. Moudden, and D.L. Donoho. Morphological component analysis : an adaptive thresholding strategy. IEEE Transactions on Image Processing, 2007, 16(1):2675 - 2681.
[38] J. Mairal, M. Elad, G. Sapiro. Sparse representation for color image restoration. IEEE Transactions on Image Processing, 2008，17(1):53-68.
[39] B. A. Olshausen, D. J Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 1996, 381:607-609.
[40] Zhu S.C. Statistical Modeling and Conceptualization of Visual Patterns [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2003, 25(6):691-712.
[41] 焦李成，孙强.多尺度变换域图象的感知与识别：进展和展望.计算机学报.2006，29(2):177-193.

要说哪一个方法被做的全面，那莫非LE莫属。如果只说LE这个方法本身，是不新的，许多年前在做mesh相关的领域就开始这莫用。但是放在黎曼几何的框架内，给出完整的几何分析的，应该是Belkin和Niyogi（LE作者）的功劳。
LE的基本思想就是用一个无向有权图来描述一个流形，然后通过用图的嵌入（graph embedding）来找低维表示。说白了，就是保持图的局部邻接关系的情况把这个图从高维空间中重新画在一个低维空间中（graph drawing）。
在至今为止的流行学习的典型方法中，LE是速度最快、效果相对来说不怎莫样的。但是LE有一个其他方法没有的特点，就是如果出现outlier情况下，它的鲁棒性（robustness）特别好。