Smoothing Surface crashing

[ilandaniele (Ilan)]

Good afternoon.

In my project i convert a vtkVolume by some given DICOM files.
Then i convert my vtkVolume on to a Surface by using vtkMarchingCubes and 6 steps more.
Code is the following (myVolume is a class that inherits from vtkVolume):

void VolumeToSurface::ConvertVolume(vtkSmartPointer vol, double r)
{
if (vol == nullptr)
return;

vol->setPicked(false);

double ratio = 1.0 - r;

bool finish = false;

vtkSmartPointer<vtkPlaneCollection> planes = vol->GetMapper()->GetClippingPlanes();

int min = -1000;
int max = -1000;

double R, G, B;


//convertToPoints(vol, 0.15, 200, 1800);

while(!finish){
	
	finish = extraerRango(min, max, vol);

	if (min < 10000) {
		Consola::append("Creando superficie con rango ");
		Consola::append(min);
		Consola::append(" - ");
		if(max < 10000)
			Consola::append(max);
		else
			Consola::append("Max");

		Consola::append("\n");

		imageData = vtkSmartPointer<vtkImageData>::New();
		surfaceExtractor = vtkSmartPointer<vtkMarchingCubes>::New();

		imageData->DeepCopy(vol->getImageData());



		surfaceExtractor->ComputeNormalsOn();
		
		Consola::append("Etapa 1 de 7: Filtrando valores...");
		if (abs(min) > abs(max)) {
			
			
			surfaceExtractor->SetValue(max, min);
			
			Consola::appendln("Salteada");
		}
		else {

			aplicarRango(vol->getImageData(), min, max, imageData);
			Consola::appendln("OK");
			surfaceExtractor->SetValue(0, 128);
		}

		//marchingcubes

		surfaceExtractor->SetInputData(imageData);

		vtkSmartPointer<vtkCallbackCommand> callback = CallbacksClass::getProgressCallBack();
		if (callback != nullptr)
		{
			CallbacksClass::setEncabezado("Generando superficie: ");
			surfaceExtractor->AddObserver(vtkCommand::ProgressEvent, callback);
		}

		Consola::append("Etapa 2 de 7: Generando superficie...");
		surfaceExtractor->Update();
		Consola::appendln("OK");
		
		//marchingcubes

		vtkNew<vtkPolyData> mesh; //reemplazar por vtkSmartPointer
		mesh->DeepCopy(surfaceExtractor->GetOutput());
		
		
		//triangulate
		vtkSmartPointer<vtkTriangleFilter> triangleFilter =
			vtkSmartPointer<vtkTriangleFilter>::New();
		triangleFilter->SetInputData(mesh);

		if (callback != nullptr)
		{
			CallbacksClass::setEncabezado("Triangulando caras: ");
			triangleFilter->AddObserver(vtkCommand::ProgressEvent, callback);
		}
		Consola::append("Etapa 3 de 7: Triangulando caras...");
		triangleFilter->Update();
		mesh->DeepCopy(triangleFilter->GetOutput());
		Consola::appendln("OK");
		//triangulate
		

		//decimator
		
		Consola::append("Etapa 4 de 7: Reduciendo caras...");

		if (ratio > 0.01) {
			vtkSmartPointer<vtkQuadricDecimation> decimator = vtkSmartPointer<vtkQuadricDecimation>::New();
			decimator->SetInputData(mesh);
			decimator->SetTargetReduction(ratio);

			if (callback != nullptr)
			{
				CallbacksClass::setEncabezado("Reduciendo caras: ");
				decimator->AddObserver(vtkCommand::ProgressEvent, callback);
			}

			decimator->Update();

			mesh->DeepCopy(decimator->GetOutput());
			Consola::appendln("OK");
		}
		else{

			Consola::appendln("Salteado");
		}

		//decimator
		

		

		//clean mesh
		
		vtkSmartPointer<vtkCleanPolyData> cleanPolyData =
			vtkSmartPointer<vtkCleanPolyData>::New();
		cleanPolyData->SetInputData(mesh);
		cleanPolyData->SetTolerance(0.001);

		if (callback != nullptr)
		{
			CallbacksClass::setEncabezado("Limpiando puntos: ");
			cleanPolyData->AddObserver(vtkCommand::ProgressEvent, callback);
		}
		Consola::append("Etapa 5 de 7: Limpiando puntos...");

		cleanPolyData->Update();
		Consola::appendln("OK");
		mesh->DeepCopy(cleanPolyData->GetOutput());
		
		
		
		//making slices

		vtkSmartPointer<vtkBox> implicitCube =
			vtkSmartPointer<vtkBox>::New();
		
		double* bnds = vol->getCustomBounds();
		double* origBounds = vol->getImageData()->GetBounds();
		//printf("BOUNDS %f\n", bnds[0]);

		double bounds[6];

		bounds[0] = bnds[0] + origBounds[1] / 2.0;
		bounds[1] = bnds[1] + origBounds[1] / 2.0;
		bounds[2] = bnds[2] + origBounds[3] / 2.0;
		bounds[3] = bnds[3] + origBounds[3] / 2.0;
		bounds[4] = bnds[4] + origBounds[5] / 2.0;
		bounds[5] = bnds[5] + origBounds[5] / 2.0;

		implicitCube->SetBounds(bounds);

		vtkSmartPointer<vtkClipPolyData> clipper =
			vtkSmartPointer<vtkClipPolyData>::New();
		clipper->SetClipFunction(implicitCube);

		clipper->SetInputData(mesh);

		clipper->InsideOutOn();
		Consola::append("Etapa 6 de 7: Realizando cortes...");
		if (callback != nullptr)
		{
			CallbacksClass::setEncabezado("Realizando cortes: ");
			clipper->AddObserver(vtkCommand::ProgressEvent, callback);
		}
		clipper->Update();
		
		Consola::appendln("OK");
		mesh->DeepCopy(clipper->GetOutput());
		

		MasterHandler::GetSurfaceInteractor()->GetMeshOperations()->moveMeshToCOSCenter(mesh);
	

		//making slices

		//remove small objects

		Consola::appendln("Etapa 7 de 7: Removiendo objetos sueltos");

		if (!vol->getIsSplitted()) {

			vtkSmartPointer<vtkPolyDataConnectivityFilter> connectivityFilter = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
			connectivityFilter->SetInputData(mesh);
			connectivityFilter->SetExtractionModeToAllRegions();
			vtkSmartPointer<vtkCallbackCommand> m_progressCallback = CallbacksClass::getProgressCallBack();

			CallbacksClass::setEncabezado("Analizando objeto: ");
			connectivityFilter->AddObserver(vtkCommand::ProgressEvent, m_progressCallback);
			Consola::appendln("...Etapa 1 de 2: Analizando objeto...");
			connectivityFilter->Update();

			// remove objects consisting of less than ratio vertexes of the biggest object
			vtkIdTypeArray* regionSizes = connectivityFilter->GetRegionSizes();
			Consola::appendln("OK");

			// find object with most vertices
			long maxSize = 0;
			for (int regions = 0; regions < connectivityFilter->GetNumberOfExtractedRegions(); regions++)
				if (regionSizes->GetValue(regions) > maxSize)
					maxSize = regionSizes->GetValue(regions);


			// append regions of sizes over the threshold
			connectivityFilter->SetExtractionModeToSpecifiedRegions();
			for (int regions = 0; regions < connectivityFilter->GetNumberOfExtractedRegions(); regions++)
				if (regionSizes->GetValue(regions) > maxSize * 0.02)
					connectivityFilter->AddSpecifiedRegion(regions);

			Consola::appendln("...Etapa 2 de 2: Removiendo objetos...");

			CallbacksClass::setEncabezado("Removiendo objetos: ");
			connectivityFilter->Update();

			
			//remove small objects
			mesh->DeepCopy(connectivityFilter->GetOutput());
			Consola::appendln("OK");
		}
		else {
			Consola::appendln("Salteada");

		}
		vtkNew<vtkPolyDataMapper> mapper;
		mapper->SetInputData(mesh);
		mapper->ScalarVisibilityOff();
		vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
		actor->SetMapper(mapper);

		
		if(max<10000)
			extraerColor(R, G, B, (max - min) / 2, vol);
		else
			extraerColor(R, G, B, (1500 - min) / 2, vol);
		actor->GetProperty()->SetColor(R, G, B);

		renwinVTS->GetRenderers()->GetFirstRenderer()->AddActor(actor);
		emit MasterHandler::GetSimpleView()->signal_renderizar();
		Consola::append("Finalizado\n");
	}
	
}

}

