图像变换的不变性和偏微分方程应用

1、图像变换的不变性和偏微分方程应用讨论具有单调和对比不变的图像变换，即形态学算子(基本思想是用具有一定形态的结构元素去量度和提取图像中的对应形状以达到对图像的分析和识别的目的，其基本元算有四个：膨胀、腐蚀、开启、闭合)，重点介绍：膨胀、腐蚀和中值算子。然后依次赋予形态学算子平移不变性、欧氏不变性和仿射不变性，并讨论它们的微分性质，导出相关的偏微分方程。6.1 形态学算子单调和对比不变的图像变换6.1.1 定义前面学过连续模型下图像空间的定义，是一族由R2R的特殊函数组成的函数空间，并记为F 。图像变换T是作用在F 上的一个算子，即T将一副图像u变换为另一幅图像Tu。图像水平集之间的变换，是对于F

2、 中所有函数，Y表示在F 所拥有的所有水平集，即 Y = clu ; uF , l0, 1 这是一个由R2的子集组成的集合族。对于图像变换T，引进算子T作用在Y上，它将一个水平集X转换为另一个水平集TX，即 T: X Y TX Y定义1：称图像变换T是单调递增的，如果对于任意两两幅图像u,vF u v Tu Tv集合算子T是单调递增的，如果对于任意X,YY X Y T(X) T(Y)定义2：图像变换T是对比不变的，如果对每一个连续对比变换g，对任意的uF ，都满足g(u)F 和 g(Tu) = T( g(u) )同时满足单调性和对比不变的图像变换被称为形态学算子。可以证明：线性算子是单调的，但

3、不是对比不变的。例1：最大值滤波是对比不变的。最大值滤波定义：其中B是包含原点的闭集，x+B = x+z;zB 。假设 ，由于x+B为闭集，zx+B，满足u(z) = a，而 u(y) u(z) , yx+B又因为对比变换g是单调递增的，所以 g( u(y) ) g( u(z) ) = g(a) yx+B即，对图像g(u)满足对比不变定义 D( gu(x) ) = g(a) = g( Du(x) )对比不变的图像变换有一特殊性质，即变换的结果使图像保留了原图像的部分灰度。一副二值图像在经过对比不变图像变换后还是一副二值图像。但线性滤波器都不具备这一特性。下面定理说明了这一性质，记R(u)为图像

4、u的值域，即 R(u) = s0,1, x, u(x) = s 其中Ru是包含R(u)的最小闭集。定理1：T是一对比不变的图像变换。那么对每一副图像u，R(Tu) Ru，特别的，如果图像u只有有限个灰度值，则Tu只取其中的部分灰度值。证明：考虑一连续单调递增函数g，满足g(s) = s，当 sRu时。否则，g(s) s。定义： g(s) = s + d(s,Ru)2其中d(s,X)表示s到X距离 。当且仅当 sRu时，有 d(s,Ru) = 0因此，当且仅当sRu时，g(s) = s，所以 g(u) = u。因为T是对比不变的，所以 Tu = T( g(u) ) = g(Tu)因此 (Tu)(

5、x)Ru 。 定义3：一个图像变换T是灰度平移不变的，如果对任意的常数C，有 T(u+C) = Tu + C如果图像变换T同时具有灰度平移不变性和对比不变性，就得到下面的结论。定理2：T是一个单调灰度平移不变算子，如果u(x)是R2上的Lipschitz函数，那么Tu(x)也是Lipschitz函数，并且Tu(x)的Lipschitz常数比u(x)的Lipschitz常数小。Lipschitz常数定义：如果函数u满足 |u(x) u(y)| k |x-y| , x, y则u为Lipschitz函数，k为u的Lipschitz常数。证明：假设u的Lipschitz常数为K。对任意的x,y,z有

