Создание треугольных сеток на сфере

Эта страница является черновиком статьи.

Два-три предложения.

Постановка задачи

Регулярная...

Генерация сетки в сферическом треугольнике

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

Метод бисекции

Назовём бисекцией операцию деления исходного треугольника на четыре треугольника нового поколения. Собственно термин «бисекция» относится к делению сторон пополам. В середины рёбер вставляются новые вершины (белые точки на рисунках), которые соединяются новыми рёбрами (пунктирные линии), образующими новые треугольники. Следующее поколение получается очередной бисекцией.

В терминах геометрии на сфере задача вставки точек в стороны треугольников решается последовательным решением обратной и прямой геодезических задач. Однако в данном случае гораздо проще использовать векторную алгебру. Пусть концы стороны заданы векторами a и b; тогда средняя точка f вычисляется как их нормированная сумма:

$\mathbf {f} ={\frac {\mathbf {a} +\mathbf {b} }{|\mathbf {a} +\mathbf {b} |}}$

Первая бисекция	Вторая бисекция
Трисекция

Метод трисекции

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

Проще всего вычислить положение центральной точки g:

$\mathbf {g} ={\frac {\mathbf {a} +\mathbf {b} +\mathbf {c} }{|\mathbf {a} +\mathbf {b} +\mathbf {c} |}}$

где a, b и c — векторы вершин исходного треугольника.

Разделить стороны на три равных отрезка сложнее. Удобное решение предлагает утилита PROJ.4 geod:

$ geod +a=радиус_Земли +lat_1=широта_1 +lon_1=долгота_1 +lat_2=широта_2 +lon_2=долгота_2 +n_S=3

Здесь параметры lat_1, lon_1, lat_2, lon_2 задают начало и конец линии, а параметр n_S определяет число отрезков. Результатом будут широты и долготы четырёх точек, лежащих на равных расстояниях вдоль линии.

Программная реализация

Приведём пример программы на Си генерации сетки в треугольнике, реализующей метод бисекции:

/* ========================================================================== *
 * triangulate
 *
 * Программа генерирует треугольную сетку в сферическом треугольнике.
 *
 * Usage: triangulate <triangle> <sections>
 *
 * Arguments:
 *    triangle - файл с координатами трёх точек
 *    sections - количество бисекций
 * ========================================================================== */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "triangulate.h"

