Drawing Using Fourier Transform

%获取图像边缘
clear
impath = 'C:\Users\phy chen\Desktop\傅里叶画图\best OW.png';%输入图片地址
msizex = 500; %更改图片大小
msizey = 500;
img = imread(impath);%读取图像
grey = rgb2gray(img);
resizeGrey = imresize(grey, [msizex, msizey]);%调整图像大小
eg = edge(resizeGrey, 'canny');%获取图像边缘
[x, y] = find(eg);%查找边缘矩阵中为零的数据并给出位置（行和列）
imshow(eg);


%将图像分解为多个连续图像，并将坐标变换为复数储存在path里面
global Path G;
Path = [];
G = rot90(eg, 3) ;%逆时针旋转数组270°
XSN = struct('path', {});%第一个为设置变量名，第二个设置为空数组
% 如果某个位置为图像边缘（即为白色），
% 则寻找通过dfs函数找到这一连续图像的其他边缘像素，
% 将这些边缘像素坐标转化为复数储存在path，XSN有多少个path证明一个图像有多少个数据
% 然后将原本的这些边缘像素归零，防止对下一图像的数据造成干扰
for i = 1:msizex
    for j = 1:msizey
        if G(i, j) ~= 0
            dfs(i, j);
            % imshow(G);
            % Path
            XSN(end + 1).path = Path;
            Path = [];
        end
    end
end

%傅里叶变换
XnSet = struct('Xn', {}, 'XnSum', {});
for k = 1:length(XSN)
    tmp = fft(XSN(k).path);%将坐标信息xn变换为复数频率向量Xk
    A = exp(-2i * pi .* ([0:length(tmp) - 1]' * [0:length(tmp) - 1]) ./ length(tmp)) ./ length(tmp);%制作一个对称矩阵，其数据为Xk变换为xn的统一前缀Ak
    for j = 1:length(A(:, 1))
        A(j, :) = tmp .* A(j, :);%Ak*Xk=xn 复数频率向量Xk变换为坐标信息xn
    end
    XnSet(end + 1).Xn = A;
    XnSet(end).XnSum = cumsum(A, 2);%对行求累计和，最后一列为坐标向量，其中第j行表示在此之前的频率向量叠加
end

%画图
figure
for i = 1:length(XnSet)
    for j=1:1:(length(XnSet(i).XnSum))/2
        clf       %清空图窗
        hold on   %添加新绘图时保留当前绘图
        plot(XnSet(i).XnSum(j, :));                %只画频率向量叠加

        if i > 1
            for k = 1:i - 1
                plot(XnSet(k).XnSum(:,end),'-r')
            end
        end

        plot(XnSet(i).XnSum(1:j,end),'-r')         %画坐标向量
        hold off

        xlim([0,500]);ylim([0,500])                %画框大小
        pause(0.01)
    end
end

function [] = dfs(x, y)
    global Path G;
    dx = [0, 0, 1, 1, 1, -1, -1, -1];
    dy = [1, -1, 0, -1, 1, 0, 1, -1];
    G(x, y) = 0;
    Path(end + 1) = x + y * 1i;
    for i = 1:8
        if G(x + dx(i), y + dy(i)) ~= 0
            dfs(x + dx(i), y + dy(i));
            Path(end + 1) = x + y * 1i;
        end

    end
end