close all; clear; clc;% 使用光流进行标签的生成
%% 视频帧的读取
npy_data = readNPY('train.npy');%% 提取标签的坐标
first_label = squeeze(npy_data(2,1,:,:));
h = fspecial("gaussian",[3,3],1);
first_label_g = imfilter(first_label, h,'replicate');%'conv'
first_label_c = edge(first_label_g,"canny");% imshow(uint8(first_label_c*255));
first_label_double = im2double(first_label_c);
first_label_bw = im2bw(first_label_double,0.5);% imshow(uint8(first_label_bw *255));[h, w]= size(first_label);
xPos =[];
yPos =[];
step =0;for i =1:h
for j =1:w
if first_label_bw(i, j)==1% xPos =[xPos, i];% 保存为行
% yPos =[yPos, j];
step = step +1;if mod(step,1)==0
xPos =[xPos; j];% 保存为列
yPos =[yPos; i];
end
end
end
end
%% 逐帧处理
first_frame = squeeze(npy_data(1,1,:,:));
first_frame = uint8(first_frame);% imshow(first_frame);[c,frame_num,img_h,img_w]= size(npy_data);
num =0;
save_npy(1,1,:,:)= first_frame;
save_npy(2,1,:,:)= first_frame;% 预留一个通道，用于保存标签
for i =2%2:frame_num
num = num +1;
currFrame = squeeze(npy_data(1,i,:,:));
currFrame = uint8(currFrame);
pyramidLayer =4;
kernelSize =3;
sigma =1.5;
iterNumMax =5;
ww =13;[u, v]= affineLKopticalFlow(first_frame, currFrame, xPos, yPos, pyramidLayer, kernelSize, sigma, iterNumMax, ww);% 显示
newFrame = repmat(currFrame,[1,1,3]);
new_xPos = xPos + u;
new_yPos = yPos + v;
newFrame1 = zeros(size(currFrame));for kk =1:length(xPos)%显示
newFrame(int16(new_yPos(kk)), int16(new_xPos(kk)),1)=255;
newFrame(:,:,2)= currFrame;
newFrame(int16(new_yPos(kk)), int16(new_xPos(kk)),3)=0;
newFrame1(int16(new_yPos(kk)), int16(new_xPos(kk)))=1;
end
save_npy(1,i,:,:)= currFrame;
save_npy(2,i,:,:)= first_frame;% writeNPY(save_npy,"test.npy");
figure(1)
imshow(newFrame)
title("Pre")% 形态学处理 离散点的连接
se = strel('disk',51);% 加大半径 可以进行填充
fc = imclose(newFrame1, se);% imwrite(fc,"fc.jpg")% bw = im2bw(fc);% fill_img = imfill(bw,"holes");% 封闭区域
figure(2)
imshow(fc)
title("fc")%将mask显示在图片上
mask_frame = repmat(currFrame,[1,1,3]);[mask_h, mask_w, channel]= size(mask_frame);for i =1:mask_h
for j =1:mask_w
if fc(i, j)==1
mask_frame(i, j,1)=255;
mask_frame(:,:,2)= currFrame;
mask_frame(i, j,3)=0;
end
end
end
figure(4)
imshow(mask_frame)
title("mask")
first_frame = currFrame;
xPos = double(int16(new_xPos));
yPos = double(int16(new_yPos));
end
3. 光流
% author : huangcan
% date :2019.4.22% comment : this code is an implementation of Local antomical affine optical flow algorithm
% reference : paper:"Fast Left Ventricle Tracking in 3D Echocardiographic Data Using Anatomical Affine Optical Flow"% parameter :% lastFrame 上一帧的图像
% currFrame 当前帧的图像
% xPos 要跟踪的点的横坐标
% yPos 要跟踪的点的纵坐标
% pyramidLayer 尺度金字塔的层数（包含原始层）
% kernelSize 高斯模糊的模板尺寸大小，3即可
% sigma 高斯模糊标准差
% iterNumMax 迭代次数上限，一般４次即可
% ww ROI框的尺寸，注意确保多尺度后的尺寸仍在图像范围内
function [u, v]= affineLKopticalFlow(lastFrame, currFrame, xPos, yPos, pyramidLayer, kernelSize, sigma, iterNumMax, ww)
u = zeros(length(xPos),1);
v = zeros(length(xPos),1);%% build image pyramid
fr1Pyramid = cell(pyramidLayer,1);
fr2Pyramid = cell(pyramidLayer,1);
fr1Pyramid{1}= lastFrame;
fr2Pyramid{1}= currFrame;
h = fspecial('gaussian',[kernelSize kernelSize], sigma);for i =2:pyramidLayer
%h = fspecial('gaussian',[33],0.83^(i-2));
temp1 = imfilter(fr1Pyramid{i-1}, h,'same');
temp2 = imfilter(fr2Pyramid{i-1}, h,'same');
fr1Pyramid{i}= imresize(temp1,1/2);
fr2Pyramid{i}= imresize(temp2,1/2);
end
%% coarse to fine
g = zeros(length(xPos),6, pyramidLayer +1);for i = pyramidLayer:-1:1
im1 = im2double(fr1Pyramid{i});
im2 = im2double(fr2Pyramid{i});%for each point, calculate I_x, I_y
Ix_m = filter2([-101;-101;-101],im1,'same');% partial on x
Iy_m = filter2([-1-1-1;000;111],im1,'same');% partial on y
xPosPyramid = xPos ./2^(i-1);
yPosPyramid = yPos ./2^(i-1);for xx =1:length(xPos)%x0 =[u; v; ux; uy; vx; vy];这个不是运动矩阵，而是参数矩阵
%matG =[1+g(xx,3, i+1) g(xx,4, i+1) g(xx,1, i+1); g(xx,5, i+1)1+g(xx,6, i+1) g(xx,2, i+1);001];%粗精度下的运动矩阵matG
matG =[10 g(xx,1, i+1);01 g(xx,2, i+1);001];
matX0 =[100;010;001]* matG;
x0 =[matX0(1,3);matX0(2,3);matX0(1,1)-1;matX0(1,2);matX0(2,1);matX0(2,2)-1];for k =1:iterNumMax
% x0为组合仿射矩阵，迭代计算，相对于fr1
x = calcAAOF(Ix_m, Iy_m, im1, im2, xPosPyramid, yPosPyramid, ww, xx, x0);% affine motion combination
matX1 =[1+x(3) x(4) x(1); x(5)1+x(6) x(2);001];
matX0 = matX1 * matX0;
x0 =[matX0(1,3); matX0(2,3); matX0(1,1)-1; matX0(1,2); matX0(2,1); matX0(2,2)-1];
end
%此时X0为该层金字塔相对于第一帧的仿射参数，matX0为相对于第一帧的仿射矩阵
if(i >=2)% 为更细精度的分析进行运动放大
iterFinal =[200;020;001]* matX0;
g(xx,:,i)=[iterFinal(1,3); iterFinal(2,3); iterFinal(1,1)-1; iterFinal(1,2); iterFinal(2,1); iterFinal(2,2)-1];else% 最终层不需要变化
iterFinal = matX0;
g(xx,:,i)=[iterFinal(1,3); iterFinal(2,3); iterFinal(1,1)-1; iterFinal(1,2); iterFinal(2,1); iterFinal(2,2)-1];
end
end
end
for xx =1:length(xPos)
u(xx)= g(xx,1,1);
v(xx)= g(xx,2,1);
end
end
%sum
function x = calcAAOF(Ix_m, Iy_m, im1, im2, xPos, yPos, ww, ix, x0)%x0 =[u v ux1 ux2 vx1 vx2];
w = floor(ww/2);if(length(xPos)>100)
weightMatPhase = ones(ww,ww);else
weightMatPhase = zeros(ww, ww);for i =1: ww
for j =1: ww
xc = int16(xPos(ix))-(ww +1)/2+ i;
yc = int16(yPos(ix))-(ww +1)/2+ j;
weightMatPhase(j,i)=0;for ii =1:length(xPos)
dis = sqrt((double(xc)- double(xPos(ii)))^2+(double(yc)- double(yPos(ii)))^2);if(dis < sqrt(2)*ww)
weightMatPhase(j, i)= weightMatPhase(j, i)+1;
end
end
end
end
end
weight = weightMatPhase(:);
xx = ix;
i = int16(xPos(xx));
j = int16(yPos(xx));%Ix = Ix_m(j-w:j+w, i-w:i+w);
Ix = iGetN(Ix_m,i,j,2*w+1);%Iy = Iy_m(j-w:j+w, i-w:i+w);
Iy = iGetN(Iy_m, i, j,2*w+1);% affine model
newX = zeros(size(Ix));
newY = zeros(size(Iy));for iSecond =-w:w %x
for jSecond =-w:w %y
iiMoveSecond = x0(1)+ x0(3)*double(iSecond)+ x0(4)*double(iSecond);
jjMoveSecond = x0(2)+ x0(5)*double(jSecond)+ x0(6)*double(jSecond);
newX(jSecond+w+1,iSecond+w+1)= double(i)+ iSecond + iiMoveSecond;
newY(jSecond+w+1,iSecond+w+1)= double(j)+ jSecond + jjMoveSecond;
end
end
% bilinear interp2
newXX = linspace(1, size(im2,2), size(im2,2));
newYY = linspace(1, size(im2,1), size(im2,1));
Xfirst = repmat(newXX,[size(im2,1)1]);
Yfirst = repmat(newYY',[1 size(im2,2)]);
newIm2 = interp2(Xfirst, Yfirst, im2, newX,newY);% calculate time difference
%It = newIm2 - im1(j-w:j+w, i-w:i+w);
It = newIm2 - iGetN(im1, i, j,2*w+1);
Ix = Ix(:);
Iy = Iy(:);
It = It(:);
tempX = linspace(-w, w,2*w+1);
tempY = linspace(-w, w,2*w +1);
xCoord = repmat(tempX,[2*w +11]);
yCoord = repmat(tempY',[12*w+1]);
xCoord = xCoord(:);
yCoord = yCoord(:);
A11 =sum(sum(weight .* Ix .* Ix));% w Ix Ix
A12 =sum(sum(weight .* Ix .* Iy));% w Ix Iy
A13 =sum(sum(xCoord .* weight .* Ix .* Ix));% x w Ix Ix
A14 =sum(sum(yCoord .* weight .* Ix .* Ix));% y w Ix Ix
A15 =sum(sum(xCoord .* weight .* Ix .* Iy));% x w Ix Iy
A16 =sum(sum(yCoord .* weight .* Ix .* Iy));% y w Ix Iy
A21 =sum(sum(weight .* Ix .* Iy));% w Ix Iy
A22 =sum(sum(weight .* Iy .* Iy));% w Iy Iy
A23 =sum(sum(xCoord .* weight .* Ix .* Iy));%x w Ix Iy
A24 =sum(sum(yCoord .* weight .* Ix .* Iy));%y w Ix Iy
A25 =sum(sum(xCoord .* weight .* Iy .* Iy));%x w Iy Iy
A26 =sum(sum(yCoord .* weight .* Iy .* Iy));%y w Iy Iy
A31 =sum(sum(xCoord .* weight .* Ix .* Ix));%x w Ix Ix
A32 =sum(sum(xCoord .* weight .* Ix .* Iy));% x w Ix Iy
A33 =sum(sum(xCoord.^2.* weight .* Ix .* Ix));%x^2 w Ix Ix
A34 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Ix));% x y Ix Ix
A35 =sum(sum(xCoord .^2.* weight .* Ix .* Iy));% x^2 w Ix Iy
A36 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Iy));% x y w Ix Iy
A41 =sum(sum(yCoord .* weight .* Ix .* Ix));% y w Ix Ix
A42 =sum(sum(yCoord .* weight .* Ix .* Iy));% y w Ix Iy
A43 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Ix));% x y w Ix Ix
A44 =sum(sum(yCoord .^2.* weight .* Ix .* Ix));%y^2 w Ix Ix
A45 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Iy));% x y w Ix Iy
A46 =sum(sum(yCoord .^2.* weight .* Ix .* Iy));% y^2 w Ix Iy
A51 =sum(sum(xCoord .* weight .* Ix .* Iy));% x^2 w Ix Iy
A52 =sum(sum(xCoord .* weight .* Iy .* Iy));% x w Iy Iy
A53 =sum(sum(xCoord .^2.* weight .* Ix .* Iy));% x^2 w Ix Iy
A54 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Iy));% x y w Ix Iy
A55 =sum(sum(xCoord .^2.* weight .* Iy .* Iy));% x^2 w Iy Iy
A56 =sum(sum(xCoord .* yCoord .* weight .* Iy .* Iy));% x y w Iy Iy
A61 =sum(sum(yCoord .* weight .* Ix .* Iy));% y w Ix Iy
A62 =sum(sum(yCoord .* weight .* Iy .* Iy));%y w Iy Iy
A63 =sum(sum(xCoord .* yCoord .* weight .* Ix .* Iy));% x y w Ix Iy
A64 =sum(sum(yCoord .^2.* weight .* Ix .* Iy));% y^2 w Ix Iy
A65 =sum(sum(xCoord .* yCoord .* weight .* Iy .* Iy));% x y w Iy Iy
A66 =sum(sum(yCoord .^2.* weight .* Iy .* Iy));% y^2 w Iy Iy
A =[A11 A12 A13 A14 A15 A16; A21 A22 A23 A24 A25 A26; A31 A32 A33 A34 A35 A36; A41 A42 A43 A44 A45 A46; A51 A52 A53 A54 A55 A56; A61 A62 A63 A64 A65 A66];
b11 =-sum(sum(weight .* Ix .* It));%w Ix It
b12 =-sum(sum(weight .* Iy .* It));%w Iy It
b13 =-sum(sum(xCoord .* weight .* Ix .* It));%x w Ix It
b14 =-sum(sum(yCoord .* weight .* Ix .* It));% y w Ix It
b15 =-sum(sum(xCoord .* weight .* Iy .* It));% x w Iy It
b16 =-sum(sum(yCoord .* weight .* Iy .* It));% y w Iy It
b =[b11;b12;b13;b14;b15;b16];%x = A \ b;
x = pinv(A)* b;
nanIdx = find(isnan(x)==1);
x(nanIdx)= zeros(size(nanIdx));
end
function r=iGetN(m,x,y,N)[h,w,c]=size(m);
halfN = floor(N/2);
n1=halfN; n2=N-halfN-1;
r=zeros(N,N,c);
xmin=max(1,x-n1);
xmax=min(w,x+n2);
ymin=max(1,y-n1);
ymax=min(h,y+n2);
pxmin=halfN-(x-xmin)+1; pxmax=halfN+(xmax-x)+1;
pymin=halfN-(y-ymin)+1; pymax=halfN+(ymax-y)+1;
r(pymin:pymax,pxmin:pxmax,:)=m(ymin:ymax,xmin:xmax,:);
end
1819
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
