文章目录
- 前言
- 结论
- 说明：
- 公式
- 测试
前言
看了一下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)(x2−x)(y2−y)(z−z1)+vold(x2,y2,z1)(x−x1)(y−y1)(z2−z)+vold(x2,y1,z2)(x−x1)(y2−y)(z−z1)+vold(x1,y2,z1)(x2−x)(y−y1)(z2−z)+vold(x2,y2,z2)(x−x1)(y−y1)(z−z1)+vold(x1,y1,z1)(x2−x)(y2−y)(z2−z)+vold(x1,y2,z2)(x2−x)(y−y1)(z−z1)+vold(x2,y1,z1)(x−x1)(y2−y)(z2−z)
测试
我在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;
}
