Репроектирование полярной с декартовой сеткой

Ответ 1

Спасибо за ваши ответы - подумав немного больше об этом, я придумал следующий код:

import numpy as np

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as mpl

from scipy.interpolate import interp1d
from scipy.ndimage import map_coordinates


def polar2cartesian(r, t, grid, x, y, order=3):

    X, Y = np.meshgrid(x, y)

    new_r = np.sqrt(X*X+Y*Y)
    new_t = np.arctan2(X, Y)

    ir = interp1d(r, np.arange(len(r)), bounds_error=False)
    it = interp1d(t, np.arange(len(t)))

    new_ir = ir(new_r.ravel())
    new_it = it(new_t.ravel())

    new_ir[new_r.ravel() > r.max()] = len(r)-1
    new_ir[new_r.ravel() < r.min()] = 0

    return map_coordinates(grid, np.array([new_ir, new_it]),
                            order=order).reshape(new_r.shape)

# Define original polar grid

nr = 10
nt = 10

r = np.linspace(1, 100, nr)
t = np.linspace(0., np.pi, nt)
z = np.random.random((nr, nt))

# Define new cartesian grid

nx = 100
ny = 200

x = np.linspace(0., 100., nx)
y = np.linspace(-100., 100., ny)

# Interpolate polar grid to cartesian grid (nearest neighbor)

fig = mpl.figure()
ax = fig.add_subplot(111)
ax.imshow(polar2cartesian(r, t, z, x, y, order=0), interpolation='nearest')
fig.savefig('test1.png')

# Interpolate polar grid to cartesian grid (cubic spline)

fig = mpl.figure()
ax = fig.add_subplot(111)
ax.imshow(polar2cartesian(r, t, z, x, y, order=3), interpolation='nearest')
fig.savefig('test2.png')

Что не является строго повторной сеткой, но отлично работает для того, что мне нужно. Просто отправьте код, если он будет полезен кому-либо еще. Не стесняйтесь предлагать улучшения!

Ответ 2

Вы можете сделать это более компактно с помощью scipy.ndimage.geometric_transform. Вот пример кода:

import numpy as N
import scipy as S
import scipy.ndimage

temperature = <whatever> 
# This is the data in your polar grid.
# The 0th and 1st axes correspond to r and θ, respectively.
# For the sake of simplicity, θ goes from 0 to 2π, 
# and r units are just its indices.

def polar2cartesian(outcoords, inputshape, origin):
    """Coordinate transform for converting a polar array to Cartesian coordinates. 
    inputshape is a tuple containing the shape of the polar array. origin is a
    tuple containing the x and y indices of where the origin should be in the
    output array."""

    xindex, yindex = outcoords
    x0, y0 = origin
    x = xindex - x0
    y = yindex - y0

    r = N.sqrt(x**2 + y**2)
    theta = N.arctan2(y, x)
    theta_index = N.round((theta + N.pi) * inputshape[1] / (2 * N.pi))

    return (r,theta_index)

temperature_cartesian = S.ndimage.geometric_transform(temperature, polar2cartesian, 
    order=0,
    output_shape = (temperature.shape[0] * 2, temperature.shape[0] * 2),
    extra_keywords = {'inputshape':temperature.shape,
        'center':(temperature.shape[0], temperature.shape[0])})

Вы можете изменить order=0 по желанию для лучшей интерполяции. Выходной массив temperature_cartesian здесь равен 2r на 2r, но вы можете указать любой размер и начало, которое вам нравится.

Ответ 3

Я пришел к этому сообщению некоторое время назад, когда пытался сделать что-то подобное, это перепрофилирование полярных данных в декартовую сетку и наоборот. Предлагаемое здесь решение отлично работает. Однако для преобразования координат требуется некоторое время. Я просто хотел поделиться другим подходом, который может сократить время обработки до 50 раз и более.

В алгоритме используется функция scipy.ndimage.interpolation.map_coordinates.

