基于双目视觉的深度图获取与视点合成

主程序文件

%% 基于双目视觉的深度图获取与视点合成
clear; close all; clc;

fprintf('=== 双目视觉深度图获取与视点合成 ===\n');

%% 系统参数初始化
params = initializeStereoSystem();

%% 加载或生成双目图像
[leftImage, rightImage, trueDisparity] = loadStereoImages(params);

%% 立体标定和校正
[rectifiedLeft, rectifiedRight, stereoParams] = stereoRectification(leftImage, rightImage, params);

%% 深度图计算
[disparityMap, depthMap] = computeDepthMap(rectifiedLeft, rectifiedRight, stereoParams, params);

%% 视点合成
[synthesizedRight, synthesizedLeft] = synthesizeNovelViews(rectifiedLeft, rectifiedRight, depthMap, stereoParams, params);

%% 质量评估和可视化
evaluateResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
               disparityMap, depthMap, trueDisparity, stereoParams, params);

fprintf('\n双目视觉处理完成！\n');

系统参数初始化

function params = initializeStereoSystem()
    % 初始化双目视觉系统参数
    
    fprintf('初始化双目视觉系统参数...\n');
    
    params = struct();
    
    % 相机参数
    params.camera.focalLength = 800;          % 焦距 (像素)
    params.camera.baseline = 0.1;             % 基线距离 (米)
    params.camera.pixelSize = 4.8e-6;         % 像素尺寸 (米)
    params.camera.imageSize = [480, 640];     % 图像尺寸 [高度, 宽度]
    
    % 立体匹配参数
    params.stereo.minDisparity = 0;           % 最小视差
    params.stereo.maxDisparity = 64;          % 最大视差
    params.stereo.windowSize = 7;             % 匹配窗口大小
    params.stereo.smoothness = 100;           % 平滑项权重
    params.stereo.consistencyCheck = true;    % 一致性检查
    
    % 深度图参数
    params.depth.maxDepth = 10;               % 最大深度 (米)
    params.depth.minDepth = 0.5;              % 最小深度 (米)
    
    % 视点合成参数
    params.viewSynthesis.interpolation = 'bilinear';  % 插值方法
    params.viewSynthesis.holeFilling = true;          % 空洞填充
    
    % 显示参数
    params.display.showSteps = true;          % 显示处理步骤
    
    fprintf('系统参数初始化完成\n');
end

双目图像加载

function [leftImage, rightImage, trueDisparity] = loadStereoImages(params)
    % 加载或生成双目图像
    
    fprintf('加载双目图像...\n');
    
    % 选择数据源
    choice = menu('选择图像来源', '使用内置测试图像', '加载自定义图像', '生成虚拟场景');
    
    switch choice
        case 1
            % 使用内置测试图像
            [leftImage, rightImage, trueDisparity] = loadTestImages(params);
        case 2
            % 加载自定义图像
            [leftImage, rightImage, trueDisparity] = loadCustomImages();
        case 3
            % 生成虚拟场景
            [leftImage, rightImage, trueDisparity] = generateVirtualScene(params);
    end
    
    % 转换为灰度图像（如果必要）
    if size(leftImage, 3) == 3
        leftImage = rgb2gray(leftImage);
        rightImage = rgb2gray(rightImage);
    end
    
    % 调整图像尺寸
    leftImage = imresize(leftImage, params.camera.imageSize);
    rightImage = imresize(rightImage, params.camera.imageSize);
    
    fprintf('双目图像加载完成: %dx%d\n', size(leftImage, 2), size(leftImage, 1));
end

function [leftImage, rightImage, trueDisparity] = loadTestImages(params)
    % 加载测试双目图像
    
    % 使用MATLAB内置图像或生成测试图像
    try
        % 尝试加载立体图像对
        leftImage = imread('left.png');
        rightImage = imread('right.png');
        trueDisparity = [];
    catch
        % 生成虚拟立体图像对
        fprintf('生成虚拟立体图像对...\n');
        [leftImage, rightImage, trueDisparity] = generateSyntheticStereoPair(params);
    end
end

function [leftImage, rightImage, trueDisparity] = generateSyntheticStereoPair(params)
    % 生成合成立体图像对
    
    imageSize = params.camera.imageSize;
    height = imageSize(1); width = imageSize(2);
    
    % 创建基础图像
    leftImage = uint8(128 * ones(height, width));
    
    % 添加不同深度的物体
    objects = {
        struct('center', [width/4, height/4], 'radius', 30, 'disparity', 20, 'intensity', 200),  % 近处物体
        struct('center', [3*width/4, height/4], 'radius', 40, 'disparity', 10, 'intensity', 150), % 中等距离
        struct('center', [width/2, 3*height/4], 'radius', 50, 'disparity', 5, 'intensity', 100),  % 远处物体
        struct('center', [width/5, 4*height/5], 'radius', 25, 'disparity', 25, 'intensity', 180), % 很近物体
    };
    
    % 生成左图像
    for i = 1:length(objects)
        obj = objects{i};
        [X, Y] = meshgrid(1:width, 1:height);
        mask = ((X - obj.center(1)).^2 + (Y - obj.center(2)).^2) <= obj.radius^2;
        leftImage(mask) = obj.intensity;
    end
    
    % 生成右图像（根据视差）
    rightImage = leftImage;
    trueDisparity = zeros(height, width);
    
    for i = 1:length(objects)
        obj = objects{i};
        [X, Y] = meshgrid(1:width, 1:height);
        mask = ((X - obj.center(1)).^2 + (Y - obj.center(2)).^2) <= obj.radius^2;
        
        % 在右图像中移动物体
        newX = X - obj.disparity;
        validMask = mask & newX >= 1 & newX <= width;
        
        % 更新右图像
        for y = 1:height
            for x = 1:width
                if validMask(y, x)
                    newXcoord = round(newX(y, x));
                    rightImage(y, newXcoord) = obj.intensity;
                    % 清除原位置（简化处理）
                    if x ~= newXcoord
                        rightImage(y, x) = 128;
                    end
                end
            end
        end
        
        % 记录真实视差
        trueDisparity(mask) = obj.disparity;
    end
    
    % 添加噪声
    leftImage = imnoise(leftImage, 'gaussian', 0, 0.001);
    rightImage = imnoise(rightImage, 'gaussian', 0, 0.001);
