Plane-Cut on segment's labelmap

[shanfl]

I need to plane-cut a labelmap of a segment with a vtkplanewidget, but the cutting position is not correct.I think there might be an issue with the conversion from voxel ijk coordinates to RAS coordinates, but I can’t pinpoint the cause.

here is code

std::vector<vtkSmartPointer<vtkOrientedImageData>> PlaneCutImage(vtkOrientedImageData* image, double planeOrigin_World[3], double planeNormal_World[3])
{
	Q_D(qSlicerSurgeryPlanPlaneCutEffect);
	if (!image) {
		return {};
	}

	qDebug() << "planeOrigin:" << planeOrigin_World[0] << " " << planeOrigin_World[1] << " " << planeOrigin_World[2];
	qDebug() << "planeNormal:" << planeNormal_World[0] << " " << planeNormal_World[1] << " " << planeNormal_World[2];


	// ijk 2 ras
	vtkNew<vtkTransform> ijkToWorldTransform;
	vtkNew<vtkMatrix4x4> ijkToWorldMatrix;
	image->GetImageToWorldMatrix(ijkToWorldMatrix);
	ijkToWorldTransform->SetMatrix(ijkToWorldMatrix);

	PrintMatrixFormatted(ijkToWorldMatrix);



	//
	int dims[3];
	double spacing[3], origin[3];
	image->GetDimensions(dims);
	image->GetSpacing(spacing);
	image->GetOrigin(origin);

	// 
	if (dims[0] <= 0 || dims[1] <= 0 || dims[2] <= 0) {
		return {};
	}



	// 
	int extent[6] = { 0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1 };
	void* wholeInputPtr = image->GetScalarPointerForExtent(extent);
	if (!wholeInputPtr) {
		wholeInputPtr = image->GetScalarPointer();
		if (!wholeInputPtr) {
			return {};
		}
	}

	int scalarType = image->GetScalarType();
	int numComponents = image->GetNumberOfScalarComponents();
	size_t scalarSize = image->GetScalarSize();

	if (scalarSize == 0) {
		return {};
	}

	size_t voxelSize = scalarSize * numComponents;

	std::vector<vtkSmartPointer<vtkOrientedImageData>> ret;

	// 
	vtkNew<vtkOrientedImageData> part1;
	vtkNew<vtkOrientedImageData> part2;
	part1->DeepCopy(image);
	part2->DeepCopy(image);



	// 
	part1->AllocateScalars(scalarType, numComponents);
	part2->AllocateScalars(scalarType, numComponents);

	//
	void* wholePart1Ptr = part1->GetScalarPointerForExtent(extent);
	void* wholePart2Ptr = part2->GetScalarPointerForExtent(extent);

	if (!wholePart1Ptr || !wholePart2Ptr) {
		wholePart1Ptr = part1->GetScalarPointer();
		wholePart2Ptr = part2->GetScalarPointer();

		if (!wholePart1Ptr || !wholePart2Ptr) {
			return {};
		}
	}


	vtkIdType totalBytes = static_cast<vtkIdType>(dims[0]) * dims[1] * dims[2] * voxelSize;
	memset(wholePart1Ptr, 0, totalBytes);
	memset(wholePart2Ptr, 0, totalBytes);

	bool sideWrite[2] = { false, false };

	// handle every voxel
	for (int z = 0; z < dims[2]; ++z) {
		for (int y = 0; y < dims[1]; ++y) {
			for (int x = 0; x < dims[0]; ++x) {

				//
				double ijkCenter[3] = { x + 0.5 ,y + 0.5 , z + 0.5 };
				double rasCenter[3];
				ijkToWorldTransform->TransformPoint(ijkCenter, rasCenter);

				if ((y == 0) && (z == 0) && (x == 0)) {
					qDebug() << "IJK (0,0,0) -> RAS:" << rasCenter[0] << rasCenter[1] << rasCenter[2];
				}



				double vector[3] = {
					rasCenter[0] - planeOrigin_World[0],
					rasCenter[1] - planeOrigin_World[1],
					rasCenter[2] - planeOrigin_World[2]
				};

				double dotProduct =
					vector[0] * planeNormal_World[0] +
					vector[1] * planeNormal_World[1] +
					vector[2] * planeNormal_World[2];


				// 
				size_t offset = (z * dims[1] * dims[0] + y * dims[0] + x) * voxelSize;

				void* inputPtr = static_cast<char*>(wholeInputPtr) + offset;
				void* part1Ptr = static_cast<char*>(wholePart1Ptr) + offset;
				void* part2Ptr = static_cast<char*>(wholePart2Ptr) + offset;

				if (dotProduct >= 0) {
					sideWrite[0] = true;
					memcpy(part1Ptr, inputPtr, voxelSize);
					
				}
				else {
					sideWrite[1] = true;
					memcpy(part2Ptr, inputPtr, voxelSize);
					
				}
			}
		}
	}

	
	part1->SetImageToWorldMatrix(ijkToWorldMatrix);
	part2->SetImageToWorldMatrix(ijkToWorldMatrix);

	//
	if (sideWrite[0]) {
		ret.emplace_back(part1);

		WriteNII("D:/part1.nii.gz", part1);
	}
	if (sideWrite[1]) {
		ret.emplace_back(part2);
		WriteNII("D:/part2.nii.gz", part2);
	}

	return ret;
}