2018-04-27 14:41:53 jayandchuxu 阅读数 342

矩阵的特征值、特征向量的概念

这里,我们讨论的是n阶的方阵A

定义

从向量的定义可知,它是方向和长度的结合体。当一个线性变换A作用在n维线性空间V中的某一非零向量x上时,便是对该向量的长度和方向进行变化。然而,存在一些向量,线性变换A并没有改变其方向,而只是改变了长度,这种向量,叫做线性变换A特征向量,它在变换中被改变的倍数,叫做它的特征值。用数学公式表示这一概念,即:

Ax=λx(1)

其中,λ的个线性变换A的某一个特征值。从公式上可以轻易发现,如果某一向量x是线性变换A的特征向量,那么与其方向相同的任意长度(不为零)向量,都是A的特征向量,并且属于同一个特征值λ。由于相同方向的特征向量具有相同特征值,我们可以同特征值来描述这一族向量(同一个特征值,可以有多个方向的特征向量)。
从公式(1)中,可以看出特征向量和特征值的计算方法:
|λEA|=0(2)
(λEA)x=0(3)

对应于同一个线性变换A,可以有多个特征向量(方向不同),但是有多个特征向量可以对应同一个特征值。 一个向量是一个方向,两个不同方向的向量就可以张成一个空间。在相同特征值的特征向量张成的空间内,任何一个向量在变换A下,都有相同的放大倍数λ

公式(2)的左侧,总可以展成如下形式的多项式:

|λEA|=λn(a11+a22++ann)λn1++(1)n|A|

所以求特征值就是求下面方程的解:

λn(a11+a22++ann)λn1++(1)n|A|=0(4)

关于从方程(4)得到的特征值,有几个比较重要的结论(参考资料1):

  1. n阶矩阵在复数范围内,一定有n个特征值(重特征值按重数计算个数)。
  2. n阶矩阵在实数范围内有多少个特征值是不一定的
  3. n阶实对称矩阵可以看成是一个特例,因为它一定有n个实特征值(重特征值按重数计算个数)。如果其中一个特征值λ=0,矩阵的秩r(A)=k,(0<k<nk是正整数),则λ=0恰为Ank重特征值。
  4. 如果 n阶矩阵A不是对称矩阵,那么,λ=0至少为Ank重特征值。

关于特征值和特征向量的理解,参考资料3写的也很好,知乎上有很多大神的回答直击要害,对问题的理解很有帮助。

作用

参考资料4中,认为矩阵的变换有三个作用:旋转拉伸投影
A是一个n×n方阵时,只涉及到旋转,拉伸。如果矩阵在实数域内可以得到n个特征值(重特征值按重数计算个数),那么利用其对应的特征向量(单位化后)组成矩阵正交Q可以使得:

A=QΛQ1

其中正交矩阵Q起到旋转作用(旋转矩阵都是正交矩阵,且行列式都为1),对角矩阵Λ起拉伸作用。
当矩阵不是方阵而是m×n时,可以对其进行SVD分解
在之前,想研究一下正交矩阵。

正交矩阵

按照定义,正交矩阵是QQT=E,它的行列式为1或者-1。

正交矩阵的性质

了解正交矩阵的性质,在很多计算方面,能够更深入了解所进行的运算的意义。
我们这里说的都是有限维欧式空间内的正交矩阵

  1. 正交矩阵的转置伴随矩阵、之间的积矩阵都是正交矩阵;
  2. .每一行(列)都是单位向量
  3. 任意两行或两列相互垂直
  4. 其行列式等于±1

假设n维欧式空间Rn中,一个正交变换Q存在一个一一对应的正交矩阵Q。所以研究正交变换的性质,可以转为研究正交矩阵。
根据正交矩阵行列式的值,将其分为两类:|Q|=1,为第一类;|Q|=1,为第二类。
第一类正交矩阵,当其左乘一个向量时,几何意义是使该量在Oxyz坐标系下旋转;
第一类正交矩阵,当其左乘一个向量时,几何意义是使该向量沿Oxyz某一轴(点)进行反射;[5]
无论是哪一类正交矩阵,其左乘向量,均不会改变向量的长度,即|Qv|=|Q||v|=|v|
所以上面将矩阵A拆成A=QΛQ1的形式后,由于Q是正交矩阵,它作用在向量v上,只是对其进行旋转或者反射,而对向量长度有影响的是矩阵Λ。其上的元素由矩阵A的特征值组成。
需要注意的是,只有当矩阵A能够在实数域内有n个特征值(重数也计入)的情况下,可以如此分解。但是,对更多的一般性矩阵A,在实数域内没有n个实特征根,或者是更一般的m×n维矩阵,此时,我们用SVD分解的方法进行研究。

SVD分解

协方差

假设有两个变量XY,他们的取值为X={x1,x2,,xn}Y={y1,y2,,yn},他们的平均值(期望)记为是E(X)=X¯=1nxiE(Y)=Y¯=1nyiXY的协方差就是研究两个变量之间的相关关系。比如当一个的取值不断增大时,另一个变量的取值如何变化,知乎上的一篇文章( 如何通俗易懂地解释「协方差」与「相关系数」的概念?)讲的很好。

根据奇异值分解

n×n阶矩阵按特征值分解相似,任一m×n阶矩阵A也可以写成类似的形式:

A=UΣVT(5)

那么得到的U是一个m×m的方阵(里面的向量是正交的,U里面的向量称为左奇异向量),Σ是一个m×n的矩阵(除了对角线的元素都是0,对角线上的元素称为奇异值),VT(V的转置)是一个n×n的矩阵,里面的向量也是正交的,V里面的向量称为右奇异向量)。

求法

利用如下公式可以计算出各矩阵:

(ATA)vi=λivi(6)

δi=λi(7)

ui=1δiAvi(8)

对方程计算,得到的v,就是的右奇异向量σ就是奇异值u就是左奇异向量

作用

当矩阵A=UΣVT左乘一个向量v时,只有Σ对向量的长度进行了拉伸(收缩),而矩阵UV都只是对其进行旋转或反射。当作用在图像上时,也只有Σ对图像进行了各个方向上的伸缩改变。

用奇异值分解图像

这里写图片描述
用奇异值的方法,将这幅图像进行分解,得形如A=UΣVT格式的矩阵。其中Σ是由矩阵奇异值由大到小排列组成的对角矩阵。
我们分别保留前10,30,100,300个奇异值,其余奇异值设为0,比较图像的变化:
奇异值保留前10
奇异值保留前10
奇异值保留前30
奇异值保留前30
奇异值保留前100
奇异值保留前100
奇异值保留前300
奇异值保留前300。

代码

import cv2
import numpy as np
from numpy import linalg as la  # 用到别名
from scipy.misc import imsave

import scipy
im = cv2.imread('lena512.bmp')
print(im.shape)
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
U, Sigma, VT = la.svd(gray)
print('矩阵U的形状:', U.shape, '    矩阵Sigma的形状:',
      Sigma.shape, '    矩阵VT的形状:', VT.shape)
