How to use the ehtim.image.Image function in ehtim

To help you get started, we’ve selected a few ehtim examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github achael / eht-imaging / ehtim / scattering / stochastic_optics.py View on Github external
Returns:
               out (Image): The ensemble-average scattered image.
            """

        # Inputs an unscattered image and an ensemble-average blurring kernel (2D array); returns the ensemble-average image
        # The pre-computed kernel can optionally be specified (ker)

        if wavelength_cm == None:
            wavelength_cm = C/im.rf*100.0 #Observing wavelength [cm]

        if ker is None:
            ker = self.Ensemble_Average_Kernel(im, wavelength_cm, use_approximate_form)

        Iim = Wrapped_Convolve((im.imvec).reshape(im.ydim, im.xdim), ker)
        out = image.Image(Iim, im.psize, im.ra, im.dec, rf=C/(wavelength_cm/100.0), source=im.source, mjd=im.mjd, pulse=im.pulse)
        if len(im.qvec):
            Qim = Wrapped_Convolve((im.qvec).reshape(im.ydim, im.xdim), ker)
            Uim = Wrapped_Convolve((im.uvec).reshape(im.ydim, im.xdim), ker)
            out.add_qu(Qim, Uim)
        if len(im.vvec):
            Vim = Wrapped_Convolve((im.vvec).reshape(im.ydim, im.xdim), ker)
            out.add_v(Vim)

        return out
github achael / eht-imaging / ehtim / imaging / pol_imager_utils.py View on Github external
outcv = unpack_poltuple(res.x, xtuple, nimage, pol_solve)
    if pol_prim == "amp_phase":
        outcut = mcv(outcv)  #change of variables
    else:
        outcut = outcv

    if np.any(np.invert(embed_mask)): 
        out = embed_pol(out, embed_mask) #embed
    else:
        out = outcut

    iimage = out[0]
    qimage = make_q_image(out, pol_prim)
    uimage = make_u_image(out, pol_prim)

    outim = image.Image(iimage.reshape(Prior.ydim, Prior.xdim), Prior.psize,
                         Prior.ra, Prior.dec, rf=Prior.rf, source=Prior.source,
                         mjd=Prior.mjd, pulse=Prior.pulse)
    outim.add_qu(qimage.reshape(Prior.ydim, Prior.xdim), uimage.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" % (chisq1(outcut), chisq2(outcut)))
    print(res.message)

    # Return Image object
    return outim
github achael / eht-imaging / ehtim / imaging / dynamical_imaging.py View on Github external
cur_len = np.sum(embed_mask_List[i])
            log_Frames[i] = embed(x[init_i:(init_i+cur_len)], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Frames[i] = np.exp(log_Frames[i])*(embed_mask_List[i].reshape((N_pixel, N_pixel)))
            init_i += cur_len

        if R_flow['alpha'] != 0.0:
            cur_len = np.sum(embed_mask_List[0]) #assumes all the priors have the same embedding
            Flow_x = embed(x[init_i:(init_i+2*cur_len-1):2],   embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow_y = embed(x[(init_i+1):(init_i+2*cur_len):2], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow = np.transpose([Flow_x.ravel(),Flow_y.ravel()]).reshape((N_pixel, N_pixel,2))
            init_i += 2*cur_len

        if stochastic_optics == True:
            EpsilonList = x[init_i:(init_i + N**2-1)]
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            im_List = [image.Image(Frames[j], Prior.psize, Prior.ra, Prior.dec, rf=Obsdata_List[j].rf, source=Prior.source, mjd=Prior.mjd) for j in range(N_frame)]
            scatt_im_List = [scattering_model.Scatter(im_List[j], Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N), ea_ker = ea_ker[j], sqrtQ=sqrtQ, Linearized_Approximation=True).imvec for j in range(N_frame)] #the list of scattered image vectors
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            init_i += len(EpsilonList)

        s1 = s2 = 0.0

        if alpha_s1 != 0.0:
            s1 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy1, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s1
        if alpha_s2 != 0.0:
            s2 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy2, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s2

        s_dynamic_grad = cm_grad = flux_grad = s_dS = s_dF = 0.0
        if R_dI['alpha'] != 0.0: s_dynamic_grad += RdI_gradient(Frames,**R_dI)*R_dI['alpha']
        if R_dt['alpha'] != 0.0: s_dynamic_grad += Rdt_gradient(Frames, B_dt, **R_dt)*R_dt['alpha']

        if alpha_dS1 != 0.0: s_dS += RdS_gradient(Frames, nprior_embed_List, embed_mask_List, entropy1, norm_reg, beam_size=beam_size, alpha_A=alpha_A)*alpha_dS1
github achael / eht-imaging / ehtim / movie.py View on Github external
Args:
               n (int): the frame number

           Returns:
               (Image): the Image object of the nth frame
        """

        if n<0 or n>=len(self.frames):
            raise Exception("n must be in the range 0 - %i"% self.nframes)

        time = self.start_hr + n * self.framedur/3600

        # Make the primary image
        imarr = self.frames[n].reshape(self.ydim, self.xdim)
        outim = ehtim.image.Image(imarr, self.psize, self.ra, self.dec, self.pa,
                                    polrep=self.polrep, pol_prim=self.pol_prim, time=time,
                                    rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)

        # Copy over the rest of the polarizations
        for pol in list(self._movdict.keys()):
            if pol==self.pol_prim: continue
            polframes = self._movdict[pol]
            if len(polframes):
                polvec = polframes[n]
                polarr = polvec.reshape(self.ydim, self.xdim).copy()
                outim.add_pol_image(polarr, pol)

        return outim
