读取VTK ImageData 和RT dose 计算DVH,储存表格

#include <vtkSphereSource.h>
#include <vtkImageMathematics.h>
#include <vtkCutter.h>
#include <vtkPlane.h>
#include <vtkStripper.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkLinearExtrusionFilter.h>
#include <vtkMetaImageWriter.h>
#include <vtkMetaImageReader.h>
#include <vtkImplicitModeller.h>
#include <vtkOutlineFilter.h>
#include <vtkImageGaussianSmooth.h>
#include <vtkAppendPolyData.h>
#include <vtkImageImport.h>
#include <vtkImageMapper.h>
#include <vtkLineSource.h>
#include <vtkFloatArray.h>
#include <vtkChartXY.h>
#include <vtkAxis.h>
#include <vtkConeSource.h>
#include <vtkContourFilter.h>
#include <vtkBrush.h>
#include <vtkPen.h>
#include <vtkPlot.h>
#include <vtkTextProperty.h>
#include <vtkTable.h>
#include <vtkFloatArray.h>
#include <vtkChartXY.h>
#include <vtkContextView.h>
#include <vtkChartLegend.h>

#include "DvhCalculate.h"
#include "stdafx.h"
#include "DataLayer.h"
#include "algorithm.h"
#include "ReadDoseData.h"
#include "DatabaseDefine.h"
#include "MedicineDatabase.h"
#include "PlanDataDefine.h"
#include "CustomInteractorStyle.h" 
#include "PluginDefine.h"
#include "LogUtil.h"
#include "RTDataModelDefine.h"

void DvhCalculate::SetDatabase(HP_TPS::IHPTPSRTDatabase* database)
{
	_database = database;
}

void  DvhCalculate::SetRTDataModel(HP_TPS::RTDataModel* rtDataModel)
{
	_rtDataModel = rtDataModel;
}