6、|u(x+z)- u(y+z)| K|x-y| u(y+z)- K|x-y|u(x+z)u(y+z)+ K|x-y|因为T单调，考虑上面关于z的函数，有 T(u(y+z)- K|x-y|)Tu(x+z)T(u(y+z)+ K|x-y|)注意到取 z=0, 有 T(u(y+z) = (Tu)(y)，用T的灰度平移不变性（将K|x-y|看做C ）得 Tu(y)- K|x-y|Tu(x)Tu(y)+ K|x-y|。 |Tu(x)- Tu(y)| K|x-y| 6.1.2 从形态学算子到集合算子记集合XW上的特征函数为1x，即1x也被认为是一个图像函数，即1xF 。借助特征函数，可从单调、对比不变的图

7、像变换（形态学算子）T衍生出一个集合变换T。定义4：令T是一个单调、对比不变的图像变换，定义T的伴随集合算子T为 X W, 1X F T(X) = c1( T(1X) )另外 T(F) = F, T(W) = W如果T作为函数是单调的，那么T作为集合变换也是单调的。因为 X Y 1X 1YT作为单调的图像变换，使单调性得以保持 T(1X) T(1Y)定理3：T是一个对比不变的单调算子。阈值函数gl(s)定义为, 如果sl，则 gl(s) = 1 ；否则gl(s) = 0 。那么T几乎处处和每一个阈值函数相交换，即 gl (Tu) = T( gl (u) )对l,x几乎处处成立。证明：定义则 g

8、el (s) 是对比变换（连续、单调递增的），且gel(s)gl，于是同样的方法，用不减函数 g el gl，可证明 T( gl(u) ) g-l(Tu)其中，g-l(s)=1 ，当sl 时； g-l(s)=0 ，当 sl时。因此，有 gl(Tu) T( gl(u ) ) g-l(Tu)我们考虑可数因而可忽略的子集R，所有的l 满足 meas( x, Tu(x) = l ) 0 对于l R ，有 g-l(Tu) = gl(Tu)几乎处处成立。这样对几乎每一个l，会得到 T( gl(u) ) = gl(Tu)几乎处处成立。 定理4：T是定义在图像函数集合F 上的单调对比不变算子，1xF 。T的伴

9、随集合算子为T，则T是单调的，并且 uF 有 T( clu) = cl(T(u) *对l, c几乎处处成立。并且 Tu(x) = sup l,xT(clu) *对x几乎处处成立。另外， T(F) = F, T(W) = W几乎处处成立。*式说明图像变换后的水平集是原图像水平集（并且是同一个 l）在伴随集合算子作用下的结果。T(F) = F 说明当l=1时，*式成立；T(W) = W 说明当l=0时，*式成立。这里涉及到水平集和最大值表示公式： cl(u) = xW ; u(x) l u(x) = sup l; x cl(u) 证明：根据T的定义，显然有 1clu = gl(u) c1( gl(

10、v) ) = clv并且T和gl几乎处处可交换。由定理3，得到： T(clu) = c1( T(1clu ) ) = c1( T( gl(u) ) ) = c1( gl(Tu) ) = cl(Tu) 对于x, l0 几乎处处成立。由 T(clu) = cl(Tu) 知 T(clu)是Tu的水平集。那么显然*式成立。令u是一个常函数0，对于l0 ，有clu= F 。 利用*式，有 cl(Tu) = T(clu ) = T(F ) 对 l0几乎处处成立。而且，由于对比不变算子T和常函数0相交换，因此Tu = 0，并且对 l 0 ，cl(Tu) = F，则有T(F) = F几乎处处成立。同理可证 T

11、(W) = W。 定理说明：如果图像变换T是单调且对比不变的，那么计算Tu可以通过一下算法实现：（1）计算u的所有水平集 cl(u) ( l0,1 )；（2）对每一个水平集 cl(u) ，用T的伴随集合算子T作用，得到 T( cl(u) )；（3）用最大值表示公式得到Tu, 整个过程如下。这种算法适用于T难以实现，而T容易计算的情况。6.1.3 从集合算子到形态学算子考虑，给定一个单调的集合算子T是否可以得到一个对比不变的单调图像变换呢？自然的思路就是令 Tu(x) = sup l,xT(clu) 定理5：令T是一个Y Y单调算子，满足 T(F) = F , T(W) = W那么，可以定义图像

12、变换 Tu(x) = sup l ,xT(clu) 对于所有的l，满足 cl(Tu) = T(cl(u)则对几乎所有的 l R g(Tu) = T(g(u)证明：对每一个l ，我们有 cl(Tu) = T(cl(u)即对R中的所有l满足meas(R) = 0。注意到uv当且仅当 clu clv ，对R的一个稠密可数子集合上的所有l ，可得T是单调的 cl(Tu) = T(cl(u) T(cl(v) = cl(Tv) Tu Tv 下面证明，T和对比变换相交换。假设g是严格增加的，设 和 ，对于 lg(+) 有clg(u) = F，因此，T(clg(u) = F；对于 lg(-) 有 clg(u)

13、 = RN，因此 T(clg(u) = RN 。T( g(u(x) ) = sup l, g(-) lg(-) , x T(clg(u) = sup g(m), x T(cgg(u) = sup g(m), x T(cmu) = g( Tu(x) )下面验证T和一般的不减对比变换g相交换。严格增加连续函数gn和hn满足 gn(s) g(s)，hn(s) g(s)对所有的s和gnghn。因此由上面结论有 T(g(u)T(gn(u) = gn(Tu) g(Tu) T(g(u)T(hn(u) = hn(Tu) g(Tu)可以推出 T(g(u) = g(Tu)。 6.1.4 应用实例：“Extrema

14、 Killer”算子“Extrema Killer”算子是一个图像光滑算法，作用是去除图像中的“峰（peak）”孤立的水平集，尤其对椒盐噪声效果显著。算法如下：（1）假设一个集合X有若干连通区域组成定义一个集合变换 Tb(X) = Xb，而（2）Extrema Killer图像变换定义为 Tbu(x) = sup l ,xTb(clu) *式定义了集合算子是Extrema Killer变换的伴随集合算子，即 cl(Tbu) = Tb(clu)噪声图像 killer算子作用后图像 （改进的Extrema Killer）6.2 平移不变的形态学算子主要描述平移不变的形态学算子以及它的伴随集合算子，

15、也就是平移不变的单调集合算子。记平移为tx，且满足（1）对于集合X： txX = x+X = x+y | yX （2）对于图像u： (tx)u(y) = u(y-x)其中x不是一个二维的点，而是表示一个二维向量。定义5：集合算子T是平移不变的，如果 tx( T(X) ) = T( txX )定义图像变换是平移不变的，如果 tx( T(u) ) = T( txu)定理6：（Matheron）令T是平移不变的单调集合算子，那么存在一个集合族B B = X, 0 T(X) 其中0是R2中的原点，T满足其中X-y = X+(-y)。相反*式也定义了一个单调、平移不变的集合算子。证明：先看*式等价于利用

16、单调性和平移不变性，得下面的等价关系第五个等价性质成立理由是：如果 BX，并且BB那么XB ，因为BB ，即 B, 0T(B) 又因为BX，则有 X, 0T(X) XB就有 。相反的，如果算子通过*式定义，显然是单调和对比不变的。 定理7：F 是图像函数空间，Y是 F 中所有水平集的集合。假设对比不变下是稳定的，并且包含了T中所有元素的特征函数。令T是T的伴随集合算子。如果xR2，集合族B x = X, xT(X) 那么uF 有：对x几乎处处成立。其中B =B 0 X,0T(X) 另外，如果T是位移不变算子，则相反的，如果一个算子通过以上的公式定义，则该算子是单调和对比不变的。证明：令 其中B

17、 x = X, xT(X) x以下证明Tu(x) = Tu(x)几乎处处相等。选择一个可数的稠密y0,1，满足ly, clTu(x) = T(clu)(x)对xRNNl成立。（对比下稳定定理4）这里Nl的Lebesge测度为0。设N=Nl ，则N的Lebesge测度也为0。为证明定理，先证明对所有的ly和所有的 xRNNl，有：（即处处相等了） Tu(x) l Tu(x) l 对任意 l,my，我们有：第五个等价关系：因为如果BB 并且BB x则 XB x。 B, 0T(B) 、 B, xT(B) 因为BX则有 X, 0T(X) 、 X, xT(X) = B x那么，如果某些BB x,那么一定

18、有Bcmu ，就也有 cmuB x，于是证明了提出的命题。 布尔代数（Boolean Algebra）中有个著名的结论：如果T是一个sup inf形式的算子，那么T也具有inf sup的形式，即此时的B 与sup inf形式中的B 是不同的。定理8：如果T是平移不变的形态学算子，那么它的伴随算子集合T可以通过一下公式来定义，证明：T满足定理7，并且可以延拓到F 中函数的所有水平集。很容易由定理7推出定理6的结果。 对于X的特征函数1x，则 ,当x + B X ,当x + B X于是， 当且仅当BB 满足 x+BX和定理6。6.3 形态学算子：膨胀和腐蚀算子6.3.1 定义定义6：X是R2中一子

19、集，t0是一尺度参数。称Dt是基于一个集合B和尺度参数t的膨胀，如果其中集合B称为结构元素。同样，Et表示结构元素B的尺度参数t的腐蚀公式说明Dt和Et是平移不变的单调集合算子，满足定理6，此时B 中只含有一个元素tB。例2：（1）如果结构元素B=x0是仅仅包含一个点的集合，那么DtX = X+tx0是一个平移算子。相应的，EtX = X - tx0也是一个平移算子。（2）如果结构元素B是一个圆心在原点半径为1的开球D(0,1)，那么DtX就是X的t-领域，即与X的距离小于t的所有点的集合。(3)令B = D(0,1)，则定理9：（1）Dt(Xc) = (EtX)c其中Xc = R2/X；（2

20、）如果结构元素为D(0,1)，那么Dt，Et是旋转不变的，即它们和旋转运算可以交换；（3）Ds+t= DsDt（ Es+t= EsEt），当且仅当结构元素B是凸的。证明：（1）（3）首先证明 任取B, (t+s)B = tB+sB当且仅当B是凸的。因为tB + sB= sx+ty | x,yB ，所以ztB+sB, x, yB，满足而(s+t)x和(s+t)y都属于(s+t)B，又因为B是凸的，所以(t+s)B = tB+sB。反过来，如果(t+s)B = tB+sB成立，那么对于任意的x,yB,存在zB，满足(s+t)z = sx+ty，也就是所以B是凸的。注意到 DtDsX = (X+sB

21、) + tB = X + sB + tB Ds+tX = X + (s+t)B根据上面的结论，Ds+t = DsDt成立当且仅当它们的结构元素B是凸的。同样的结论也适用于腐蚀算子。 定义7：u是一副图像，称Dt是基于结构元素B和尺度参数t的膨胀变换，如果类似，基于结构元素B和尺度参数t的腐蚀算子Et被定义为上面的定义说明Dtu(x)具有inf sup的形式，即而B 只含有tB一个元素，所以Dt是一个平移不变的形态学算子，Et亦然。定理10：对于图像的膨胀和腐蚀算子，如果B是关于0对称的，那么 -Et(-u) = Dt(u)证明：第三个等号是由于B的对称性。原图 对黑色的膨胀（对背景色的腐蚀）

22、结构元素是圆盘 I = imread(star.bmp);subplot(1,2,1),imshow(I);J = I;w,h = size(I);r = 10;for i= r+1:w-r for j= r+1:h-r min = 256; for x = -r:r for y = -r:r if sqrt(x*x+y*y)r & I(i+x,j+y)min min = I(i+x,j+y); end end end J(i,j) = min; endendsubplot(1,2,2),imshow(J);I = imread(star.bmp);subplot(1,2,1),imshow(

23、I);J = I;w,h = size(I);r = 10;for i= r+1:w-r for j= r+1:h-r max = 0; for x = -r:r for y = -r:r if sqrt(x*x+y*y)max max = I(i+x,j+y); end end end J(i,j) = max; endendsubplot(1,2,2),imshow(J); 原图膨胀腐蚀I = imread(girl.bmp);subplot(1,3,1),imshow(I);se = strel(disk,4);D = imdilate(I,se);subplot(1,3,2),imsh

24、ow(D);E = imerode(I,se);subplot(1,3,3),imshow(E);从前面膨胀和腐蚀的结果图像可以看出，单独的膨胀和腐蚀都不可能成为一个好的滤波器，因为膨胀是图像变亮，从而增加了许多白色区域（腐蚀增加了黑色区域）。但经过一些组合可以产生很好的滤波效果。 开运算：先对图像进行腐蚀然后再膨胀； 闭运算：先对图像进行膨胀然后再腐蚀。例如： T = Dt 。Et 。Et 。Dt T= Et 。Dt 。Dt 。Et原图 腐蚀 膨胀、腐蚀膨胀、膨胀、腐蚀 腐蚀、膨胀、膨胀、腐蚀原图膨胀 腐蚀、膨胀腐蚀、腐蚀、膨胀 膨胀、腐蚀、腐蚀、膨胀集合的膨胀和腐蚀算子与图像的膨胀和腐蚀变

25、换存在以下关系：（1）令图像膨胀变换的结构元素为 -tb，bB ，（2）令集合膨胀算子的结构元素为B，由于此时 Dt(clu) = cl(Dtu)所以，上述的集合膨胀算子是图像膨胀变换的伴随集合算子。所以有： Dt(clu) = cl(Dtu)6.3.2 偏微分方程和膨胀（腐蚀）算子 记 其中 表示内积，当B = D(0,1)时，| |B就是Euchlid范数。定理11：（Lax formula）如果u(t,x) = Dtu0(x)，并且结构元素B是凸的，那么u(t,x)满足 ut = | Du |-B其中u对x两次可微。对应的，如果u(t,x) = Etu0(x)，那么u(t,x)满足 ut

26、 = - | Du |-B其中u对x两次可微。证明：先看在t=0时的性质。假设u0在x处是C2的，已知 u(t,x) = Dtu0(x)，u(0,x) = u0(x)所以既然u0在x处是可微的，那么两边同除以h，并令h0，得到（t=0时成立）下证对于任意尺度t时具有同样的关系：由于B是凸的，因此Dt+h= DtDh= DhDt，所以 u(t+h,x) - u(t,x) = Dh(u(t)(x) - u(t)(x)两边同除以h，并令h0，由上面的结果且用u(t)代替u0，就得到一般的结果。 6.3.3 膨胀和腐蚀算子的离散算法膨胀与腐蚀算子具有类似的性质，下面只给出膨胀算子的伪代码。已知图像u的

27、分辨率为mn，ui,j（i=1,m,j=1n）表示在i,j处的灰度值，取结构元素为B，t为参数，Vi,j表示运算的结果：Dilation(u,B,t)FOR x0 = 1 to m FOR y0 = 1 to n FOREACH(tB中的像素点p) 得到p的坐标(px,py) 将ux0+px,y0+py的灰度值存入数组array中 求array中的最大值max_array 令Vx0,y0= max_array END FOREACH ENDEND6.4 形态学算子：中值算子6.4.1 定义首先假设图像u的所有水平集都是Lebesgue可测的。权函数k(y)满足集合B的k-测度为定义8：令X是R

28、2的一个可测子集，k是权函数。称X的中值集（k加权）并用medkX来表为 medkX= x, |X-x|k1/2 定理12：算子 medk:YY是单调的算子并且满足集合连续性质，即，如果 (Xl)lR是一个递减的可测集合，并且满足 Xl = mlXm, lR那么 medk(Xl) = mlmedk(Xm)证明：由medk的定义知其是单调的， XY X-xY-x 所以 medk(Xl) mlmedk(Xm)反过来，令 x mlmedk(Xm)由medk的定义可得 m , |Xm - x|k 12因为Xm是一个单调减的集合序列，同时具有有限的测度，根据实变函数中的Lebesgue极限收敛定理 |X

29、m - x|k |Xl - x|k并且|Xl - x|k1/2，再由medk的定义，有xmedkXl即 medk(Xl)ml medk(Xm) 其中 可测，且其中 可测，且则有若定义9：图像的加权中值滤波器基于一个结构元素B，其定义为显然，medk(u)是个sup inf型的，是一个平移不变的形态学算子。如果在记号上不区分图像中值滤波器和集合中值算子，都记为medk。 medk作为一个单调的集合算子，可以用最大值表现公式扩展到一个图像的变换T（定理4） Tu(x) = sup u(x) | xmedk(X) 同时满足 cl(Tu) = medk(clu)由于medk作为一个单调、平移不变的集合

30、算子，所以T是一个单调的、对比不变的变换，根据定理7其中B = B,0medk(B) 。所以有：定理13：medk是图像加权中值滤波器的伴随集合算子。6.4.2 中值算子的离散算法已知图像u的分辨率为mn，ui,j（i=1,m,j=1n）表示在i,j处的灰度值，取结构元素为B，t为参数，Vi,j表示运算的结果：Median(u,B,t)NB = 结构元素B的像素数目FOR x0 = 1 to m FOR y0 = 1 to n FOREACH(B中的像素点p) 得到p的坐标(px,py) 将ux0+px,y0+py的灰度值存入数组array中 对array进行排序，排在第（NB/2+1）的值m

31、ed_array 令vx0,y0=med_array END FOREACH ENDEND噪声图像及32：32：224水平线中值滤波后图像及32：32：224水平线Matlab源码：Image = imread(lena.bmp);subplot(2,3,1),imshow(Image);I = imnoise(Image,salt & pepper,0.05);subplot(2,3,2),imshow(I);L_s_I = Level_Set_Line(I,32,32,256);subplot(2,3,3),imshow(uint8(L_s_I);M = medfilt2(I,3,3);s

32、ubplot(2,3,5),imshow(M);L_s = Level_Set_Line(M,32,32,256);subplot(2,3,6),imshow(uint8(L_s);6.5 欧氏不变的形态学算子6.5.1 定义和微分性质 记 Ha = Tx1+ax22(0)其中Tx1+ax22表示T(u)，u = x1+ax22。如果T是单调的，那么H也是单调的。定义 Thx = hTx = hH0 Thx1+ax22(0) = hTx1+hax22(0) = hHah用D(0,M)表示圆心在0，半径为M的圆。定理14：令B 是一族R2有界子集族（ BB ，BD(0,M) ）并且是各向同性的（

33、 BB ，RBD(0,M)，R是一个旋转变换 ）。令和相关的带参数h的变换那么，对于任意C2的函数u，有 (Thu)(x) = u(x) + hTx(0)| Du |(x) + O(h2)证明：由于T( u-u(x) ) = Tu - u(x)，不失一般性，令u(x)=0。因为T是平移不变和旋转不变的，选择图像支撑集W的左下角定位在x，两个轴的方向(i,j)定义为所以，对y = (y1, y2)，Talyor展开有 u(y) = py1+ O(|y|2)其中p = |Du|(0)0，这是因为 ，u=x1+ax22。从而可以有，如果y = (y1,y2)hD(0,M)，那么 py1 - O(h2

34、) u(y) py1+ O(h2) |y|2 |(hM)2| O(|y2|)O(|h2|)M2 O(h2)根据Th的单调性，又因为BB ,BD(0,M)，故 Thu(0) = Thpy1(0) + O(h2)再由前面的定义 Thu(0) = Thpy1(0) + O(h2) = ph( Ty(0) ) + O(h2)利用p = | Du(0) |，u(x) = 0可得结论。 定理15：令B 是一族R2有界子集族（ BB ，BD(0,M) ）并且是各向同性的（ BB ，RBD(0,M)，R是一个保距变换 ）。令和相关的带参数t的变换设H(0)=Tx(0)=0，则对R2上的每一个C3函数u，有（i

35、）对每一个紧集K x, Du(x)0 Thu(x) = u(x) + h|Du(x)| H(1/2hcurv(u) + Ox(h3)其中D(u)表示微分算子，|Ox(h3)|Ckh3 ,常数Ck依赖于u和K，且（ii）对每一个紧集 KR2且 K x,Du(x) = 0 有 | Thu(x) u(x) | M2 h2 | D2u(x) | + Ox(h3)其中 Ox(h3) Ckh3,常数Ck依赖于u和K。（i）结论，刻画各向同性的形态学算子的微分性质，和第五章定理1类似，对Th有下面的结论：令h0，n，并且nh2 = t，那么是下面偏微分方程初值问题在t时刻的解其中g是一个单调函数。6.5.2

36、 中值滤波的微分性质定义9中的结构元素B = D(0,h)，k测度是B上的平均测度，即 k = lD(o,h) ph2所以这里选择了特殊的结构元素B，所以此时med是个欧氏不变的形态学算子。下面讨论一下它的性质：定理16：令u是R2上的一个C3函数，则(i) medD(0,h)u(x) = u(x) + 1/6 curv(u) |Du|(x) h2 + O(h3)其中O(h3) Ckh3，对 x,Du(x)0 上的每一个紧子集；(ii) | medD(0,h)u(x) u(x) | |D2u(x)| h2 + O(h3)其中O(h3)Ckh3，对R2上的每一个紧子集。定理17： medD(0,

37、1)y1 + hy22(0) = h3 + O(h3)证明：当h比较小时，抛物线p(h)把D(0,1)分割成两块相交的区域，其中p(h)的方程是 x+hy2 = m。0yxy=1y=-1m(h)D(0,1)P(h)ABCDE由抛物线和OY轴组成的区域，其代数面积为（OY轴左侧面积为负，右侧面积为正）抛物线顶点位置为(m(h),0)决定这两块区域的面积。如果抛物线p(h)把D(0,1)分割为两块面积相等的区域当且仅当 2m(h) 2h3 = 2area(ABE)不计算area(ABE)的值，只需给出一个估计所以 |2m(h) - 2h3| |m(h) - h|3最后可以得到 m(h) = h3

38、+ O(h3) 对于x + hy2 = m截圆面积一半，其实是h0，此时m(h)0令t = m-h3，则|2t| (t - 2h3)3 |2th| (th - 23)(th - 23)2h2h0右边0，即th02th 8h227 即 th3 854 得到t与O(h3)同阶。即得 m-h3 = O(h3)定理16的证明：根据定理17，与medD(0,h)相关的H(0) = 0。所以可以使用定理13，由定理17得到 H(h) = h3 + O(h3)所以（i）（ii）可以由定理15的（i）（ii）得到。 6.5.3 中值曲率驱动（mean curvatrue motion）方程由定理16，令h0，

39、n，并且nh2 = t，那么是下面偏微分方程初值问题在t时刻的解：当C=1时，就是中值曲率驱动（MCM）方程。原图 MCM步长迭代2次步长迭代5次 10次 均值滤波步长迭代5次 步长迭代2次I = imread(other.bmp);subplot(1,2,1),imshow(I);w,h = size(I);I = double(I);t = 0.1; N = 10;for k = 1:N Dx = zeros(w,h);Dy = zeros(w,h); Dxx = zeros(w,h);Dyy = zeros(w,h);Dxy = zeros(w,h); for i = 2:w-1 for

40、 j = 2:h-1 Dx(i,j)= I(i+1,j)-I(i,j); Dy(i,j)= I(i,j+1)-I(i,j); Dxx(i,j)= -2*I(i,j)+I(i-1,j)+I(i+1,j); Dyy(i,j)= -2*I(i,j)+I(i,j-1)+I(i,j+1); Dxy(i,j)= Dx(i,j+1)-Dx(i,j); end end for i = 1:w for j = 1:h if Dx(i,j)*Dx(i,j)+Dy(i,j)*Dy(i,j)= 0 I(i,j) = I(i,j) +t*(Dxx(i,j)*Dy(i,j)*Dy(i,j)-2*Dxy(i,j)*Dx(i,j)*Dy(i,j)+Dyy(i,j)*Dx(i,j)*Dx(i,j)/(Dx(i,j)*Dx(i,j)+Dy(i,j)*Dy(i,j); end end endendsubplot(1,2,2), imshow(uint8(I);从试验