3D骨架提取算法的Python实现和C++实现,包括骨架提取以及写入树文件中:
(1)Python
import itk
# python -m pip install --upgrade pip
# python -m pip install itk-thickness3d
def OutputTXT(filePath, list_centerLinePoint):
with open(filePath, "w+") as fout:
fout.writelines(
[" ".join([str(i) for i in centerLinePoint]) + "\n" for centerLinePoint in list_centerLinePoint])
def Collect26Neighbors(image, p_index):
offset = [
[-1, -1, 1], [0, -1, 1], [1, -1, 1], [-1, -1, 0], [0, -1, 0], [1, -1, 0], [-1, -1, -1], [0, -1, -1],
[1, -1, -1],
[-1, 0, 1], [0, 0, 1], [1, 0, 1], [-1, 0, 0], [1, 0, 0], [-1, 0, -1], [0, 0, -1], [1, 0, -1],
[-1, 1, 1], [0, 1, 1], [1, 1, 1], [-1, 1, 0], [0, 1, 0], [1, 1, 0], [-1, 1, -1], [0, 1, -1], [1, 1, -1]]
Np = [1 if image.GetPixel(tuple([p_index[i] + offset[l][i] for i in range(3)])) == 1 else 0 for l in range(26)]
return Np
def ReturnId(p_index, list_centerLinePoint):
# [(id, x, y, z, parentId, leftChildId, rightChildId),...]
id = None
for i in range(len(list_centerLinePoint)):
if list_centerLinePoint[i][1:4] == list(p_index):
id = list_centerLinePoint[i][0]
break
return id
def ReturnIndex(j, p_index):
offset = [
[-1, -1, 1], [0, -1, 1], [1, -1, 1], [-1, -1, 0], [0, -1, 0], [1, -1, 0], [-1, -1, -1], [0, -1, -1],
[1, -1, -1],
[-1, 0, 1], [0, 0, 1], [1, 0, 1], [-1, 0, 0], [1, 0, 0], [-1, 0, -1], [0, 0, -1], [1, 0, -1],
[-1, 1, 1], [0, 1, 1], [1, 1, 1], [-1, 1, 0], [0, 1, 0], [1, 1, 0], [-1, 1, -1], [0, 1, -1], [1, 1, -1]]
re_index = [p_index[i] + offset[j][i] for i in range(3)]
return re_index
def ReturnJx(Np, count=1):
ct = 0
j = []
for i in range(26):
if Np[i] == 1:
j.append(i)
ct += 1
if count == ct:
break
return j[0] if count == 1 else j
def MarkAllPoints(image):
size = image.GetLargestPossibleRegion().GetSize()
id_ini = 0
list_centerline_vec = []
for i in range(size[0] - 1):
for j in range(size[1] - 1):
for k in range(size[2] - 1):
index = (i, j, k)
value = image.GetPixel(index)
if value == 0:
continue
else:
list_centerline_vec.append([id_ini, i, j, k]) # [(id, x, y, z),...]
id_ini = id_ini + 1
return list_centerline_vec
def GenerateTree(image):
centerline_vec = MarkAllPoints(image)
for i in range(len(centerline_vec)):
p_index = centerline_vec[i][1:4]
Np = Collect26Neighbors(image, p_index)
sum = 0
for j in range(26):
sum += Np[j]
if sum == 0:
centerline_vec[i].extend([-1, -1, -1])
# 判断是端点的情况
if sum == 1:
j = ReturnJx(Np, count=1)
if j <= 12:
parent_index = ReturnIndex(j, p_index)
centerline_vec[i].extend([ReturnId(parent_index, centerline_vec), -1, -1])
# 如果只有一个childId一律设置成左childId
if j >= 13:
leftChild_index = ReturnIndex(j, p_index)
centerline_vec[i].extend([-1, ReturnId(leftChild_index, centerline_vec), -1])
if sum == 2:
j = ReturnJx(Np, count=2)
parent_index = ReturnIndex(j[0], p_index)
leftChild_index = ReturnIndex(j[1], p_index)
centerline_vec[i].extend(
[ReturnId(parent_index, centerline_vec), ReturnId(leftChild_index, centerline_vec), -1])
if sum == 3:
j = ReturnJx(Np, count=3)
parent_index = ReturnIndex(j[0], p_index)
leftChild_index = ReturnIndex(j[1], p_index)
rightChild_index = ReturnIndex(j[2], p_index)
centerline_vec[i].extend([ReturnId(parent_index, centerline_vec), ReturnId(leftChild_index, centerline_vec),
ReturnId(rightChild_index, centerline_vec)])
return centerline_vec
if __name__ == '__main__':
input_filename = "3.nii.gz"
# input_filename = "3_mask.nii.gz"
img = itk.imread(input_filename)
thining_map = itk.BinaryThinningImageFilter3D.New(img)
# print(dir(thining_map))
itk.imwrite(thining_map, "3_thining.nii.gz")
image = thining_map.GetOutput()
centerline_vec_final = GenerateTree(image)
path = "CenterlinePoints.txt"
OutputTXT(path, centerline_vec_final)
(2)C++
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkBinaryThinningImageFilter3D.h"
#include "itkTestingMacros.h"
// PNG对应
#include <itkPNGImageIOFactory.h>
// BMP对应
#include <itkBMPImageIOFactory.h>
// JPG对应
#include <itkJPEGImageIOFactory.h>
// mhd格式图像
#include <itkMetaImageIOFactory.h>
// nii格式图片
#include <itkNiftiImageIOFactory.h>
#include <iostream>
#include <cmath>
#include <ctime>
#include <fstream>
#include <string>
const unsigned int Dimension = 3;
using PixelType = unsigned char;
using ImageType = itk::Image<PixelType, Dimension>;
typedef itk::Image<PixelType, Dimension> ImageType;
typedef itk::ImageFileReader<ImageType> ReaderType;
typedef itk::ImageFileWriter<ImageType> WriterType;
typedef itk::Point<PixelType, 3> PointType;
//下面是生成中心线树的代码
class CenterlinePoint {
public:
int x;
int y;
int z;
int id;
int parentId;
int liftChildId;
int rightChildId;
};
//输出坐标TXT文件
void OutputTXT(const char* path, std::vector<CenterlinePoint> centerline_vec_final) {
std::ofstream fout(path);
if (fout) {
for (int i = 0; i < centerline_vec_final.size(); i++) {
fout << centerline_vec_final[i].id << " " << centerline_vec_final[i].x << " " << centerline_vec_final[i].y << " " << centerline_vec_final[i].z << " " << centerline_vec_final[i].parentId << " " << centerline_vec_final[i].liftChildId << " " << centerline_vec_final[i].rightChildId << std::endl;
}
}
fout.close();
}
//对p的26邻接体素编号,并且改成0,1赋值(0:背景,1:前景)
void Collect26Neighbors(ImageType::Pointer image, int p_index[3], int *Np) {
int offset[26][3] = {
{-1,-1,1},{0,-1,1},{1,-1,1},{-1,-1,0},{0,-1,0},{1,-1,0},{-1,-1,-1},{0,-1,-1},{1,-1,-1},
{-1,0,1},{0,0,1},{1,0,1},{-1,0,0},{1,0,0},{-1,0,-1},{0,0,-1},{1,0,-1},
{-1,1,1},{0,1,1},{1,1,1},{-1,1,0},{0,1,0},{1,1,0},{-1,1,-1},{0,1,-1},{1,1,-1} };
ImageType::IndexType index;
ImageType::PixelType value;
for (int l = 0; l < 26; l++) {
index[0] = p_index[0] + offset[l][0];
index[1] = p_index[1] + offset[l][1];
index[2] = p_index[2] + offset[l][2];
value = image->GetPixel(index);
std::cout << int(value) << std::endl;
if (value == 1) { Np[l] = 1; }
else { Np[l] = 0; }
}
}
int ReturnId(int p_index[3], std::vector<CenterlinePoint> centerline_vec) {
int id = -1;
for (int i = 0; i < centerline_vec.size(); i++) {
if (centerline_vec[i].x == p_index[0] && centerline_vec[i].y == p_index[1] && centerline_vec[i].z == p_index[2]) {
id = centerline_vec[i].id;
}
}
return id;
}
void ReturnIndex(int re_index[3], int j, int p_index[3]) {
int offset[26][3] = {
{-1,-1,1},{0,-1,1},{1,-1,1},{-1,-1,0},{0,-1,0},{1,-1,0},{-1,-1,-1},{0,-1,-1},{1,-1,-1},
{-1,0,1},{0,0,1},{1,0,1},{-1,0,0},{1,0,0},{-1,0,-1},{0,0,-1},{1,0,-1},
{-1,1,1},{0,1,1},{1,1,1},{-1,1,0},{0,1,0},{1,1,0},{-1,1,-1},{0,1,-1},{1,1,-1} };
re_index[0] = p_index[0] + offset[j][0];
re_index[1] = p_index[1] + offset[j][1];
re_index[2] = p_index[2] + offset[j][2];
}
int ReturnJ1(int Np[26]) {
int j;
for (int i = 0; i < 26; i++) {
if (Np[i] == 1) {
j = i;
}
}
return j;
}
void ReturnJ2(int *j, int Np[26]) {
int count = 0;
for (int i = 0; i < 26; i++) {
if (Np[i] == 1) {
j[count++] = i;
}
if (count == 2)
break;
}
}
void ReturnJ3(int *j, int Np[26])
{
int count = 0;
for (int i = 0; i < 26; i++) {
if (Np[i] == 1) {
j[count++] = i;
}
if (count == 3)
break;
}
}
std::vector<CenterlinePoint> MarkAllPoints(ImageType::Pointer image)
{
std::vector<CenterlinePoint> centerline_vec;
ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
ImageType::IndexType index;
ImageType::PixelType value;
int id_ini = 0;
for (int i = 1; i < size[0] - 1; i++) {
for (int j = 1; j < size[1] - 1; j++) {
for (int k = 1; k < size[2] - 1; k++) {
index[0] = i;
index[1] = j;
index[2] = k;
value = image->GetPixel(index);
if (value == 0) { continue; }
else {
CenterlinePoint cpt;
cpt.x = i;
cpt.y = j;
cpt.z = k;
cpt.id = id_ini;
id_ini = id_ini + 1;
centerline_vec.push_back(cpt);
}
}
}
}
return centerline_vec;
}
std::vector<CenterlinePoint> GenerateTree(ImageType::Pointer image) {
std::vector<CenterlinePoint> centerline_vec;
centerline_vec = MarkAllPoints(image);
for (int i = 0; i < centerline_vec.size(); i++) {
int p_index[3];
p_index[0] = centerline_vec[i].x;
p_index[1] = centerline_vec[i].y;
p_index[2] = centerline_vec[i].z;
int Np[26];
Collect26Neighbors(image, p_index, Np);
int sum = 0;
for (int j = 0; j < 26; j++) {
sum = sum + Np[j];
}
if (sum == 0) {
centerline_vec[i].parentId = -1;
centerline_vec[i].liftChildId = -1;
centerline_vec[i].rightChildId = -1;
}
//判断是端点的情况
if (sum == 1) {
int j = ReturnJ1(Np);
if (j <= 12) {
int parent_index[3];
ReturnIndex(parent_index, j, p_index);
centerline_vec[i].parentId = ReturnId(parent_index, centerline_vec);
centerline_vec[i].liftChildId = -1;
centerline_vec[i].rightChildId = -1;
}
//如果只有一个childId一律设置成左childId
if (j >= 13) {
int liftChild_index[3];
ReturnIndex(liftChild_index, j, p_index);
centerline_vec[i].parentId = -1;
centerline_vec[i].liftChildId = ReturnId(liftChild_index, centerline_vec);
centerline_vec[i].rightChildId = -1;
}
}
if (sum == 2) {
int j[2];
ReturnJ2(j, Np);
int parent_index[3];
ReturnIndex(parent_index, j[0], p_index);
int liftChild_index[3];
ReturnIndex(liftChild_index, j[1], p_index);
centerline_vec[i].parentId = ReturnId(parent_index, centerline_vec);
centerline_vec[i].liftChildId = ReturnId(liftChild_index, centerline_vec);
centerline_vec[i].rightChildId = -1;
}
if (sum == 3) {
int j[3];
ReturnJ3(j, Np);
int parent_index[3];
ReturnIndex(parent_index, j[0], p_index);
int liftChild_index[3];
ReturnIndex(liftChild_index, j[1], p_index);
int rightChild_index[3];
ReturnIndex(rightChild_index, j[2], p_index);
centerline_vec[i].parentId = ReturnId(parent_index, centerline_vec);
centerline_vec[i].liftChildId = ReturnId(liftChild_index, centerline_vec);
centerline_vec[i].rightChildId = ReturnId(rightChild_index, centerline_vec);
}
}
return centerline_vec;
}
int main(int argc, char * argv[])
{
// Verify the number of parameters in the command line
if (argc == 2)
{
std::cerr << "Missing parameters." << std::endl;
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " inputImageFile outputImageFile" << std::endl;
return EXIT_FAILURE;
}
const char * infilename = "3.nii";
const char * outfilename = "3_thining.nii";
// Read image
//itk::ObjectFactoryBase::UnRegisterAllFactories();
//itk::NiftiImageIOFactory::RegisterOneFactory();
itk::ObjectFactoryBase::RegisterFactory(itk::NiftiImageIOFactory::New());
// itk::ObjectFactoryBase::RegisterFactory(itk::MetaImageIOFactory::New()); mhd
// itk::PNGImageIOFactory::New()或itk::BMPImageIOFactory::New() // 或itk::JPGImageIOFactory::New()
using ReaderType = itk::ImageFileReader<ImageType>;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infilename);
ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update());
// Define the thinning filter
using ThinningFilterType = itk::BinaryThinningImageFilter3D<ImageType, ImageType>;
ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
ITK_EXERCISE_BASIC_OBJECT_METHODS(thinningFilter, BinaryThinningImageFilter3D, ImageToImageFilter);
thinningFilter->SetInput(reader->GetOutput());
ITK_TRY_EXPECT_NO_EXCEPTION(thinningFilter->Update());
ImageType::Pointer image = thinningFilter->GetOutput();
// 以二叉树的形式保存结果,本人非CS出身,如有问题欢迎指正。
std::vector<CenterlinePoint> centerline_vec_final;
centerline_vec_final = GenerateTree(image);
const char* path = "CenterlinePoints.txt";
OutputTXT(path, centerline_vec_final);
// output to file
// itk::ObjectFactoryBase::RegisterFactory(itk::MetaImageIOFactory::New()); mhd
// itk::PNGImageIOFactory::New()或itk::BMPImageIOFactory::New() // 或itk::JPGImageIOFactory::New()
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetInput(thinningFilter->GetOutput());
writer->SetFileName(outfilename);
ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update());
std::cout << "Test finished.";
return EXIT_SUCCESS;
}
如有问题,欢迎指正