Rust教程 – 学习天文图像的多尺度处理

最近，人们投入了大量精力开发新颖的图像处理技术。其中许多技术都源自于傅里叶和小波变换等数字信号处理方法。

这些技术不仅使得各种图像处理技术如降噪、锐化和动态范围扩展成为可能，而且还使得计算机视觉中使用的许多技术如边缘检测、目标检测等成为可能。

多尺度分析是相对较新的技术之一，已经在广泛的应用中得到采用，特别是在天文图像和数据处理应用中。这种基于小波变换的技术允许我们将数据分解为多个信号，所有这些信号加在一起形成最终信号。

然后，我们可以对这些单独的子信号进行处理或分析工作，从而能够进行有针对性的操作，而不影响其他子信号。

在本文中，我们首先将通过一种特定的算法来探讨这种技术的本质。然后，我们将进一步看看如何在Rust编程语言中实现我们在第一部分讨论的内容，并重新创建文章前半部分中所看到的示例。

（本文视频讲解：java567.com）

阅读之前：

第一部分的先决条件：

所描述的技术源自于“小波变换”的概念。你不需要对它了解得很深，但基本的了解会帮助你更好地理解材料。

由于本文集中于图像处理和分析，了解数字格式中像素是如何工作的基本知识是有帮助的，但不是必需的。

第二部分的先决条件：

在这里，我们专注于使用Rust编程语言实现算法，而不深入探讨语言本身的细节。因此，需要熟悉编写Rust程序并且能够阅读crate文档。

如果这不适用于你，你仍然可以阅读第一部分并学习这种技术，然后也许你会想尝试在你选择的语言中使用它。如果你不熟悉Rust，我强烈建议你学习基础知识。这里有一个可以让你入门的交互式Rust课程。

目录

1. 第一部分：理解多尺度处理技术和算法
   1. 什么是多尺度图像处理
   2. À Trous小波变换
   3. 缩放函数
   4. 每个尺度的卷积像素
   5. 处理边界条件
   6. 计算任意给定图像的最大可能尺度
   7. 总结说明
2. 第二部分：如何在 Rust 中实现 À Trous 变换
   1. 先决条件
   2. À Trous 变换
   3. 迭代器和 À Trous 变换
   4. 卷积
   5. 实现迭代器
   6. 重组
   7. 使用 À Trous 变换
   8. 进一步阅读
3. 总结

第一部分：理解多尺度处理技术和算法

当我们谈论对某些数据进行多尺度处理或分析时，我们指的是什么？嗯，我们通常是指将输入数据分解为多个信号，每个信号代表特定尺度的信息。

在图像分析中，尺度简单地指的是我们在任何给定时间内所观察到的结构的大小。它忽略了当前尺度之外的所有东西，无论是更小还是更大的。

什么是多尺度图像处理？

对于图像，“尺度”通常指图像中各种结构或细节的像素大小。通过以下示例，你将能够直观地理解：

梅西耶33号，又称三角座星系

假设我们的朴素理解是正确的，我们可以至少得到以下3个尺度的图像：

• 非常小的结构，通常是单个像素的大小。这一层，当与图像的其余部分分离时，主要只包含噪声和一些尖锐的星星。
• 较小的结构，通常几个像素大小。这一层，当分离时，将包含所有的星星和星系臂中非常细微的细节。
• 大型和非常大型的结构，通常是100个像素以上。这一层，当分离时，将包含星系中心的一般大小和形状。

现在问题变成了，为什么我们需要一开始就做所有这些呢？

答案很简单：它使我们能够对图像进行有针对性的增强和改变。

例如，对整体图像进行降噪通常会导致星系的清晰度降低。但由于我们已经将图像分解为多个尺度，我们可以轻松地将降噪应用于前几层，因为大多数易于去除的随机噪声仅存在于较低尺度的层中。

然后，我们将降噪后的低尺度层与未修改的大尺度层重新组合，我们就得到了一个输出，它在降噪的同时保持了质量。

噪声的另一个奇特之处在于，它几乎总是存在于这些层中的一个，使得降噪过程既简单又非破坏性。

如果你更喜欢以视觉方式学习，让我们通过使用上面使用的图像来实践。我们将使用以下这张图像的灰度版本，我也添加了随机高斯噪声：

梅西耶33号 AKA 三角座星系，转换为灰度并添加高斯噪声