bool VolumeToSurface::extraerRango(int & min, int & max, vtkSmartPointer vol)
{
QVector X_function = *vol->GetOpacityFunction_X();
QVector Y_function = *vol->GetOpacityFunction_Y();

int cant = NUMERODEPUNTOS_SCALAR;

int i = 0;

while (i < cant && X_function[i] < max) {
	i++;
}

while (i < cant && Y_function[i] < 0.9) {
	i++;
}

if (i < cant)
	min = X_function[i];
else {
	min = 50000;
	max = 50000;
	return true;
}

while (i < cant && Y_function[i] >= 0.9) {
	i++;
}

if (i < cant) {
	max = X_function[i];
	return false;
}
else {
	max = 50000;
	return true;
}

}

void VolumeToSurface::aplicarRango(vtkSmartPointer input, int min, int max, vtkSmartPointer output)
{

emit  MasterHandler::GetSimpleView()->visibleProgressBar(true);
int* dims = output->GetDimensions();
int recorridas = 0;
int numCmp = dims[0] * dims[1] * dims[2];

short* input_ptr = static_cast<short*>(input->GetScalarPointer());
short* output_ptr = static_cast<short*>(output->GetScalarPointer());

int dim_x = dims[0];
int dim_y = dims[1];
int dim_z = dims[2];

int update = dims[2] / 20.0;
int iternum = 0;

for (int k = 0; k < dim_z; k++) {
	for (int i = 0; i < dim_x; i++) {
		for (int j = 0; j < dim_y; j++) {

			//double val = input->GetScalarComponentAsDouble(i, j, k, 0);
			double val = GetScalarValue(input_ptr, i, j, k, dim_x, dim_y);

			if (val < min || val > max) {
				/**
				if (neighbourInRange(input, i, j, k, dims, min, max)) {
					val = 255;
				}
				else/**/
					val = 0;
				
			}
			else {
				val = 255;
			}
			//output->SetScalarComponentFromDouble(i, j, k, 0, val);
			SetScalarValue(output_ptr, i, j, k, dim_x, dim_y, val);
			recorridas++;
		}
		
		//progress->setValue((int)(x*100.0f / (double)dim_x));
		if (iternum == 0)
		{
			emit MasterHandler::GetSimpleView()->progressChanged(recorridas * 100.0f / (double)numCmp);
			emit MasterHandler::GetSimpleView()->signal_WriteLabel(QString(
				QString("Filtrando valores: ")
				+ QString::number((int)(recorridas * 100.0f / (double)numCmp))
				+ QString("%")
			));
		}
		iternum = (iternum + 1) % update;
		
	}
}
emit  MasterHandler::GetSimpleView()->visibleProgressBar(false);

}

