3D图像双线性插值

news/2024/5/22 8:22:28

文章目录

  • 前言
  • 结论
    • 说明:
    • 公式
  • 测试

前言

看了一下2d图像的双线性插值的理论,基本上都是在原图上找到对应的浮点坐标 p f p_f pf后,将以 p f p_f pf外围的4个点进行计算。计算的方法类似于二维直线方程的理论,但是写成了权重的方式。我看了一下,权重的方式还蛮好理解的。因为要用到,不自量力类推了一下3d图像插值,如果有错误,麻烦看到的大佬指出,感谢。

结论

在这里插入图片描述

说明:

v ( x , y , z ) v(x,y,z) v(x,y,z): ( x , y , z ) (x,y,z) (x,y,z)点处的value
v n e w v_{new} vnew: 新值
v o l d v_{old} vold: 旧值

公式

v n e w ( x , y , z ) =     v o l d ( x 1 , y 1 , z 2 ) ( x 2 − x ) ( y 2 − y ) ( z − z 1 ) + v o l d ( x 2 , y 2 , z 1 ) ( x − x 1 ) ( y − y 1 ) ( z 2 − z ) + v o l d ( x 2 , y 1 , z 2 ) ( x − x 1 ) ( y 2 − y ) ( z − z 1 ) + v o l d ( x 1 , y 2 , z 1 ) ( x 2 − x ) ( y − y 1 ) ( z 2 − z ) + v o l d ( x 2 , y 2 , z 2 ) ( x − x 1 ) ( y − y 1 ) ( z − z 1 ) + v o l d ( x 1 , y 1 , z 1 ) ( x 2 − x ) ( y 2 − y ) ( z 2 − z ) + v o l d ( x 1 , y 2 , z 2 ) ( x 2 − x ) ( y − y 1 ) ( z − z 1 ) + v o l d ( x 2 , y 1 , z 1 ) ( x − x 1 ) ( y 2 − y ) ( z 2 − z ) v_{new}(x,y,z)=\\ \ \ \ v_{old}(x1,y1,z2)(x2-x)(y2-y)(z-z1) \\+v_{old}(x2,y2,z1)(x-x1)(y-y1)(z2-z) \\+v_{old}(x2,y1,z2)(x-x1)(y2-y)(z-z1) \\+v_{old}(x1,y2,z1)(x2-x)(y-y1)(z2-z) \\+v_{old}(x2,y2,z2)(x-x1)(y-y1)(z-z1) \\+v_{old}(x1,y1,z1)(x2-x)(y2-y)(z2-z) \\+v_{old}(x1,y2,z2)(x2-x)(y-y1)(z-z1) \\+v_{old}(x2,y1,z1)(x-x1)(y2-y)(z2-z) vnew(x,y,z)=   vold(x1,y1,z2)(x2x)(y2y)(zz1)+vold(x2,y2,z1)(xx1)(yy1)(z2z)+vold(x2,y1,z2)(xx1)(y2y)(zz1)+vold(x1,y2,z1)(x2x)(yy1)(z2z)+vold(x2,y2,z2)(xx1)(yy1)(zz1)+vold(x1,y1,z1)(x2x)(y2y)(z2z)+vold(x1,y2,z2)(x2x)(yy1)(zz1)+vold(x2,y1,z1)(xx1)(y2y)(z2z)

测试

我在3d图像上测试了一下放大1.5倍,效果看起来还行
在这里插入图片描述

// 测试 3d 双线性插值

#include<iostream>
#include<itkImage.h>
#include<itkImageFileReader.h>
#include<itkImageFileWriter.h>
#include<itkNiftiImageIO.h>
#include<itkPNGImageIO.h>

using namespace std;

using PixelType = short;
const int Dimension = 3;
using ImageType = itk::Image<PixelType, Dimension>;
using ImagePointerType = ImageType::Pointer;


template<typename image_type, typename image_pointer>
void readData(const std::string& file_path, image_pointer& out_image) {
	using imageIOType = itk::NiftiImageIO;
	//using imageIOType = itk::PNGImageIO;
	using readerType = itk::ImageFileReader<image_type>;

	auto reader = readerType::New();
	auto imageIO = imageIOType::New();

	reader->SetImageIO(imageIO);
	reader->SetFileName(file_path);
	try {
		reader->Update();
	}
	catch (const itk::ExceptionObject& e) {
		cout << e.what() << endl;
	}

	out_image = reader->GetOutput();
}

template<typename image_type, typename image_pointer>
void writeData(const image_pointer& write_image, const std::string& write_path) {
	using writerType = itk::ImageFileWriter<image_type>;
	auto writer = writerType::New();

	using ImageIOType = itk::NiftiImageIO;
	// using ImageIOType = itk::PNGImageIO;
	auto ImageIO = ImageIOType::New();
	writer->SetImageIO(ImageIO);

	writer->SetInput(write_image);
	writer->SetFileName(write_path);
	writer->Update();
}

