沿极坐标系的图像信息

Question

我有一组png图像，希望能用Python和相关工具进行处理。每个图像都代表一个已知尺寸的物体。

在每个图像中，都有物体的某个特定特征在一个特定的像素/物理位置，而每个图像的位置都不相同。

我想在给定的图像上施加一个极坐标系，原点位于这个特征的位置。

然后，我想获得以下信息： - 对于给定的极角，随径向位置变化的图像强度函数 - 当值在所有极角上平均时，随径向位置变化的图像强度函数。

我精通Python编程以及NumPy和SciPy中许多函数的使用，但是对于图像分析完全是新手。

如果您能就解决此问题采用可能的方法给我任何建议，我将不胜感激。

谢谢。

- cytochrome

---

根据你的描述，这似乎都是数组数学问题，而不是图像处理问题。你说你熟悉Python中的数组数学，那么你想要了解什么？也就是说，你所描述的并不是一个特定于图像处理的问题，比如你不是在尝试寻找对象、边缘、模糊或锐化等等。你想知道如何将图像导入Numpy中吗（这将是我解决这个问题的第一步）？ - tom10

Answer 1

您所描述的并不完全是传统意义上的图像处理，但使用numpy等工具很容易实现。

以下是一个相当大的示例，执行了您提到的一些操作，以指导您朝着正确的方向前进...请注意，示例图像都显示了原点在图像中心的结果，但函数接受一个原点参数，因此您可以直接为您的目的进行适应。

import numpy as np
import scipy as sp
import scipy.ndimage

import Image

import matplotlib.pyplot as plt

def main():
    im = Image.open('mri_demo.png')
    im = im.convert('RGB')
    data = np.array(im)

    plot_polar_image(data, origin=None)
    plot_directional_intensity(data, origin=None)

    plt.show()

def plot_directional_intensity(data, origin=None):
    """Makes a cicular histogram showing average intensity binned by direction
    from "origin" for each band in "data" (a 3D numpy array). "origin" defaults
    to the center of the image."""
    def intensity_rose(theta, band, color):
        theta, band = theta.flatten(), band.flatten()
        intensities, theta_bins = bin_by(band, theta)
        mean_intensity = map(np.mean, intensities)
        width = np.diff(theta_bins)[0]
        plt.bar(theta_bins, mean_intensity, width=width, color=color)
        plt.xlabel(color + ' Band')
        plt.yticks([])

    # Make cartesian coordinates for the pixel indicies
    # (The origin defaults to the center of the image)
    x, y = index_coords(data, origin)

    # Convert the pixel indices into polar coords.
    r, theta = cart2polar(x, y)

    # Unpack bands of the image
    red, green, blue = data.T

    # Plot...
    plt.figure()

    plt.subplot(2,2,1, projection='polar')
    intensity_rose(theta, red, 'Red')

    plt.subplot(2,2,2, projection='polar')
    intensity_rose(theta, green, 'Green')

    plt.subplot(2,1,2, projection='polar')
    intensity_rose(theta, blue, 'Blue')

    plt.suptitle('Average intensity as a function of direction')

def plot_polar_image(data, origin=None):
    """Plots an image reprojected into polar coordinages with the origin
    at "origin" (a tuple of (x0, y0), defaults to the center of the image)"""
    polar_grid, r, theta = reproject_image_into_polar(data, origin)
    plt.figure()
    plt.imshow(polar_grid, extent=(theta.min(), theta.max(), r.max(), r.min()))
    plt.axis('auto')
    plt.ylim(plt.ylim()[::-1])
    plt.xlabel('Theta Coordinate (radians)')
    plt.ylabel('R Coordinate (pixels)')
    plt.title('Image in Polar Coordinates')

def index_coords(data, origin=None):
    """Creates x & y coords for the indicies in a numpy array "data".
    "origin" defaults to the center of the image. Specify origin=(0,0)
    to set the origin to the lower left corner of the image."""
    ny, nx = data.shape[:2]
    if origin is None:
        origin_x, origin_y = nx // 2, ny // 2
    else:
        origin_x, origin_y = origin
    x, y = np.meshgrid(np.arange(nx), np.arange(ny))
    x -= origin_x
    y -= origin_y
    return x, y

def cart2polar(x, y):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    return r, theta

def polar2cart(r, theta):
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return x, y


def bin_by(x, y, nbins=30):
    """Bin x by y, given paired observations of x & y.
    Returns the binned "x" values and the left edges of the bins."""
    bins = np.linspace(y.min(), y.max(), nbins+1)
    # To avoid extra bin for the max value
    bins[-1] += 1 

    indicies = np.digitize(y, bins)

    output = []
    for i in xrange(1, len(bins)):
        output.append(x[indicies==i])

    # Just return the left edges of the bins
    bins = bins[:-1]

    return output, bins

