如何根据 contourf 或 pcolormesh 图上的旋转数据绘制特定轮廓

Question

我有一个存储在 numpy 数组中的图像。我创建了一个函数来将数据旋转一个角度 theta。为了执行旋转，该函数将图像 (i,j) 的索引坐标转换为 (x,y) 并应用旋转矩阵。然后该函数返回旋转后的 (X, Y) 坐标的网格。

我想在同一坐标系上叠加非旋转图像和旋转图像，并提取特定的垂直和水平剖面。我无法正确导航旋转的图像，因为它只能使用 map_coordinates 函数（据我所知）使用“ij”进行导航。

设置和功能定义：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pyplot as plt

def rotate_image(arr, dpi, theta_degrees = 0.0, pivot_point = [0,0]):

  theta_radians = (np.pi/180.0)* theta_degrees
  c = round(np.cos(theta_radians), 3)
  s = round(np.sin(theta_radians), 3)

  rotation_matrix = np.array([[c, -s, 0],
                              [s, c, 0],
                              [0, 0,  1]])
  #print(rotation_matrix)

  width, height = arr.shape
  pivot_point_xy = np.array([(25.4 / dpi[0])* pivot_point[0], (25.4/dpi[1])*pivot_point[1]])
  pivot_shift_vector = np.array([[pivot_point_xy[0]],
                                 [pivot_point_xy[1]],
                                 [0]])
  
  x = (25.4 / dpi[0]) * np.array(range(width)) #convert pixels to mm units
  y = (25.4 / dpi[1]) * np.array(range(height))#convert pixels to mm units
  
  XX , YY = np.meshgrid(x,y)
  ZZ = arr
  coordinates = np.stack([XX,YY,ZZ])
  #shift to rotation point, apply rotation, shift back to original coordinates
  coordinates_reshape = np.reshape(coordinates, (3,-1))
  translated_coordinates = coordinates_reshape - pivot_shift_vector
  rotated_coordinates = np.matmul(rotation_matrix, translated_coordinates)
  final_coordinates = rotated_coordinates + pivot_shift_vector
  final_coordinates_reshaped = np.reshape(final_coordinates, (3, width, height))
  
  return final_coordinates_reshaped

示例图：

img = np.arange(1,26).reshape((5,5))

rotated_img_0 = rotate_image(img, theta_degrees= 0, dpi =[1,1], pivot_point = [2.5,2.5])
rotated_img_1 = rotate_image(img, theta_degrees= 45, dpi =[1,1], pivot_point = [2.5,2.5])

# plot
fig, ax = plt.subplots(2, 1, figsize = (10,20))

ax[0].pcolormesh(*rotated_img_0, vmin=0, vmax=rotated_img_0[2].max())
ax[0].pcolormesh(*rotated_img_1, vmin=0, vmax=rotated_img_1[2].max(), alpha = 0.7)
ax[0].hlines(60, rotated_img_1[0].min(), rotated_img_1[0].max() , color = 'black')

ax[1].contourf(*rotated_img_0, vmin=0, vmax=rotated_img_0[2].max())
ax[1].contourf(*rotated_img_1, vmin=0, vmax=rotated_img_1[2].max(), alpha = 0.7)
ax[1].hlines(60, rotated_img_1[0].min(), rotated_img_1[0].max() , color = 'black')

plt.show()

我试图从 scipy 改编这里概述的 interpolate2d 方法，但它不适用于旋转数据：https://docs.scipy.org/doc//scipy-0.17.0/reference/generated/scipy.interpolate.interp2d.html

Map_coordinates 也适用于使用“ij”坐标的非旋转数据。 i,j 的简单切片也可以满足我的目的。

我希望能够在相同的 xy 坐标下从每个图表中提取相同的配置文件。

Answer 1

虽然这不是一个直接的答案，但我决定最好解决这个问题。我重写了 rotate_image 函数，以便可以使用简单的数组切片来使用 map_coordinates 函数提取配置文件。

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage

def rotate_image(img, theta_degrees = 0.0, pivot_point = np.array([[0], 
[0]])):

    width , height = img.shape

    #compute the 2D rotation matrix
    theta_radians = (np.pi/180.0)* theta_degrees
    c = round(np.cos(theta_radians), 3)
    s = round(np.sin(theta_radians), 3)

    rotation_matrix = np.array([[c, -s],
                          [s, c]])

    #create a sampling point cloud using meshgrid
    X, Y = range(width), range(height)
    XX , YY = np.meshgrid(X, Y)
    coordinates = np.stack([XX, YY])
    coordinates_reshape = np.reshape(coordinates, (2,-1))

    #rotate the image around the chosen pivot point
    translated_coordinates = coordinates_reshape - pivot_point
    rotated_coordinates = np.matmul(rotation_matrix, 
      translated_coordinates)
    final_coordinates = rotated_coordinates + pivot_point
    final_coordinates_reshaped = np.reshape(final_coordinates, (2, width, 
              height))

    #use scipy map_coordinates function to resample the image at the new 
    #coordinates
    rotated_image = ndimage.map_coordinates(img, 
                     final_coordinates_reshaped, mode = 'constant')

    return rotated_image

img = np.random.rand(100, 100)*100

rotated_img = rotate_image(img, 20, np.array([[50],[50]]))

#graph Results
fig, ax = plt.subplots(2,1)
ax[0].imshow(rotated_img)
ax[0].hlines(50, 0, 99)
ax[1].plot(rotated_img[50])