Есть ли итеративный способ расчета радиусов вдоль линии сканирования? - PullRequest
Новая
       8

Есть ли итеративный способ расчета радиусов вдоль линии сканирования?

1 голос
/ 15 сентября 2009

Я обрабатываю ряд точек, которые имеют одинаковое значение Y, но разные значения X. Я прохожу точки, увеличивая X на единицу. Например, у меня может быть Y = 50, а X - целые числа от -30 до 30. Часть моего алгоритма заключается в нахождении расстояния до начала координат в каждой точке и последующей обработке.

После профилирования я обнаружил, что вызов sqrt при расчете расстояния занимает значительное количество моего времени. Есть ли итеративный способ расчета расстояния?

Другими словами:

Я хочу эффективно рассчитать: r[n] = sqrt(x[n]*x[n] + y*y)). Я могу сохранить информацию из предыдущей итерации. Каждая итерация изменяется при увеличении x, поэтому x[n] = x[n-1] + 1. Я не могу использовать функции sqrt или trig, потому что они работают слишком медленно, за исключением начала каждой строки развертки.

Я могу использовать аппроксимации при условии, что они достаточно хороши (ошибка менее 0,1%), а введенные ошибки гладкие (я не могу связать предварительно рассчитанную таблицу аппроксимаций).

Дополнительная информация: х и у всегда целые числа от -150 до 150

Я собираюсь попробовать пару идей завтра и отметить лучший ответ, основанный на том, который является самым быстрым.

Результаты

Я сделал несколько таймингов

  • Формула расстояния: 16 мс / итерация
  • Интерпретирующее решение Пита: 8 мс / итерация
  • решение по предварительным расчетам wrang-wrang: 8 мс / итерация

Я надеялся, что тест примет решение между двумя, потому что мне нравятся оба ответа. Я собираюсь пойти с Питом, потому что он использует меньше памяти.

Ответы [ 6 ]

4 голосов
/ 15 сентября 2009

Просто чтобы почувствовать это, для вашего диапазона y = 50, x = 0 дает r = 50, а y = 50, x = +/- 30 дает r ~ = 58.3. Вы хотите хорошее приближение для +/- 0,1% или +/- 0,05 абсолютного значения. Точность намного ниже, чем у большинства библиотек.

Два приблизительных подхода - вы вычисляете r на основе интерполяции из предыдущего значения или используете несколько членов подходящего ряда.

Интерполяция от предыдущего r

r = (x 2 + y 2 ) 1/2

др / дх = 1/2. 2х. (x 2 + y 2 ) -1 / 2 = x / r 

    double r = 50;

    for ( int x = 0; x <= 30; ++x ) {

        double r_true = Math.sqrt ( 50*50 + x*x );

        System.out.printf ( "x: %d r_true: %f r_approx: %f error: %f%%\n", x, r, r_true, 100 * Math.abs ( r_true - r ) / r );

        r = r + ( x + 0.5 ) / r; 
    }

Дает: 

x: 0 r_true: 50.000000 r_approx: 50.000000 error: 0.000000%
x: 1 r_true: 50.010000 r_approx: 50.009999 error: 0.000002%
....
x: 29 r_true: 57.825065 r_approx: 57.801384 error: 0.040953%
x: 30 r_true: 58.335225 r_approx: 58.309519 error: 0.044065%

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

Усеченная серия

Ряд Тейлора для sqrt (1 + x) для x около нуля равен

sqrt (1 + x) = 1 + 1/2 x - 1/8 x 2 ... + (- 1/2) n + 1 x п

Используя r = y sqrt (1 + (x / y) 2 ), вы ищете термин t = (- 1/2) n + 1 0.36 n с магнитудой меньше 0,001, log (0,002)> n log (0,18) или n> 3,6, поэтому принятие условий для x ^ 4 должно быть в порядке.

2 голосов
/ 15 сентября 2009 
Y=10000
Y2=Y*Y
for x=0..Y2 do
  D[x]=sqrt(Y2+x*x)

