def compute_operation(file1, file2, operation, output, display,
                      args_z1z2, args_bbox, args_keystitle, args_geometry):
    with fits.open(file1) as hdulist:
        image_header1 = hdulist[0].header
        image1 = hdulist[0].data.astype(np.float)
    naxis1 = image_header1['naxis1']
    naxis2 = image_header1['naxis2']
    if display == 'all':
        ximshow_file(file1.name,
                     args_z1z2=args_z1z2, args_bbox=args_bbox,
                     args_keystitle=args_keystitle,
                     args_geometry=args_geometry,
                     debugplot=12)
    with fits.open(file2) as hdulist:
        image_header2 = hdulist[0].header
        image2 = hdulist[0].data.astype(np.float)
    naxis1_ = image_header2['naxis1']
    naxis2_ = image_header2['naxis2']
    if display == 'all':
        ximshow_file(file2.name,
                     args_z1z2=args_z1z2, args_bbox=args_bbox,
                     args_keystitle=args_keystitle,
                     args_geometry=args_geometry,
                     debugplot=12)
    if naxis1 != naxis1_:
        raise ValueError("NAXIS1 values are different.")
    if naxis2 != naxis2_:
        raise ValueError("NAXIS2 values are different.")
    if operation == "+":
        solution = image1 + image2
    elif operation == "-":
        solution = image1 - image2
    elif operation == "*":
        solution = image1 * image2
    elif operation == "/":
        solution = image1 / image2
    else:
        raise ValueError("Unexpected operation=" + str(operation))
    hdu = fits.PrimaryHDU(solution.astype(np.float), image_header1)
    hdu.writeto(output, overwrite=True)
    if display in ['all', 'result']:
        ximshow_file(output.name,
                     args_z1z2=args_z1z2, args_bbox=args_bbox,
                     args_keystitle=args_keystitle,
                     args_geometry=args_geometry,
                     debugplot=12)