Простой класс для построения линий уровня двумерной сеточной функции

При обработке данных в предметных областях, связанных с научной деятельностью, часто возникает необходимость в построении и визуализации функции двух независимых переменных. Типичным примером является необходимость визуального представления результатов решения двумерных дифференциальных уравнений в частных производных, получаемых в виде так называемых сеточных функций.

Предлагается простой класс для построения линий уровня (изолиний) функции: Z=F(X,Y) в виде линий на плоскости X-Y, удовлетворяющих уравнениям Z=const (где const — набор заданных значений).

Предполагается, что функция Z задана в виде массива z[J,K] на произвольной сетке с четырехугольными ячейками. Сетка задается двумя массивами x[J,K], y[J,K], где J и K размеры сетки.

Значения функции определены в углах четырехугольной ячейки. В каждой ячейке проверяется прохождение рассчитываемой линии уровня через ее грани и, при условии, что линия проходит через ячейку, вычисляются координаты пересечения линии уровня с гранями. Внутри ячейки линия проводится прямолинейным отрезком.

Исходный текст снабжен подробными комментариями.

Файл LinesLevels.cs:

using System.Collections.Generic;
using System.Linq;
using System.Windows;
namespace WpfLinesLevels
{
    public class LinesOfLevels
    {
        private int J, K;
        private double[,] X;
        private double[,] Y;
        private double[,] Z;
        // Список изолиний
        public List<LineLevel> Lines { get; set; }
        /// <summary>
        /// Подготовка 
        /// </summary>
        /// <param name="_levels">Массив уровней</param>
        /// <param name="_x">Координаты X области</param>
        /// <param name="_y">Координаты Y области</param>
        /// <param name="_z">Сеточная функция</param>
        public LinesOfLevels(double[] _levels, double[,] _x, double[,] _y, double[,] _z)
        {
            Lines = new List<LineLevel>(_levels.Count());
            foreach (double l in _levels)
            {
                Lines.Add(new LineLevel(l));
            }
            X = _x;
            Y = _y;
            Z = _z;
            J = X.GetLength(0);
            K = X.GetLength(1);
        }
        /// <summary>
        /// Расчет изолиний.
        /// </summary>
        public void Calculate()
        {
            for (int j = 0; j < J - 1; j++)
                for (int k = 0; k < K - 1; k++)
                {
                    Ceil ir = new Ceil(j, k, X, Y, Z);
                    for (int l = 0; l < Lines.Count(); l++)
                        ir.AddIntoLineLevel(Lines[l]);
                }
        }
    }
    /// <summary>
    /// Одна изолиния
    /// </summary>
    public class LineLevel
    {
        // Список точек изолинии в виде пар точек
        // принадлежащих одной четырехугольной ячеейке
        public List<PairOfPoints> Pairs { get; set; }
        // Уровень изолинии
        public double Level { get; set; }
        public LineLevel(double _level)
        {
            Level = _level;
            Pairs = new List<PairOfPoints>();
        }
    }
    /// <summary>
    /// Пара точек изолинии, принадлежащая одной ячейке
    /// </summary>
    public class PairOfPoints
    {
        public List<Point> Points { get; set; }
        public PairOfPoints() { Points = new List<Point>(); }
    }
    /// <summary>
    /// Угол ячейки.
    /// Индексы для определения одного угла четырехугольной ячейки
    /// </summary>
    internal struct Dot
    {
        internal int j { get; set; }
        internal int k { get; set; }
        internal Dot(int _j, int _k) { j = _j; k = _k; }
    }
    /// <summary>
    /// Четырехугольная ячейка сетки. Определяет текущую ячейку.
    /// Рассчитывает отрезки изолиний в ячейке
    /// </summary>
    internal class Ceil
    {
        // Углы ячейки
        private Dot[] d = new Dot[4];
        // Координатные точки углов
        private Point[] r = new Point[4];
        // Массивы координат всей области
        private double[,] X;
        private double[,] Y;
        // Массив сеточной функции
        private double[,] Z;
        /// <summary>
        /// Определение ячейки 
        /// Определяется левым нижним углом. Циклы перебора индексов должны быть на 1 меньше размерностей J,K массивов
        /// </summary>
        /// <param name="_j">j - индекс левого нижнего угла</param>
        /// <param name="_k">k - индекс левого нижнего угла</param>
        /// <param name="_x">Массив X[J,K]</param>
        /// <param name="_y">Массив Y[J,K]</param>
        /// <param name="_z">Массив сеточной функции Z[J,K]</param>
        internal Ceil(int _j, int _k, double[,] _x, double[,] _y, double[,] _z)
        {
            d[0] = new Dot(_j, _k);
            d[1] = new Dot(_j + 1, _k);
            d[2] = new Dot(_j + 1, _k + 1);
            d[3] = new Dot(_j, _k + 1);
            X = _x;
            Y = _y;
            Z = _z;
            r[0] = dotPoint(d[0]);
            r[1] = dotPoint(d[1]);
            r[2] = dotPoint(d[2]);
            r[3] = dotPoint(d[3]);
        }
        /// <summary>
        /// Определение координатной точки Point угла
        /// </summary>
        /// <param name="_d">Угол, заданный стуктурой Dot</param>
        /// <returns></returns>
        private Point dotPoint(Dot _d) { return new Point(X[_d.j, _d.k], Y[_d.j, _d.k]); }
        /// <summary>
        /// Определение функции в заданном углу
        /// </summary>
        /// <param name="_d">Угол, заданный стуктурой Dot</param>
        /// <returns></returns>
        private double dotZ(Dot _d) { return Z[_d.j, _d.k]; }
        /// <summary>
        /// Определение пары точек, через которые проходит линия уровня
        /// Точки на границах ячейки определяются линейной интераоляцией.
        /// </summary>
        /// <param name="_l">Значение уровня функции</param>
        /// <returns></returns>
        private PairOfPoints ByLevel(double _l)
        {
            PairOfPoints p = new PairOfPoints();
            // Ребро 0
            if ((dotZ(d[0]) >= _l && dotZ(d[1]) < _l) || (dotZ(d[1]) > _l && dotZ(d[0]) <= _l))
            {
                double t = (_l - dotZ(d[1])) / (dotZ(d[0]) - dotZ(d[1]));
                double x = r[0].X * t + r[1].X * (1 - t);
                double y = r[0].Y * t + r[1].Y * (1 - t);
                p.Points.Add(new Point(x, y));
            }
            // Ребро 1
            if ((dotZ(d[1]) >= _l && dotZ(d[2]) < _l) || (dotZ(d[2]) > _l && dotZ(d[1]) <= _l))
            {
                double t = (_l - dotZ(d[2])) / (dotZ(d[1]) - dotZ(d[2]));
                double x = r[1].X * t + r[2].X * (1 - t);
                double y = r[1].Y * t + r[2].Y * (1 - t);
                p.Points.Add(new Point(x, y));
                if (p.Points.Count == 2) return p;
            }
            // Ребро 2
            if ((dotZ(d[2]) >= _l && dotZ(d[3]) < _l) || (dotZ(d[3]) > _l && dotZ(d[2]) <= _l))
            {
                double t = (_l - dotZ(d[3])) / (dotZ(d[2]) - dotZ(d[3]));
                double x = r[2].X * t + r[3].X * (1 - t);
                double y = r[2].Y * t + r[3].Y * (1 - t);
                p.Points.Add(new Point(x, y));
                if (p.Points.Count == 2) return p;
            }
            // Ребро 3
            if ((dotZ(d[3]) >= _l && dotZ(d[0]) < _l) || (dotZ(d[0]) > _l && dotZ(d[3]) <= _l))
            {
                double t = (_l - dotZ(d[0])) / (dotZ(d[3]) - dotZ(d[0]));
                double x = r[3].X * t + r[0].X * (1 - t);
                double y = r[3].Y * t + r[0].Y * (1 - t);
                p.Points.Add(new Point(x, y));
            }
            return p;
        }
        /// <summary>
        /// Добавление пары точек в линию уравня
        /// </summary>
        /// <param name="_lL">Линия уровня</param>
        internal void AddIntoLineLevel(LineLevel _lL)
        {
            PairOfPoints lp = ByLevel(_lL.Level);
            if (lp.Points.Count > 0) _lL.Pairs.Add(lp);
        }
    }
}