norm(x,y)=
  if (y==0) x
  else if (x>y) norm(y,x) 
  else {
     s=Y/y
     D[round(x*s)]/s
  }

Если ваши координаты гладкие, то идею можно расширить с помощью линейной интерполяции. Для большей точности увеличьте Y.

Идея состоит в том, что s * (x, y) находится на линии y = Y, для которой вы предварительно рассчитали расстояния. Получите расстояние, затем разделите его на s.

Полагаю, вам действительно нужно нужно расстояние, а не его квадрат.

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

Под линейной интерполяцией я подразумеваю изменить D[round(x)] на: 

f=floor(x)
a=x-f
D[f]*(1-a)+D[f+1]*a
1 голос
/ 15 сентября 2009

Ну, всегда есть попытка оптимизировать ваш sqrt, самый быстрый, который я видел, это старое автомобильное землетрясение 3 sqrt:

http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

Тем не менее, поскольку sqrt нелинейный, вы не сможете выполнить простую линейную интерполяцию по вашей линии, чтобы получить свой результат. Лучше всего использовать поиск по таблицам, поскольку это обеспечит вам быстрый доступ к данным. И, поскольку вы, кажется, итерируете целые числа, поиск в таблице должен быть чрезвычайно точным.

1 голос
/ 15 сентября 2009

Это на самом деле не отвечает на ваш вопрос, но может помочь ...

Первые вопросы, которые я хотел бы задать:

  • «Мне вообще нужен sqrt?».
  • «Если нет, то как я могу уменьшить количество sqrts?»
  • тогда ваш: «Могу ли я заменить оставшиеся sqrts умным расчетом?»

Итак, я бы начал с:

  • Вам нужен точный радиус, или радиус-квадрат будет приемлемым? Есть быстрое приближение к sqrt, но, вероятно, недостаточно точное для вашей спецификации.
  • Можете ли вы обработать изображение, используя зеркальные квадранты или восьмые? Обрабатывая все пиксели с одинаковым значением радиуса в пакете, вы можете уменьшить количество вычислений в 8 раз.
  • Можете ли вы предварительно рассчитать значения радиуса? Вам нужна только таблица, которая составляет четверть (или, возможно, восьмую) от размера изображения, которое вы обрабатываете, и таблицу нужно будет только предварительно рассчитать один раз, а затем повторно использовать для многих прогонов алгоритма.

Так что умная математика может быть не самым быстрым решением.

0 голосов
/ 15 сентября 2009

Это что-то вроде HAKMEM :

ПУНКТ 149 (Минский): КРУГОВОЙ АЛГОРИТМ Вот элегантный способ рисовать почти круги на точечно-графическом дисплее: 

NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y + epsilon * NEW(!) X

Это делает очень круглый эллипс по центру в начале координат с его размером определяется начальная точка. эпсилон определяет угловой скорость циркулирующей точки и слегка влияет на эксцентриситет. Если эпсилон - это степень 2, тогда мы не даже нужно умножение, не говоря уже о квадратные корни, синусы и косинусы! «круг» будет идеально стабильным потому что очки скоро станут периодично.

Алгоритм окружности был изобретен ошибка, когда я пытался спасти один зарегистрироваться на дисплее взломать! Бен Герли был удивительный взлом, используя только около шести или семи инструкций, и это было великое чудо. Но это было в основном линейно-ориентированный. Это случилось для меня это было бы интересно есть кривые, и я пытался получить кривая отображения взломать с минимальным инструкции.

0 голосов
/ 15 сентября 2009

Что ж, вы можете начать с x = 0 (вам нужно только вычислить n> = 0 и превратить эти результаты в соответствующее n <0). После этого я бы взглянул на использование производной от sqrt (a ^ 2 + b ^ 2) (или соответствующего греха), чтобы воспользоваться константой dx. </p>

Если это не совсем точно, могу я указать, что это очень хорошая работа для SIMD, которая предоставит вам обратную операцию с квадратным корнем как для SSE, так и для VMX (и модель шейдера 2).

...