Рисование сеточных графиков трехмерных функций и изолиний к ним

Статья представляет собой нечто вроде “практического руководства” для построения весьма интересных трехмерных графиков функций вида z=f(x,y), с примером реализации на C#.

Рисовать будем в компоненте Chart ^[1], разумеется, можно использовать любой canvas-подобный элемент. При большом количестве точек будут наблюдаться задержки в рисовании, но этот компонент легко использовать и для демонстрации он вполне сгодится.

Строим график

Для этого нам понадобится сетка или поле, каждый узел которой будет иметь координаты x, y и содержать значение z. Эту сетку можно заполнять любым способом, раз мы строим график функции, то будем заполнять ее значениями этой функции (z=f(x,y)). Шаг между двумя последовательными значениями x или y выберем равный одному, для более простой и наглядной демонстрации, размерность сетки N x N. Значения в узлах сетки достаточно считать один раз. Для преобразования трехмерных координат в двухмерные и для поворота графика будем использовать матрицу поворота ^[2] и углы Эйлера ^[3].

Обозначим углы α, β и γ углами a, b и c соответственно и будем использовать их значения в градусах. Запишем матрицу через коэффициенты:

   | l1, l2, l3 |	
 M=| m1, m2, m3 |
   | n1, n2, n3 |
l1 = cos(a) * cos(c) - cos(b) * sin(a) * sin(c)      
m1 = sin(a) * cos(c) + cos(b) * cos(a) * sin(c)       
n1 = sin(b) * sin(c)       
l2 = -cos(a) * sin(c) + cos(b) * sin(a) * cos(c)        
m2 = -sin(a) * sin(c) + cos(b) * cos(a) * cos(c)      
n2 = sin(b) * cos(c)
l3 = sin(b) * sin(a)
m3 = -sin(b) * cos(a)       
n3 = cos(b)

Запишем итоговое преобразование из трехмерных координат в двухмерные:
X=l₁x+l₂y+l₃z
Y=m₁x+m₂y+m₃z
x,y,z –это “внутренние” координаты для графика, X,Y – итоговые экранные координаты. Коэффициенты n1,n2 и n3 нам для этой задачи не понадобятся.

double[,] a; //массив размерностью N x N
…
double X, Y;
//Рисуем вертикальные линии при постоянном значении x
for (int x = 0; x < N; x++)
{
	for (int y = 0; y < N; y++)
	{
		X = l1() * (x - N / 2.0) + l2() * (y - N / 2.0) + l3() * a[x, y];
		Y = m1() * (x - N / 2.0) + m2() * (y - N / 2.0) + m3() * a[x, y];
		chart1.Series[n].Points.AddXY(X, Y);
	}
n++;
}
//Рисуем горизонтальные линии при постоянном значении y
for (int y = 0; y < N; y++)
{
	for (int x = 0; x < N; x++)
	{
		X = l1() * (x - N / 2.0) + l2() * (y - N / 2.0) + l3() * a[x, y];
		Y = m1() * (x - N / 2.0) + m2() * (y - N / 2.0) + m3() * a[x, y];
		chart1.Series[n].Points.AddXY(X, Y);
	}
	n++;
}

Добавим события, при которых будут меняться углы и будет происходить вращение графика. Например, это будет движение курсора при зажатой левой кнопке мышки и изменение углов b и с:

        bool onmove = false;
        Point startpos;
        …
        private void chart1_MouseDown(object sender, MouseEventArgs e)
        {
            if (e.Button == MouseButtons.Left)
            {
                onmove = true;
                startpos = e.Location;
            }
        }
        private void chart1_MouseMove(object sender, MouseEventArgs e)
        {
            if (onmove)
            {
                if ((startpos.Y - e.Y) < 0) b--;
                if ((startpos.Y - e.Y) > 0) b++;
                if ((startpos.X - e.X) < 0) c--;
                if ((startpos.X - e.X) > 0) c++;

                if (b > 359) b = 0;
                if (c > 359) c = 0;
                if (b < 0) b = 359;
                if (c < 0) c = 359;

                drawscene();
            }
        }
        private void chart1_MouseUp(object sender, MouseEventArgs e)
        {
            if (e.Button == MouseButtons.Left) onmove = false;
        }

