python - 在一个图中同时绘制两个横截面强度

我有一个形状数组（512,512）。
看起来像（（行= x，列= y，密度= z =数组的数量）

[[0.012825 0.020408 0.022976 ... 0.015938 0.02165  0.024357]
 [0.036332 0.031904 0.025462 ... 0.031095 0.019812 0.024523]
 [0.015831 0.027392 0.031939 ... 0.016249 0.01697  0.028686]
 ...
 [0.024545 0.011895 0.022235 ... 0.033226 0.03223  0.030235]]

我已经将其绘制为2D密度图。我的目标是找到圆的中心，并在一个图中绘制垂直和水平的横截面。

现在，我很难找到圆的中心并将两个横截面合并为一个图形。

请帮忙。

这是我的代码：

    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    import scipy.ndimage
    data = pd.read_csv('D:/BFP.csv', header=None)


    # create data
    data = np.array(data)
    print(data)

    #plot data
    side = np.linspace(-1.5,1.5,512)
    x,y = np.meshgrid(side,side)
    z = [[data[i][j] for i in range(len(data[0]))]for j in range(len(data))]


    #-- Extract the line...
    # Make a line with "num" points...
    x0, y0 = 270, 0  # These are in _pixel_ coordinates!!
    x1, y1 = 270, 500
    num = 512
    x_, y_ = np.linspace(x0, x1, num), np.linspace(y0, y1, num)

    # Extract the values along the line, using cubic interpolation
    zi = scipy.ndimage.map_coordinates(z, np.vstack((x_,y_)))


    #-- Plot...
    fig, axes = plt.subplots(nrows=2)
    axes[0].imshow(z,origin='lower')
    axes[0].plot([x0, x1], [y0, y1], 'ro-')
    #axes[0].axis('image')

    axes[1].plot(zi)
    plt.savefig('D:/vertical.png')
    plt.show()

图片在这里：

最佳答案

我无法帮助您找到圆心，但是您可以通过在网格中创建3个轴来创建横截面的漂亮可视化。通常，我会为此使用GridSpec，但是imhsow倾向于弄乱轴的相对大小以保持正方形像素。值得庆幸的是，AxesGrid toolkit可以提供帮助。

该代码的基础受this matplotlib example的启发。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.stats import multivariate_normal
import scipy

fig, main_ax = plt.subplots(figsize=(5, 5))
divider = make_axes_locatable(main_ax)
top_ax = divider.append_axes("top", 1.05, pad=0.1, sharex=main_ax)
right_ax = divider.append_axes("right", 1.05, pad=0.1, sharey=main_ax)

# make some labels invisible
top_ax.xaxis.set_tick_params(labelbottom=False)
right_ax.yaxis.set_tick_params(labelleft=False)

main_ax.set_xlabel('dim 1')
main_ax.set_ylabel('dim 2')
top_ax.set_ylabel('Z profile')
right_ax.set_xlabel('Z profile')

x, y = np.mgrid[-1:1:.01, -1:1:.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
rv = multivariate_normal([-0.2, 0.2], [[1, 1.5], [0.25, 0.25]])
z = rv.pdf(pos)
z_max = z.max()

cur_x = 110
cur_y = 40

main_ax.imshow(z, origin='lower')
main_ax.autoscale(enable=False)
right_ax.autoscale(enable=False)
top_ax.autoscale(enable=False)
right_ax.set_xlim(right=z_max)
top_ax.set_ylim(top=z_max)
v_line = main_ax.axvline(cur_x, color='r')
h_line = main_ax.axhline(cur_y, color='g')
v_prof, = right_ax.plot(z[:,int(cur_x)],np.arange(x.shape[1]), 'r-')
h_prof, = top_ax.plot(np.arange(x.shape[0]),z[int(cur_y),:], 'g-')

plt.show()

只是为了好玩，您甚至可以使其互动

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.stats import multivariate_normal
import scipy

fig, main_ax = plt.subplots(figsize=(5, 5))
divider = make_axes_locatable(main_ax)
top_ax = divider.append_axes("top", 1.05, pad=0.1, sharex=main_ax)
right_ax = divider.append_axes("right", 1.05, pad=0.1, sharey=main_ax)

# make some labels invisible
top_ax.xaxis.set_tick_params(labelbottom=False)
right_ax.yaxis.set_tick_params(labelleft=False)

main_ax.set_xlabel('dim 1')
main_ax.set_ylabel('dim 2')
top_ax.set_ylabel('Z profile')
right_ax.set_xlabel('Z profile')

x, y = np.mgrid[-1:1:.01, -1:1:.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
rv = multivariate_normal([-0.2, 0.2], [[1, 1.5], [0.25, 0.25]])
z = rv.pdf(pos)
z_max = z.max()

main_ax.imshow(z, origin='lower')
main_ax.autoscale(enable=False)
right_ax.autoscale(enable=False)
top_ax.autoscale(enable=False)
right_ax.set_xlim(right=z_max)
top_ax.set_ylim(top=z_max)
v_line = main_ax.axvline(np.nan, color='r')
h_line = main_ax.axhline(np.nan, color='g')
v_prof, = right_ax.plot(np.zeros(x.shape[1]),np.arange(x.shape[1]), 'r-')
h_prof, = top_ax.plot(np.arange(x.shape[0]),np.zeros(x.shape[0]), 'g-')

def on_move(event):
    if event.inaxes is main_ax:
        cur_x = event.xdata
        cur_y = event.ydata

        v_line.set_xdata([cur_x,cur_x])
        h_line.set_ydata([cur_y,cur_y])
        v_prof.set_xdata(z[:,int(cur_x)])
        h_prof.set_ydata(z[int(cur_y),:])

        fig.canvas.draw_idle()

fig.canvas.mpl_connect('motion_notify_event', on_move)

plt.show()

注意：滞后只是由于gif转换，更新在我的机器上更加流畅