end

function [leftImage, rightImage, trueDisparity] = loadCustomImages()
    % 加载自定义图像
    
    [leftFile, leftPath] = uigetfile({'*.jpg;*.png;*.bmp;*.tif', 'Image Files'}, '选择左图像');
    [rightFile, rightPath] = uigetfile({'*.jpg;*.png;*.bmp;*.tif', 'Image Files'}, '选择右图像');
    
    if isequal(leftFile, 0) || isequal(rightFile, 0)
        error('必须选择左右图像文件');
    end
    
    leftImage = imread(fullfile(leftPath, leftFile));
    rightImage = imread(fullfile(rightPath, rightFile));
    trueDisparity = [];
    
    fprintf('自定义图像加载成功\n');
end

function [leftImage, rightImage, trueDisparity] = generateVirtualScene(params)
    % 生成虚拟3D场景
    
    fprintf('生成虚拟3D场景...\n');
    
    imageSize = params.camera.imageSize;
    height = imageSize(1); width = imageSize(2);
    
    % 创建虚拟场景（平面 + 球体）
    leftImage = uint8(zeros(height, width));
    rightImage = uint8(zeros(height, width));
    trueDisparity = zeros(height, width);
    
    % 相机参数
    f = params.camera.focalLength;
    B = params.camera.baseline;
    
    % 生成场景深度
    [X, Y] = meshgrid(1:width, 1:height);
    centerX = width / 2;
    centerY = height / 2;
    
    % 创建深度图（平面 + 球体）
    depthMap = params.depth.maxDepth * ones(height, width);
    
    % 添加倾斜平面
    planeDepth = params.depth.maxDepth - (params.depth.maxDepth - params.depth.minDepth) * ...
                 (Y / height);
    depthMap = min(depthMap, planeDepth);
    
    % 添加球体
    sphereCenter = [centerX, centerY * 0.7];
    sphereRadius = min(width, height) * 0.2;
    sphereDepth = params.depth.minDepth * 1.5;
    
    sphereMask = ((X - sphereCenter(1)).^2 + (Y - sphereCenter(2)).^2) <= sphereRadius^2;
    depthMap(sphereMask) = sphereDepth;
    
    % 添加立方体
    cubeLeft = centerX * 0.6;
    cubeRight = centerX * 0.8;
    cubeTop = centerY * 1.2;
    cubeBottom = centerY * 1.4;
    cubeDepth = params.depth.minDepth * 1.2;
    
    cubeMask = (X >= cubeLeft & X <= cubeRight & Y >= cubeTop & Y <= cubeBottom);
    depthMap(cubeMask) = cubeDepth;
    
    % 生成纹理
    texture = generateTexture(imageSize);
    
    % 生成左右图像
    for y = 1:height
        for x = 1:width
            Z = depthMap(y, x);
            disparity = (f * B) / Z;  % 视差公式
            
            % 左图像像素
            leftImage(y, x) = texture(y, x);
            
            % 右图像对应像素
            x_right = round(x - disparity);
            if x_right >= 1 && x_right <= width
                rightImage(y, x_right) = texture(y, x);
            end
            
            trueDisparity(y, x) = disparity;
        end
    end
    
    % 填充右图像的空白区域
    for y = 1:height
        for x = 1:width
            if rightImage(y, x) == 0
                % 使用邻近像素填充
                x_search = max(1, min(width, x));
                rightImage(y, x) = rightImage(y, x_search);
            end
        end
    end
    
    % 添加噪声
    leftImage = imnoise(leftImage, 'gaussian', 0, 0.005);
    rightImage = imnoise(rightImage, 'gaussian', 0, 0.005);
    
    fprintf('虚拟场景生成完成\n');
end

function texture = generateTexture(imageSize)
    % 生成纹理
    
    height = imageSize(1); width = imageSize(2);
    
    % 使用多种频率的正弦波合成纹理
    texture = zeros(height, width);
    
    frequencies = [0.01, 0.02, 0.05, 0.1];
    amplitudes = [0.3, 0.2, 0.1, 0.05];
    
    for i = 1:length(frequencies)
        [X, Y] = meshgrid(1:width, 1:height);
        pattern = sin(2 * pi * frequencies(i) * X) .* cos(2 * pi * frequencies(i) * Y);
        texture = texture + amplitudes(i) * pattern;
    end
    
    % 归一化到0-255
    texture = (texture - min(texture(:))) / (max(texture(:)) - min(texture(:)));
    texture = uint8(255 * texture);
end

立体标定和校正

function [rectifiedLeft, rectifiedRight, stereoParams] = stereoRectification(leftImage, rightImage, params)
    % 立体标定和图像校正
    
    fprintf('进行立体校正...\n');
    
    % 如果没有标定数据，使用虚拟标定参数
    stereoParams = generateVirtualStereoParams(params);
    
    % 使用立体校正
    [rectifiedLeft, rectifiedRight] = rectifyStereoImages(leftImage, rightImage, stereoParams);
    
    % 可视化校正结果
    if params.display.showSteps
        visualizeRectification(leftImage, rightImage, rectifiedLeft, rectifiedRight);
    end
    
    fprintf('立体校正完成\n');
end