Проверим на тестовой функции z=x²+y²:

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

Вид на функцию под углом.
Возможно, стоит по-другому изменять углы, при вращении используя все 3 угла и использовать сетку с шагом меньше 1, но мы упростили ситуацию.

Строим изолинии

Воспользуемся алгоритмом Движущихся квадратов ^[4] (“Marching squares”), в Википедии дано достаточно подробное описание. Гугл также выдает очень хорошую статью ^[5], где описывается этот алгоритм и реализовывается на C#.
Суть кратко:
1. Требуется найти начальную позицию — откуда пойдет изолиния.
2. Затем сравнить значения в соседних узлах сетки, которые образуют квадрат с вершинами: (x_i,y_j), (x_i-1,y_j), (x_i,y_j-1), (x_i-1,y_j-1).

3. Выбрать дальнейшее направление обхода, в зависимости от полученного значения на шаге 2, и перейти в следующую точку, в которой опять проделать шаг 2.
4. Выполнять до тех пор, пока не попадем обратно в начальную позицию или не достигнем края сетки.
Всего может быть 16 вариантов:

Не будем писать весь код заново, возьмем часть кода из той статьи ^[5] и изменим под нашу задачу:

        enum dir
        {
            None,
            Up,
            Left,
            Down,
            Right
        }
        dir prevStep;
        dir nextStep;
        bool border;
        int startx, starty;
        void findstartpos()
        {
            for (int y = 0; y < N; y++)
                for (int x = 0; x < N; x++)
                    if (arr[x, y] < Z)
                        {
                            startx = x;
                            starty = y;
                            return;
                        }
        }
        bool check(int x, int y)
        {
            if (x == N - 1 || y == N - 1 || x == 0 || y == 0) border = true;
            if (x < 0 || y < 0 || x >= N || y >= N) return false;
            if (arr[x, y] < Z) return true;
            return false;
        }
        void step(int x, int y)
        {
            bool ul = check(x - 1, y - 1);
            bool ur = check(x, y - 1);
            bool dl = check(x - 1, y);
            bool dr = check(x, y);
            prevStep = nextStep;
            int state = 0;
            if (ul) state |= 1;
            if (ur) state |= 2;
            if (dl) state |= 4;
            if (dr) state |= 8;
            switch (state)
            {
                case 1: nextStep = dir.Down; break;
                case 2: nextStep = dir.Right; break;
                case 3: nextStep = dir.Right; break;
                case 4: nextStep = dir.Left; break;
                case 5: nextStep = dir.Down; break;
                case 6:
                    if (prevStep == dir.Down)
                    {
                        nextStep = dir.Left;
                    }
                    else
                    {
                        nextStep = dir.Right;
                    }
                    break;
                case 7: nextStep = dir.Right; break;
                case 8: nextStep = dir.Up; break;
                case 9:
                    if (prevStep == dir.Right)
                    {
                        nextStep = dir.Down;
                    }
                    else
                    {
                        nextStep = dir.Up;
                    }
                    break;
                case 10: nextStep = dir.Up; break;
                case 11: nextStep = dir.Up; break;
                case 12: nextStep = dir.Left; break;
                case 13: nextStep = dir.Down; break;
                case 14: nextStep = dir.Left; break;
                default:
                    nextStep = dir.None;
                    break;
            }
        }

Попробуем опять же на нашей тестовой функции z=x²+y²:

