nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

QuatCore: улучшаем программу аффинного алгоритма

Чтобы было, на чём отлаживать новый компилятор, сразу хочу переписать программу аффинного алгоритма (наша самая крупная программа для QuatCore на данный момент), с учётом новых возможностей, которые представит модуль "непосредственных значений" QuatCoreImm.

К примеру, первые 3 строчки кода:
	SP		StackAdr
	NOP		0		;избежать Write Hazard: только здесь в SP запишется новое значение!
	SP		[SP]


Можно заменить одной:
	SP		Stack


Ещё и выкинуть из сегмента данных строку
StackAdr	dw	Stack	;адрес стека (он недоступен напрямую!)


Мелочь, а приятно. Таких мест не так уж много, но кое-чего есть...


Нахождение наиболее отдалённой точки

		;состояние регистров неизвестно (за искл. PC и SP)
		FindMDD3	proc
		;перво-наперво, находим МДД3
		;как точка, сумма квадратов расстояний до которой максимальна (т.е самая отдалённая)
		;для каждой точки считаем сумму квадр. расстояний до всех остальных, в т.ч до самой себя (это 0, зато код упрощается)
		;внешний цикл (по строкам) - перем. j
		;внутр. цикл - по индексу перем (X или Y) - перем k
		;самый-самый внутр. цикл (по столбцам) - перем. i
				Y		Points2D
				[SP+1]	0	;максимальная отдалённость, инициализируем нулём
				;рег. X будет хранить индекс точки с макс. отд. её иниц. не надо-заведомо на одной из итераций запишется
					
				j		3
		@@j_loop:	k		1	;также от 0 до 3, чтобы все расстояния просуммировать
				[SP]		0	;здесь будет храниться текущий максимум
		@@k_loop:	i		3	;от 0 до 1, т.е значения X и Y
		;а теперь непосредственно прибавляем к акк. ([X+2i+k]-[X+2j+k])^2
		@@i_loop:	Acc		[Y+2j+k]	;загрузили одно значение
				SUB		[Y+2i+k]	;вычитаем второе
				SQRD2		Acc			;возводим в квадрат
				ADD		[SP]		;прибавляем к пред. значению
				[SP]		Acc
				iLOOP		@@i_loop		;теперь то же самое для второй координаты
				kLOOP		@@k_loop
		;хорошо, нашли сумму расстояний от точки j до всех остальных
		;она лежит в Acc и в [SP]
		;нужно сравнить с самым большим, и если больше - заменить
		;самое большое значение мы храним в [SP+1],
		;а его индекс - в [SP+2]
				SUB		[SP+1]	;можно и "пожертвовать" значением в Acc,
		;т.к в [SP] лежит точно такое же!
				JL		@@skip
		;дистанция оказалась больше - надо обновить максимум и его индекс
				[SP+1]	[SP]	;двухпортовая память-это зашибись!
				X		j
		@@skip:	jLOOP		@@j_loop
		;по окончании этих циклов, у нас в X лежит индекс точки, которая является МДД3
		;осталось эту точку поменять местами с точкой под индексом 0...
		;здесь i=j=k=0 (все циклы завершились!)
				i		X
				X		Points2D
				Z		Points2D
				CALL		SwapPoints	;потёрли текущий максимум (лежал в [SP])-и хрен с ним
		
		;на этом первая часть марлезонского балета закончена.
		FindMDD3	endp


Как ни странно, здесь мы ничего улучшить не можем, и так всё относительно неплохо. Если здесь и есть откровенно "дурацкий" код - это последние строки, где мы перекладываем значения из регистра в регистр, перед вызовом процедуры SwapPoints. Эта процедура чуть более универсальна, чтобы и при сортировке против часовой стрелки её можно было вызывать, поэтому хочешь-не хочешь, а единую канву пришлось выдержать. Что-то я плясал с бубном вокруг этих трёх регистров, X,Y,Z, ничего шибко лучше не получилось.

