% initial design variables
function [x, vmin, vmax, numVars] = microDensityInit(nelx,nely)
    x = ones(nely, nelx);
    R = sqrt(1.2/pi);

    % init 2
    for i = 0:nelx-1
        for j = 0:nely-1
            if sqrt((i - nelx / 2 - 0.5)^2 + (j - nely / 2 - 0.5)^2) < min(nelx, nely) *(R/2)
                x(j+1, i+1) = 0;
            end
        end
    end
    
    vmin = zeros(nely*nelx,1);
    vmax = ones(nely*nelx,1);
    numVars = nelx*nely;
    
%     a = [-0.4, 0.4];
%     b = [-0.4, -0.4];
%     c = [0.4, 0];
%     [cx,cy] = meshgrid((0:nelx)/nelx - 1/2, (0:nely)/nely - 1/2);
%     for i = 1:nelx
%         for j = 1:nely
%             p = [cx(j,i), cy(j,i)];
%             if isPointInTriangle(a,b,c,p)
%                 x(j,i) = 0;
%             end
%         end
%     end
%     x = int64(x);
%     imshow(x);
end

% is point p in triangle abc
function re = isPointInTriangle(a,b,c,p)
    ab = b-a;
    ac = c-a;
    ap = p-a;
    
    u = (dot(ab,ab)*dot(ap,ac) - dot(ac,ab)*dot(ap,ab)) / ...
        (dot(ac,ac)*dot(ab,ab) - dot(ac,ab)*dot(ab,ac));
    v = (dot(ac,ac)*dot(ap,ab) - dot(ac,ab)*dot(ap,ac)) / ...
        (dot(ac,ac)*dot(ab,ab) - dot(ac,ab)*dot(ab,ac));
    
    if (u<1 || abs(u-1)<1e-5) && (u>0 || abs(u)<1e-5) && ...
       (v<1 || abs(v-1)<1e-5) && (v>0 || abs(v)<1e-5) && ...
       (u+v<1 || abs(u+v-1)<1e-5)
        re = 1;
    else 
        re = 0;
    end
end