def overlay_chunk(src_chunk, dest_chunk):
    # Courtesy of Peter B. (peter.be...ton@gmail.com)
    chunk = numpy.zeros((Level.chunk_size_Z, Level.chunk_size_X, Level.color_depth), dtype=numpy.uint32)
    """

    blend( (sL,sA), (dL,dA) ) = t1 + t2 + t3

    t1 = sL * sA *   dA  / (255^2)
    t2 = dL * dA * (!sA) / (255^2)
    t3 = sL *      (!dA) /  255

    """

    src_channel_R  = src_chunk[:,:,0].astype(numpy.uint32)
    dest_channel_R = dest_chunk[:,:,0].astype(numpy.uint32)
    src_channel_G  = src_chunk[:,:,1].astype(numpy.uint32)
    dest_channel_G = dest_chunk[:,:,1].astype(numpy.uint32)
    src_channel_B  = src_chunk[:,:,2].astype(numpy.uint32)
    dest_channel_B = dest_chunk[:,:,2].astype(numpy.uint32)
    src_channel_A  = src_chunk[:,:,3].astype(numpy.uint32)
    dest_channel_A = dest_chunk[:,:,3].astype(numpy.uint32)

    for i, (s,d) in enumerate(zip(
        (src_channel_R, src_channel_G, src_channel_B),
        (dest_channel_R, dest_channel_G, dest_channel_B)
    )):
        chunk[:,:,i] = (
            #t1 =  sL * sA *   dA  / (255^2)
            (s * src_channel_A * dest_channel_A / 65025) +
            #t2 = dL * dA * (!sA) / (255^2)
            (d * dest_channel_A * (255 - src_channel_A) / 65025) +
            #t3 = sL *      (!dA) /  255
            (s * (255 - dest_channel_A) / 255)
        )
    chunk[:,:,3] = (src_channel_A * dest_channel_A) / 255 + \
        (src_channel_A * (255 - dest_channel_A)) / 255 + \
        (dest_channel_A * (255 - src_channel_A) / 255)

    return chunk.reshape(Level.chunk_size_Z, Level.chunk_size_X, Level.color_depth).astype(numpy.uint8)
dataset/ETHPy150Open aheadley/pynemap/pynemap.py/overlay_chunk/original