mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-27 11:27:23 +08:00
Add surface to mesh test
Fix forward function call
This commit is contained in:
D Fournier
@@ -9,6 +9,7 @@ from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
import scipy.interpolate as interpolation
|
||||
import time
|
||||
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
|
||||
@@ -39,24 +40,33 @@ xr = np.linspace(-99., 99., 40)
|
||||
yr = np.linspace(-49., 49., 20)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
|
||||
d_iter = 1
|
||||
lrl = np.zeros(d_iter)
|
||||
|
||||
|
||||
sclx = 100.
|
||||
dx = 5
|
||||
dx = np.asarray([15., 10., 5., 2.5])
|
||||
|
||||
d_iter = len(dx)
|
||||
l1_r = np.zeros(d_iter)
|
||||
l2_r = np.zeros(d_iter)
|
||||
linf_r = np.zeros(d_iter)
|
||||
timer = np.zeros(d_iter)
|
||||
mcell = np.zeros(d_iter)
|
||||
#%% Loop through decreasing meshes and measure the residual
|
||||
# Create mesh using simpeg and write out in GIF format
|
||||
|
||||
for ii in range(d_iter):
|
||||
|
||||
|
||||
nc = int(sclx/dx)
|
||||
nc = int(sclx/dx[ii])
|
||||
|
||||
hxind = [(dx, 2*nc)]
|
||||
hyind = [(dx, nc)]
|
||||
hzind = [(dx, nc)]
|
||||
hxind = [(dx[ii], 2*nc)]
|
||||
hyind = [(dx[ii], nc)]
|
||||
hzind = [(dx[ii], nc)]
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
mcell[ii] = mesh.nC
|
||||
|
||||
actv = PF.Magnetics.getActiveTopo(mesh,topo,'N')
|
||||
|
||||
# Drape observations on topo + offset
|
||||
@@ -74,15 +84,13 @@ for ii in range(d_iter):
|
||||
xn = mesh.vectorNx
|
||||
yn = mesh.vectorNy
|
||||
zn = mesh.vectorNz
|
||||
|
||||
mcell = mesh.nC
|
||||
|
||||
print 'Mesh size: ' + str(mcell)
|
||||
|
||||
print 'Mesh size: ' + str(mcell[ii])
|
||||
|
||||
#%% Create model
|
||||
chibkg = 0.
|
||||
chiblk = 0.01
|
||||
model = np.ones(mcell)*chibkg
|
||||
model = np.ones(mcell[ii])*chibkg
|
||||
|
||||
# Do a three sphere problem for more frequencies
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, 0., 0., -sclx/3, R)
|
||||
@@ -98,7 +106,13 @@ for ii in range(d_iter):
|
||||
Utils.writeUBCTensorModel('Model.sus',mesh,model)
|
||||
actv = np.ones(mesh.nC)
|
||||
#%% Forward mode ldata
|
||||
|
||||
start_time = time.time()
|
||||
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
timer[ii] = (time.time() - start_time)
|
||||
|
||||
#fwr_tmi = d[0:ndata]
|
||||
#fwr_y = d[ndata:2*ndata]
|
||||
#fwr_z = d[2*ndata:]
|
||||
@@ -130,8 +144,9 @@ for ii in range(d_iter):
|
||||
#r_By = fwr_y - bya
|
||||
#r_Bz = fwr_z - bza
|
||||
|
||||
lrl[ii] = sum( r_tmi**2 ) **0.5
|
||||
|
||||
l2_r[ii] = np.sum( r_tmi**2 ) **0.5
|
||||
l1_r[ii] = np.sum( np.abs( r_tmi ) )
|
||||
linf_r[ii] = np.max( np.abs( r_tmi ) )
|
||||
|
||||
#%% Write predicted to file
|
||||
|
||||
@@ -139,34 +154,34 @@ PF.Magnetics.writeUBCobs('Obsloc.loc',B,M,rxLoc,d,np.ones(len(d)))
|
||||
|
||||
#%% Plot results
|
||||
print 'Residual between analytical sphere and integral forward'
|
||||
print "dx \t nc \t l1 \t l2 \t linf \t Runtime"
|
||||
for ii in range(d_iter):
|
||||
nc = 3**(ii+1)
|
||||
|
||||
print "||r||= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
|
||||
print str(dx[ii]) + "\t" + str(mcell[ii]) + "\t" + str(l1_r[ii]) + "\t" + str(l2_r[ii]) + "\t" + str(linf_r[ii]) + "\t" + str(timer[ii])
|
||||
|
||||
#%% Plot fields
|
||||
plt.figure(1)
|
||||
#ax = plt.subplot(221)
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(b_tmi,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.04)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(b_tmi,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(b_tmi,X.shape), s=20)
|
||||
ax.set_title('Analytical')
|
||||
|
||||
#%% Plot the forward solution from integral
|
||||
plt.figure(2)
|
||||
#ax = plt.subplot(222)
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
plt.colorbar(fraction=0.04)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
|
||||
ax.set_title('Numerical')
|
||||
|
||||
#%% Plot residual data
|
||||
plt.figure(3)
|
||||
#ax = plt.subplot(212)
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(r_tmi,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.04)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(r_tmi,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(r_tmi,X.shape), s=20)
|
||||
ax.set_title('Sphere Ana Bx')
|
||||
@@ -1,19 +1,18 @@
|
||||
import os
|
||||
|
||||
home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
#%%
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
import scipy.interpolate as interpolation
|
||||
import scipy.signal as sn
|
||||
import time
|
||||
|
||||
from interpFFT import interpFFT
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
|
||||
import os
|
||||
|
||||
home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
plt.close('all')
|
||||
|
||||
topofile = 'Gaussian.topo'
|
||||
@@ -78,7 +77,7 @@ mcell = mesh.nC
|
||||
print 'Mesh size: ' + str(mcell)
|
||||
|
||||
#%% Create model
|
||||
chibkg = 0.000
|
||||
chibkg = 0.0001
|
||||
chiblk = 0.01
|
||||
model = np.ones(mcell)*chibkg
|
||||
|
||||
@@ -92,9 +91,12 @@ model[sph_ind] = 4.*chiblk
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, sclx/2., 0., -sclx/2.5, R/2.5)
|
||||
model[sph_ind] = 2.5*chiblk
|
||||
|
||||
# Zero out
|
||||
model[actv==0] = -100
|
||||
Utils.writeUBCTensorMesh('Mesh.msh',mesh)
|
||||
Utils.writeUBCTensorModel('Model.sus',mesh,model)
|
||||
actv = np.ones(mesh.nC)
|
||||
Utils.writeUBCTensorModel('nullcell.dat',mesh,actv)
|
||||
#actv = np.ones(mesh.nC)
|
||||
#%% Forward mode ldata
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
#fwr_tmi = d[0:ndata]
|
||||
@@ -103,7 +105,7 @@ d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
|
||||
#%% Compute data on a line
|
||||
xx = np.linspace(-101., 99., 200)
|
||||
xx = np.linspace(xr.min(), xr.max(), 200)
|
||||
yy = np.zeros(len(xx))
|
||||
zz = F(xx,0.) + zoffset
|
||||
|
||||
@@ -118,20 +120,21 @@ l2_r = np.zeros((d_iter,5))
|
||||
linf_r = np.zeros((d_iter,5))
|
||||
timer = np.zeros((d_iter,5))
|
||||
|
||||
d2d = np.reshape(d, (len(yr),len(xr)))
|
||||
#%% Try different interpolation schemes
|
||||
for ii in range(d_iter):
|
||||
|
||||
indx = ii+1
|
||||
|
||||
d2d = np.reshape(d, (len(yr),len(xr)))
|
||||
d2d = mkvc(d2d[::indx,::indx].T)
|
||||
|
||||
dsub = d2d[::indx,::indx]
|
||||
|
||||
xsub = xr[::indx]
|
||||
ysub = yr[::indx]
|
||||
|
||||
# Nearest Neighbourg
|
||||
start_time = time.time()
|
||||
F = interpolation.NearestNDInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],d2d)
|
||||
F = interpolation.NearestNDInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T))
|
||||
d_i2d_nnb = mkvc( F(0.,xx) )
|
||||
l1_r[ii,0] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5
|
||||
l2_r[ii,0] = np.sum( (d_line - d_i2d_nnb)**2. )
|
||||
@@ -140,7 +143,7 @@ for ii in range(d_iter):
|
||||
|
||||
# Linear interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,d2d)
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T))
|
||||
d_i2d_lin = mkvc( F(0.,xx) )
|
||||
l1_r[ii,1] = np.sum( np.abs(d_line - d_i2d_lin) )**0.5
|
||||
l2_r[ii,1] = np.sum( (d_line - d_i2d_lin)**2. )
|
||||
@@ -150,7 +153,7 @@ for ii in range(d_iter):
|
||||
|
||||
# Cubic interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,d2d,kind='cubic')
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='cubic')
|
||||
d_i2d_cub = mkvc( F(0.,xx) )
|
||||
l1_r[ii,2] = np.sum( np.abs(d_line - d_i2d_cub) )**0.5
|
||||
l2_r[ii,2] = np.sum( (d_line - d_i2d_cub)**2. )
|
||||
@@ -160,95 +163,48 @@ for ii in range(d_iter):
|
||||
|
||||
# Quintic interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,d2d,kind='quintic')
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='quintic')
|
||||
d_i2d_qui = mkvc( F(0.,xx) )
|
||||
l1_r[ii,3] = np.sum( np.abs(d_line - d_i2d_qui) )**0.5
|
||||
l2_r[ii,3] = np.sum( (d_line - d_i2d_qui)**2. )
|
||||
linf_r[ii,3] = np.max( np.abs(d_line - d_i2d_qui) )
|
||||
timer[ii,3] = (time.time() - start_time)
|
||||
|
||||
# CloughTocher interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.CloughTocher2DInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T))
|
||||
d_i2d_CTI = mkvc( F(0.,xx) )
|
||||
l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_CTI) )**0.5
|
||||
l2_r[ii,4] = np.sum( (d_line - d_i2d_CTI)**2. )
|
||||
linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_CTI) )
|
||||
timer[ii,4] = (time.time() - start_time)
|
||||
|
||||
|
||||
|
||||
|
||||
#%% FFT interpolation
|
||||
|
||||
d2d = np.reshape(d, (len(yr),len(xr)))
|
||||
|
||||
# Add padding values by reflection (2**n)
|
||||
lenx = np.floor( np.log2( 2*len(xr) ) )
|
||||
npadx = int(np.floor( ( 2**lenx - len(xr) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npadx*2+len(xr),2) != 0:
|
||||
oddx = 1
|
||||
|
||||
else:
|
||||
oddx = 0
|
||||
|
||||
tap0 = sn.hamming(npadx*2)
|
||||
tapl = sp.spdiags(tap0[0:npadx],0,npadx,npadx)
|
||||
tapr = sp.spdiags(tap0[npadx:],0,npadx,npadx+oddx)
|
||||
|
||||
# Mirror the 2d data over the half lenght and apply 0-taper
|
||||
d2dpad = np.hstack([np.fliplr(d2d[:,0:npadx]) * tapl, d2d, np.fliplr(d2d[:,-npadx:]) * tapr])
|
||||
|
||||
# Repeat along the second dimension
|
||||
leny = np.floor( np.log2( 2*len(yr) ) )
|
||||
npady = int(np.floor( ( 2**leny - len(yr) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npady*2+len(yr),2) != 0:
|
||||
oddy = 1
|
||||
|
||||
else:
|
||||
oddy = 0
|
||||
|
||||
tap0 = sn.hamming(npady*2)
|
||||
tapu = sp.spdiags(tap0[0:npady],0,npady,npady)
|
||||
tapd = sp.spdiags(tap0[npady:],0,npady+oddy,npady)
|
||||
|
||||
d2dpad = np.vstack([tapu*np.flipud(d2dpad[0:npady,:]), d2dpad, tapd*np.flipud(d2dpad[-npady:,:])])
|
||||
|
||||
# Compute FFT
|
||||
FFTd2d = np.fft.fft2(d2dpad)
|
||||
|
||||
# Compute IFFT at a given location
|
||||
# Do an FFT shift
|
||||
#FFTshift = np.fft.fftshift(FFTd2d)
|
||||
FFTd2d = np.hstack([FFTd2d[:,0:31],np.zeros((32,64)),FFTd2d[:,31:]])
|
||||
FFTd2d = np.vstack([FFTd2d[0:15,:],np.zeros((32,128)),FFTd2d[15:,:]])
|
||||
# Pad with zeros
|
||||
#temp = np.zeros((FFTd2d.shape[0]*2,FFTd2d.shape[1]*2))
|
||||
|
||||
|
||||
# Compute inverse FFT
|
||||
IFFTd2d = np.fft.ifft2(FFTd2d)*FFTd2d.size/d2dpad.size
|
||||
|
||||
|
||||
|
||||
# Extract core
|
||||
#d2d_out = np.real(IFFTd2d[npady*2:-(npady*2+oddy+1),npadx*2:-(npadx*2+oddx+1)])
|
||||
d2d_out = np.real(IFFTd2d[:-1,:-1])
|
||||
d2d_out = d2d_out[npady*2:-(npady*2+oddy),npadx*2:-(npadx*2+oddx)]
|
||||
|
||||
# Create new distance vector
|
||||
XX = np.linspace(np.min(xr),np.max(xr),d2d_out.shape[1])
|
||||
YY = np.linspace(np.min(yr),np.max(yr),d2d_out.shape[0])
|
||||
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(XX,YY,d2d_out)
|
||||
d_i2d_fft = mkvc( F(xx,0.) )
|
||||
l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5
|
||||
l2_r[ii,4] = np.sum( (d_line - d_i2d_nnb)**2. )
|
||||
linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_nnb) )
|
||||
timer[ii,4] = (time.time() - start_time)
|
||||
|
||||
print("--- FFT completed in %s seconds ---" % (time.time() - start_time))
|
||||
|
||||
plt.figure()
|
||||
ax = plt.subplot()
|
||||
plt.imshow(d2d_out, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
plt.colorbar(fraction=0.04)
|
||||
#==============================================================================
|
||||
# #%% FFT interpolation
|
||||
# d2d_out = interpFFT(xsub,ysub,dsub)
|
||||
#
|
||||
# # Create new distance vector
|
||||
# XX = np.linspace(np.min(xsub),np.max(xsub),d2d_out.shape[1])
|
||||
# YY = np.linspace(np.min(ysub),np.max(ysub),d2d_out.shape[0])
|
||||
#
|
||||
# start_time = time.time()
|
||||
# F = interpolation.interp2d(XX,YY,d2d_out)
|
||||
# d_i2d_fft = mkvc( F(xx,0.) )
|
||||
# l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5
|
||||
# l2_r[ii,4] = np.sum( (d_line - d_i2d_nnb)**2. )
|
||||
# linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_nnb) )
|
||||
# timer[ii,4] = (time.time() - start_time)
|
||||
#
|
||||
# print("--- FFT completed in %s seconds ---" % (time.time() - start_time))
|
||||
#
|
||||
#
|
||||
# plt.figure()
|
||||
# ax = plt.subplot()
|
||||
# plt.imshow(d2d_out, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
# plt.colorbar(fraction=0.04)
|
||||
#==============================================================================
|
||||
|
||||
|
||||
#%% Write predicted to file
|
||||
@@ -287,8 +243,12 @@ plt.plot(xx,d_i2d_lin,c='b')
|
||||
plt.plot(xx,d_i2d_cub,c='g')
|
||||
plt.plot(xx,d_i2d_qui,c='m')
|
||||
plt.plot(xx,d_i2d_nnb,c='k')
|
||||
plt.plot(xx,d_i2d_fft,c='c',marker='o')
|
||||
plt.plot(xr[::indx],np.zeros(len(xr[::indx])),c='r',linewidth=0.,marker='o')
|
||||
plt.plot(xx,d_i2d_CTI,c='c')
|
||||
|
||||
# Plot interpolation from true value on line
|
||||
F = interpolation.interp1d(xx,d_line)
|
||||
dtrue = F(xr[::indx])
|
||||
plt.plot(xr[::indx],dtrue,c='r',linewidth=0.,marker='o')
|
||||
ax.set_title('Analytical')
|
||||
|
||||
#%% Write result to file
|
||||
|
||||
@@ -61,8 +61,10 @@ for ii in range(d_iter):
|
||||
model = np.ones(mcell)*chibkg
|
||||
model[sph_ind] = chiblk
|
||||
|
||||
actv = np.ones(mcell)
|
||||
|
||||
#%% Forward mode ldata
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,'xyz')
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'xyz')
|
||||
fwr_x = d[0:ndata]
|
||||
fwr_y = d[ndata:2*ndata]
|
||||
fwr_z = d[2*ndata:]
|
||||
|
||||
@@ -0,0 +1,85 @@
|
||||
def interpFFT(x,y,m):
|
||||
"""
|
||||
Load in a 2D grid and resample
|
||||
|
||||
OUTPUT:
|
||||
m_out
|
||||
|
||||
|
||||
"""
|
||||
|
||||
from SimPEG import np, sp
|
||||
import scipy.signal as sn
|
||||
|
||||
|
||||
# Add padding values by reflection (2**n)
|
||||
lenx = np.round( np.log2( 2*len(x) ) )
|
||||
npadx = int(np.floor( ( 2**lenx - len(x) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npadx*2+len(x),2) != 0:
|
||||
oddx = 1
|
||||
|
||||
else:
|
||||
oddx = 0
|
||||
|
||||
tap0 = sn.hamming(npadx*2)
|
||||
tapl = sp.spdiags(tap0[0:npadx],0,npadx,npadx)
|
||||
tapr = sp.spdiags(tap0[npadx:],0,npadx,npadx+oddx)
|
||||
|
||||
# Mirror the 2d data over the half lenght and apply 0-taper
|
||||
mpad = np.hstack([np.fliplr(m[:,0:npadx]) * tapl, m, np.fliplr(m[:,-npadx:]) * tapr])
|
||||
|
||||
# Repeat along the second dimension
|
||||
leny = np.round( np.log2( 2*len(y) ) )
|
||||
npady = int(np.floor( ( 2**leny - len(y) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npady*2+len(y),2) != 0:
|
||||
oddy = 1
|
||||
|
||||
else:
|
||||
oddy = 0
|
||||
|
||||
tap0 = sn.hamming(npady*2)
|
||||
tapu = sp.spdiags(tap0[0:npady],0,npady,npady)
|
||||
tapd = sp.spdiags(tap0[npady:],0,npady+oddy,npady)
|
||||
|
||||
mpad = np.vstack([tapu*np.flipud(mpad[0:npady,:]), mpad, tapd*np.flipud(mpad[-npady:,:])])
|
||||
|
||||
# Compute FFT
|
||||
FFTm = np.fft.fft2(mpad)
|
||||
|
||||
# Do an FFT shift
|
||||
FFTshift = np.fft.fftshift(FFTm)
|
||||
|
||||
# Pad high frequencies with zeros to increase the sampling rate
|
||||
py = int(FFTm.shape[0]/2)
|
||||
px = int(FFTm.shape[1]/2)
|
||||
|
||||
FFTshift = np.hstack([np.zeros((FFTshift.shape[0],px)),FFTshift,np.zeros((FFTshift.shape[0],px))])
|
||||
FFTshift = np.vstack([np.zeros((py,FFTshift.shape[1])),FFTshift,np.zeros((py,FFTshift.shape[1]))])
|
||||
|
||||
# Inverse shift
|
||||
FFTm = np.fft.ifftshift(FFTshift)
|
||||
|
||||
|
||||
# Compute inverse FFT
|
||||
IFFTm = np.fft.ifft2(FFTm)*FFTm.size/mpad.size
|
||||
|
||||
|
||||
m_out = np.real(IFFTm)
|
||||
# Extract core
|
||||
#m_out = np.real(IFFTm[npady*2:-(npady*2+oddy+1),npadx*2:-(npadx*2+oddx+1)])
|
||||
|
||||
|
||||
|
||||
m_out = m_out[npady*2:-(npady+oddy)*2,npadx*2:-(npadx+oddx)*2]
|
||||
|
||||
if np.mod(m.shape[0],2) != 0:
|
||||
m_out = m_out[:-1,:]
|
||||
|
||||
if np.mod(m.shape[1],2) != 0:
|
||||
m_out = m_out[:,:-1]
|
||||
|
||||
return m_out
|
||||
Reference in New Issue
Block a user