Find and merge groups of polygons in the surface that meet the 
        following criteria:
        
            * Are coplanars.
            * Are contiguous.
            * The result is convex.
            
        This method is very useful at reducing the number the items and,
        therefore, the shadowing time computing. Before override this
        instance, it is saved and can be restored with ``.restore()``
        
        :param plot: If True, generates the before and after 
            visualizations for the surface. Use it to check the results.
        :type plot: bool
        :returns: None
        
        .. warning:: This method do not check if the merged polygons are 
            actually convex. The convex hull of the union is directly 
            calculated. For this reason, it is very important to visualy
            check the solution.

def melt(self, plot=False):
        """
        Find and merge groups of polygons in the surface that meet the 
        following criteria:
        
            * Are coplanars.
            * Are contiguous.
            * The result is convex.
            
        This method is very useful at reducing the number the items and,
        therefore, the shadowing time computing. Before override this
        instance, it is saved and can be restored with ``.restore()``
        
        :param plot: If True, generates the before and after 
            visualizations for the surface. Use it to check the results.
        :type plot: bool
        :returns: None
        
        .. warning:: This method do not check if the merged polygons are 
            actually convex. The convex hull of the union is directly 
            calculated. For this reason, it is very important to visualy
            check the solution.
        """
        from pyny3d.utils import bool2index
        from scipy.spatial import ConvexHull
        
        # First, coplanarity
        ## Normalize parametric equations
        para = [poly.get_parametric() for poly in self]
        para = np.array([p/np.linalg.norm(p) for p in para])
        n = para.shape[0]
        ## Coincidences
        cop = []
        for i, plane in enumerate(para[:-1]):
            indexes = np.zeros((n-i-1, 4))
            for c in range(4):
                indexes[:, c] = np.isclose(para[i+1:, c], plane[c])
            pos = bool2index(indexes.sum(axis=1)==4)+i+1
            if pos.shape[0] > 0:
                cop.append(np.hstack((i, pos)))
                para[pos, :] = np.nan
        
        # Second, contiguity
        substituted = []
        cop_cont = []
        for i, group in enumerate(cop):
            polygons = [self[i] for i in group]
            if Surface.contiguous(polygons):
                cop_cont.append(polygons)
                substituted.append(group)
                        
        if len(substituted) != 0:
            self.save()
            if plot: self.plot()
            substituted = sum(substituted)
                    
            # Hull        
            merged = []
            for polygons in cop_cont:
                points = np.concatenate([polygon.points 
                                         for polygon in polygons])
                hull = ConvexHull(points[:, :2])
                merged.append(Polygon(points[hull.vertices]))
    
            # Final substitution
            new_surface = [self[i] for i in range(len(self.polygons))
                                   if i not in substituted]
            new_surface += merged
            self.polygons = new_surface
            self.sorted_areas = None
            
            if plot: self.plot()