Сортировка против часовой стрелки
	;состояние регистров к данному моменту:
	;X = Y =Z = Points2D
	;i = индекс наиболее отдалённой точки  (не понадобится)
	;j=k=0
	;Inv неизвестен
	;C = [Fx0]  (не понадобится)
	;Acc = разность полусумм квадратов (не понадобится)
	SortCCW		proc
	;Далее: нужно расставить оставшиеся точки против часовой стрелки
	;Первый шаг к этому: переместить начало координат в МДД3, который сейчас по нулевому индексу
			X		Fx1	;чтобы индекс от 2 до 0 соотв. точкам (Fx1,Fy1) ... (Fx3, Fy3)
			Z		Fx2	;чтобы иметь сдвинутую на 1 адресацию
	;потом вернём как ни в чём не бывало, у нас фикс. точка, все правила арифм. соблюдаются						
			CALL		ShiftOrigin
	;вот теперь пора сортировать!
			j		1
			i		1
	;Делаем "сдвинутую" адресацию [Fx2 + j], т.е j=1 указывает на 3-ю точку, j=0 - на 2-ю.
	;но при этом, i = 1 указывает на 2-ю, i=0 - на 1-ю
	;а нулевую пока вообще не рассматриваем - она на своём месте!
	;хотим на позицию j+2 поставить "самую большую", для этого надо выбрать между j+1,...1
	;здесь заведомо k=0, так что X+2i / X+2j не нужны!
	@@loop:		C		[Z+2j+k]
			MUL		[X+2i+1]
			C		[Z+2j+1]
			FMS		[X+2i+k]	;нахождение "векторного произведения"
			JL		@@skip
		;меняем местами, как обычно, точки с коорд. i и с коорд. j
			CALL 		SwapPoints
	@@skip:		iLOOP		@@loop
			jLOOP		@@loop
					
			CALL		ShiftOrigin		
			;а теперь ещё надо два последних поменять местами, для совместимости с алг. сопровождения
			;здесь i=j=k=0
			X		Fx3
			CALL		SwapPoints
	SortCCW		endp


И здесь тоже ничего не попишешь. Для работы с двумерными массивами у нас всё-таки очень неплохо всё предусмотрено.

Да и процедуры SwapPoints и ShiftOrigin особенно не поменяешь:
;меняет местами точку с базой Z, инд. j  и базой X, инд. i
;изменяет регистр k (по окончании выходит 0), и C
SwapPoints	proc
		k		1
		NOP		0	;хотим избежать warning'ов
@@swap:		C		[Z+2j+k]
		[Z+2j+k]	[X+2i+k]
		[X+2i+k]	C
		kLOOP		@@swap
		JMP		[--SP]
SwapPoints	endp	

	;отражает точки относительно наиболее отдалённой. Повторное применение этой операции возврщает точки на место!
	;состояние регистров к данному моменту:
	;X = Points2D+2 = Fx1 (т.е указывает на точку с индексом 1)
	;Y = Points2D (точка с индексом 0)
	;Z неважен (вообще Fx2, точка с индексом 2)
	;i неважен (вообще это индекс наиболее отдалённой точки, а потом 0)
	;j,k неважны (вообще j=k=0)
	;Inv неизвестен
	;C неважен (вообще [Fx0] при первом вызове, одна из коорд. при втором) 
	;Acc неважен (вобще разность полусумм квадратов при первом вызове, результат сравн. при втором)
ShiftOrigin proc
	;надо из всех точек вычесть коорд. наиболее отдалённой (с индексом 0)
	;кроме неё самой, разумеется
		k		1
@@k_loop:	i		2	;как всегда, чтобы Y и X
@@i_loop:	Acc		[Y+k]
		SUB		[X+2i+k]
		[X+2i+k]	Acc
		iLOOP		@@i_loop
		kLOOP		@@k_loop
		JMP		[--SP]
ShiftOrigin endp


Нахождение матрицы аффинного преобразования
	;состояние регистров к этому моменту
	;Y = Points2D
	;X = Fx3
	;Z = Fx2
	;i=j=k=0
	;Inv неизвестен
	;C, Acc - неважно
	Compute4PointAffine	proc
				X		AffineMat
				Z		Matrix ;вообще AfTransf, но это промежуточная матрица, её занесём сюда (пока не наступило сопровождение, эти ячейки напрочь не нужны!)
	;по сути, Z = X * Y, где все операнды - матрицы.
				k		1	;номер строки результата (и строки AffineMat)
		@@k_loop:	j		2	;номер столбца результата (и столбца Points2D)
		@@j_loop:	i		3	;номер столбца AffineMat и строки Points2D
				ZAcc		RoundZero	;обнулить до 1/2 мл. разр
		@@i_loop:	C		[X+4j+i]
				FMA		[Y+2i+k]
				iLOOP		@@i_loop
	;очередной результат готов...
				[Z+2j+k]	Acc
				jLOOP		@@j_loop
				kLOOP		@@k_loop
	;вот, какбе и всё...	
	Compute4PointAffine	endp


