%--------------------------
% @Author: Jingqiao Hu
% @Date: 2021-12-16 17:19:07
% @LastEditTime: 2021-12-17 15:15:04

% assume length of microx = microx, i.e. dx = 1
%--------------------------
function [bnodes] = generateCVT(nelx, nely, microx, rho)

%     e1 = 1; nu1 = 0.3;

    globalx = nelx * microx;
    globaly = nely * microx;

%     [~, fixeddofs, F, ~] = designDomain(optDesign, globalx, globaly, vF);

%     rho = top88(globalx, globaly, v0, e1, nu1, F, fixeddofs);
%     save rho rho;
%     load("data/mat/rho.mat");

    % generate quadtree for seeds of CVT
    min_size = microx / 4;
    max_size = microx / 2;
    [blocks, bnodes] = quadtree1(rho, min_size, max_size);
    hold on; scatter(bnodes(:,1), bnodes(:,2),8,'filled','k');

    % transform to [-nelx/2, -nely,2] & [nelx/2, nely/2];
    bnodes(:,1) = bnodes(:,1) - globalx / 2;
    bnodes(:,2) = bnodes(:,2) - globaly / 2;
    write_seeds(bnodes);

    figure; voronoi(bnodes(:,1), bnodes(:,2)); axis equal;
%     bnodes(100,1) = bnodes(100,1) + 0.8;
%     figure(3); voronoi(bnodes(:,1), bnodes(:,2)); 
end

function write_seeds(bnodes)  
    fileID = fopen('data/seeds.cin','w');
    fprintf(fileID,'%f %f\n',bnodes');
    fclose(fileID);
end

function [blocks, bidx] = quadtree1(rho, min_size, max_size)

%     S = qtdecomp(rho, 0.5, [min_size, max_size]); % split-version
    S = qtdecomp(rho, @density_filter, 0.45, [min_size, max_size]); % density-version
    blocks = zeros(size(S));

    show_division(S);

    dim = max_size;
    while dim >= min_size
        numblocks = length(find(S == dim));

        if (numblocks > 0)
            vals = qtgetblk(rho, S, dim);

            valmodes = zeros(size(vals));
            for blknum = 1:size(vals,3)
                valmodes(:,:,blknum) = mode(vals(:,:,blknum),'all');
            end

            blocks = qtsetblk(blocks, S, dim, valmodes);
        end
        dim = dim / 2;        
    end
%     figure(2); imshow(1-blocks); axis equal
    
    bidx = blocks_idx(rho, S, min_size, max_size);
end

function dosplit = density_filter(A, threshold, dims)

    meanVals = mean(mean(A,1),2);
    dosplit = meanVals > threshold;
    
    dosplit = (dosplit & (size(A,1) > dims(1))) | (size(A,2) > dims(2));
end

function bidx = blocks_idx(rho, S, min_size, max_size)

    bidx = [];

    dim = min_size;
    while dim <= max_size
        numblocks = length(find(S == dim));

        if (numblocks > 0)
            [~, r, c] = qtgetblk(rho, S, dim);
            r = r + dim/2; % y
            c = c + dim/2; % x           

            if isempty(bidx)
                bidx = [bidx; c, r]; % [x, y]
%                 hold on; scatter(c, r, 'filled');
            else
                rc = setdiff([c,r], bidx, "rows");
                bidx = [bidx; rc];
%                 hold on; scatter(rc(:,2), rc(:,1), 'filled');
            end
        end
        dim = dim * 2;        
    end
    bidx(:,2) = size(rho,1) - bidx(:,2);
end

function show_division(S)
    blocks = repmat(uint8(0), size(S));

    for dim = [512 256 128 64 32 16 8 4 2 1]
        numblocks = length(find(S == dim));

        if (numblocks > 0)
            values = repmat(uint8(1), [dim dim numblocks]);
            values(2:dim, 2:dim, :) = 0;
            blocks = qtsetblk(blocks, S, dim, values);
        end
    end

    blocks(end, 1:end) = 1;
    blocks(1:end, end) = 1;

    figure(1); clf; imagesc(flipud(1-blocks)); axis equal; axis off;
end