#%%
import numpy as np
import cv2

def merge_intervals(intervals):
    # Merge overlapping intervals.
    if not intervals:
        return []
    intervals.sort(key=lambda x: x[0])
    merged = [intervals[0]]
    for current in intervals[1:]:
        last = merged[-1]
        if current[0] <= last[1]:
            merged[-1] = (last[0], max(last[1], current[1]))
        else:
            merged.append(current)
    return merged

def get_allowed_intervals(neighbor_values, delta):
    # Compute forbidden intervals for this channel from neighbor values.
    forbidden = []
    for c in neighbor_values:
        low = max(0, c - delta)
        high = min(1, c + delta)
        forbidden.append((low, high))
    # Merge overlapping forbidden intervals.
    merged = merge_intervals(forbidden)
    allowed = []
    # Allowed region before the first forbidden interval.
    if merged:
        if merged[0][0] > 0:
            allowed.append((0, merged[0][0]))
        # Allowed regions between forbidden intervals.
        for i in range(1, len(merged)):
            allowed.append((merged[i-1][1], merged[i][0]))
        # Allowed region after the last forbidden interval.
        if merged[-1][1] < 1:
            allowed.append((merged[-1][1], 1))
    else:
        allowed.append((0, 1))
    return allowed

def sample_from_allowed(allowed):
    # Choose a random value from the allowed intervals with probability proportional to interval length.
    lengths = [end - start for start, end in allowed]
    total_length = sum(lengths)
    if total_length <= 0:
        return None
    point = np.random.rand() * total_length
    cumulative = 0
    for (start, end), length in zip(allowed, lengths):
        cumulative += length
        if point <= cumulative:
            return np.random.uniform(start, end)
    return allowed[-1][0]

def sample_channel(neighbor_values, delta):
    # Compute allowed intervals and sample a value.
    allowed = get_allowed_intervals(neighbor_values, delta)
    if not allowed:
        return None
    return sample_from_allowed(allowed)

def generate_valid_color(neighbors, r):
    """
    Generate a random color that satisfies, for each neighbor,
    |R - R_neighbor| >= sqrt(r), and similarly for G and B.
    neighbors: list of neighbor colors (each is an array of shape (3,))
    r: threshold for squared difference.
    """
    delta = np.sqrt(r)
    color = np.zeros(3)
    for channel in range(3):
        neighbor_channel = [n[channel] for n in neighbors]
        sampled = sample_channel(neighbor_channel, delta)
        # If the allowed intervals are empty, use a fallback.
        if sampled is None:
            sampled = np.random.rand()
        color[channel] = sampled
    return color

def generate_color_pattern(h, w, d, r):
    """
    Generate a dense color pattern.
    h, w: height and width of the output image in pixels.
    d: patch size in pixels (each patch is d x d).
    r: randomness regulation factor. For each neighbor the condition is:
       (|channel_new - channel_neighbor|)^2 >= r.
    """
    pattern = np.zeros((h, w, 3), dtype=np.float32)
    num_patches_x = h // d
    num_patches_y = w // d

    for i in range(num_patches_x):
        for j in range(num_patches_y):
            # Collect neighbor patch colors (using mean because patches have uniform color).
            neighbors = []
            if i > 0:
                patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
                neighbors.append(np.mean(patch_above, axis=(0, 1)))
            if j > 0:
                patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
                neighbors.append(np.mean(patch_left, axis=(0, 1)))
                
            if len(neighbors) == 0:
                # For the first patch use a random color.
                new_color = np.random.rand(3)
            else:
                new_color = generate_valid_color(neighbors, r)
                
            pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
            print('Patch ({}, {})'.format(i, j))
    return pattern

def main():
    # Parameters for the pattern.
    # For example, using 100x100 pixels image with patches of 3x3 pixels.
    height = 2000
    width = 2000
    patch_size = 3
    # Set r to a value that gives sufficient difference.
    # Note: r is the threshold for squared difference per channel, and the condition is
    # |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
    r = 1000

    pattern = generate_color_pattern(height, width, patch_size, r)
    # Scale to 0-255 for display and saving.
    pattern_display = (pattern * 255).astype(np.uint8)

    # Display the pattern with OpenCV.
    cv2.imshow('Pattern', pattern_display)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

    # Save the pattern as an image file.
    cv2.imwrite('color_pattern.png', pattern_display)

if __name__ == '__main__':
    main()