def grain(f, labels, measurement, option='image'):
    """
    y = grain(f, labels, measurement, option="image")

    Gray-scale statistics for each labeled region.

    Computes gray-scale statistics of each grain in the image. The
    grains regions are specified by the labeled image `labels` and the
    gray-scale information is specified by the image `f`. The
    statistics to compute is specified by the parameter measurement,
    which has the same options as in function stats. The
    parameter option defines the output format:

    'image'
        if the output is an uint16 image where each label value is changed to
        the measurement value
    'data'
        a double column vector. In this case, the first element (index 0) is
        the measurement of region 1. The region with label zero is not measure
        as it is normally the background.

    Parameters
    ----------
    f : ndarray (gray-scale)
    labels : ndarray of integer dtype
        Labeled image. As usual, label 0 is backaground.
    measurement : {'max', 'min', 'median', 'mean', 'sum', 'std', 'std1' }
        Which measure to compute.
    option : {'image' [default], 'data' }
        Output format: if 'image', results as a gray-scale mosaic image
        (default): if 'data', results a column vector of measurements.

    Returns
    -------
    y: Gray-scale (uint8 or uint16) image.
        Or an array with gray-scale statistics per region
    """
    """
        - Examples
            #
            #   example 1
            1#
            f=to_uint8([range(6),range(6),range(6)])
            labels=labelflat(f)
            grain(labels,f,'sum','data')
            grain(labels,f,'sum')
            #
            #   example 2
            #
            f=readgray('astablet.tif')
            g=gradm(f)
            marker=regmin(close(g))
            ws=cwatershed(g,marker,sebox(),'regions')
            g=grain(ws,f,'mean')
            show(f)
            show(g)
    """
    from numpy import zeros_like, asarray, zeros, float32
    measurement = measurement.lower()
    is_data = option.lower() == 'data'
    if is_data:
        y = []
    else:
        type = f.dtype
        if measurement in ('mean', 'std'):
            type = float32
        y = zeros(labels.shape)
    fdata = f.ravel()
    labelsravel = labels.ravel()
    for obj_id in range(labelsravel.max()):
        pixels = fdata[labelsravel == obj_id + 1]
        if measurement == 'max':
            val = pixels.max()
        elif measurement == 'min':
            val = pixels.min()
        elif measurement == 'sum':
            val = pixels.sum()
        elif measurement == 'mean':
            val = pixels.mean()
        elif measurement == 'std':
            val = pixels.std()
        elif measurement == 'std1':
            print("'std1' is not implemented")
        else:
            print(
                "measurement should be 'max', 'min', 'mean', 'sum', 'std', 'std1'."
                )
        if is_data:
            y.append(val)
        else:
            y += val * (labels == obj_id + 1)
    return asarray(y)