int main(int argc, char *argv[])
{
  s_tri *tri;
  char buf[1024], basename[1024], fname[1024], *ptr;
  int i, j, k, n, zoom, n_tri;
  int m, m0, m1;
  double x[3], y[2], lon, lat;
  FILE *fp;

  if (argc < 3) {
    puts("Usage: triangulate <triangle> <bisections>");
    exit(EXIT_SUCCESS);
  }
  zoom = atoi(argv[2]);

  strcpy(basename, argv[1]);
  if ((ptr = strrchr(basename, '.')) != NULL)
    *ptr = '\0';

  /* рассчитать число треугольников всех уровней */
  n = 1;     /* число отрезков разбиения базового треугольника */
  n_tri = 1; /* полная сумма треугольников */
  for (k = 0; k < zoom; k++) {
    n <<= 1;        /* бисекция */
    n_tri += n * n; /* добавить треугольники к сумме */
  }

  /* выделить память под треугольники */
  if ((tri = calloc(n_tri, sizeof(s_tri))) == NULL) {
    printf("can't allocate %d triangles\n", n_tri);
    exit(EXIT_FAILURE);
  }

  /* прочесть координаты точек базового треугольника */
  if ((fp = fopen(argv[1], "r")) == NULL) {
    printf("can't open %s\n", argv[1]);
    exit(EXIT_FAILURE);
  }
  for (i = 0; i < 3; i++) {
    if (fgets(buf, 1024, fp) == NULL) {
      printf("error reading %s\n", argv[1]);
      exit(EXIT_FAILURE);
    }
    sscanf(buf, "%lf %lf", &lon, &lat);
    y[0] = Radians(lon);
    y[1] = Radians(lat);
    SpherToCart(y, tri[0].pt[i]);
  }
  fclose(fp);

  /* начать вывод в файл точек */
  strcpy(fname, basename);
  strcat(fname, "_p.MIF");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  fputs("Version   300\n", fp);
  fputs("Charset \"WindowsCyrillic\"\n", fp);
  fputs("Delimiter \",\"\n", fp);
  fputs("CoordSys Earth Projection 1, 0\n", fp);
  fputs("Columns 1\n", fp);
  fputs("  id Decimal(10, 0)\n", fp);
  fputs("Data\n", fp);
  fputs("\n", fp);
  for (i = 0; i < 3; i++) {
    CartToSpher(tri[0].pt[i], y);
    fprintf(fp, "Point %f %f\n", Degrees(y[0]), Degrees(y[1]));
    fputs("    Symbol (49,0,12) \n", fp);
  }

  tri[0].full = 1;

  n = 1;     /* число отрезков разбиения базового треугольника */
  m0 = 0;    /* номер 1-го треугольника предыдущего уровня */
  m1 = 1;    /* номер 1-го треугольника следующего уровня */
  for (k = 0; k < zoom; k++) {
    n <<= 1;        /* бисекция */
    for (m = 0; m < m1 - m0; m++) {

      InsertPoint(tri[m0 + m].pt[1], tri[m0 + m].pt[2], tri[m1 + m*4].pt[0]);
      InsertPoint(tri[m0 + m].pt[2], tri[m0 + m].pt[0], tri[m1 + m*4].pt[1]);
      InsertPoint(tri[m0 + m].pt[0], tri[m0 + m].pt[1], tri[m1 + m*4].pt[2]);
      tri[m1 + m*4].full = !tri[m0 + m].full;

      memcpy(tri[m1 + m*4 + 1].pt[0], tri[m0 + m].pt[0], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 1].pt[1], tri[m1 + m*4].pt[2], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 1].pt[2], tri[m1 + m*4].pt[1], sizeof(double) * 3);
      tri[m1 + m*4 + 1].full = tri[m0 + m].full;

      memcpy(tri[m1 + m*4 + 2].pt[0], tri[m1 + m*4].pt[2], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 2].pt[1], tri[m0 + m].pt[1], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 2].pt[2], tri[m1 + m*4].pt[0], sizeof(double) * 3);
      tri[m1 + m*4 + 2].full = tri[m0 + m].full;

      memcpy(tri[m1 + m*4 + 3].pt[0], tri[m1 + m*4].pt[1], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 3].pt[1], tri[m1 + m*4].pt[0], sizeof(double) * 3);
      memcpy(tri[m1 + m*4 + 3].pt[2], tri[m0 + m].pt[2], sizeof(double) * 3);
      tri[m1 + m*4 + 3].full = tri[m0 + m].full;

      if (!tri[m0 + m].full)
	continue;

      /* вывести точки */
      for (j = 0; j < 3; j++) {
	CartToSpher(tri[m1+m*4].pt[j], y);
	fprintf(fp, "Point %f %f\n", Degrees(y[0]), Degrees(y[1]));
	fputs("    Symbol (49,0,12) \n", fp);
      }
      i += 3;
    }
    m0 = m1;
    m1 += n * n;
  }
  fclose(fp);

  /* создать атрибуты точек */
  strcpy(fname, basename);
  strcat(fname, "_p.MID");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  for (i = 0; i < (n + 1) * (n + 2) / 2; i++)
    fprintf(fp, "%d\n", i + 1);
  fclose(fp);

  /*
   * Вывод треугольников
   */

  /* вывести треугольники в файл полигонов */
  strcpy(fname, basename);
  strcat(fname, "_a.MIF");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  fputs("Version   300\n", fp);
  fputs("Charset \"WindowsCyrillic\"\n", fp);
  fputs("Delimiter \",\"\n", fp);
  fputs("CoordSys Earth Projection 1, 0\n", fp);
  fputs("Columns 1\n", fp);
  fputs("  id Decimal(10, 0)\n", fp);
  fputs("Data\n", fp);
  fputs("\n", fp);
  n = 1;
  m0 = 0;
  m1 = 1;
  for (k = 0; k < zoom; k++) {
    n <<= 1;
    m0 = m1;
    m1 += n * n;
  }
  for (m = 0; m < m1 - m0; m++) {
    fputs("Region  1\n", fp);
    fputs("  4\n", fp);
    for (j = 0; j <= 3; j++) {
      i = j % 3;
      CartToSpher(tri[m0+m].pt[i], y);
      fprintf(fp, "%f %f\n", Degrees(y[0]), Degrees(y[1]));
    }
    fputs("    Pen (1,2,0) \n", fp);
    fputs("    Brush (1,0,16777215)\n", fp);
    for (j = 0; j < 3; j++)
      x[j] = tri[m0 + m].pt[0][j]
	+ tri[m0 + m].pt[1][j] + tri[m0 + m].pt[2][j];
    CartToSpher(x, y);
    fprintf(fp, "    Center %f %f\n", Degrees(y[0]), Degrees(y[1]));
  }
  fclose(fp);

  /* создать атрибуты полигонов */
  strcpy(fname, basename);
  strcat(fname, "_a.MID");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  n = 1;
  m0 = 0;
  m1 = 1;
  for (k = 0; k < zoom; k++) {
    n <<= 1;
    m0 = m1;
    m1 += n * n;
  }
  for (m = 0; m < m1 - m0; m++)
    fprintf(fp, "%d\n", m + 1);
  fclose(fp);

  /*
   * Вывод линий
   */

  /* вывести рёбра в файл линий */
  strcpy(fname, basename);
  strcat(fname, "_l.MIF");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  fputs("Version   300\n", fp);
  fputs("Charset \"WindowsCyrillic\"\n", fp);
  fputs("Delimiter \",\"\n", fp);
  fputs("CoordSys Earth Projection 1, 0\n", fp);
  fputs("Columns 1\n", fp);
  fputs("  id Decimal(10, 0)\n", fp);
  fputs("Data\n", fp);
  fputs("\n", fp);
  n = 1;
  m0 = 0;
  m1 = 1;
  for (k = 0; k < zoom; k++) {
    n <<= 1;
    m0 = m1;
    m1 += n * n;
  }
  for (m = 0; m < m1 - m0; m++) {
    if (!tri[m0 + m].full)
      continue;
    CartToSpher(tri[m0+m].pt[0], y);
    for (j = 0; j < 3; j++) {
      fputs("Pline 2\n", fp);
      fprintf(fp, "%f %f\n", Degrees(y[0]), Degrees(y[1]));
      i = (j + 1) % 3;
      CartToSpher(tri[m0+m].pt[i], y);
      fprintf(fp, "%f %f\n", Degrees(y[0]), Degrees(y[1]));
      fputs("    Pen (1,2,0) \n", fp);
    }
  }

  /* создать атрибуты линий */
  strcpy(fname, basename);
  strcat(fname, "_l.MID");
  if ((fp = fopen(fname, "w")) == NULL) {
    printf("can't create %s\n", fname);
    exit(EXIT_FAILURE);
  }
  j = 0;
  n = 1;
  m0 = 0;
  m1 = 1;
  for (k = 0; k < zoom; k++) {
    n <<= 1;
    m0 = m1;
    m1 += n * n;
  }
  for (m = 0; m < m1 - m0; m++)
    if (tri[m0 + m].full)
      for (i = 0; i < 3; i++)
	fprintf(fp, "%d\n", ++j);
  fclose(fp);

  free(tri);
  return 0;
}

