一个简单的影像切片工具,生成xyz格式

本文介绍了一个使用C#编写的简单遥感影像切片工具,针对TIFF格式,提供基础界面和初步切片逻辑。通过MapTile方法处理切片,着重于代码优化和多线程应用。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

使用C#调用gdal写了个简单的遥感影像切xyz工具,只能切tiff影像。界面没设计,代码没优化,也没上多线程,有兴趣的自己拿去改。

 

 

全部代码如下,切片逻辑参见 MapTile 这个方法 

using System;
using System.Collections.Generic;
using System.Drawing;
using System.Windows.Forms;
using GDAL = OSGeo.GDAL;
using OGR = OSGeo.OGR;
using System.IO;

namespace MapTileTool
{
    public partial class MainForm : Form
    {
        public MainForm()
        {
            InitializeComponent();
        }

        private void MainForm_Load(object sender, EventArgs e)
        {
           
        }

        private void HandleGeoTiffSelect(object sender, EventArgs e)
        {
            string title = "选择遥感影像";
            string filter = "tif文件(*.tif,*.tiff)|*.tif;*.tiff";
            string fileName = FileChooserSelect(title, filter);
            if (null != fileName) {
                this.tifTextBox.Text = fileName;
            }
        }
        private void HandleGeoTiffRest(object sender, EventArgs e)
        {
            this.tifTextBox.Text = ""; ;
        }

        private void HandleOutputSelect(object sender, EventArgs e)
        {
            string outputPath = FoloderChooserSelect("选择输出文件夹");
            if (null != outputPath) {
                this.outputTextBox.Text = outputPath;
            }
        }

        private void HandleOutputRest(object sender, EventArgs e)
        {
            this.outputTextBox.Text = "";
        }

        private void ExecuteMapTile(object sender, EventArgs e)
        {
            string geoTiffPath = tifTextBox.Text;
            string outputPath = outputTextBox.Text;
            string zMin = zMinTextBox.Text;
            string zMax = zMaxTextBox.Text;

            if (string.IsNullOrEmpty(geoTiffPath)
                || string.IsNullOrEmpty(outputPath)
                || string.IsNullOrEmpty(zMin)
                || string.IsNullOrEmpty(zMax)) {
                MessageBox.Show("信息填写有误");
                return;
            }

            try
            {
                int zMinLevel = Convert.ToInt32(zMin);
                int zMaxLevel = Convert.ToInt32(zMax);
                executeBtn.Enabled = false;
                consoleLabel.Text = "正在切片,请稍候。。。";
                MapTile(geoTiffPath, outputPath, zMinLevel, zMaxLevel);
                consoleLabel.Text = "恭喜,影像切片已完成";
                executeBtn.Enabled = true;
            }
            catch (Exception ex) {
                Console.WriteLine(ex.Message);
                consoleLabel.Text = "切图异常:" + ex.Message;
                executeBtn.Enabled = true;
            }

        }

