function chr_sim_marginals(bT1, bT2, xplot, yplot,minRunLength, gapTolerance, extendLength, tol)
% SCATTERWITHMARGINALS
% I wanted to have theoretical barcodes in the marginals of the plot, this
% adds them.... still quite manual but works
% Also colors red for hotspot regions
%
% bT1, bT2 : structures with .rawBarcode and .rawBitmask for chr1 and chr2
% minRunLength, gapTolerance, extendLength, tol : parameters for extend_chr_sim

%% Figure
figure('Color','w');
scatterSize = [0.2 0.2 0.6 0.6]; % [x y width height]
gap = 0; rightWidth = 0.05; topHeight = 0.07;

% Main scatter axes
axScatter = axes('Position', scatterSize);
hold(axScatter,'on');

% Scatter using extend_chr_sim
extend_chr_sim(xplot, yplot, 'minRunLength', minRunLength, 'gapTolerance', gapTolerance, ...
    'extendLength', extendLength, 'tol', tol);
xmin=min( 1:length(bT1.rawBarcode(bT1.rawBitmask)));
xmax=max(1:length(bT1.rawBarcode(bT1.rawBitmask)));
ymin=min( 1:length(bT2.rawBarcode(bT2.rawBitmask)));
ymax=max(1:length(bT2.rawBarcode(bT2.rawBitmask)));

xlim(axScatter, [xmin xmax]);  % define X limits explicitly
ylim(axScatter, [ymin ymax]);  % define Y limits explicitly
axis(axScatter, 'manual');     % turn off autoscaling

xlabel(axScatter,'Chr1')
ylabel(axScatter,'Chr2')

%% Extract all plotted points from extend_chr_sim
allX = [];
allY = [];

% Get line objects (plotted matches/extensions)
lines = findobj(axScatter, 'Type', 'Line');
for k = 1:numel(lines)
    xData = get(lines(k), 'XData');
    yData = get(lines(k), 'YData');
    allX = [allX, xData(:)'];  % flatten and append
    allY = [allY, yData(:)'];
end

% Get scatter objects (individual points)
scatters = findobj(axScatter, 'Type', 'Scatter');
for k = 1:numel(scatters)
    xData = get(scatters(k), 'XData');
    yData = get(scatters(k), 'YData');
    allX = [allX, xData(:)'];
    allY = [allY, yData(:)'];
end

% Now allX and allY contain *everything* extend_chr_sim plotted

% Get unique coordinates
uniqueX = unique(allX);
uniqueY = unique(allY);

% Define tolerance for parallel-line collapsing
tol = 3;
%% For X
highlightX = [];
for i = 1:length(uniqueX)
    xVal = uniqueX(i);

    % indices of all points at this X
    idxPts = find(allX == xVal);
    yVals = allY(idxPts);

    % Sort Y positions
    yVals = sort(yVals);

    % Count only separated Y's
    if ~isempty(yVals)
        keep = yVals(1);
        for j = 2:length(yVals)
            if abs(yVals(j) - keep(end)) > tol
                keep(end+1) = yVals(j); %#ok<AGROW>
            end
        end

        % If more than 1 "separated" point, mark this X for highlight
        if numel(keep) > 1
            highlightX(end+1) = xVal; %#ok<AGROW>
        end
    end
end

%% For Y
highlightY = [];
for i = 1:length(uniqueY)
    yVal = uniqueY(i);

    % indices of all points at this Y
    idxPts = find(allY == yVal);
    xVals = allX(idxPts);

    % Sort X positions
    xVals = sort(xVals);

    % Count only separated X's
    if ~isempty(xVals)
        keep = xVals(1);
        for j = 2:length(xVals)
            if abs(xVals(j) - keep(end)) > tol
                keep(end+1) = xVals(j); %#ok<AGROW>
            end
        end

        % If more than 1 "separated" point, mark this Y for highlight
        if numel(keep) > 1
            highlightY(end+1) = yVal; %#ok<AGROW>
        end
    end
end


%% Right marginal (chr2)
axRight = axes('Position', [scatterSize(1)+scatterSize(3)+gap scatterSize(2) rightWidth scatterSize(4)]);
plot(axRight, bT2.rawBarcode(bT2.rawBitmask), 1:length(bT2.rawBarcode(bT2.rawBitmask)),'k');
ylim(axRight, axScatter.YLim); % use main scatter Y limits
xlim(axRight, [min(get(axRight,'XLim'))-0.1*range(get(axRight,'XLim')) max(get(axRight,'XLim'))+0.1*range(get(axRight,'XLim'))]); % keep current X limits if needed
set(axRight,'YDir','normal','YTick',[],'XTick',[],'XAxisLocation','top')
%xlabel(axRight,'Intensity')
%axis(axRight,'tight')
axRight.YTick = axScatter.YTick;
axRight.YTickLabel = axScatter.YTickLabel;

% Highlight high-coverage Y pixels
% I dont like that these are literally lines... also too red for me
hold(axRight,'on');
for yy = highlightY
    line([min(get(axRight,'XLim')) max(get(axRight,'XLim'))],[yy yy],'Color','r','LineWidth',1.5);
    alpha(.5)
end
plot(axRight, bT2.rawBarcode(bT2.rawBitmask), 1:length(bT2.rawBarcode(bT2.rawBitmask)),'k');
%% Top marginal (chr1)
axTop = axes('Position', [scatterSize(1) scatterSize(2)+scatterSize(4)+gap scatterSize(3) topHeight]);
plot(1:length(bT1.rawBarcode(bT1.rawBitmask)), bT1.rawBarcode(bT1.rawBitmask),'k');
xlim(axTop, axScatter.XLim); % use main scatter X limits
ylim(axTop, [min(get(axTop,'YLim'))-0.1*range(get(axTop,'YLim')) max(get(axTop,'YLim'))+0.1*range(get(axTop,'YLim'))]); % keep current Y limits
set(axTop,'XDir','normal','XTick',[],'YTick',[],'YAxisLocation','right')
%ylabel(axTop,'Intensity')
%axis(axTop,'tight')
axTop.XTick = axScatter.XTick;
axTop.XTickLabel = axScatter.XTickLabel;

% Highlight high-coverage X pixels
hold(axTop,'on');
for xx = highlightX
    line([xx xx],[min(get(axTop,'YLim')) max(get(axTop,'YLim'))],'Color','r','LineWidth',1.5);
    alpha(.5)
end
plot(1:length(bT1.rawBarcode(bT1.rawBitmask)), bT1.rawBarcode(bT1.rawBitmask),'k');
end
