How can I unsubdivide regions of flat faces with Python?
Edit Mode bmesh script.


Tested on subdivided and poked cube.


Tags all edges that have a face angle > EPS
Gets all triangular faces made up of only edges not tagged and joins them. (To Quad)
Feeds all verts that have no linked edges tagged into bmesh unsubdivide operator. 


Can edit properties of both to quads, eg max face angle, and iterations for unsubdivide

Test script, run in edit mode. Needs no selection.

import bmesh
import bpy
from math import radians

EPS = 1e-5 # edge angle
ITER = 10 # unsubdiv iters
AST = radians(40)
AFT = radians(40)        

me = bpy.context.object.data
bm = bmesh.from_edit_mesh(me)

for e in bm.edges:
    if len(e.link_faces) &lt; 2:
        e.tag = True
    else:
        e.tag = e.calc_face_angle() &gt; EPS


ret = bmesh.ops.join_triangles(
        bm,
        faces=[f for f in bm.faces if all(not e.tag for e in f.edges)],
        angle_face_threshold=AFT, 
        angle_shape_threshold=AST,
        )

bmesh.ops.unsubdivide(bm,
        verts=[v for v in bm.verts if not any(e.tag for e in v.link_edges)],
        iterations=ITER,
        )

bmesh.update_edit_mesh(me) 


Alternatively Dissolve edges with 0 face angle.

In edge select mode select an edge with no face angle, then select similar by face angle



Hit X

and dissolve edges.

A script



A quick edit mode bmesh script to do this.

import bmesh
import bpy

EPS = 1e-5


me = bpy.context.object.data
bm = bmesh.from_edit_mesh(me)


bmesh.ops.dissolve_edges(bm, 
        edges = [e for e in bm.edges if e.calc_face_angle() &lt; EPS],
        )

bmesh.update_edit_mesh(me)        


Not sure how it will go on 5M ++ faces.  Could change slightly to work on selection.

        edges = [e for e in bm.edges if e.calc_face_angle() &lt; EPS and e.select]


This will ofcourse create a number of ngons, that can be grid filled, and fed into the unsubdivided.