def METHOD_NAME():
    options, flags = gs.parser()
    # TODO: allow using existing slices with a flag
    input_ = options["input"]
    output_basename = options["output"]
    slices_basename = options["slices"]
    surface = options["surface"]
    divide = flags["d"]
    use_slices = flags["s"]
    check_raster3d_exists(input_)
    if surface:
        check_raster2d_exists(surface)
    size = None
    if options["size"]:
        size = int(options["size"])
        if size % 2 == 0:
            gs.fatal(
                _("Please provide an odd number for the moving" " window size, not %d")
                % size
            )
    # options to be passed to r3.to.rast
    additional_options = {}
    for option in ["multiply", "add"]:
        if options[option]:
            additional_options[option] = options[option]
    additional_options["type"] = "CELL"
    basename_sep = "_"
    output_pattern = "^%s%s[0-9]+$" % (output_basename, basename_sep)
    slices_pattern = "^%s%s[0-9]+$" % (slices_basename, basename_sep)
    # TODO: slices test seems to be broken
    overwrite = gs.overwrite()
    rasters = list_in_current_mapset(output_pattern)
    if not overwrite and rasters:
        basename_in_use_message(output_basename, rasters)
    rasters = list_in_current_mapset(slices_pattern)
    if not use_slices and not overwrite and rasters:
        basename_in_use_message(slices_basename, rasters)
    # TODO: list of categories would be better
    info = gs.parse_command("r3.info", map=input_, flags="r")
    map_min = int(float(info["min"]))
    map_max = int(float(info["max"]))
    gs.verbose(_("Categories from %d to %d") % (map_min, map_max))
    rasters = []
    try:
        # TODO: we should switch to 3D region for the following computations?
        # or should we require resolutions to be the same?
        if not use_slices:
            gs.run_command(
                "r3.to.rast", input=input_, output=slices_basename, **additional_options
            )
        # it's clear we created just the ones in the current mapset
        # and we need to make the command line as short as possible
        rasters = list_in_current_mapset(slices_pattern)
        if size:
            base_expr = "int({m}[{a},{b}] == $num)"
            size = int((size - 1) / 2)  # window size to max index
            expr_list = []
            for raster in rasters:
                for j in range(-size, size + 1):
                    for i in range(-size, size + 1):
                        expr_list.append(base_expr.format(m=raster, a=i, b=j))
            expr = " + ".join(expr_list)
            expr = "{out}{sep}$num = {inp}".format(
                out=output_basename, sep=basename_sep, inp=expr
            )
        else:
            if surface:
                # TODO: < or <= ?
                base_expr = "if({d} < {s}, int({m} == $num), 0)"
                expr_list = []
                for depth, raster in enumerate(rasters):
                    expr_list.append(base_expr.format(m=raster, d=depth + 1, s=surface))
                expr = " + ".join(expr_list)
                if divide:
                    expr = "({e}) / {otype}({s})".format(
                        e=expr, otype="float", s=surface
                    )
                expr = "{out}{sep}$num = {inp}".format(
                    out=output_basename, sep=basename_sep, inp=expr
                )
            else:
                expr = "{out}{sep}$num = int({inp} == $num)".format(
                    out=output_basename,
                    sep=basename_sep,
                    inp=" == $num) + int(".join(rasters),
                )
        # TODO: define behavior when the map is empty or has just one class
        try:
            # TODO: this uses uncommitted experimental code for trunk
            from grass.script.parallel import ModuleCallList, execute_by_module
            call = ModuleCallList()
            parallel = True
        except ImportError:
            call = gs  # fall back to grass.script
            parallel = False
        for i in range(map_min, map_max + 1):
            call.mapcalc(expr, num=i)
        if parallel:
            execute_by_module(call, nprocs=4)
    except CalledModuleError as error:
        remove(rasters)
        gs.fatal(_("Module %s failed. Check the above error messages.") % error.cmd)