И тут ничего улучшить не представляется возможным, и так по-моему максимально информативна каждая строка.

Нахождение крена

	;состояние регистров к этому моменту
	;X = AffineMat (константная матрица 3х4, чтобы посчитать преобр)
	;Y = Points2D
	;Z = Matrix (матрица 2х2 и вектор 1х2, пока не AfTransf)
	;i=j=k=0
	;Inv неизвестен
	;C, Acc - пофиг наверное
	FindRoll	proc
		;нам понадобятся [Z]+[Z+3] и [Z+1]-[Z+2]
			Y		QuatY 	;в этом кватернионе построим,используя Y/Z значения как временные
						;(потом они занулятся)
			;а можно Y вместо X, если надо
			i		1
			j		3
		@@sico:	Acc		[Z+i]
			Inv		i
			PM		[Z+i^j]
			[Y+i]		Acc
			iLOOP		@@sico
	
	;теперь наших друзей надо отнормировать, причём динамический диапазон очень велик
	;максимум 13 итераций, давайте столько и делать...
	;здесь у нас появляется беззнаковое число, без него получалось бы очень грустно	
			CALL NormSiCo
	;здесь у нас i=j=k=Inv=0
	
	;теперь синус-косинус превращаем в кватернион, то есть в синус-косинус половинного угла
	;если co>0, то
	;QuatX = (1+abs(co))/2,
	;QuatY = si/2
	;в противном случае
	;QuatX = si/2,
	;QuatY = (1+abs(co))/2
	
	;сейчас X = AffineMat (матрица 3х4, уже неактуальна)
	;Y = QuatY (компоненты co / si),
	;Z = Matrix (матрица 2х2 и вектор 1х2, к ним ещё вернёмся)
	;i=j=k=Inv=0
				
;	ABSP1D2	[Z+2j+k]	;самая укуренная команда-взять модуль B, прибавить единицу, и поделить на 2!
;замена для ABSP1D2:
			ABS		[Y+k]
			X		OneHalf				
			DIV2		Acc
			ADD		[X+k]
;конец замены
	;флаг S (sign) сейчас показывает знак co, что для нас очень полезно
	;но от адресов [Y+S] и [Y+~S] мы пока что избавились. Придётся их "сэмулировать"
			JGE		@@skip
			i		1		;выходит i=S
		@@skip:	X		QuatA
			[X+i]		Acc
			j		1
			DIV2		[Y+1]
			Y		QuatA
			[X+i^j]		Acc			;по сути, Y+(~S)
			CALL		NormSiCo	;подготовили первые 2 компонента кватерниона
					
	;сейчас X=Y=QuatA (2 компоненты кватерниона),
	;Z = Matrix
	;теперь помножаем две матрицы, одна в явном виде - Matrix (Z)
	;вторая - задана неявно, в виде co / si (в QuatY)
	;результат идет в AfTransf
	;здесь i=1, k неизвестен, j=0, Inv=0
			Y		QuatY
			X		AfTransf	;сюда результат складируем
			i		1	;номер строки результата, и строки co/si
	@@i_loop:	k		1	;номер столбца результата, и столбца AfTransf
	@@k_loop:	j		1	;номер столбца co/si и строки AfTransf
			ZAcc		RoundZero	;обнулить до 1/2 мл. разр
	@@j_loop:	C		[Y+i^j]	
			FMPM		[Z+2j+k]
			jLOOP		@@j_loop
			[X+2i+k]	Acc
			kLOOP		@@k_loop
			iLOOP		@@i_loop

	;здесь у нас i=j=k=0
	;теперь у нас должна получиться матрица, включающая в себя масштаб и ракурс (ужатие вдоль некоторой оси), но не крен...
	;ещё кватернион подчистим - по Y,Z может лежать всякий мусор
			[Y+k]		0
			[Y+1]		0
	FindRoll	endp


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

