Аппроксимация

[Igors]

Добрый день

Дано: на вход поступают точки одна за одной (напр QPointF). Гарантируется что x каждой новой точки больше предыдущей (напр x - время).

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

Спасибо

[Akon]

А в чем конкретно (формула) выражается точность аппроксимации?

[ssoft]

Не все условия задачи указаны.
Как конкретно нужно аппроксимировать.
- линейно
- кусочно-линейно
- сплайнами
- полиномами
- и т.д.
- др.

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

[Disa]

Доброе утро.

А в чем конкретно (формула) выражается точность аппроксимации?

Если для нахождения минимума этой функции использовать метод градиентного спуска, то это не так важно. (итоговая формулуа, что для МНК с полиномом, что для сигмоиды (например) одинаковая)

y = h(x), h - гипотеза (например, h = t0 + t1x + t2x^2 + .... tnx^n)

минимизируем J(t) = 1/m СУММ(1..m)(g(h(x), y)), где g(h) - функция точности (cost function), m - число точек. Берем метод градиентного спуска:

repeat {
t0 = t0 - alpha dJ(t) / d(t0)
t1 = t1 - alpha dJ(t) / d(t1)
...
tn = tn - ...
}

(dJ(t) / d(tn) - частная производная J(t) по tn, alpha - параметр метода)

Дальше в зависимости от вида гипотезы и функции цены (точности) можно делать различные оптимизации и изменять или заменять метод градиентного спуска, чем-то другим (решая СЛАУ, просто перемножаю матрицы по формулам, в лоб пересчитывать производные и прочее)

[ssoft]

Я, наверное, понял, что имелось ввиду.
Скорее всего нужна не аппроксимация, а загрубление.

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

Управление "точностью" несколькими параметрами - угол, отклонение и длина.

1. Берем точку p0.
    Строим векторы к следующей p1 и последующей p2.
2. Если дистанция от точки (p0, p1) больше заданной, помечаем точку p1, как используемую, делаем ее текущей (заменяем p0 на p1, p1 на p2, p2 на p3) и переходим к 1.
3. Если угол между векторами (p0, p2) (p0, p1) больше заданного, помечаем точку p1, как используемую, делаем ее текущей и переходим к пункту 1.
4. Если отклонение точки p1 до линии (p0,p2) больше заданной, помечаем точку p1, как используемую, делаем ее текущей и переходим к 1.
5 ... любые другие критерии, аналогично предыдущим трем пунктам.
6. Игнорируем точку p1, заменяем p1 на p2, p2 на p3. Переходим к 1.

До тех пор пока не кончатся точки.

[Igors]

А в чем конкретно (формула) выражается точность аппроксимации?

Было бы конечно прекрасно если бы кто-то давал формулы, и наше дело сводилось к их реализации. Но увы, этого почему-то не происходит  

Я, наверное, понял, что имелось ввиду.
Скорее всего нужна не аппроксимация, а загрубление.

А "загрубление" это уже не аппроксимация?  Ладно, бог с ними, терминами, по сути совершенно верно. Скорость здесь слишком критична, нет возможности даже накопить все исходные, не говоря уже о градиенте, гипотезе, МНК и.т.п. Чем меньше выходной контейнер - тем лучше т.к. он сохраняется, грубо говоря, для каждого пикселя

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

Управление "точностью" несколькими параметрами - угол, отклонение и длина.

1. Берем точку p0.
    Строим векторы к следующей p1 и последующей p2.
2. Если дистанция от точки (p0, p1) больше заданной, помечаем точку p1, как используемую, делаем ее текущей (заменяем p0 на p1, p1 на p2, p2 на p3) и переходим к 1.
3. Если угол между векторами (p0, p2) (p0, p1) больше заданного, помечаем точку p1, как используемую, делаем ее текущей и переходим к пункту 1.
4. Если отклонение точки p1 до линии (p0,p2) больше заданной, помечаем точку p1, как используемую, делаем ее текущей и переходим к 1.
5 ... любые другие критерии, аналогично предыдущим трем пунктам.
6. Игнорируем точку p1, заменяем p1 на p2, p2 на p3. Переходим к 1.

До тех пор пока не кончатся точки.