        /**
         * 切片处理逻辑
         */
        private void MapTile(string geoTiffPath, string outputPath, int zMin, int zMax) {
            try
            {
                OGR.Ogr.RegisterAll();
                GDAL.Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
            }
            catch (Exception ex)
            {
                Console.WriteLine("Gdal驱动注册失败:" + ex.Message);
            }

            OSGeo.GDAL.Dataset dataset = GDAL.Gdal.Open(geoTiffPath, OSGeo.GDAL.Access.GA_ReadOnly);
            int width = dataset.RasterXSize;
            int height = dataset.RasterYSize;
            int bandCount = dataset.RasterCount;
            double[] transform = new double[6];
            dataset.GetGeoTransform(transform);
            double lonMin = transform[0];
            double lonMax = transform[0] + (width * transform[1]);
            double latMin = transform[3] + (height * transform[5]);
            double latMax = transform[3];

            if (zMin < 6)
            {
                zMin = 6;
            }

            if (zMax > 18)
            {
                zMax = 18;
            }

            for (int z = zMin; z <= zMax; z++)
            {
                int tileRowMin = GetXyTileByZ(lonMin, latMax, z)[0];
                int tileRowMax = GetXyTileByZ(lonMax, latMin, z)[0];
                int tileColMin = GetXyTileByZ(lonMin, latMax, z)[1];
                int tileColMax = GetXyTileByZ(lonMax, latMin, z)[1];
                double tempLonMin = GetLngLatByXyz(tileRowMin, tileColMin, z)[0];
                double tempLonMax = GetLngLatByXyz(tileRowMin + 1, tileColMin, z)[0];
                double tempLatMin = GetLngLatByXyz(tileRowMin, tileColMin + 1, z)[1];
                double tempLatMax = GetLngLatByXyz(tileRowMin, tileColMin, z)[1];

                // 获取X轴方向分辨率
                double xResolution = (tempLonMax - tempLonMin) / 256;
                // 获取Y轴方向分辨率
                double yResolution = (tempLatMax - tempLatMin) / 256;

                for (int x = tileRowMin; x <= tileRowMax; x++)
                {
                    for (int y = tileColMin; y <= tileColMax; y++)
                    {
                        double tileLonMin = tempLonMin + (x - tileRowMin) * xResolution * 256;
                        double tileLatMax = tempLatMax - (y - tileColMin) * yResolution * 256;

                        try
                        {
                            double xTileMax = tileLonMin + 256 * xResolution;
                            double yTileMin = tileLatMax - 256 * yResolution;
                            List<int> positionMin = GetPosition(tileLonMin, tileLatMax, transform);
                            List<int> positionMax = GetPosition(xTileMax, yTileMin, transform);
                            Boolean contain = positionMin[0] < width && positionMin[0] >= 0
                                    && positionMin[1] < height && positionMin[1] >= 0
                                    && positionMax[0] < width && positionMax[0] >= 0
                                    && positionMax[1] < height && positionMax[1] >= 0;

                            Bitmap bitmap = new Bitmap(256, 256);
                            if (contain)
                            {
                                int xDistant = positionMax[0] - positionMin[0] + 1;
                                int yDistant = positionMax[1] - positionMin[1] + 1;
                                bitmap = new Bitmap(xDistant, yDistant);
                                if (bandCount == 1)
                                {
                                    float[] gs = new float[xDistant * yDistant];
                                    dataset.GetRasterBand(1).ReadRaster(positionMin[0], positionMin[1], xDistant, yDistant, gs, xDistant, yDistant, 0, 0);
                                    for (int i = 0; i < xDistant; i++)
                                    {
                                        for (int j = 0; j < yDistant; j++)
                                        {
                                            float g = gs[xDistant * j + i];
                                            if (g < 0)
                                            {
                                                continue;
                                            }
                                            int gray = (int)(255 * g);
                                            Color color = Color.FromArgb(gray, gray, gray);
                                            bitmap.SetPixel(i, j, color);
                                        }
                                    }

                                }
                                else
                                {
                                    int[] rs = new int[xDistant * yDistant];
                                    int[] gs = new int[xDistant * yDistant];
                                    int[] bs = new int[xDistant * xDistant];
                                    dataset.GetRasterBand(1).ReadRaster(positionMin[0], positionMin[1], xDistant, yDistant, rs, xDistant, yDistant, 0, 0);
                                    dataset.GetRasterBand(2).ReadRaster(positionMin[0], positionMin[1], xDistant, yDistant, gs, xDistant, yDistant, 0, 0);
                                    dataset.GetRasterBand(3).ReadRaster(positionMin[0], positionMin[1], xDistant, yDistant, bs, xDistant, yDistant, 0, 0);

                                    for (int i = 0; i < xDistant; i++)
                                    {
                                        for (int j = 0; j < yDistant; j++)
                                        {
                                            int r = rs[xDistant * j + i];
                                            int g = gs[xDistant * j + i];
                                            int b = bs[xDistant * j + i];
                                            if (r < 0 || r > 255 || g < 0 || g > 255 || b < 0 || b > 255)
                                            {
                                                continue;
                                            }

                                            Color color = Color.FromArgb(r, g, b);
                                            bitmap.SetPixel(i, j, color);
                                        }
                                    }
                                }
                            }
                            else
                            {
                                for (int i = 0; i < 256; i++)
                                {
                                    for (int j = 0; j < 256; j++)
                                    {
                                        double tileLon = tileLonMin + i * xResolution;
                                        double tileLat = tileLatMax - j * yResolution;
                                        List<int> position = GetPosition(tileLon, tileLat, transform);
                                        int xPosition = position[0];
                                        int yPosition = position[1];

                                        if (xPosition < 0 || yPosition < 0 || xPosition >= width || yPosition >= height)
                                        {
                                            continue;
                                        }
                                        Color color;
                                        if (bandCount == 1)
                                        {
                                            float[] values = new float[1];
                                            dataset.GetRasterBand(1).ReadRaster(xPosition, yPosition, 1, 1, values, 1, 1, 0, 0);
                                            if (values[0] < 0)
                                            {
                                                continue;
                                            }
                                            int gray = (int)(255 * values[0]);
                                            color = Color.FromArgb(gray, gray, gray);
                                        }
                                        else
                                        {
                                            int[] rs = new int[1];
                                            int[] gs = new int[1];
                                            int[] bs = new int[1];
                                            dataset.GetRasterBand(1).ReadRaster(xPosition, yPosition, 1, 1, rs, 1, 1, 0, 0 );
                                            dataset.GetRasterBand(2).ReadRaster(xPosition, yPosition, 1, 1, gs, 1, 1, 0, 0);
                                            dataset.GetRasterBand(3).ReadRaster(xPosition, yPosition, 1, 1, bs, 1, 1, 0, 0);
                                            int r = rs[0];
                                            int g = gs[0];
                                            int b = bs[0];
                                            if (r < 0 || r > 255 || g < 0 || g > 255 || b < 0 || b > 255)
                                            {
                                                continue;
                                            }

                                            color = Color.FromArgb(r, g, b);
                                        }
                                        bitmap.SetPixel(i, j, color);
                                    }
                                }
                            }

                            try
                            {
                                String layerPath = outputPath + "/" + z + "/" + x;
                                if (!System.IO.Directory.Exists(layerPath))
                                {
                                    System.IO.Directory.CreateDirectory(layerPath);
                                }

                                String imagePath = layerPath + "/" + y + ".png";
                                bitmap.Save(imagePath);
                            }
                            catch (Exception e)
                            {
                                Console.WriteLine("切图异常:" + e.Message);
                                consoleLabel.Text = "切图异常:" + e.Message;
                            }
                        }
                        catch (Exception e)
                        {
                            Console.WriteLine("切图异常:" + e.Message);
                            consoleLabel.Text = "切图异常:" + e.Message;
                        }
                    }
                }
            }

        }


