3D骨架提取

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;
}

如有问题,欢迎指正

  • 5
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 9
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值