基于matlab的大鼠旷场试验二值化追踪

大鼠旷场试验视频的处理搞不定，没有找到能用的旷场试验软件，不会写代码，扒拉别人的代码删删改改拼的，Matlab版本R2021a

主体代码参考：MATLAB科研图像处理——动物运动轨迹追踪 - 知乎

包括位置信息、移动速度、移动距离、穿越次数的输出，实现方法比较简陋，也没有做批处理和整合，不过根据视频拍摄的不同应该也是需要自己调的。没有设置运动的阈值，所以动物即便在原地待着也会有位移，运动的结果会偏高不少，不过嘛，系统误差也算误差？或者有老哥自己整个阈值处理也行。

直立次数、修饰行为需要自己看视频慢慢统计，区域停留时间这个有位置信息了我就懒得再弄了。

当前是对白色大鼠黑色背景的追踪，需要追踪其他颜色比如黑色小鼠白色背景的话调二值化那里的小于大于号。

要是有老哥愿意把这个玩意整合一下就行了。

% 本脚本演示如何在开放场地测试期间跟踪个体小鼠的位置
% 原作者：Ethan Zhao，2022年8月
% 原教程：https://zhuanlan.zhihu.com/p/550378280
% 修改者：Myron Ziyan，2023年7月

% 清除之前的数据和图形窗口
clear; close all;

% 加载视频文件并获取基本信息
rawVideo = VideoReader("空白-3.mp4"); % 视频文件名，如有需要，请更新文件名
frameNum = rawVideo.NumFrames; % 视频总帧数
frameRate = rawVideo.FrameRate; % 视频帧率

% 定义感兴趣区域 (ROI)
Width = 600; % ROI宽度
Height = 600; % ROI高度
xPos = 45; % ROI在x方向的起始位置 MATLAB以1开始,ImageJ以0开始
yPos = 428; % ROI在y方向的起始位置
xRange = xPos:(xPos + Width - 1); % ROI在x方向上的范围
yRange = yPos:(yPos + Height - 1); % ROI在y方向上的范围

%% 跟踪小鼠位置
bwThd = 35; % 二值化阈值
sizeThd = 200; % 过滤小噪声的面积阈值
timeStep = 1; % 帧间隔

% 创建存储小鼠位置信息的数组
mousePos = zeros(floor(frameNum / timeStep), 2);

% 遍历每一帧图像并进行处理
for ii = 1:size(mousePos, 1)
    currentFrame = 1 + timeStep * (ii - 1); % 当前帧数
    img = read(rawVideo, currentFrame); % 读取当前帧图像

    % 对ROI进行二值化处理
    img = rgb2gray(img(yRange, xRange, :)); % 提取ROI中的图像
    imgBW = imfill(img > bwThd, 'holes'); % 二值化处理，大于号或小于号决定选择的二值化部分，当前为选择白色部分
    imgBW = bwareaopen(imgBW, sizeThd); % 过滤小噪声

    % 选择最大的对象并记录位置
    tmpStats = regionprops(imgBW, 'Area', 'Centroid'); % 获取区域属性
    [~, index] = max([tmpStats.Area]); % 找到最大的区域
    mousePos(ii, :) = tmpStats(index).Centroid; % 记录小鼠位置
    disp(['第', num2str(currentFrame), '帧已完成。']); % 显示进度
end

% 将存储的小鼠位置信息保存为.mat文件，文件名由TimeStep数值决定。
save(['mousePos-timeStep', num2str(timeStep), '.mat'], 'mousePos', 'xRange', 'yRange', 'timeStep');

% 显示x和y的最大值和最小值
maxX = max(mousePos(:, 1));
minX = min(mousePos(:, 1));
maxY = max(mousePos(:, 2));
minY = min(mousePos(:, 2));

disp('X的最大值: ' + string(maxX));
disp('X的最小值: ' + string(minX));
disp('Y的最大值: ' + string(maxY));
disp('Y的最小值: ' + string(minY));

% 将位置数据保存为CSV文件
csvFileName = sprintf('%s mouse_position.csv', erase(rawVideo.Name, '.mp4'));
csvwrite(csvFileName, mousePos);

第一块是二值化和追踪，二值化阈值需要自己调整，我用的imageJ来看的，按理说用ps或其他软件应该都可以，找个能限制灰度通道直方图的软件，自己慢慢调就行。

