sensor fusion for attitude estimation.

[xiahaa]

I have uploaded the sensor fusion algorithm. This implementation is based on the following paper: Mahony, Robert, Tarek Hamel, and Jean-Michel Pflimlin. "Nonlinear complementary filters on the special orthogonal group." IEEE Transactions on automatic control 53.5 (2008): 1203-1217.

the special orthogonal group is nothing but the rotation matrix (I hope you still remember what we learned in the first course: representation of rotation and pose).

If you have some interesting topices you would like to learn on Friday, just write to me and we will see whether or not we can prepare something in Friday. Otherwise, there will be no new things this week.

[xiahaa]

I hope you can also finish this. video1 video2.

I will later add a docs to help you.

[hannaLyu]

Okey. and now I have some problems about 3 picture stitching. How can I get matching pairs? cause when I stitch 2 picture, I use hamming distance to get data about two figure.

[xiahaa]

Okey. and now I have some problems about 3 picture stitching. How can I get matching pairs? cause when I stitch 2 picture, I use hamming distance to get data about two figure.

One way is that you stich two images firstly, and then use the stiched image and image 3 to do a further stich. This is the easiest way. Another way is this tutorial.

[xiahaa]

cylindrical projection. page 8.

More about image stiching

So what you need to do is a mapping function that can transform a pixel to its corresponding pixel on a cylindrical plane.

[hannaLyu]

Now i try to calculate matchingpairs for each picture,But why i have some inf value in workspace of matching pairs? https://github.com/xiaohanlyu/image-stitching/blob/master/src/stitching/multistit.m

[hannaLyu]

% Compute the output limits for each transform for i = 1:numel(tforms) [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(i,2)], [1 imageSize(i,1)]); end

what's meaning of outputlimits function?cause i do not find it by using help in matlab

[xiahaa]

% Compute the output limits for each transform for i = 1:numel(tforms) [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(i,2)], [1 imageSize(i,1)]); end

what's meaning of outputlimits function?cause i do not find it by using help in matlab

Have you installed the computer vision toolbox?

[xiahaa]

Now i try to calculate matchingpairs for each picture,But why i have some inf value in workspace of matching pairs? https://github.com/xiaohanlyu/image-stitching/blob/master/src/stitching/multistit.m

            [val, id] = sort(despxor,'ascend');
            if val(1) < 0.7 * val(2)
                matching_pairs(i,:) = [i, id(1)];
            else
                matching_pairs(i,:) = [i, inf];
            end

When the best match and second best match are too similar (based on matching score, if the best one is not 0.7 smaller than the second best one), we will treat this feature as feature that cannot be matched and assign an inf.

[xiahaa]

For multiple images stiching: see ex9. For IMU aided image stiching: see ex10. For open project on visual/inertial sensor fusion for image stiching: see ex11.

[xiahaa]

% Compute the output limits for each transform for i = 1:numel(tforms) [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(i,2)], [1 imageSize(i,1)]); end

what's meaning of outputlimits function?cause i do not find it by using help in matlab

我也找不到。 你可以用imtransform的第二个第三个输出参数来替代。

[hannaLyu]

现在我遇到了一个问题。使用Imtransform之前要计算tform，但是计算tfrom的时候要求变换矩阵最后一列是0，0，1，用什么办法可以转换？

[xiahaa]

