%pollen data from: http://lib.stat.cmu.edu/datasets/
%dimensions of pollen data: ridge, nub, crack, weight, density
%Paper: Uncovering Clusters in Crowded Parallel Coordinates Visualizations
%Bresenham algorithm from wikipedia
%Purpose: highlight the most dense regions of the parallel coordinates
%dataset to indentify structures lying withing
%output: a discrete raster file i.e. image of the density plot, but doesn't
%tell you which points are in the regions of the highest density
%Author: Ankur Kumar (ak364)

function main()
%for the financial dataset
% load 'C:\Documents and Settings\labuser\Desktop\data\aapl.csv';
% aapl = ddt(flipud(aapl(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\amzn.csv';
% amzn = ddt(flipud(amzn(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\ba.csv';
% ba = ddt(flipud(ba(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\bidu.csv';
% bidu = ddt(flipud(bidu(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\bsc.csv';
% bsc = ddt(flipud(bsc(:,5)));
% %dow jones needs to be scaled for better viewing
% load 'C:\Documents and Settings\labuser\Desktop\data\dji.csv';
% dji = ddt(flipud(dji(:,5)/(max(dji(:,5))-min(dji(:,5)))));
% load 'C:\Documents and Settings\labuser\Desktop\data\goog.csv';
% goog = ddt(flipud(goog(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\grmn.csv';
% grmn = ddt(flipud(grmn(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\gs.csv';
% gs = ddt(flipud(gs(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\isrg.csv';
% isrg = ddt(flipud(isrg(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\leh.csv';
% leh = ddt(flipud(leh(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\lmt.csv';
% lmt = ddt(flipud(lmt(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\mer.csv';
% mer = ddt(flipud(mer(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\nok.csv';
% nok = ddt(flipud(nok(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\ptr.csv';
% ptr = ddt(flipud(ptr(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\rimm.csv';
% rimm = ddt(flipud(rimm(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\slb.csv';
% slb = ddt(flipud(slb(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\usdyen.csv';
% usdyen = ddt(flipud(usdyen(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\vlo.csv';
% vlo = ddt(flipud(vlo(:,5)));
% load 'C:\Documents and Settings\labuser\Desktop\data\xom.csv';
% xom = ddt(flipud(xom(:,5)));
% %end of loading closing price data.. about 5 sectors
% pollen = [aapl amzn ba bidu bsc dji goog grmn gs isrg leh lmt mer nok ptr rimm slb usdyen vlo xom];

%for the pollen dataset
load pollen.csv
pollen = pollen(:,[1:5]);


m = length(pollen(:,1)); %number of entries in dataset
n = length(pollen(1,:)); %number of dimensions in dataset
map = jet(m); %color map
figure(1)
hold on
%plot the entire dataset in parallel coordinates
for i=1:m
    plot(pollen(i,:),'color',map(i,:));
end
%the axes for parallel coords
axisx = [1;n];
axisy = [0;0];
label = [{'r'},{'n'},{'c'},{'w'},{'d'}];
for i = 1:n
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-40;40]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:n]);
set(gca, 'xticklabel', label);
hold off

%from VisionX the Matlab plot's pixel resolution is 435x343 (LxW)
L=435;W=343;
G = zeros(L,W); %raster representation matrix, density matrix
maxs = max(pollen); %max_j for all n-dimensions
mins = min(pollen); %min_j for all n-dimensions
A = zeros(m,n); %auxiliary matrix
%fill auxiliary matrix, but need integers only
%normalize all the datavalues from [1,L-1]
for j=1:n
    for i=1:m
        A(i,j) = 1+(((L-1)*(pollen(i,j)-mins(j)))/(maxs(j)-mins(j)));
    end
end
%make everything integers
A = cast(A,'int32');
%compute frequencies and smooth the frequencies using mean filter for
%better density plot
for i=2:n
    F = zeros(L,L); %frequency matrix
    for r=1:m
        b = A(r,i-1);
        c = A(r,i);
        F(b,c) = F(b,c) + 1; %increment frequency
    end
    F = meanfilter(F); %apply smoothing filter (mean filter) for density
    nc = F.nc; %for retreive.m
    nr = F.nr; %for retreive.m
    F = F.F;
    %assign values to the density matrix, G
    for j=1:L
        for k=1:L
            if(F(j,k)~=0)
                %find coordinates of line segment pixels
                ux = (i-1)*W/n;
                uy = j;
                vx = i*W/n;
                vy = k;
                %use Bresenham algorithm to compute all pixel locations
                %joining pixels (ux,uy) to (vx,vy)
                coords = bresenham(ux,uy,vx,vy);
                for index=1:length(coords(:,1))
                    s=coords(index,1);
                    t=coords(index,2);
                    if(G(t,s)<F(j,k))
                        G(t,s) = F(j,k); %update corresponding pixel intensity
                    end
                end
            end
        end
    end
end

%IPC plot density plot, use hot colormap for better viewing
scale = 30;
% I = min(255, (255*G*scale)/max(max(G))); %grayscale pollen data
I = (G*scale)/max(max(G));
figure(2)
colormap('hot');
image(I);%isolate the structure in the data
save img.mat A nc nr F %save values for retreive.m



%square wave smoothing filter, mean filter 3x3
%apply (1/9)*[1 1 1; 1 1 1; 1 1 1] smoothing filter at every value
function ans = meanfilter(X)
ans.nr=[];ans.nc=[];
limy = length(X(:,1));
limx = length(X(1,:));
ans.F = zeros(limx,limy);
%zero pad X for easy calculation
Y = zeros(limx+2,limy+2);
%embed X in Y
Y(2:1+limx,2:1+limy)=X;
for i=2:limx+1
    for j=2:limy+1
        %calculate mean 
        sum = Y(i-1,j-1)+Y(i,j-1)+Y(i+1,j-1) + Y(i-1,j)+Y(i,j)+Y(i+1,j) + Y(i-1,j+1)+Y(i,j+1)+Y(i+1,j+1);
        ans.F(i-1,j-1) = sum/9; 
        if(sum/9 > 8/9) %help in retreive.m
            ans.nr=[ans.nr;i-1];
            ans.nc=[ans.nc;j-1];
        end
    end
end
ans;
    
%Bresenham line algorithm to calculate pixel location values between two points 
%source: Wikipedia
function loc = bresenham(x0,y0,x1,y1)
steep = abs(y1 - y0) > abs(x1 - x0);
if (steep~=0)
    %swap x0 y0
    temp = x0;
    x0 = y0;
    y0 = temp;
    %swap x1 y1
    temp = x1;
    x1 = y1;
    y1 = temp;
end
if (x0>x1)
    %swap x0 x1
    temp = x0;
    x0 = x1;
    x1 = temp;
    %swap y0 y1
    temp = y0;
    y0 = y1;
    y1 = temp;
end
deltax = x1-x0;
deltay = abs(y1-y0);
error = -deltax/2;
ystep = 0;
y = y0;
if(y0<y1)
    ystep = 1;
else
    ystep = -1;
end
loc = []; %the locations
for x=x0:x1
    if(steep~=0)
        loc = [loc;y x];
    else
        loc = [loc;x y];
    end
    error = error + deltay;
    if (error>=0)
        y = y + ystep;
        error = error - deltax;
    end
end
loc = cast(loc,'int32');

%calculate derivative f'(n) = (f(n+1) - f(n-1))/2
%used for financial dataset
function tarr=ddt(arr)
tarr = []; %temporary array
for i=2:length(arr)-1
    tarr=[tarr; (arr(i+1)-arr(i-1))/2];
end