use python to randomly inset each face, bmesh?
You can't loop through bm.faces, because you are creating more faces in the loop, which are also iterated over. Instead make a copy of the faces

faces_copy = [f for f in bm.faces]


Now, your code would give us the following result.



Notable problems
- Insets, even at 0 depth.
- Your example image shows and inset with 0 depth and then a random extrusion.

Let's implement that. Luckily, after a inset_region operation, the face index stays the same so we can simply extrude it by insetting another time and moving the face's vertices along the face normal.

bmesh.ops.inset_region(bm, faces = [face], thickness=0.4, depth=0)
bmesh.ops.inset_region(bm, faces = [face], thickness=0, depth=0)
bmesh.ops.translate(bm, verts = face.verts, vec = face.normal)




Now lets add the randomness:

import random
#[...]

bmesh.ops.inset_region(bm, faces = [face], thickness=random.uniform(0.1, 0.9), depth=0)
bmesh.ops.inset_region(bm, faces = [face], thickness=0, depth=0)
bmesh.ops.translate(bm, verts = face.verts, vec = face.normal * random.uniform(0.1, 4))


Iterations

We need to catch the newly "extruded" faces. Luckily inset_region returns geometry. We can create a new face list like this:

geom = bmesh.ops.inset_region(bm, faces = [face], thickness=0, depth=0)
new_faces = geom['faces'] # list of bmesh faces


Create an empty list for the new faces. After the insetting and extrusion, we'll add the active face and the extruded side faces.

new_faces.extend(geom['faces'])
new_faces.append(face)


Complete code (with two iterations):

import bpy
import bmesh
import random
from mathutils import Vector

obj = bpy.context.active_object

bm = bmesh.new()
bm.from_mesh(obj.data,face_normals=True)

# faces, which could be extruded
faces_copy = [f for f in bm.faces]
new_faces = []
# number of iterations goes here
for i in range(0, 2):
    for face in faces_copy:
        do_inset = random.randint(0,1)
        if do_inset:
            bmesh.ops.inset_region(bm, faces = [face], thickness=random.uniform(0.1, 0.9), depth=0)
            geom = bmesh.ops.inset_region(bm, faces = [face], thickness=0, depth=0)
            bmesh.ops.translate(bm, verts = face.verts, vec = face.normal * random.uniform(0.1, 4))
            new_faces.extend(geom['faces'])
            new_faces.append(face)
    faces_copy = new_faces
    new_faces = []

bm.to_mesh(obj.data)
bm.free()
obj.data.update()




The concept is pretty cool, but it already show, that we need to implement a method to take the faces size into account before using bmesh.ops.translate.