void createNewImage(
	ImagePointerType& new_image,
	const ImageType::SpacingType& spacing,
	const ImageType::DirectionType& direction,
	const ImageType::PointType& origin,
	const ImageType::SizeType& size)
{
	new_image = ImageType::New();
	ImageType::IndexType start;
	start[0] = 0; // x轴起始值
	start[1] = 0; // y轴起始值
	start[2] = 0; // z轴起始值
	ImageType::RegionType region;
	region.SetSize(size);
	region.SetIndex(start);
	new_image->SetRegions(region);
	new_image->Allocate();
	new_image->SetSpacing(spacing);
	new_image->SetDirection(direction);
	new_image->SetOrigin(origin);
	PixelType* regionBuffer = new_image->GetBufferPointer();
	for (int x = 0; x < size[0]; x++) {
		for (int y = 0; y < size[1]; y++) {
			for (int z = 0; z < size[2]; z++) {
				int index = size[0] * size[1] * z +
					size[0] * y + x;
				regionBuffer[index] = 0;
			}
		}
	}
}


void bilinearInterpolation(const ImagePointerType& input_image,
	ImagePointerType& out_image)
{
	// 取出各自的size
	ImageType::SizeType inputSize = input_image->GetBufferedRegion().GetSize();
	ImageType::SizeType outSize = out_image->GetBufferedRegion().GetSize();

	PixelType* inputBuffer = input_image->GetBufferPointer();
	PixelType* outBuffer = out_image->GetBufferPointer();

	int XInput = inputSize[0], YInput = inputSize[1], ZInput = inputSize[2];
	int XOut = outSize[0], YOut = outSize[1], ZOut = outSize[2];
	float deltax = XInput * 1.0 / XOut, deltay = YInput * 1.0 / YOut, deltaz = ZInput * 1.0 / ZOut;
	cout << deltax << ", " << deltay << ", " << deltaz << endl;
	// 目标: 求得outImage(xi,yi,zi)下的像素值
	for (int xi = 0; xi < XOut; ++xi) {
		for (int yi = 0; yi < YOut; ++yi) {
			for (int zi = 0; zi < ZOut; ++zi) {
				// 计算对应到inputImage的浮点坐标
				float x = xi * deltax, y = yi * deltay, z = zi * deltaz;
				// 计算(x1,y1,z1)
				int x1 = floor(x), y1 = floor(y), z1 = floor(z);
				// 计算(x2,y2,z2)
				int x2 = x1 + 1, y2 = y1 + 1, z2 = z1 + 1;

				int index = XOut * YOut * zi + XOut * yi + xi;

				if (x2 >= XInput || y2 >= YInput || z2 >= ZInput) {
					int index_x1y1z1 = XInput * YInput * z1 + XInput * y1 + x1;
					outBuffer[index] = inputBuffer[index_x1y1z1];
					continue;
				}
				// 执行插值计算
				// x1,y1,z2
				int index_x1y1z2 = XInput * YInput * z2 + XInput * y1 + x1;
				PixelType v_x1y1z2 = inputBuffer[index_x1y1z2];
				// cout << "here1" << endl;
				// x2,y2,z1
				int index_x2y2z1 = XInput * YInput * z1 + XInput * y2 + x2;
				PixelType v_x2y2z1 = inputBuffer[index_x2y2z1];
				// cout << "here2" << endl;
				// x2,y1,z2
				int index_x2y1z2 = XInput * YInput * z2 + XInput * y1 + x2;
				PixelType v_x2y1z2 = inputBuffer[index_x2y1z2];
				// cout << "here3" << endl;
				// x1, y2, z1
				int index_x1y2z1 = XInput * YInput * z1 + XInput * y2 + x1;
				PixelType v_x1y2z1 = inputBuffer[index_x1y2z1];
				// cout << "here4" << endl;
				// x2, y2, z2
				int index_x2y2z2 = XInput * YInput * z2 + XInput * y2 + x2;
				PixelType v_x2y2z2 = inputBuffer[index_x2y2z2];
				// cout << "here5" << endl;				
				// x1, y1, z1
				int index_x1y1z1 = XInput * YInput * z1 + XInput * y1 + x1;
				PixelType v_x1y1z1 = inputBuffer[index_x1y1z1];
				// cout << "here6" << endl;
				// x1, y2, z2
				int index_x1y2z2 = XInput * YInput * z2 + XInput * y2 + x1;
				PixelType v_x1y2z2 = inputBuffer[index_x1y2z2];
				// cout << "here7" << endl;
				// x2, y1, z1
				int index_x2y1z1 = XInput * YInput * z1 + XInput * y1 + x2;
				PixelType v_x2y1z1 = inputBuffer[index_x2y1z1];
				// cout << "here8" << endl;

				PixelType value = v_x1y1z2 * (x2 - x) * (y2 - y) * (z - z1)
					+ v_x2y2z1 * (x - x1) * (y - y1) * (z2 - z)
					+ v_x2y1z2 * (x - x1) * (y2 - y) * (z - z1)
					+ v_x1y2z1 * (x2 - x) * (y - y1) * (z2 - z)
					+ v_x2y2z2 * (x - x1) * (y - y1) * (z - z1)
					+ v_x1y1z1 * (x2 - x) * (y2 - y) * (z2 - z)
					+ v_x1y2z2 * (x2 - x) * (y - y1) * (z - z1)
					+ v_x2y1z1 * (x - x1) * (y2 - y) * (z2 - z);
				// 赋值							
				outBuffer[index] = value;			
			}
		}
	}
	
}