        /**
         * 获取行列号
         *
         * @param lng
         * @param lat
         * @param z
         */
        public static List<int> GetXyTileByZ(double lng, double lat, int z)
        {
            List<int> xy = new List<int>(2);
            int x = (int)Math.Floor((lng + 180) / 360 * Math.Pow(2, z));
            int y = (int)Math.Floor((1 - Math.Log(Math.Tan(Math.PI*lat/180) + 1 / Math.Cos( Math.PI*lat/180)) / Math.PI) / 2 * (Math.Pow(2, z)));
            xy.Add(x);
            xy.Add(y);
            return xy;
        }

        /**
         * 根据经纬度获取影像像素位置 
         */
        public static List<int> GetPosition(double lng, double lat, double[] transform)
        {
            double dTemp = transform[1] * transform[5] - transform[2] * transform[4];
            int xPix = (int)((transform[5] * (lng - transform[0]) - transform[2] * (lat - transform[3])) / dTemp);
            int yPix = (int)((transform[1] * (lat - transform[3]) - transform[4] * (lng - transform[0])) / dTemp);
            List<int> xy = new List<int>(2);
            xy.Add(xPix);
            xy.Add(yPix);
            return xy;
        }

        /**
         * 根据切片序列号获取经纬度坐标
         * @return
         */
        public static List<Double> GetLngLatByXyz(int x, int y, int z)
        {
            List<Double> xy = new List<Double>(2);
            double n = Math.Pow(2, z);
            double lng = x / n * 360.0 - 180.0;
            double lat = Math.Atan(Math.Sinh(Math.PI * (1 - 2 * y / n)));
            lat = lat * 180.0 / Math.PI;
            xy.Add(lng);
            xy.Add(lat);
            return xy;
        }


        /**
        * 选择文件控件
        */
        private static string FileChooserSelect(string title, string filter)
        {
            string fileName = "";
            OpenFileDialog fileDialog = new OpenFileDialog();
            fileDialog.Multiselect = true;
            fileDialog.Title = title;
            fileDialog.Filter = filter;
            if (fileDialog.ShowDialog() == DialogResult.OK)
            {
                fileName = fileDialog.FileName;
            }

            return fileName;
        }

        /**
         * 选择文件夹控件
         */
        private static string FoloderChooserSelect(string title)
        {
            string folderPath = "";
            FolderBrowserDialog folderDialog = new FolderBrowserDialog();
            folderDialog.Description = title;
            if (folderDialog.ShowDialog() == DialogResult.OK)
            {
                folderPath = folderDialog.SelectedPath;
            }

            return folderPath;
        }
    }

}
 