void VolumeToSurface::extraerColor(double & R, double & G, double & B, int index, vtkSmartPointer vol)
{

QVector<double> Color_function = *vol->GetColorFunction();



int i = (index + 1000) / STEP_COLOR;

R = Color_function[i * 3 + 0];
G = Color_function[i * 3 + 1];
B = Color_function[i * 3 + 2];

}

It generates the surface successfully. However on part 5 of cleanPolyData, it takes a little bit of a long time.

Next if i press smooth mesh button on my app, it runs the next code:

void MeshOperations::smoothMesh( vtkSmartPointer actor, unsigned int iterations )
{

vtkSmartPointer<vtkPolyData> mesh = (vtkPolyData*)actor->GetMapper()->GetInput();
//cout << "Mesh smoothing with " << nbrOfSmoothingIterations << " iterations." << endl;
Consola::appendln("Suavizando objeto");

vtkSmartPointer<vtkSmoothPolyDataFilter> smoother = vtkSmartPointer<vtkSmoothPolyDataFilter>::New();
smoother->SetInputData( mesh );
smoother->SetNumberOfIterations( int(iterations) );
smoother->SetFeatureAngle(45);
smoother->SetRelaxationFactor(0.05);

vtkSmartPointer<vtkCallbackCommand> m_progressCallback = CallbacksClass::getProgressCallBack();

{
	CallbacksClass::setEncabezado("Suavizando objeto: ");
    smoother->AddObserver(vtkCommand::ProgressEvent, m_progressCallback);
}
smoother->Update();

vtkSmartPointer<vtkPolyDataMapper> polyMapper = (vtkPolyDataMapper*)actor->GetMapper();

//lock mutex
MasterHandler::GetThread()->lock_Actor();
emit MasterHandler::GetSimpleView()->signal_RemoveActor(actor);

//lock mutex
MasterHandler::GetThread()->lock_Actor();
mesh->DeepCopy( smoother->GetOutput() );

polyMapper->SetInputData(mesh);

polyMapper->ScalarVisibilityOff();

actor->SetMapper(polyMapper);

MasterHandler::GetSurfaceInteractor()->AddActor(actor);

//unlock mutex
MasterHandler::GetThread()->unlock_Actor();
Consola::appendln("Suavizando objeto: Finalizado");

}