对该图像执行基于尺度的层分离，我们得到以下结果。请注意，为了展示目的，结果已经重新调整到一个范围，以便可以将其视为图像。实际的变换生成的像素值在单独查看时没有意义，但是本教程中描述的所有技术和计算仍然可以安全地应用而不需要重新调整。重新组合过程会自动给我们返回正确的范围：

9级 À Trous 分解。从左上到右下，我们有以下像素尺度的图像：1、2、4、8、16、32、64、128、256（2的幂）

1. 第一层和第二层包含噪声和星星。在这个特定的示例中，噪声与星星混在一起。但是使用第一层和第二层，我们可以轻松地针对不在第二层中的区域进行定位，因为我们可以确信这些区域是第一层中存在噪声的地方。
2. 对于第三层，我们仍然可以看到星星残余的亮度。但是如果你仔细观察，你会发现星系臂开始非常微弱地出现。
3. 从第四层开始，我们看到星系在不同尺度和细节级别上的图像，完全不包含星星。我们从更细微的细节（相对较小的尺度细节）开始，并逐渐转移到越来越大的尺度样本。最后，我们只能看到星系曾经所在的一个模糊的形状。

从这里开始，我们可以选择性地将降噪应用于前两层。然后，我们重新组合所有层以创建以下图像，该图像几乎没有噪声，同时保留了星星和星系臂中相同数量的细节：

梅西耶33号 AKA 三角座星系，重新组合所有层但对像素尺度1和2的层应用了降噪处理

在其最基本的形式中，多尺度分析涉及将您的源图像（通常称为“信号”）分解为多个“信号” – 每个信号包含源信号中特定尺度的数据。

在这里谈论图像信号时，尺度指的是我们从源图像创建层时采取的相邻像素之间的距离。

在实践中，这种技术被用作各种天文数据分析和图像处理的第一步。

例如，你可以使用该技术更容易地检测星星的位置，同时忽略比以其他方式更容易的大型结构。

À Trous 小波变换

我之前向你展示的所有内容，以及你将在本教程中看到的所有内容，都是使用小波分解和重组，使用 À Trous 算法进行的离散小波变换。

这个算法多年来一直被用于各种应用。但是它最近在天文图像处理应用中变得尤为重要，其中图像中的不同对象和信号可以完全基于结构尺度进行分离。

以下是该算法的工作原理：

1. 我们从源图像输入和分解为n级别的开始。
2. 对于每个级别n：
   • 我们使用我们的缩放函数对图像进行卷积（稍后我们会看到这是什么），其中相邻像素被认为相距 2n 个单位，给出结果 resultn。这就是“À Trous”名称的由来，字面上翻译为“带有孔洞”。
   • 然后，使用 input - resultn 计算出层输出 outputn。
   • 然后更新 input 等于 resultn。这也称为残差数据，它作为下一层的源数据。
3. 重复上述步骤，直到所有级别完成。
4. 最后，我们有9个小波层和1个残差层。所有10个层都需要用于重组。

对于更多理解这个算法的数学方法，我鼓励你阅读关于 à trous 算法的这里 。

重组过程非常简单：我们只需要将所有10层相加即可。我们可以选择对任何层应用正或负的 偏差，这是一个在重组期间乘以该层像素值的因子。您可以使用它来增强或减弱特定层的特征。

缩放函数

缩放函数是特定的卷积核，可以帮助我们更好地表示特定尺度的数据，根据我们的用例。以下是3种最常用的缩放函数：

梅西耶33号 AKA 三角座星系的图像的第3级分解可视化的 À Trous 算法中最常用的3种缩放函数。

• B3样条是一个非常平滑的核。它主要用于单独处理大型结构。如果我们想要锐化我们的星系，我们将使用这个核心。
• 低尺度是一个非常尖锐的核心，最适合处理小尺度结构。
• 线性插值核心给了我们两全其美，因此在需要处理小尺度和大尺度结构时使用。这是我们之前所有示例中使用的方法。

在每个尺度上的卷积像素

在算法中我提到，在每个尺度上，图像中的像素被认为相距2n个单位。让我们试着通过以下可视化更好地理解这一点：

考虑以下8像素乘以8像素的图像。每个像素都标有从1到64的索引，这是它们的索引。

一个8x8像素图像的代表性像素网格