function stereoParams = generateVirtualStereoParams(params)
    % 生成虚拟立体标定参数
    
    imageSize = params.camera.imageSize;
    
    % 相机内参矩阵
    focalLength = params.camera.focalLength;
    principalPoint = [imageSize(2)/2, imageSize(1)/2];
    
    K1 = [focalLength, 0, principalPoint(1);
          0, focalLength, principalPoint(2);
          0, 0, 1];
    
    K2 = K1;
    
    % 畸变系数（假设无畸变）
    distortion1 = [0, 0, 0, 0];
    distortion2 = [0, 0, 0, 0];
    
    % 旋转矩阵（假设相机平行）
    R = eye(3);
    
    % 平移向量（沿X轴）
    T = [params.camera.baseline / params.camera.pixelSize; 0; 0];
    
    % 创建stereoParameters对象
    stereoParams = stereoParameters(K1, distortion1, K2, distortion2, R, T, imageSize);
end

function [rectifiedLeft, rectifiedRight] = rectifyStereoImages(leftImage, rightImage, stereoParams)
    % 立体图像校正
    
    % 使用MATLAB立体校正函数
    [rectifiedLeft, rectifiedRight] = rectifyStereoImagesImpl(leftImage, rightImage, stereoParams);
end

function [rectifiedLeft, rectifiedRight] = rectifyStereoImagesImpl(leftImage, rightImage, stereoParams)
    % 立体图像校正实现
    
    % 简化实现：假设图像已经校正
    % 在实际应用中应该使用完整的校正变换
    
    rectifiedLeft = leftImage;
    rectifiedRight = rightImage;
    
    % 这里可以添加实际的校正代码
    % [rectifiedLeft, rectifiedRight] = rectifyStereoImages(leftImage, rightImage, stereoParams);
end

function visualizeRectification(originalLeft, originalRight, rectifiedLeft, rectifiedRight)
    % 可视化校正结果
    
    figure('Position', [100, 100, 1200, 800]);
    
    subplot(2,2,1);
    imshow(originalLeft);
    title('原始左图像');
    
    subplot(2,2,2);
    imshow(originalRight);
    title('原始右图像');
    
    subplot(2,2,3);
    imshow(rectifiedLeft);
    title('校正后左图像');
    
    subplot(2,2,4);
    imshow(rectifiedRight);
    title('校正后右图像');
    
    % 显示极线对齐
    figure('Position', [100, 100, 1000, 600]);
    
    % 创建合成图像显示极线对齐
    compositeImage = [rectifiedLeft, rectifiedRight];
    imshow(compositeImage);
    title('校正后图像极线对齐检查');
    hold on;
    
    % 绘制几条水平线检查极线对齐
    height = size(rectifiedLeft, 1);
    for y = 1:floor(height/5):height
        plot([1, size(compositeImage, 2)], [y, y], 'r-', 'LineWidth', 1);
    end
    
    sgtitle('立体校正结果', 'FontSize', 14, 'FontWeight', 'bold');
end

深度图计算

function [disparityMap, depthMap] = computeDepthMap(leftImage, rightImage, stereoParams, params)
    % 计算深度图
    
    fprintf('计算深度图...\n');
    
    % 选择立体匹配方法
    method = menu('选择立体匹配方法', '块匹配', '半全局匹配', '两者比较');
    
    switch method
        case 1
            % 块匹配
            disparityMap = blockMatching(leftImage, rightImage, params);
        case 2
            % 半全局匹配
            disparityMap = semiGlobalMatching(leftImage, rightImage, params);
        case 3
            % 两者比较
            disparityMapBM = blockMatching(leftImage, rightImage, params);
            disparityMapSGM = semiGlobalMatching(leftImage, rightImage, params);
            
            % 比较结果
            compareMatchingMethods(leftImage, disparityMapBM, disparityMapSGM, params);
            disparityMap = disparityMapSGM;  % 默认使用SGM
    end
    
    % 后处理视差图
    disparityMap = postProcessDisparity(disparityMap, leftImage, params);
    
    % 转换为深度图
    depthMap = disparityToDepth(disparityMap, stereoParams, params);
    
    % 可视化结果
    if params.display.showSteps
        visualizeDepthComputation(leftImage, disparityMap, depthMap, params);
    end
    
    fprintf('深度图计算完成\n');
end

function disparityMap = blockMatching(leftImage, rightImage, params)
    % 块匹配立体匹配
    
    fprintf('  使用块匹配算法...\n');
    
    height = size(leftImage, 1);
    width = size(leftImage, 2);
    
    minDisp = params.stereo.minDisparity;
    maxDisp = params.stereo.maxDisparity;
    windowSize = params.stereo.windowSize;
    halfWindow = floor(windowSize / 2);
    
    disparityMap = zeros(height, width);
    
    % 为每个像素计算视差
    for y = halfWindow+1:height-halfWindow
        for x = halfWindow+1:width-halfWindow
            bestDisparity = 0;
            bestCost = inf;
            
            % 提取左图像块
            leftPatch = double(leftImage(y-halfWindow:y+halfWindow, x-halfWindow:x+halfWindow));
            
            % 在视差范围内搜索
            for d = minDisp:maxDisp
                x_right = x - d;
                if x_right < halfWindow+1 || x_right > width-halfWindow
                    continue;
                end
                
                % 提取右图像块
                rightPatch = double(rightImage(y-halfWindow:y+halfWindow, x_right-halfWindow:x_right+halfWindow));
                
                % 计算匹配代价（SAD）
                cost = sum(abs(leftPatch(:) - rightPatch(:)));
                
                if cost < bestCost
                    bestCost = cost;
                    bestDisparity = d;
                end
            end
            
            disparityMap(y, x) = bestDisparity;
        end
        
        if mod(y, 50) == 0
            fprintf('    处理进度: %.1f%%\n', 100 * y / height);
        end
    end
    
    % 填充边界
    disparityMap = fillDisparityBoundaries(disparityMap, halfWindow);
end