se = np.eye(512, dtype=np.float64)
n = 512
i = 0
k = 30  # 保留特征值数目
# 改变特征值
while i < n:
    if i > k - 1:
        Sigma[i] = 0
    se[i, i] = Sigma[i]
    i += 1

svt = np.dot(se, VT)
usvt = np.dot(U, svt)
imsave('USVT_Sigma=30m.bmp', usvt)

参考资料

  1. 秦川, 李小飞. 方阵的秩与特征值的关系[J]. 课程教育研究:学法教法研究, 2015(27):120-120.
  2. 如何通俗易懂地解释「协方差」与「相关系数」的概念?
  3. 如何理解矩阵特征值?马同学的回答
  4. 矩阵的特征值分解与奇异值分解的几何意义
  5. 杜美华, 孙建英. 正交变换的几何意义及其应用[J]. 哈尔滨师范大学自然科学学报, 2014, 30(3):36-39.
2018-05-04 09:49:58 weixin_42026802 阅读数 728

   数学不行就先拿别人的看看吧,比较全:

滤波:FName1为原始数据文件名,FName2为滤波后数据文件名。procedure TForm.Smooth(FName1,FName2:string);var  L1 : longint;  I, Icount, Fp1, Fp2: integer;  pre:array [1..14] of Byte;  Last:array [1..3] of Byte;  X, Y ,Z: buffer;  F:file of byte;begin  assignfile(F,FName1);  reset(F);  datalength:=filesize(F);  closefile(f);  Fp1:=fileopen(Fname1,0);  Icount:=fileread(Fp1,X,14);  L1:=Icount;  assignfile(F,fname2);  rewrite(F);  closefile(F);  Fp2:=fileopen(Fname2,1);  for I:=1 to 14 do    pre[I]:=X[I];  Y[1]:=X[1];  Y[2]:=(X[2]+2*X[1])div 3;  Y[3]:=(x[3]+2*x[2]+3*x[1])div 6;  Y[4]:=(x[4]+2*x[3]+3*x[2]+4*x[1])div 10;  Y[5]:=(x[5]+2*x[4]+3*x[3]+4*x[2]+5*x[1])div 15;  Y[6]:=(x[6]+2*x[5]+3*x[4]+4*x[3]+5*x[2]+6*x[1])div 21;  Y[7]:=(x[7]+2*x[6]+3*x[5]+4*x[4]+5*x[3]+6*x[2]+7*x[1])div 28;  Y[8]:=(x[8]+2*x[7]+3*x[6]+4*x[5]+5*x[4]+6*x[3]+7*x[2]+8*x[1])div 36;  Y[9]:=(x[9]+2*x[8]+3*x[7]+4*x[6]+5*x[5]+6*x[4]+7*x[3]+8*x[2]+7*x[1])div 43;  Y[10]:=(x[10]+2*x[9]+3*x[8]+4*x[7]+5*x[6]+6*x[5]+7*x[4]+8*x[3]+7*x[2]+6*x[1])div 49;  Y[11]:=(x[11]+2*x[10]+3*x[9]+4*x[8]+5*x[7]+6*x[6]+7*x[5]+8*x[4]+7*x[3]+6*x[2]+5*x[1])div 54;  Y[12]:=(x[12]+2*x[11]+3*x[10]+4*x[9]+5*x[8]+6*x[7]+7*x[6]+8*x[5]+7*x[4]+6*x[3]+5*x[2]+4*x[1])div 58;  Y[13]:=(x[13]+2*x[12]+3*x[11]+4*x[10]+5*x[9]+6*x[8]+7*x[7]+8*x[6]+7*x[5]+6*x[4]+5*x[3]+4*x[2]+3*x[1])div 61;  Y[14]:=(x[14]+2*x[13]+3*x[12]+4*x[11]+5*x[10]+6*x[9]+7*x[8]+8*x[7]+7*x[6]+6*x[5]+5*x[4]+4*x[3]+3*x[2]+2*x[1])div 63;  Filewrite(fp2,Y,14);  last[1]:=Y[12];  last[2]:=Y[13];  last[3]:=Y[14];  while L1<dataLength do    begin     Icount:=fileread(Fp1,X,N);     L1:=L1+Icount;     Y[1]:=(x[1]+2*pre[14]+3*pre[13]+4*pre[12]+5*pre[11]+6*pre[10]+7*pre[9]+8*pre[8]+7*pre[7]+6*pre[6]+5*pre[6]+4*pre[4]+3*pre[3]+2*pre[2]+1*pre[1])div 64;     Y[2]:=(x[2]+2*X[1]+3*pre[14]+4*pre[13]+5*pre[12]+6*pre[11]+7*pre[10]+8*pre[9]+7*pre[8]+6*pre[7]+5*pre[6]+4*pre[6]+3*pre[4]+2*pre[3]+1*pre[2])div 64;     Y[3]:=(X[3]+2*x[2]+3*X[1]+4*pre[14]+5*pre[13]+6*pre[12]+7*pre[11]+8*pre[10]+7*pre[9]+6*pre[8]+5*pre[7]+4*pre[6]+3*pre[6]+2*pre[4]+1*pre[3])div 64;     Y[4]:=(X[4]+2*X[3]+3*x[2]+4*X[1]+5*pre[14]+6*pre[13]+7*pre[12]+8*pre[11]+7*pre[10]+6*pre[9]+5*pre[8]+4*pre[7]+3*pre[6]+2*pre[5]+1*pre[4])div 64;     Y[5]:=(X[5]+2*X[4]+3*X[3]+4*x[2]+5*X[1]+6*pre[14]+7*pre[13]+8*pre[12]+7*pre[11]+6*pre[10]+5*pre[9]+4*pre[8]+3*pre[7]+2*pre[6]+1*pre[5])div 64;     Y[6]:=(X[6]+2*X[5]+3*X[4]+4*X[3]+5*x[2]+6*X[1]+7*pre[14]+8*pre[13]+7*pre[12]+6*pre[11]+5*pre[10]+4*pre[9]+3*pre[8]+2*pre[7]+1*pre[6])div 64;     Y[7]:=(X[7]+2*X[6]+3*X[5]+4*X[4]+5*X[3]+6*x[2]+7*X[1]+8*pre[14]+7*pre[13]+6*pre[12]+5*pre[11]+4*pre[10]+3*pre[9]+2*pre[8]+1*pre[7])div 64;     Y[8]:=(X[8]+2*X[7]+3*X[6]+4*X[5]+5*X[4]+6*X[3]+7*x[2]+8*X[1]+7*pre[14]+6*pre[13]+5*pre[12]+4*pre[11]+3*pre[10]+2*pre[9]+1*pre[8])div 64;     Y[9]:=(X[9]+2*X[8]+3*X[7]+4*X[6]+5*X[5]+6*X[4]+7*X[3]+8*x[2]+7*X[1]+6*pre[14]+5*pre[13]+4*pre[12]+3*pre[11]+2*pre[10]+1*pre[9])div 64;     Y[10]:=(X[10]+2*X[9]+3*X[8]+4*X[7]+5*X[6]+6*X[5]+7*X[4]+8*X[3]+7*x[2]+6*X[1]+5*pre[14]+4*pre[13]+3*pre[12]+2*pre[11]+1*pre[10])div 64;     Y[11]:=(X[11]+2*X[10]+3*X[9]+4*X[8]+5*X[7]+6*X[6]+7*X[5]+8*X[4]+7*X[3]+6*x[2]+5*X[1]+4*pre[14]+3*pre[13]+2*pre[12]+1*pre[11])div 64;     Y[12]:=(X[12]+2*X[11]+3*X[10]+4*X[9]+5*X[8]+6*X[7]+7*X[6]+8*X[5]+7*X[4]+6*X[3]+5*x[2]+4*X[1]+3*pre[14]+2*pre[13]+1*pre[12])div 64;     Y[13]:=(X[13]+2*X[12]+3*X[11]+4*X[10]+5*X[9]+6*X[8]+7*X[7]+8*X[6]+7*X[5]+6*X[4]+5*X[3]+4*x[2]+3*X[1]+2*pre[14]+1*pre[13])div 64;     Y[14]:=(X[14]+2*X[13]+3*X[12]+4*X[11]+5*X[10]+6*X[9]+7*X[8]+8*X[7]+7*X[6]+6*X[5]+5*X[4]+4*X[3]+3*x[2]+2*X[1]+1*pre[14])div 64;     Z[1]:=(Y[1]+LAST[3]+LAST[2]+LAST[1]) div 4;     Z[2]:=(Y[2]+Y[1]+LAST[3]+LAST[2]) div 4;     Z[3]:=(Y[3]+Y[2]+Y[1]+LAST[3]) div 4;     for I :=15 to ICount do      Y[I]:=(X[I]+2*X[I-1]+3*X[I-2]+4*X[I-3]+5*X[I-4]+6*X[I-5]+7*X[I-6]+8*X[I-7]+7*X[I-8]+6*X[I-9]+5*X[I-10]+4*X[I-11]+3*x[I-12]+2*X[I-13]+1*X[I-14]) div 64;     for I:= 4 to Icount do      Z[I]:=(Y[I]+Y[I-1]+Y[I-2]+Y[I-3]) div 4;     for I:=1 to 14 do      pre[I]:=X[Icount-I+1];      last[1]:=Y[Icount-2];      Last[2]:=Y[Icount-1];      Last[3]:=Y[Icount];      filewrite(fp2,Z,Icount);    end;      fileclose(fp2);      fileclose(fp1);end;