Это первое что приходит в голову. но эта схема не ловит плавный переход. Допустим на каждом шаге угол меньше критического и ф-ция монотонно убывает (или возрастает). Тогда вектор (p0, p1) будет "разворачивться" по направлению роста. Напр на шаге 100 мы получим вектор (p0. p100), 98 точек выброшено - гуд. Но угол между (p0, p100) и (p0, p2) может быть ого-го. Залатать это "длиной"  - неясно как.

[ssoft]

Это первое что приходит в голову. но эта схема не ловит плавный переход. Допустим на каждом шаге угол меньше критического и ф-ция монотонно убывает (или возрастает). Тогда вектор (p0, p1) будет "разворачивться" по направлению роста. Напр на шаге 100 мы получим вектор (p0. p100), 98 точек выброшено - гуд. Но угол между (p0, p100) и (p0, p2) может быть ого-го. Залатать это "длиной"  - неясно как.

Да и правда, потомучто я наврал, не специально, конечно.
Сначала смотрим отклонение p1 относительно линии (p0, p2). Если p1 можно проигнорировать, то смотрим до какой точки, т.е. p2 и последующие анализируем относительно линии (p0, p1). Там достаточно все просто, нужно порисовать чуть-чуть. Алгоритм очень быстро работает, причем является однопроходным, я использую его для загрубления графики.

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

Код:

#ifndef QEX_GRAPHICS_TOLLERANCE_H
#define QEX_GRAPHICS_TOLLERANCE_H

#include "config.h"

class WIN_EXPORT QexGraphicsTollerance
{
private:
	double absolute_distance;
	double relative_distance;

public:
	QexGraphicsTollerance ();

	QexGraphicsTollerance & setAbsoluteDistance ( double distance );
	double absoluteDistance () const;

	QexGraphicsTollerance & setRelativeDistance ( double distance );
	double relativeDistance () const;

	bool operator == ( const QexGraphicsTollerance & other ) const;
	bool operator != ( const QexGraphicsTollerance & other ) const;
};

#endif

Код:

#ifndef QEX_GRAPHICS_CURVE_FILTER_H
#define QEX_GRAPHICS_CURVE_FILTER_H

#include "QexGraphicsTollerance.h"

typedef QVector< int > glIndexArray;
typedef QVector< QexPoint3Df > glPointArray;

class WIN_EXPORT QexGraphicsCurveFilter
{
private:
	QexGraphicsTollerance filter_tollerance;

public:
	QexGraphicsCurveFilter ( const QexGraphicsTollerance & tollerance );

	void setTollerance ( const QexGraphicsTollerance & tollerance );
	const QexGraphicsTollerance & tollerance () const;

	glIndexArray filter ( const glPointArray & points ) const;
	glIndexArray filter ( const glPointArray & points, glIndexArray indices ) const;

	glIndexArray filter ( int point_count, const QexPoint3Df * points ) const;
	glIndexArray filter ( const QexPoint3Df * points, int index_count, const int * indices = 0 ) const;

	GLfloat length ( const QexPoint3Df * points, int index_count, const int * indices ) const;
};

#endif // QEX_GRAPHICS_CURVE_FILTER_H

Код:

#include "QexGraphicsCurveFilter.h"

QexGraphicsCurveFilter::QexGraphicsCurveFilter ( const QexGraphicsTollerance & tollerance )
	: filter_tollerance( tollerance )
{
};

void QexGraphicsCurveFilter::setTollerance ( const QexGraphicsTollerance & tollerance )
{
	filter_tollerance = tollerance;
};

const QexGraphicsTollerance & QexGraphicsCurveFilter::tollerance () const
{
	return filter_tollerance;
};

glIndexArray QexGraphicsCurveFilter::filter ( const glPointArray & points ) const
{
	return filter( points.count(), points.constData() );
};

glIndexArray QexGraphicsCurveFilter::filter ( const glPointArray & points, glIndexArray indices ) const
{
	return filter( points.constData(), indices.count(), indices.constData() );
};

glIndexArray QexGraphicsCurveFilter::filter ( int point_count, const QexPoint3Df * points ) const
{
	return filter( points, point_count, 0 );
};

