|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
||
以下是对散射图像处理各核心步骤的详细说明：
【1. 扣除探测器本底噪声】
– 目的：去除探测器自身固有的电子噪声、暗电流等非信号干扰，确保后续数据仅反映样品的散射信号。
– 原理：探测器在无X射线照射时，各像素会产生一定的基线噪声，需通过预设值或暗场图像扣除。
– 操作细节：
– 对于特定探测器（如线站常用的greateyes探测器`），其本底噪声数值通常为固定值（约3100），可直接在软件中输入该值进行全局扣除；
– 若噪声不稳定，需采集暗场图像（关闭X射线时的探测器读数），逐像素扣除对应噪声值。
– 关键影响：若本底噪声未完全扣除，会导致散射强度基线偏高，掩盖样品的弱信号。
【2. 光斑中心确定与量化】
– 目的：确定X射线入射光斑在探测器上的中心位置，为后续坐标转换和角度校准提供基准点。
– 原理：入射光斑在探测器上形成的散射中心是对称分布的原点，其坐标是散射矢量`q`计算的参考零点。
– 操作细节：
– 光斑中心确定：入射光斑通常表现为中心区域的高强度亮斑，可通过计算图像重心，初步定位；
– 基于对称性精确量化：利用散射图案的对称性（各向同性材料为圆形对称），通过选取强度差异大的行列，计算其对称轴，精确得到中心坐标（x0, y0）。
– 关键影响：中心坐标偏差会导致后续q值计算的角度误差，尤其对高角度（大q）区域影响显著。
function tiff_image_analysis(file_path)% TIFF图像分析与相似度计算脚本% 用法: tiff_image_analysis('path/to/your/image.tif')% 检查输入参数if nargin 1[file_name, file_dir] = uigetfile({'*.tif;*.tiff', 'TIFF Files (*.tif, *.tiff)'});if isequal(file_name, 0) || isequal(file_dir, 0)fprintf('未选择文件n');return;endfile_path = fullfile(file_dir, file_name);endtry% 读取TIF图像img = imread(file_path);% 转换为double类型img_double = double(img);% 输出图像信息[rows, cols] = size(img_double);fprintf('图像尺寸: %d x %dn', rows, cols);% 处理图像数据filtered_img = img_double;filtered_img(filtered_img 0) = NaN;% 计算目标值并在±5的误差范围内查找匹配像素target_value = exp(0.5 * (log(max(max(filtered_img))) + log(mean(mean(filtered_img)))));tolerance = input('请输入误差范围 (默认5): ') || 5; % 固定误差范围±5[m, n] = find(abs(filtered_img - target_value)% 计算中心位置X = [m, n];n_points = length(X(:, 1));y = ones(n_points, 1);b = [rand(1)*1000, rand(1), rand(1)];fun = @(a, X) X(:, 1).^2 + X(:, 2).^2 + a(1)*X(:, 1) + a(2)*X(:, 2) + a(3);[a, ~, ~] = nlinfit(X, y, fun, b);init_center_y = round(-a(1)/2);init_center_x = round(-a(2)/2);fprintf('初始中心位置: (%d, %d)n', init_center_y, init_center_x);% 设置分析参数row_center = input('请输入误差范围 (默认2500): ') ; % 中心行col_center = input('请输入误差范围 (默认3000): ') ; % 中心列half_width = input('请输入误差范围 (默认50): ') ; % 半宽% 计算行方向平均曲线rows_range = row_center - half_width : row_center + half_width;row_average_curve = mean(img_double(rows_range, init_center_x - 300 : init_center_x + 300), 1);% 计算行方向相似度calculate_similarity = @(guess_center) calculate_dtw_similarity(row_average_curve, guess_center);similarity = arrayfun(calculate_similarity, 1:length(row_average_curve));[~, best_row_center] = max(similarity);best_row_center = best_row_center + init_center_x - 300;fprintf('行方向最佳中心位置: %dn', best_row_center);% 计算列方向平均曲线cols_range = col_center - half_width : col_center + half_width;col_average_curve = mean(img_double(init_center_y - 300 : init_center_y + 300, cols_range), 2);% 计算列方向相似度calculate_similarity = @(guess_center) calculate_dtw_similarity(col_average_curve', guess_center);similarity = arrayfun(calculate_similarity, 1:length(col_average_curve));[~, best_col_center] = max(similarity);best_col_center = best_col_center + init_center_y - 300;fprintf('列方向最佳中心位置: %dn', best_col_center);% 绘制行方向曲线figure;x = 1:length(row_average_curve);plot(x, row_average_curve, 'LineWidth', 2);hold on;xline(best_row_center - (init_center_x - 300), 'r--', 'LineWidth', 1.5);title('行方向平均曲线与最佳中心点');xlabel('X轴');ylabel('Y轴');grid on;% 绘制列方向曲线figure;x = 1:length(col_average_curve);plot(x, col_average_curve, 'LineWidth', 2);hold on;xline(best_col_center - (init_center_y - 300), 'r--', 'LineWidth', 1.5);title('列方向平均曲线与最佳中心点');xlabel('X轴');ylabel('Y轴');grid on;% 显示图像figure;imagesc(img_double);colormap(jet);axis image;title('原始图像');% 绘制最佳中心点的整行和整列红线hold on;line([1, cols], [best_col_center, best_col_center], 'Color', 'r', 'LineWidth', 2);line([best_row_center, best_row_center], [1, rows], 'Color', 'r', 'LineWidth', 2);% 以最佳中心点为中心，半径500画绿色圆theta = linspace(0, 2*pi, 100);radius = 500;x_circle = best_row_center + radius * cos(theta);y_circle = best_col_center + radius * sin(theta);plot(x_circle, y_circle, 'g', 'LineWidth', 2);fprintf('分析完成n');catch efprintf('错误: %sn', e.message);endendfunction sim = calculate_dtw_similarity(curve, guess_center)% 使用DTW算法计算相似度% 输入: curve - 待分析的曲线% guess_center - 猜测的中心点位置% 输出: sim - 相似度值% 分割数据data1 = fliplr(curve(1:guess_center));data2 = curve(guess_center:end);% 计算DTW距离dtw_distance = dtw(data1, data2);% 转换为相似度 (距离越小，相似度越大)sim = 1 / (1 + dtw_distance);end
【3. 能量与样品到探测器距离校准】
–目的：校准X射线光子能量（λ）和样品到探测器靶面的距离（SDD），这两个参数是计算散射矢量q的核心物理量。
–原理：散射矢量q = 4πsinθ/λ（θ为半散射角），而sinθ ≈ tanθ = r/SDD（r为探测器上某点到中心的距离），因此q的准确性直接依赖于λ和SDD的精度。
–操作细节：
– 使用聚苯乙烯纳米球体（PS300）等标准样品：其散射图案中存在已知的特征峰（由球间有序排列或单分散粒径决定），通过测量特征峰对应的r值，结合标准样品的理论峰位q_std，反推当前实验的λ和D；
– 例如：PS300的平均粒径已知，其散射峰的q值可通过公式计算，对比实验测得的峰位r，即可校准D和λ（q = 4πr/(λD)）。
– 关键影响：λ或D的误差会导致q值整体偏移，影响后续结构参数（如粒径、间距）的计算准确性。
A = double(imread('250612_PM055741MS730_PS300_10s_410eV_120 mm_1.tif')-3180);cenx = 3480; ceny = 4096-855;phi1=0; phi2=10;rmax=4000;sampletodetector=150;E = 410; lambda = 1239/E;pixelsize=0.015;[results,Fq,FI,Li,Q,I_ps300]=caliLE(A,cenx,ceny,phi1,phi2,rmax,lambda,sampletodetector,pixelsize);disp(Li)function [results,Fq,FI,Li,Q,I_ps300]= caliLE(A,cenx,ceny,phi1,phi2,rmax,lambda,sampletodetector,pixelsize)[q,IR] = cake2qintensity(A,cenx,ceny,phi1,phi2,rmax,lambda,sampletodetector,pixelsize);[Fq,~] = cutuseless(q,IR,0.08,0.3);Q = Fq; R = 150; I_ps300 = (3*(sin(Q*R)-(Q*R).*cos(Q*R))./(Q*R).^3).^2;%此处通常需要替换为实测的PS300的散射曲线[~,ps300peak_indices] = findpeaks(I_ps300, 'MinPeakDistance', 1);peak_ps300 = Q(ps300peak_indices);[L_grid, E_grid] = ndgrid(145:0.5:155, 350:0.5:440);LE = [L_grid(:), E_grid(:)];% 使用arrayfun替代for循环results = arrayfun(@(i) evaluate_LE_pair(i, LE, IR, 4000, 0.015, peak_ps300), 1:length(LE), 'UniformOutput', false);% 提取结果并找到最佳参数[~, top1_indices] = maxk(cellfun(@(x) x.corr_dtw, results), 1);Li = LE(top1_indices,:);% 重新计算并绘制最佳结果[q,IR] = cake2qintensity(A,3480,4096-855,0,10,4000,1239/Li(1,2),Li(1,1),0.015);[Fq,FI] = cutuseless(q,IR,0.08,0.3);figure, loglog(Fq,FI/50000000); hold on; loglog(Q,I_ps300);end% 辅助函数定义function result = evaluate_LE_pair(i, LE, IR, rmax, pixelsize, peak_ps300)q = 4*pi*sin(atan((1:rmax)*pixelsize/LE(i,1))/2)/(1239/LE(i,2));[Fq,FI] = cutuseless(q,IR,0.08,0.3);[~, peak_indices] = findpeaks(FI, 'MinPeakDistance', 40);result.peak_F0 = Fq(peak_indices);result.corr_dtw = vector_correlation(result.peak_F0, peak_ps300, 'dtw');endfunction [q,inten] = cake2qintensity(A,cenx,ceny,phi1,phi2,rmax,lambda,sampletodetector,pixelsize)A(A[lx,la]=size(A);interpolatedA=@(xx,yy) interp2(1:la,1:lx,A,xx,yy);r = 1:rmax;theta=phi1:phi2;theta=theta+180;theta_rad = theta*pi/180;x = r'*cos(theta_rad)+cenx;y = r'*sin(theta_rad)+ceny;Si=interpolatedA(x,y);Si=Si';inten = nanmean(Si, 1);q = 4*pi*sin(atan((1:rmax)*pixelsize/sampletodetector)/2)/lambda;endfunction [Fq,FI] = cutuseless(q,inten,q1,q2)index = find(q >= q1 & qFq = q(index);FI = inten(index);end
【4. 坐标变换（像素→倒易空间）】
– 目的：将探测器上的二维像素坐标（x, y）转换为倒易空间的散射矢量q坐标（单位：nm⁻¹），实现从“空间位置”到“结构特征尺度”的映射。
– 原理：倒易空间中，q的大小反映散射对应的真实空间尺度（d = 2π/q，d为结构特征尺寸），例如q=1 nm⁻¹对应真实空间尺度约6.28 nm。
– 操作细节：
– 首先将像素坐标（x, y）转换为相对中心的距离：r = √[(x – x0)² + (y – y0)²]（单位：像素）；
– 结合像素尺寸（s，单位：μm/像素），将r转换为物理距离：r_phys = r × s（单位：μm）；
– 利用校准后的λ（单位：nm）和D（单位：nm），计算q值：q = 4πr_phys / (λ × D)（单位：nm⁻¹）；
– 最终得到二维散射图中每个像素对应的(q_x, q_y)或极坐标(q, φ)（φ为方位角）。
– 关键影响：坐标变换是连接实验数据与微观结构的桥梁，转换误差会直接导致结构尺寸计算错误。
A = double(imread('250612_PM055741MS730_PS300_10s_410eV_120 mm_1.tif') - 3180);cenx = 3488; % 中心列（注意：MATLAB中X对应列，Y对应行）ceny = 3268; % 中心行phi1 = 175;phi2 = 185;rmax = 4000;sampletodetector = 170; % 样品到探测器距离（单位：mm）E = 410; % 能量（单位：eV）lambda = 1239 / E; % 波长（单位：nm）pixelsize = 15 / 1000; % 像素尺寸（单位：mm）% 计算q空间坐标（单位：nm⁻¹）[psizex, psizey] = size(A); % psizex=行数（Y方向），psizey=列数（X方向）X = (1 - cenx) : (psizey - cenx); % X方向像素偏移（列）qX = 4 * pi * sin(0.5 * atan((X * pixelsize) / sampletodetector)) / lambda;Y = (1 - ceny) : (psizex - ceny); % Y方向像素偏移（行）qY = 4 * pi * sin(0.5 * atan((Y * pixelsize) / sampletodetector)) / lambda;% 优化图像显示figure('Color', 'w'); % 新建白色背景窗口imagesc(qX, qY, A, [0, 5000]); % 按q空间坐标显示axis equal; % 关键：保证qx和qy方向比例一致（按原图比例）colormap('jet'); % 使用美观的配色方案（可替换为'viridis'等）colorbar; % 添加颜色条% 标注单位和标题xlabel('qx (nm^{-1})'); % 标注x轴单位ylabel('qy (nm^{-1})'); % 标注y轴单位set(gca, 'FontName', 'Times New Roman', 'FontSize', 10); % 统一字体和大小xlim([min(qX), max(qX)]); % x轴范围限制为qX的实际范围ylim([min(qY), max(qY)]); % y轴范围限制为qY的实际范围axis tight; % 进一步收紧坐标轴，消除冗余空间
【5. 目标区域转为一维信号（积分处理）】
– 目的：将二维散射图转换为一维I(q)-q曲线（或I(φ)-φ曲线），简化数据形式，突出结构特征。
– 原理：二维散射图包含空间各方向的散射信息，通过积分可提取不同维度的特征（如径向的尺寸分布、方位角的各向异性）。
-操作细节：
– 径向积分（各向同性结构）：
对二维散射图以中心为原点，按不同`q`值进行环形积分（同一`q`值的所有方位角像素强度平均），得到`I(q)-q`曲线。适用于球形颗粒、无规分布的纳米结构等各向同性体系，可反映颗粒尺寸、分布、间距等信息；
– 方位角积分（各向异性结构）：
固定某一q值，沿方位角φ（0~360°）积分，得到I(φ)-φ曲线。适用于纤维、片层、取向纳米结构等各向异性体系，可反映结构的取向程度（如峰的半高宽越小，取向越有序）。
– 关键影响：积分区域的选择（如排除边缘噪声、限定q范围）会影响一维曲线的信噪比，需结合样品特性调整。
[file, path] = uigetfile({'*.tif;*.tiff;*.jpg;*.jpeg;*.png', '图像文件 (*.tif, *.jpg, *.png)'; '*.*', '所有文件 (*.*)'}, '选择图像文件');if isequal(file, 0) || isequal(path, 0)disp('用户取消了选择');return;endimg_path = fullfile(path, file);% 获取其他参数输入offset = input('请输入偏移值: ');cenx = input('请输入中心点x坐标: ');ceny = input('请输入中心点y坐标: ');phi1 = input('请输入起始角度: ');phi2 = input('请输入结束角度: ');rmax = input('请输入最大半径: ');sampletodetector = input('请输入样品到探测器的距离: ');E = input('请输入能量值: ');pixelsize = input('请输入像素尺寸: ')/1000;% 读取图像并处理tryA = double(imread(img_path)) - offset;catch MEerror('图像读取错误: %s', ME.message);end% 计算波长lambda = 1239.842 / E;% 处理数据[q, IR] = cake2qintensity(A, cenx, ceny, phi1, phi2, rmax, lambda, sampletodetector, pixelsize);% 裁剪无用数据q_min = input('请输入q的最小值: ');q_max = input('请输入q的最大值: ');[Fq, FI] = cutuseless(q, IR, q_min, q_max);% 绘制结果figure;semilogy(Fq, FI / 1000000);xlabel('q (Å^{-1})');ylabel('Intensity (arb. units)');title('q-Intensity Plot');grid on;function [q, inten] = cake2qintensity(A, cenx, ceny, phi1, phi2, rmax, lambda, sampletodetector, pixelsize)% 将小于1的值设为NaNA(A 1) = NaN;% 获取图像尺寸[lx, la] = size(A);% 创建插值函数interpolatedA = @(xx, yy) interp2(1:la, 1:lx, A, xx, yy);% 创建半径和角度向量r = 1:rmax;theta = phi1:phi2;theta = theta + 180;theta_rad = theta * pi / 180;% 计算坐标x = r' * cos(theta_rad) + cenx;y = r' * sin(theta_rad) + ceny;% 插值并计算强度Si = interpolatedA(x, y);Si = Si';inten = nanmean(Si, 1);% 计算q值q = 4 * pi * sin(atan((1:rmax) * pixelsize / sampletodetector) / 2) / lambda;endfunction [Fq, FI] = cutuseless(q, inten, q1, q2)% 找到在指定q范围内的数据点index = find(q >= q1 & qFq = q(index);FI = inten(index);end
【6. 背景散射扣除】
– 目的：去除样品以外的散射信号（如样品支撑物、溶剂、空气等），仅保留样品自身的散射信息。
– 原理：背景散射会叠加在样品信号上，尤其在低`q`区域（小角度）可能掩盖样品的弱信号（如大尺寸结构的散射）。
– 操作细节：
– 样品支撑部分的窗口散射**：如石英毛细管、聚合物薄膜等支撑物本身会产生散射，需单独采集支撑物的散射图（无样品时），在相同实验条件下（能量、距离、曝光时间）扣除；
– 溶剂散射：对于溶液样品，需采集纯溶剂（如缓冲液）的散射图，按样品浓度比例扣除（若样品浓度低，溶剂散射可能占主导）；
– 扣除方式：通常采用像素级减法（样品图 – 背景图×校正系数），校正系数需根据曝光时间、样品厚度等参数调整。
– 关键影响：背景扣除不足会导致信号基线偏高，过度扣除则可能引入负强度，需通过对比空白样的散射特征优化。
【7. 数据拼接】
– 目的：扩展散射矢量q的覆盖范围，获取更完整的结构信息（从大尺寸到小尺寸）。
– 原理：单次实验的q范围受探测器尺寸、样品-探测器距离D限制：
– 硬X射线（高能，λ小）穿透能力强，可通过调整D（近距测小q，远距测大q）覆盖0.05~3 nm⁻¹；
– 软X射线（低能，λ大）穿透能力弱，D调整范围有限，单次q范围通常仅为硬X射线的1/5（如0.05~0.5 nm⁻¹），需多次实验拼接。
– 操作细节：
– 对同一样品，在不同D或不同能量下采集多组数据，确保q范围有重叠区域；
– 以重叠区域的强度为基准，对各组数据进行归一化（如按曝光时间、光束强度校准），拼接为一条完整的I(q)-q曲线。
– 关键影响：拼接误差会导致曲线在重叠区域出现跳变，需通过插值或拟合平滑过渡。
本文源自微信公众号：同步辐射共振软X射线散射线站
原文标题：《获取二维散射图像后的数据预处理教程》
原文链接：https://mp.weixin.qq.com/s/6AM4V6DMZKP02S2QYGRsBA
本转载仅出于分享优质测试干货，旨在传递更多观点，并不代表赞同其全部观点或证实其内容的真实性。文章中所包含的图片、音频、视频等素材的版权均归原作者所有。如有侵权请告知删除。