unit hhx_effectex;interfaceuses hhx_Effects;var mxEmbossColor:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:        (( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0,-1,-1,-1, 0, 0),        ( 0, 0, 0, 1, 0, 0, 0),        ( 0, 0, 1, 1, 1, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:0;        FilterName:'彩色浮雕';);var mxEmbossLight:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:192;        FilterName:'高亮浮雕';);var mxEmbossMedium:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:192;        FilterName:'中值浮雕';);var mxEmbossDark:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:128;        FilterName:'黑色浮雕';);var mxEdgeEnhance:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0,-2,16,-2, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:4;        Bias:0;        FilterName:'边缘增强 (线性锐化)';);var mxBlurBartlett:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx7;        Matrix:      (( 1, 2, 3, 4, 3, 2, 1),      ( 2, 4, 6, 8, 6, 4, 2),      ( 3, 6, 9,12, 9, 6, 3),      ( 4, 8,12,16,12, 8, 4),      ( 3, 6, 9,12, 9, 6, 3),      ( 2, 4, 6, 8, 6, 4, 2),      ( 1, 2, 3, 4, 3, 2, 1));        Divisor:256;        Bias:0;        FilterName:'Bartlett模糊 (线性模糊)';);var mxBlurGaussian:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx7;        Matrix:     (( 1, 4, 8, 10, 8, 4, 1),      ( 4,12,25,29,25,12, 4),      ( 8,25,49,58,49,25, 8),      (10,29,58,67,58,29,10),      ( 8,25,49,58,49,25, 8),      ( 4,12,25,29,25,12, 4),      ( 1, 4, 8, 10, 8, 4, 1));      Divisor:999;      Bias:0;      FilterName:'高斯模糊 (线性模糊)';);var mxNegative:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:255;      FilterName:'反相 (线性效果)';);var mxAverage:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:9;      Bias:0;      FilterName:'平均值滤波 (线性模糊)';);var mxBlur:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 2, 4, 2, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:16;      Bias:0;      FilterName:'模糊 Blur';);var mxBlurSoftly:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 3,16, 3, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:32;      Bias:0;      FilterName:'轻度模糊 Blur softly';);var mxBlurMore:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx5;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 1, 4, 6, 4, 1, 0),      ( 0, 2, 6, 8, 6, 2, 0),      ( 0, 1, 4, 6, 4, 1, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:64;      Bias:0;      FilterName:'进一步模糊 Blur more';);var mxPrewitt:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1,-2, 1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Prewitt边缘 (线性边缘检测)';);var mxTraceContour:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0,-1,32,-1, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:240;      FilterName:'轮廓描绘 (线性边缘检测)';);//      FilterName:'Trace contour (Edge detect linear)';);var mxSharpen:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,16,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:8;      Bias:0;      FilterName:'锐化 (线性锐化)';);var mxSharpenMore:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,12,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:4;      Bias:0;      FilterName:'进一步锐化 (线性锐化)';);var mxSharpenLess:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,24,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:16;      Bias:0;      FilterName:'轻度锐化 (线性锐化)';);var mxUnSharpMask:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0,-2,16,-2, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:4;      Bias:0;      FilterName:'非锐化模板 (线性锐化)';);//      FilterName:'非锐化模板Unsharp mask (Sharpen linear)';);

var mxEdgesStrong:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 3,-16,3, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'强化边缘 (线性边缘检测)';);var mxEdgesWeak:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 1,-4, 1, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'边缘弱化 (线性边缘检测)';);var mxEtch:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-6, 2,-6, 0, 0),      ( 0, 0,-1,32,-1, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:240;      FilterName:'腐蚀';);//      FilterName:'Etch (Effects linear)';);var mxLaplacianHV:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0,-1, 4,-1, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Laplacian水平/纵向边缘 (线性边缘检测)';);//      FilterName:'Laplacian horz./vert. (Edge detect linear)';);var mxLaplacianOmni:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1, 8,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Laplacian所有方向边缘?(线性边缘检测)';);//      FilterName:'Laplacian omnidir? (Edge detect linear)';);var mxSharpenDirectional:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-3,-3,-3, 0, 0),      ( 0, 0, 0,16, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:10;      Bias:0;      FilterName:'有方向的锐化';);//      FilterName:'Sharpen directional (Sharpen linear)';);var mxSobelPass:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Sobel边缘 (线性边缘检测)';);//      FilterName:'Sobel pass (Edge detect linear)';);//******************************************************************************////  Two-Pass filters////******************************************************************************var nmxLaplacianInvert:TMultiPassGraphicFilter=      (FilterType:ftMultiPass;       Filters:(@mxLaplacianOmni,@mxNegative,@mxZero,@mxZero);       Functions:(gaNone,gaNone,gaNone);       FilterName:'Laplacian Negative');implementationend.

