Plan-Cut on segment's labelmap

Support
1

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

2026-01-09_214810
2026-01-09_214810722×394 61.1 KB

2026-01-09_214856
2026-01-09_214856719×392 42.6 KB

2

May you consider this tentative change?

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

→ 

double dotProduct =
  rasCenter[0] * planeOrigin_World[0] +
  rasCenter[1] * planeOrigin_World[1] +
  rasCenter[2] * planeOrigin_World[2];

I have not test run anything.

3

I doubt the above modification would solve anything after second thought.

This line seems to be considered:

double ijkCenter[3] = { x + 0.5 ,y + 0.5 , z + 0.5 };

x, y and z are positive integers of an IJK coordinate. It should perhaps be transformed to RAS first, and add half of the spacing on each axis to get the center RAS coordinate of each voxel.