def getGriddata(x,y,z,extend):
''' data x,y,z and boundbox to print '''
(xmin,xmax,ymin,ymax)=extend
grid_y, grid_x = np.mgrid[xmin:xmax:(xmax-xmin)*10j, ymin:ymax:(ymax-ymin)*10j]
points=[]
for i in range(x.shape[0]):
points.append([y[i],x[i]])
values=z
# see http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
from scipy.interpolate import griddata
# grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
# grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='cubic')
return grid_z2
python类griddata()的实例源码
def computeInterpolatedICImg(self,roi=None):
"""Interpolates ICs back onto 2D image.
Uses ``scipy.interpolate.griddata``, see also http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
If ``roi`` is specified, will only interpolate nodes of this ROI.
Keyword Args:
roi (pyfrp.subclasses.pyfrp_ROI.ROI): A PyFRAP ROI.
Returns:
tuple: Tuple containing:
* X (numpy.ndarray): Meshgrid x-coordinates.
* Y (numpy.ndarray): Meshgrid y-coordinates.
* interpIC (numpy.ndarray): Interpolated ICs.
"""
X,Y,interpIC=self.computeInterpolatedSolutionToImg(self.IC,roi=roi)
return X,Y,interpIC
def interpolate_apertures(self, aperture_centers, aperture_means):
"""
This function ...
:param aperture_centers:
:param aperture_means:
:return:
"""
# Inform the user
log.info("Interpolating between the mean values of each aperture to fill the sky frame ...")
x_values = np.array([center.x for center in aperture_centers])
y_values = np.array([center.y for center in aperture_centers])
x_ticks = np.arange(0, self.frame.xsize, 1)
y_ticks = np.arange(0, self.frame.ysize, 1)
z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks)
# Set the sky frame
self.sky = Frame(z_grid)
# -----------------------------------------------------------------
def interpolate_apertures(self, aperture_centers, aperture_means):
"""
This function ...
:param aperture_centers:
:param aperture_means:
:return:
"""
# Inform the user
log.info("Interpolating between the mean values of each aperture to fill the sky frame ...")
x_values = np.array([center.x for center in aperture_centers])
y_values = np.array([center.y for center in aperture_centers])
x_ticks = np.arange(0, self.frame.xsize, 1)
y_ticks = np.arange(0, self.frame.ysize, 1)
z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks)
# Set the sky frame
self.sky = Frame(z_grid)
# -----------------------------------------------------------------
def map_data_cosinetolinear(self,values_on_cosine_grid,Ny,a,b):
"""
Map data on a cosine grid to a linear grid
"""
ycells = np.linspace(0, Ny, Ny)
ylin = np.linspace(a, b, Ny)
ycos = 0.5*(b+a) - 0.5*(b-a)*np.cos((ycells*np.pi)/(Ny-1))
#print(ycos.shape,values_on_cosine_grid.shape)
#plt.plot(ycos,values_on_cosine_grid,'x-',label='cosinetolinear Before')
values_on_linear_grid = interp.griddata(ycos, values_on_cosine_grid,
ylin, method='cubic',
fill_value=values_on_cosine_grid[-1])
#values_on_linear_grid = interp2.map_coordinates(values_on_cosine_grid,ycos,output=ylin)
#plt.plot(ylin,values_on_linear_grid,'o-',alpha=0.4,label='cosinetolinear After')
#plt.legend()
#plt.show()
return values_on_linear_grid
def interpolate_in_depth(self,z,zi,Var):
shp = Var.shape
ndepths = z.size
shp1= [shp[0],ndepths]
im, zm = np.meshgrid(np.arange(shp[0]),z)
imi, zmi = np.meshgrid(np.arange(shp[0]),zi)
im = im.T.ravel()
zm= zm.T.ravel()
imi = imi.T.ravel()
zmi= zmi.T.ravel()
Vari = interpolate.griddata((im,zm),Var.ravel(),(imi,zmi)).reshape(shp[0],zi.size)
return Vari
def atoms_to_density_map(atoms, voxelSZ):
(x, y, z) = atoms[:,1:4].T.copy()
(x_min, x_max) = (x.min(), x.max())
(y_min, y_max) = (y.min(), y.max())
(z_min, z_max) = (z.min(), z.max())
grid_len = max([x_max - x_min, y_max - y_min, z_max - z_min])
R = np.int(np.ceil(grid_len / voxelSZ))
if R % 2 == 0:
R += 1
msg = "Length of particle (voxels), %d"%(R)
logging.info(msg)
elec_den = atoms[:,0].copy()
#x = (x-x_min)/voxelSZ
#y = (y-y_min)/voxelSZ
#z = (z-z_min)/voxelSZ
x = (x-0.5*(x_max+x_min-grid_len))/voxelSZ
y = (y-0.5*(y_max+y_min-grid_len))/voxelSZ
z = (z-0.5*(z_max+z_min-grid_len))/voxelSZ
bins = np.arange(R+1)
all_bins = np.vstack((bins,bins,bins))
coords = np.asarray([x,y,z]).T
#(h, h_edges) = np.histogramdd(coords, bins=all_bins, weights=elec_den)
#return h
#return griddata(coords, elec_den, np.mgrid[0:R,0:R,0:R].T, method='linear', fill_value=0.).T
integ = np.floor(coords)
frac = coords - integ
ix = integ[:,0]; iy = integ[:,1]; iz = integ[:,2]
fx = frac[:,0]; fy = frac[:,1]; fz = frac[:,2]
cx = 1. - fx; cy = 1. - fy; cz = 1. - fz
h_total = np.histogramdd(np.asarray([ix,iy,iz]).T, weights=elec_den*cx*cy*cz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix,iy,iz+1]).T, weights=elec_den*cx*cy*fz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix,iy+1,iz]).T, weights=elec_den*cx*fy*cz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix,iy+1,iz+1]).T, weights=elec_den*cx*fy*fz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix+1,iy,iz]).T, weights=elec_den*fx*cy*cz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix+1,iy,iz+1]).T, weights=elec_den*fx*cy*fz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix+1,iy+1,iz]).T, weights=elec_den*fx*fy*cz, bins=all_bins)[0]
h_total += np.histogramdd(np.asarray([ix+1,iy+1,iz+1]).T, weights=elec_den*fx*fy*fz, bins=all_bins)[0]
return h_total
def get_ddd_by_energy(self, energy, points):
""" TODO: documentation
"""
try:
from scipy import interpolate
except ImportError as e:
logger.error("Please install scipy to be able to use spline-based interpolation")
raise e
ev_point = np.array([points, [energy] * len(points)])
return interpolate.griddata(self.points, self.ddd_list, np.transpose(ev_point), method='linear')
def plot_interpolated(self, aperture_centers, aperture_means):
"""
This function ...
:param aperture_centers:
:param aperture_means:
:return:
"""
x_values = np.array([center.x for center in aperture_centers])
y_values = np.array([center.y for center in aperture_centers])
x_ticks = np.arange(0, self.frame.xsize, 1)
y_ticks = np.arange(0, self.frame.ysize, 1)
z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks)
self.sky = Frame(z_grid)
from matplotlib.backends import backend_agg as agg
from matplotlib import cm
# plot
#fig = Figure() # create the figure
fig = plt.figure()
agg.FigureCanvasAgg(fig) # attach the rasterizer
ax = fig.add_subplot(1, 1, 1) # make axes to plot on
ax.set_title("Interpolated Contour Plot of Experimental Data")
ax.set_xlabel("X")
ax.set_ylabel("Y")
cmap = cm.get_cmap("hot") # get the "hot" color map
contourset = ax.contourf(x_ticks, y_ticks, z_grid, 10, cmap=cmap)
cbar = fig.colorbar(contourset)
cbar.set_ticks([0, 100])
fig.axes[-1].set_ylabel("Z") # last axes instance is the colorbar
plt.show()
# -----------------------------------------------------------------
def plot_interpolated(self, aperture_centers, aperture_means):
"""
This function ...
:param aperture_centers:
:param aperture_means:
:return:
"""
x_values = np.array([center.x for center in aperture_centers])
y_values = np.array([center.y for center in aperture_centers])
x_ticks = np.arange(0, self.frame.xsize, 1)
y_ticks = np.arange(0, self.frame.ysize, 1)
z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks)
self.sky = Frame(z_grid)
from matplotlib.backends import backend_agg as agg
from matplotlib import cm
# plot
#fig = Figure() # create the figure
fig = plt.figure()
agg.FigureCanvasAgg(fig) # attach the rasterizer
ax = fig.add_subplot(1, 1, 1) # make axes to plot on
ax.set_title("Interpolated Contour Plot of Experimental Data")
ax.set_xlabel("X")
ax.set_ylabel("Y")
cmap = cm.get_cmap("hot") # get the "hot" color map
contourset = ax.contourf(x_ticks, y_ticks, z_grid, 10, cmap=cmap)
cbar = fig.colorbar(contourset)
cbar.set_ticks([0, 100])
fig.axes[-1].set_ylabel("Z") # last axes instance is the colorbar
plt.show()
# -----------------------------------------------------------------
def warp_image(im, flow):
"""
Use optical flow to warp image to the next
:param im: image to warp
:param flow: optical flow
:return: warped image
"""
from scipy import interpolate
image_height = im.shape[0]
image_width = im.shape[1]
flow_height = flow.shape[0]
flow_width = flow.shape[1]
n = image_height * image_width
(iy, ix) = np.mgrid[0:image_height, 0:image_width]
(fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
fx += flow[:,:,0]
fy += flow[:,:,1]
mask = np.logical_or(fx <0 , fx > flow_width)
mask = np.logical_or(mask, fy < 0)
mask = np.logical_or(mask, fy > flow_height)
fx = np.minimum(np.maximum(fx, 0), flow_width)
fy = np.minimum(np.maximum(fy, 0), flow_height)
points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
warp = np.zeros((image_height, image_width, im.shape[2]))
for i in range(im.shape[2]):
channel = im[:, :, i]
plt.imshow(channel, cmap='gray')
values = channel.reshape(n, 1)
new_channel = interpolate.griddata(points, values, xi, method='cubic')
new_channel = np.reshape(new_channel, [flow_height, flow_width])
new_channel[mask] = 1
warp[:, :, i] = new_channel.astype(np.uint8)
return warp.astype(np.uint8)
def map_data_lineartocosine(self,values_on_linear_grid,Ny,a,b):
"""
Map data on a linear grid to a cosine grid
"""
ycells = np.linspace(0, Ny, Ny)
ylin = np.linspace(a, b, Ny)
ycos = 0.5*(b+a) - 0.5*(b-a)*np.cos((ycells*np.pi)/(Ny-1))
plt.plot(ylin,values_on_linear_grid,'o-',alpha=0.4,label='lineartocosine Before')
values_on_cosine_grid = interp.griddata(ylin, values_on_linear_grid,
ycos, method='cubic',
fill_value=values_on_linear_grid[-1])
plt.plot(ycos,values_on_cosine_grid,'x-',label='lineartocosine After')
plt.legend()
plt.show()
return values_on_cosine_grid
def interpolate_xml_array(data, low_res_coords, full_res_size):
"""Interpolate arbitrary size dataset to a full sized grid."""
from scipy.interpolate import griddata
grid_x, grid_y = np.mgrid[0:full_res_size[0], 0:full_res_size[1]]
x, y = low_res_coords
return griddata(np.vstack((np.asarray(y), np.asarray(x))).T,
data,
(grid_x, grid_y),
method='linear')
def cars_profile(filename, doy, latitude, longitude, depth):
"""
For now only the nearest value
For now only for one position, not an array of positions
longitude 0-360
"""
assert np.size(doy) == 1
assert np.size(latitude) == 1
assert np.size(longitude) == 1
#assert np.size(depth) == 1
assert (longitude >= 0) & (longitude <= 360)
assert depth >= 0
nc = netCDF4.Dataset(filename)
t = 2 * np.pi * doy/366
# Improve this. Optimize to get only necessary Z
Z = slice(0, nc['depth'].size)
I = np.absolute(nc['lat'][:] - latitude).argmin()
J = np.absolute(nc['lon'][:] - longitude).argmin()
# Not efficient, but it works
assert (nc['depth'][:64] == nc['depth_ann'][:]).all()
assert (nc['depth'][:55] == nc['depth_semiann'][:]).all()
value = nc['mean'][:, I, J]
value[:64] += nc['an_cos'][Z, I, J] * np.cos(t) + \
nc['an_sin'][:, I, J] * np.sin(t)
value[:55] += nc['sa_cos'][Z, I, J] * np.cos(2*t) + \
nc['sa_sin'][:, I, J] * np.sin(2*t)
value = value
output = {'depth': np.asanyarray(depth)}
from scipy.interpolate import griddata
output['value'] = griddata(nc['depth'][Z], value[Z], depth)
for v in ['std_dev']:
output[v] = griddata(nc['depth'][Z], nc[v][Z, I, J], depth)
return output
def interpolate_coarser2finer_2D(self,x,y,xi,yi,Var,I):
aux = np.array([Var[i[0],i[1]] for i in I]) #T in a.ann_i sites
Vari = interpolate.griddata((x,y),aux,(xi,yi),method='linear')
return Vari
def interpolate_coarser2finer_2D_depth(self,x,y,xi,yi,Var,I):
aux = np.array([Var[:,i[0],i[1]] for i in I]) #T in a.ann_i sites
d = aux.shape[1]
Vari= np.zeros((xi.shape[0],d))
for i in range(d):
Vari[:,i] = interpolate.griddata((x,y),aux[:,i],(xi,yi),method='linear')
return Vari
def interpola(self,x,y,xi,yi,var,ndepths):
Vari = np.zeros((xi.shape[0],ndepths))
for k in range(ndepths):
Vari[:,k] = interpolate.griddata((x,y),var[k,:,:].ravel(),(xi,yi),method='nearest')
Vari[Vari==Vari.min()]=Vari.max()
return Vari
def coupe(X, Y, Z, x, y):
"""
Calcule une coupe et renvoie le Z qui correspond
"""
X=X.flatten()
Y=Y.flatten()
Z=Z.flatten()
z = griddata((X, Y), Z, (x, y))
return z
def warp_image(im, flow):
"""
Use optical flow to warp image to the next
:param im: image to warp
:param flow: optical flow
:return: warped image
"""
from scipy import interpolate
image_height = im.shape[0]
image_width = im.shape[1]
flow_height = flow.shape[0]
flow_width = flow.shape[1]
n = image_height * image_width
(iy, ix) = np.mgrid[0:image_height, 0:image_width]
(fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
fx += flow[:,:,0]
fy += flow[:,:,1]
mask = np.logical_or(fx <0 , fx > flow_width)
mask = np.logical_or(mask, fy < 0)
mask = np.logical_or(mask, fy > flow_height)
fx = np.minimum(np.maximum(fx, 0), flow_width)
fy = np.minimum(np.maximum(fy, 0), flow_height)
points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
warp = np.zeros((image_height, image_width, im.shape[2]))
for i in range(im.shape[2]):
channel = im[:, :, i]
plt.imshow(channel, cmap='gray')
values = channel.reshape(n, 1)
new_channel = interpolate.griddata(points, values, xi, method='cubic')
new_channel = np.reshape(new_channel, [flow_height, flow_width])
new_channel[mask] = 1
warp[:, :, i] = new_channel.astype(np.uint8)
return warp.astype(np.uint8)
def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3):
self.A = A
self.h = h
self.Asp = 2*self.h**2/self.A
self.rho = rho
# Input lift and drag data
self.n = len(alpha)
self.alpha = alpha
self.CL = CL
self.CD = CD
self.k_spline = k_spline
# -------- Spline interpolation -----------------------------
if len(self.alpha.shape) == 1:
self.interpolationMethod = 'spline'
self.nrControlParameters = 1
self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline)
self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline)
elif len(self.alpha.shape) == 2:
self.interpolationMethod = 'griddata'
self.nrControlParameters = 2
self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing)
self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)
def CLCD(self, alpha):
if self.interpolationMethod == 'spline':
CL = interpolate.splev(alpha, self.CL_spline)
CD = interpolate.splev(alpha, self.CD_spline)
elif self.interpolationMethod == 'Rbf':
CL = self.CL_RbfModel(alpha[0], alpha[1])
CD = self.CD_RbfModel(alpha[0], alpha[1])
elif self.interpolationMethod == 'griddata':
CL = interpolate.griddata(self.alpha, self.CL, alpha)
CD = interpolate.griddata(self.alpha, self.CD, alpha)
if self.nrControlParameters == 2:
return float(CL), float(CD)
else:
return CL, CD
def mesh2grid(v, mesh):
""" Interpolates from an unstructured coordinates (mesh) to a structured
coordinates (grid)
"""
x = mesh[:,0]
z = mesh[:,1]
lx = x.max() - x.min()
lz = z.max() - z.min()
nn = v.size
nx = np.around(np.sqrt(nn*lx/lz))
nz = np.around(np.sqrt(nn*lz/lx))
dx = lx/nx
dz = lz/nz
# construct structured grid
x = np.linspace(x.min(), x.max(), nx)
z = np.linspace(z.min(), z.max(), nz)
X, Z = np.meshgrid(x, z)
grid = stack(X.flatten(), Z.flatten())
# interpolate to structured grid
V = _interp.griddata(mesh, v, grid, 'linear')
# workaround edge issues
if np.any(np.isnan(V)):
W = _interp.griddata(mesh, v, grid, 'nearest')
for i in np.where(np.isnan(V)):
V[i] = W[i]
V = np.reshape(V, (nz, nx))
return V, grid
def grid2mesh(V, grid, mesh):
""" Interpolates from structured coordinates (grid) to unstructured
coordinates (mesh)
"""
return _interp.griddata(grid, V.flatten(), mesh, 'linear')
def dic(self):
r""" Returns the corrected Deviance Information Criterion (DIC) for all chains loaded into ChainConsumer.
If a chain does not have a posterior, this method will return `None` for that chain. **Note that
the DIC metric is only valid on posterior surfaces which closely resemble multivariate normals!**
Formally, we follow Liddle (2007) and first define *Bayesian complexity* as
.. math::
p_D = \bar{D}(\theta) - D(\bar{\theta}),
where :math:`D(\theta) = -2\ln(P(\theta)) + C` is the deviance, where :math:`P` is the posterior
and :math:`C` a constant. From here the DIC is defined as
.. math::
DIC \equiv D(\bar{\theta}) + 2p_D = \bar{D}(\theta) + p_D.
Returns
-------
list[float]
A list of all the DIC values - one per chain, in the order in which the chains were added.
References
----------
[1] Andrew R. Liddle, "Information criteria for astrophysical model selection", MNRAS (2007)
"""
dics = []
dics_bool = []
for i, chain in enumerate(self.parent.chains):
p = chain.posterior
if p is None:
dics_bool.append(False)
self._logger.warn("You need to set the posterior for chain %s to get the DIC" % chain.name)
else:
dics_bool.append(True)
num_params = chain.chain.shape[1]
means = np.array([np.average(chain.chain[:, ii], weights=chain.weights) for ii in range(num_params)])
d = -2 * p
d_of_mean = griddata(chain.chain, d, means, method='nearest')[0]
mean_d = np.average(d, weights=chain.weights)
p_d = mean_d - d_of_mean
dic = mean_d + p_d
dics.append(dic)
if len(dics) > 0:
dics -= np.min(dics)
dics_fin = []
i = 0
for b in dics_bool:
if not b:
dics_fin.append(None)
else:
dics_fin.append(dics[i])
i += 1
return dics_fin
timeslice.py 文件源码
项目：algorithm-reference-library
作者: SKA-ScienceDataProcessor
项目源码
文件源码
阅读 26
收藏 0
点赞 0
评论 0
def predict_timeslice_single(vis: Visibility, model: Image, predict=predict_2d_base, **kwargs) -> Visibility:
""" Predict using a single time slices.
This fits a single plane and corrects the image geometry.
:param vis: Visibility to be predicted
:param model: model image
:return: resulting visibility (in place works)
"""
log.debug("predict_timeslice_single: predicting using single time slice")
inchan, inpol, ny, nx = model.shape
vis.data['vis'] *= 0.0
if not isinstance(vis, Visibility):
avis = coalesce_visibility(vis, **kwargs)
else:
avis = vis
# Fit and remove best fitting plane for this slice
avis, p, q = fit_uvwplane(avis, remove=False)
# Calculate nominal and distorted coordinate systems. We will convert the model
# from nominal to distorted before predicting.
workimage = copy_image(model)
# Use griddata to do the conversion. This could be improved. Only cubic is possible in griddata.
# The interpolation is ok for invert since the image is smooth but for clean images the
# interpolation is particularly poor, leading to speckle in the residual image.
lnominal, mnominal, ldistorted, mdistorted = lm_distortion(model, -p, -q)
for chan in range(inchan):
for pol in range(inpol):
workimage.data[chan, pol, ...] = \
griddata((mnominal.flatten(), lnominal.flatten()),
values=workimage.data[chan, pol, ...].flatten(),
xi=(mdistorted.flatten(), ldistorted.flatten()),
method='cubic',
fill_value=0.0,
rescale=True).reshape(workimage.data[chan, pol, ...].shape)
vis = predict(vis, workimage, **kwargs)
return vis
timeslice.py 文件源码
项目：algorithm-reference-library
作者: SKA-ScienceDataProcessor
项目源码
文件源码
阅读 22
收藏 0
点赞 0
评论 0
def invert_timeslice_single(vis: Visibility, im: Image, dopsf, normalize=True, **kwargs) -> (Image, numpy.ndarray):
"""Process single time slice
Extracted for re-use in parallel version
:param vis: Visibility to be inverted
:param im: image template (not changed)
:param dopsf: Make the psf instead of the dirty image
:param normalize: Normalize by the sum of weights (True)
"""
inchan, inpol, ny, nx = im.shape
if not isinstance(vis, Visibility):
avis = coalesce_visibility(vis, **kwargs)
else:
avis = vis
log.debug("invert_timeslice_single: inverting using single time slice")
avis, p, q = fit_uvwplane(avis, remove=False)
workimage, sumwt = invert_2d_base(avis, im, dopsf, normalize=normalize, **kwargs)
finalimage = create_empty_image_like(im)
# Use griddata to do the conversion. This could be improved. Only cubic is possible in griddata.
# The interpolation is ok for invert since the image is smooth.
# Calculate nominal and distorted coordinates. The image is in distorted coordinates so we
# need to convert back to nominal
lnominal, mnominal, ldistorted, mdistorted = lm_distortion(workimage, -p, -q)
for chan in range(inchan):
for pol in range(inpol):
finalimage.data[chan, pol, ...] = \
griddata((mdistorted.flatten(), ldistorted.flatten()),
values=workimage.data[chan, pol, ...].flatten(),
method='cubic',
xi=(mnominal.flatten(), lnominal.flatten()),
fill_value=0.0,
rescale=True).reshape(finalimage.data[chan, pol, ...].shape)
return finalimage, sumwt
def get_ddd_grid(self, energy_list, n):
""" TODO: documentation
"""
energy = []
dist = []
data = []
ddd_e = self.ddd.keys()
ddd_e = sorted(ddd_e)
for e in energy_list:
idx = np.where((np.array(ddd_e) >= e))[0][0] - 1
d_lower = self.ddd[ddd_e[idx]]
d_upper = self.ddd[ddd_e[idx + 1]]
lower_idx = np.where(max(d_lower[1, :]) == d_lower[1, :])[0][0]
upper_idx = np.where(max(d_upper[1, :]) == d_upper[1, :])[0][0]
offset = 1 / (ddd_e[idx + 1] - ddd_e[idx]) * (e - ddd_e[idx + 1])
x_offset = (d_upper[0, upper_idx] - d_lower[0, lower_idx]) * offset
y_offset = 1 + (1 - d_upper[1, upper_idx] / d_lower[1, lower_idx]) * offset
depth = d_upper[0, :] + x_offset
ddd = d_upper[1, :] * y_offset
xi = np.linspace(0, depth[-1], n)
spl = RegularInterpolator(x=depth, y=ddd)
data.extend(spl(xi))
dist.extend(xi)
energy.extend([e] * n)
out = [dist, energy, data]
return np.reshape(np.transpose(out), (len(energy_list), n, 3))
# TODO why it is not used ?
# ddd_list = []
# energy = []
# dist = []
# point = []
# for e in energy_list:
# dist.extend(np.linspace(0, self.get_dist(e), n))
# energy.extend([e] * n)
# point.append(dist)
# point.append(energy)
# data = interpolate.griddata(self.points,
# self.ddd_list,
# np.transpose(point),
# method='linear')
# out = [dist, energy, data]
# return np.reshape(np.transpose(out), (len(energy_list), n, 3))
def computeInterpolatedSolutionToImg(self,vals,roi=None,method='linear',res=None):
"""Interpolates solution back onto 2D image.
Uses ``scipy.interpolate.griddata``, see also http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
If ``roi`` is specified, will only interpolate nodes of this ROI.
For more details about interpolation methods, check out
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html .
Keyword Args:
vals (numpy.ndarray): Solution to be interpolated.
roi (pyfrp.subclasses.pyfrp_ROI.ROI): A PyFRAP ROI.
method (str): Interpolation method.
fillVal (float): Value applied outside of ROI.
res (int): Resolution of resulting images in pixels.
Returns:
tuple: Tuple containing:
* X (numpy.ndarray): Meshgrid x-coordinates.
* Y (numpy.ndarray): Meshgrid y-coordinates.
* interpIC (numpy.ndarray): Interpolated solution.
"""
#Get image resolution and center of geometry
if res==None:
res=self.ICimg.shape[0]
#Build Empty Img
X,Y=np.meshgrid(np.arange(res),np.arange(res))
#Get cellcenters
x,y,z=self.mesh.getCellCenters()
##print x
##print y
##print z
##print roi.meshIdx
#print max(roi.meshIdx)
#print len(x)
if roi!=None:
xInt=x[roi.meshIdx]
yInt=y[roi.meshIdx]
val=vals[roi.meshIdx]
else:
xInt=x
yInt=y
val=vals
interpIC=interp.griddata((xInt,yInt),val,(X,Y),method=method)
return X,Y,interpIC
def findClosestVector(point, arr_shape=None, pixel_shape=None, xyorig=None):
'''
Finds the closest vector of array coordinates (x, y) from an input vector of pixel coordinates (x, y).
Parameters:
point : tuple
Original point of interest in pixel units, order of (x,y)
arr_shape : tuple
Shape of data array in (x,y) order
pixel_shape : tuple
Shape of image in pixels in (x,y) order
xyorig : str
Indicates the origin point of coordinates. Set to "relative" switches to an array coordinate
system relative to galaxy center. Default is absolute array coordinates (x=0, y=0) = upper left corner
Returns:
minind : tuple
A tuple of array coordinates in x, y order
'''
# set as numpy arrays
arr_shape = np.array(arr_shape, dtype=int)
pixel_shape = np.array(pixel_shape, dtype=int)
# compute midpoints
xmid, ymid = arr_shape/2
xpixmid, ypixmid = pixel_shape/2
# default absolute array coordinates
xcoords = np.array([0, arr_shape[0]], dtype=int)
ycoords = np.array([0, arr_shape[1]], dtype=int)
# split x,y coords and pixel coords
x1, x2 = xcoords
y1, y2 = ycoords
xpix, ypix = pixel_shape
# build interpolates between array coordinates and pixel coordinates
points = [[x1, y1], [x1, y2], [xmid, ymid], [x2, y1], [x2, y2]]
values = [[0, ypix], [0, 0], [xpixmid, ypixmid], [xpix, ypix], [xpix, 0]] # full image
# values = [[xpixmid-xmid, ypixmid+ymid], [xpixmid-xmid, ypixmid-ymid], [xpixmid, ypixmid], [xpixmid+xmid, ypixmid+ymid], [xpixmid+xmid, ypixmid-ymid]] # pixels based on arr_shape
#values = [[xpixmid-x2, ypixmid+y2], [xpixmid-x2, ypixmid-y2], [xpixmid, ypixmid], [xpixmid+x2, ypixmid+y2], [xpixmid+x2, ypixmid-y2]] # pixels based on arr_shape
# make 2d array of array indices in absolute or (our) relative coordindates
arrinds = np.mgrid[x1:x2, y1:y2].swapaxes(0, 2).swapaxes(0, 1)
# interpolate a new 2d pixel coordinate array
final = griddata(points, values, arrinds)
# find minimum array vector closest to input coordinate point
diff = np.abs(point - final)
prod = diff[:, :, 0]*diff[:, :, 1]
minind = np.unravel_index(prod.argmin(), arr_shape)
# toggle relative array coordinates
if xyorig in ['relative', 'center']:
minind = np.array(minind, dtype=int)
xmin = minind[0] - xmid
ymin = ymid - minind[1]
minind = (xmin, ymin)
return minind
def interpolate_missing_data(Y):
"""
Interpolate any missing data using nearest neighbor interpolation.
Missing data is identified as entries with values NaN
Input:
Y np.ndarray (3D)
movie, raw data in 3D format (d1 x d2 x T)
Outputs:
Y np.ndarray (3D)
movie, data with interpolated entries (d1 x d2 x T)
coor list
list of interpolated coordinates
"""
coor=[];
if np.any(np.isnan(Y)):
raise Exception('The algorithm has not been tested with missing values (NaNs). Remove NaNs and rerun the algorithm.')
# need to
for idx,row in enumerate(Y):
nans=np.where(np.isnan(row))[0]
n_nans=np.where(~np.isnan(row))[0]
coor.append((idx,nans))
Y[idx,nans]=np.interp(nans, n_nans, row[n_nans])
# mis_data = np.isnan(Y)
# coor = mis_data.nonzero()
# ok_data = ~mis_data
# coor_ok = ok_data.nonzero()
# Yvals=[np.where(np.isnan(Y)) for row in Y]
#
# Yvals = griddata(coor_ok,Y[coor_ok],coor,method='nearest')
# un_t = np.unique(coor[-1])
# coorx = []
# coory = []
# Yvals = []
# for i, unv in enumerate(un_t):
# tm = np.where(coor[-1]==unv)
# coorx.append(coor[0][tm].tolist())
# coory.append(coor[1][tm].tolist())
# Yt = Y[:,:,unv]
# ok = ~np.isnan(Yt)
# coor_ok = ok.nonzero()
# ytemp = griddata(coor_ok,Yt[coor_ok],(coor[0][tm],coor[1][tm]),method='nearest')
# Yvals.append(ytemp)
return Y, coor
#%%