2016-04-21 20:44:32 yangdashi888 阅读数 2502

6.5 矩阵的运算及其运算规则

一、矩阵的加法与减法


  1、运算规则
  设矩阵
  则
     
  简言之,两个矩阵相加减,即它们相同位置的元素相加减!
  注意:只有对于两个行数、列数分别相等的矩阵(即同型矩阵),加减法运算才有意义,即加减运算是可行的.

  2、 运算性质 (假设运算都是可行的)
  满足交换律和结合律
  交换律 
  结合律 

二、矩阵与数的乘法


  1、 运算规则
  乘矩阵A,就是将数乘矩阵A中的每一个元素,记为
  特别地,称称为的负矩阵.
  2、 运算性质
  满足结合律和分配律
  结合律: (λμ)A=λ(μA) ; (λ+μ)A =λA+μA
  分配律: λ (A+B)=λA+λB

  典型例题
  例6.5.1 已知两个矩阵
  满足矩阵方程,求未知矩阵
   由已知条件知
    
    

三、矩阵与矩阵的乘法


  1、 运算规则
  设,则A与B的乘积是这样一个矩阵:
  (1) 行数与(左矩阵)A相同,列数与(右矩阵)B相同,即
  (2) C的第行第列的元素由A的第行元素与B的第列元素对应相乘,再取乘积之和.

  典型例题
  例6.5.2 设矩阵
  计算
   的矩阵.设它为
    

    
  想一想:设列矩阵,行矩阵的行数和列数分别是多少呢
  是3×3的矩阵,是1×1的矩阵,即只有一个元素.

  课堂练习
  1、设,求
  2、在第1道练习题中,两个矩阵相乘的顺序是A在左边,B在右边,称为A左乘B或B右乘A.如果交换顺序,让B在左边,A在右边,即A右乘B,运算还能进行吗?请算算试试看.并由此思考:两个矩阵应当满足什么条件,才能够做乘法运算.
  3、设列矩阵,行矩阵,求,比较两个计算结果,能得出什么结论吗?
  4、设三阶方阵,三阶单位阵为,试求,并将计算结果与A比较,看有什么样的结论.

  解:
  第1题
  第2题
  对于
  求是有意义的,而是无意义的.

  结论1 只有在下列情况下,两个矩阵的乘法才有意义,或说乘法运算是可行的:左矩阵的列数=右矩阵的行数.
  第3题
  矩阵,的矩阵.
         
           
    结论2 在矩阵的乘法中,必须注意相乘的顺序.即使在均有意义时,也未必有=成立.可见矩阵乘法不满足交换律.
  第4题
  计算得:
  结论3 方阵A和它同阶的单位阵作乘积,结果仍为A,即
  单位阵在矩阵乘法中的作用相当于数1在我们普通乘法中的作用.

  典型例题
  例6.5.3 设,试计算
   
      
      
    
      
      
    结论4 两个非零矩阵的乘积可以是零矩阵.由此若,不能得出的结论.

  例6.5.4 利用矩阵的乘法,三元线性方程组
  可以写成矩阵的形式
  若记系数、未知量和常数项构成的三个矩阵分别为
  则线性方程组又可以简写为矩阵方程的形式:

  2、 运算性质(假设运算都是可行的)
  (1) 结合律 
  (2) 分配律 (左分配律);
         (右分配律).
  (3) 
   3、 方阵的幂
 
定义:设A是方阵,是一个正整数,规定
显然,记号表示个A的连乘积.

四、矩阵的转置


  1、 定义
 
定义:将矩阵A的行换成同序号的列所得到的新矩阵称为矩阵A的转置矩阵,记作
  例如,矩阵的转置矩阵为
  2、运算性质(假设运算都是可行的)
  (1) 
  (2) 
  (3) 
  (4) 是常数.

  典型例题
  例6.5.5  利用矩阵
  验证运算性质:
   
  而
    
  所以
   

 
定义:如果方阵满足,即,则称A为对称矩阵
  对称矩阵的特点是:它的元素以主对角线为对称轴对应相等.

五、方阵的行列式


  1、定义
 
定义:由方阵A的元素所构成的行列式(各元素的位置不变),称为方阵A的行列式,记作

  2 、运算性质
  (1) (行列式的性质)
  (2) ,特别地:
  (3) 是常数,A的阶数为n)
  思考:设A为阶方阵,那么的行列式与A的行列式之间的关系为什么不是,而是

  不妨自行设计一个二阶方阵,计算一下
  例如,则
  于是,而
  思考:,有几种方法可以求
    方法一:先求矩阵乘法,得到一个二阶方阵,再求其行列式.
    方法二:先分别求行列式,再取它们的乘积.

3、逆矩阵
     逆矩阵: 设A是数域上的一个n阶方阵,若在相同数域上存在另一个n阶矩阵B,使得: AB=BA=E。 则我们称B是A的逆矩阵,而A则被称为可逆矩阵。其用A^(-1)表示。
   逆矩阵的性质为:
    1 矩阵A可逆的充要条件是A的行列式不等于0。
   2 可逆矩阵一定是方阵
3 如果矩阵A是可逆的,A的逆矩阵是唯一的。
4 可逆矩阵也被称为非奇异矩阵、满秩矩阵。
5 两个可逆矩阵的乘积依然可逆。
6 可逆矩阵的转置矩阵也可逆。
7 矩阵可逆当且仅当它是满秩矩阵

 其中的方阵:就是行和列数相等的矩阵。
行列式:其就是求一个矩阵的数值,一般使用上三角或者下三角。其如下:


满秩矩阵:一个判断线性方程是否有解的关键条件,其是其说的就是要线性无关,其可分为列满秩和行满秩。其概念是:
满秩矩阵(non-singular matrix): 设A是n阶矩阵, 若r(A) = n, 则称A为满秩矩阵。但满秩不局限于n阶矩阵。若矩阵秩等于行数,称为行满秩;若矩阵秩等于列数,称为列满秩。既是行满秩又是列满秩则为n阶矩阵即n阶方阵。其中的R(A)指的就是求其秩。其中的行满秩就是行向量线性无关;列满秩就是列向量线性无关;其示意图如下:



   3、 列(行)满秩阵与方程关系:
 
这也是为什么相机标定的时候每个视场只能提前四个有用的点,这时因为第五个点肯定跟前面的四个点其中一个线性相关,意思就是(x,y)(X,Y)这个肯定跟前面的一个点成倍数关系,这样子在求秩的时候,有成倍数的行会被化成零,则这就成了降秩了,就不是满秩了。








2018-04-01 09:56:30 u013921430 阅读数 10311

 

fishing-panhttps://blog.csdn.net/u013921430转载请注明出处】

前言

       Hessian Matrix(海森矩阵)在图像处理中有广泛的应用,比如边缘检测、特征点检测等。而海森矩阵本身也包含了大量的数学知识,例如泰勒展开、多元函数求导、矩阵、特征值等。写这篇博客的目的,就是想从原理和实现上讲一讲Hessian Matrix,肯定有不足的地方,希望大家批评指正。

