diff --git a/docs/changes.md b/docs/changes.md
index 84957c05..720684fd 100644
--- a/docs/changes.md
+++ b/docs/changes.md
@@ -17,6 +17,7 @@
 - add `smooth_data()` to smooth/diffuse data attributes
 - add `shapes.Tubes()`
 - add `utils.Minimizer()` class
+- add `CellCenters(Points)` class
 
 
 ## Breaking changes
@@ -25,6 +26,7 @@
 - improvements to `pointcloud.pca_ellipse()` and bug fixes
 - change `clone2d(scale=...)` to `clone2d(size=...)`
 - remove `shapes.StreamLines()` becoming `object.compute_streamlines()`
+- split `mesh.decimate()` into `mesh.decimate()`, `mesh.decimate_pro()` and `mesh.decimate_binned()` as per #992
 
 
 ### Renaming
diff --git a/vedo/core.py b/vedo/core.py
index 0b5416a8..ccc27c60 100644
--- a/vedo/core.py
+++ b/vedo/core.py
@@ -622,8 +622,11 @@ def cell_centers(self):
 
         Examples:
             - [delaunay2d.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delaunay2d.py)
+        
+        See also: `CellCenters()`.
         """
         vcen = vtk.new("CellCenters")
+        vcen.CopyArraysOff()
         vcen.SetInputData(self.dataset)
         vcen.Update()
         return utils.vtk2numpy(vcen.GetOutput().GetPoints().GetData())
@@ -941,7 +944,7 @@ def resample_data_from(self, source, tol=None, categorical=False):
         m1.pointdata["xvalues"] = np.power(pts[:,0], 3)
         m1.celldata["yvalues"]  = np.power(ces[:,1], 3)
         m2 = Mesh(dataurl+'bunny.obj')
-        m2.resample_arrays_from(m1)
+        m2.resample_data_from(m1)
         # print(m2.pointdata["xvalues"])
         show(m1, m2 , N=2, axes=1)
         ```
@@ -1053,7 +1056,7 @@ def interpolate_data_from(
         interpolator.SetSourceData(points)
         interpolator.SetKernel(kern)
         interpolator.SetLocator(locator)
-        interpolator.PassFieldArraysOff()
+        interpolator.PassFieldArraysOn()
         interpolator.SetNullPointsStrategy(null_strategy)
         interpolator.SetNullValue(null_value)
         interpolator.SetValidPointsMaskArrayName("ValidPointMask")
diff --git a/vedo/mesh.py b/vedo/mesh.py
index fa19b2de..a06040c0 100644
--- a/vedo/mesh.py
+++ b/vedo/mesh.py
@@ -959,24 +959,35 @@ def subdivide(self, n=1, method=0, mel=None):
         )
         return self
 
-    def decimate(self, fraction=0.5, n=None, method="quadric", boundaries=False):
+
+    def decimate(self, fraction=0.5, n=None, preserve_volume=True, regularize=True):
         """
         Downsample the number of vertices in a mesh to `fraction`.
 
+        This filter preserves the `pointdata` of the input dataset.
+
         Arguments:
             fraction : (float)
                 the desired target of reduction.
             n : (int)
-                the desired number of final points (`fraction` is recalculated based on it).
-            method : (str)
-                can be either 'quadric' or 'pro'. In the first case triagulation
-                will look like more regular, irrespective of the mesh original curvature.
-                In the second case triangles are more irregular but mesh is more precise on more
-                curved regions.
-            boundaries : (bool)
-                in "pro" mode decide whether to leave boundaries untouched or not
+                the desired number of final points
+                (`fraction` is recalculated based on it).
+            preserve_volume : (bool)
+                Decide whether to activate volume preservation which greatly
+                reduces errors in triangle normal direction.
+            regularize : (bool)
+                regularize the point finding algorithm so as to have better quality
+                mesh elements at the cost of a slightly lower precision on the
+                geometry potentially (mostly at sharp edges).
+                Can be useful for decimating meshes that have been triangulated on noisy data.
 