#define CF_INDEX( i ) ( indices ? indices[ ( i ) ] : ( i ) )

GLfloat QexGraphicsCurveFilter::length ( const QexPoint3Df * points, int index_count, const int * indices ) const
{
	GLfloat result( 0 );
	for ( int i = 1; i < index_count; ++i )
	{
		result += QexVector3Df::abs( QexVector3Df( points[ CF_INDEX( i ) ], points[ CF_INDEX( i - 1 ) ] ) );
	}
	return result;
}

glIndexArray QexGraphicsCurveFilter::filter ( const QexPoint3Df * points, int index_count, const int * indices ) const
{
	if ( !index_count )
		return glIndexArray();

	GLfloat ref_abs_dist= qMax< GLfloat >( tollerance().absoluteDistance(), 0 );
	GLfloat ref_rel_dist = qMax< GLfloat >( tollerance().relativeDistance(), 0 );

	if ( ref_rel_dist )
		ref_rel_dist *= length( points, index_count, indices );

	GLfloat ref_distance = ref_abs_dist && ref_rel_dist
		? qMin( ref_abs_dist, ref_rel_dist )
		: qMax( ref_abs_dist, ref_rel_dist );

	glIndexArray ret;
	int index_first = CF_INDEX( 0 );
	ret.append( index_first ); // всегда добавляется первый

	int index;
	for ( int i = 1; i < index_count; ++i )
	{
		index = CF_INDEX( i );

		if ( ref_distance )
		{
			index_first = CF_INDEX( i - 1 );
			QexVector3Df vector( points[ index_first ], points[ index ] );

			// проверяем index_next
			int index_next;
			for ( int j = i + 1; j < index_count; ++j )
			{
				index_next = CF_INDEX( j );

				GLfloat diviation = QexVector3Df::abs( QexVector3Df( points[ index ], points[ index_next ] ) );
				if ( diviation < ref_distance )
					continue;

				QexVector3Df vector_next( points[ index_first ], points[ index_next ] );
				diviation = QexVector3Df::abs( vector ) * sin( QexVector3Df::angle( vector, vector_next ) );
				if ( diviation >= ref_distance )
				{
					i = j - 1;
					ret.append( index );
					break;
				}
			}
		}
		else
		{
			ret.append( index );
		}
	}

	 // добавляется последний, если еще не добавлен
	index = CF_INDEX( index_count - 1 );
	if ( ret.last() != index || ret.count() == 1 )
		ret.append( index );

	return ret;
};

#undef CF_INDEX

Возможно, критерии загрубления немного другие.

[Igors]

Код:

diviation = QexVector3Df::abs( vector ) * sin( QexVector3Df::angle( vector, vector_next ) );

Если QexVector3Df::abs возвращает длину вектора, то лучше так

Код

C++ (Qt)
diviation = fabs(QexVector3Df::cross( vector, vector_next )) / QexVector3Df::abs( vector_next );
 

Возможно, критерии загрубления немного другие.

Да. У Вас все норм, но требуется разумный ref_dist. У меня возможна (и желательна) ситуевина когда аппроксимируется 2-3 отрезками и даже одним, независимо от длины.

В любом случае спасибо за код 

[Igors]

Вот нашел сильный алгоритм (аттач)

[_Vitaliy_]

Сейчас тоже работаю с похожей тематикой. Рекомендую посмотреть в сторону МНК или если можно интерполировать, то сплайн Акимы или кубический сплайн.

[Igors]

Рекомендую посмотреть в сторону МНК или если можно интерполировать, то сплайн Акимы или кубический сплайн.

Все-таки лучше прочитать тему, а потом (может быть) "рекомендовать" 

Скорость здесь слишком критична, нет возможности даже накопить все исходные, не говоря уже о градиенте, гипотезе, МНК и.т.п.

[Disa]

А что за книга из которой рисунок?

[Igors]

А что за книга из которой рисунок?

Это статья 2000 года (Locovic, Veach) - видимо "первоисточник" придумал скромный Локович. Потом десятки перепечаток в которых почему-то график вверх ногами (наверное не заметили). Как всегда возня вокруг GPU и.т.п.

[Disa]

Спасибо!)