泰勒展开及海森矩阵

      将一个一元函数f(x)在x0处进行泰勒展开,可以得到以下公式。

       其中余项为皮亚诺余项

       其中二阶导数的部分映射到二维以及多维空间就是Hessian Matrix。在二维图像中,假设图像像素值关于坐标(x, y)的函数是f(x, y),那么将f(x+dx,y+dy)在f(x0, y0)处展开,得到如下式子;

     如果将这个式子用矩阵表示,并且舍去余项,则式子会是下面这个样子。

      上面等式右边的第三项中的第二个矩阵就是二维空间中的海森矩阵了;从而有了一个结论,海森矩阵就是空间中一点处的二阶导数。进而推广开来,多维空间中的海森矩阵可以表示为。

 

海森矩阵的意义

     众所周知,二阶导数表示的导数的变化规律,如果函数是一条曲线,且曲线存在二阶导数,那么二阶导数表示的是曲线的曲率,曲率越大,曲线越是弯曲。以此类推,多维空间中的一个点的二阶导数就表示该点梯度下降的快慢。以二维图像为例,一阶导数是图像灰度变化即灰度梯度,二阶导数就是灰度梯度变化程度,二阶导数越大灰度变化越不具有线性性(这里有一点绕口了,意思就是这里灰度梯度改变越大,不是线性的梯度)。

      但是在二维图像中,海森矩阵是二维正定矩阵,有两个特征值和对应的两个特征向量。两个特征值表示出了图像在两个特征向量所指方向上图像变化的各向异性。如果利用特征向量与特征值构成一个椭圆,那么这个椭圆就标注出了图像变化的各向异性。那么在二维图像中,什么样的结构最具各项同性,又是什么样的结构的各向异性更强呢?很显然,圆具有最强的各项同性,线性越强的结构越具有各向异性。如下图;

注:图中箭头的方向不一定正确,我只是随意标注

       且特征值应该具有如下特性;

 

λ1

λ2

图像特征

-High

-High

斑点结构(前景为亮)

+High

+High

斑点结构(前景为暗)

Low

-High

线性结构(前景为亮)

Low

+High

线性结构(前景为暗) 

海森矩阵的应用

       海森矩阵的应用很多,我在此不一一列举。上文中提到了矩阵的特征值与特征向量构成的椭圆表现出了图像的各向异性,这种各项异性图像处理就得到了应用。以二维图像为例,图像中的点性结构具有各项同性,而线性结构具有各向异性。因此我们可以利用海森矩阵对图像中的线性结构进行增强,滤去点状的结构和噪声点。同样,也可以用于找出图像中的点状结构,滤除其他信息。

       当然,我们在使用海森矩阵时,不需要把图像进行泰勒展开,我们只需要直接求取矩阵中的元素即可。一般,对数字图像进行二阶求导使用的是以下方法;

      但是这种方法鲁棒性很差,容易受到图像中局部信号的干扰,甚至可以说,这种求导方式是存在争议的。因为这一点的二阶导数也可以采用如下方法表示;

                                              

       除了上面两种表示方法外,二阶导数还可以采用其他方式表示,而且往往不同的方法求得的值不同,因为这种方法只把包含其自身在内的三个点的信息囊括进去,信息量不足。因此,提倡大家换一种方法。根据线性尺度空间理论(LOG),对一个函数求导,等于函数与高斯函数导数的卷积。式子如下;

       由于高斯模板可以将周围一矩形范围内所有的点的信息都包含进来,这样就不会有误差。所以利用图像求取hessian矩阵中的元素时,将图像与高斯函数的二阶导数做卷积即可。在此,为大家提供高斯函数的二阶偏导。

      在编写程序时,我们只需要事先将图像分别于三个模板进行卷积,生成三种偏导数的“图”,然后每次根据需要索引对应位置的偏导数即可。

代码

        需要说明的是,我在代码中运用了一些OpenCV的函数;

///---------------------------///
//---海森矩阵二维图像增强
//---不用先生---2018.03.27

#include <iostream>
#include <vector>
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include <map>

#define STEP 6
#define ABS(X) ((X)>0? X:(-(X)))
#define PI 3.1415926

using namespace std;
using namespace cv;



int main()
{
	Mat srcImage = imread("G:\\博客\\图像处理\\hessian\\hessian_matrix\\Multiscale vessel enhancement filtering1.bmp");

	if (srcImage.empty())
	{
		cout << "图像未被读入";
		system("pause");
		return 0;
	}
	if (srcImage.channels() != 1)
	{
		cvtColor(srcImage, srcImage, CV_RGB2GRAY);
	}
	


	int width = srcImage.cols;
	int height = srcImage.rows;

	Mat outImage(height, width, CV_8UC1,Scalar::all(0));
	int W = 5;
	float sigma = 01;
	Mat xxGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
	Mat xyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
	Mat yyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));

        //构建高斯二阶偏导数模板
	for (int i = -W; i <= W;i++)
	{
		for (int j = -W; j <= W; j++)
		{
			xxGauKernel.at<float>(i + W, j + W) = (1 - (i*i) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
			yyGauKernel.at<float>(i + W, j + W) = (1 - (j*j) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
			xyGauKernel.at<float>(i + W, j + W) = ((i*j))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * PI*pow(sigma, 6)));
		}
	}


	for (int i = 0; i < (2 * W + 1); i++)
	{
		for (int j = 0; j < (2 * W + 1); j++)
		{
			cout << xxGauKernel.at<float>(i, j) << "  ";
		}
		cout << endl;
	}

	Mat xxDerivae(height, width, CV_32FC1, Scalar::all(0));
	Mat yyDerivae(height, width, CV_32FC1, Scalar::all(0));
	Mat xyDerivae(height, width, CV_32FC1, Scalar::all(0));
        //图像与高斯二阶偏导数模板进行卷积
	filter2D(srcImage, xxDerivae, xxDerivae.depth(), xxGauKernel);
	filter2D(srcImage, yyDerivae, yyDerivae.depth(), yyGauKernel);
	filter2D(srcImage, xyDerivae, xyDerivae.depth(), xyGauKernel);


	for (int h = 0; h < height; h++)
	{
		for (int w = 0; w < width; w++)
		{
			
			
				//map<int, float> best_step;
				
			/*	int HLx = h - STEP; if (HLx < 0){ HLx = 0; }
				int HUx = h + STEP; if (HUx >= height){ HUx = height - 1; }
				int WLy = w - STEP; if (WLy < 0){ WLy = 0; }
				int WUy = w + STEP; if (WUy >= width){ WUy = width - 1; }


				float fxx = srcImage.at<uchar>(h, WUy) + srcImage.at<uchar>(h, WLy) - 2 * srcImage.at<uchar>(h, w);
				float fyy = srcImage.at<uchar>(HLx, w) + srcImage.at<uchar>(HUx, w) - 2 * srcImage.at<uchar>(h, w);
				float fxy = 0.25*(srcImage.at<uchar>(HUx, WUy) + srcImage.at<uchar>(HLx, WLy) - srcImage.at<uchar>(HUx, WLy) - srcImage.at<uchar>(HLx, WUy));*/


			float fxx = xxDerivae.at<float>(h, w);
			float fyy = yyDerivae.at<float>(h, w);
			float fxy = xyDerivae.at<float>(h, w);


			float myArray[2][2] = { { fxx, fxy }, { fxy, fyy } };          //构建矩阵,求取特征值

			Mat Array(2, 2, CV_32FC1, myArray);
			Mat eValue;
			Mat eVector;

			eigen(Array, eValue, eVector);                               //矩阵是降序排列的
			float a1 = eValue.at<float>(0, 0);
			float a2 = eValue.at<float>(1, 0);

			if ((a1>0) && (ABS(a1)>(1+ ABS(a2))))             //根据特征向量判断线性结构
			{


				outImage.at<uchar>(h, w) =  pow((ABS(a1) - ABS(a2)), 4);
				//outImage.at<uchar>(h, w) = pow((ABS(a1) / ABS(a2))*(ABS(a1) - ABS(a2)), 1.5);
				
				
			}


				
		}

	}

	
