if ~exist('soccer.mp4', 'file')
    error('Video file soccer.mp4 not found.');
end
video = VideoReader('soccer.mp4');
frameCount = floor(video.Duration * video.FrameRate); 
firstFrame = readFrame(video);
hFig = figure;
imshow(firstFrame);
title('Drag to select ROI (double-click to confirm)');
h = imrect;
roi = wait(h); 
roi = round(roi); 
close(hFig);
firstFrame = imcrop(firstFrame, roi);
hFig = figure; 
imshow(firstFrame);
title('Click on the yellow player''s color');
[xY, yY] = ginput(1);
yellowColor = squeeze(firstFrame(round(yY), round(xY), :));
title('Click on the green player''s color');
[xG, yG] = ginput(1);
greenColor = squeeze(firstFrame(round(yG), round(xG), :));
title('Click on the ball''s center and a point on its radius');
[xc, yc] = ginput(2);
ballRadius = sqrt((xc(1)-xc(2))^2 + (yc(1)-yc(2))^2);
cropSize = round(2 * ballRadius);
x1 = round(xc(1) - cropSize/2);
y1 = round(yc(1) - cropSize/2);
ballTemplate = imcrop(firstFrame, [x1, y1, cropSize, cropSize]);
title('Drag to select axis (goalpost) area (double-click to confirm)');
h = imrect;
axisRect = wait(h);
axisRect = round(axisRect);
axisTemplate = imcrop(firstFrame, axisRect);
axisCenterX = axisRect(1) + axisRect(3)/2; 
axisCenterY0 = axisRect(2) + axisRect(4)/2; 
close(hFig);
processedFrameCount = ceil(frameCount / 20);
greenCenters = zeros(processedFrameCount, 10, 2);
yellowCenter = zeros(processedFrameCount, 2);
ballCenter = zeros(processedFrameCount, 2);
axisCenterY = zeros(processedFrameCount, 1);
greenCentersCorrected = zeros(processedFrameCount, 10, 2);
yellowCenterCorrected = zeros(processedFrameCount, 2);
ballCenterCorrected = zeros(processedFrameCount, 2);
hWait = waitbar(0, 'Tracking in progress...');
outputVideo = VideoWriter('tracking_output_cropped_axis.avi');
open(outputVideo);
searchWidthAxis = axisRect(3) + 20;
frameIdx = 1;
for f = 1:20:frameCount
    video.CurrentTime = (f - 1) / video.FrameRate;
    frame = readFrame(video);
    frame = imcrop(frame, roi);
    labFrame = rgb2lab(frame);
    greenColorLAB = rgb2lab(reshape(greenColor, [1,1,3]));
    distGreen = sqrt(sum((labFrame - greenColorLAB).^2, 3));
    maskGreen = distGreen < 15;
    maskGreen = bwareaopen(maskGreen, 50);
    statsG = regionprops(maskGreen, 'Centroid', 'Area');
    if length(statsG) < 10
        hFig = figure;
        imshow(frame);
        title(sprintf('Frame %d: Green players missing. Click 10 players.', f));
        [xg, yg] = ginput(10);
        greenCenters(frameIdx, :, 1) = xg;
        greenCenters(frameIdx, :, 2) = yg;
        close(hFig);
    else
        [~, idx] = maxk([statsG.Area], 10);
        for i = 1:10
            greenCenters(frameIdx, i, :) = statsG(idx(i)).Centroid;
        end
    end
    yellowColorLAB = rgb2lab(reshape(yellowColor, [1,1,3]));
    distYellow = sqrt(sum((labFrame - yellowColorLAB).^2, 3));
    maskYellow = distYellow < 15;
    maskYellow = bwareaopen(maskYellow, 50);
    statsY = regionprops(maskYellow, 'Centroid', 'Area');
    if isempty(statsY)
        hFig = figure;
        imshow(frame);
        title(sprintf('Frame %d: Yellow player missing. Click directly.', f));
        [xy, yy] = ginput(1);
        yellowCenter(frameIdx, :) = [xy, yy];
        close(hFig);
    else
        [~, maxIdx] = max([statsY.Area]);
        yellowCenter(frameIdx, :) = statsY(maxIdx).Centroid;
    end
    grayFrame = rgb2gray(frame);
    edges = edge(grayFrame, 'Canny');
    [centers, radii] = imfindcircles(edges, [round(ballRadius - 5), round(ballRadius + 5)], 'Sensitivity', 0.9);
    if ~isempty(centers)
        ballCenter(frameIdx, :) = centers(1, :);
    else
        [y, x] = find(edges);
        points = [x, y];
        numCircles = 1;
        try
            idx = kmeans(points, numCircles);
            clusterPoints = points(idx == 1, :);
            if size(clusterPoints, 1) >= 3
                p1 = clusterPoints(1, :); p2 = clusterPoints(2, :); p3 = clusterPoints(3, :);
                midAB = [(p1(1) + p2(1))/2, (p1(2) + p2(2))/2];
                if p2(1) - p1(1) == 0
                    slopeAB = inf;
                else
                    slopeAB = (p2(2) - p1(2)) / (p2(1) - p1(1));
                end
                if slopeAB == 0
                    slopePerpAB = inf;
                else
                    slopePerpAB = -1 / slopeAB;
                end
                midBC = [(p2(1) + p3(1))/2, (p2(2) + p3(2))/2];
                if p3(1) - p2(1) == 0
                    slopeBC = inf;
                else
                    slopeBC = (p3(2) - p3(2)) / (p3(1) - p2(1));
                end
                if slopeBC == 0
                    slopePerpBC = inf;
                else
                    slopePerpBC = -1 / slopeBC;
                end
                if isinf(slopePerpAB)
                    x0 = midAB(1);
                    y0 = slopePerpBC * (x0 - midBC(1)) + midBC(2);
                elseif isinf(slopePerpBC)
                    x0 = midBC(1);
                    y0 = slopePerpAB * (x0 - midAB(1)) + midAB(2);
                else
                    A = [-slopePerpAB, 1; -slopePerpBC, 1];
                    B = [midAB(2) - slopePerpAB * midAB(1); midBC(2) - slopePerpBC * midBC(1)];
                    sol = A \ B;
                    x0 = sol(1);
                    y0 = sol(2);
                end
                ballCenter(frameIdx, :) = [x0, y0];
            else
                ballCenter(frameIdx, :) = prevBallCenter;
            end
        catch
            ballCenter(frameIdx, :) = prevBallCenter;
        end
    end
    prevBallCenter = ballCenter(frameIdx, :);
    grayAxisTemplate = rgb2gray(axisTemplate);
    xStartAxis = max(1, round(axisCenterX - searchWidthAxis/2));
    xEndAxis = min(size(frame,2), round(axisCenterX + searchWidthAxis/2));
    grayFrameStrip = grayFrame(:, xStartAxis:xEndAxis);
    corrMapAxis = normxcorr2(grayAxisTemplate, grayFrameStrip);
    [maxCorrAxis, maxIdxAxis] = max(corrMapAxis(:));
    [ypeakAxis, ~] = ind2sub(size(corrMapAxis), maxIdxAxis);
    yCenterAxis = ypeakAxis - (size(grayAxisTemplate,1) - 1)/2;
    if frameIdx == 1 || maxCorrAxis >= 0.3
        axisCenterY(frameIdx) = yCenterAxis;
        prevAxisCenterY = yCenterAxis;
    else
        axisCenterY(frameIdx) = prevAxisCenterY;
    end
    yCorrection = axisCenterY(frameIdx) - axisCenterY0;
    for i = 1:10
        greenCentersCorrected(frameIdx, i, 1) = greenCenters(frameIdx, i, 1);
        greenCentersCorrected(frameIdx, i, 2) = greenCenters(frameIdx, i, 2) - yCorrection;
    end
    yellowCenterCorrected(frameIdx, 1) = yellowCenter(frameIdx, 1);
    yellowCenterCorrected(frameIdx, 2) = yellowCenter(frameIdx, 2) - yCorrection;
    ballCenterCorrected(frameIdx, 1) = ballCenter(frameIdx, 1);
    ballCenterCorrected(frameIdx, 2) = ballCenter(frameIdx, 2) - yCorrection;
    waitbar(frameIdx / processedFrameCount, hWait, sprintf('%d / %d frames processed', frameIdx, processedFrameCount));
    markedFrame = frame;
    greenPos = squeeze(greenCentersCorrected(frameIdx, :, :));
    markedFrame = insertMarker(markedFrame, greenPos, 'o', 'Color', 'green', 'Size', 5);
    markedFrame = insertMarker(markedFrame, yellowCenterCorrected(frameIdx, :), 'o', 'Color', 'yellow', 'Size', 5);
    markedFrame = insertMarker(markedFrame, ballCenterCorrected(frameIdx, :), 'o', 'Color', 'red', 'Size', 5);
    markedFrame = insertShape(markedFrame, 'Rectangle', axisRect + [0, axisCenterY(frameIdx) - axisCenterY0, 0, 0], 'Color', 'blue', 'LineWidth', 2);
    writeVideo(outputVideo, markedFrame);
    frameIdx = frameIdx + 1;
