How to use the ehtim.image function in ehtim

rf = float(file.readline().split()[2]) * 1e9
    xdim = file.readline().split()
    xdim_p = int(xdim[2])
    psize_x = float(xdim[4])*RADPERAS/xdim_p
    ydim = file.readline().split()
    ydim_p = int(ydim[2])
    psize_y = float(ydim[4])*RADPERAS/ydim_p
    file.close()

    if psize_x != psize_y:
        raise Exception("Pixel dimensions in x and y are inconsistent!")

    # Load the data, convert to list format, make object
    datatable = np.loadtxt(filename, dtype=float)
    image = datatable[:,2].reshape(ydim_p, xdim_p)
    outim = ehtim.image.Image(image, psize_x, ra, dec, rf=rf, source=src, mjd=mjd, time=time, pulse=pulse,
                              polrep='stokes',pol_prim='I')

    # Look for Stokes Q and U
    qimage = uimage = vimage = np.zeros(image.shape)
    if datatable.shape[1] == 6:
        qimage = datatable[:,3].reshape(ydim_p, xdim_p)
        uimage = datatable[:,4].reshape(ydim_p, xdim_p)
        vimage = datatable[:,5].reshape(ydim_p, xdim_p)
    elif datatable.shape[1] == 5:
        qimage = datatable[:,3].reshape(ydim_p, xdim_p)
        uimage = datatable[:,4].reshape(ydim_p, xdim_p)

    if np.any((qimage != 0) + (uimage != 0)) and np.any((vimage != 0)):
        #print('Loaded Stokes I, Q, U, and V Images')
        outim.add_qu(qimage, uimage)
        outim.add_v(vimage)

#skip headline
    TableMOD = np.genfromtxt(nameF, skip_header=4)
    ScaleR = 1.
    FluxConst = 1.
    Flux = FluxConst*TableMOD[:,0]
    xPS = ScaleR*TableMOD[:,1]*np.cos(np.pi/2.-(np.pi/180.)*TableMOD[:,2])*(1.e3)*RADPERUAS #to radians
    yPS = ScaleR*TableMOD[:,1]*np.sin(np.pi/2.-(np.pi/180.)*TableMOD[:,2])*(1.e3)*RADPERUAS #to radians
    NumbPoints = np.shape(yPS)[0]

    #set image parameters
    if fov==0:
        MaxR = np.amax(TableMOD[:,1]) #in mas
        fov = 1.*MaxR*(1.e3)*RADPERUAS

    image0 = np.zeros((int(npix),int(npix)))
    im = image.Image(image0, fov/npix, 0., 0., rf=86e9)

    beamMaj = beamPar[0]
    if beamMaj==0:
        beamMaj = 4.*fov/npix

    beamMin = beamPar[1]
    if beamMin==0:
        beamMin = 4.*fov/npix

    beamTh = beamPar[2]

    sigma_maj = beamMaj / (2. * np.sqrt(2. * np.log(2.)))
    sigma_min = beamMin / (2. * np.sqrt(2. * np.log(2.)))
    cth = np.cos(beamTh)
    sth = np.sin(beamTh)
    xfov = im.xdim * im.psize

# load in the data
obs = eh.obsdata.load_uvfits(obsname)

# split the observations based upon the time
obs_List = sw.splitObs(obs)

############## reconstruct movie with no warp field ##############

# initialize the mean and the image covariance for the prior. 
# this can be a single image to be the same mean and covariance for each 
# time, or different for each time by appending an image/matrix for each timestep

# initialize mean
meanImg = []
emptyprior = eh.image.make_square(obs, NPIX, fov)
gaussprior = emptyprior.add_gauss(flux, (fwhm, fwhm, 0, 0, 0))
meanImg.append(gaussprior.copy())

# initialize covariance
imCov = []
imCov.append( sw.gaussImgCovariance_2(meanImg[0], powerDropoff=2.0, frac=1./2.) )

# make the covariance matrix that says how much variation there should be between frames in time 
noiseCov_img = np.eye(npixels)*variance_img_diff

# initialize the flowbasis and get the initTheta which says how to specify no motion for the specified flow basis
init_x, init_y, flowbasis_x, flowbasis_y, initTheta = sw.affineMotionBasis_noTranslation(meanImg[0])

