數位訊號處理程式設計與應用期末專題

專題背景

數字影像相關法（Digital Image Correlation,DIC） 是一種強大且相對低成本的技術，用於以微米尺度解析度測量材料和結構的變形場。該技術的基本原理是通過拍攝材料表面的斑點圖案連續序列圖片，在變形過程中匹配兩幅連續圖像中相似區域，以確定位移場。DIC的應用範圍廣泛，包括材料科學、土木工程、機械工程、生醫工程等領域，本專題將運用MATLAB中的影像處理工具箱，通過DIC方法來測量材料的2D位移場，並進一步處理這些數據以計算並繪製2D的變形測量值。

2d Image cross Correlation的原理:

數字影像相關法是一種非接觸式的測量技術，它基於連續圖像之間的匹配來推斷物體的變形情況。當物體變形時，表面上的斑點圖案會發生位移和變形，DIC通過跟蹤這些斑點圖案的位置變化，來計算物體的位移場。操作方式為通過拍攝物體的連續圖像序列，得到一系列待比對的圖像。然後，DIC算法會將相鄰的兩幅圖像進行比對，通過計算圖像中斑點圖案的相似度，找到相似的圖像區域並確定其位移，並透過對多個圖像對進行連續比對，重建出物體的位移場。對於透明材料，我們還可以利用3D DIC測量完整的三維位移場，例如測量細胞對可變形基板施加的力，3D DIC 通過在不同角度下拍攝物體的圖像，利用立體視覺技術推斷出物體表面各點的三維位移，從而實現對物體變形場的全面測量。

DIC代碼的獲取與應用

本專題中，我們將運用MATLAB中的影像處理工具箱，設計一個2d digital image cross correlation 的模型程式碼，並根據一組圖像確定2D位移場。同時，我們將對這些數據進行進一步的處理和分析，以計算並繪製2D的變形測量值。

dic

資料集選用

資料來源為Ncorr的開源資料Rigid Body Translation & Rotation，包含兩組數據：一組是平移（在x方向上預定位移為-4.25個像素，在y方向上預定位移為-2.75個像素），另一組是旋轉（預定旋轉角度為5度），這兩組數據的參數為：375x375 Image Resolution。

理論模式的建立情形

Cross Correlation理論模式: （Cross-correlation）是一種在訊號處理和統計分析中常用的方法，用於衡量兩個訊號之間的相似度。它通常用於在兩個訊號中找到相同模式的位置或檢測兩個訊號之間的時間延遲，也能用於衡量兩個影像之間的相似度，在一個影像中找到另一個影像的位置或檢測兩個影像之間的空間偏移。

假設有兩個離散時間訊號： - 訊號1：x = [x₁, x₂, x₃, …, xₙ] - 訊號2：y = [y₁, y₂, y₃, …, yₘ]

其中，n 和 m 分別代表兩個訊號的長度。

公式:

Cross-correlation(k) = Σ(x[i] * y[i+k]), for i = 1 to n

這裡的 k 是一個偏移參數，它表示訊號2相對於訊號1的延遲量，通過改變 k 的值，可以計算出不同延遲量下的互相關。

計算步驟如下： 1. 對訊號1和訊號2進行標準化，確保其均值為0並具有相同的變異數。 2. 對於每個 k 值，計算訊號1和訊號2之間的相關性，使用上述公式。 3. 可以使用不同的方法來計算互相關，例如迴圈計算或使用內建函數。 4. 根據互相關的結果，可以找到訊號2在訊號1中的最佳匹配位置或找到兩個訊號之間的時間延遲。

影像表示: 假設有兩個二維影像： - 影像1：I₁（大小為 m × n） - 影像2：I₂（大小為 p × q）

Cross-correlation(k, l) = Σ[ I₁(i, j) * I₂(i+k, j+l) ]

這裡的 k 和 l 是偏移參數，它們表示影像2相對於影像1的水平和垂直偏移量。通過改變 k 和 l 的值，可以計算出不同偏移量下的互相關。

互相關在二維影像上的計算步驟如下： 1. 對影像1和影像2進行前處理，例如灰度化、標準化等，以確保它們具有相似的特徵。 2. 對於每個偏移量 (k, l)，計算影像1和影像2之間的互相關，使用上述公式。 3. 根據互相關的結果，可以找到影像2在影像1中的最佳匹配位置或找到兩個影像之間的空間偏移。

應用範例:

互相關在二維影像處理中有廣泛的應用，包括: 特徵匹配：用於在一個影像中尋找與另一個影像相似的特徵、物體檢測：用於在一個影像中檢測特定物體的位置、影像配準：用於對齊兩個影像，以進行影像融合或影像拼接、影像壓縮或拉伸：基於影像相似性計算壓縮或拉伸。

dic

軟體程式的發展情形

程式解析

導入數據: 首先，使用load函數導入名為rotation_data.mat的數據，該數據包含了參考圖像和當前圖像。

load('rotation_data.mat');
imshow(ref);
imshow(cur);

current

reference

相關分析: 在數位影像相關分析中，我們首先使用normxcorr2函數計算相關係數。

c = normxcorr2(ref, cur);
surf(c);
shading flat;

最大相關值、相關係數及座標

找出最大相關值和座標: 我們可以使用以下代碼找到相關係數矩陣中的最大值和其對應的座標。 最大相關值：0.69493 最大相關座標：(374, 376)

[maxValue, maxIndices] = max(c(:));
[maxRows, maxCols] = find(c == maxValue);
disp(['最大相關值：', num2str(maxValue)]);
disp(['最大相關座標：(', num2str(maxRows), ', ', num2str(maxCols), ')']);

在參考圖像上繪製矩形框: 我們可以在參考圖像上使用drawrectangle函數繪製一個矩形框，該矩形框表示位移場的估計。

imshow(ref);
drawrectangle(gca, 'Position', [xoffSet + 1, yoffSet + 1, size(cur, 2), size(cur, 1)], 'FaceAlpha', 0);

繪製矩形框

網格生成: 生成一個網格，用於在參考圖像和當前圖像上進行位移場計算。

grid_size = 8;
grid_spacing = 25;

x_min = 52;  % 指定 x 座標的最小值
x_max = 329;  % 指定 x 座標的最大值
y_min = 52;  % 指定 y 座標的最小值
y_max = 329;  % 指定 y 座標的最大值

% 生成網格點的 x 座標
ref_grid_x = 1:grid_spacing:size(ref, 2);
ref_grid_x = ref_grid_x(ref_grid_x >= x_min & ref_grid_x + grid_size <= x_max);

% 生成網格點的 y 座標
ref_grid_y = 1:grid_spacing:size(ref, 1);
ref_grid_y = ref_grid_y(ref_grid_y >= y_min & ref_grid_y + grid_size <= y_max);


ref_grid_x = 1:grid_spacing:size(ref, 2);
ref_grid_y = 1:grid_spacing:size(ref, 1);

進行位移場計算: 將在生成的網格點上進行位移場計算。

cur_grid_x = zeros(size(ref_grid_x));
cur_grid_y = zeros(size(ref_grid_y));

for i = 1:numel(ref_grid_x)
    ref_grid_point = [ref_grid_x(i), ref_grid_y(i)];
    ref_grid_rect = ref(ref_grid_point(2):(ref_grid_point(2) + grid_size - 1), ref_grid_point(1):(ref_grid_point(1) + grid_size - 1));
    c = normxcorr2(ref_grid_rect(:,:,1), cur(:,:,1));
    [ypeak, xpeak] = find(c == max(c(:)));
    cur_grid_topleft = [xpeak - size(ref_grid_rect, 2), ypeak - size(ref_grid_rect, 1)];
    ref_grid_center = [ref_grid_x(i) + grid_size / 2, ref_grid_y(i) + grid_size / 2];
    cur_grid_center = [ref_grid_point(1) + cur_grid_topleft(1) + grid_size / 2, ref_grid_point(2) + cur_grid_topleft(2) + grid_size / 2];
    displacement_vector = cur_grid_center - ref_grid_center;
    cur_grid_x(i) = displacement_vector(1);
    cur_grid_y(i) = displacement_vector(2);
end

繪製位移場

繪製位移場: 使用quiver函數繪製位移場。

figure;
quiver(ref_grid_x(:), ref_grid_y(:), cur_grid_x(:), cur_grid_y(:), 'Color', 'g');

x displacement

y displacement

使用smooth2a函數:

smoothed_cur_grid_x = smooth2a(cur_grid_x, 3, 3);
smoothed_cur_grid_y = smooth2a(cur_grid_y, 3, 3);

figure;
contourf(ref_grid_y, ref_grid_x, smoothed_cur_grid_x, 10);
colorbar;

figure;
contourf(ref_grid_y, ref_grid_x, smoothed_cur_grid_y, 10);
colorbar;

x displacement

y displacement

其他變量計算和繪圖: 除了位移場之外，還可以計算其他與變形相關的量，如Jacobian、拉伸、應變等。