end
close(hWait);
close(outputVideo);
data = [];
for idx = 1:processedFrameCount
    row = [];
    for i = 1:10
        row = [row, squeeze(greenCentersCorrected(idx, i, :))'];
    end
    row = [row, yellowCenterCorrected(idx, :), ballCenterCorrected(idx, :)];
    data = [data; row];
end
headers = {};
for i = 1:10
    headers{end+1} = sprintf('Green%d_X', i);
    headers{end+1} = sprintf('Green%d_Y', i);
end
headers{end+1} = 'Yellow_X';
headers{end+1} = 'Yellow_Y';
headers{end+1} = 'Ball_X';
headers{end+1} = 'Ball_Y';
writecell([headers; num2cell(data)], 'tracking_results_cropped_axis.xlsx');
prompt = sprintf('Enter frame number to visualize (1 to %d): ', processedFrameCount);
frameToShowIdx = input(prompt);
while ~isnumeric(frameToShowIdx) || frameToShowIdx < 1 || frameToShowIdx > processedFrameCount || mod(frameToShowIdx, 1) ~= 0
    disp('Invalid input. Enter an integer between 1 and processedFrameCount.');
    frameToShowIdx = input(prompt);
end
frameToShow = 1 + (frameToShowIdx - 1) * 20;
video.CurrentTime = (frameToShow - 1) / video.FrameRate;
frame = readFrame(video);
frame = imcrop(frame, roi);
markedFrame = frame;
greenPos = squeeze(greenCentersCorrected(frameToShowIdx, :, :));
markedFrame = insertMarker(markedFrame, greenPos, 'o', 'Color', 'green', 'Size', 5);
markedFrame = insertMarker(markedFrame, yellowCenterCorrected(frameToShowIdx, :), 'o', 'Color', 'yellow', 'Size', 5);
markedFrame = insertMarker(markedFrame, ballCenterCorrected(frameToShowIdx, :), 'o', 'Color', 'red', 'Size', 5);
markedFrame = insertShape(markedFrame, 'Rectangle', axisRect + [0, axisCenterY(frameToShowIdx) - axisCenterY0, 0, 0], 'Color', 'blue', 'LineWidth', 2);
figure;
imshow(markedFrame);
title(sprintf('Tracking result for frame %d (ROI cropped, axis corrected)', frameToShow));
disp('Completed! Saved to tracking_results_cropped_axis.xlsx');
    

data = readmatrix('tracking_results_cropped_axis.xlsx', 'Range', 'A2:X1716');

framecount = size(data, 1);
if framecount == 0
    error('No data found in the file.');
end

averagedistances = zeros(framecount, 1);
stddistances = zeros(framecount, 1);
ballcoords = zeros(framecount, 2);

for frameidx = 1:framecount
    greenpoints = reshape(data(frameidx, 1:20), [2, 10])'; 
    ballpos = data(frameidx, 23:24);

    if any(isnan(ballpos)) || any(isnan(greenpoints(:)))
        averagedistances(frameidx) = NaN;
        stddistances(frameidx) = NaN;
        ballcoords(frameidx, :) = [NaN NaN];
        continue;
    end

    ballcoords(frameidx, :) = ballpos;

    distances = [];
    for i = 1:9
        for j = i+1:10
            dist = norm(greenpoints(i,:) - greenpoints(j,:));
            distances(end+1) = dist;
        end
    end

    averagedistances(frameidx) = mean(distances);
    stddistances(frameidx) = std(distances);
end

validframes = ~isnan(averagedistances) & ~isnan(stddistances) & ~isnan(ballcoords(:,1)) & ~isnan(ballcoords(:,2));
if sum(validframes) < 3
    error('Not enough valid data points.');
end

ballcoords = ballcoords(validframes, :);
averagedistances = averagedistances(validframes);
stddistances = stddistances(validframes);
validcount = sum(validframes);

normalizedstd = (stddistances - min(stddistances)) / (max(stddistances) - min(stddistances));
colormapdata = parula(256);
coloridx = round(normalizedstd * 255) + 1;

gridcountx = 50;
gridcounty = 50;
xlimits = linspace(min(ballcoords(:,1)), max(ballcoords(:,1)), gridcountx + 1);
ylimits = linspace(min(ballcoords(:,2)), max(ballcoords(:,2)), gridcounty + 1);

figure;
subplot(1,2,1);
scatter3(ballcoords(:,1), ballcoords(:,2), averagedistances, 20, colormapdata(coloridx,:), 'filled');
xlabel('Ball X (pixels)');
ylabel('Ball Y (pixels)');
zlabel('Average Distance (pixels)');
title('3D Scatter Plot');
grid on;
view(45, 30);
hold on;
daspect([1 1 0.3]);

gridvalues = nan(gridcounty, gridcountx);
for i = 1:gridcountx
    for j = 1:gridcounty
        inx = ballcoords(:,1) >= xlimits(i) & ballcoords(:,1) < xlimits(i+1);
        iny = ballcoords(:,2) >= ylimits(j) & ballcoords(:,2) < ylimits(j+1);
        incell = inx & iny;
        if any(incell)
            gridvalues(j,i) = mean(normalizedstd(incell));
        end
    end
end

filledgrid = gridvalues;
for iter = 1:2
    tempgrid = filledgrid;
    for i = 1:gridcountx
        for j = 1:gridcounty
            if isnan(filledgrid(j,i))
                neighbors = [];
                for dx = -1:1
                    for dy = -1:1
                        if dx == 0 && dy == 0
                            continue;
                        end
                        xidx = i + dx;
                        yidx = j + dy;
                        if xidx >= 1 && xidx <= gridcountx && yidx >= 1 && yidx <= gridcounty
                            if ~isnan(filledgrid(yidx, xidx))
                                neighbors(end+1) = filledgrid(yidx, xidx);
                            end
                        end
                    end
                end
                if ~isempty(neighbors)
                    tempgrid(j,i) = mean(neighbors);
                end
            end
        end
    end
    filledgrid = tempgrid;
end

subplot(1,2,2);
hold on;
graycolor = [0.8 0.8 0.8];
for i = 1:gridcountx
    for j = 1:gridcounty
        xpatch = [xlimits(i), xlimits(i+1), xlimits(i+1), xlimits(i)];
        ypatch = [ylimits(j), ylimits(j), ylimits(j+1), ylimits(j+1)];
        if isnan(filledgrid(j,i))
            patch(xpatch, ypatch, graycolor, 'EdgeColor', 'none');
        else
            idx = round(filledgrid(j,i)*255) + 1;
            patch(xpatch, ypatch, colormapdata(idx, :), 'EdgeColor', 'none');
        end
    end
end

scatter(ballcoords(:,1), ballcoords(:,2), 12, 'k', 'filled');

xlabel('Ball X (pixels)');
ylabel('Ball Y (pixels)');
title('Smoothed Normalized Std Dev Map');
axis equal;
grid on;
colormap(parula);
clim([0 1]);
colorbar.Label.String = 'Normalized Std Dev';

xlim([min(xlimits), max(xlimits)]);
ylim([min(ylimits), max(ylimits)]);



%
%선수-선수(위), 선수-공(아래)
%



try
    data = readmatrix('tracking_results_cropped_axis.xlsx', 'Range', 'A2:X1716');
catch e
    error('Error reading Excel file: %s', e.message);
end

numFrames = size(data, 1);
if numFrames == 0
    error('No data read from Excel file. Check file content or range.');
end
avgDistances = zeros(numFrames, 1);
stdDistances = zeros(numFrames, 1);
ballPositions = zeros(numFrames, 2); 

for idx = 1:numFrames
    greenPos = reshape(data(idx, 1:20), [2, 10])';  
    ballPos = data(idx, 23:24);  
    if any(isnan(ballPos)) || any(isnan(greenPos(:)))
        warning('Frame %d: NaN values detected. Skipping distance calculation.', idx);
        continue;
    end
    ballPositions(idx, :) = ballPos;
    
    distances = sqrt(sum((greenPos - ballPos).^2, 2)); 
    avgDistances(idx) = mean(distances);
    stdDistances(idx) = std(distances);
end


validIdx = ~isnan(avgDistances) & ~isnan(stdDistances) & ~isnan(ballPositions(:,1)) & ~isnan(ballPositions(:,2));
if sum(validIdx) < 3
    error('Too few valid data points (%d) for fitting. Need at least 3 non-NaN points.', sum(validIdx));
end
ballPositions = ballPositions(validIdx, :);
avgDistances = avgDistances(validIdx);
stdDistances = stdDistances(validIdx);
numValidFrames = sum(validIdx);

normStd = (stdDistances - min(stdDistances)) / (max(stdDistances) - min(stdDistances));
colors = parula(256);
colorIndices = round(normStd * 255) + 1;

figure;
subplot(1,2,1);
scatter3(ballPositions(:,1), ballPositions(:,2), avgDistances, ...
         20, colors(colorIndices,:), 'filled');
hold on;
colormap(parula);
c = colorbar;
c.Label.String = 'Standard Deviation of Distances (pixels)';
xlabel('Ball X Position (pixels)');
ylabel('Ball Y Position (pixels)');
zlabel('Average Distance to Ball (pixels)');
title('3D Scatter of Player Distances');
grid on;
view(45, 30);
axis equal;

subplot(1,2,2);

gridCountX = 50;
gridCountY = 50;
xEdges = linspace(min(ballPositions(:,1)), max(ballPositions(:,1)), gridCountX+1);
yEdges = linspace(min(ballPositions(:,2)), max(ballPositions(:,2)), gridCountY+1);

gridValues = nan(gridCountY, gridCountX);

for i = 1:gridCountX
    for j = 1:gridCountY
        inX = ballPositions(:,1) >= xEdges(i) & ballPositions(:,1) < xEdges(i+1);
        inY = ballPositions(:,2) >= yEdges(j) & ballPositions(:,2) < yEdges(j+1);
        inCell = inX & inY;
        if any(inCell)
            gridValues(j,i) = mean(normStd(inCell));
        end
    end
end

for iter = 1:2
    tempGrid = gridValues;
    for i = 1:gridCountX
        for j = 1:gridCountY
            if isnan(gridValues(j,i))
                neighbors = [];
                for dx = -1:1
                    for dy = -1:1
                        if dx == 0 && dy == 0, continue; end
                        xIdx = i + dx;
                        yIdx = j + dy;
                        if xIdx >= 1 && xIdx <= gridCountX && yIdx >= 1 && yIdx <= gridCountY
                            if ~isnan(gridValues(yIdx, xIdx))
                                neighbors(end+1) = gridValues(yIdx, xIdx);
                            end
                        end
                    end
                end
                if ~isempty(neighbors)
                    tempGrid(j,i) = mean(neighbors);
                end
            end
        end
    end
    gridValues = tempGrid;
end

hold on;
grayColor = [0.8 0.8 0.8];
for i = 1:gridCountX
    for j = 1:gridCountY
        xPatch = [xEdges(i), xEdges(i+1), xEdges(i+1), xEdges(i)];
        yPatch = [yEdges(j), yEdges(j), yEdges(j+1), yEdges(j+1)];
        if isnan(gridValues(j,i))
            patch(xPatch, yPatch, grayColor, 'EdgeColor', 'none');
        else
            colorIdx = round(gridValues(j,i) * 255) + 1;
            patch(xPatch, yPatch, colors(colorIdx,:), 'EdgeColor', 'none');
        end
    end
end

scatter(ballPositions(:,1), ballPositions(:,2), 10, 'k', 'filled');
xlabel('Ball X Position (pixels)');
ylabel('Ball Y Position (pixels)');
title('Smoothed Normalized Std Dev Map');
axis equal;
grid on;
colormap(parula);
clim([0 1]);
colorbar.Label.String = 'Normalized Std Dev';

xlim([min(xEdges), max(xEdges)]);
ylim([min(yEdges), max(yEdges)]);