Merge pull request ITISFoundation#115 from dyollb/feature/topology_fi…

…x_manifold

Feature/topology fix manifold

Core/Graph.h

@@ -7,6 +7,7 @@
 
 #include <array>
 #include <cmath>
+#include <type_traits>
 
 namespace XCore
 {
@@ -50,7 +51,7 @@ namespace XCore
 
 	public:
 		using idx_type = TIndex;
-		using idx_value = typename idx_type::IndexValueType;
+		using idx_value = typename std::decay<decltype(std::declval<idx_type&>()[0])>::type;
 		using idx_list_type = StackVector<idx_type, Pow(3, D) - 1>;
 
 		CUniformGridGraph(const std::array<idx_value, D>& dims, const std::array<double, D>& spacing)

Core/TopologyInvariants.h

@@ -30,7 +30,6 @@ bool EulerInvariant(const TNeighborhood& neighbors, const TLabel label)
 	// ==> If we mark the center as BG, then P, F, E and V are '0'.
 	static const int c = 27 / 2;
 
-//fprintf(stderr, "neighbors[c]: %d vs %d", static_cast<int>(neighbors[c]), static_cast<int>(label));
 	assert(neighbors[c] == label);
 
 	// Voxels (parallelepipeds) - count center
@@ -211,4 +210,28 @@ bool CCInvariant(TNeighborhood neighbors, const TLabel label)
 	return (cc_before == cc_after);
 }
 
+/** \brief Will voxels be manifold after removing center voxel?
+ */
+template<typename TNeighborhood, typename TLabel>
+bool NonmanifoldRemove(const TNeighborhood& neighbors, const TLabel label)
+{
+	static const int c = 27 / 2;
+
+	// test for non-manifold edges around center voxel
+	return (neighbors[c-1]==label && neighbors[c-3]==label && neighbors[c-1-3]!=label)
+		|| (neighbors[c+1]==label && neighbors[c-3]==label && neighbors[c+1-3]!=label)
+		|| (neighbors[c-1]==label && neighbors[c+3]==label && neighbors[c-1+3]!=label)
+		|| (neighbors[c+1]==label && neighbors[c+3]==label && neighbors[c+1+3]!=label)
+		//
+		|| (neighbors[c-1]==label && neighbors[c-9]==label && neighbors[c-1-9]!=label)
+		|| (neighbors[c+1]==label && neighbors[c-9]==label && neighbors[c+1-9]!=label)
+		|| (neighbors[c-1]==label && neighbors[c+9]==label && neighbors[c-1+9]!=label)
+		|| (neighbors[c+1]==label && neighbors[c+9]==label && neighbors[c+1+9]!=label)
+		//
+		|| (neighbors[c-3]==label && neighbors[c-9]==label && neighbors[c-3-9]!=label)
+		|| (neighbors[c+3]==label && neighbors[c-9]==label && neighbors[c+3-9]!=label)
+		|| (neighbors[c-3]==label && neighbors[c+9]==label && neighbors[c-3+9]!=label)
+		|| (neighbors[c+3]==label && neighbors[c+9]==label && neighbors[c+3+9]!=label);
+}
+
 } // namespace topology

Core/itkFixTopologyCarveOutside.h

@@ -81,6 +81,9 @@ class /*ITK_TEMPLATE_EXPORT*/ FixTopologyCarveOutside : public ImageToImageFilte
 	itkSetMacro(OutsideValue, InputImagePixelType);
 	itkGetConstMacro(OutsideValue, InputImagePixelType);
 
+	itkSetMacro(EnforceManifold, bool);
+	itkGetConstMacro(EnforceManifold, bool);
+
 	/** Get Skeleton by thinning image. */
 	OutputImageType* GetThinning();
 
@@ -116,6 +119,7 @@ class /*ITK_TEMPLATE_EXPORT*/ FixTopologyCarveOutside : public ImageToImageFilte
 private:
 	InputImagePixelType m_InsideValue = 1;
 	InputImagePixelType m_OutsideValue = 0;
+	bool m_EnforceManifold = true;
 
 }; // end of FixTopologyCarveOutside class
 

Core/itkFixTopologyCarveOutside.hxx

@@ -154,10 +154,13 @@ void FixTopologyCarveOutside<TInputImage, TOutputImage>::ComputeThinImage()
 			bool can_carve = false;
 			if (EulerInvariant(ot.GetNeighborhood(), m_InsideValue))
 			{
-				// Check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood)
-				if (CCInvariant(ot.GetNeighborhood(), m_InsideValue))
+				if (!m_EnforceManifold || !NonmanifoldRemove(ot.GetNeighborhood(), m_InsideValue))
 				{
-					can_carve = true;
+					// Check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood)
+					if (CCInvariant(ot.GetNeighborhood(), m_InsideValue))
+					{
+						can_carve = true;
+					}
 				}
 			}
 