<think>嗯,用户想了解如何从DICOM医学影像文件中提取三维数组并生成点云。首先,我需要回忆一下DICOM文件的结构和处理方法。DICOM是医学影像的标准格式,通常存储的是CT、MRI等扫描数据,每个文件可能代表一个切片,多个切片组合成三维数据。 首先,提取三维数组的话,可能需要用到专门的库,比如Python中的pydicom。每个DICOM文件读取后,像素数据存储在PixelData元素中,转换成numpy数组。然后,将多个切片堆叠起来形成三维数组。记得调整方向,因为不同设备的扫描方向可能不同,可能需要调整坐标轴。 接下来是生成点云。三维数组中的每个体素对应一个三维坐标,但直接转换会导致数据量太大,可能需要降采样或者阈值处理。比如,使用Marching Cubes算法提取等值面,或者根据阈值选择特定灰度值的体素,转换为点坐标。这需要用到像VTK或者SimpleITK这样的工具,或者Python中的scikit-image。 然后,用户提到参考中的OpenGL顶点、索引、法线等信息。生成点云后,可能需要进一步处理顶点和法线,但点云通常只需要顶点位置。如果要渲染,可能需要IBO和VBO,但用户的问题可能只需要点云数据,所以可能不需要索引和法线。 工具方面,Python的pydicom处理读取,SimpleITK或VTK处理三维数据,PyVista或Open3D可视化和导出点云。需要注意内存管理,因为医学影像数据量大,容易导致内存不足,可能需要分块处理或者优化。 用户可能还关心具体的代码示例,比如如何读取DICOM系列,转换为numpy数组,应用阈值生成点坐标,保存为PLY或XYZ格式。此外,可能需要考虑坐标系的转换,比如患者坐标系到世界坐标系的转换,这在DICOM中有相关的元数据,比如Image Position Patient和Image Orientation Patient,需要应用到每个切片的位置上,以确保点云的空间正确性。 还要注意,直接阈值法可能不够精确,可能需要更高级的分割算法,比如区域生长、水平集等,但这可能超出用户当前的需求。用户可能只需要基础的方法,所以保持回答简洁,提供主要步骤和工具即可。 总结步骤:读取DICOM文件系列,转换为三维数组,处理坐标转换,应用阈值或算法生成点云坐标,导出为点云格式。需要使用的库和工具,以及可能的优化方法,比如降采样。</think>### 从DICOM文件提取三维数组生成点云的方法 #### 1. **DICOM数据读取与三维数组构建** - **工具选择**:使用Python的`pydicom`库或`SimpleITK`读取DICOM文件。医学影像通常由多层二维切片构成,需将多个切片按顺序堆叠为三维数组。 - **代码示例**: ```python import pydicom import numpy as np import os # 读取DICOM序列(假设文件按顺序存放在文件夹中) dicom_files = sorted([os.path.join(dicom_dir, f) for f in os.listdir(dicom_dir)]) slices = [pydicom.dcmread(f) for f in dicom_files] volume = np.stack([s.pixel_array for s in slices], axis=2) # 三维数组(z,y,x) ``` - **注意事项**: - 需处理DICOM元数据中的`Image Position (Patient)`和`Image Orientation (Patient)`字段,确保三维空间坐标系正确[^2]。 #### 2. **三维数组到点云转换** - **阈值提取**:根据灰度值筛选感兴趣区域(如骨骼、器官),例如: ```python threshold = 200 # 根据实际影像调整阈值 coordinates = np.argwhere(volume > threshold) # 获取体素坐标(z,y,x) ``` - **坐标转换**:将体素坐标转换为真实物理坐标(毫米单位),需结合DICOM的`Pixel Spacing`和`Slice Thickness`参数: ```python pixel_spacing = slices[0].PixelSpacing # (y方向间距, x方向间距) slice_thickness = slices[0].SliceThickness physical_coords = coordinates * np.array([slice_thickness, pixel_spacing[0], pixel_spacing[1]]) ``` - **优化方法**:若数据量过大,可通过降采样减少点云密度,例如每间隔2个体素采样一次。 #### 3. **点云导出与可视化** - **格式导出**:将坐标保存为PLY或XYZ格式: ```python with open("output.ply", "w") as f: f.write("ply\nformat ascii 1.0\n") f.write(f"element vertex {len(physical_coords)}\n") f.write("property float x\nproperty float y\nproperty float z\nend_header\n") for pt in physical_coords: f.write(f"{pt[2]} {pt[1]} {pt[0]}\n") # 注意坐标轴方向转换 ``` - **可视化工具**:使用`CloudCompare`、`ParaView`或Python库`pyntcloud`进行点云查看和分析。 #### 4. **高级方法(等值面提取)** - **Marching Cubes算法**:通过提取等值面生成表面网格(更适合三维重建): ```python from skimage.measure import marching_cubes verts, faces, normals, _ = marching_cubes(volume, level=threshold) # verts即为顶点坐标(需转换为物理坐标) ``` #### 工具与库总结 | 步骤 | 工具/库 | 用途 | |--------------|-------------------|-----------------------------| | DICOM解析 | pydicom, SimpleITK| 读取元数据和像素阵列 | | 三维处理 | NumPy, VTK | 数组操作、坐标转换 | | 点云生成 | Open3D, pyntcloud| 点云构建、降采样、导出 | | 可视化 | CloudCompare | 交互式查看点云 |
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值