- PVSM.RU - https://www.pvsm.ru -
Статья представляет собой нечто вроде “практического руководства” для построения весьма интересных трехмерных графиков функций вида 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=l1x+l2y+l3z
Y=m1x+m2y+m3z
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=x2+y2:
Вид сверху на функцию где видно ту самую сетку и показаны значения в узлах.
Вид на функцию под углом.
Возможно, стоит по-другому изменять углы, при вращении используя все 3 угла и использовать сетку с шагом меньше 1, но мы упростили ситуацию.
Воспользуемся алгоритмом Движущихся квадратов [4] (“Marching squares”), в Википедии дано достаточно подробное описание. Гугл также выдает очень хорошую статью [5], где описывается этот алгоритм и реализовывается на C#.
Суть кратко:
1. Требуется найти начальную позицию — откуда пойдет изолиния.
2. Затем сравнить значения в соседних узлах сетки, которые образуют квадрат с вершинами: (xi,yj), (xi-1,yj), (xi,yj-1), (xi-1,yj-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=x2+y2:
Как видно на картинке алгоритм довольно успешно справился и отделил точки где значение функции выше 5, но немного правее и выше. Изолиния получилась угловатой, поэтому проведем интерполяцию. Смысл интерполяции в том, что мы на основании значениях z в соседних узлах сетки вычисляем более близкое значение x, или y, к реальной изолинии, что делает изолинию более правдоподобной.
Будем использовать формулу линейной интерполяции:
x=(Z-f(xi-1,yj)/(f(xi,yj)-f(xi-1,yj))+xi-1,
y=(Z-f(xi,yj-1)/(f(xi,yj)-f(xi,yj-1))+yj-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
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/programmirovanie/5162
Ссылки в тексте:
[1] Chart: http://msdn.microsoft.com/en-us/library/system.windows.forms.datavisualization.charting.chart.aspx
[2] матрицу поворота: http://ru.wikipedia.org/wiki/Матрица_поворота
[3] углы Эйлера: http://ru.wikipedia.org/wiki/Углы_Эйлера
[4] Движущихся квадратов: http://ru.wikipedia.org/wiki/Marching_squares
[5] очень хорошую статью: http://devblog.phillipspiess.com/2010/02/23/better-know-an-algorithm-1-marching-squares/
[6] тут: http://narod.ru/disk/45760351001.6ac3b489e8c920a3091f7610505ed0b8/masq-demo.zip.html
[7] тут: http://rghost.ru/37520611
[8] ru.wikipedia.org/wiki/Интерполяция: http://ru.wikipedia.org/wiki/Интерполяция
Нажмите здесь для печати.