多尺度tophat的实现问题

各位好，最近在用tophat作图像增强，看了下原理，还蛮简单的，但是在实现多尺度tophat上面遇到些问题，还望高手指点一下。

单尺度tophat和多尺度tophat我都做了，但是可能因为我的实现上面有些问题，导致我的多尺度tophat的增强效果几乎和单尺度的没多大差别，有时候甚至还不如单尺度的，而且运行时间上明显增加，我参考的文献是【Peihe Tang, Hao Liu, Kaiqiong Sun. Enhancement of coronary angiogram by estimation of local background】，附件是我的代码，顺便贴出来一下。

这个我已经折腾蛮久了，真心希望有高手指点一下哈，谢谢了。

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

% Single scale enhancement for pre-observation.

se = strel('disk', 11);

res = imtophat(Iori,se); res=imadjust(res,stretchlim(res),[0 1]);

figure;

subplot(1,3,1); imshow(Iori); title('Original Image');

subplot(1,3,2); imshow(Open_res); title('Background Estimation using Single Scale');

subplot(1,3,3); imshow(res); title('Top-hat Result using Single Scale');

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

% Multiscale enhancement.

[nrow,ncol]=size(Iori);

% Construct different Structure Elements(SE).

se1=strel('disk',1); se2=strel('disk',3); se3=strel('disk',5);

se4=strel('disk',7); se5=strel('disk',9); se6=strel('disk',11);

% Opening operators.

Img=double(Iori);

Open0=Img;

Open1=imopen(Img,se1); Open2=imopen(Img,se2); Open3=imopen(Img,se3);

Open4=imopen(Img,se4); Open5=imopen(Img,se5); Open6=imopen(Img,se6);

% % Reconstruct

% Open1=imreconstruct(Open1,Img); Open2=imreconstruct(Open2,Img); Open3=imreconstruct(Open3,Img);

% Open4=imreconstruct(Open4,Img); Open5=imreconstruct(Open5,Img); Open6=imreconstruct(Open6,Img);

% % Difference between two adjacent opening results.

DI0=Open0-Open1; DI1=Open1-Open2; DI2=Open2-Open3;

DI3=Open3-Open4; DI4=Open4-Open5; DI5=Open5-Open6;

% DI0=Open0.*Open1; DI1=Open1.*Open2; DI2=Open2.*Open3;

% DI3=Open3.*Open4; DI4=Open4.*Open5; DI5=Open5.*Open6;

% DI0=Open1-Open0; DI1=Open2-Open0; DI2=Open3-Open0;

% DI3=Open4-Open0; DI4=Open5-Open0; DI5=Open6-Open0;

% search max abs (DD) and order num (ind)

% % chosen 1 : only one scale according to maximum difference

% -----------------------------------------------------------

bInd=ones(nrow,ncol);

DD=DI0;

ind0 = bInd;

KK=DI1>=DD; ind1 = bInd.*double(KK) + double(KK); bInd = bInd + 1; DD=max(DD,DI1);

KK=DI2>=DD; ind2 = bInd.*double(KK) + double(KK); bInd = bInd + 1; DD=max(DD,DI2);

KK=DI3>=DD; ind3 = bInd.*double(KK) + double(KK); bInd = bInd + 1; DD=max(DD,DI3);

KK=DI4>=DD; ind4 = bInd.*double(KK) + double(KK); bInd = bInd + 1; DD=max(DD,DI4);

KK=DI5>=DD; ind5 = bInd.*double(KK) + double(KK); bInd = bInd + 1; DD=max(DD,DI5);

ind=max(ind0,ind1);ind=max(ind,ind2); ind=max(ind,ind3); ind=max(ind,ind4); ind=max(ind,ind5);

%%

%Multi-scale opening for background estimation and enhancement

bk=zeros(nrow,ncol); % background

kk=zeros(nrow,ncol);

for i=1:nrow

for j=1:ncol

if ind(i,j)==0 bk(i,j) = Open0(i,j); end

if ind(i,j)==1 bk(i,j) = Open1(i,j); end

if ind(i,j)==2 bk(i,j) = Open2(i,j); end

if ind(i,j)==3 bk(i,j) = Open3(i,j); end

if ind(i,j)==4 bk(i,j) = Open4(i,j); end

if ind(i,j)==5 bk(i,j) = Open5(i,j); end

if ind(i,j)==6 bk(i,j) = Open6(i,j); end

end

end

figure;

subplot(1,3,1);imshow(Iori); title('Original Image');

subplot(1,3,2);imshow(bk,[]);title('Multiscale background estimation');

kk=Img-bk;

I2=uint8(kk-1);

I3=imadjust(I2,stretchlim(I2),[0 1]);

subplot(1,3,3);imshow(I3);title('Multiscale top-hat result');

% imwrite(I3,'img002enResult005.bmp');