function disparityMap = semiGlobalMatching(leftImage, rightImage, params)
    % 半全局立体匹配
    
    fprintf('  使用半全局匹配算法...\n');
    
    % 简化实现 - 使用MATLAB的disparity函数
    try
        disparityMap = disparity(leftImage, rightImage, ...
            'Method', 'SemiGlobal', ...
            'DisparityRange', [params.stereo.minDisparity, params.stereo.maxDisparity], ...
            'BlockSize', params.stereo.windowSize, ...
            'ContrastThreshold', 0.5, ...
            'UniquenessThreshold', 15);
    catch
        % 如果MATLAB函数不可用，使用简化实现
        fprintf('    MATLAB disparity函数不可用，使用简化SGM实现\n');
        disparityMap = simplifiedSGM(leftImage, rightImage, params);
    end
end

function disparityMap = simplifiedSGM(leftImage, rightImage, params)
    % 简化的半全局匹配实现
    
    height = size(leftImage, 1);
    width = size(leftImage, 2);
    
    minDisp = params.stereo.minDisparity;
    maxDisp = params.stereo.maxDisparity;
    numDisparities = maxDisp - minDisp + 1;
    
    % 初始化代价立方体
    costVolume = inf(height, width, numDisparities);
    
    % 计算匹配代价（Census变换 + 汉明距离）
    for d = 1:numDisparities
        disparity = minDisp + d - 1;
        
        % 计算Census变换
        leftCensus = censusTransform(leftImage);
        rightCensus = censusTransform(rightImage);
        
        % 计算汉明距离
        for y = 2:height-1
            for x = max(2, disparity+2):width-1
                if x - disparity >= 2
                    % 提取Census描述子
                    leftDesc = leftCensus(y-1:y+1, x-1:x+1);
                    rightDesc = rightCensus(y-1:y+1, x-1-disparity:x+1-disparity);
                    
                    % 计算汉明距离
                    hammingDist = sum(xor(leftDesc(:), rightDesc(:)));
                    costVolume(y, x, d) = hammingDist;
                end
            end
        end
    end
    
    % 路径聚合（简化 - 只考虑4个方向）
    directions = [0, 1; 1, 0; 1, 1; 1, -1];  % 右, 下, 右下, 左下
    aggregatedCost = costVolume;
    
    for dirIdx = 1:size(directions, 1)
        dy = directions(dirIdx, 1);
        dx = directions(dirIdx, 2);
        
        for d = 1:numDisparities
            % 沿路径聚合代价
            for y = (1+abs(dy)):(height-abs(dy))
                for x = (1+abs(dx)):(width-abs(dx))
                    if costVolume(y, x, d) < inf
                        prevCost = aggregatedCost(y-dy, x-dx, d);
                        minPrevCost = min(aggregatedCost(y-dy, x-dx, :));
                        
                        newCost = costVolume(y, x, d) + ...
                                 min(prevCost, minPrevCost + params.stereo.smoothness) - minPrevCost;
                        
                        aggregatedCost(y, x, d) = min(aggregatedCost(y, x, d), newCost);
                    end
                end
            end
        end
    end
    
    % 赢家通吃
    disparityMap = zeros(height, width);
    for y = 1:height
        for x = 1:width
            [~, bestD] = min(aggregatedCost(y, x, :));
            disparityMap(y, x) = minDisp + bestD - 1;
        end
    end
end

function census = censusTransform(image)
    % Census变换
    
    [height, width] = size(image);
    census = false(height, width);
    
    for y = 2:height-1
        for x = 2:width-1
            center = image(y, x);
            neighborhood = image(y-1:y+1, x-1:x+1);
            census(y, x) = neighborhood > center;
        end
    end
end

function disparityMap = postProcessDisparity(disparityMap, leftImage, params)
    % 视差图后处理
    
    fprintf('  视差图后处理...\n');
    
    % 中值滤波去噪
    disparityMap = medfilt2(disparityMap, [3, 3]);
    
    % 一致性检查（如果启用）
    if params.stereo.consistencyCheck
        disparityMap = consistencyCheck(disparityMap, leftImage);
    end
    
    % 空洞填充
    disparityMap = fillDisparityHoles(disparityMap);
end

function disparityMap = consistencyCheck(disparityMap, image)
    % 左右一致性检查
    
    % 简化实现 - 基于纹理的一致性检查
    [gradientMag, ~] = imgradient(image);
    textureMask = gradientMag > 10;  % 纹理区域
    
    % 在低纹理区域进行平滑
    disparityMap(~textureMask) = medfilt2(disparityMap, [5, 5]);
end

function disparityMap = fillDisparityHoles(disparityMap)
    % 填充视差图空洞
    
    % 使用形态学操作填充小空洞
    se = strel('disk', 2);
    disparityMap = imclose(disparityMap, se);
    
    % 使用区域生长填充大空洞
    mask = disparityMap > 0;
    disparityMap = regionfill(disparityMap, ~mask);
end

function depthMap = disparityToDepth(disparityMap, stereoParams, params)
    % 将视差图转换为深度图
    
    focalLength = params.camera.focalLength;
    baseline = params.camera.baseline;
    
    % 避免除零
    disparityMap(disparityMap == 0) = 0.1;
    
    % 深度计算: Z = (f * B) / d
    depthMap = (focalLength * baseline) ./ double(disparityMap);
    
    % 限制深度范围
    depthMap(depthMap > params.depth.maxDepth) = params.depth.maxDepth;
    depthMap(depthMap < params.depth.minDepth) = params.depth.minDepth;
end