OneHalf	dw		16384 ;для нахождения кватерниона крена.


и затем обращались к нему через регистр X:

	X		OneHalf				
	DIV2		Acc
	ADD		[X+k]


Деление на два запихали посередине, чтобы избежать Hazard'а, когда в регистр X ещё не успело записаться новое значение, а мы уже обращаемся через него.

Выкидываем эту константу из сегмента данных, и прибавляем напрямую:
	ADD		16384


О да, целых два слова сэкономили, одно в коде, другое в данных...

Нахождение масштаба

	;Состояние регистров к этому моменту:
	;X=AfTransf,
	;Y=QuatY,
	;Z=Matrix,
	;i=j=k=Inv=0,
	;C=Acc=Txx
	FindScale	proc
	;1. Вычисляем величины |X-Z| (т.е |Txx - Tyy| ) и |2Y| (например, |Txy+Tyx|, хотя можно 2|Txy| или 2 |Tyx| )
	;важно, что в данный момент Inv=0
			[SP]		0	;манхэттенская метрика, сумма модулей (метрика городских кварталов)
			[SP+1]	0	;чебышевская, максимум от модулей 
			;текущее значение будем хранить в регистре C - дешево и сердито!
			i		1
			j		3					
	@@loop:		Acc		[X+i]	;+Tyx (i=1), -Tyy (i=0)
			PM		[X+i^j]	;Txy (i=1), Txx (i=0)   ;поменяли местами чтобы проверить модуль
			ABS		Acc		;|Txy+Tyx| (i=1), |Txx-Tyy| (i=0)
			C		Acc		;сохранили это значение
			ADD		[SP]	
			[SP]		Acc	;обновили манхэттенскую
			Acc		C
			SUB		[SP+1]
			JL		@@skipCheb ;лучше подошел бы JGE и поменять местами C и [SP+1]. Разница при нуле-на 1 больше операций.
			[SP+1]	C	;обновили Чебышевскую	
	@@skipCheb:	Inv		1		;так что на следующей итерации будет "-"
			iLOOP		@@loop
	;теперь в [SP] лежит манхэттенская метрика, а в [SP+1] - чебышевская
	;осталось посчитать выражение целиком - сумму Txx+Tyy с весом 1/2, а эти - с весами 0,332 и 0,168
	;или сначала всё сложить - а потом поделить пополам...
	;можно, если ввести команду UDIV2 - беззнаковое деление
	;по сути, SHR1 
	
	;кстати, у нас до сих пор j=3 :)
	;i=k=0
			Acc		[X+i]	;добавили Txx
			ADD		[X+i^j]	;и Tyy
	;теперь наши метрики берём с весами
			Z		MetricW		
			C		[SP]
			FMA		[Z+i]
			C		[SP+1]
			FMA		[Z+1]
				
	;ага, сделали
	;продемонстрируем результат
			JO		@@endLoop
	@@shl1:		ADD		Acc
			j++		0
			JNO		@@shl1		
	;ага, сделали
	@@endLoop:	[SP+1]		Acc			;финт ушами: заносим значение 32767, ибо нефиг! (другими способами его сложно раздобыть)
			DIV2		UAC ;получается отрицательное значение от -16384 до -1 (число -1 при "делении на 2" так и останется -1)
			SUB		[SP+1] ;вычитает 32767, что даст отрицательное значение от -49151 до -32766. Но в UAC отобразится от 16385 до 32768.
					;если теперь найти обратную величину, это должно получиться в диапазоне 32768..65532
					;а если мы вычтем ещё половинку, может выйти от 16384 до 32768, и диапазон 32768..65536.
					;так что ну её нахрен, эту половинку?
					;или можем уже на этапе Ньютона этот случай проверить.
			DIV2S		1 ;это для нужд округления. 					
			[SP]		UAC							
			OUT		j	;для отладки
	FindScale	endp


И тут можно слегка сэкономить. Мы ввели две константы, веса, с которыми нужно взять Манхэттенскую и Чебышевскую метрику, чтобы из них получить что-то напоминающее Эвклидову (см. как быстро оценить длину вектора):