第二块是输出轨迹运动图的，没啥用，毕竟有了位置信息就能用其他方式作图了。

load('mousePos-timeStep1.mat'); % 加载保存的小鼠位置数据，如有需要，请更新文件名

figureID = 9; % 设置图形窗口的ID
figure(figureID);
set(gca, 'YDir','reverse'); % 将y轴方向设置为反向，使图像与原始图像方向一致
xlim([0 length(xRange)]); % 设置x轴的显示范围
ylim([0 length(yRange)]); % 设置y轴的显示范围
set(gca,'visible','off'); % 关闭坐标轴显示

for ii = 1:size(mousePos,1)-1
    tmpX = mousePos(ii:ii+1,1);
    tmpY = mousePos(ii:ii+1,2);

    figure(figureID); hold on; 
    plot(tmpX,tmpY,'r-'); % 在图形窗口中画出当前帧的小鼠位置轨迹（用红色线段表示）
end

% 保存生成的图形为"轨迹图.png"
outputFileName = '轨迹图.png';
saveas(figureID, outputFileName);

第三块是热图绘制，老实说我对这个热图绘制的结果很不满意，但不知道怎么搞别的热图绘制。

%% plot trace heatmap
% 加载之前保存的小鼠位置数据（存储了每帧小鼠的位置信息）
load('mousePos-timeStep1.mat');% 如有需要，请更新文件名

% 设置高斯模糊参数
blurSig = 8; 

% 创建热图（用于展示小鼠位置的密度分布）
% heatField 是一个二维矩阵，其大小为 length(yRange) × length(xRange)，初始值全部为0
heatField = zeros(length(yRange), length(xRange));

% 遍历每一帧小鼠的位置信息并构建热图
for ii = 1:size(mousePos, 1) - 1
    % 获取当前帧和下一帧的小鼠位置
    tmpFrame = mousePos(ii, :);
    nextFrame = mousePos(ii + 1, :);

    % 将当前帧和下一帧的x、y坐标分别组成两个一维数组x和y
    x = [tmpFrame(1), nextFrame(1)];
    y = [tmpFrame(2), nextFrame(2)];

    % 根据当前帧和下一帧的坐标，计算在两帧之间的路径上的一系列坐标点
    nPoints = ceil(sqrt((x(2) - x(1)).^2 + (y(2) - y(1)).^2)) + 1;
    xvalues = round(linspace(x(1), x(2), nPoints));
    yvalues = round(linspace(y(1), y(2), nPoints));

    % 将这一系列坐标点转换成线性索引，用于更新热图中的对应位置的数值
    tmpIdx = sub2ind(size(heatField), yvalues, xvalues);

    % 将沿路径的坐标点对应的位置在热图中的数值加1，用于路径集成
    heatField(tmpIdx) = heatField(tmpIdx) + 1; 
end

% 对构建好的热图进行高斯模糊处理，以平滑显示效果
heatFieldBlur = imgaussfilt(heatField, blurSig);

% 创建图形窗口并绘制热图（以颜色表示小鼠位置的密度分布）
figure(10);
imagesc(heatFieldBlur); % 在图形窗口中绘制热图
axis image; % 设置坐标轴纵横比例相等，保持热图像素显示不变形
colormap jet; % 设置颜色映射为jet colormap
colorbar; % 添加颜色条作为图例
outputFileName = '轨迹热图.png';

第四块是移动速度和移动距离的计算

% 加载保存的小鼠位置数据
load('mousePos-timeStep1.mat'); % 如有需要，请更新文件名

% 视频帧率和像素与实际距离的转换参数
frameRate = 30; % 视频帧率（假设每秒30帧）
pixelToDistance = 80 / 600; % 像素与实际距离的转换参数（80cm对应ROI宽度600像素）

% 计算每帧的时间间隔（假设每帧为1秒）
timeInterval = 1 / frameRate; % 每帧的时间间隔（秒）

% 计算每帧移动的像素值
deltaX = diff(mousePos(:, 1)); % X方向上的位移（像素）
deltaY = diff(mousePos(:, 2)); % Y方向上的位移（像素）

% 计算每帧移动的距离（厘米）
distancePerFrame = sqrt(deltaX.^2 + deltaY.^2) * pixelToDistance;

% 计算移动距离和平均速度
totalDistance = sum(distancePerFrame); % 移动总距离（厘米）
averageSpeed = totalDistance / (size(mousePos, 1) - 1) / timeInterval; % 平均速度（厘米/秒）