function compareMatchingMethods(leftImage, disparityBM, disparitySGM, params)
    % 比较不同匹配方法的结果
    
    figure('Position', [100, 100, 1200, 800]);
    
    subplot(2,3,1);
    imshow(leftImage);
    title('左图像');
    
    subplot(2,3,2);
    imagesc(disparityBM);
    colorbar;
    title('块匹配视差图');
    axis image;
    
    subplot(2,3,3);
    imagesc(disparitySGM);
    colorbar;
    title('半全局匹配视差图');
    axis image;
    
    subplot(2,3,5);
    histogram(disparityBM(:), 50);
    xlabel('视差值');
    ylabel('频次');
    title('块匹配视差分布');
    grid on;
    
    subplot(2,3,6);
    histogram(disparitySGM(:), 50);
    xlabel('视差值');
    ylabel('频次');
    title('半全局匹配视差分布');
    grid on;
    
    sgtitle('立体匹配方法比较', 'FontSize', 14, 'FontWeight', 'bold');
end

function visualizeDepthComputation(leftImage, disparityMap, depthMap, params)
    % 可视化深度计算过程
    
    figure('Position', [100, 100, 1200, 600]);
    
    subplot(2,3,1);
    imshow(leftImage);
    title('参考图像（左）');
    
    subplot(2,3,2);
    imagesc(disparityMap);
    colorbar;
    title('视差图');
    axis image;
    
    subplot(2,3,3);
    imagesc(depthMap);
    colorbar;
    title('深度图');
    axis image;
    
    subplot(2,3,4);
    mesh(depthMap);
    title('深度图3D显示');
    xlabel('X'); ylabel('Y'); zlabel('深度 (m)');
    
    subplot(2,3,5);
    histogram(disparityMap(:), 50);
    xlabel('视差值');
    ylabel('频次');
    title('视差分布');
    grid on;
    
    subplot(2,3,6);
    histogram(depthMap(:), 50);
    xlabel('深度值 (m)');
    ylabel('频次');
    title('深度分布');
    grid on;
    
    sgtitle('深度图计算过程', 'FontSize', 14, 'FontWeight', 'bold');
end

视点合成

function [synthesizedRight, synthesizedLeft] = synthesizeNovelViews(leftImage, rightImage, depthMap, stereoParams, params)
    % 合成新的视点
    
    fprintf('合成新视点...\n');
    
    % 从左图像合成右视点
    synthesizedRight = synthesizeView(leftImage, depthMap, stereoParams, 'left2right', params);
    
    % 从右图像合成左视点
    synthesizedLeft = synthesizeView(rightImage, depthMap, stereoParams, 'right2left', params);
    
    % 可视化合成结果
    if params.display.showSteps
        visualizeViewSynthesis(leftImage, rightImage, synthesizedLeft, synthesizedRight, depthMap, params);
    end
    
    fprintf('视点合成完成\n');
end

function synthesizedView = synthesizeView(referenceImage, depthMap, stereoParams, direction, params)
    % 合成指定方向的视点
    
    [height, width] = size(referenceImage);
    synthesizedView = zeros(height, width, 'like', referenceImage);
    
    focalLength = params.camera.focalLength;
    baseline = params.camera.baseline;
    
    % 创建权重图用于混合
    weightMap = zeros(height, width);
    
    for y = 1:height
        for x = 1:width
            Z = depthMap(y, x);
            
            if Z > params.depth.minDepth && Z < params.depth.maxDepth
                % 计算视差
                disparity = (focalLength * baseline) / Z;
                
                if strcmp(direction, 'left2right')
                    % 从左到右：像素向左移动
                    newX = round(x - disparity);
                else
                    % 从右到左：像素向右移动
                    newX = round(x + disparity);
                end
                
                % 检查边界
                if newX >= 1 && newX <= width
                    % 前向映射
                    synthesizedView(y, newX) = referenceImage(y, x);
                    weightMap(y, newX) = weightMap(y, newX) + 1;
                end
            end
        end
    end
    
    % 处理重叠（平均）
    overlapMask = weightMap > 1;
    synthesizedView(overlapMask) = synthesizedView(overlapMask) ./ weightMap(overlapMask);
    
    % 填充空洞
    if params.viewSynthesis.holeFilling
        synthesizedView = fillHoles(synthesizedView, weightMap, params);
    end
    
    % 转换为原始数据类型
    synthesizedView = cast(synthesizedView, class(referenceImage));
end

function filledImage = fillHoles(image, weightMap, params)
    % 填充图像中的空洞
    
    filledImage = image;
    holeMask = weightMap == 0;
    
    % 使用图像修复填充空洞
    if any(holeMask(:))
        switch params.viewSynthesis.interpolation
            case 'bilinear'
                % 双线性插值
                [X, Y] = meshgrid(1:size(image, 2), 1:size(image, 1));
                validMask = ~holeMask;
                
                if any(validMask(:))
                    filledImage(holeMask) = griddata(X(validMask), Y(validMask), ...
                                                     double(image(validMask)), ...
                                                     X(holeMask), Y(holeMask), 'linear');
                end
                
            case 'inpainting'
                % 使用图像修复算法
                try
                    filledImage = inpaintExemplar(image, holeMask);
                catch
                    % 如果修复失败，使用最近邻填充
                    filledImage = fillHolesNearest(image, holeMask);
                end
        end
    end
end

function filledImage = fillHolesNearest(image, holeMask)
    % 使用最近邻填充空洞
    
    filledImage = image;
    
    % 找到空洞边界
    boundaryMask = bwperim(holeMask);
    
    % 为每个空洞像素分配最近的边界像素值
    [holeY, holeX] = find(holeMask);
    [boundaryY, boundaryX] = find(boundaryMask);
    
    for i = 1:length(holeY)
        % 找到最近的边界像素
        distances = sqrt((boundaryX - holeX(i)).^2 + (boundaryY - holeY(i)).^2);
        [~, nearestIdx] = min(distances);
        filledImage(holeY(i), holeX(i)) = image(boundaryY(nearestIdx), boundaryX(nearestIdx));
    end
end