MetricW:
ManhW	dw		11010	;вес, с которым взять Манхэттенскую метрику ;
ChebW	dw		21758	;вес Чебышевской метрики					;
;с такими весами получаем недурственное приближение к Эвклидовой метрике


В более старой версии программы, обращение к этим весам шло внутри цикла, но потом мы решили сделать Loop Unrolling (благо, как было 4 команды, так и осталось!):

	Z		MetricW		
	C		[SP]
	FMA		[Z+i]
	C		[SP+1]
	FMA		[Z+1]


Подставляем константы напрямую - и выигрываем одно слово кода и 2 слова данных:
	C		[SP]
	FMA		11010	;ManhW
	C		[SP+1]
	FMA		21758	;ChebW


Нахождение вектора

	;состояние регистров:
	;X=AfTransf
	;Y=QuatY
	;Z=MetricW
	;i=k=0
	;j=Exp
	;Inv=1
	;Acc=Scale
	;C = 2*Acc (побоч эффект DIV2S)
	FindVector	proc
				Z		ZeroInv
				k		3
				C		[Z+i]	;наше нулевое приближение к обр. величине
		@@Newton:	Acc		0		;занести двойку
				UFMS		[SP]	;Acc = 2 - a * x
				MULU		UAC		;Acc = Acc * x = x * (2- a*x), т.е итерация метода Ньютона
				C		UAC
				kLOOP		@@Newton			
	;золотой ключик почти у нас в кармане
	;если на входе было 16384, то правильный результат 65536 не влезет в 16 бит, сократившись до нуля. 
	;это мы должны заметить по переполнению, только вот не соображу, оно должно ПОЯВИТЬСЯ, или исчезнуть?
				Z		Exp	
				JGE		@@SkipOflo
				j++		0 			;раз не вмещается в мантиссу, значит добавим экспоненту, а в мантиссу нужно 32768 теперь.Ща найдём...
				ZAcc		MinusOne 	;вместо бесхозной минус двойки
		@@SkipOflo:	[Z+i]		j						
				[Z+1]		UAC
	;и ещё осталось найти Ty, Tz
				i		1
				Y		Rx
				X		Ty
		@@vector:	MULSU		[Y+i]
				[X+i]		Acc
				iLOOP		@@vector
	FindVector	endp


И здесь одна константа затесалась, нулевое приближение для метода Ньютона,

ZeroInv	dw		43648 ;нулевое приближение для метода Ньютона нахождения обратной величины


И, как водится, нам приходится сначала инициализировать регистр, а только потом через него обращаться в память:
	Z		ZeroInv
	k		3
	C		[Z+i]	;наше нулевое приближение к обр. величине


причём инициализируем всегда чуть заранее, чтобы Hazard'ов не возникало, и это ухудшает читаемость кода, т.к между объявлением и использованием появляются интервалы. Что ж, мы знаем, что надо делать:
	C		43648	;наше нулевое приближение к обр. величине


И снова сэкономили одно слово данных и одно слово кода.


Раньше эта программа занимала 189 слов данных и 204 слова кода, теперь на 5 слов данных и на 5 слов кода меньше. Не так уж много. Такова уж специфика этой "вычислительной части" - там мы практически всё время работаем с массивами, и поэтому один раз инициализировать базовые регистры и затем обращаться через хитрючие адреса типа [X+2i+k] выходит очень эффективно. "Одиночных констант", как мы видим, раз-два и обчёлся. Потому и жил себе спокойно со старым модулем Imm.

Но в обработке изображений их, похоже, будет куда больше: это и ширина и высота картинки (на первых порах 1024х720), и пороговое значение яркости, и константы для видеопроцессора, задающие кадровые и строчные синхроимпульсы, и куча другого... Да и с доступом к памяти мы, похоже, стояли на пределе - ещё усложнения - и добираться до некоторых значений стало бы совсем туго.
Tags: ПЛИС, программки, работа, странные девайсы
Subscribe

  • Ремонт лыжных мостиков

    Вернулся с сегодняшнего субботника. Очень продуктивно: отремонтировали все ТРИ мостика! Правда, для этого надо было разделиться, благо народу…

  • Гетто-байк

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

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments