<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def merge2images(hdu1, hdu2, debugplot):
 """Merge 2 EMIR images, averaging the common region. Parameters hdu1 : HDUList object Input image #1. hdu2 : HDUList object Input image #2. debugplot : int Determines whether intermediate computations and/or plots are displayed. The valid codes are defined in numina.array.display.pause_debugplot. Returns ------- image_merged : HDUList object Output image. """


    # check image dimensions
    image_header = hdu1[0].header
    image_header_ = hdu2[0].header
    #
    naxis = image_header['naxis']
    naxis_ = image_header_['naxis']
    if naxis != naxis_:
        raise ValueError('Incompatible NAXIS values: {}, {}'.format(
            naxis, naxis_
        ))
    naxis1 = image_header['naxis1']
    naxis1_ = image_header_['naxis1']
    if naxis1 != naxis1_:
        raise ValueError('Incompatible NAXIS1 values: {}, {}'.format(
            naxis1, naxis1_
        ))
    if naxis == 1:
        naxis2 = 1
    elif naxis == 2:
        naxis2 = image_header['naxis2']
        naxis2_ = image_header_['naxis2']
        if naxis2 != naxis2_:
            raise ValueError('Incompatible NAXIS2 values: {}, {}'.format(
                naxis2, naxis2_
            ))
    else:
        raise ValueError('Unexpected NAXIS value: {}'.format(naxis))

    # initialize output array
    image2d_merged = np.zeros((naxis2, naxis1))

    # main loop
    data1 = hdu1[0].data
    data2 = hdu2[0].data
    for i in range(naxis2):
        sp1 = data1[i, :]
        sp2 = data2[i, :]
        jmin1, jmax1 = find_pix_borders(sp1, sought_value=0)
        jmin2, jmax2 = find_pix_borders(sp2, sought_value=0)
        image2d_merged[i, :] = sp1 + sp2
        useful1 = (jmin1 != -1) and (jmax1 != naxis1)
        useful2 = (jmin2 != -1) and (jmax2 != naxis1)
        if useful1 and useful2:
            jmineff = max(jmin1, jmin2)
            jmaxeff = min(jmax1, jmax2)
            image2d_merged[i, jmineff:(jmaxeff+1)] /= 2

    # return result
    image_merged = fits.PrimaryHDU(data=image2d_merged.astype(np.float32),
                                   header=image_header)

    return image_merged