function visualizeViewSynthesis(originalLeft, originalRight, synthesizedLeft, synthesizedRight, depthMap, params)
    % 可视化视点合成结果
    
    figure('Position', [100, 100, 1400, 900]);
    
    subplot(2,3,1);
    imshow(originalLeft);
    title('原始左图像');
    
    subplot(2,3,2);
    imshow(originalRight);
    title('原始右图像');
    
    subplot(2,3,4);
    imshow(synthesizedLeft);
    title('合成左图像');
    
    subplot(2,3,5);
    imshow(synthesizedRight);
    title('合成右图像');
    
    subplot(2,3,3);
    imagesc(depthMap);
    colorbar;
    title('用于合成的深度图');
    axis image;
    
    subplot(2,3,6);
    % 显示差异图
    diffRight = double(originalRight) - double(synthesizedRight);
    imagesc(abs(diffRight));
    colorbar;
    title('合成右图像差异');
    axis image;
    
    sgtitle('视点合成结果', 'FontSize', 14, 'FontWeight', 'bold');
    
    % 计算质量指标
    calculateQualityMetrics(originalLeft, originalRight, synthesizedLeft, synthesizedRight);
end

function calculateQualityMetrics(originalLeft, originalRight, synthesizedLeft, synthesizedRight)
    % 计算合成图像的质量指标
    
    fprintf('\n=== 合成图像质量评估 ===\n');
    
    % PSNR（峰值信噪比）
    psnrLeft = psnr(synthesizedLeft, originalLeft);
    psnrRight = psnr(synthesizedRight, originalRight);
    
    fprintf('PSNR - 合成左图像: %.2f dB\n', psnrLeft);
    fprintf('PSNR - 合成右图像: %.2f dB\n', psnrRight);
    
    % SSIM（结构相似性）
    ssimLeft = ssim(synthesizedLeft, originalLeft);
    ssimRight = ssim(synthesizedRight, originalRight);
    
    fprintf('SSIM - 合成左图像: %.4f\n', ssimLeft);
    fprintf('SSIM - 合成右图像: %.4f\n', ssimRight);
    
    % MSE（均方误差）
    mseLeft = immse(synthesizedLeft, originalLeft);
    mseRight = immse(synthesizedRight, originalRight);
    
    fprintf('MSE - 合成左图像: %.4f\n', mseLeft);
    fprintf('MSE - 合成右图像: %.4f\n', mseRight);
end

结果评估和可视化

function evaluateResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                        disparityMap, depthMap, trueDisparity, stereoParams, params)
    % 综合结果评估和可视化
    
    fprintf('进行综合结果评估...\n');
    
    % 创建综合结果图
    createComprehensiveResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                              disparityMap, depthMap, trueDisparity, params);
    
    % 深度图精度评估（如果有真实视差）
    if ~isempty(trueDisparity)
        evaluateDepthAccuracy(disparityMap, depthMap, trueDisparity, stereoParams, params);
    end
    
    % 生成报告
    generateReport(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                  disparityMap, depthMap, params);
    
    fprintf('结果评估完成\n');
end

function createComprehensiveResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                                   disparityMap, depthMap, trueDisparity, params)
    % 创建综合结果可视化
    
    figure('Position', [50, 50, 1600, 1000]);
    
    % 原始图像和合成图像对比
    subplot(3,4,1);
    imshow(rectifiedLeft);
    title('校正左图像');
    
    subplot(3,4,2);
    imshow(rectifiedRight);
    title('校正右图像');
    
    subplot(3,4,5);
    imshow(synthesizedLeft);
    title('合成左图像');
    
    subplot(3,4,6);
    imshow(synthesizedRight);
    title('合成右图像');
    
    % 视差图和深度图
    subplot(3,4,3);
    imagesc(disparityMap);
    colorbar;
    title('计算视差图');
    axis image;
    
    subplot(3,4,4);
    imagesc(depthMap);
    colorbar;
    title('计算深度图');
    axis image;
    
    % 3D点云显示
    subplot(3,4,7);
    displayPointCloud(rectifiedLeft, depthMap, params);
    title('3D点云重建');
    
    % 深度图剖面
    subplot(3,4,8);
    plotDepthProfile(depthMap, params);
    title('深度剖面分析');
    
    % 误差分析（如果有真实视差）
    if ~isempty(trueDisparity)
        subplot(3,4,9);
        errorMap = abs(disparityMap - trueDisparity);
        imagesc(errorMap);
        colorbar;
        title('视差误差图');
        axis image;
        
        subplot(3,4,10);
        histogram(errorMap(:), 50);
        xlabel('视差误差');
        ylabel('频次');
        title('视差误差分布');
        grid on;
    end
    
    % 合成质量热图
    subplot(3,4,11);
    qualityMap = calculateQualityMap(rectifiedRight, synthesizedRight);
    imagesc(qualityMap);
    colorbar;
    title('合成质量热图');
    axis image;
    
    % 系统性能总结
    subplot(3,4,12);
    displayPerformanceSummary(rectifiedRight, synthesizedRight, disparityMap, trueDisparity);
    
    sgtitle('双目视觉系统完整结果分析', 'FontSize', 16, 'FontWeight', 'bold');
end

function displayPointCloud(image, depthMap, params)
    % 显示3D点云
    
    [height, width] = size(image);
    
    % 创建点云
    points3D = zeros(height * width, 3);
    colors = zeros(height * width, 3);
    
    focalLength = params.camera.focalLength;
    cx = width / 2;
    cy = height / 2;
    
    idx = 1;
    for y = 1:height
        for x = 1:width
            Z = depthMap(y, x);
            if Z > params.depth.minDepth && Z < params.depth.maxDepth
                % 从图像坐标转换到3D坐标
                X = (x - cx) * Z / focalLength;
                Y = (y - cy) * Z / focalLength;
                
                points3D(idx, :) = [X, Y, Z];
                colors(idx, :) = double(image(y, x)) / 255;
                idx = idx + 1;
            end
        end
    end
    
    % 移除未使用的点
    points3D = points3D(1:idx-1, :);
    colors = colors(1:idx-1, :);
    
    % 显示点云
    scatter3(points3D(:,1), points3D(:,3), -points3D(:,2), 10, colors, 'filled');
    xlabel('X (m)'); ylabel('Z (m)'); zlabel('Y (m)');
    axis equal;
    grid on;
    view(45, 30);