//----------做一个闭操作
	Mat element = getStructuringElement(MORPH_RECT, Size(3, 2));
	morphologyEx(outImage, outImage, MORPH_CLOSE, element);
	
	imwrite("temp.bmp", outImage);

	imshow("[原始图]", outImage);
	waitKey(0);


	system("pause");
	return 0;
}

输出结果

      左侧是原图,中间是增强后的结果,右侧是将增强后的结果做了二值化的结果图;

                                                                                               

结果分析

      从结果来看,图像中的大部分的线性结构都被增强了,但是有一些细微的结构并未被增强太多,而且有些粗的结构中出现了空洞,其实这都与求导窗口的大小有关,求导窗口太小,很多粗的结构会出现中空的现象,因为中心区域被认为是点结构了;求导窗口太大,就容易出现细微结构丢失的情况。此外,高斯模板的方差选取也影响了偏导数的大小。其实,这种方式是使用一个方差为s 的高斯核的二阶导数生成一个探测核,用于测量导数方向范围内(-s,s)内外区域之间的对比度。

      但是同一图像中,线性结构的粗细肯定是不同的,同样的窗口大小是无法全部适用的,针对上面的问题,有人提出了多模板的方法。即对一个点用多种尺度的高斯模板进行卷积,然后选择各向异性最强的结果作为该点的输出。值得一试。

      回过头来看,根据海森矩阵,还可以确定一张图像中的角点部分,即前面表格中提到的两个特征值的绝对值都较大的情况。其实这就是Harris 角点检测的主要思想。

最后

      这篇博客准备了一段时间了,主要是公式的原因,我想使用markdown编辑器以及LaTex公式编辑器编辑公式,但是没搞定,还需要琢磨一下,等我琢磨好了,就把这篇文章中的图片格式的公式换掉。

 

参考文献:

     Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A. (1998) Multiscale vessel enhancement filtering. In: Wells W.M., Colchester A., Delp S. (eds) Medical Image Computing and Computer-Assisted Intervention — MICCAI’98. MICCAI 1998. Lecture Notes in Computer Science, vol 1496. Springer, Berlin, Heidelberg

      

 

2015-05-27 13:13:26 lxw907304340 阅读数 2336


图像掩膜:

用选定的图像、图形或物体,对处理的图像(全部或局部)进行遮挡,来控制图像处理的区域或处理过程。用于覆盖的特定图像或物体称为掩模或模板。光学图像处理中,掩模可以是胶片、滤光片等。数字图像处理中,掩模为二维矩阵数组,有时也用多值图像。数字图像处理中,图像掩模主要用于:①提取感兴趣区,用预先制作的感兴趣区掩模与待处理图像相乘,得到感兴趣区图像,感兴趣区内图像值保持不变,而区外图像值都为0。②屏蔽作用,用掩模对图像上某些区域作屏蔽,使其不参加处理或不参加处理参数的计算,或仅对屏蔽区作处理或统计。③结构特征提取,用相似性变量或图像匹配方法检测和提取图像中与掩模相似的结构特征。④特殊形状图像的制作。

 

图像处理(image processing)

用计算机对图像进行分析,以达到所需结果的技术。又称影像处理。图像处理一般指数字图像处理。数字图像是指用工业相机、摄像机、扫描仪等设备经过拍摄得到的一个大的二维数组,该数组的元素称为像素,其值称为灰度值。图像处理技术的一般包括图像压缩,增强和复原,匹配、描述和识别3个部分。 常见的系统有康耐视系统、图智能系统等,目前是正在逐渐兴起的技术。

21世纪是一个充满信息的时代,图像作为人类感知世界的视觉基础,是人类获取信息、表达信息和传递信息的重要手段。数字图像处理[9],即用计算机对图像进行处理,其发展历史并不长。数字图像处理技术源于20世纪20年代,当时通过海底电缆从英国伦敦到美国纽约传输了一幅照片,采用了数字压缩技术。首先数字图像处理技术可以帮助人们更客观、准确地认识世界,人的视觉系统可以帮助人类从外界获取3/4以上的信息,而图像、图形又是所有视觉信息的载体,尽管人眼的鉴别力很高,可以识别上千种颜色,但很多情况下,图像对于人眼来说是模糊的甚至是不可见的,通过图象增强技术,可以使模糊甚至不可见的图像变得清晰明亮。

在计算机中,按照颜色和灰度的多少可以将图像分为二值图像、灰度图像、索引图像和真彩色RGB图像四种基本类型。大多数图像处理软件都支持这四种类型的图像。

 

1.二值图像

一幅二值图像的二维矩阵仅由01两个值构成,“0”代表黑色,“1”代白色。由于每一像素(矩阵中每一元素)取值仅有01两种可能,所以计算机中二值图像的数据类型通常为1个二进制位。二值图像通常用于文字、线条图的扫描识别(OCR)和掩膜图像的存储。

 

2.灰度图像

灰度图像矩阵元素的取值范围通常为[0255]。因此其数据类型一般为8位无符号整数的(int8),这就是人们经常提到的256灰度图像。“0”表示纯黑色,“255”表示纯白色,中间的数字从小到大表示由黑到白的过渡色。在某些软件中,灰度图像也可以用双精度数据类型(double)表示,像素的值域为[01]0代表黑色,1代表白色,01之间的小数表示不同的灰度等级。二值图像可以看成是灰度图像的一个特例。

 

3.索引图像