在这个例子中，我们只关注一个中心像素的卷积操作，假设是像素编号28。

尺度 0： 在尺度0上，2n的值变为1。这意味着在卷积中，我们将考虑距离我们目标中心像素1个单位的像素。下面突出显示了这些像素：

在尺度0上，突出显示了像素编号28的卷积像素

尺度 1： 这是有趣的地方。在尺度1上，2n的值变为2。这意味着在卷积中，我们将直接跳到距离目标像素2个位置的像素：

在尺度1上，突出显示了像素编号28的卷积像素

正如你所看到的，我们通过跳过2n - 1个相邻像素并选择第2n个像素，在计算目标像素的值时创建了“孔洞”。这就是算法的基础。

这个过程对图像中的每个像素都会重复进行，就像常规的卷积过程一样。每次，我们都会考虑像素之间的距离随着尺度增加而增加。

让我们再看一个尺度。

尺度 2： 这里更加有趣。在尺度2上，2n的值变为4。这意味着在卷积中，我们将直接跳到距离目标像素4个位置的像素：

在尺度2上，突出显示了像素编号28的卷积像素

等等？为什么我们选择像素1、4、8、25和57？1和4之间只相距3个位置，25只相距2个位置，8和57甚至与目标像素不对角对齐。发生了什么？

处理边界条件

正如我们提到的，这个过程适用于图像中的所有像素，我们还需要考虑卷积像素的位置是否位于图像之外。

这不是这个算法独有的概念。在卷积过程中，这被称为边界条件或处理边界像素。有各种各样的技术来处理这个问题，它们都涉及到虚拟地扩展图像，以使得我们似乎根本没有遇到边界。

一些技术包括：

• 根据需要扩展，通过复制最后一行/列的值
• 在所有边缘和角落镜像图像
• 在边缘周围包裹图像。

在我们的例子中，我们采用了“镜像”技术。当实现这样的算法时，我们不需要实际创建扩展的图像。任何边界处理都可以使用基本的数学公式来实现。

我们的扩展图像，在尺度2中选择了正确的像素，如下所示：

使用镜像技术在所有边缘和角落扩展的源图像。所有淡化的区域代表了扩展区域。

再次强调，扩展仅仅是逻辑上的，完全使用公式计算，而不是实际扩展源图像然后检查。我们可以很容易地看到，在放置了镜像图像后，我们仍然遵循选择相距2n个位置的像素的基本规则。

计算任何给定图像的最大可能尺度

仔细思考一下，你会发现图像可以分解的最大层数可以通过计算图像宽度或高度（较低者）的log2，并且丢弃小数部分来计算。

对于我们的5x5图像，log2(5) ~= 2.32。如果我们丢弃小数部分，那么剩下的就是2层。同样，对于一个1000x1000像素的图像，log21000 ~= 9.96，这意味着我们可以将一个1000x1000像素的图像分解成最多9层。这简单地意味着我们的“孔洞”不能大于图像的宽度或高度。

即使使用了上面所述的镜像扩展，如果孔洞大于图像的宽度，它们仍将位于扩展区域之外，特别是对于角落或边界像素，这将使得在该尺度上执行卷积变得不可能。

总结说明

再仔细考虑一下这些示例和可视化，你可以清楚地看到为什么这个算法起作用，以及它是如何根据结构的尺寸将图像中的结构分离出来的。逐渐增加的孔洞尺寸使得在任何给定层中只保留比孔洞本身大的结构。

使用这个算法的一个重要优点是计算成本。由于这不涉及傅立叶或小波变换，

相对而言计算成本很低。然而，内存成本确实较高。但往往这是一个很好的权衡。

与其他离散小波变换算法相比，当源图像的大小在整个过程中保持不变时，这个算法的另一个优点是。这里没有发生降采样或上采样，使得这个算法成为理解和实现最容易的算法之一。

这个算法在几乎所有的天文图像处理软件中都得到了使用，比如PixInsight等等。

这个算法也以其他名称而闻名，如Stationary Wavelet Transform和Starlet Transform。

第二部分：如何在 Rust 中实现 À Trous 变换

现在我将向你展示如何在 Rust 中实现这个算法。

在本教程中，我将假设你对 Rust 及其基本概念（如数据类型、迭代器和特性）相当熟悉，并且能够轻松编写使用这些概念的程序。

我还将假设你了解在这个上下文中卷积和卷积核的含义。