iSeg/SurfaceViewerWidget.cpp

@@ -28,18 +28,18 @@
 #include <vtkPointData.h>
 #include <vtkUnsignedShortArray.h>
 
+#include <vtkDecimatePro.h>
 #include <vtkEventQtSlotConnect.h>
 #include <vtkInteractorStyleTrackballCamera.h>
 #include <vtkLookupTable.h>
+#include <vtkPolyDataConnectivityFilter.h>
 #include <vtkPolyDataMapper.h>
 #include <vtkPropPicker.h>
 #include <vtkProperty.h>
 #include <vtkQuadricLODActor.h>
 #include <vtkRenderWindow.h>
 #include <vtkRenderer.h>
 
-#include <vtkPolyDataConnectivityFilter.h>
-
 #include <vtkAutoInit.h>
 #ifdef ISEG_VTK_OPENGL2
 VTK_MODULE_INIT(vtkRenderingOpenGL2);
@@ -77,8 +77,7 @@ void transform_slices_vtk(const std::vector<TIn*>& slices, size_t slice_size, TA
 	}
 }
 
-enum eActions
-{
+enum eActions {
 	kSelectTissue,
 	kGotoSlice,
 	kMarkPoint
@@ -111,6 +110,21 @@ SurfaceViewerWidget::SurfaceViewerWidget(SlicesHandler* hand3D1, eInputType inpu
 	lb_connectivity_count = new QLabel("");
 	lb_connectivity_count->setToolTip("Number of connected regions");
 
+	reduction = new QLineEdit(QString::number(0));
+	reduction->setValidator(new QIntValidator(0, 100));
+	reduction->setToolTip("Reduction of number triangles in percent.");
+
+	// set default reduction depending number of voxels.
+	const auto N = static_cast<float>(hand3D->return_area()) * hand3D->num_slices();
+	if (N > 1e8)
+	{
+		reduction->setText(QString::number(80));
+	}
+	else if (N > 1e6)
+	{
+		reduction->setText(QString::number(50));
+	}
+
 	// layout
 	auto vbox = new QVBoxLayout;
 	vbox->addWidget(vtkWidget);
@@ -120,6 +134,10 @@ SurfaceViewerWidget::SurfaceViewerWidget(SlicesHandler* hand3D1, eInputType inpu
 	transparency_hbox->addWidget(sl_trans);
 	vbox->addLayout(transparency_hbox);
 
+	auto reduction_hbox = new QHBoxLayout;
+	reduction_hbox->addWidget(reduction);
+	vbox->addLayout(reduction_hbox);
+
 	if (input_type == kSource)
 	{
 		auto lb_thresh = new QLabel("Contour iso-value");
@@ -145,9 +163,10 @@ SurfaceViewerWidget::SurfaceViewerWidget(SlicesHandler* hand3D1, eInputType inpu
 	setLayout(vbox);
 
 	// connections
-	QObject::connect(sl_trans, SIGNAL(sliderReleased()), this, SLOT(transp_changed()));
-	QObject::connect(bt_update, SIGNAL(clicked()), this, SLOT(reload()));
-	QObject::connect(bt_connectivity, SIGNAL(clicked()), this, SLOT(split_surface()));
+	connect(sl_trans, SIGNAL(sliderReleased()), this, SLOT(transp_changed()));
+	connect(bt_update, SIGNAL(clicked()), this, SLOT(reload()));
+	connect(bt_connectivity, SIGNAL(clicked()), this, SLOT(split_surface()));
+	connect(reduction, SIGNAL(editingFinished()), this, SLOT(reduction_changed));
 
 	// setup vtk scene
 	ren3D = vtkSmartPointer<vtkRenderer>::New();
@@ -177,20 +196,24 @@ SurfaceViewerWidget::SurfaceViewerWidget(SlicesHandler* hand3D1, eInputType inpu
 	input = vtkSmartPointer<vtkImageData>::New();
 	discreteCubes = vtkSmartPointer<vtkDiscreteFlyingEdges3D>::New();
 	cubes = vtkSmartPointer<vtkFlyingEdges3D>::New();
+	decimate = vtkSmartPointer<vtkDecimatePro>::New();
+	decimate->PreserveTopologyOn();
+	decimate->BoundaryVertexDeletionOff();
+	decimate->SplittingOff();
+	decimate->SetTargetReduction(reduction->text().toDouble() / 100.0);
 
 	load();
 
 	vtkWidget->GetRenderWindow()->Render();
 }
 
-SurfaceViewerWidget::~SurfaceViewerWidget() 
+SurfaceViewerWidget::~SurfaceViewerWidget()
 {
-
 }
 
 bool SurfaceViewerWidget::isOpenGLSupported()
 {
-	// todo: check e.g. via some helper process 
+	// todo: check e.g. via some helper process
 	return true;
 }
 
@@ -278,7 +301,9 @@ void SurfaceViewerWidget::load()
 		cubes->SetInputData(input);
 		cubes->SetValue(0, range[0] + 0.01 * (range[1] - range[0]) * sl_thresh->value());
 
-		mapper->SetInputConnection(cubes->GetOutputPort());
+		decimate->SetInputConnection(cubes->GetOutputPort());
+
+		mapper->SetInputConnection(decimate->GetOutputPort());
 		mapper->ScalarVisibilityOff();
 	}
 	else
@@ -290,8 +315,9 @@ void SurfaceViewerWidget::load()
 		// merge duplicate points (check if necessary)
 		// connectivity filter & set random colors
 		// mapper set input to connectivity output
+		decimate->SetInputConnection(discreteCubes->GetOutputPort());
 
-		mapper->SetInputConnection(discreteCubes->GetOutputPort());
+		mapper->SetInputConnection(decimate->GetOutputPort());
 		if (input_type == kTarget)
 		{
 			mapper->ScalarVisibilityOff();
@@ -322,7 +348,7 @@ void SurfaceViewerWidget::split_surface()
 
 	auto num_regions = connectivity->GetNumberOfExtractedRegions();
 	ISEG_INFO("Number of disconnected regions: " << num_regions);
-	
+
 	// attach lookuptable
 	auto lut = vtkSmartPointer<vtkLookupTable>::New();
 	lut->SetNumberOfTableValues(num_regions);
@@ -336,7 +362,7 @@ void SurfaceViewerWidget::split_surface()
 
 	auto output = vtkSmartPointer<vtkPolyData>::New();
 	output->ShallowCopy(connectivity->GetOutput());
-	
+
 	mapper->SetInputData(output);
 	mapper->ScalarVisibilityOn();
 	mapper->SetColorModeToMapScalars();
@@ -555,6 +581,13 @@ void SurfaceViewerWidget::thresh_changed()
 	}
 }
 
+void SurfaceViewerWidget::reduction_changed()
+{
+	decimate->SetTargetReduction(reduction->text().toDouble() / 100.0);
+
+	vtkWidget->GetRenderWindow()->Render();
+}
+
 int SurfaceViewerWidget::get_picked_tissue() const
 {
 	double* worldPosition = picker->GetPickPosition();
@@ -580,4 +613,4 @@ int SurfaceViewerWidget::get_picked_tissue() const
 	return -1;
 }
 
-}// namespace iseg
+} // namespace iseg

iSeg/SurfaceViewerWidget.h

@@ -24,6 +24,7 @@ class QVTKWidget;
 class QVTKInteractor;
 class QSlider;
 class QLabel;
+class QLineEdit;
 class QPushButton;
 class QCheckBox;
 
@@ -32,6 +33,7 @@ class vtkInteractorStyleTrackballCamera;
 class vtkImageData;
 class vtkFlyingEdges3D;
 class vtkDiscreteFlyingEdges3D;
+class vtkDecimatePro;
 class vtkPolyDataMapper;
 class vtkRenderer;
 class vtkEventQtSlotConnect;
@@ -75,6 +77,7 @@ public slots:
 protected slots:
 	void transp_changed();
 	void thresh_changed();
+	void reduction_changed();
 	void popup(vtkObject* obj, unsigned long, void* client_data, void*, vtkCommand* command);
 	void select_action(QAction*);
 
@@ -90,6 +93,7 @@ protected slots:
 	QVTKWidget* vtkWidget;
 	QSlider* sl_trans;
 	QSlider* sl_thresh;
+	QLineEdit* reduction;
 	QPushButton* bt_update;
 	QPushButton* bt_connectivity;
 	QLabel* lb_connectivity_count;
@@ -102,6 +106,7 @@ protected slots:
 	vtkSmartPointer<vtkInteractorStyleTrackballCamera> style;
 	vtkSmartPointer<vtkDiscreteFlyingEdges3D> discreteCubes;
 	vtkSmartPointer<vtkFlyingEdges3D> cubes;
+	vtkSmartPointer<vtkDecimatePro> decimate;
 	vtkSmartPointer<vtkPolyDataMapper> mapper;
 	vtkSmartPointer<vtkActor> actor;
 	vtkSmartPointer<vtkLookupTable> lut;