索引图像的文件结构比较复杂,除了存放图像的二维矩阵外,还包括一个称之为颜色索引矩阵MAP的二维数组。MAP的大小由存放图像的矩阵元素值域决定,如矩阵元素值域为[0255],则MAP矩阵的大小为2563,用MAP=[RGB]表示。MAP中每一行的三个元素分别指定该行对应颜色的红、绿、蓝单色值,MAP中每一行对应图像矩阵像素的一个灰度值,如某一像素的灰度值为64,则该像素就与MAP中的第64行建立了映射关系,该像素在屏幕上的实际颜色由第64行的[RGB]组合决定。也就是说,图像在屏幕上显示时,每一像素的颜色由存放在矩阵中该像素的灰度值作为索引通过检索颜色索引矩阵MAP得到。索引图像的数据类型一般为8位无符号整形(int8),相应索引矩阵MAP的大小为2563,因此一般索引图像只能同时显示256种颜色,但通过改变索引矩阵,颜色的类型可以调整。索引图像的数据类型也可采用双精度浮点型(double)。索引图像一般用于存放色彩要求比较简单的图像,如Windows中色彩构成比较简单的壁纸多采用索引图像存放,如果图像的色彩比较复杂,就要用到RGB真彩色图像。

 

4.RGB彩色图像

RGB图像与索引图像一样都可以用来表示彩色图像。与索引图像一样,它分别用红(R)、绿(G)、蓝(B)三原色的组合来表示每个像素的颜色。但与索引图像不同的是,RGB图像每一个像素的颜色值(由RGB三原色表示)直接存放在图像矩阵中,由于每一像素的颜色需由RGB三个分量来表示,MN分别表示图像的行列数,三个M x N的二维矩阵分别表示各个像素的RGB三个颜色分量。RGB图像的数据类型一般为8位无符号整形,通常用于表示和存放真彩色图像,当然也可以存放灰度图像。

 

5.数字化图像数据有两种存储方式:<1>位图存储(Bitmap)<2>矢量存储(Vector)

我们平常是以图像分辨率(即像素点)和颜色数来描述数字图象的。例如一张分辨率为640*480,16位色的数字图片,就由2^16=65536种颜色的307200(=640*480)个素点组成。

<1>位图图像:

位图方式是将图像的每一个象素点转换为一个数据,当图像是单色(只有黑白二色)时,8个象素点的数据只占据一个字节(一个字节就是8个二进制数,1个二进制数存放象素点);16色(区别于前段“16位色”)的图像每两个象素点用一个字节存储;256色图像每一个象素点用一个字节存储。这样就能够精确地描述各种不同颜色模式的图像图面。位图图像弥补了矢量式图像的缺陷,它能够制作出色彩和色调变化丰富的图像,可以逼真地表现自然界的景象,同时也可以很容易地在不同软件之间交换文件,这就是位图图像的优点;而其缺点则是它无法制作真正的3D图像,并且图像缩放和旋转时会产生失真的现象,同时文件较大,对内存和硬盘空间容量的需求也较高。位图方式就是将图像的每一像素点转换为一个数据。如果用1位数据来记录,那么它只能代表2种颜色(2^1=2);如果以8位来记录,便可以表现出256种颜色或色调(2^8=256),因此使用的位元素越多所能表现的色彩也越多。通常我们使用的颜色有16色、256色、增强16位和真彩色24位。一般所说的真彩色是指24(2^24)的位图存储模式适合于内容复杂的图像和真实照片。但随着分辨率以及颜色数的提高,图像所占用的磁盘空间也就相当大;另外由于在放大图像的过程中,其图像势必要变得模糊而失真,放大后的图像像素点实际上变成了像素“方格”。 用数码相机和扫描仪获取的图像都属于位图。

 

<2>矢量图像:

矢量图像存储的是图像信息的轮廓部分,而不是图像的每一个象素点。例如,一个圆形图案只要存储圆心的坐标位置和半径长度,以及圆的边线和内部的颜色即可。该存储方式的缺点是经常耗费大量的时间做一些复杂的分析演算工作,图像的显示速度较慢;但图像缩放不会失真;图像的存储空间也要小得多。所以,矢量图比较适合存储各种图表和工程

 

6数字化

 

通过取样和量化过程将一个以自然形式存在的图像变换为适合计算机处理的数字形式。图像在计算机内部被表示为一个数字矩阵,矩阵中每一元素称为像素。图像数字化需要专门的设备,常见的有各种电子的和光学的扫描设备,还有机电扫描设备和手工操作的数字化仪。

 

7.图像编码

 

对图像信息编码,以满足传输和存储的要求。编码能压缩图像的信息量,但图像质量几乎不变。为此,可以采用模拟处理技术,再通过模-数转换得到编码,不过多数是采用数字编码技术。编码方法有对图像逐点进行加工的方法,也有对图像施加某种变换或基于区域、特征进行编码的方法。脉码调制、微分脉码调制、预测码和各种变换都是常用的编码技术。

 

8.图像压缩

 

由数字化得到的一幅图像的数据量十分巨大,一幅典型的数字图像通常由500×5001000×1000个像素组成。如果是动态图像,是其数据量更大。因此图像压缩对于图像的存储和传输都十分必要。

有两类压缩算法,即无损压缩和有损压缩。最常用的无损压缩算法取空间或时间上相邻像素值的差,再进行编码。游程码就是这类压缩码的例子。有损压缩算法大都采用图像交换的途径,例如对图像进行快速傅里叶变换或离散的余弦变换。著名的、已作为图像压缩国际标准的JPEGMPEG均属于有损压缩算法。前者用于静态图像,后者用于动态图像。它们已由芯片实现。

 

9.增强复原

 

图像增强的目标是改进图片的质量,例如增加对比度,去掉模糊和噪声,修正几何畸变等;图像复原是在假定已知模糊或噪声的模型时,试图估计原图像的一种技术。

图像增强按所用方法可分成频率域法和空间域法。前者把图像看成一种二维信号,对其进行基于二维傅里叶变换的信号增强。采用低通滤波(即只让低频信号通过)法,可去掉图中的噪声;采用高通滤波法,则可增强边缘等高频信号,使模糊的图片变得清晰。具有代表性的空间域算法有局部求平均值法和中值滤波(取局部邻域中的中间像素值)法等,它们可用于去除或减弱噪声。

早期的数字图像复原亦来自频率域的概念。现代采取的是一种代数的方法,即通过解一个大的方程组来复原理想的图片。

以提高图像质量为目的的图像增强和复原对于一些难以得到的图片或者在拍摄条件十分恶劣情况下得到的图片都有广泛的应用。例如从太空中拍摄到的地球或其他星球的照片,用电子显微镜或X光拍摄的生物医疗图片等。

图像增强:

使图像清晰或将其转换为更适合人或机器分析的形式。与图像复原不同,图像增强并不要求忠实地反映原始图像。相反,含有某种失真(例如突出轮廓线)的图像可能比无失真的原始图像更为清晰。常用的图像增强方法有:①灰度等级直方图处理:使加工后的图像在某一灰度范围内有更好的对比度;②干扰抑制:通过低通滤波、多图像平均、施行某类空间域算子等处理,抑制叠加在图像上的随机性干扰;③边缘锐化:通过高通滤波、差分运算或某种变换,使图形的轮廓线增强;④伪彩色处理:将黑白图像转换为彩色图像,从而使人们易于分析和检测图像包含的信息。

图像复原:

