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

Question

我有一个形状数组（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]]

我已经将它绘制成二维密度图。我的目标是找到圆心并在一个图中画出一个垂直和水平的横截面。

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

请帮忙。

这是我的代码：

    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()

图片在这里：

Answer 1

我无法帮助您找到圆心，但您可以通过在网格中创建 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，在我的机器上更新更流畅