tests

Thiago Franco de Moraes
1 parent b8a2a2a8
Showing 3 changed files with 50 additions and 20 deletions Show diff stats
invesalius/data/smooth_cy.pyx
invesalius/data/surface.py
invesalius/data/surface_process.py
@@ -118,7 +118,7 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
     cdef np.ndarray[DTYPE8_t, ndim=3] _mask = np.zeros_like(image)
     cdef np.ndarray[DTYPEF64_t, ndim=3] aux = np.zeros_like(out)
  
-    cdef int i, x, y, z, S;
+    cdef int i, x, y, z, S
     cdef DTYPEF64_t H, v, cn
     cdef DTYPEF64_t diff=0.0
     cdef DTYPEF64_t dt=1/6.0
@@ -152,31 +152,31 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
                     S += 1
  
     for i in xrange(n):
-        replicate(out, aux);
-        diff = 0.0;
+        replicate(out, aux)
+        diff = 0.0
  
         for z in xrange(dz):
             for y in xrange(dy):
                 for x in xrange(dx):
                     if mask[z, y, x]:
-                        H = calculate_H(aux, z, y, x);
-                        v = aux[z, y, x] + dt*H;
+                        H = calculate_H(aux, z, y, x)
+                        v = aux[z, y, x] + dt*H
  
                         if image[z, y, x]:
-                            if v > 0:
-                                out[z, y, x] = v
-                            else:
-                                out[z, y, x] = 0.0001
-                        else:
                             if v < 0:
+                                out[z, y, x] = 0.0
+                            else:
                                 out[z, y, x] = v
+                        else:
+                            if v > 0:
+                                out[z, y, x] = 0.0
                             else:
-                                out[z, y, x] = -0.0001
+                                out[z, y, x] = v
  
                         diff += (out[z, y, x] - aux[z, y, x])*(out[z, y, x] - aux[z, y, x])
  
-        cn = sqrt((1.0/S) * diff);
+        cn = sqrt((1.0/S) * diff)
         print "%d - CN: %.28f - diff: %.28f\n" % (i, cn, diff)
  
         if cn <= E:
-            break;
+            break
@@ -472,11 +472,9 @@ class SurfaceManager():
             iteractions = 10
             bsize = 4
             mask_tfile = tempfile.mktemp()
+
             mmask = np.memmap(mask_tfile, shape=mask.matrix.shape, dtype='float64', mode='w+')
-            smooth_cy.smooth(mask.matrix, iteractions, bsize, mmask)
-            mmask[0, :, :] = -1
-            mmask[:, 0, :] = -1
-            mmask[:, :, 0] = -1
+            smooth_cy.smooth((mask.matrix > 127).astype('uint8'), iteractions, bsize, mmask)
  
         else:
             mask_tfile = mask.temp_file
@@ -592,6 +590,38 @@ class SurfaceManager():
             #  polydata.SetSource(None)
             #  polydata.DebugOn()
  
+        elif algorithm == 'Whitaker':
+            normals = vtk.vtkPolyDataNormals()
+            normals_ref = weakref.ref(normals)
+            normals_ref().AddObserver("ProgressEvent", lambda obj,evt:
+                                      UpdateProgress(normals_ref(), _("Creating 3D surface...")))
+            normals.SetInput(polydata)
+            normals.ReleaseDataFlagOn()
+            #normals.SetFeatureAngle(80)
+            #normals.AutoOrientNormalsOn()
+            normals.ComputeCellNormalsOn()
+            normals.GetOutput().ReleaseDataFlagOn()
+            normals.Update()
+            del polydata
+            polydata = normals.GetOutput()
+            polydata.SetSource(None)
+            del normals
+
+            clean = vtk.vtkCleanPolyData()
+            clean.ReleaseDataFlagOn()
+            clean.GetOutput().ReleaseDataFlagOn()
+            clean_ref = weakref.ref(clean)
+            clean_ref().AddObserver("ProgressEvent", lambda obj,evt:
+                            UpdateProgress(clean_ref(), _("Creating 3D surface...")))
+            clean.SetInput(polydata)
+            clean.PointMergingOn()
+            clean.Update()
+
+            del polydata
+            polydata = clean.GetOutput()
+            polydata.SetSource(None)
+            del clean
+
         else:
             #smoother = vtk.vtkWindowedSincPolyDataFilter()
             smoother = vtk.vtkSmoothPolyDataFilter()
@@ -158,9 +158,9 @@ class SurfaceProcess(multiprocessing.Process):
         else:
             contour.SetValue(0, self.min_value) # initial threshold
             contour.SetValue(1, self.max_value) # final threshold
-        contour.ComputeScalarsOn()
-        contour.ComputeGradientsOn()
-        contour.ComputeNormalsOn()
+        #  contour.ComputeScalarsOn()
+        #  contour.ComputeGradientsOn()
+        #  contour.ComputeNormalsOn()
         contour.ReleaseDataFlagOn()
         contour.Update()
         #contour.AddObserver("ProgressEvent", lambda obj,evt:
...	...	@@ -118,7 +118,7 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
118	118	cdef np.ndarray[DTYPE8_t, ndim=3] _mask = np.zeros_like(image)
119	119	cdef np.ndarray[DTYPEF64_t, ndim=3] aux = np.zeros_like(out)
120	120
121		- cdef int i, x, y, z, S;
	121	+ cdef int i, x, y, z, S
