图像极坐标变换及超声图像上的应用-EW帮帮网

为什么需要极坐标变换？

在图像处理和计算机视觉中，我们最常见的坐标体系是笛卡尔坐标系（Cartesian Coordinates），它用水平和垂直方向的像素位置来描述图像。然而，在某些成像模式中，原始数据并不是以直角坐标采样的，而是以极坐标系（Polar Coordinates）方式获取的。医学超声成像就是一个典型例子：探头发出的超声波沿着扇形范围逐线扫描，每条扫描线记录从探头中心出发、沿特定角度方向的回波信息。这种采样方式所"拍摄"的图像一般是一个扇形区域，而我们采集的数据构成一个了一个矩阵，因此我们需要将采样后的数据进行合适的变换以呈现出"拍摄"图像的真实形状。

极坐标与笛卡尔坐标的关系

如上图所示，极坐标与笛卡尔坐标共享原点。任意点的笛卡尔坐标为该点到x轴与y轴的方向向量的投影决定。任意点的极坐标表示了该点到原点的距离以及该点从x轴上绕原点逆时针旋转一定角度。

已知极坐标 $(r,\theta)$ ，根据几何关系其对应的笛卡尔为 $\left\{\begin{matrix} x=rcos\theta\\ y=rsin\theta \end{matrix}\right.$ 。已知笛卡尔坐标 $(x,y)$ ，根据几何关系其对应的极坐标 $\left\{\begin{matrix} r=\sqrt{x^{2}+y^{2}}\\ \theta = atan2(y,x) \end{matrix}\right.$ 。

超声图像扫描的基本原理

医学超声成像是利用高频声波（通常 2–20 MHz）在组织中的传播特性来获取图像的技术。探头内的压电晶体既可以发射超声脉冲，又能接收从组织界面反射回来的回波信号。常见的超声扫描模式有：

线性扫描（Linear Scan）：超声波束沿平行方向发射，形成矩形视野（常用于浅表组织）。
扇形扫描（Sector Scan）：以探头表面一点为中心，超声波束沿不同角度方向扫描，形成扇形视野。我们所处理的超声图像即为扇形扫描方式获取。

在扇形扫描中，探头发射的每一条扫描线（scan line）对应一个固定的角度 $\theta$ ，沿这条方向上会采样多个深度点（depth samples），记录回波强度。因此，原始超声数据其实存储在极坐标系下：

半径 $r$ ：回波的深度（采样点序号 × 声速/采样频率）
角度 $\theta$ ：扫描线的角度位置

如果探头有 $N$ 条扫描线，每条线有 $M$ 个采样点，那么原始数据就是一个 $M\times N$ 的矩阵：

列方向对应扫描角度变化
行方向对应深度变化

为了将 $M\times N$ 扫描线矩阵转换为真实的超声图像(即扇形图像)，可以采样极坐标进行映射。

确定扫描几何参数，包括：总扫描角度，单个扫描线对应的扫描角度，扫描线的扫描深度。
这里我们的总扫描角度为 $\theta$ ，扫描线为 $N$ 条，则单个扫描线对应的扫描角度为 $\frac{\theta}{N}$ ，同时每条扫描线的扫描深度为 $M$ 。这就构成了基本的扫描信息。
极坐标到笛卡尔坐标转换
根据基本扫描几何参数，任意扫描矩阵上的点 $(m,n)$ 对应一个极坐标点 $(m, \frac{n \theta}{N})$ 。然后将其转换为笛卡尔坐标为 $(m*cos(\frac{n\theta}{N}),m*sin(\frac{n\theta}{N}))$ 。
绘制扇形图像
对于扫描线上的每一个采样点，我们总是可以转换到笛卡尔坐标系下。将所有的扫描线点均转换到笛卡尔坐标系下，即可实现扇形图像的绘制。

在绘制扇形图像时，有一些注意点如下：

我们总是从目标图像坐标反推原始图像坐标，而不是从原始图像坐标计算目标图像坐标。这样做可以避免映射过程中产生目标图像上的空隙问题（如小图像变换到大图像产生像素空隙），同时也方便我们在原始图像上进行非整数坐标上的插值运算。
我们需要将原点设置合适位置以绘制一个上下对称的扇形区域，如将 $(0,\frac{N}{2})$ 作为绘制原点。
目标图像尺寸由扫描扇形的角度决定，一般情况下我们保持目标图像的宽度为采样点 $N$ ，目标图像的高度由扫描角度决定，一般为 $2Nsin(\frac{\theta}{2})$ 。

从零构造超声图像变换

// scanangle: 超声图像扫描角度
// src.cols: 每条超声扫描线像素个数
// src.rows: 超声扫描线数量
void UltrasonicMap(cv::Mat& src, cv::Mat& dst, double scanangle)
{
	double half_scanangle = scanangle * .5;

	// 根据传入参数计算变换后图像尺寸
	// 基本规则：保持超声扫描线长度不变，即512分辨率扫描线被变换成扇形后长度仍然为512
	double proj_width = src.cols;
	double proj_height = src.cols * sin(half_scanangle) * 2.;

	// 计算变换后图像的整数尺寸（向上取整以容纳所有数据点）
	cv::Point dst_size(cvCeil(proj_width), cvCeil(proj_height));
	
	// 计算变换后极坐标中心点
	cv::Point2d circle_center(0, proj_height * .5);

	// 将第一条扫描线变换后目标图像上的向量（以circle_center为中心）
	double begin_angle = -scanangle * .5;
	cv::Point2d first_vec(1., tan(begin_angle));
	double vec_norm = sqrt(first_vec.x * first_vec.x +
		first_vec.y * first_vec.y);
	first_vec.x /= vec_norm;
	first_vec.y /= vec_norm;

	
	// 遍历变换后图像(目标图像)，计算目标图像上每一个像素对应的原始图像位置，并在原始图像上进行插值运算
	dst.create(dst_size, CV_32FC1);
	dst.setTo(cv::Scalar(0));
	for (int y = 0; y < dst_size.y; ++y)
	{
		float* datadst = dst.ptr<float>(y);
		for (int x = 0; x < dst_size.x; ++x)
		{
			// 确定目标图像点所在的扫描线向量
			cv::Point2d current_vec(x - circle_center.x, y - circle_center.y);
			double dr = sqrt(current_vec.x * current_vec.x + current_vec.y * current_vec.y);
			current_vec.x /= dr; current_vec.y /= dr;

			// 仅对符合条件的扫描线进行插值运算，条件包括：
			// 1 相对于第一条扫描线逆时针旋转的扫描线
			// 2 扫描线旋转角度范围小于scanangle
			// 3 扫描线像素个数小于src.cols - 1
			double rotate_direc = current_vec.x * first_vec.y - current_vec.y * first_vec.x;
			double rotate_radian = current_vec.x * first_vec.x + current_vec.y * first_vec.y;
			if (rotate_direc > 0 || rotate_radian < cos(scanangle) || dr >= src.cols - 1)
				continue;

			// 求目标图像点(x,y)在原图像上的坐标(x_src, y_src)，原图坐标为浮点数据，需要插值获取
			double theta = acos(rotate_radian);
			double y_src = theta / scanangle * (src.rows - 1);
			double x_src = dr;


			// 使用双线性插值：
			// 获取原图坐标(x_src, y_src)邻域上的4个整数坐标的值
			// val_00, val_01
			// val_10, val_11
			int x_floor = cvFloor(x_src);
			int y_floor = cvFloor(y_src);
			
			uchar* datasrc0 = src.ptr<uchar>(y_floor);
			uchar* datasrc1 = src.ptr<uchar>(y_floor + 1);

			uchar v00 = datasrc0[x_floor];
			uchar v01 = datasrc0[x_floor + 1];
			uchar v10 = datasrc1[x_floor];
			uchar v11 = datasrc1[x_floor + 1];

			// 获取水平方向(w0)与垂直方向(w1)线性插值权值
			double w0 = x_src - x_floor;  
			double w1 = y_src - y_floor;

			// 首先在水平方向上进行两次线性插值(也称作X方向上插值)
			double v0 = v00 * (1. - w0) + v01 * w0;
			double v1 = v10 * (1. - w0) + v11 * w0;

			// 使用水平方向上的插值v0,v1在垂直方向上(Y方向上)进行一次线性插值，
			// 获得(x_src, y_src)坐标上的插值，作为目标图像(x,y)坐标上的填充值
			double v = v0 * (1. - w1) + v1 * w1;

			// 填充目标图像
			datadst[x] = v;

		}
	}
}

int main() {

    // 读取超声扫描线图像
	cv::Mat img_rect = cv::imread("raw0.bmp", cv::IMREAD_GRAYSCALE);
	if (img_rect.empty())
		return 0;

	cv::Mat img_sector;
	double scanangle = CV_PI / 6. * 2;  // 扫描角度（单位为弧度）
	AntiAliasing(img_rect, scanangle);
	UltrasonicMap(img_rect, img_sector, scanangle);

	cv::imwrite("raw0_sector.bmp", img_sector);
	return 0;
}

以上代码实现了超声图像扫描线到扇形区域的变换，函数UltrasonicMap定义了几个关键的超声扫描参数：

scanangle: 超声图像扫描角度( $\theta$ )
src.cols: 每条超声扫描线像素个数( $M$ )
src.rows: 超声扫描线数量( $N$ )

变换后的扇形区域关于水平轴( $x$ 轴)对称，即在 $x$ 轴上下区域分别放置 $\frac{\theta}{2}$ 的扇形图像，其扇形中心点位于图像最左边且上下居中位置。

我们遍历目标图像上的所有点，分别计算出该点在原图像上的对应点，然后使用线性插值计算出对应点的灰度值，从而完成了超声图像的极坐标映射。

对于目标图像上的某些点，映射到原图后超出了有效范围，如大于最大扫描角度或者大于最大扫描像素，我们则不进行变换。

使用OpenCV进行超声图像变换

在上一章节中，我们从零手动构造了超声图像的变换，在变换过程中使用了双线性插值。OpenCV提供了cv::remap函数可以简化变换工作，同时该函数提供了多种插值方式以获得更佳的效果。

void cv::remap(
    InputArray  src,        // 输入图像
    OutputArray dst,        // 输出图像
    InputArray  map1,       // 映射表1
    InputArray  map2,       // 映射表2（可选）
    int         interpolation, // 插值方式
    int         borderMode = BORDER_CONSTANT, // 边界模式
    const Scalar& borderValue = Scalar()      // 边界填充值
);

src：输入图像，可以是灰度或彩色（任意通道）
dst：输出图像，大小由 map1 决定，类型与通道数与src保持一致
map1,map2：映射表，告诉 OpenCV输出图像中每个像素应该去源图像的哪个位置取值。
可使用单表模式，即map1类型为CV_32FC2，表示对应像素在原图像上的x,y坐标；
或双表模式，map1与map2的类型均为CV_32FC1，map1表示x坐标，map2表示y坐标。
interpolation：插值方式，具体可参考深入理解图像插值：从原理到应用
borderMode：当映射坐标超出 src 范围时，边界数据的处理方式

基于UltrasonicMap函数进行适当改写，可利用cv::remap实现同样的超声变换效果。我们首先构造map1，map2映射表，该表存储了目标图像对应位置在原图上的x,y坐标，然后调用remap实现重映射。

// scanangle: 超声图像扫描角度
// src.cols: 每条超声扫描线像素个数
// src.rows: 超声扫描线数量
void UltrasonicMapUsingOpenCV(cv::Mat& src, cv::Mat& dst, double scanangle)
{

	double half_scanangle = scanangle * .5;

	// 根据传入参数计算变换后图像尺寸
	// 基本规则：保持超声扫描线长度不变，即512分辨率扫描线被变换成扇形后长度仍然为512
	double proj_width = src.cols;
	double proj_height = src.cols * sin(half_scanangle) * 2.;

	// 计算变换后图像的整数尺寸（向上取整以容纳所有数据点）
	cv::Point dst_size(cvCeil(proj_width), cvCeil(proj_height));
	
	// 计算变换后极坐标中心点
	cv::Point2d circle_center(0, proj_height * .5);

	// 将第一条扫描线变换后目标图像上的向量（以circle_center为中心）
	double begin_angle = -scanangle * .5;
	cv::Point2d first_vec(1., tan(begin_angle));
	double vec_norm = sqrt(first_vec.x * first_vec.x +
		first_vec.y * first_vec.y);
	first_vec.x /= vec_norm;
	first_vec.y /= vec_norm;


	// 遍历变换后图像(目标图像)，计算目标图像上每一个像素对应的原始图像位置，
	// 填充到map_x,map_y矩阵中。如果没有对应坐标，则使用-1代替！
	cv::Mat map_x(dst_size.y, dst_size.x, CV_32FC1);
	cv::Mat map_y(dst_size.y, dst_size.x, CV_32FC1);
	map_x.setTo(cv::Scalar(-1.));
	map_y.setTo(cv::Scalar(-1.));
	for (int y = 0; y < dst_size.y; ++y)
	{
		float* datamapx = map_x.ptr<float>(y);
		float* datamapy = map_y.ptr<float>(y);
		for (int x = 0; x < dst_size.x; ++x)
		{
			// 确定目标图像点所在的扫描线向量
			cv::Point2d current_vec(x - circle_center.x, y - circle_center.y);
			double dr = sqrt(current_vec.x * current_vec.x + current_vec.y * current_vec.y);
			current_vec.x /= dr; current_vec.y /= dr;

			// 仅对符合条件的扫描线进行插值运算，条件包括：
			// 1 相对于第一条扫描线逆时针旋转的扫描线
			// 2 扫描线旋转角度范围小于scanangle
			// 3 扫描线像素个数小于src.cols - 1
			double rotate_direc = current_vec.x * first_vec.y - current_vec.y * first_vec.x;
			double rotate_radian = current_vec.x * first_vec.x + current_vec.y * first_vec.y;
			if (rotate_direc > 0 || rotate_radian < cos(scanangle) || dr >= src.cols - 1)
				continue;

			// 求目标图像点(x,y)在原图像上的坐标(x_src, y_src)
			double theta = acos(rotate_radian);
			double y_src = theta / scanangle * (src.rows - 1);
			double x_src = dr;

			// 填充映射表
			datamapx[x] = x_src;
			datamapy[x] = y_src;
		}
	}

	// 使用remap
	cv::remap(src, dst, map_x, map_y, cv::INTER_LINEAR);
}

UltrasonicMap与UltrasonicMapUsingOpenCV效果一致

进一步效果优化

观察超声图像在原点附近的图像质量，我们发现在原点附近的扇形超声图像存在混叠问题。混叠是采样率不足导致高频信号错误映射为低频信号的现象，在图像中表现为锯齿、摩尔纹等。这是什么原因引起的呢?

我们在深入理解图像插值：从原理到应用中曾说过当对大图像进行缩小时，常规插值处理可能产生混叠，从而引入了cv::INTER_AREA的插值方式，通过求一个区域的均值以避免缩小图像时的混叠问题。这里在原点附近，我们将512条扫描线上对应点压缩到几个或者几十个像素，使用cv::INTER_LINEAR插值显然会产生混叠问题。

然而，随着扫描深度增加，同样512条扫描线映射的长度会逐渐增加，甚至会大于512。这样，我们就没有办法选择一个统一的方式进行插值。甚至在同样缩小的区域，不同部位所使用的区域面积都存在差异。因此，我们需要进行自定义的反混叠处理。

在扫描图像的列方向上，对同一深度的扫描点应用统一的重采样处理。如使用cv::resize(roi, roi_resized, size, 0, 0, cv::INTER_AREA)函数实现重采样。
然后再恢复到原始尺寸。如使用cv::resize(roi_resized, roi, cv::Size(src_t.cols, 1), 0, 0, cv::INTER_LINEAR)恢复到原始尺寸。
不同的采样深度，使用不同的重采样尺寸，使得重采样后图像数量级与变换后图像数量级基本一致。
当采样深度大于某个阈值时，变换后图像尺寸大于变换前图像尺寸(或数量级基本一致)，不会产生混叠问题，则大于该阈值的采样深度点不再进行重采样处理。

void AntiAliasing(cv::Mat& src, double scanangle)
{
	// 计算扇形半径与最大弧长
	double scan_resolutions = src.cols;
	double scan_arc_max = scan_resolutions * scanangle;

	// 计算弧长等于扫描线一半时所对应的子扇形半径
	double half_scan_lines = src.rows * .5;
	double radius_threshold = half_scan_lines / scanangle;

	// 由于重采样基于roi处理(数据浅拷贝)，而我们期望roi数据是连续内存区域，
	// 通过转置实现后续的roi均在连续内存区域上进行，提升处理效率！
	cv::Mat src_t;
	cv::transpose(src, src_t);

	for (int row = 0; row < radius_threshold; ++row)
	{
		// 定义重采样图像roi
		cv::Rect roi_rect(0, row, src_t.cols, 1);
		cv::Mat roi = src_t(roi_rect);

		// 计算重采样后图像尺寸
		cv::Mat roi_resized;
		cv::Size size(row * scanangle + 1, 1);

		// 重采样
		cv::resize(roi, roi_resized, size, 0, 0, cv::INTER_AREA);

		// 恢复到原始尺寸
		cv::resize(roi_resized, roi, cv::Size(src_t.cols, 1), 0, 0, cv::INTER_LINEAR);
	}

	
	// 二次转置还原数据
	cv::transpose(src_t, src);

}

在进行超声图像变换前调用AntiAliasing可避免原点附近的混叠问题，使用如下：

int main() {

	cv::Mat img_rect = cv::imread("raw0.bmp", cv::IMREAD_GRAYSCALE);
	if (img_rect.empty())
		return 0;

	cv::Mat img_sector;
	double scanangle = CV_PI / 6. * 2;  // 扫描角度（单位为弧度）
	AntiAliasing(img_rect, scanangle);
	UltrasonicMap(img_rect, img_sector, scanangle);

	return 0;
}

总结

本文首先讲解了极坐标与笛卡尔坐标变换的基本关系，以及极坐标变换在超声图像扫描上的应用。然后提供了两种实现方案：手工实现的双线性插值变换和使用OpenCV的remap函数优化版本。针对变换过程中原点附近的混叠问题，构造了基于区域可变重采样的反混叠优化方法。通过坐标变换和图像处理技术，实现了超声扫描数据到扇形图像的精准可视化转换。

图像极坐标变换及超声图像上的应用

为什么需要极坐标变换？

极坐标与笛卡尔坐标的关系

超声图像扫描的基本原理

从零构造超声图像变换

使用OpenCV进行超声图像变换

进一步效果优化

总结

网站公告

今日签到

热门文章

最新发布