Then, when it comes to an end, it suddenly stops and throws an error:
In vtkAbstractArray.h in line 180 GetNumberOfValues() says:
Read access violation on this was 0xFFFFFFFFFFFFFFC7

I know it has to do that it did not find some value on the array, but i do not understand where it is corrupted. This used to work on version VTK 8.1 but now on VTK 8.2 or VTK 9.0.0 is not working anymore. Does anyone know if any of this used vtk classes is not working properly anymore?
I tried debugging and it passed all the lines, i suspect it has to do with the smoother->Update() line.
Best regards.

[lassoan (Andras Lasso)]

It is useful that you have shared your source code, but it is so much that I don’t think anybody would have the time to review it or create an executable that can be tested in a debugger. Could you provide a minimal example that reproduces the problem?

[ilandaniele (Ilan)]

Thanks for the answer.
Sure.
As i said before, i load in my application a couple of dicom files, which generates a vtkVolume object, like the next one:

Then, i select the vtkVolume and transform it to a surface with one of my interactions, it follows 7 steps:

1. Filtering values (to set a threshold of values that i want to transform to surface, based on to the Hounsfield coefficients)
2. Generate the surface using vtkMarchingCubes
3. Triangulate the faces using vtkTriangleFilter
4. Reduce faces using vtkQuadricDecimation (starts taking a little bit of a longer time)
5. Cleaning the mesh using vtkCleanPolyData (i use a progress bar to see the progress of the update, dont know why but it suddenly stops at 50% and takes a long time to go)
6. Clipping the data so i only use surface that is of my interest, using vtkBox and vtkClipPolyData
7. Removing small objects from the surface, using vtkPolyDataConnectivityFilter (I suspect the error has to do with this part, since it uses a vtkIdTypeArray, which is the class giving me errors).

And it ends up with a surface like the next one:

Which really needs a smoothing.
Next, i press on my smooth interaction and it uses vtkSmoothPolyDataFilter to smooth (have tried using vtkWindowedSincPolyDataFilter but it gave me the same crash result).

It reaches 97% on the progress bar, and then it crashes and shows the next exception (sorry it is in spanish):

which means a read acces violation was produced

The strange thing is, that if i follow the same steps with a smaller volume or a “clipped before” volume, it works properly and does not crash.
Let me know if you need any further example.

[lassoan (Andras Lasso)]

This data processing pipeline should work well. We have been using these filters extensively without any problems for many years, so most likely the issue is how exactly you use these filters or something specific to your data or environment. If you can provide minimal source code and data reproduces the issue then we can help you investigate this issue further.

[ilandaniele (Ilan)]

Sure, thank you for the help.
I am uploading you here the files that are needed. I tried to make it minimal, but still some features have to be added so you can check if they affect the performance of the application.BonesViewer.rar (218.8 KB)
Take into account that i did the minimilization as best as i could since it is a big project, let me know if something fails or a part of code is missing, and i will fix it asap so you can check it.
I tried commenting steps 3 to 6 and it worked fine, but still i don’t know why it fails.
Thanks!

[ilandaniele (Ilan)]

I forgot to pass you the dataset that i use, so you can recreate the whole scenario. Since some features were taken out, maybe if i pass you the whole project it would be easier to recreate the crashing. Let me know if you need it also.
Dataset: https://www.mediafire.com/file/s7w34avj6or7u7h/S20.rar/file
Best regards.