1
$\begingroup$

enter image description here

At the Top "right", there is the 2D-density plot of the recorded data (actual), fewer in number. Recorded data has been sampled a on the 8 arms of a regular octagon. These 8 arms are placed along x and y axes and along the diagonals. On the top-left is the interpolated data (thin-plate) of the actual data, i.e. larger in number. Compared to the Right (actual) data-set we can observe that "Interpolated" data is showing some kind of "Quadrupole" kind of (which also looks like wings of a "Butterfly") artifact which is not desired and expected. I was expecting kind circular kind of "Interpolated" density plot as can be shown at the feet of the bottom figures.

I have tried both my own written-code and in-built python module to perform "Thin-Plate-spline" interpolation and both the results are showing same "Quadrupole" kind of artifact in the interpolated data (see the fig. at top-left corner). Codes are attached below with both in-bult scipy.interpolate() module and my own code as well. In my own code I have solved the equation as described below,

https://elonen.iki.fi/code/tpsdemo/

My question is why is "Quadrupole" or "Butterfly" sort of artifact is arising instead of a rather circular or octagonal pattren (which is seen in actual data, Top-Right corner Figure) in "Thin Plate Spline" interpolation ?

How shall I be able to reduce that ? For example,The top-left figure is the actual data and top-right one is the corresponding "Interpolated" data which is different kind of pattern, right ?. This implies that I might be making some errors in the process of interpolation, but where am I committing the mistake is not clear to me.

Now,I have tried 3 different codes to have the interpolated data. And, in each case I found the artifact to be present there in the interpolated data.

2 of these 3 codes have used "in-built"

scipy.interpolate()-thin_plate spline whereas I also ran my own code of "Thin-Plate-Spline". I was hoping if there is any mistake in my codes and thus I posted my codes along with the question. I was hoping that it might aid reviewer to understand where does the onset of "Quadrupole" or "Butterfly" kind of "artifact" lie. I would like to request if any one could kindly response and guide me about what has to be done. Kindly, please have your indispensable expertise.

CODE - 1 => Density-Plot

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100 
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
plt.plot(x, y, 'ko')
plt.imshow(Ei, extent=[0, 1, 1, 0], cmap='gist_earth')
plt.colorbar()
###########################################################
plt.show()

CODE - 2 => 3D-Surface Plot

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

#############################################################
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5
#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100 
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
#p = axes.plot_surface(xi,yi,Ei, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
p = axes.plot_surface(xi, yi, Ei,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
###########################################################
plt.show()

CODE - 3 => My own code - 3D

import numpy as np
import math as mt


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5

#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis


#########################################################################################################
print(" ")
################################################################

#data = np.loadtxt('data1.dat',unpack=True)
#x = data[0,:] 
#y = data[1,:]
#E1 = data[2,:]

################################################################

data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,16:])/100 
y = (data[1,16:])/100
E1 = data1[0,16:]

################################################################
################################################################
n1 = len(x)
#print(" n1 : ", n1)
r = np.stack((x, y), axis=-1)
print(r)


N = n1 + 3

o = np.zeros(3)


#################################################################

#################################################################

V1 = np.hstack((E1, o))
V = V1.reshape(N, 1)

##########################################################################################################


C = np.ones(n1)

P = np.stack((C, x, y), axis=-1)
P_T = P.transpose()  # Transpose of P

#print(" P : ")
#print(P)

#print(" P_T : ")
#print(P_T)


#########################################################################################################

O1 = np.zeros(9)
O = O1.reshape((3,3))

#print(" O : ")
#print(O)

#########################################################################################################

K = np.zeros((n1,n1))

for i in range(n1):
    for j in range(n1):

        if(i == j):
            K[i][j] = 0.0
        else:
            R = r[i] - r[j]
            s = np.sqrt( (R[1]*R[1]) + (R[0]*R[0]) )
#           print(s)
            K[i][j] =  (s*s) * (np.log(s)) 

#print(" U[i][j]")
#print(K)

##########################################################################################################

#L1 = np.block([K, P])
L1 = np.hstack((K, P))
#L2 = np.block([P_T, O])
L2 = np.hstack((P_T, O))
#L = np.block([[L1],[L2]])
L = np.vstack((L1,L2))
#print(" L : ")
#print(L)

LI = np.linalg.inv(L)

###########################################################################################################

W = np.dot(LI, V)   # Weights

print(" ")
#print(" W : ")
#print(W)

###########################################################################################################

def f(x, y, r, N, W):

    f = W[N-3] + (W[N-2])*x + (W[N-1])*y
 
    for i in range(N-3):

        r1 = np.sqrt(((r[i][0] - x)**2) + ((r[i][1] - y)**2))
#       r1 = (((r[i][0] - x)**2) + ((r[i][1] - y)**2))
        if(r1==0.0):
            U1 = 50     #IMPORATANT
        else:
            U1 = (r1*r1) * (np.log(r1))

                
        f = f + W[i]*U1
    
    if(f < 0.0):
        f = 0.0

    return f


#f9 = f(35, 35, r, N, W)  # value of f9 is in good approximation with the value achieved from simulation
#print(f9)

## upto this point CODE is working well


print("_______________________________________________________________________________________")
print(" ")


#F = f(xv, yv, x, y, N, W) # This is not working 



x1 = np.linspace(-200, 200, 1000)
y1 = np.linspace(-200, 200, 1000)
xv, yv = np.meshgrid(x1,y1)

f_xy = np.zeros((1000, 1000))

for i in range(1000):
    for j in range(1000):
        f_xy[i][j] = f(x1[i], y1[j], r, N, W)


#############################################################################################################################

p = axes.plot_surface(xv,yv,f_xy, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
##p = axes.plot_surface(xv, yv, f_xy,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
plt.show()
$\endgroup$
6
  • $\begingroup$ What exactly is the question? You show three pieces of code, but it's not clear what we're supposed to see from them. The post could do with a more concise question. $\endgroup$ May 14 at 3:38
  • $\begingroup$ @WolfgangBangerth , Sir my question is essentially why is "Butterfly" sort of artifact arising in "Thin Plate Spline" ? And, how shall I be able to reduce that ? For example,The top-left figure is the actual data and top-right one is the corresponding "Interpolated" data which is possessing the "Artifact". Now, I have tried 3 different codes to have the interpolated data. And, in each case I found the artifact to be present there in the interpolated data. 2 of these 3 codes have used "in-built" scipy.interpolate()-thin_plate spline whereas I also ran my own code of "Thin-Plate-Spline". $\endgroup$ May 14 at 5:56
  • $\begingroup$ @WolfgangBangerth , I was wondering if there is any mistake in my codes and thus I posted my codes along with the question. I was hoping that it might aid reviewer to understand where does the onset of "Butterfly" kind of "artifact" lie. Thanks a lot. Sir, Is it clear now ? Kindly let me know $\endgroup$ May 14 at 5:58
  • $\begingroup$ Please edit your question to make it clearer. Incorporate your comments here into your edited question. While you're there, you may want to define what "butterfly artifact" means -- it's not a term I am familiar with, and I suspect many others aren't either. $\endgroup$ May 16 at 21:28
  • $\begingroup$ @WolfgangBangerth , Sir, I have edited the question. Kindly let me know if it is not in proper order. Thanks a lot. $\endgroup$ yesterday

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.