2017-05-31 46 views
0

我無法使用起始點和停止點的座標切片圖像。到目前爲止,我有以下代碼在兩個aribtray點上切片圖像

hdulist = fits.open(filename) 
hdr= hdulist[0].header 

import numpy as np 
import scipy 
import matplotlib.pyplot as plt 
from matplotlib.pyplot import figure, show 
from astropy.io import fits 
from scipy import interpolate 

data=hdulist[0].data 

#Make a line with "num" points 
D, B = input('Enter the coordinates of the starting point:').split(',') 
E, C = input("Enter the coordinates of the stopping point: ").split(',') 
x0= float(D) 
x1= float(E) 
y0= float(B) 
y1= float(C) 

x = np.arange(data.shape[1]) 
y = np.arange(data.shape[0]) 
#length = int((np.hypot(x1-x0, y1-y0))) (can be used instead of num_points) 
num_points = 1000 
xvalues = np.linspace(x0, x1, num_points) 
yvalues = np.linspace(y0, y1, num_points) 
f = scipy.interpolate.interp2d(x, y, data) #default is linear 


# Extract the values along the line 
profile = f(xvalues, yvalues) #this gives me a 2D array, I think it needs to be 1D 
#c = profile.flatten() 
print(profile.shape) 

'輪廓'不是線性的,而是立方體。有沒有辦法讓我的輪廓線性化,這樣我就可以在起點和終點之間的點上切分圖像?我只需要製作'輪廓'1D而不是2D。

我要繪製這樣的:

import numpy as np 
from numpy import random 
from matplotlib.pyplot import figure, show 

vels = np.linspace(0, 530, len(profile)) 
fig = figure() 
frame = fig.add_subplot(1,1,1) 
frame.plot(vels, profile) 
frame.set_ylabel('y-axis') 
frame.set_xlabel('x-axis') 
frame.grid(True) 
show() 
print(vels.shape) 
print(profile.shape) 
print(len(profile)) 

我的代碼不工作,因爲我得到的陰謀沒有顯示線路的片,但立方體的一個切片。

回答

0

interp2D的文檔看來,似乎網格是從插值構建的。就直覺而言,在我看來,你需要該網格的diagonal。建立快速的實驗與您的代碼的改編:

import numpy as np 
import scipy 
import matplotlib.pyplot as plt 
from matplotlib.pyplot import figure, show 
from scipy import interpolate 

X = np.arange(-5, 5, 0.25) 
Y = np.arange(-5, 5, 0.25) 
X, Y = np.meshgrid(X, Y) 
R = np.sqrt(X**2 + Y**2) 
data = np.sin(R) 

x0 = 2 
x1 = 23 
y0 = 1 
y1 = 36 

rx0 = 4 
rx1 = 2 
ry0 = 7 
ry1 = 32 

x = np.arange(data.shape[1]) 
y = np.arange(data.shape[0]) 
num_points = 1000 
xvalues = np.linspace(x0, x1, num_points) 
yvalues = np.linspace(y0, y1, num_points) 
f = scipy.interpolate.interp2d(x, y, data) #default is linear 


# Extract the values along the line 
profile = f(xvalues, yvalues) 

xvalues2 = np.linspace(rx0, rx1, num_points) 
yvalues2 = np.linspace(ry0, ry1, num_points) 
profile2 = f(xvalues2, yvalues2) 

plt.subplot(121) 
plt.imshow(data.T, origin="lower", interpolation="nearest") 
plt.scatter([x0, x1], [y0, y1]) 
plt.plot([x0, x1], [y0, y1]) 
plt.scatter([rx0, rx1], [ry0, ry1], c="r") 
plt.plot([rx0, rx1], [ry0, ry1], c="r") 
# plt.show() 

diag = np.diag(profile) 
diag2 = np.diag(profile2) 
plt.subplot(122) 
plt.plot(np.arange(diag.shape[0]), diag) 
plt.plot(np.arange(diag2.shape[0]), diag2, c="r") 
plt.show() 

這將返回以下:

Slicing a 2D image

注:我沒考慮到的座標到2D圖(這就是爲什麼這兩個線看起來相同的大小)。

+0

非常感謝,我的代碼現在按照我想要的方式工作! – Thomas