% 显示计算结果
disp(['移动总距离：', num2str(totalDistance), '厘米']);
disp(['平均速度：', num2str(averageSpeed), '厘米/秒']);

第五块是区域穿越计数，现在是1/1/1分了九块，想调整自己调

% 从保存的.mat文件中加载小鼠位置数据
load('mousePos-timeStep1.mat'); % 如有需要，请更新文件名

% 获取帧数
numFrames = size(mousePos, 1);

% 计算四条线的Y轴位置
yPos1 = 200.00; % 第一条线的Y轴位置，需要根据浮点精度增加.000
yPos2 = 400.00; % 第二条线的Y轴位置

% 计算四条线的X轴位置
xPos1 = 200.000; % 第一条线的X轴位置
xPos2 = 400.000; % 第二条线的X轴位置

% 初始化八个方向的穿越计数器
crossings_y1_from_above = 0;
crossings_y1_from_below = 0;
crossings_y2_from_above = 0;
crossings_y2_from_below = 0;
crossings_x1_from_left = 0;
crossings_x1_from_right = 0;
crossings_x2_from_left = 0;
crossings_x2_from_right = 0;

% 遍历每一帧来统计穿越次数
for frame = 2:numFrames
    % 获取当前帧小鼠的位置
    current_pos = mousePos(frame, :);
    % 获取上一帧小鼠的位置
    prev_pos = mousePos(frame - 1, :);
    
    % 检测是否从y 1/3上方穿越
if prev_pos(2) <= yPos1 && current_pos(2) > yPos1
    crossings_y1_from_above = crossings_y1_from_above + 1; % 从上方穿越次数加一
end

% 检测是否从y 1/3下方穿越
if prev_pos(2) >= yPos1 && current_pos(2) < yPos1
    crossings_y1_from_below = crossings_y1_from_below + 1; % 从下方穿越次数加一
end

% 检测是否从y 2/3上方穿越
if prev_pos(2) <= yPos2 && current_pos(2) > yPos2
    crossings_y2_from_above = crossings_y2_from_above + 1; % 从上方穿越次数加一
end

% 检测是否从y 2/3下方穿越
if prev_pos(2) >= yPos2 && current_pos(2) < yPos2
    crossings_y2_from_below = crossings_y2_from_below + 1; % 从下方穿越次数加一
end

% 检测是否从x 1/3左方穿越
if prev_pos(1) <= xPos1 && current_pos(1) > xPos1
    crossings_x1_from_left = crossings_x1_from_left + 1; % 从左方穿越次数加一
end

% 检测是否从x 1/3右方穿越
if prev_pos(1) >= xPos1 && current_pos(1) < xPos1
    crossings_x1_from_right = crossings_x1_from_right + 1; % 从右方穿越次数加一
end

% 检测是否从x 2/3左方穿越
if prev_pos(1) <= xPos2 && current_pos(1) > xPos2
    crossings_x2_from_left = crossings_x2_from_left + 1; % 从左方穿越次数加一
end

% 检测是否从x 2/3右方穿越
if prev_pos(1) >= xPos2 && current_pos(1) < xPos2
    crossings_x2_from_right = crossings_x2_from_right + 1; % 从右方穿越次数加一
end
end

% 计算总的八个方向的穿越次数
total_crossings = crossings_y1_from_above + crossings_y1_from_below + ...
                  crossings_y2_from_above + crossings_y2_from_below + ...
                  crossings_x1_from_left + crossings_x1_from_right + ...
                  crossings_x2_from_left + crossings_x2_from_right;

% 显示总的八个方向的穿越次数
disp('八个方向的总穿越次数:');
disp('Y1 From Above: ' + string(crossings_y1_from_above));
disp('Y1 From Below: ' + string(crossings_y1_from_below));
disp('Y2 From Above: ' + string(crossings_y2_from_above));
disp('Y2 From Below: ' + string(crossings_y2_from_below));
disp('X1 From Left: ' + string(crossings_x1_from_left));
disp('X1 From Right: ' + string(crossings_x1_from_right));
disp('X2 From Left: ' + string(crossings_x2_from_left));
disp('X2 From Right: ' + string(crossings_x2_from_right));

% 显示总的穿越次数
disp('所有方向的总穿越次数:');
disp(total_crossings);

% 将总穿越次数保存到CSV文件
csvwrite('total_crossings.csv', total_crossings);