# run StarWarps to find the distribution of the image at each timestep
expVal_t, expVal_t_t, expVal_tm1_t, loglikelihood, apxImgs = sw.computeSuffStatistics(
    meanImg, imCov, obs_List, noiseCov_img, initTheta, init_x, init_y,

def padNewFOV(im, fov_arcseconds):
    
    oldfov = im.psize * im.xdim
    newfov = fov_arcseconds * ehtim.RADPERUAS
    tnpixels = np.ceil(im.xdim * newfov/oldfov).astype('int')
    
    origimg = np.reshape(im.imvec, [im.xdim, im.xdim])
    padimg = np.pad(origimg, ((0,tnpixels-im.xdim), (0,tnpixels-im.xdim)), 'constant')
    padimg = np.roll(padimg, np.floor((tnpixels-im.xdim)/2.).astype('int'), axis=0)
    padimg = np.roll(padimg, np.floor((tnpixels-im.xdim)/2.).astype('int'), axis=1)
    
    return image.Image(padimg.reshape((tnpixels, tnpixels)), im.psize, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd, pulse=im.pulse)

def get_frame(j):
            if type(Unscattered_Movie) == movie.Movie:
                im = image.Image(Unscattered_Movie.frames[j].reshape((N,N)), psize, ra, dec, rf, pulse, source, mjd)
                if len(Unscattered_Movie.qframes) > 0:
                    im.add_qu(Unscattered_Movie.qframes[j].reshape((N,N)), Unscattered_Movie.uframes[j].reshape((N,N)))

                return im
            elif type(Unscattered_Movie) == list:
                return Unscattered_Movie[j]
            else:
                return Unscattered_Movie

def objfunc_scattering(self, minvec):
        """Current stochastic optics objective function.
        """
        N = self.prior_next.xdim

        imvec       = minvec[:N**2]
        EpsilonList = minvec[N**2:]
        if self.transform_next == 'log':
            imvec = np.exp(imvec)

        IM = ehtim.image.Image(imvec.reshape(N,N), self.prior_next.psize,
                               self.prior_next.ra, self.prior_next.dec, self.prior_next.pa,
                               rf=self.obs_next.rf, source=self.prior_next.source, mjd=self.prior_next.mjd)

        #the scattered image vector
        scatt_im = self.scattering_model.Scatter(IM, Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N),
                                                 ea_ker = self._ea_ker, sqrtQ=self._sqrtQ,
                                                 Linearized_Approximation=True).imvec

        # Calculate the chi^2 using the scattered image
        datterm = 0.
        chi2_term_dict = self.make_chisq_dict(scatt_im)
        for dname in sorted(self.dat_term_next.keys()):
            datterm += self.dat_term_next[dname] * (chi2_term_dict[dname] - 1.)

        # Calculate the entropy using the unscattered image
        regterm = 0

tstart = time.time()
    if grads:
        res = opt.minimize(objfunc, xinit, method='L-BFGS-B', jac=objgrad,
                       options=optdict, callback=plotcur)
    else:
        res = opt.minimize(objfunc, xinit, method='L-BFGS-B',
                       options=optdict, callback=plotcur)

    tstop = time.time()

    # Format output
    out = res.x
    if logim: out = np.exp(res.x)
    if np.any(np.invert(embed_mask)): out = embed(out, embed_mask)

    outim = image.Image(out.reshape(Prior.ydim, Prior.xdim), Prior.psize,
                     Prior.ra, Prior.dec, rf=Prior.rf, source=Prior.source,
                     mjd=Prior.mjd, pulse=Prior.pulse)

    if len(Prior.qvec):
        print("Preserving image complex polarization fractions!")
        qvec = Prior.qvec * out / Prior.imvec
        uvec = Prior.uvec * out / Prior.imvec
        outim.add_qu(qvec.reshape(Prior.ydim, Prior.xdim), uvec.reshape(Prior.ydim, Prior.xdim))

    # Print stats
    print("time: %f s" % (tstop - tstart))
    print("J: %f" % res.fun)
    print("Final Chi^2_1: %f Chi^2_2: %f  Chi^2_3: %f" % (chisq1(out[embed_mask]), chisq2(out[embed_mask]), chisq3(out[embed_mask])))
    print(res.message)

    # Return Image object

# Observe the average image
eht = eh.array.load_txt('../arrays/VLBA_86GHz.txt') # Load the observing array
#observing parameters:
tint_sec = 60
tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 0.5e9
#create the observation
obs = scatt.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True)

# Generate an image prior
npix = 55 #This must be odd
prior_fwhm = 300*eh.RADPERUAS # Gaussian size in microarcssec
gaussprior = eh.image.Image(np.zeros((npix,npix)), fov/npix, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
gaussprior = gaussprior.add_gauss(total_flux,(prior_fwhm, prior_fwhm, 0, 0, 0))

# Now try imaging
# First image with no scattering mitigation
imgr = eh.imager.Imager(obs, gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr.make_image_I()
imgr.out_last().display()

# Now try deblurring before imaging
imgr_deblur = eh.imager.Imager(sm.Deblur_obs(obs), gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr_deblur.make_image_I()
imgr_deblur.out_last().display()

# Now image using stochastic optics
imgr_so = eh.imager.Imager(obs, gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr_so.make_image_I_stochastic_optics()