-        .. note:: Setting `fraction=0.1` leaves 10% of the original number of vertices
+        Note:
+            Setting `fraction=0.1` leaves 10% of the original number of vertices.
+            Internally the VTK class
+            [vtkQuadricDecimation](https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html)
+            is used for this operation.
+        
+        See also: `decimate_binned()` and `decimate_pro()`.
         """
         poly = self.dataset
         if n:  # N = desired number of points
@@ -985,22 +996,17 @@ def decimate(self, fraction=0.5, n=None, method="quadric", boundaries=False):
             if fraction >= 1:
                 return self
 
-        if "quad" in method:
-            decimate = vtk.new("QuadricDecimation")
-            # decimate.SetVolumePreservation(True)
+        decimate = vtk.new("QuadricDecimation")
+        decimate.SetVolumePreservation(preserve_volume)
+        decimate.SetRegularize(regularize)
+        decimate.MapPointDataOn()
 
-        else:
-            decimate = vtk.new("DecimatePro")
-            decimate.PreserveTopologyOn()
-            if boundaries:
-                decimate.BoundaryVertexDeletionOff()
-            else:
-                decimate.BoundaryVertexDeletionOn()
-        decimate.SetInputData(poly)
         decimate.SetTargetReduction(1 - fraction)
+        decimate.SetInputData(poly)
         decimate.Update()
 
         self._update(decimate.GetOutput())
+        self.metadata["decimate_actual_fraction"] = 1 - decimate.GetActualReduction()
 
         self.pipeline = OperationNode(
             "decimate",
@@ -1008,6 +1014,133 @@ def decimate(self, fraction=0.5, n=None, method="quadric", boundaries=False):
             comment=f"#pts {self.dataset.GetNumberOfPoints()}",
         )
         return self
+    
+    def decimate_pro(
+            self,
+            fraction=0.5,
+            n=None,
+            preserve_topology=True,
+            preserve_boundaries=True,
+            splitting=False,
+            splitting_angle=75,
+            feature_angle=0,
+            inflection_point_ratio=10,
+        ):
+        """
+        Downsample the number of vertices in a mesh to `fraction`.
+
+        This filter preserves the `pointdata` of the input dataset.
+
+        Arguments:
+            fraction : (float)
+                The desired target of reduction.
+                Setting `fraction=0.1` leaves 10% of the original number of vertices.
+            n : (int)
+                the desired number of final points (`fraction` is recalculated based on it).
+            preserve_topology : (bool)
+                If on, mesh splitting and hole elimination will not occur.
+                This may limit the maximum reduction that may be achieved.
+            preserve_boundaries : (bool)
+                Turn on/off the deletion of vertices on the boundary of a mesh.
+                Control whether mesh boundaries are preserved during decimation.
+            feature_angle : (float)
+                Specify the angle that defines a feature.
+                This angle is used to define what an edge is
+                (i.e., if the surface normal between two adjacent triangles
+                is >= FeatureAngle, an edge exists).
+            splitting : (bool)
+                Turn on/off the splitting of the mesh at corners,
+                along edges, at non-manifold points, or anywhere else a split is required.
+                Turning splitting off will better preserve the original topology of the mesh,
+                but you may not obtain the requested reduction.
+            splitting_angle : (float)
+                Specify the angle that defines a sharp edge.
+                This angle is used to control the splitting of the mesh.
+                A split line exists when the surface normals between two edge connected triangles
+                are >= `splitting_angle`.
+            inflection_point_ratio : (float)
+                An inflection point occurs when the ratio of reduction error between two iterations
+                is greater than or equal to the `inflection_point_ratio` value.
+
+        Note:
+            Setting `fraction=0.1` leaves 10% of the original number of vertices
+        
+        See also:
+            `decimate()` and `decimate_binned()`.
+        """
+        poly = self.dataset
+        if n:  # N = desired number of points
+            npt = poly.GetNumberOfPoints()
+            fraction = n / npt
+            if fraction >= 1:
+                return self
+
+        decimate = vtk.new("DecimatePro")
+        decimate.SetPreserveTopology(preserve_topology)
+        decimate.SetBoundaryVertexDeletion(preserve_boundaries)
+        if feature_angle:
+            decimate.SetFeatureAngle(feature_angle)
+        decimate.SetSplitting(splitting)
+        decimate.SetSplitAngle(splitting_angle)
+        decimate.SetInflectionPointRatio(inflection_point_ratio)
+
+        decimate.SetTargetReduction(1 - fraction)
+        decimate.SetInputData(poly)
+        decimate.Update()
+        self._update(decimate.GetOutput())
+
+        self.pipeline = OperationNode(
+            "decimate_pro",
+            parents=[self],
+            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
+        )
+        return self
+    
+    def decimate_binned(self, divisions=(), use_clustering=False):
+        """
+        Downsample the number of vertices in a mesh.
+        
+        This filter preserves the `celldata` of the input dataset,
+        if `use_clustering=True` also the `pointdata` will be preserved in the result.
+
+        Arguments:
+            divisions : (list)
+                number of divisions along x, y and z axes.
+            auto_adjust : (bool)
+                if True, the number of divisions is automatically adjusted to
+                create more uniform cells.
+            use_clustering : (bool)
+                use [vtkQuadricClustering](https://vtk.org/doc/nightly/html/classvtkQuadricClustering.html)
+                instead of 
+                [vtkBinnedDecimation](https://vtk.org/doc/nightly/html/classvtkBinnedDecimation.html).
+        
+        See also: `decimate()` and `decimate_pro()`.
+        """
+        if use_clustering:
+            decimate = vtk.new("QuadricClustering")
+            decimate.CopyCellDataOn()
+        else:
+            decimate = vtk.new("BinnedDecimation")
+            decimate.ProducePointDataOn()
+            decimate.ProduceCellDataOn()
+
+        decimate.SetInputData(self.dataset)
+
+        if len(divisions) == 0:
+            decimate.SetAutoAdjustNumberOfDivisions(1)
+        else:
+            decimate.SetAutoAdjustNumberOfDivisions(0)
+            decimate.SetNumberOfDivisions(divisions)
+        decimate.Update()
+
+        self._update(decimate.GetOutput())
+        self.metadata["decimate_binned_divisions"] = decimate.GetNumberOfDivisions()
+        self.pipeline = OperationNode(
+            "decimate_binned",
+            parents=[self],
+            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
+        )
+        return self
 
     def delete_cells(self, ids):
         """
diff --git a/vedo/pointcloud.py b/vedo/pointcloud.py
index b7fbd3ce..397e6389 100644
--- a/vedo/pointcloud.py
+++ b/vedo/pointcloud.py
@@ -24,6 +24,7 @@
 __all__ = [
     "Points",
     "Point",
+    "CellCenters",
     "merge",
     "delaunay2d",  # deprecated, use .generate_delaunay2d()
     "fit_line",
@@ -3579,6 +3580,26 @@ def visible_points(self, area=(), tol=None, invert=False):
             svp.SelectInvisibleOn()
         svp.Update()
 
-        m = Points(svp.GetOutput())#.point_size(5)
+        m = Points(svp.GetOutput())
         m.name = "VisiblePoints"
         return m
+
+####################################################
+class CellCenters(Points):
+    def __init__(self, pcloud):
+        """
+        Generate `Points` at the center of the cells of any type of object.
+
+        Check out also `cell_centers()`.
+        """
+        vcen = vtk.new("CellCenters")
+        vcen.CopyArraysOn()
+        vcen.VertexCellsOn()
+        # vcen.ConvertGhostCellsToGhostPointsOn()
+        try:
+            vcen.SetInputData(pcloud.dataset)
+        except AttributeError:
+            vcen.SetInputData(pcloud)
+        vcen.Update()
+        super().__init__(vcen.GetOutput())
+        self.name = "CellCenters"
diff --git a/vedo/version.py b/vedo/version.py
index fcf17404..5fbefcbb 100644
--- a/vedo/version.py
+++ b/vedo/version.py
@@ -1 +1 @@
-_version = '2023.5.0+dev16'
+_version = '2023.5.0+dev17'
diff --git a/vedo/vtkclasses.py b/vedo/vtkclasses.py
index f5a5bbce..dc76edf8 100644
--- a/vedo/vtkclasses.py
+++ b/vedo/vtkclasses.py
@@ -246,6 +246,7 @@
     "vtk3DLinearGridCrinkleExtractor",
     "vtkAppendFilter",
     "vtkAppendPolyData",
+    "vtkBinnedDecimation",
     "vtkCellCenters",
     "vtkCellDataToPointData",
     "vtkCenterOfMass",
@@ -272,6 +273,7 @@
     "vtkPointDataToCellData",
     "vtkPolyDataNormals",
     "vtkProbeFilter",
+    "vtkQuadricClustering",
     "vtkQuadricDecimation",
     "vtkResampleWithDataSet",
     "vtkReverseSense",