肺结节 最长径计算 切缘球的计算和生成

1、计算网格对应的mask

2、mask外扩

3、重新生成网格

代码

头文件

#ifndef YCUTTINGEDGEBALL_H
#define YCUTTINGEDGEBALL_H
#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkSmartPointer.h>
#include <vector>
#include <array>
double CalculateLongestDiameter(vtkPolyData* polyData,
                                double* p1,double* p2);
vtkImageData *fillBall(vtkImageData *mask, int len, int center[3]);
std::vector<std::array<int,3>> insideBall(
    vtkImageData *mask,int len, int center[3]);

class YCuttingEdgeBall
{
public:
    YCuttingEdgeBall();
    void caculateBall(vtkPolyData* poly);
    vtkSmartPointer<vtkPolyData> extend(double len);
    double voxelSize=0.3;
    vtkSmartPointer<vtkPolyData> polyData;
    vtkSmartPointer<vtkPolyData> polyData_dia;
    vtkSmartPointer<vtkPolyData> polyData_20;
    double longestdiameter=0;
    double p1[3];
    double p2[3];
};

#endif // YCUTTINGEDGEBALL_H

实现代码文件

#include "ycuttingedgeball.h"
#include <vtkMath.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkImageData.h>
#include <vtkDiscreteMarchingCubes.h>
#include <vtkKdTree.h>
#include <cmath>

// 计算两点之间的距离
static inline double CalculateDistance(const double* point1, const double* point2) {
    return std::sqrt(vtkMath::Distance2BetweenPoints(point1, point2));
}

// 计算mesh的最长径
double CalculateLongestDiameter(vtkPolyData* polyData,double* p1,double* p2) {
    double longestDiameter2 = 0.0;
    vtkPoints* points = polyData->GetPoints();
    int numPoints = points->GetNumberOfPoints();

    int idx1=0;
    int idx2=0;

    for (int i = 0; i < numPoints; ++i) {
        double point1[3];
        points->GetPoint(i, point1);

        for (int j = i + 1; j < numPoints; ++j) {
            double point2[3];
            points->GetPoint(j, point2);

            double distance2 = vtkMath::Distance2BetweenPoints(point1, point2);
            if (distance2 > longestDiameter2) {
                idx1=i;
                idx2=j;
                longestDiameter2 = distance2;
            }
        }
    }
    points->GetPoint(idx1, p1);
    points->GetPoint(idx2, p2);
    double longestDiameter=std::sqrt(longestDiameter2);
    return longestDiameter;
}
vtkImageData *fillBall(vtkImageData *mask, int len, int center[])
{
    int dims[3];
    mask->GetDimensions(dims);
    int64_t r2=len*len;
    int ib[3];
    ib[0]=center[0]>len?center[0]-len:0;
    ib[1]=center[1]>len?center[1]-len:0;
    ib[2]=center[2]>len?center[2]-len:0;
    int imax[3];
    imax[0]=center[0]+len+1<dims[0]?center[0]+len+1:dims[0];
    imax[1]=center[1]+len+1<dims[1]?center[1]+len+1:dims[1];
    imax[2]=center[2]+len+1<dims[2]?center[2]+len+1:dims[2];
    for(auto i=ib[2];i<imax[2];i++)
    {
        for(auto j=ib[1];j<imax[1];j++)
        {
            uint8_t* pdata=(uint8_t*)mask->GetScalarPointer(0,j,i);
            for(auto k=ib[0];k<imax[0];k++)
            {
                int64_t x=k-center[0];
                int64_t y=j-center[1];
                int64_t z=i-center[2];
                int64_t len2= x*x+y*y+z*z;
                if(len2<=r2)
                {
                    pdata[k]=1;
                }
            }
        }
    }
    return mask;
}
std::vector<std::array<int,3>> insideBall(vtkImageData *mask, int len, int center[])
{
    std::vector<std::array<int,3>> points;
    int dims[3];
    mask->GetDimensions(dims);
    int64_t r2=len*len;
    int ib[3];
    ib[0]=center[0]>len?center[0]-len:0;
    ib[1]=center[1]>len?center[1]-len:0;
    ib[2]=center[2]>len?center[2]-len:0;
    int imax[3];
    imax[0]=center[0]+len+1<dims[0]?center[0]+len+1:dims[0];
    imax[1]=center[1]+len+1<dims[1]?center[1]+len+1:dims[1];
    imax[2]=center[2]+len+1<dims[2]?center[2]+len+1:dims[2];
    for(auto i=ib[2];i<imax[2];i++)
    {
        for(auto j=ib[1];j<imax[1];j++)
        {
            uint8_t* pdata=(uint8_t*)mask->GetScalarPointer(0,j,i);
            for(auto k=ib[0];k<imax[0];k++)
            {
                int64_t x=k-center[0];
                int64_t y=j-center[1];
                int64_t z=i-center[2];
                int64_t len2= x*x+y*y+z*z;
                if(len2<=r2)
                {
                    if(!pdata[k])
                    {
                        points.emplace_back(std::array<int,3>{k,j,i});
                    }

                }
            }
        }
    }
    return points;
}
YCuttingEdgeBall::YCuttingEdgeBall() {}