int main()
{
	const string inputImagePath = "G:/blood_vessel2023/images/train/62.nii.gz";
	ImagePointerType inputImage;
	readData<ImageType, ImagePointerType>(inputImagePath, inputImage);
	ImageType::SizeType inputSize = inputImage->GetBufferedRegion().GetSize();
	cout << inputSize << endl;	

	ImageType::SizeType outSize; // 1.5倍输入大小
	outSize[0] = inputSize[0] * 1.5; outSize[1] = inputSize[1] * 1.5; outSize[2] = inputSize[2] * 1.5;
	cout << outSize << endl;
	// char c = getchar();

	ImagePointerType outImage;
	createNewImage(outImage,
		inputImage->GetSpacing(), inputImage->GetDirection(), inputImage->GetOrigin(),
		outSize);

	// 插值
	bilinearInterpolation(inputImage, outImage);
	writeData<ImageType, ImagePointerType>(outImage, "bilinear.nii.gz");

	return 0;
}


http://lihuaxi.xjx100.cn/news/1175492.html

相关文章

数据结构 第四章:串

文章目录 一、串的定义和实现1.1串的定义和基本操作1.1.1串的定义1.1.2串的基本操作1.1.3小结 1.2串的存储结构1.2.1顺序存储1.2.2链式存储1.2.3基于顺序存储实现基本操作1.2.4小结 二、串的模式匹配2.1什么是字符串的模式匹配2.2朴素模式匹配算法2.3KMP算法2.4求next数组2.5KM…

Rinne Loves Graph

Rinne Loves Graph (nowcoder.com) 链接&#xff1a;登录—专业IT笔试面试备考平台_牛客网 来源&#xff1a;牛客网 题目描述 Island 发生了一场暴乱&#xff01;现在 Rinne 要和 Setsuna 立马到地上世界去。 众所周知&#xff1a;Island 是有一些奇怪的城镇和道路构成的…

如何用海外代理辅助对接 ChatGPT

许多朋友问我有没有好用的海外代理。说实话&#xff0c;真的好用的并不多。 最近我了解到了一家还不错的海外代理&#xff0c;叫做 IPIDEA&#xff0c;我已经使用了一段时间了&#xff0c;觉得质量挺不错。 你可能知道&#xff0c;我最近在进行一些 ChatGPT 相关的研究&#xf…

平衡二叉树的插入,删除以及平衡调整。

一&#xff0c;平衡二叉树插入失衡情况及解决方案 由于各种的插入导致的不平衡&#xff0c;每次调整都是最小不平衡子树。 LL&#xff1a;由于在结点A的 左孩子的左子树 插入结点导致失衡。 右单旋&#xff1a;①将A的 左孩子B 向右上旋转 代替A成为根节点       ②将A结…

springboot+vue广场舞团系统(java项目源码+文档)

风定落花生&#xff0c;歌声逐流水&#xff0c;大家好我是风歌&#xff0c;混迹在java圈的辛苦码农。今天要和大家聊的是一款基于springboot的广场团舞系统。项目源码以及部署相关请联系风歌&#xff0c;文末附上联系信息 。 &#x1f495;&#x1f495;作者&#xff1a;风歌&a…

带你深入了解RecyclerView的缓存机制

RecyclerView是Android中常用的列表显示控件&#xff0c;它具有强大的灵活性和性能优势。其中&#xff0c;RecyclerView的缓存机制是其实现高效滚动和快速显示大量数据的重要因素之一。在本文中&#xff0c;我们将详细解析RecyclerView的缓存机制&#xff0c;并附上相关示例代码…

springboot+vue大学生租房系统(java项目源码+文档)

风定落花生&#xff0c;歌声逐流水&#xff0c;大家好我是风歌&#xff0c;混迹在java圈的辛苦码农。今天要和大家聊的是一款基于springboot的大学生租房系统。项目源码以及部署相关请联系风歌&#xff0c;文末附上联系信息 。 &#x1f495;&#x1f495;作者&#xff1a;风歌…

挂耳式耳机哪个牌子好?这次推荐准没错!

在校园里&#xff0c;上网课、复习考证、刷视频……耳机总是我的第一选择&#xff0c;既不会打扰别人又能有效学习&#xff0c;普通的入耳式耳机戴久了总是会耳道痛&#xff0c;甚至出现耳道发炎问题&#xff0c;耳机也总是很难清洁。但是生活中我又离不开耳机&#xff0c;所以…