先决条件

我们需要一些依赖项。在我们开始之前，让我们快速创建一个新项目：

cargo new --lib atrous-rs
cd atrous-rs

现在让我们添加我们需要的所有依赖项。我们实际上只需要 2 个：

cargo add image ndarray

image 是一个 Rust 库，我们将使用它来处理所有标准格式和编码的图像。它还可以帮助我们在各种格式之间进行转换，并提供对像素数据的轻松访问。

ndarray 是一个 Rust 库，它可以帮助你创建、操作和处理二维、三维或多维数组。我们可以使用嵌套的向量，但在这种情况下使用像 ndarray 这样的项目更好，因为我们需要对单个值以及它们的相邻值执行许多操作。不仅使用 ndarray 更容易，而且对许多操作和 CPU 类型也进行了性能优化。

虽然我将涵盖我们从这些 crate 中使用的基本函数/特性/方法/数据类型，但我不打算为它们详细说明。我建议你阅读文档。

我们实际上将直接跳到算法实现，并稍后回来看如何使用它。

À Trous 变换

创建一个新文件来保存我们的实现。让我们把它命名为 transform.rs。

首先添加以下结构体，它将保存我们执行变换所需的信息：

// transform.rs

use ndarray::Array2;

pub struct ATrousTransform {
	input: Array2<f32>, // `Array2<f32>` 是一个二维数组，每个值都是 `f32` 类型。这将保存我们输入图像的像素数据。
    levels: usize, // 将图像分解为的级别或尺度的数量
    current_level: usize, // 我们需要生成的当前级别。这保存了我们迭代器的状态。
    width: usize, // 输入图像的宽度
    height: usize, // 输入图像的高度
}

我们还需要一种轻松创建此结构的方法。在我们的情况下，我们希望能够直接从输入图像创建它。此外，输入图像可以是任何受支持的格式和编码，但我们希望将图像转换为我们期望的一致的颜色类型以执行计算，因此我们还需要将图像转换为我们期望的格式。

使用 Rust 中的“构造函数”模式将所有这些逻辑提取到一起是有帮助的。让我们实现它：

// transform.rs

use image::GenericImageView;

impl ATrousTransform {
	pub fn new(input: &image::DynamicImage, levels: usize) -> Self {
    	let (width, height) = input.dimensions();
        let (width, height) = (width as usize, height as usize);
        
        // 创建一个新的二维数组，用于保存所有输入像素数据的适当大小的每个维度。方法 `zeros` 接受一个“形状”参数，这是一个元组 (行数，列数)。
    	let mut data = Array2::<f32>::zeros((height, width));
        
        // 将图像转换为灰度图像，其中每个像素值都是 `f32` 类型。循环遍历输入图像中的所有像素以及其二维位置。
        for (x, y, pixel) in input.to_luma32f().enumerate_pixels() {
        	// 将像素值放置在我们的数据数组中适当的位置。使用 `[[]]` 语法提供二维索引，例如 `[[行索引，列索引]]`。
        	data[[y as usize, x as usize]] = pixel.0[0];
        }
        
        Self {
        	input: data,
            levels,
            current_level: 0,
            width,
            height
		}
    }
}

这将负责将图像转换为灰度图像，并将像素值转换为 f32。如果你还不知道，对于具有浮点像素值的图像，值始终被归一化。这意味着它们始终在 0 和 1 之间——0 表示黑色，1 表示白色。

迭代器和 À Trous 变换

在继续之前，让我们思考一下算法。我们需要能够生成逐渐增加的尺度的图像，直到达到我们需要的最大级别。

我们希望我们库的使用者能够访问所有这些尺度，并能够对它们进行操作，也可以轻松地在完成后重新组合它们。他们需要能够过滤层以忽略特定尺度的结构，操纵或“映射”它们以更改其特征，在它们上执行操作，甚至存储每个图像（如果需要的话）。

这听起来很像迭代器！迭代器为我们提供了诸如 filter、skip、take、map、for_each 等方法，所有这些方法都是我们在重新组合之前使用的所有方法。

迭代器的一个额外优点是，它允许您在转移到下一个之前完全完成每个层的处理。如果你不确定这是为什么，我建议你阅读更多关于在 Rust 中使用迭代器处理一系列项目的内容。

我们将为我们的 ATrousTransform 类型实现 Iterator 特性，它应该为每次迭代产生一个波纹层作为输出。

