def combine(method, images, masks=None, dtype=None,
            region=None, zeros=None, scales=None, weights=None,
            include_variance=True, datamodel=None, method_name=""):
    """Combine HDUList objects using algorithm 'method'"""

    import numina.datamodel

    nimages = len(images)

    # processing masks
    if masks is None:
        intl_masks = [None for _ in images]
    elif isinstance(masks, Extension):
        intl_masks = [img[masks.name].data for img in images]
    elif isinstance(masks, collections.Sequence):
        intl_masks = [mask[0].data for mask in masks]
    else:
        raise TypeError('mask in invalid')

    if len(intl_masks) != nimages:
        raise TypeError('len(masks) != len(images)')

    # Processing region
    if region is None:
        region = Ellipsis # meaning arr[...]

    # processing zeros, scales and weights
    num_values = dict(zeros=None, scales=None, weights=None)
    arg_values = dict(zeros=zeros, scales=scales, weights=weights)

    for key in num_values:
        arg = arg_values[key] # locals()[key]
        value_arg, arg_is_func = _inspect_method(arg)
        if arg_is_func:
            num_values[key] = [value_arg(img, mask, region) for img, mask in zip(images, intl_masks)]
        else:
            num_values[key] = value_arg
        if num_values[key] is not None:
            # Intl combine works only with np arrays
            num_values[key] = numpy.array(num_values[key])

    arrays = [hdu[0].data for hdu in images]
    if masks is None:
        fmasks = None
    else:
        fmasks = intl_masks

    out = C.generic_combine(
        method, arrays, masks=fmasks, dtype=dtype, out=None,
        zeros=num_values['zeros'], scales=num_values['scales'],
        weights=num_values['weights'])

    # Build HDUList
    headers = [img[0].header for img in images]

    base_header = headers[0].copy()
    hdu1 = fits.PrimaryHDU(data=out[0], header=base_header)
    list_of_hdu = [hdu1]

    #_logger.debug('update result header')
    prolog = ""
    if prolog:
        #_logger.debug('write prolog')
        hdu1.header['history'] = prolog

    if method_name != "":
        #_logger.info("Combined %d images using '%s'", nimages, method.__name__)
        hdu1.header['history'] = "Combined %d images using '%s'" % (nimages, method)
    else:
        hdu1.header['history'] = "Combined %d images" % (nimages,)

    hdu1.header['history'] = 'Combination time {}'.format(datetime.datetime.utcnow().isoformat())

    if datamodel is None:
        datamodel = numina.datamodel.DataModel()

    for idx, img in enumerate(images, start=1):
        hdu1.header['history'] = "Image{} {}".format(idx, datamodel.get_imgid(img))

    prevnum = base_header.get('NUM-NCOM', 1)

    hdu1.header['NUM-NCOM'] = prevnum * nimages
    hdu1.header['UUID'] = str(uuid.uuid1())

    # Headers of last image
    # hdu1.header['TSUTC2'] = headers[-1]['TSUTC2']

    if include_variance:
        varhdu = fits.ImageHDU(out[1], name='VARIANCE')
        list_of_hdu.append(varhdu)

    num = fits.ImageHDU(out[2].astype('uint8'), name='MAP')
    list_of_hdu.append(num)

    result = fits.HDUList(list_of_hdu)
    return result
Combine HDUList objects using algorithm 'method