[uxx, uxy] = gradient(cur_grid_x, grid_spacing, grid_spacing);
[uyx, uyy] = gradient(cur_grid_y, grid_spacing, grid_spacing);

Jacobian = uxx .* uyy - uxy .* uyx;
left_stretch = sqrt(uxx.^2 + uyx.^2);
right_stretch = sqrt(uyy.^2 + uxy.^2);
strain_invariant = uxx + uyy;
rotation_tensor = atan2(uyx, uxx) * 180 / pi;

繪製位移梯度和其他量的圖形:

figure;
quiver(ref_grid_x(:), -ref_grid_y(:), smoothed_cur_grid_x(:), -smoothed_cur_grid_y(:), 'Color', 'g');
hold on;
quiver(ref_grid_x(:), -ref_grid_y(:), uxx(:), -uyx(:), 'Color', 'r');
quiver(ref_grid_x(:), -ref_grid_y(:), uxy(:), -uyy(:), 'Color', 'b');
title('位移梯度');
legend('位移', 'uxx', 'uyx', 'Location', 'best');
hold off;

位移梯度

繪製變量圖形: 使用contourf函數繪製變量的等高線圖。

figure;
contourf(ref_grid_y, ref_grid_x, Jacobian, 10);
colorbar;
title('Jacobian');

繪製左拉伸、右拉伸、應變和旋轉張量的等高線圖:

figure;
subplot(2, 2, 1);
contourf(ref_grid_y, ref_grid_x, left_stretch, 10);
colorbar;
title('Left Stretch');

subplot(2, 2, 2);
contourf(ref_grid_y, ref_grid_x, right_stretch, 10);
colorbar;
title('Right Stretch');

subplot(2, 2, 3);
contourf(ref_grid_y, ref_grid_x, strain_invariant, 10);
colorbar;
title('Strain Invariant');

figure;
contourf(ref_grid_y, ref_grid_x, rotation_tensor, 10);
colorbar;
title('Rotation Tensor');

變量圖形

旋轉張量

完整程式碼

load('rotation_data.mat');
imshow(ref);
imshow(cur);

% cursubimage = cur(200:209,100:109);
% imshow(cursubimage)

c = normxcorr2(ref, cur);
surf(c)
shading flat

[ypeak, xpeak] = find(c == max(c(:)));
yoffSet = ypeak - size(cur, 1);
xoffSet = xpeak - size(cur, 2);

[maxValue, maxIndices] = max(c(:));
[maxRows, maxCols] = find(c == maxValue);

disp(['最大相關值：', num2str(maxValue)]);
disp(['最大相關座標：(', num2str(maxRows), ', ', num2str(maxCols), ')']);

imshow(ref)
drawrectangle(gca, 'Position', [xoffSet + 1, yoffSet + 1, size(cur, 2), size(cur, 1)], 'FaceAlpha', 0);

ref_point = [100, 200];
rect_width = 10;
rect_height = 10;
ref_rect = ref(ref_point(2):(ref_point(2) + rect_height - 1), ref_point(1):(ref_point(1) + rect_width - 1));
c = normxcorr2(ref_rect, ref);
[ypeak, xpeak] = find(c == max(c(:)));
ref_rect_topleft = [xpeak - size(ref_rect, 2), ypeak - size(ref_rect, 1)];
rectangle = [ref_rect_topleft(1), ref_rect_topleft(2), rect_width, rect_height];
ref = insertShape(ref, 'Rectangle', rectangle, 'LineWidth', 2, 'Color', [1, 0, 0]);
imshow(cur);

grid_size = 8;
grid_spacing = 25;

x_min = 52;  % 指定 x 座標的最小值
x_max = 329;  % 指定 x 座標的最大值
y_min = 52;  % 指定 y 座標的最小值
y_max = 329;  % 指定 y 座標的最大值

% 生成網格點的 x 座標
ref_grid_x = 1:grid_spacing:size(ref, 2);
ref_grid_x = ref_grid_x(ref_grid_x >= x_min & ref_grid_x + grid_size <= x_max);

% 生成網格點的 y 座標
ref_grid_y = 1:grid_spacing:size(ref, 1);
ref_grid_y = ref_grid_y(ref_grid_y >= y_min & ref_grid_y + grid_size <= y_max);


ref_grid_x = 1:grid_spacing:size(ref, 2);
ref_grid_y = 1:grid_spacing:size(ref, 1);


% 檢查網格點是否超出圖像範圍
ref_grid_x = ref_grid_x(ref_grid_x >= x_min & ref_grid_x + grid_size <= x_max);
ref_grid_y = ref_grid_y(ref_grid_y >= y_min & ref_grid_y + grid_size <= y_max);