我们将首先实现算法的最内部部分，然后再从那里构建出来。因此，我们首先需要一种方法，即在我们的循环中的第一步中，使用缩放函数对输入数据缓冲区进行卷积，并确保相邻像素的间隔为 2n 位置，这是我们循环的第一步。

卷积

在我们做任何其他事情之前，我们需要定义我们的卷积核。创建一个名为 kernel.rs 的新文件，并将其添加到 lib.rs 中，其中包含以下内容：

// kernel.rs

#[derive(Copy, Clone)]
pub struct LinearInterpolationKernel {
    values: [[f32; 3]; 3]
}

impl Default for LinearInterpolationKernel {
    fn default() -> Self {
        Self {
            values: [
                [1. / 16., 1. / 8., 1. / 16.],
                [1. / 8., 1. / 4., 1. / 8.],
                [1. / 16., 1. / 8., 1. / 16.],
            ]
        }
    }
}

我们使用结构体而不是常量数组的数组来定义它，因为我们需要在它上面定义一些与索引处理相关的小方法。我们稍后会回到这一点。

创建另一个文件 convolve.rs。这里将放置处理单个像素的卷积的所有代码。我们将定义一个 Convolution 特性，该特性将定义执行每个像素上的卷积所需的方法。

// convolve.rs

pub trait Convolution {
    fn compute_pixel_index(
        &self,
        distance: usize,
        kernel_index: [isize; 2],
        target_pixel_index: [usize; 2]
    ) -> [usize; 2];
    
    fn compute_convoluted_pixel(
        &self, 
        distance: usize, 
        index: [usize; 2]
    ) -> f32;
}

你可能会问为什么我们需要一个特性而不是一个简单的 impl 块。在本文中，我们只使用灰度图像，但你可能也想扩展它以实现 RGB 或其他颜色模式。

现在，你需要为你的 ATrousTransform 结构实现此特性：

// convolve.rs

impl Convolution for ATrousTransform {
    fn compute_pixel_index(
        &self, 
        distance: usize, 
        kernel_index: [isize; 2], 
        target_pixel_index: [usize; 2]
    ) -> [usize; 2] {
        let [kernel_index_x, kernel_index_y] = kernel_index;
        
        // 通过将它们的相对位置与孔的大小相乘来计算相邻像素的实际距离。
        let x_distance = kernel_index_x * distance as isize;
        let y_distance = kernel_index_y * distance as isize;
        
        let [x, y] = target_pixel_index;

        // 根据当前像素的索引计算 2D 图像中相邻像素的索引。
        let mut x = x as isize + x_distance;
        let mut y = y as isize + y_distance;

        // 如果 x 索引超出边界，将 x 视为最近的边界位置
        if x < 0 {
            x = 0;
        } else if x > self.width as isize - 1 {
            x = self.width as isize - 1;
        }

        // 如果 y 索引超出边界，将 y 视为最近的边界位置
        if y < 0 {
            y = 0;
        } else if y > self.height as isize - 1 {
            y = self.height as isize - 1;
        }

        // 像素的最终 2D 索引。
        [y as usize, x as usize]
    }

    fn compute_convoluted_pixel(
        &self, 
        distance: usize, 
        [x, y]: [usize; 2]
    ) -> f32 {
        // 创建一个新变量来保存当前像素的卷积结果。
        let mut pixels_sum = 0.0;

        let kernel = LinearInterpolationKernel::default();

        // 遍历从中心像素到相对位置的像素的相对位置以执行卷积。
        for kernel_index_x in -1..=1 {
            for kernel_index_y in -1..=1 {
                // 获取映射到内核中当前位置的计算像素位置
                let pixel_index = self.compute_pixel_index(
                	distance,
                    [kernel_index_x, kernel_index_y],
                    [x, y]
                );
                
                // 获取该位置的乘法因子（内核值）来自内核。
                let kernel_value = kernel.value_from_relative_index(
                	kernel_index_x,
                    kernel_index_y
                );

                // 将像素值与内核缩放因子相乘并将其添加到像素总和中。
                pixels_sum += kernel_value * self.input[pixel_index];
            }
        }

        // 从卷积过程中返回计算像素的值。
        pixels_sum
    }
}