Для демонстрации работы класса предлагается небольшое тестовое приложение WPF, которое строит линии уровня для функции вида: z = x^2 + y^2 на сетке 10 на 10.

Файл MainWindow.xaml:

<Window x:Class="WpfLinesLevels.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
        xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
        xmlns:local="clr-namespace:WpfLinesLevels"
        mc:Ignorable="d"
        Title="Изолинии Z = X^2+Y^2" Height="590" Width="525">
    <Grid Name="grid1">        
    </Grid>
</Window>

И файл кода MainWindow.xaml.cs:

using System.Linq;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Media;
using System.Windows.Shapes;

namespace WpfLinesLevels
{
    /// <summary>
    /// Логика взаимодействия для MainWindow.xaml
    /// </summary>
    public partial class MainWindow : Window
    {
        private double Xmax;
        private double Xmin;
        private double Ymax;
        private double Ymin;
        private double xSt;
        private double ySt;

        public MainWindow()
        {
            InitializeComponent();       
            // Определение уровней, которые будут отображаться
            double[] levels = { 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
            double[,] X = new double[10, 10];
            double[,] Y = new double[10, 10];
            double[,] Z = new double[10, 10];

            // Переменные для пересчета физических координат в экранные
            Xmax = 10;
            Xmin = 0;
            Ymax = 10;
            Ymin = 0;
            xSt = 525 / (Xmax - Xmin);
            ySt = 525 / (Ymax - Ymin);

            // Определение массивов координат и функции
            for (int k = 0; k < 10; k++)
                for (int j = 0; j < 10; j++)
                {
                    X[j, k] = j;
                    Y[j, k] = k;
                    Z[j, k] = j * j + k * k;
                }
            // Создание изолиний
            LinesOfLevels lol = new LinesOfLevels(levels, X, Y, Z);
            // Их расчет
            lol.Calculate();
            // Построение
            DrowLevelLine(lol, X, Y);
        }
        /// <summary>
        /// Метод построения изолиний
        /// </summary>
        /// <param name="lL">Расчитанный объект с изолиниями</param>
        /// <param name="x">массив X координат</param>
        /// <param name="y">массив Y координат</param>
        private void DrowLevelLine(LinesOfLevels lL, double[,] x, double[,] y)
        {
            Canvas can = new Canvas();
            foreach (LineLevel l in lL.Lines)
            {
                foreach (PairOfPoints pp in l.Pairs)
                {

                    if (pp.Points.Count() == 2)
                    {
                        Line pl = new Line();
                        pl.Stroke = new SolidColorBrush(Colors.BlueViolet);
                        pl.X1 = xCalc(pp.Points[0].X);
                        pl.X2 = xCalc(pp.Points[1].X);
                        pl.Y1 = yCalc(pp.Points[0].Y);
                        pl.Y2 = yCalc(pp.Points[1].Y);
                        can.Children.Add(pl);
                    }
                }
            }
            can.Margin = new Thickness(10, 10, 10, 10);
            can.VerticalAlignment = VerticalAlignment.Stretch;
            can.HorizontalAlignment = HorizontalAlignment.Stretch;
            grid1.Children.Add(can);
        }
        /// <summary>
        /// Пересчет физической координаты X в экранную
        /// </summary>
        /// <param name="_x">Физическая кордината X</param>
        /// <returns>Экранная координата X</returns>
        private double xCalc(double _x)
        {
            return xSt * (_x - Xmin);
        }
        /// <summary>
        /// Пересчет физической координаты Y в экранную
        /// </summary>
        /// <param name="_x">Физическая кордината Y</param>
        /// <returns>Экранная координата Y</returns>
        private double yCalc(double _y)
        {
            return ySt * (Ymax - _y);
        }
    }
}

Результат работы тестового примера представлен на рисунке:

Спасибо за внимание!