122	122	cdef DTYPEF64_t H, v, cn
123	123	cdef DTYPEF64_t diff=0.0
124	124	cdef DTYPEF64_t dt=1/6.0
...	...	@@ -152,31 +152,31 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
152	152	S += 1
153	153
154	154	for i in xrange(n):
155		- replicate(out, aux);
156		- diff = 0.0;
	155	+ replicate(out, aux)
	156	+ diff = 0.0
157	157
158	158	for z in xrange(dz):
159	159	for y in xrange(dy):
160	160	for x in xrange(dx):
161	161	if mask[z, y, x]:
162		- H = calculate_H(aux, z, y, x);
163		- v = aux[z, y, x] + dt*H;
	162	+ H = calculate_H(aux, z, y, x)
	163	+ v = aux[z, y, x] + dt*H
164	164
165	165	if image[z, y, x]:
166		- if v > 0:
167		- out[z, y, x] = v
168		- else:
169		- out[z, y, x] = 0.0001
170		- else:
171	166	if v < 0:
	167	+ out[z, y, x] = 0.0
	168	+ else:
172	169	out[z, y, x] = v
	170	+ else:
	171	+ if v > 0:
	172	+ out[z, y, x] = 0.0
173	173	else:
174		- out[z, y, x] = -0.0001
	174	+ out[z, y, x] = v
175	175
176	176	diff += (out[z, y, x] - aux[z, y, x])*(out[z, y, x] - aux[z, y, x])
177	177
178		- cn = sqrt((1.0/S) * diff);
	178	+ cn = sqrt((1.0/S) * diff)
179	179	print "%d - CN: %.28f - diff: %.28f\n" % (i, cn, diff)
180	180
181	181	if cn <= E:
182		- break;
	182	+ break
...	...
...	...	@@ -472,11 +472,9 @@ class SurfaceManager():
472	472	iteractions = 10
473	473	bsize = 4
474	474	mask_tfile = tempfile.mktemp()
	475	+
475	476	mmask = np.memmap(mask_tfile, shape=mask.matrix.shape, dtype='float64', mode='w+')
476		- smooth_cy.smooth(mask.matrix, iteractions, bsize, mmask)
477		- mmask[0, :, :] = -1
478		- mmask[:, 0, :] = -1
479		- mmask[:, :, 0] = -1
	477	+ smooth_cy.smooth((mask.matrix > 127).astype('uint8'), iteractions, bsize, mmask)
480	478
481	479	else:
482	480	mask_tfile = mask.temp_file
...	...	@@ -592,6 +590,38 @@ class SurfaceManager():
592	590	# polydata.SetSource(None)
593	591	# polydata.DebugOn()
594	592
	593	+ elif algorithm == 'Whitaker':
	594	+ normals = vtk.vtkPolyDataNormals()
	595	+ normals_ref = weakref.ref(normals)
	596	+ normals_ref().AddObserver("ProgressEvent", lambda obj,evt:
	597	+ UpdateProgress(normals_ref(), _("Creating 3D surface...")))
	598	+ normals.SetInput(polydata)
	599	+ normals.ReleaseDataFlagOn()
	600	+ #normals.SetFeatureAngle(80)
	601	+ #normals.AutoOrientNormalsOn()
	602	+ normals.ComputeCellNormalsOn()
	603	+ normals.GetOutput().ReleaseDataFlagOn()
	604	+ normals.Update()
	605	+ del polydata
	606	+ polydata = normals.GetOutput()
	607	+ polydata.SetSource(None)
	608	+ del normals
	609	+
	610	+ clean = vtk.vtkCleanPolyData()
	611	+ clean.ReleaseDataFlagOn()
	612	+ clean.GetOutput().ReleaseDataFlagOn()
	613	+ clean_ref = weakref.ref(clean)
	614	+ clean_ref().AddObserver("ProgressEvent", lambda obj,evt:
	615	+ UpdateProgress(clean_ref(), _("Creating 3D surface...")))
	616	+ clean.SetInput(polydata)
	617	+ clean.PointMergingOn()
	618	+ clean.Update()
	619	+
	620	+ del polydata
	621	+ polydata = clean.GetOutput()
	622	+ polydata.SetSource(None)
	623	+ del clean
	624	+
595	625	else:
596	626	#smoother = vtk.vtkWindowedSincPolyDataFilter()
597	627	smoother = vtk.vtkSmoothPolyDataFilter()
...	...
...	...	@@ -158,9 +158,9 @@ class SurfaceProcess(multiprocessing.Process):
158	158	else:
159	159	contour.SetValue(0, self.min_value) # initial threshold
160	160	contour.SetValue(1, self.max_value) # final threshold
161		- contour.ComputeScalarsOn()
162		- contour.ComputeGradientsOn()
163		- contour.ComputeNormalsOn()
	161	+ # contour.ComputeScalarsOn()
	162	+ # contour.ComputeGradientsOn()
	163	+ # contour.ComputeNormalsOn()
164	164	contour.ReleaseDataFlagOn()
165	165	contour.Update()
166	166	#contour.AddObserver("ProgressEvent", lambda obj,evt:
...	...