H./H(3,3) 你要选择projective2d T = maketform('projective',H');

[hannaLyu]

明天是不是没有lecture了？

[hannaLyu]

所以在图片拼接的过程中，进行的是投影变换而不是放射变换？因为一张图片要转移到另一张图所在的坐标系中？我这样理解对吗

[xiahaa]

所以在图片拼接的过程中，进行的是投影变换而不是放射变换？因为一张图片要转移到另一张图所在的坐标系中？我这样理解对吗

前面对，后面不对。2D变换有欧式，相似，仿射，投影，都可以做坐标变换。但是变换能力的强弱和变换的自由对有关，变换自由度越高，变换就越多。 2D欧式变换2个自由对，相似变换3个，仿射是6个，投影是8个。

[xiahaa]

明天是不是没有lecture了？

没有。 有问题就答疑，没有就过了。 下周开始benchmark，做得快的话，看看能不能找个相关的题目结合一下，带你发一篇conference paper。

[hannaLyu]

好的，那我下午就不去了，尽力在明天把三图拼接搞定

[xiahaa]

好的，那我下午就不去了，尽力在明天把三图拼接搞定

👌

[hannaLyu]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

[xiahaa]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

% check answer
I1 = WarpAndBlend(seq.H{1},seq.imgs{1},seq.imgs{2});
subplot(2,2,1);
imshow(I1);
subplot(2,2,2);
I2 = WarpAndBlend(seq.H{2},seq.imgs{1},seq.imgs{3});
imshow(I2);
subplot(2,2,3);
I3 = WarpAndBlend(seq.H{3},seq.imgs{2},seq.imgs{3});
imshow(I3);

这三个两图拼接结果对不对？

[xiahaa]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

H没错，Matlab里面x沿着column方向，算imtransform的时候按照常规定义，x是沿着图像的width方向，那么是不是你要做点啥才行？

2. 哪个inverse函数。

[hannaLyu]

不好意思打错了，是invert invtform = invert(tform) returns the inverse of the geometric transformation tform.

[xiahaa]

不好意思打错了，是invert invtform = invert(tform) returns the inverse of the geometric transformation tform.

哥，maketform(inv(H))不就结了么。/(ㄒoㄒ)/~~

[hannaLyu]

不好意思打错了，是invert invtform = invert(tform) returns the inverse of the geometric transformation tform.

哥，maketform(inv(H))不就结了么。/(ㄒoㄒ)/~~

哈哈我晕了

[hannaLyu]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

H没错，Matlab里面x沿着column方向，算imtransform的时候按照常规定义，x是沿着图像的width方向，那么是不是你要做点啥才行？

1. 哪个inverse函数。

我没理解你的意思。x指的是什么？

[xiahaa]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

H没错，Matlab里面x沿着column方向，算imtransform的时候按照常规定义，x是沿着图像的width方向，那么是不是你要做点啥才行？

1. 哪个inverse函数。

我没理解你的意思。x指的是什么？

MATLAB里面矩阵，（1，2），1是沿着那个方向算的？教过你一个名词叫column-major。 图像变换的时候，x是row-major。 不能再提示咯 😯

[hannaLyu]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

H没错，Matlab里面x沿着column方向，算imtransform的时候按照常规定义，x是沿着图像的width方向，那么是不是你要做点啥才行？

1. 哪个inverse函数。

我没理解你的意思。x指的是什么？

MATLAB里面矩阵，（1，2），1是沿着那个方向算的？教过你一个名词叫column-major。 图像变换的时候，x是row-major。 不能再提示咯

明白了

[xiahaa]

1.我上传了未完成的三图拼接。计算H矩阵检查答案的时候发现出来的图片不对，或许是H矩阵算错了，但我找不出来bug。 2.有没有哪个函数可以代替Inverse的？

H没错，Matlab里面x沿着column方向，算imtransform的时候按照常规定义，x是沿着图像的width方向，那么是不是你要做点啥才行？

1. 哪个inverse函数。

我没理解你的意思。x指的是什么？

MATLAB里面矩阵，（1，2），1是沿着那个方向算的？教过你一个名词叫column-major。 图像变换的时候，x是row-major。 不能再提示咯

明白了

[hannaLyu]

我没看懂你写的blend..直接用会报错，我尝试用tutorial 里面的方式来做，错误提示我的数据类型不对。这种情况是哪里出了问题？ 另外，我的recompute H部分算的对吗？

[xiahaa]

function Out = WarpAndBlend(H,ImL,ImR)
    % 生成一个tform
    T = maketform('projective',H');
    % 用生成的tform变换IML，但是只用输出的x, y range。
    [~,XData,YData]=imtransform(ImL,T,'FillValues',zeros(size(ImL,3),1));

   % 上下界
    XData=[floor(XData(1)) ceil(XData(2))];
    YData=[floor(YData(1)) ceil(YData(2))];
    % 和图2比一下，得到最大的bounding box
    nX = [min(0,XData(1)) max(size(ImR,2),XData(2))];
    nY = [min(0,YData(1)) max(size(ImR,1),YData(2))];
    % 变换，此时指定x y的范围，这样就可以把图变道bounding box的尺寸
    [ImLH]=imtransform(ImL,T,'XData',nX,'YData',nY,'FillValues',zeros(size(ImL,3),1));
    % 变换，此时指定x y的范围，这样就可以把图变道bounding box的尺寸 
  [ImRH]=imtransform(ImR,maketform('affine',eye(3)),'XData',nX,'YData',nY,'FillValues',zeros(size(ImL,3),1));

    % mask这部分没啥用了，旧的代码
    Mask1=ones(size(ImL,1),size(ImL,2));
    Mask2=ones(size(ImR,1),size(ImR,2));
    Mask1=imtransform(Mask1,T,'XData',nX,'YData',nY,'FillValues',0)>0;
    Mask2=imtransform(Mask2,maketform('affine',eye(3)),'XData',nX,'YData',nY,'FillValues',0)>0;
    Maskavg = Mask1 & Mask2;
    Out = uint8(zeros(size(ImLH)));

    % 生成一个距离图，图像boarder上是1，然后调用bwdist，会把这个距离扩散开来，比如第二行就会是2，你可以imshow一下这个dist1
    dist1 = zeros(size(ImL,1),size(ImL,2));%dist1(round(size(ImL,1)*0.5),round(size(ImL,2)*0.5)) = 1;
    dist1(1,:) = 1;dist1(end,:) = 1;dist1(:,1) = 1;dist1(:,end) = 1;
    dist1 = bwdist(dist1,'euclidean');%maxdist = max(dist1(:));dist1 = (maxdist+1) - dist1;

    dist2 = zeros(size(ImR,1),size(ImR,2));%dist2(round(size(ImR,1)*0.5),round(size(ImR,2)*0.5)) = 1;
    dist2(1,:) = 1;dist2(end,:) = 1;dist2(:,1) = 1;dist2(:,end) = 1;
    dist2 = bwdist(dist2,'euclidean');%maxdist = max(dist2(:));dist2 = (maxdist+1) - dist2;
   % 把dist变到bounding box
    dist1t=imtransform(dist1,T,'XData',nX,'YData',nY,'FillValues',0);
    dist2t=imtransform(dist2,maketform('affine',eye(3)),'XData',nX,'YData',nY,'FillValues',0);

%     blender = vision.AlphaBlender('Operation', 'Binary mask', ...
%         'MaskSource', 'Input port');
% 
%     Out = step(blender, ImLH, ImRH, Mask2);

    % 根据距离加权两个图的intensity，原则就是离中心越近，weight越高。
    if size(ImL,3) == 1
        Out = combine(ImLH, ImRH, dist1t, dist2t);
    else
        Out(:,:,1) = combine(ImLH(:,:,1), ImRH(:,:,1), dist1t, dist2t);
        Out(:,:,2) = combine(ImLH(:,:,2), ImRH(:,:,2), dist1t, dist2t);
        Out(:,:,3) = combine(ImLH(:,:,3), ImRH(:,:,3), dist1t, dist2t);
    end
end

function imblend = poissonBlend(source,target,mask)
    m = size(source,1);n = size(source,2);
    N = m*n;
    A = spdiags([-4.*ones(N,1) ones(N,1) ones(N,1) ones(N,1) ones(N,1)],[0,1,-1,m,-m],N,N);
    ii = m:m:N-m;
    jj = ii+1;
    A((ii-1)*N+jj) = 0;
    ii = m:m:N-m;
    A((ii.*N+ii)) = 0;
    % gradient of source
    ds = A * source(:);
    % gradient of target
    dt = A * target(:);

end

function im12 = combine(im1, im2, dist1,dist2)
    weight1 = dist1./(max(dist1(:)));
    weight2 = dist2./(max(dist2(:)));
    weights = weight1 + weight2;

    im12 = weight1.*im2double(im1) + weight2.*im2double(im2);
    im12 = im12 ./ weights;

    im12 = uint8(255/(max(im12(:))-min(im12(:))).*(im12-min(im12(:))));
end

[xiahaa]

你这一步就不太对了 ··· % compute out put limits for each transform seq.tform=cell(numImages,1); seq.T=cell(numImages,1); for i = 1:numel(seq.H) H=seq.H{i}; img=seq.imgs{i}; tform = maketform('projective',H'); [B,xlim(i,:),ylim(i,:)] = imtransform(img,tform); tforms(i)=tform; end ··· 你这里的H分别是1到2,1到3,2到3. 你这样做算不出来limit的，你都做成统一的基。 1到2,1到3这种。

[hannaLyu]

也就是说其实只需要两个tform？

[xiahaa]

也就是说其实只需要两个tform？

是的。 但是感觉你并没有整明白。 你可以先假设所有图都用1来当base，那你要找的是i到1的H，然后算出limit，然后找到中心，然后在选中心图做base，重新算H。 举例： 1当base，你要找2->1, 3->1, .... 你现在是图少，且有匹配所以你可以直接算3->1，假如5->1没有匹配，直接算，算不出来H，怎么办？算5->4->1。

[hannaLyu]

我最终还是向命运屈服。。选择了方案一。 你周五下午有空吗？我想找你再问一下Opion2,就是寻找center image的这个流程，感觉我其实不太明白自己要做什么。而且我下周开始的工作肯定也还有问题，我攒着一块问你

[xiahaa]

我最终还是向命运屈服。。选择了方案一。 你周五下午有空吗？我想找你再问一下Opion2,就是寻找center image的这个流程，感觉我其实不太明白自己要做什么。而且我下周开始的工作肯定也还有问题，我攒着一块问你

1做出来了么？

这就跟拼积木一样，你选了一个做起点，是不是把后面的都往第一个靠。每靠一个，图是不是就大了一些。全拼完是不是就是一个大图？那么离最后大图的中心位置最近的那个积木是不是就能找到？ 找到之后，把它当成起点，再拼一下就好了。

[hannaLyu]

我上传了1，就是把同样的工作重复了两遍。 我现在卡住的地方是，首先我找到了center image，然后通过重新算的H矩阵把每一个图和center image拼起来，接下来要怎么把他们合在一起呢？用blending吗？其实我不太明白blending 的原理

[xiahaa]

我上传了1，就是把同样的工作重复了两遍。 我现在卡住的地方是，首先我找到了center image，然后通过重新算的H矩阵把每一个图和center image拼起来，接下来要怎么把他们合在一起呢？用blending吗？其实我不太明白blending 的原理

还是用warpandblend。 warp就是变换，把图转一转，对上。 blend就是A+B，直接加的话，会有锯齿啊，不连续，那么blend就通过算一些量来加权的加，这样锯齿和不连续就不明显。blend不影响任何拼接，只影响最后成图的视觉效果。

[hannaLyu]

让我明天再试一下

[hannaLyu]

我好像搞明白H矩阵怎么算了。但还是那个问题，我不会拼三张图。。 用warpandblend函数只能拼接已知H变换的两张图，对吗？可我没有拼接完之后的H变换矩阵

[hannaLyu]

有没有四张图的资源？我想试一下。我尝试用自己拍的，但都是海，效果不是很好

[xiahaa]

用warpandblend函数只能拼接已知H变换的两张图，对吗？可我没有拼接完之后的H变换矩阵

咋没有？好好想想。 if 1 is the base, you need H_1^2, H_1^3; now you find 2 should be the base, then you need H_2^1, H_2^3 so H_2^1 = inv(H_1^2) H_2^3 = H_1^3*H_2^1.

[xiahaa]

有没有四张图的资源？我想试一下。我尝试用自己拍的，但都是海，效果不是很好

海没有足够的可区别feature啊。 https://github.com/yihui-he/panorama/tree/master/imgs 这里面有多图。

[xiahaa]

我好像搞明白H矩阵怎么算了。但还是那个问题，我不会拼三张图。。 用warpandblend函数只能拼接已知H变换的两张图，对吗？可我没有拼接完之后的H变换矩阵

知道你啥意思了。 你选了1作为base，把2往1拼，拼完了之后得到12 此时12和3拼，你不知道3到12的H，是这意思吧？ 3到12的H就是3到1的H。因为1是base。

你这样，你这个问题本质就是不会拼三张图，和怎么算H并没啥太大关系。 你跳过计算center image的步骤，就用1当做center image，这样你不用重新算H，就直接拼3张图，先拼1+2=12，然后用12和3拼，用3到1的H当做H。

[hannaLyu]

我刚开始就是这么想的，但出来的结果很奇怪。。我把修改过的threeimg上传了

[xiahaa]

我刚开始就是这么想的，但出来的结果很奇怪。。我把修改过的threeimg上传了

warpandblend的问题，不好意思。不太支持多图。 我改了一个临时的。 你看看吧。 你把warpImage1 2 3 显示出来，这三个就是最后变换的图，要把这三个blend起来就可以了。

[hannaLyu]

warpimage只能显示一半？ 按你说的blend了一下，但是只有一个角，所以判断不出来

[xiahaa]

warpimage只能显示一半？ 按你说的blend了一下，但是只有一个角，所以判断不出来

blend的不对吧？把结果贴出来看看。

[xiahaa]

你提交的代码的做法不太对。 warpimage 1,2,3是已经做过homography变换的了，理论上你把三张加起来就是结果。 更好的做法你可以参照这个：

for i = 1:numImages

    I = readimage(buildingScene, i);

    % Transform I into the panorama.
    warpedImage = imtransform(I, tforms(i), 'XData',xLimits,'YData',yLimits,'FillValues',zeros(size(I,3),1));

%%%%%%%%%%%%%%%%%%%%% Alpha blend %%%%%%%%%%%%%%%%%%%%%%%%%
%     % Generate a binary mask.
%     mask = imtransform(ones(size(I,1),size(I,2)), tforms(i), 'XData',xLimits,'YData',yLimits,'FillValues',0)>0;
%     % Overlay the warpedImage onto the panorama.
%     panorama = step(blender, panorama, warpedImage, mask);
%%%%%%%%%%%%%%%%%%%%% sinple blend using distance %%%%%%%%%%%%%%%%%%%%%%%%%
    dist1 = zeros(size(I,1),size(I,2));%dist1(round(size(I,1)*0.5),round(size(I,2)*0.5)) = 1;
    dist1(1:margin,:) = 1;dist1(end-margin:end,:) = 1;dist1(:,1:margin) = 1;dist1(:,end-margin:end) = 1;
    dist1 = bwdist(dist1,'euclidean');%maxdist = max(dist1(:));dist1 = (maxdist+1) - dist1;
    dist1t=imtransform(dist1,tforms(i),'XData',xLimits,'YData',yLimits,'FillValues',0) + 1e-3;
    [panorama,weights] = simpleBlend(warpedImage, panorama, weights, dist1t);

%%%%%%%%%%%%%%%%%%%%% laplacian blend using distance %%%%%%%%%%%%%%%%%%%%%%%%%
%     mask = ones(size(I,1),size(I,2))*1;
%     mask = imtransform(mask, tforms(i), 'XData',xLimits,'YData',yLimits,'FillValues',0);
%     mask((rgb2gray(panorama) ~= 0)&mask) = 1;
%     panorama = laplacianBlend(warpedImage, panorama, mask);
end
weights = weights + 1e-3;
panorama = panorama ./ weights;
figure
imshow(panorama)

function [panorama,weights] = simpleBlend(im1, panorama, weights, dist)
    weight = dist./(max(dist(:)));
    weights = weights + weight;
    if size(im1,3) == 1
        panorama = combine(im1, panorama, weight);
    else
        panorama(:,:,1) = combine(im1(:,:,1), panorama(:,:,1), weight);
        panorama(:,:,2) = combine(im1(:,:,2), panorama(:,:,2), weight);
        panorama(:,:,3) = combine(im1(:,:,3), panorama(:,:,3), weight);
    end
end

function im12 = combine(im1, im2, weight)
    im12 = im2;
    im12 = im12 + weight.*im2double(im1);
end

[hannaLyu]

我终于成功看到结果了。参数用了你的stitching_multiple. 我有两个问题： 1.dist1(1:margin,:) = 1;dist1(end-margin:end,:) = 1;dist1(:,1:margin) = 1;dist1(:,end-margin:end) = 1; 为什么要加margin？ 2.stitch之后的图为什么有明有暗？

[xiahaa]

我终于成功看到结果了。参数用了你的stitching_multiple. 我有两个问题： 1.dist1(1:margin,:) = 1;dist1(end-margin:end,:) = 1;dist1(:,1:margin) = 1;dist1(:,end-margin:end) = 1; 为什么要加margin？ 2.stitch之后的图为什么有明有暗？

1. 不要margin也可以。
2. 贴个图出来看看~应该不会出现有明有暗的情况。