% Step 1: Read the star image
imageFileName = '4.png'; % Replace with your star image file
starImage = imread(imageFileName);

% Convert the image to grayscale
grayStarImage = rgb2gray(starImage);

% Step 2: Detect edges using Canny edge detector
edges = edge(grayStarImage, 'canny');

% Step 3: Detect line segments using Hough transform
[H, T, R] = hough(edges);
P = houghpeaks(H, 100, 'threshold', ceil(0.5 * max(H(:)))); % Increase threshold
lines = houghlines(edges, T, R, P, 'FillGap', 10, 'MinLength', 15); % Increase FillGap and MinLength

% Visualize the detected line segments (optional)
figure;
imshow(starImage);
hold on;
for k = 1:length(lines)
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:, 1), xy(:, 2), 'LineWidth', 2, 'Color', 'red');
end
title('Detected Line Segments');
hold off;

% Step 4: Find intersection points of line segments more robustly
intersectionPoints = [];
for i = 1:length(lines)
    for j = i+1:length(lines)
        % Line 1
        x1 = lines(i).point1(1);
        y1 = lines(i).point1(2);
        x2 = lines(i).point2(1);
        y2 = lines(i).point2(2);
        
        % Line 2
        x3 = lines(j).point1(1);
        y3 = lines(j).point1(2);
        x4 = lines(j).point2(1);
        y4 = lines(j).point2(2);
        
        % Calculate intersection
        denom = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4);
        if denom ~= 0
            px = ((x1*y2 - y1*x2)*(x3 - x4) - (x1 - x2)*(x3*y4 - y3*x4)) / denom;
            py = ((x1*y2 - y1*x2)*(y3 - y4) - (y1 - y2)*(x3*y4 - y3*x4)) / denom;
            if px >= 1 && px <= size(starImage, 2) && py >= 1 && py <= size(starImage, 1)
                intersectionPoints = [intersectionPoints; px py];
            end
        end
    end
end

% Remove outliers using the interquartile range method
if ~isempty(intersectionPoints)
    Q1 = prctile(intersectionPoints, 25);
    Q3 = prctile(intersectionPoints, 75);
    IQR = Q3 - Q1;
    inlierIdx = all(intersectionPoints >= (Q1 - 1.5 * IQR) & intersectionPoints <= (Q3 + 1.5 * IQR), 2);
    filteredPoints = intersectionPoints(inlierIdx, :);
else
    filteredPoints = intersectionPoints;
end

% Use the median of the filtered points to find the central point
center = median(filteredPoints, 1);

% Convert center to integer values for indexing
center = round(center);

% Step 5: Mark the same center on the original star image
figure;
imshow(starImage);
hold on;
plot(center(1), center(2), 'r+', 'MarkerSize', 10, 'LineWidth', 2);
title('Star Image with Center Marked using Hough Transform');
hold off;