end

function plotDepthProfile(depthMap, params)
    % 绘制深度剖面
    
    [height, width] = size(depthMap);
    
    % 选择中心行
    centerRow = round(height / 2);
    depthProfile = depthMap(centerRow, :);
    
    plot(1:width, depthProfile, 'b-', 'LineWidth', 2);
    xlabel('像素位置');
    ylabel('深度 (m)');
    title(sprintf('第%d行深度剖面', centerRow));
    grid on;
    ylim([params.depth.minDepth, params.depth.maxDepth]);
end

function qualityMap = calculateQualityMap(original, synthesized)
    % 计算合成质量热图
    
    if size(original, 3) == 3
        original = rgb2gray(original);
        synthesized = rgb2gray(synthesized);
    end
    
    % 使用局部SSIM
    qualityMap = zeros(size(original));
    windowSize = 11;
    halfWindow = floor(windowSize / 2);
    
    for y = halfWindow+1:size(original,1)-halfWindow
        for x = halfWindow+1:size(original,2)-halfWindow
            origPatch = original(y-halfWindow:y+halfWindow, x-halfWindow:x+halfWindow);
            synthPatch = synthesized(y-halfWindow:y+halfWindow, x-halfWindow:x+halfWindow);
            
            qualityMap(y, x) = ssim(synthPatch, origPatch);
        end
    end
    
    % 填充边界
    qualityMap = fillQualityBoundaries(qualityMap, halfWindow);
end

function qualityMap = fillQualityBoundaries(qualityMap, halfWindow)
    % 填充质量图的边界
    
    [height, width] = size(qualityMap);
    
    % 填充上边界
    for y = 1:halfWindow
        qualityMap(y, :) = qualityMap(halfWindow+1, :);
    end
    
    % 填充下边界
    for y = height-halfWindow+1:height
        qualityMap(y, :) = qualityMap(height-halfWindow, :);
    end
    
    % 填充左边界
    for x = 1:halfWindow
        qualityMap(:, x) = qualityMap(:, halfWindow+1);
    end
    
    % 填充右边界
    for x = width-halfWindow+1:width
        qualityMap(:, x) = qualityMap(:, width-halfWindow);
    end
end

function displayPerformanceSummary(original, synthesized, disparityMap, trueDisparity)
    % 显示性能总结
    
    % 清空子图
    cla;
    
    % 计算性能指标
    psnrVal = psnr(synthesized, original);
    ssimVal = ssim(synthesized, original);
    mseVal = immse(synthesized, original);
    
    % 视差图统计
    disparityStats = [min(disparityMap(:)), max(disparityMap(:)), mean(disparityMap(:)), std(disparityMap(:))];
    
    text(0.1, 0.9, '系统性能总结:', 'FontSize', 12, 'FontWeight', 'bold');
    text(0.1, 0.8, sprintf('PSNR: %.2f dB', psnrVal), 'FontSize', 10);
    text(0.1, 0.7, sprintf('SSIM: %.4f', ssimVal), 'FontSize', 10);
    text(0.1, 0.6, sprintf('MSE: %.4f', mseVal), 'FontSize', 10);
    text(0.1, 0.5, sprintf('视差范围: [%d, %d]', disparityStats(1), disparityStats(2)), 'FontSize', 10);
    text(0.1, 0.4, sprintf('平均视差: %.2f', disparityStats(3)), 'FontSize', 10);
    
    if ~isempty(trueDisparity)
        disparityError = mean(abs(disparityMap(:) - trueDisparity(:)));
        text(0.1, 0.3, sprintf('视差误差: %.2f', disparityError), 'FontSize', 10);
    end
    
    axis off;
    xlim([0, 1]);
    ylim([0, 1]);
end

function evaluateDepthAccuracy(disparityMap, depthMap, trueDisparity, stereoParams, params)
    % 评估深度图精度
    
    fprintf('\n=== 深度图精度评估 ===\n');
    
    % 计算视差误差
    disparityError = abs(disparityMap - trueDisparity);
    validMask = trueDisparity > 0;
    
    meanDisparityError = mean(disparityError(validMask));
    stdDisparityError = std(disparityError(validMask));
    maxDisparityError = max(disparityError(validMask));
    
    fprintf('视差误差统计:\n');
    fprintf('  均值: %.3f 像素\n', meanDisparityError);
    fprintf('  标准差: %.3f 像素\n', stdDisparityError);
    fprintf('  最大值: %.3f 像素\n', maxDisparityError);
    
    % 计算深度误差
    trueDepth = (params.camera.focalLength * params.camera.baseline) ./ double(trueDisparity);
    depthError = abs(depthMap - trueDepth);
    
    meanDepthError = mean(depthError(validMask));
    stdDepthError = std(depthError(validMask));
    
    fprintf('深度误差统计:\n');
    fprintf('  均值: %.3f 米\n', meanDepthError);
    fprintf('  标准差: %.3f 米\n', stdDepthError);
    
    % 可视化误差分析
    visualizeErrorAnalysis(disparityError, depthError, validMask, params);
end