[ref_grid_x, ref_grid_y] = meshgrid(ref_grid_x, ref_grid_y);
cur_grid_x = zeros(size(ref_grid_x));
cur_grid_y = zeros(size(ref_grid_y));

imshow(ref);
imshow(cur);
hold on;
plot(ref_grid_x(:), ref_grid_y(:), 'ro', 'MarkerSize', 5);

for i = 1:numel(ref_grid_x)
    ref_grid_point = [ref_grid_x(i), ref_grid_y(i)];
    ref_grid_rect = ref(ref_grid_point(2):(ref_grid_point(2) + grid_size - 1), ref_grid_point(1):(ref_grid_point(1) + grid_size - 1));
    c = normxcorr2(ref_grid_rect(:,:,1), cur(:,:,1));
    [ypeak, xpeak] = find(c == max(c(:)));
    cur_grid_topleft = [xpeak - size(ref_grid_rect, 2), ypeak - size(ref_grid_rect, 1)];
    ref_grid_center = [ref_grid_x(i) + grid_size / 2, ref_grid_y(i) + grid_size / 2];
    cur_grid_center = [ref_grid_point(1) + cur_grid_topleft(1) + grid_size / 2, ref_grid_point(2) + cur_grid_topleft(2) + grid_size / 2];
    displacement_vector = cur_grid_center - ref_grid_center;
    cur_grid_x(i) = displacement_vector(1);
    cur_grid_y(i) = displacement_vector(2);
end

plot(cur_grid_x(:), cur_grid_y(:), 'go', 'MarkerSize', 5);
quiver(ref_grid_x(:), ref_grid_y(:), cur_grid_x(:), cur_grid_y(:), 'Color', 'g');
hold off;

figure;
contourf(ref_grid_y, ref_grid_x, cur_grid_x, 10);
colorbar;

figure;
contourf(ref_grid_y, ref_grid_x, cur_grid_y, 10);
colorbar;

smoothed_cur_grid_x = smooth2a(cur_grid_x, 3, 3);
smoothed_cur_grid_y = smooth2a(cur_grid_y, 3, 3);

figure;
contourf(ref_grid_y, ref_grid_x, smoothed_cur_grid_x, 10);
colorbar;

figure;
contourf(ref_grid_y, ref_grid_x, smoothed_cur_grid_y, 10);
colorbar;

[uxx, uxy] = gradient(cur_grid_x, grid_spacing, grid_spacing);
[uyx, uyy] = gradient(cur_grid_y, grid_spacing, grid_spacing);

% 計算其他變量和繪製圖形
Jacobian = uxx .* uyy - uxy .* uyx;
left_stretch = sqrt(uxx.^2 + uyx.^2);
right_stretch = sqrt(uyy.^2 + uxy.^2);
strain_invariant = uxx + uyy;

rotation_tensor = atan2(uyx, uxx) * 180 / pi; % 轉換為角度制


[uxx, uxy] = gradient(smoothed_cur_grid_x, grid_spacing, grid_spacing);
[uyx, uyy] = gradient(smoothed_cur_grid_y, grid_spacing, grid_spacing);



% 繪製位移梯度和其他量的圖形
figure;
quiver(ref_grid_x(:), -ref_grid_y(:), smoothed_cur_grid_x(:), -smoothed_cur_grid_y(:), 'Color', 'g');
hold on;
quiver(ref_grid_x(:), -ref_grid_y(:), uxx(:), -uyx(:), 'Color', 'r');
quiver(ref_grid_x(:), -ref_grid_y(:), uxy(:), -uyy(:), 'Color', 'b');
title('位移梯度');
legend('位移', 'uxx', 'uyx', 'Location', 'best');
hold off;

figure;
subplot(2, 2, 1);
contourf(ref_grid_y, ref_grid_x, Jacobian, 10);
colorbar;
title('Jacobian');

subplot(2, 2, 2);
contourf(ref_grid_y, ref_grid_x, left_stretch, 10);
colorbar;
title('Left Stretch');

subplot(2, 2, 3);
contourf(ref_grid_y, ref_grid_x, right_stretch, 10);
colorbar;
title('Right Stretch');

subplot(2, 2, 4);
contourf(ref_grid_y, ref_grid_x, strain_invariant, 10);
colorbar;
title('Strain Invariant');

figure;
contourf(ref_grid_y, ref_grid_x, rotation_tensor, 10);
colorbar;
title('Rotation Tensor');