Как видно на картинке алгоритм довольно успешно справился и отделил точки где значение функции выше 5, но немного правее и выше. Изолиния получилась угловатой, поэтому проведем интерполяцию. Смысл интерполяции в том, что мы на основании значениях z в соседних узлах сетки вычисляем более близкое значение x, или y, к реальной изолинии, что делает изолинию более правдоподобной.
Будем использовать формулу линейной интерполяции:
x=(Z-f(x_i-1,y_j)/(f(x_i,y_j)-f(x_i-1,y_j))+x_i-1,
y=(Z-f(x_i,y_j-1)/(f(x_i,y_j)-f(x_i,y_j-1))+y_j-1,
где Z- значение, на котором нужно провести изолинию.
Интерполирование по координате x, или по координате y, выбирается в зависимости от направления движения для предыдущего и следующего шагов.
Напишем для этого такой, не очень хороший код:

            ...
            List<PointF> res;
            ...
            //Изменение координаты x или y
            int dx = 0, dy = 0;
            switch (prevStep)
           {
               case dir.Down: dy = 1; break;
               case dir.Left: dx = 1;break;
               case dir.Up: dy = -1; break;
               case dir.Right: dx = -1; break;
               default: break;
            }
            ...
            double X = x0 + x;
            double Y = y0 + y;
            if (ip) //ip - interpolation
            {
                //Интерполируем при неизменном направлении
                if (dx != 0 && prevStep == nextStep) 
                       Y = y0 + y + (Z - a[x, y - 1]) / (a[x, y] - a[x, y - 1]) - 1;
                if (dy != 0 && prevStep == nextStep) 
                       X = x0 + x + (Z - a[x - 1, y]) / (a[x, y] - a[x - 1, y]) - 1;
                //Интерполируем при изменении направления
                if (nextStep == dir.Down && prevStep == dir.Left) 
                       Y = y0 + y + (Z - a[x, y - 1]) / (a[x, y] - a[x, y - 1]) - 1;
                if (nextStep == dir.Left && prevStep == dir.Down) 
                       X = x0 + x + (Z - a[x - 1, y]) / (a[x, y] - a[x - 1, y]) - 1;
                if (nextStep == dir.Up && prevStep == dir.Right) 
                       X = x0 + x + (Z - a[x - 1, y]) / (a[x, y] - a[x - 1, y]) - 1;
                if (nextStep == dir.Up && prevStep == dir.Left) 
                       X = x0 + x + (Z - a[x - 1, y]) / (a[x, y] - a[x - 1, y]) - 1;
                if (nextStep == dir.Right && prevStep == dir.Up) 
                       Y = y0 + y + (Z - a[x, y - 1]) / (a[x, y] - a[x, y - 1]) - 1;
                if (nextStep == dir.Right && prevStep == dir.Down) 
                       X = x0 + x + (Z - a[x - 1, y]) / (a[x, y] - a[x - 1, y]) - 1;
                //Исключаем "ненужные" точки
                if (!(nextStep == dir.Down && prevStep == dir.Right) && 
                       !(nextStep == dir.Left && prevStep == dir.Up))
                             res.Add(new PointF((float)X, (float)Y));
            }
            ...

И получаем результат, который нас уже более устраивает:


Пример «ненужных точек», которые мы исключили.
Из возможных улучшений следует отметить, что наша реализация будет строить все изолинии только для монотонно возрастающих (для убывающих следует поменять знак на “>”) функций. Для функций периодически возрастающих и убывающих следует изменить, например функцию поиска начальной позиции и выполнять её несколько раз.

Пример работы итоговой программы.
На практике это можно применять для анализа карт, или для поиска контура на изображении.

Исходный код программы можно скачать тут ^[6] или тут ^[7].
Еще раз ссылки, которые использовались:
ru.wikipedia.org/wiki/Матрица_поворота ^[2]
ru.wikipedia.org/wiki/Углы_Эйлера ^[3]
ru.wikipedia.org/wiki/Интерполяция ^[8]
ru.wikipedia.org/wiki/Marching_squares ^[4]
devblog.phillipspiess.com/2010/02/23/better-know-an-algorithm-1-marching-squares/ ^[5]

Автор: sobak