Particle Systems

import bpy

def quantize(original_obj, numdiv_shortest_side):

    # make a plane
    bpy.ops.mesh.primitive_plane_add(radius=0.5)
    plane = bpy.context.active_object

    # add particle system to original object
    ps = original_obj.modifiers.new("grid particles", type='PARTICLE_SYSTEM')
    psettings = ps.particle_system.settings

    # particle settings
    psettings.distribution = 'GRID'
    psettings.emit_from = 'FACE'
    psettings.physics_type = 'NO'
    psettings.grid_resolution = numdiv_shortest_side
    psettings.use_render_emitter = True
    psettings.show_unborn = True
    psettings.use_scale_dupli = True
    psettings.particle_size = 1.0
    psettings.render_type = 'OBJECT'
    psettings.dupli_object = plane

    plane.select = False
    
    original_obj.select = True
    bpy.context.scene.objects.active = original_obj
    bpy.ops.object.duplicates_make_real()
    original_obj.hide = True
    original_obj.select = False
    
    bpy.ops.object.make_local(type='SELECT_OBDATA')

    # find first of the new objects, to join the rest onto
    joiner = None
    for obj in bpy.data.objects:
        if not obj in {plane, original_obj}:
            joiner = obj
            bpy.context.scene.objects.active = obj
            break
    
    bpy.ops.object.delete(use_global=False)

    
    # join, and remove doubles
    bpy.ops.object.join()
    # bpy.ops.mesh.remove_doubles()


numdiv_shortest_side = 21
original_obj = bpy.context.active_object
quantize(original_obj, numdiv_shortest_side)