void YCuttingEdgeBall::caculateBall(vtkPolyData *poly)
{
    polyData=poly;
    longestdiameter=CalculateLongestDiameter(poly,p1,p2);
#pragma omp parallel // starts a new team
    {
#pragma omp sections // divides the team into sections
        {
#pragma omp section
            { polyData_20=extend(20);;
            }
#pragma omp section
            {polyData_dia=extend(longestdiameter);}
        }
    }
}

vtkSmartPointer<vtkPolyData> YCuttingEdgeBall::extend(double len)
{
    vtkSmartPointer<vtkPolyData> polyex;
    double bound[6];
    polyData->ComputeBounds();
    polyData->GetBounds(bound);
    double origin[3];
    int dims[3];
    int ilen=len/voxelSize;
    int offset=ilen+20;
    for(auto i=0;i<3;i++)
    {
        dims[i]=(bound[2*i+1]-bound[2*i])/voxelSize+offset*2+2;
        origin[i]=bound[2*i]-offset*voxelSize;
    }

    vtkNew<vtkImageData> mask;
    mask->SetDimensions(dims);
    mask->SetSpacing(voxelSize,voxelSize,voxelSize);
    mask->SetOrigin(origin);
    mask->AllocateScalars(VTK_UNSIGNED_CHAR,1);
    memset(mask->GetScalarPointer(),0,dims[0]*dims[1]*dims[2]);

    int ijk[3];
    for(int j=0;j<3;j++)
    {
        ijk[j]=(p1[j]-bound[2*j])/voxelSize+offset;
    }

    fillBall(mask,ilen,ijk);
    for(int j=0;j<3;j++)
    {
        ijk[j]=(p2[j]-bound[2*j])/voxelSize+offset;
    }
    fillBall(mask,ilen,ijk);
    std::vector<std::array<int,3>> vpoints;
    {
        for(int j=0;j<3;j++)
        {
            double p = (p2[j] + p1[j]) / 2;
            ijk[j]=(p-bound[2*j])/voxelSize+offset;
        }
        int dis = std::floor(CalculateDistance(p1, p2)/2/voxelSize)+1;
        int r = ilen + dis;
        vpoints =insideBall(mask,r,ijk);
    }
    vtkNew<vtkKdTree> dkdtree;
    // 构建vtkKdTree
    vtkKdTree* kdTree = dkdtree;
    auto points=polyData->GetPoints();
    int numpoint = points->GetNumberOfPoints();
    vtkNew<vtkPoints> newpoint;
    newpoint->SetNumberOfPoints(numpoint);
    for (int i = 0; i < numpoint; i++)
    {
        double pt[3];
        points->GetPoint(i,pt);
        pt[0] = (pt[0] - origin[0]) / voxelSize;
        pt[1] = (pt[1] - origin[1]) / voxelSize;
        pt[2] = (pt[2] - origin[2]) / voxelSize;
        newpoint->SetPoint(i, pt);
    }

    kdTree->BuildLocatorFromPoints(newpoint);
    vtkIdTypeArray* idarray = kdTree->BuildMapForDuplicatePoints(1);
    idarray->UnRegister(nullptr);
    double len2 = len * len / (voxelSize * voxelSize);
    uint8_t* pdata = (uint8_t*)mask->GetScalarPointer();
    if (vpoints.size() > 500000)
    {
#pragma omp parallel for num_threads(4)
        for (int i = 0; i < vpoints.size(); i++)
        {
            auto& pt = vpoints[i];
            double p[3];
            p[0] = pt[0];// *precision + origin[0];
            p[1] = pt[1];// *precision + origin[1];
            p[2] = pt[2];// *precision + origin[2];
            double dist2; // 平方距离
            // 使用vtkKdTree查找最近点
            kdTree->FindClosestPoint(p[0], p[1], p[2], dist2);
            if (dist2 < len2)
            {
                //pdata[pt[0]+pt[1]*dims[0]+pt[2]*dims[0]*dims[1]]=1;
                mask->SetScalarComponentFromDouble(pt[0], pt[1], pt[2], 0, 1);
            }
        }
    }
    else
    {
        for (auto pt : vpoints)
        {
            double p[3];
            p[0] = pt[0];// * precision + origin[0];
            p[1] = pt[1];// * precision + origin[1];
            p[2] = pt[2];// * precision + origin[2];
            double dist2; // 平方距离
            // 使用vtkKdTree查找最近点
            kdTree->FindClosestPoint(p[0], p[1], p[2], dist2);
            if (dist2 < len2)
            {
                //pdata[pt[0]+pt[1]*dims[0]+pt[2]*dims[0]*dims[1]]=1;
                mask->SetScalarComponentFromDouble(pt[0], pt[1], pt[2], 0, 1);
            }
        }
    }
    vtkSmartPointer<vtkDiscreteMarchingCubes> marchingCubes = vtkSmartPointer<vtkDiscreteMarchingCubes>::New();
    marchingCubes->SetInputData(mask);
    marchingCubes->SetValue(0, 1);   //设置等值线
    marchingCubes->ComputeNormalsOn();   //将法向量计算开启
    marchingCubes->ComputeScalarsOff();  //如果计算Scalars 会导致mesh法线错误
    marchingCubes->Update();
    polyex=marchingCubes->GetOutput();

    return polyex;
}