我们需要进行计算，以根据内核中的相对位置从中心像素计算出每个像素的位置，同时确保最终像素索引也考虑了“孔大小”。正如你可能注意到的，当计算索引时，你还要处理边界条件。

我建议你花些时间仔细阅读代码和注释。

实现迭代器

终于到了为你的 ATrousTransform 实现 Iterator 特性的时候了：

// transform.rs

impl Iterator for ATrousTransform {
    // 我们的输出也是一个图像，以及每次迭代的当前级别。当前级别是一个 `Option`，用于表示中间层生成后的最终残差层。
    type Item = (Array2::<f32>, Option<usize>);

    fn next(&mut self) -> Option<Self::Item> {
        let pixel_scale = self.current_level;
        self.current_level += 1;

        // 我们已经生成了所有的层。返回 None 以退出迭代器。
        if pixel_scale > self.levels {
            return None;
        }

        // 我们已经生成了所有中间层，返回残差层。
        if pixel_scale == self.levels {
            return Some((self.input.clone(), None))
        }

        let (width, height) = (self.width, self.height);

        // 像素之间的卷积距离（也称为“孔”的大小）。
        let distance = 2_usize.pow(pixel_scale as u32);
        
        // 创建一个新的缓冲区来保存此层的计算数据。
        let mut current_data = Array2::<f32>::zeros((height, width));

        // 遍历 2D 图像中的每个像素位置
        for x in 0..width {
            for y in 0..height {
            	// 将当前层中的当前像素设置为在输入数据中的当前像素上进行卷积的结果。
                current_data[[y, x]] = self.compute_convoluted_pixel(
                	distance, 
                    [x, y]
                );
            }
        }

        // 通过从先前层中减去当前计算的像素来创建当前层
        let final_data = self.input.clone() - &current_data;
        
        // 将输入层设置为等于当前计算的层，以便在下一次迭代中将其用作“上一层”。这也是每个层的残差数据。
        self.input = current_data;

        // 返回当前层数据以及当前级别信息。
        Some((final_data, Some(self.current_level)))
    }
}

我要指出这里有很多优化性能的潜力，但这超出了本文的范围。

最后我们将看看如何将所有这些层组合起来并重构我们的输入图像。

重组

正如我之前所说的，重构通过 A Trous 变换分解的图像就像将所有层相加一样简单。

我们将为此定义一个特性。一旦查看了实现，为什么我们需要一个特性就会变得清楚。

创建一个名为 recompose.rs 的新文件，其中包含以下内容：

// recompose.rs

use image::{DynamicImage, ImageBuffer, Luma};
use ndarray::Array2;

pub trait RecomposableLayers: Iterator<Item = (Array2<f32>, Option<usize>)> {
    fn recompose_into_image(
        self,
        width: usize,
        height: usize,
    ) -> DynamicImage
        where
            Self: Sized,
    {
        // 创建一个结果缓冲区来保存我们输出图像的像素数据。
        let mut result = Array2::<f32>::zeros((height, width));

        // 对于每一层，将层数据添加到结果缓冲区的当前值中。
        for layer in self {
            result += &layer.0;
        }

        // 计算最终数据中的最小和最大像素强度值，以便我们可以执行“重新缩放”，将所有像素值规范化为 0 到 1 的范围，这是浮点 32 图像所期望的。
        let min_pixel = result.iter().copied().reduce(f32::min).unwrap();
        let max_pixel = result.iter().copied().reduce(f32::max).unwrap();

        // 创建一个新的 `ImageBuffer`，这是 `image` 包提供的一种类型，用于作为图像的像素数据缓冲区。在这里，我们正在创建一个新的 `Luma` ImageBuffer，其像素值类型为 `u16`。Luma 只是指灰度。
        let mut result_img: ImageBuffer<Luma<u16>, Vec<u16>> =
            ImageBuffer::new(width as u32, height as u32);
        
        // 预先计算缩放计算的分母，以便我们不必为每次迭代重复这个过程。
        let rescale_ratio = max_pixel - min_pixel;

        // 遍历 `ImageBuffer` 中的所有像素，并根据 `result` 缓冲区中的数据填充它，然后重新缩放值。
        for (x, y, pixel) in result_img.enumerate_pixels_mut() {
            let intensity = result[(y as usize, x as usize)];

            *pixel =
                Luma([((intensity - min_pixel) / rescale_ratio * u16::MAX as f32) as u16]);
        }