vtkSmartPointer<vtkTable> DvhCalculate::DvhChartCalculate(int planID) {
	//由planId获取plan
	HP_TPS::Plan plan;
	_database->GetPlanDatabase()->GetPlanById(planID, plan);
	//打印日志查看planId和plan.id一致性
	_logUtil.LogError(("plan.id:" + std::to_string(plan.id)).c_str());
	// 查询勾画集下的所有 ROI
	std::vector<HP_TPS::ROI> roiVec;
	_database->GetStructureDatabase()->QueryROIByStructureSetId(plan.referenceStructureSetId, roiVec);
	_logUtil.LogError(("referenceStructureSetId:" + std::to_string(plan.referenceStructureSetId)).c_str());
	// 用于存储满足条件的 roi.name
	std::vector<std::string> realroiNames;
	std::vector<std::vector<int>> realroiColors;
	std::vector<std::vector<int>> realroiID;
	// 遍历 roiVec 并存储非空的 roi.name
	if (!roiVec.empty()) {  //如果StructureSetId是非空就计算真实ROI
		for (const auto& roi : roiVec) {
			std::string roiName = roi.name;
			if (int(roi.isEmpty) == 0) {
				realroiNames.push_back(roiName); // 将满足条件的 roi.name 存储到 realroiNames 中
				realroiID.push_back({ roi.id });
				realroiColors.push_back({ roi.color.r, roi.color.g, roi.color.b });
				_logUtil.LogError(("realroiNames:" + roiName).c_str());
			}
		}
	}


	// 读取所有的 ROIImageData,为下一步计算 ROI 的 DVH 做准备
	std::unordered_map<int, vtkSmartPointer<vtkImageData>> RoisImageData = _rtDataModel->imageDataSource.imageDataROIMap;
	int dimX, dimY, dimZ;
	if (!RoisImageData.empty()) {
		if (RoisImageData.find(realroiID[0][0]) == RoisImageData.end()) {
			_logUtil.LogError(("Error: ROI not found in imageDataROIs: " + realroiNames[0]).c_str());
		}
		else {
			vtkSmartPointer<vtkImageData> ROIImageData = RoisImageData[realroiID[0][0]];//  保护
//获取ROI的尺寸
			int dimensions[3];
			ROIImageData->GetDimensions(dimensions);
			dimX = dimensions[0];
			dimY = dimensions[1];
			dimZ = dimensions[2];
			_logUtil.LogError(("ROIXsize:" + std::to_string(dimX)).c_str());
			_logUtil.LogError(("ROIYsize:" + std::to_string(dimY)).c_str());
			_logUtil.LogError(("ROIZsize:" + std::to_string(dimZ)).c_str());
		}
	}
	else {
		_logUtil.LogError("Error: ROIImageData is null.");
	}

	//积分微分的DVH表
	vtkSmartPointer<vtkTable> tableInt = vtkSmartPointer<vtkTable>::New();
	vtkSmartPointer<vtkTable> tablediff = vtkSmartPointer<vtkTable>::New();
	//
	vtkSmartPointer<vtkFloatArray> doseArrayInt = vtkSmartPointer<vtkFloatArray>::New();
	doseArrayInt->SetName("Dose");
	tableInt->AddColumn(doseArrayInt);
	//
	vtkSmartPointer<vtkFloatArray> intDvhArray = vtkSmartPointer<vtkFloatArray>::New();
	intDvhArray->SetName("Cumulative DVH");
	tableInt->AddColumn(intDvhArray);
	// 读取总剂量并转为三维矩阵
	//Readdosedata& readdosedata = Readdosedata::getInstance();
	//保存plan下面的剂量个数
	int number;
	if (!plan.bbcBeamDoseVec.empty()) {
		if (!plan.bbcBeamDoseVec[0].doseVec.empty()) {
			number = plan.bbcBeamDoseVec[0].doseVec.size();//  保护
		}
		else
		{
			_logUtil.LogError("Error: bbcBeamDoseVec[0].doseVec is null.");
		}
	}
	else {
		_logUtil.LogError("Error: bbcBeamDoseVec is null.");
	}
	//读取剂量下的地址
	std::string Filepath1;
	if (!number == 0) {
		Filepath1 = plan.bbcBeamDoseVec[0].doseVec[number - 1].doseGridFilepath;//  保护
		_logUtil.LogError(("doseFilepath:" + Filepath1).c_str());
	}
	else {
		_logUtil.LogError("Error: bbcBeamDoseVec is null.");
	}

	//判断剂量类型是不是总剂量值。
	int dose_type = plan.bbcBeamDoseVec[0].doseVec[number - 1].doseType;
	_logUtil.LogError(("dose_type:" + std::to_string(dose_type)).c_str());
	//赋值剂量地址给实际地址进行读取
	std::string Filepath;
	if (dose_type == 5) {
		Filepath = Filepath1;
	}
	else {
		_logUtil.LogError("Error: dose_type is not 5.");
	}

	//std::string Filepath = "D:\\mcCard_theta_90_RbeDose_equivalentDose_TIME_30.bin";  //测试的dose地址
	// 获取剂量三维矩阵 doseMatrix
	Readdosedata reader;
	reader.ReadDosebin(Filepath, dimX, dimY, dimZ);
	std::vector<std::vector<std::vector<double>>> doseMatrix = reader.getDoseMatrix();
	if (!doseMatrix.empty()) {
		_logUtil.LogError(("doseMatrix.size:" + std::to_string(doseMatrix.size())).c_str());
		_logUtil.LogError(("doseMatrix.size:" + std::to_string(doseMatrix[0].size())).c_str());
		_logUtil.LogError(("doseMatrix.size:" + std::to_string(doseMatrix[0][0].size())).c_str());
	}
	else {
		_logUtil.LogError("Error: Filepath is null.");
	}

	// 获取剂量矩阵的最大值和最小值
	double minValue = std::numeric_limits<double>::max();
	double maxValue = std::numeric_limits<double>::min();

	if (!doseMatrix.empty()) {
		for (const auto& slice : doseMatrix) {
			for (const auto& row : slice) {
				for (const auto& value : row) {
					maxValue = std::max(maxValue, value);
					minValue = std::min(minValue, value);
				}
			}
		}
		//记录剂量的最大值最小值。
		_logUtil.LogError(("maxValue:" + std::to_string(maxValue)).c_str());
		_logUtil.LogError(("minValue:" + std::to_string(minValue)).c_str());
	}
	else {
		_logUtil.LogError("Error: doseMatrix is null.");
	}

	int XofDVH = static_cast<int>((maxValue + 1));
	// 确认读取的非空 ROI 的 vtkImageData 并计算 DVH
	if (!doseMatrix.empty() && !realroiNames.empty()) {
		for (size_t i = 0; i < realroiID.size(); ++i) {
			if (RoisImageData.find(realroiID[i][0]) != RoisImageData.end()) {
				vtkSmartPointer<vtkImageData> roiImage = RoisImageData[realroiID[i][0]];
				if (!roiImage) {
					_logUtil.LogError(("No vtkImageData found for ROI: " + realroiNames[i]).c_str());
				}
				else {
					// 尝试获取 ROI 的判定值
					unsigned char* pixelsSourceROI = static_cast<unsigned char*>(roiImage->GetScalarPointer(0, 0, 0));
					int pixel = static_cast<int>(*pixelsSourceROI);
					std::string msgStr(std::to_string(pixel));
					_logUtil.LogInfo(msgStr.c_str());
				}

				// 计算体积直方图(DVH)
				std::vector<double> roiDVH(XofDVH, 0.0);
				double roiVolume = 0;
				for (int k = 0; k < dimZ; k++) {
					for (int j = 0; j < dimY; j++) {
						for (int i = 0; i < dimX; i++) {
							unsigned char* pixelsSourceROI = static_cast<unsigned char*>(roiImage->GetScalarPointer(i, j, k));
							if (static_cast<int>(*pixelsSourceROI) == 0) {
								roiVolume++;
								int dose = static_cast<int>(doseMatrix[i][j][k]);
								roiDVH[dose]++;
							}
						}
					}
				}
				_logUtil.LogError(("realroiNames:" + realroiNames[i] + "创建roiDVH成功.").c_str());
				//_logUtil.LogError(("realroiNames:" + roiDVH[0]).c_str());
				// 计算百分比
				for (double& bin : roiDVH) {
					bin = (bin / roiVolume) * 100.0;
				}

				// 计算积分 DVH
				std::vector<double> roiIntDVH(XofDVH, 0.0);
				double cumulativeVolume = 0.0;
				for (int i = static_cast<int>(maxValue); i >= 0; --i) {
					cumulativeVolume += roiDVH[i];
					roiIntDVH[i] = cumulativeVolume;
				}
				_logUtil.LogError(("realroiNames:" + realroiNames[i] + "创建roiIntDVH成功.").c_str());
				// 积分 DVH 的 table
				for (size_t i = 0; i < roiDVH.size(); ++i) {
					vtkIdType id = tableInt->InsertNextBlankRow();
					tableInt->SetValue(id, 0, minValue + i);
					tableInt->SetValue(id, 1, roiIntDVH[i]);
				}

				//std::ofstream outFile("roiDVH_" + roiName + ".csv");
				//outFile << "Dose,BinValue\n";
				//for (size_t i = 0; i < roiIntDVH.size(); ++i) {
				//	outFile << (minValue + i) << "," << roiIntDVH[i] << "\n";
				//}
				//outFile.close();
				_logUtil.LogError(("realroiNames:" + realroiNames[i] + "创建tableInt成功.").c_str());

			}
		}
	}
	else {
		_logUtil.LogError("Error: doseMatrix or realroiNames is null.");
	}

	// 输出积分 DVH 表格
	std::string info = "Integral DVH Table:";

	vtkSmartPointer<vtkTable> tableIntdraw = vtkSmartPointer<vtkTable>::New();
	// ROI 的数量
	int numROIs = realroiNames.size();
	_logUtil.LogInfo(("ROI 的数量:" + std::to_string(numROIs)).c_str());
	// 每个 ROI 的 DVH 数量
	int numDVHs = XofDVH;
	_logUtil.LogInfo(("每个 ROI 的 DVH 数量:" + std::to_string(numDVHs)).c_str());

	// 添加 Dose 列
	vtkSmartPointer<vtkFloatArray> doseArray = vtkSmartPointer<vtkFloatArray>::New();
	doseArray->SetName("Dose");
	tableIntdraw->AddColumn(doseArray);

	if (tableInt) {
		// 为每个 ROI 添加一个新的列到 tableIntdraw
		for (int roiIndex = 0; roiIndex < numROIs; ++roiIndex) {
			// 构造每个 ROI 的列名,例如 "ROI_1_DVH"
			std::string columnName = "ROI_" + std::to_string(roiIndex) + "_DVH";

			// 创建一个新的 vtkFloatArray 来存储该 ROI 的 DVH 数据
			vtkSmartPointer<vtkFloatArray> roiDvhArray = vtkSmartPointer<vtkFloatArray>::New();

			// 设置该 vtkFloatArray 的名称为构造的列名
			roiDvhArray->SetName(columnName.c_str());

			// 将该 vtkFloatArray 添加到 tableIntdraw 作为一列
			tableIntdraw->AddColumn(roiDvhArray);
		}

		// 填充 tableIntdraw 的数据
		for (int i = 0; i < numDVHs; ++i) {
			// 在 tableIntdraw 中插入一行空白行,并返回该行的 ID
			vtkIdType id = tableIntdraw->InsertNextBlankRow();

			// 从 m_IntegralDVHTable 中获取对应行的剂量值(第 0 列)
			double doseValue = tableInt->GetValue(i, 0).ToDouble();

			// 将剂量值设置到 tableIntdraw 的第 0 列的当前行
			tableIntdraw->SetValue(id, 0, vtkVariant(doseValue));

			// 遍历每个 ROI 并设置相应的 DVH 值
			for (int roiIndex = 0; roiIndex < numROIs; ++roiIndex) {
				// 从 m_IntegralDVHTable 中获取相应的 ROI 的 DVH 值
				// 假设 m_IntegralDVHTable 的列排列方式为: 第0列 - 剂量值, 第1列 - ROI1 DVH, 第2列 - ROI2 DVH, ...
				double roiDvhValue = tableInt->GetValue(roiIndex * numDVHs + i, 1).ToDouble();

				// 将该 ROI 的 DVH 值设置到 tableIntdraw 的相应列的当前行
				tableIntdraw->SetValue(id, roiIndex + 1, vtkVariant(roiDvhValue));
			}
		}
		_logUtil.LogInfo("创建tableIntdraw成功");
	}
	else {
		_logUtil.LogError("Error: tableInt is null.");
	}
	return tableIntdraw;
}

  • 5
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

大孔隆

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值