function plotSleepStateScatter(sleepState)

sheepNames = fieldnames(sleepState);
figure
for p = 1:numel(sheepNames)
    
    R1_NomFrom = 6.5;
    R1_NomTo = 9;
    R1_DeNomFrom = 0.5;
    R1_DeNomTo = 9;
    R2_NomFrom = 0.5;
    R2_NomTo = 20;
    R2_DeNomFrom = 0.5;
    R2_DeNomTo = 100;


    % wakeData = sleepState.(sheepName).allFFT(sleepState.(sheepName).allUsed, :);
    % 
    % x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
    % y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
    % subplot(2, 5, p)
    % dscatter(x, y)
    x = [];
    y = [];

    indexesWake = find(sleepState.(sheepNames{p}).isWake);
    indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
    indexesNREM = find(sleepState.(sheepNames{p}).isNREM);
    indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
    wakeData = sleepState.(sheepNames{p}).allFFT(indexesNREM(1:500), :);
    wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexesWake(1:500), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
    x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
    y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
    subplot(2, 5, p)
    % dscatter(x, y)



    gridx1 = linspace(0, 1, 100);
    gridx2 = linspace(0, 1, 100);
    [x1,x2] = meshgrid(gridx1, gridx2);
    x1 = x1(:);
    x2 = x2(:);
    xi = [x1 x2];

    [f,xi] = ksdensity([x, y], xi); 
    dataImage = [];
    [a,b,xI] = unique(xi(:, 1));
    [a,b,yI] = unique(xi(:, 2));
    for i = 1:numel(f)
        dataImage(yI(i), xI(i)) = f(i);
    end
    imagesc(dataImage)


    %%%%%%%%%%%%%%%%%%%


    indexes = find(sleepState.(sheepNames{p}).isWake);
    indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
    x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
    y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
    % dscatter(x, y)



    gridx1 = linspace(0, 1, 100);
    gridx2 = linspace(0, 1, 100);
    [x1,x2] = meshgrid(gridx1, gridx2);
    x1 = x1(:);
    x2 = x2(:);
    xi = [x1 x2];

    [f,xi] = ksdensity([x, y], xi); 
    dataImage = [];
    [a,b,xI] = unique(xi(:, 1));
    [a,b,yI] = unique(xi(:, 2));
    for i = 1:numel(f)
        dataImage(yI(i), xI(i)) = f(i);
    end

    %imagesc(dataImage)
    hold on
    bw = dataImage>10;
    B = bwboundaries(bw);
    for k = 1:length(B)
       boundary = B{k};
       plot(boundary(:,2), boundary(:,1), 'k', 'LineWidth', 2)
    end


    indexes = find(sleepState.(sheepNames{p}).isNREM);
    indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
    x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
    y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
    % dscatter(x, y)


    gridx1 = linspace(0, 1, 100);
    gridx2 = linspace(0, 1, 100);
    [x1,x2] = meshgrid(gridx1, gridx2);
    x1 = x1(:);
    x2 = x2(:);
    xi = [x1 x2];

    [f,xi] = ksdensity([x, y], xi); 
    dataImage = [];
    [a,b,xI] = unique(xi(:, 1));
    [a,b,yI] = unique(xi(:, 2));
    for i = 1:numel(f)
        dataImage(yI(i), xI(i)) = f(i);
    end

    %imagesc(dataImage)
    hold on
    bw = dataImage>10;
    B = bwboundaries(bw);
    for k = 1:length(B)
       boundary = B{k};
       plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)
    end
    caxis([1 30])
    set(gca,'YDir','normal')
    xlim([0 50])
    ylim([30 100])


end

