1. 题目描述

断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1m m的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。

假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。

现有某管道的相继100张平行切片图象，记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、…、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。为简化起见，假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为1。

取坐标系的Z轴垂直于切片，第1张切片为平面Z=0，第100张切片为平面Z=99。Z=z切片图象中象素的坐标依它们在文件中出现的前后次序为

（-256，-256，z），（-256，-255，z），…（-256，255，z），

（-255，-256，z），（-255，-255，z），…（-255，255，z），

……

（ 255，-256，z），（ 255，-255，z），…（255，255，z）。

试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。

第2页是100张平行切片图象中的6张，全部图象请从网上（http://mcm.edu.cn）下载。

关于BMP图象格式可参考：

1.《Visual C++数字图象处理》第12页2.3.1节。何斌等编著，人民邮电出版社，2001年4月。

2.http://www.dcs.ed.ac.uk/home/mxr/gfx/2d/BMP.txt

2.题目分析

对题目进行分析可知，题目的重点为求每一张图片的最大内切圆

求最大内切圆的matlab代码：

function [max_num, max_pos]=GetMax(fig)
bw_fig = fig;
black_dot = zeros(length(find(bw_fig==0)), 2);
board = [];
pos = 1;
for posi=1:length(bw_fig(:,1))
    for posj=1:length(bw_fig(:,1))
        if bw_fig(posi, posj)==0
            black_dot(pos, 1) = posi;
            black_dot(pos, 2) = posj;
            pos  = pos + 1;
        end
    end
end
edge_pos = edge(bw_fig);
pos = 1;
for posi=2:length(bw_fig(:,1))-1
    for posj=2:length(bw_fig(:,1))-1
        if edge_pos(posi, posj)==1
            board(pos, 1) = posi;
            board(pos, 2) = posj;
            pos = pos + 1;
        end
    end
end
max_num = 0;
max_pos = [1, 1];
for posi=1:length(black_dot(:,1))
    pre_min = inf;
    for posj=1:length(board(:,1))
        tmp = sqrt((board(posj,1)-black_dot(posi,1))^2+(board(posj,2)-black_dot(posi,2))^2);
        if tmp <=0.0001
            tmp = inf;
        end
        if tmp < pre_min
            pre_min = tmp;
        end
    end
    if pre_min > max_num
        max_num = pre_min;
        max_pos = black_dot(posi, :);
    end
end

返回值max_num表示内切圆半径，max_pos表示内切圆圆心坐标。

写脚本文件对图片进行处理，画出圆，需要修改文件路径。

clear;
close all;
clc;
dir_path = 'data/1.bmp';
fig = cell(1, 1); %RGB图像
result = cell(1, 1); %RGB图像
for pos=1:1
    i = imread(dir_path); %读取图像
    fig{pos} = i;
end

for pos=1:length(fig)
    result{pos} = fig{pos}*255;
end

dis = 10;
max_num = zeros(length(fig), 1);
max_pos = zeros(length(fig), 2);
for pos=1:1
    [max_num(pos), max_pos(pos, :)] = GetMax(fig{pos});
    result{pos}(max_pos(pos, 1), max_pos(pos, 2)) = 0.5;
    for posi=1:length(result{pos}(:,1))
        for posj=1:length(result{pos}(1,:))
            tmp = sqrt((posi-max_pos(pos, 1))^2+(posj-max_pos(pos, 2))^2);
            if abs(tmp-max_num(pos))<=dis
                result{pos}(posi, posj) = 150;
            end
        end
    end
end

for pos=1:1
    imwrite(result{pos}, "result/"+num2str(pos)+".bmp");
end

3.效果

内切圆

4.组合处理

在题目中，共有100张图片，通过逐一修改图片的路径获取所有内切圆是不现实的，需要统一处理。可以使用matlab读取读取文件，获取所有图片，进而统一处理。

100张图片统一处理

5.完整代码

clear;
close all;
clc;
dir_path = 'data';
name = dir(fullfile(dir_path, '*.bmp')); %获取所有bmp文件
fig = cell(1, length(name)); %RGB图像
result = cell(1, length(name)); %RGB图像
for pos=1:length(name)
    i = imread(dir_path+"/"+name(pos).name); %读取图像
    fig{pos} = i;
end

for pos=1:length(fig)
    result{pos} = fig{pos}*255;
end

dis = 5;
max_num = zeros(length(fig), 1);
max_pos = zeros(length(fig), 2);
for pos=1:100
    [max_num(pos), max_pos(pos, :)] = GetMax(fig{pos});
    result{pos}(max_pos(pos, 1), max_pos(pos, 2)) = 0.5;
    for posi=1:length(result{pos}(:,1))
        for posj=1:length(result{pos}(1,:))
            tmp = sqrt((posi-max_pos(pos, 1))^2+(posj-max_pos(pos, 2))^2);
            if abs(tmp-max_num(pos))<=dis
                result{pos}(posi, posj) = 150;
            end
        end
    end
end

for pos=1:100
    imwrite(result{pos}, "result/"+num2str(pos)+".bmp");
end