Julia中散点数据插值生成平滑曲线及求导问题求助

首先修正你原始代码中的全局变量bug——color_means作为全局变量会在每次调用calculate_std_over_windows_gray时累积数据，导致后续标准差计算错误。修正后的代码如下：

using Images, Plots, Colors, Statistics, TestImages

img10 = testimage("water_512.tiff")

function average_color_gray(window)
    return mean(window)
end

function color_std_gray(values)
    return std(values)
end

# 将color_means移到函数内部，避免全局累积
function calculate_std_over_windows_gray(img, window_size)
    color_means = []
    h, w = size(img)
    for i in 1:window_size:h-window_size+1
        for j in 1:window_size:w-window_size+1
            window = img[i:i+window_size-1, j:j+window_size-1]
            mean_color = average_color_gray(window)
            push!(color_means, mean_color)
        end
    end
    std_dev = color_std_gray(color_means)
    return std_dev
end

window_sizes = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]

std_devs = []
for size in window_sizes
    if size < min(size(img10)...)
        std_dev = calculate_std_over_windows_gray(img10, size)
        push!(std_devs, std_dev)
    end
end

std_devss = std_devs[.!isnan.(std_devs)]
x = 1:length(std_devss)  # 对应你的1:9
ss = log.(1 ./ std_devss)

方法1：用Interpolations.jl做三次样条插值（平滑+求导） #

你之前报错大概率是因为数据类型不匹配（比如用整数x构造插值）或版本差异。以下是兼容最新版本的代码：

using Interpolations

# 将x转换为浮点数数组，避免整数类型导致的方法匹配问题
x_float = Float64.(x)
y_float = Float64.(ss)

# 构造三次样条插值器（平滑）
itp = interpolate((x_float,), y_float, Gridded(Interpolations.Cubic(Line(OnGrid()))))
# 生成高密度x点用于绘制平滑曲线
x_smooth = range(first(x_float), last(x_float), length=100)
y_smooth = itp.(x_smooth)

# 求导：Interpolations.jl支持直接生成导数插值器
itp_deriv = derivative(itp)
y_deriv = itp_deriv.(x_smooth)

# 绘图展示
plot(x_smooth, y_smooth, label="平滑曲线", linewidth=2)
scatter!(x, ss, label="原始散点", markersize=5)
plot!(x_smooth, y_deriv, label="导数曲线", linewidth=2, linestyle=:dash, color=:red)

方法2：用Loess.jl做局部回归平滑（适合S型非线性曲线） #

如果插值过度拟合噪声，局部回归（Loess）是更优的平滑选择：

using Loess

# 构造Loess模型，span参数控制平滑程度（0-1，值越大越平滑）
model = loess(x_float, y_float, span=0.7)
# 生成平滑曲线
y_smooth_loess = predict(model, x_smooth)

# 求导：Loess模型的导数需要手动计算（或用数值差分）
# 数值差分求导
y_deriv_loess = diff(y_smooth_loess) ./ diff(x_smooth)
# 补全长度使其与x_smooth匹配
push!(y_deriv_loess, y_deriv_loess[end])

# 绘图
plot(x_smooth, y_smooth_loess, label="Loess平滑曲线", linewidth=2)
scatter!(x, ss, label="原始散点", markersize=5)
plot!(x_smooth, y_deriv_loess, label="导数曲线", linewidth=2, linestyle=:dash, color=:green)

方法3：用Dierckx.jl做B样条插值（专业平滑曲线） #

如果之前用Dierckx报错，可能是未正确传递参数，以下是正确用法：

using Dierckx

# 构造B样条插值器，s参数控制平滑度（0为插值所有点，越大越平滑）
spl = Spline1D(x_float, y_float, s=0.5)
y_smooth_dierckx = spl(x_smooth)

# 求导：直接调用导数方法
y_deriv_dierckx = derivative(spl, x_smooth)

# 绘图
plot(x_smooth, y_smooth_dierckx, label="B样条平滑曲线", linewidth=2)
scatter!(x, ss, label="原始散点", markersize=5)
plot!(x_smooth, y_deriv_dierckx, label="导数曲线", linewidth=2, linestyle=:dash, color=:blue)

常见报错排查 #

• 方法不匹配：确保x和y都是Float64类型，避免整数数组直接传入插值函数。
• 函数不存在：检查包的版本，比如Interpolations.jl在v0.14+版本中构造插值器的语法有变化，建议使用] add Interpolations@latest更新到最新版。
• ScatteredInterpolation.jl不适用：你的数据是规则网格点（x是连续的1:9），不是散乱点，因此不需要用散乱点插值包，用上述规则插值或平滑方法即可。

内容的提问来源于stack exchange，提问作者Maxim Tiburziano