function visualizeErrorAnalysis(disparityError, depthError, validMask, params)
    % 可视化误差分析
    
    figure('Position', [100, 100, 1200, 600]);
    
    subplot(2,3,1);
    imagesc(disparityError);
    colorbar;
    title('视差误差图');
    axis image;
    
    subplot(2,3,2);
    histogram(disparityError(validMask), 50);
    xlabel('视差误差 (像素)');
    ylabel('频次');
    title('视差误差分布');
    grid on;
    
    subplot(2,3,3);
    cdfplot(disparityError(validMask));
    xlabel('视差误差 (像素)');
    ylabel('累积概率');
    title('视差误差CDF');
    grid on;
    
    subplot(2,3,4);
    imagesc(depthError);
    colorbar;
    title('深度误差图');
    axis image;
    
    subplot(2,3,5);
    histogram(depthError(validMask), 50);
    xlabel('深度误差 (米)');
    ylabel('频次');
    title('深度误差分布');
    grid on;
    
    subplot(2,3,6);
    cdfplot(depthError(validMask));
    xlabel('深度误差 (米)');
    ylabel('累积概率');
    title('深度误差CDF');
    grid on;
    
    sgtitle('深度图精度分析', 'FontSize', 14, 'FontWeight', 'bold');
end

function generateReport(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                       disparityMap, depthMap, params)
    % 生成综合报告
    
    fprintf('\n=== 系统综合报告 ===\n');
    fprintf('图像尺寸: %dx%d\n', size(rectifiedLeft, 2), size(rectifiedLeft, 1));
    fprintf('视差范围: %d - %d 像素\n', params.stereo.minDisparity, params.stereo.maxDisparity);
    fprintf('深度范围: %.1f - %.1f 米\n', params.depth.minDepth, params.depth.maxDepth);
    fprintf('相机基线: %.3f 米\n', params.camera.baseline);
    fprintf('焦距: %.1f 像素\n', params.camera.focalLength);
    
    % 计算处理时间估计
    imageSize = size(rectifiedLeft);
    estimatedTime = estimateProcessingTime(imageSize, params);
    fprintf('估计处理时间: %.2f 秒\n', estimatedTime);
    
    % 保存结果
    saveResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                disparityMap, depthMap, params);
end

function estimatedTime = estimateProcessingTime(imageSize, params)
    % 估计处理时间
    
    numPixels = imageSize(1) * imageSize(2);
    disparityRange = params.stereo.maxDisparity - params.stereo.minDisparity + 1;
    
    % 简化的时间复杂度模型
    baseTime = 1e-6;  % 每个像素的基本处理时间
    estimatedTime = numPixels * disparityRange * baseTime;
end

function saveResults(rectifiedLeft, rectifiedRight, synthesizedLeft, synthesizedRight, ...
                    disparityMap, depthMap, params)
    % 保存处理结果
    
    timestamp = datestr(now, 'yyyymmdd_HHMMSS');
    folderName = sprintf('stereo_results_%s', timestamp);
    
    if ~exist(folderName, 'dir')
        mkdir(folderName);
    end
    
    % 保存图像
    imwrite(rectifiedLeft, fullfile(folderName, 'rectified_left.png'));
    imwrite(rectifiedRight, fullfile(folderName, 'rectified_right.png'));
    imwrite(synthesizedLeft, fullfile(folderName, 'synthesized_left.png'));
    imwrite(synthesizedRight, fullfile(folderName, 'synthesized_right.png'));
    
    % 保存深度数据
    save(fullfile(folderName, 'depth_data.mat'), ...
         'disparityMap', 'depthMap', 'params');
    
    % 保存深度图为图像
    depthImage = mat2gray(depthMap);
    imwrite(depthImage, fullfile(folderName, 'depth_map.png'));
    
    % 保存参数
    writeParametersToFile(params, fullfile(folderName, 'parameters.txt'));
    
    fprintf('结果已保存到文件夹: %s\n', folderName);
end

function writeParametersToFile(params, filename)
    % 将参数写入文件
    
    fid = fopen(filename, 'w');
    if fid == -1
        return;
    end
    
    fprintf(fid, '双目视觉系统参数\n');
    fprintf(fid, '==================\n\n');
    
    fprintf(fid, '相机参数:\n');
    fprintf(fid, '  焦距: %.1f 像素\n', params.camera.focalLength);
    fprintf(fid, '  基线: %.3f 米\n', params.camera.baseline);
    fprintf(fid, '  像素尺寸: %.2e 米\n', params.camera.pixelSize);
    fprintf(fid, '  图像尺寸: %dx%d\n\n', params.camera.imageSize(2), params.camera.imageSize(1));
    
    fprintf(fid, '立体匹配参数:\n');
    fprintf(fid, '  最小视差: %d\n', params.stereo.minDisparity);
    fprintf(fid, '  最大视差: %d\n', params.stereo.maxDisparity);
    fprintf(fid, '  窗口大小: %d\n', params.stereo.windowSize);
    fprintf(fid, '  平滑权重: %d\n\n', params.stereo.smoothness);
    
    fprintf(fid, '深度图参数:\n');
    fprintf(fid, '  最小深度: %.1f 米\n', params.depth.minDepth);
    fprintf(fid, '  最大深度: %.1f 米\n', params.depth.maxDepth);
    
    fclose(fid);
end

参考代码 基于双目的深度图获取,根据图像合成深度图左右视点 www.youwenfan.com/contentcni/65569.html

系统特点

1. 完整流程

• 图像获取: 支持测试图像、自定义图像和虚拟场景
• 立体校正: 极线对齐和图像校正
• 深度计算: 多种立体匹配算法
• 视点合成: 基于深度的新视点生成

2. 核心技术

• 立体匹配: 块匹配和半全局匹配
• 深度计算: 精确的视差-深度转换
• 视点合成: 前向映射和空洞填充
• 质量评估: PSNR、SSIM等客观指标

3. 高级功能

• 3D点云重建: 从深度图生成3D模型
• 误差分析: 完整的精度评估
• 性能优化: 处理时间估计和优化建议
• 结果导出: 完整的数据保存和报告生成