/*
 * Вычисление среднего между двумя положениями
 *
 * Аргументы исходные:
 *     a, b - трёхмерные векторы двух положений
 *
 * Аргументы определяемые:
 *     c - вектор середины
 */
void InsertPoint(double a[], double b[], double c[])
{
  double r;
  int j;

  for (j = 0; j < 3; j++)
    c[j] = a[j] + b[j];
  r = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
  if (r == 0.)
    return;
  for (j = 0; j < 3; j++)
    c[j] /= r;

  return;
}

/*
 * Преобразование сферических координат в вектор
 *
 * Аргументы исходные:
 *     y - {долгота, широта}
 *
 * Аргументы определяемые:
 *     x - вектор {x, y, z}
 */
void SpherToCart(double y[], double x[])
{
  double p;

  p = cos(y[1]);
  x[2] = sin(y[1]);
  x[1] = p * sin(y[0]);
  x[0] = p * cos(y[0]);

  return;
}

/*
 * Преобразование вектора в сферические координаты
 *
 * Аргументы исходные:
 *     x - {x, y, z}
 *
 * Аргументы определяемые:
 *     y - {долгота, широта}
 *
 * Возвращает:
 *     длину вектора
 */
double CartToSpher(double x[], double y[])
{
  double p;

  p = hypot(x[0], x[1]);
  y[0] = atan2(x[1], x[0]);
  y[1] = atan2(x[2], p);

  return hypot(p, x[2]);
}