Посмотрим на небольшой пример:

import numpy as np

# Auxiliary function to map polar data to a cartesian plane
def polar_to_cart(polar_data, theta_step, range_step, x, y, order=3):

    from scipy.ndimage.interpolation import map_coordinates as mp

    # "x" and "y" are numpy arrays with the desired cartesian coordinates
    # we make a meshgrid with them
    X, Y = np.meshgrid(x, y)

    # Now that we have the X and Y coordinates of each point in the output plane
    # we can calculate their corresponding theta and range
    Tc = np.degrees(np.arctan2(Y, X)).ravel()
    Rc = (np.sqrt(X**2 + Y**2)).ravel()

    # Negative angles are corrected
    Tc[Tc < 0] = 360 + Tc[Tc < 0]

    # Using the known theta and range steps, the coordinates are mapped to
    # those of the data grid
    Tc = Tc / theta_step
    Rc = Rc / range_step

    # An array of polar coordinates is created stacking the previous arrays
    coords = np.vstack((Ac, Sc))

    # To avoid holes in the 360º - 0º boundary, the last column of the data
    # copied in the begining
    polar_data = np.vstack((polar_data, polar_data[-1,:]))

    # The data is mapped to the new coordinates
    # Values outside range are substituted with nans
    cart_data = mp(polar_data, coords, order=order, mode='constant', cval=np.nan)

    # The data is reshaped and returned
    return(cart_data.reshape(len(y), len(x)).T)

polar_data = ... # Here a 2D array of data is assumed, with shape thetas x ranges

# We create the x and y axes of the output cartesian data
x = y = np.arange(-100000, 100000, 1000)

# We call the mapping function assuming 1 degree of theta step and 500 meters of
# range step. The default order of 3 is used.
cart_data = polar_to_cart(polar_data, 1, 500, x, y)

Я надеюсь, что это поможет кому-то в той же ситуации, что и я.

Ответ 4

Есть ли пакеты Python, которые могут это сделать?

Да! В настоящее время существует, по крайней мере, один пакет Python, который имеет функцию для перепрограммирования матрицы из декартовых координат в полярные: abel.tools.polar.reproject_image_into_polar(), который является частью пакета PyAbel.

(Иниго Эрнаес Коррес прав, scipy.ndimage.interpolation.map_coordinates - это самый быстрый способ, который мы до сих пор использовали для перепрограммирования из декартовых координат в полярные).

PyAbel можно установить из PyPi, введя в командной строке следующее:

pip install pyabel

Затем в python вы можете использовать следующий код для повторного проецирования изображения в полярные координаты:

import abel
abel.tools.polar.reproject_image_into_polar(MyImage)

[В зависимости от приложения вы можете рассмотреть передачу аргумента jacobian=True, который перечерчивает интенсивность матрицы, чтобы принять во внимание растяжку сетки (меняя "размеры бункера" ), которая имеет место, когда вы превращаются из декартовых в полярные.)

Вот полный пример:

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

CartImage = abel.tools.analytical.sample_image(501)[201:-200, 201:-200]

PolarImage, r_grid, theta_grid = abel.tools.polar.reproject_image_into_polar(CartImage)

fig, axs = plt.subplots(1,2, figsize=(7,3.5))
axs[0].imshow(CartImage , aspect='auto', origin='lower')
axs[1].imshow(PolarImage, aspect='auto', origin='lower', 
              extent=(np.min(theta_grid), np.max(theta_grid), np.min(r_grid), np.max(r_grid)))

axs[0].set_title('Cartesian')
axs[0].set_xlabel('x')
axs[0].set_ylabel('y')

axs[1].set_title('Polar')
axs[1].set_xlabel('Theta')
axs[1].set_ylabel('r')

plt.tight_layout()
plt.show()

Примечание: есть еще одно хорошее обсуждение (о повторном отображении цветных изображений в полярные координаты) на SO: информация о изображении вдоль полярной системы координат