github achael / eht-imaging / ehtim / imaging / dynamical_imaging.py View on Github external
plotcur(res.x, final=True)

    # Print stats
    print ("time: %f s" % (tstop - tstart))
    print ("J: %f" % res.fun)
    print (res.message)

    #Note: the global variables are *not* released to avoid recalculation
    if processes != -1:
        pool.close()

    #Return Frames

    outim = [image.Image(Frames[i].reshape(Prior.ydim, Prior.xdim), Prior.psize,
                         Prior.ra, Prior.dec, rf=Obsdata_List[i].rf, source=Prior.source,
                         mjd=Prior.mjd, pulse=Prior.pulse) for i in range(N_frame)]

    if R_flow['alpha'] == 0.0 and stochastic_optics == False:
        return outim
    else:
        return {'Frames':outim, 'Flow':Flow, 'EpsilonList':EpsilonList }
github achael / eht-imaging / ehtim / imaging / starwarps.py View on Github external
def applyAppxWarp(im, theta, centerTheta, init_x, init_y, flowbasis_x, flowbasis_y, initTheta, method1='phase', method2='phase'):
    
    centerIm, dImg_dTheta = calAppxWarpTerms(im, centerTheta, init_x, init_y, flowbasis_x, flowbasis_y, initTheta, method1=method1, method2=method2)
    out = centerIm.imvec + np.dot(dImg_dTheta, theta - centerTheta)
    outim = image.Image(np.reshape(out, (im.ydim, im.xdim)), im.psize, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd, pulse=im.pulse)
    return outim
github achael / eht-imaging / ehtim / image.py View on Github external
if len(self.qvec)==0 or len(self.uvec)==0:
                    rlvec = []
                    lrvec = []
                else:
                    rlvec = (self.qvec + 1j*self.uvec)
                    lrvec = (self.qvec - 1j*self.uvec)

                imdict = {'RR':rrvec,'LL':llvec,'RL':rlvec,'LR':lrvec}

        # Assemble the new image
        imvec = imdict[pol_prim_out]
        if len(imvec)==0:
            raise Exception("for switch_polrep to %s with pol_prim_out=%s, \n"%(polrep_out,pol_prim_out) +
                            "output image is not defined")

        newim = Image(imvec.reshape(self.ydim,self.xdim), self.psize, self.ra, self.dec, self.pa,
                      polrep=polrep_out, pol_prim=pol_prim_out, time=self.time,
                      rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)

        # Add in any other polarizations
        for pol in list(imdict.keys()):
            if pol==newim.pol_prim: continue
            polvec = imdict[pol]
            if len(polvec):
                newim.add_pol_image(polvec.reshape(self.ydim,self.xdim), pol)

        return newim
github achael / eht-imaging / ehtim / image.py View on Github external
sigma = fwhm_i / (2. * np.sqrt(2. * np.log(2.)))
        sigmap = sigma / self.psize

        # Define a convolution function
        def blur(imarr, sigma):
            if np.any(np.imag(imarr)!=0):
                return blur(np.real(imarr),sigma) + 1j*blur(np.imag(imarr),sigma)
            imarr_blur = filt.gaussian_filter(imarr, (sigma, sigma))
            return imarr_blur

        # Blur the primary image
        imarr = self.imvec.reshape(self.ydim, self.xdim)
        imarr = blur(imarr,sigmap)
        outim = Image(imarr, self.psize, self.ra, self.dec, self.pa,
                      polrep=self.polrep, pol_prim=self.pol_prim, time=self.time,
                      rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)

        # Blur all polarizations and copy over
        for pol in list(self._imdict.keys()):
            if pol==self.pol_prim: continue
            polvec = self._imdict[pol]
            if len(polvec):
                polarr = polvec.reshape(self.ydim, self.xdim)
                if fwhm_pol:
                    print ("Blurring polarization")
                    sigma = fwhm_pol / (2. * np.sqrt(2. * np.log(2.)))
                    sigmap = sigma / self.psize
                    polarr = blur(polarr, sigmap)
                outim.add_pol_image(polarr, pol)
github achael / eht-imaging / ehtim / image.py View on Github external
def im_new(imvec):
            imarr_new = np.array([[im_new_val(imvec, x_idx , y_idx)
                                   for x_idx in np.arange(0, -xdim_new, -1)]
                                   for y_idx in np.arange(0, -ydim_new, -1)])
            return imarr_new

        # Compute new primary image vector
        imarr_new = im_new(self.imvec)

        # Normalize
        scaling = np.sum(self.imvec) / np.sum(imarr_new)
        imarr_new *= scaling

        # Make new image
        outim = Image(imarr_new, psize_new, self.ra, self.dec, self.pa,
                      polrep=self.polrep, pol_prim=self.pol_prim, time=self.time,
                      rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)

        # Interpolate all polarizations and copy over
        for pol in list(self._imdict.keys()):
            if pol==self.pol_prim: continue
            polvec = self._imdict[pol]
            if len(polvec):
                polarr_new = im_new(polvec)
                polarr_new *= scaling
                outim.add_pol_image(polarr_new, pol)

        return outim