В архиве Sph_tri.zip этот код находится в файле triangulate.c. Файл triangulate.h содержит необходимые объявления:

#define Radians(x) (x / 57.29577951308232)
#define Degrees(x) (x * 57.29577951308232)

/* структура треугольника */
typedef struct {
  double pt[3][3]; /* координаты вершин */
  int full;        /* 0 - "пустой", 1 - "полный" */
} s_tri;

void InsertPoint(double a[],double b[],double c[]);
double CartToSpher(double x[],double y[]);
void SpherToCart(double y[],double x[]);
int main(int argc,char *argv[]);

Создадим исполняемый модуль triangulate компилятором gcc:

$ gcc -o triangulate triangulate.c -lm

Готовый исполняемый модуль для MS Windows triangulate.exe можно найти в архиве Sph_tri-win32.zip.

Сферические многогранники

Сферический многогранник — разбиение сферы дугами больших окружностей на замкнутые области, называемые сферическими многоугольниками. Способы разбиения сферы ничем не ограничены. Однако регулярные построения обычно основаны на пространственной симметрии тетраэдра, октаэдра или икосаэдра.

Нас интересуют способы разбиения сферы на равносторонние или близкие к равносторонним треугольники. В центре такого треугольника построенная на его основе сетка будет близка к регулярной. Наибольшие искажения формы сетки будут вблизи углов базового треугольника. Несложный анализ показывает, что с позиции сохранения формы выгоднее всего опираться на симметрию икосаэдра. Подходящие многогранники — правильные и полуправильные, образованные треугольниками либо пяти- и шестиугольниками, которые разбиваются на почти равносторонние треугольники. Рассмотрим лишь несколько многогранников, удовлетворяющих этим требованиям.

Икосаэдр	Пентакисдодекаэдр
Икосододекаэдр	Усечённый икосаэдр	Курносый додекаэдр

Икосаэдр — классическое Платоново тело, сложенное двадцатью равносторонними треугольниками.

Пентакисдодекаэдр. Сплошные линии — рёбра классического додекаэдра. Пунктирные линии, разделяя каждую грань додэкаэдра на треугольники, превращают его в пентакисдодэкаэдр.

Икосододекаэдр образован жёлтыми треугольниками и красными пятиугольниками. Разделив пятиугольники на треугольники, как показано пунктирными линиями, получим основу для построения сетки. Полученная фигура не имеет самостоятельного значения, поскольку получается после бисекции граней икосаэдра.

Усечённый икосаэдр, он же «футбольный мяч», образован жёлтыми шестиугольниками и красными пятиугольниками. Разбив все грани на треугольники, получим основу для построения сетки. Полученная фигура структурно совпадает с результатом трисекции граней икосаэдра, хотя геометрия немного отличается: из трёх частей, на которые разделено каждое ребро икосаэдра, средняя несколько длиннее крайних.

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

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

Ссылки

• Задачи на сфере: обратная геодезическая задача
• Задачи на сфере: прямая геодезическая задача
• man_geod – PROJ.4
• Полуправильный многогранник
• Spherical polyhedron