figure
for p = 1:numel(sheepNames)
    subplot(2, 5, p)
    x = [];
    y = [];

    indexes = find(sleepState.(sheepNames{p}).isWake);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:250), :);
    x = [x; sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [y; sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    %dscatter(x, y)


    indexes = find(sleepState.(sheepNames{p}).isWakeStar);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:250), :);
    x = [x; sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [y; sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    %dscatter(x, y)
    %scatter(x, y, 'm', '.')

    indexes = find(sleepState.(sheepNames{p}).isNREM);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
    x = [x; sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [y; sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    %dscatter(x, y)
    %scatter(x, y, 'g', '.')

    indexes = find(sleepState.(sheepNames{p}).isNREMStar);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:100), :);
    x = [x; sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [y; sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    dscatter(x, y)
    %scatter(x, y, 'y', '.')

    xlim([0 .5])
    ylim([.3 1])
end

figure
    for p = 1:numel(sheepNames)
    subplot(2, 5, p)
    x = [];
    y = [];

    indexes = find(sleepState.(sheepNames{p}).isWake);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:250), :);
    x = [x; sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [y; sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    scatter(x, y, 'k', '.')


    indexes = find(sleepState.(sheepNames{p}).isWakeStar);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:250), :);
    x = [sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    scatter(x, y, 'm', '.')

    indexes = find(sleepState.(sheepNames{p}).isNREM);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
    x = [sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    scatter(x, y, 'g', '.')

    indexes = find(sleepState.(sheepNames{p}).isNREMStar);
    wakeData = sleepState.(sheepNames{p}).allFFT(indexes(1:100), :);
    x = [sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2)];
    y = [sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2)];
    hold on
    scatter(x, y, 'y', '.')

    xlim([0 .5])
    ylim([.3 1])
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nSheep = [2, 5, 6, 7, 9];
hdSheep = [1, 3, 4, 8];
figure
%%%%%%%%%%%%%%%%
x = [];
y = [];
wakeData = [];
for s = 1:5
    p = nSheep(s);
    indexesWake = find(sleepState.(sheepNames{p}).isWake);
    indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
    indexesNREM = find(sleepState.(sheepNames{p}).isNREM);
    indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
    wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexesNREM(1:500), :);
    wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexesWake(1:500), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
subplot(1, 2, 1)

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];
[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
totalDensity = sum(sum(dataImage));
dataImage = dataImage/totalDensity;
imagesc(dataImage)


%%%%%%%%%%%%%%%%%%%
wakeData = [];
for s = 1:5
    p = nSheep(s);
    indexes = find(sleepState.(sheepNames{p}).isWake);
    indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
    wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
    wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
dataImage = dataImage/totalDensity;

hold on
bw = dataImage>0.001;
B = bwboundaries(bw);
for k = 1:length(B)
   boundary = B{k};
   plot(boundary(:,2), boundary(:,1), 'k', 'LineWidth', 2)
end
set(gca,'YDir','normal')

wakeData = [];
for s = 1:5
p = nSheep(s);
indexes = find(sleepState.(sheepNames{p}).isNREM);
indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
dataImage = dataImage/totalDensity;


hold on
bw = dataImage>0.001;
B = bwboundaries(bw);
for k = 1:length(B)
   boundary = B{k};
   plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)
end
caxis([0 0.002])
set(gca,'YDir','normal')
xlim([0 50])
ylim([30 100])

%%%%%%%%%%%
%%%%%%%%%%%%%%%%
x = [];
y = [];
wakeData = [];
for s = 1:4
p = hdSheep(s);
indexesWake = find(sleepState.(sheepNames{p}).isWake);
indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
indexesNREM = find(sleepState.(sheepNames{p}).isNREM);
indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
wakeData(end+1:end+500, :)  = sleepState.(sheepNames{p}).allFFT(indexesNREM(1:500), :);
wakeData(end+1:end+500, :) = sleepState.(sheepNames{p}).allFFT(indexesWake(1:500), :);
wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);
subplot(1, 2, 2)

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
totalDensity = sum(sum(dataImage));
dataImage = dataImage/totalDensity;
imagesc(dataImage)


%%%%%%%%%%%%%%%%%%%
wakeData = [];
for s = 1:4
p = hdSheep(s);
indexes = find(sleepState.(sheepNames{p}).isWake);
indexesWakeStar = find(sleepState.(sheepNames{p}).isWakeStar); 
wakeData(end+1:end+500, :)  = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesWakeStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
dataImage = dataImage/totalDensity;


hold on
bw = dataImage>0.001;
B = bwboundaries(bw);
for k = 1:length(B)
   boundary = B{k};
   plot(boundary(:,2), boundary(:,1), 'k', 'LineWidth', 2)
end
set(gca,'YDir','normal')
wakeData = [];
for s = 1:4
p = hdSheep(s);
indexes = find(sleepState.(sheepNames{p}).isNREM);
indexesNREMStar = find(sleepState.(sheepNames{p}).isNREMStar); 
wakeData(end+1:end+500, :)  = sleepState.(sheepNames{p}).allFFT(indexes(1:500), :);
wakeData(end+1:end+100, :) = sleepState.(sheepNames{p}).allFFT(indexesNREMStar(1:100), :);
end

x = sum(wakeData(:, (R1_NomFrom*2 + 1):(R1_NomTo*2 + 1)), 2)./sum(wakeData(:, (R1_DeNomFrom*2 + 1):(R1_DeNomTo*2 + 1)), 2);
y = sum(wakeData(:, (R2_NomFrom*2 + 1):(R2_NomTo*2 + 1)), 2)./sum(wakeData(:, (R2_DeNomFrom*2 + 1):(R2_DeNomTo*2 + 1)), 2);

gridx1 = linspace(0, 1, 100);
gridx2 = linspace(0, 1, 100);
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

[f,xi] = ksdensity([x, y], xi); 
dataImage = [];
[a,b,xI] = unique(xi(:, 1));
[a,b,yI] = unique(xi(:, 2));
for i = 1:numel(f)
    dataImage(yI(i), xI(i)) = f(i);
end
dataImage = dataImage/totalDensity;

hold on
bw = dataImage>0.001;
B = bwboundaries(bw);
for k = 1:length(B)
   boundary = B{k};
   plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)
end
caxis([0 0.002])
set(gca,'YDir','normal')
xlim([0 50])
ylim([30 100])