def reproject_image_into_polar(data, origin=None):
    """Reprojects a 3D numpy array ("data") into a polar coordinate system.
    "origin" is a tuple of (x0, y0) and defaults to the center of the image."""
    ny, nx = data.shape[:2]
    if origin is None:
        origin = (nx//2, ny//2)

    # Determine that the min and max r and theta coords will be...
    x, y = index_coords(data, origin=origin)
    r, theta = cart2polar(x, y)

    # Make a regular (in polar space) grid based on the min and max r & theta
    r_i = np.linspace(r.min(), r.max(), nx)
    theta_i = np.linspace(theta.min(), theta.max(), ny)
    theta_grid, r_grid = np.meshgrid(theta_i, r_i)

    # Project the r and theta grid back into pixel coordinates
    xi, yi = polar2cart(r_grid, theta_grid)
    xi += origin[0] # We need to shift the origin back to 
    yi += origin[1] # back to the lower-left corner...
    xi, yi = xi.flatten(), yi.flatten()
    coords = np.vstack((xi, yi)) # (map_coordinates requires a 2xn array)

    # Reproject each band individually and the restack
    # (uses less memory than reprojection the 3-dimensional array in one step)
    bands = []
    for band in data.T:
        zi = sp.ndimage.map_coordinates(band, coords, order=1)
        bands.append(zi.reshape((nx, ny)))
    output = np.dstack(bands)
    return output, r_i, theta_i

if __name__ == '__main__':
    main()

原始图像：

投影到极坐标系：

从图像中心方向的强度函数：

- Joe Kington

---

@user458738，您应该通过点击其左侧的复选标记来接受他的答案。 - Justin Peel

不错的回答。不过，x轴标签似乎是反过来的。在极坐标图像中，颈部大约显示在pi/2弧度处，但在原始图像中则为-pi/2弧度。 - Matt Hancock

@JustinPeel 我知道这个答案有点老了.. 但我有几个问题.. 如何从转换后的图像（如第二张图片）计算出极坐标到笛卡尔坐标的反投影，以及 ＃基于最小和最大r＆theta创建规则（在极坐标空间中）网格的用途是什么 - Minato

非常感谢Joe Kinghorn提供的出色答案，节省了我很多时间。如果有人发现有用的话，这个reproject_image_into_polar函数的简单“灰度”版本现在已经并入了PyAbel包中：https://github.com/PyAbel/PyAbel/blob/master/abel/tools/polar.py - DanHickstein

Answer 2

以下是我使用Scipy的 geometric_transform 方法的操作：

topolar.py

import numpy as np
from scipy.ndimage.interpolation import geometric_transform

def topolar(img, order=1):
    """
    Transform img to its polar coordinate representation.

    order: int, default 1
        Specify the spline interpolation order. 
        High orders may be slow for large images.
    """
    # max_radius is the length of the diagonal 
    # from a corner to the mid-point of img.
    max_radius = 0.5*np.linalg.norm( img.shape )

    def transform(coords):
        # Put coord[1] in the interval, [-pi, pi]
        theta = 2*np.pi*coords[1] / (img.shape[1] - 1.)

        # Then map it to the interval [0, max_radius].
        #radius = float(img.shape[0]-coords[0]) / img.shape[0] * max_radius
        radius = max_radius * coords[0] / img.shape[0]

        i = 0.5*img.shape[0] - radius*np.sin(theta)
        j = radius*np.cos(theta) + 0.5*img.shape[1]
        return i,j

    polar = geometric_transform(img, transform, order=order)

    rads = max_radius * np.linspace(0,1,img.shape[0])
    angs = np.linspace(0, 2*np.pi, img.shape[1])

    return polar, (rads, angs)

以下是一些测试用法：

testpolar.py

from topolar import topolar
from skimage.data import chelsea

import matplotlib.pyplot as plt

img = chelsea()[...,0] / 255.
pol, (rads,angs) = topolar(img)

fig,ax = plt.subplots(2,1,figsize=(6,8))

ax[0].imshow(img, cmap=plt.cm.gray, interpolation='bicubic')

ax[1].imshow(pol, cmap=plt.cm.gray, interpolation='bicubic')

ax[1].set_ylabel("Radius in pixels")
ax[1].set_yticks(range(0, img.shape[0]+1, 50))
ax[1].set_yticklabels(rads[::50].round().astype(int))

ax[1].set_xlabel("Angle in degrees")
ax[1].set_xticks(range(0, img.shape[1]+1, 50))
ax[1].set_xticklabels((angs[::50]*180/3.14159).round().astype(int))

plt.show()

...并输出：

- Matt Hancock