除去或减少在获得图像过程中因各种原因产生的退化。这类原因可能是光学系统的像差或离焦、摄像系统与被摄物之间的相对运动、电子或光学系统的噪声和介于摄像系统与被摄像物间的大气湍流等。图像复原常用二种方法。当不知道图像本身的性质时,可以建立退化源的数学模型,然后施行复原算法除去或减少退化源的影响。当有了关于图像本身的先验知识时,可以建立原始图像的模型,然后在观测到的退化图像中通过检测原始图像而复原图像。

图像分割将图像划分为一些互不重叠的区域,每一区域是像素的一个连续集。通常采用把像素分入特定区域的区域法和寻求区域之间边界的境界法。区域法根据被分割对象与背景的对比度进行阈值运算,将对象从背景中分割出来。有时用固定的阈值不能得到满意的分割,可根据局部的对比度调整阈值,这称为自适应阈值。境界法利用各种边缘检测技术,即根据图像边缘处具有很大的梯度值进行检测。这两种方法都可以利用图像的纹理特性实现图像分割。

 

10.形态学

 

形态学一词通常指生物学的一个分支,它用于处理动物和植物的形状和结构。在数学形态学的语境中也使用该词来作为提取图像分量的一种工具,这些分量在表示和描述区域形状(如边界,骨骼和凸壳)时是很有用的。此外,我们还很关注用于预处理和后处理的形态学技术,如形态学滤波、细化和裁剪。

数学形态学的基本运算:

数学形态学的基本运算有4个:腐蚀、膨胀、开启和闭合。数学形态学方法利用一个称作结构元素的”探针”收集图像的信息,当探针在图像中不断移动时,便可考察图像各个部分之间的相互关系,从而了解图像的结构特征。在连续空间中,灰度图像的腐蚀、膨胀、开启和闭合运算分别表述如下。

腐蚀

腐蚀“收缩”或“细化”二值图像中的对象。收缩的方式和程度由一个结构元素控制。数学上,AB腐蚀,记为AΘB,定义为:

换言之,AB腐蚀是所有结构元素的原点位置的集合,其中平移的BA的背景并不叠加。

膨胀

膨胀是在二值图像中“加长”或“变粗”的操作。这种特殊的方式和变粗的程度由一个称为结构元素的集合控制。结构元素通常用01的矩阵表示。数学上,膨胀定义为集合运算。AB膨胀,记为AB,定义为:

其中,Φ为空集,B为结构元素。总之,AB膨胀是所有结构元素原点位置组成的集合,其中映射并平移后的B至少与A的某些部分重叠。这种在膨胀过程中对结构元素的平移类似于空间卷积。

膨胀满足交换律,即AB=BA。在图像处理中,我们习惯令AB的第一个操作数为图像,而第二个操作数为结构元素,结构元素往往比图像小得多。

膨胀满足结合律,即A(BC)=(AB)C。假设一个结构元素B可以表示为两个结构元素B1B2的膨胀,即B=B1B2,则AB=A(B1B2)=(AB1)B2,换言之,用B膨胀A等同于用B1先膨胀A,再用B2膨胀前面的结果。我们称B能够分解成B1B2两个结构元素。结合律很重要,因为计算膨胀所需要的时间正比于结构元素中的非零像素的个数。通过结合律,分解结构元素,然后再分别用子结构元素进行膨胀操作往往会实现很客观的速度的增长。

 

开启

AB的形态学开

运算可以记做A?B,这种运算是AB腐蚀后再用B来膨胀腐蚀结果,即:

开运算的数学公式为:

其中

,∪{·}指大括号中所有集合的并集。该公式的简单几何解释为:A?BBA内完全匹配的平移的并集。形态学开运算完全删除了不能包含结构元素的对象区域,平滑了对象的轮廓,断开了狭窄的连接,去掉了细小的突出部分。

 

闭合

AB形态学闭运算记做A·B,它是先膨胀后腐蚀的结果:

从几何学

上讲,A·B是所有不与A重叠的B的平移的并集。想开运算一样,形态学闭运算会平滑对象的轮廓。然后,与开运算不同的是,闭运算一般会将狭窄的缺口连接起来形成细长的弯口,并填充比结构元素小的洞。

基于这些基本运算可以推导和组合成各种数学形态学实用算法,用它们可以进行图像形状和结构的分析及处理,包括图像分割、特征提取、边界检测、图像降噪、图像增强和恢复等。

 

11.图像分析

 

从图像中抽取某些有用的度量、数据或信息。目的是得到某种数值结果,而不是产生另一个图像。图像分析的内容和模式识别、人工智能的研究领域有交叉,但图像分析与典型的模式识别有所区别。图像分析不限于把图像中的特定区域按固定数目的类别加以分类,它主要是提供关于被分析图像的一种描述。为此,既要利用模式识别技术,又要利用关于图像内容的知识库,即人工智能中关于知识表达方面的内容。图像分析需要用图像分割方法抽取出图像的特征,然后对图像进行符号化的描述。这种描述不仅能对图像中是否存在某一特定对象作出回答,还能对图像内容作出详细描述。

图像处理的各个内容是互相有联系的。一个实用的图像处理系统往往结合应用几种图像处理技术才能得到所需要的结果。图像数字化是将一个图像变换为适合计算机处理的形式的第一步。图像编码技术可用以传输和存储图像。图像增强和复原可以是图像处理的最后目的,也可以是为进一步的处理作准备。通过图像分割得出的图像特征可以作为最后结果,也可以作为下一步图像分析的基础。

图像匹配、描述和识别对图像进行比较和配准,通过分制提取图像的特征及相互关系,得到图像符号化的描述,再把它同模型比较,以确定其分类。图像匹配试图建立两张图片之间的几何对应关系,度量其类似或不同的程度。匹配用于图片之间或图片与地图之间的配准,例如检测不同时间所拍图片之间景物的变化,找出运动物体的轨迹。

从图像中抽取某些有用的度量、数据或信息称为图像分析。图像分析的基本步骤是把图像分割成一些互不重叠的区域,每一区域是像素的一个连续集,度量它们的性质和关系,最后把得到的图像关系结构和描述景物分类的模型进行比较,以确定其类型。识别或分类的基础是图像的相似度。一种简单的相似度可用区域特征空间中的距离来定义。另一种基于像素值的相似度量是图像函数的相关性。最后一种定义在关系结构上的相似度称为结构相似度。

以图片分析和理解为目的的分割、描述和识别将用于各种自动化的系统,如字符和图形识别、用机器人进行产品的装配和检验、自动军事目标识别和跟踪、指纹识别、X光照片和血样的自动处理等。在这类应用中,往往需综合应用模式识别和计算机视觉等技术,图像处理更多的是作为前置处理而出现的。

多媒体应用的掀起,对图像压缩技术的应用起了很大的推动作用。图像,包括录像带一类动态图像将转为数字图像,并和文字、声音、图形一起存储在计算机内,显示在计算机的屏幕上。它的应用将扩展到教育、培训和娱乐等新的领域。

 

 

MATLAB图像处理

阅读数 5485

没有更多推荐了,返回首页