1. 框架

2. opticalFlow_label

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:hfor j = 1:wif first_label_bw(i, j) == 1
%             xPos = [xPos, i];       % 保存为行
%             yPos = [yPos, j];step = step + 1;if mod(step, 1) == 0xPos = [xPos; j];   % 保存为列yPos = [yPos; i];endendend
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_numnum = 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_hfor j = 1:mask_wif fc(i, j) == 1mask_frame(i, j, 1) = 255;mask_frame(:, :, 2) = currFrame;mask_frame(i, j, 3) = 0;end endendfigure(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', [3 3], 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 fineg = zeros(length(xPos), 6, pyramidLayer + 1);for i = pyramidLayer:-1:1im1 = im2double(fr1Pyramid{i});im2 = im2double(fr2Pyramid{i});% for each point, calculate I_x, I_yIx_m = filter2([-1 0 1; -1 0 1; -1 0 1],im1, 'same'); % partial on xIy_m = filter2([-1 -1 -1; 0 0 0; 1 1 1],im1,  'same'); % partial on yxPosPyramid = 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); 0 0 1];%粗精度下的运动矩阵matGmatG = [1 0 g(xx, 1, i+1); 0 1 g(xx, 2, i+1); 0 0 1];matX0 = [1 0 0;0 1 0;0 0 1] * 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为组合仿射矩阵，迭代计算，相对于fr1x = calcAAOF(Ix_m, Iy_m, im1, im2, xPosPyramid, yPosPyramid, ww, xx, x0);% affine motion combinationmatX1 = [1+x(3) x(4) x(1); x(5) 1+x(6) x(2); 0 0 1];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 = [2 0 0;0 2 0;0 0 1] * 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];endendendfor xx = 1:length(xPos)u(xx) = g(xx, 1, 1);v(xx) = g(xx, 2, 1);
endend% 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);
elseweightMatPhase = zeros(ww, ww);for i = 1 : wwfor j = 1 : wwxc = 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;endendendend
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 %xfor jSecond = -w:w %yiiMoveSecond = 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 + 1 1]);
yCoord = repmat(tempY', [1 2*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 IyA = [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 Itb = [b11;b12;b13;b14;b15;b16];%x = A \ b;
x = pinv(A) * b;
nanIdx = find(isnan(x) == 1);
x(nanIdx) = zeros(size(nanIdx));
endfunction 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