Сумма площадей прямоугольников с учетом их пересечения

double CalcSquare(List<Rect> rects)
{
    // 1. У нас есть два отсортированных массива с границами (по одному для каждой координаты).
    double[] xlimits = rects.Select(p => p.TopLeft.X).Concat(rects.Select(p => p.BottomRight.X)).Distinct().OrderBy(x => x).ToArray();
    double[] ylimits = rects.Select(p => p.TopLeft.Y).Concat(rects.Select(p => p.BottomRight.Y)).Distinct().OrderBy(x => x).ToArray();
    
    // 2. Для каждого прямоугольника "закрашиваем" клетки матрицы, которые лежат внутри него.
    bool[,] matrix = new bool[xlimits.Length - 1, ylimits.Length - 1];
    foreach(Rect rect in rects)
    {
        int imin = Array.BinarySearch(xlimits, rect.TopLeft.X);
        int imax = Array.BinarySearch(xlimits, rect.BottomRight.X);
        int jmin = Array.BinarySearch(ylimits, rect.TopLeft.Y);
        int jmax = Array.BinarySearch(ylimits, rect.BottomRight.Y);
        
        for(int i = imin; i < imax; i++)
            for(int j = jmin; j < jmax; j++)
                matrix[i, j] = true;
    }
    
    // 3. Суммируем площади закрашенных клеток.
    double square = 0;
    for(int i = 0; i < xlimits.Length - 1; i++)
    {
        double width = xlimits[i+1] - xlimits[i];
        double height = 0;
        for(int j = 0; j < ylimits.Length - 1; j++)
            if(matrix[i,j])
                height += ylimits[j+1] - ylimits[j];
        square += width * height;
    }
    return square;
}
 
 
public struct Rect
{
    public Point TopLeft;
    public Point BottomRight;
}
 
public struct Point
{
    public double X;
    public double Y;
}

Tuple<double[], double[]> GetLimits(List<Rect> rects)
{
    double[] xlimits = new double[rects.Count * 2];
    double[] ylimits = new double[rects.Count * 2];
    
    int i = 0;
    foreach(Rect rect in rects)
    {
        xlimits[i] = rect.TopLeft.X;
        ylimits[i++] = rect.TopLeft.Y;
        xlimits[i] = rect.BottomRight.X;
        ylimits[i++] = rect.BottomRight.Y;
    }
    Array.Sort(xlimits);
    Array.Sort(ylimits);
    
    int i1 = 1, j1 = 1;
    for(int k = 1; k < xlimits.Length; k++)
    {
        if(xlimits[k] != xlimits[k-1])
            xlimits[i1++] = xlimits[k];
        if(ylimits[k] != ylimits[k-1])
            ylimits[j1++] = ylimits[k];
    }
    Array.Resize(ref xlimits, i1);
    Array.Resize(ref ylimits, j1);
 
    return Tuple.Create(xlimits, ylimits);
}

Dictionary<double,int> xdict = xlimits.Select((x,i) => new {x,i}).ToDictionary(p => p.x, p => p.i);