Вращение кватернионов имеет странное поведение (Haskell OpenGL) - PullRequest
4 голосов
/ 11 июля 2020

Я следил за Haskell учебником OpenGL . Вращения в трехмерном пространстве заинтриговали меня, поэтому я начал изучать углы Эйлера и, наконец, кватернионы. эти два документа: в основном этот и этот .

Моя функция работает нормально, когда я выполняю вращение только на одной оси, но когда я это делаю например, по X и Y куб начинает случайным образом go вперед и блокируется при вращении.

Видео куба, выполняющего вращение по XY .

Когда я устанавливаю три оси (X, Y, Z), масштабирование увеличивается еще больше (но без этой странной блокирующей штуки): видео .

Вот код моей программы:

Вот главный файл, который создает окно, устанавливает функцию ожидания и выводит результат поворота на угол A на экране, где A увеличивается на 0,05 для каждого кадра.

module Main (main) where
import Core
import Utils
import Data.IORef
import Graphics.UI.GLUT
import Graphics.Rendering.OpenGL

main :: IO ()
main = do
    createAWindow "177013"
    mainLoop

createAWindow :: [Char] -> IO ()
createAWindow windowName = do
    (procName, _args) <- getArgsAndInitialize
    createWindow windowName
    initialDisplayMode $= [DoubleBuffered]
    angle <- newIORef 0.0
    delta <- newIORef 0.05
    displayCallback $= (start angle)
    reshapeCallback $= Just reshape
    keyboardMouseCallback $= Just keyboardMouse
    idleCallback $= Just (idle angle delta)

reshape :: ReshapeCallback
reshape size = do
             viewport $= (Position 0 0, size)
             postRedisplay Nothing


keyboardMouse :: KeyboardMouseCallback
keyboardMouse _ _ _ _ = return ()

idle :: IORef GLfloat -> IORef GLfloat -> IdleCallback
idle angle delta = do
           d <- get delta
           a <- get angle
           angle $~! (+d)
           postRedisplay Nothing
start :: IORef GLfloat -> DisplayCallback
start angle = do
            clear [ColorBuffer]
            loadIdentity
            a <- get angle
            let c = rotate3f (0, 0, 0) [X,Y,Z] a $ cube3f 0.2 -- here I'm rotating on X, Y and Z axis
            draw3f Quads c CCyan
            flush
            swapBuffers
                where

Вот файл ядра, в котором ion функция определена (с некоторыми другими). Я добавил несколько комментариев, так как это, вероятно, код низкого качества haskell.

module Core (draw3f, vertex3f, rotate3f, translate3f, rotate3d, Colors(..), Axes(..)) where

import Control.Lens
import Graphics.Rendering.OpenGL

data Axes = X | Y | Z
            deriving Eq
data Colors = CRed | CGreen | CBlue | CYellow | CWhite | CMagenta | CCyan | CBlack | CNone | CPreset
              deriving Eq


rotate3f :: (GLfloat, GLfloat, GLfloat) -> [Axes] -> GLfloat -> [(GLfloat, GLfloat, GLfloat)] -> [(GLfloat, GLfloat, GLfloat)]
rotate3f _ _ _ [] = []
rotate3f _ [] _ _ = []
rotate3f o axes a p = let p' = translate3f p u -- translation if I don't want to rotate it by the origin
                          q = cos a' : ((\x -> if x `elem` axes then sin a' else 0) <$> [X,Y,Z]) -- if the axe is set then its related component is equal to sin theta/2, otherwise it will be 0
                          q' = q !! 0 : (negate <$> (tail q)) -- quaternion inversion
                      in translate3f ((rotate q q') <$> p') [(0,0,0),o] -- rotate and translate again to put the object where it belongs
                          where
                              a' = (a * (pi / 180)) / 2 -- convert to radians and divide by 2 as all q components takes theta/2
                              u :: [(GLfloat, GLfloat, GLfloat)]
                              u = [o,(0,0,0)]
                              rotate :: [GLfloat] -> [GLfloat] -> (GLfloat, GLfloat, GLfloat) -> (GLfloat, GLfloat, GLfloat)
                              rotate q q' (x,y,z) = let p = [0,x,y,z]
                                                        qmul q1 q2 = [(q1 !! 0) * (q2 !! 0) - (q1 !! 1) * (q2 !! 1) - (q1 !! 2) * (q2 !! 2) - (q1 !! 3) * (q2 !! 3),
                                                                      (q1 !! 0) * (q2 !! 1) + (q1 !! 1) * (q2 !! 0) + (q1 !! 2) * (q2 !! 3) - (q1 !! 3) * (q2 !! 2),
                                                                      (q1 !! 0) * (q2 !! 2) - (q1 !! 1) * (q2 !! 3) + (q1 !! 2) * (q2 !! 0) + (q1 !! 3) * (q2 !! 1),
                                                                      (q1 !! 0) * (q2 !! 3) + (q1 !! 1) * (q2 !! 2) - (q1 !! 2) * (q2 !! 1) + (q1 !! 3) * (q2 !! 0)]
                                                        p' = qmul (qmul q p) q'
                                                    in (p' !! 1, p' !! 2, p' !! 3)


                    
translate3f :: [(GLfloat, GLfloat, GLfloat)] -> [(GLfloat, GLfloat, GLfloat)] -> [(GLfloat, GLfloat, GLfloat)]
translate3f p [(ax,ay,az),(bx,by,bz)] = map (\(x,y,z) -> (x + (bx - ax), y + (by - ay), z + (bz - az))) p



draw3f :: PrimitiveMode -> [(GLfloat, GLfloat, GLfloat)] -> Colors -> IO()
draw3f shape points color = renderPrimitive shape $ mapM_ (\(x,y,z) -> vertex3f x y z color) points

vertex3f :: GLfloat -> GLfloat -> GLfloat -> Colors -> IO()
vertex3f x y z c = do
                 if c /= CPreset
                    then color $ Color3 (c' ^. _1) (c' ^. _2) ((c' ^. _3) :: GLfloat)
                 else return ()
                 vertex $ Vertex3 x y z
                     where
                         c' :: (GLfloat, GLfloat, GLfloat)
                         c' = case c of CRed -> (1,0,0)
                                        CGreen -> (0,1,0)
                                        CBlue -> (0,0,1)
                                        CYellow -> (1,1,0)
                                        CMagenta -> (1,0,1)
                                        CCyan -> (0,1,1)
                                        CBlack -> (0,0,0)
                                        _ -> (1,1,1)

А вот файл utils, где есть только определение куба из Haskell Учебник по OpenGL

module Utils (cube3f) where

import Core
import Graphics.UI.GLUT
import Graphics.Rendering.OpenGL

cube3f :: GLfloat -> [(GLfloat, GLfloat, GLfloat)]
cube3f w = [( w, w, w), ( w, w,-w), ( w,-w,-w), ( w,-w, w),
            ( w, w, w), ( w, w,-w), (-w, w,-w), (-w, w, w),
            ( w, w, w), ( w,-w, w), (-w,-w, w), (-w, w, w),
            (-w, w, w), (-w, w,-w), (-w,-w,-w), (-w,-w, w),
            ( w,-w, w), ( w,-w,-w), (-w,-w,-w), (-w,-w, w),
            ( w, w,-w), ( w,-w,-w), (-w,-w,-w), (-w, w,-w)]

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

Вращение на 90 ° точки (1, 2, 3) по оси X вокруг точки (0, 0, 0) (начало координат) дает: (0.99999994,-3.0,2.0)

То же вращение, но по осям X и Y дает: (5.4999995,-0.99999994,-0.49999988)

Снова то же вращение, но по осям X, Y и Z дает: (5.9999995,1.9999999,3.9999995)

1 Ответ

2 голосов
/ 12 июля 2020

Второй документ о вращении кватернионами, на который вы указываете, имеет следующее предложение:

«(x̂, ŷ, ẑ) - единичный вектор, определяющий ось вращения». .

Итак, кватернион должен быть нормализован, сумма компонентов в квадрате равна 1.

Так, например, если у вас задействованы все 3 оси, это должно быть (cos θ / 2, r3 sin θ / 2, r3 sin θ / 2, r3 * sin θ / 2), где r3 - величина, обратная квадрату root числа 3. Вот как я бы объяснил, что вращение результаты, которые вы упомянули в конце своего сообщения, не сохраняют длину вектора, когда задействовано несколько осей.

Таким образом, критически важной частью является эта строка в функции rotate3f:

q = cos a' : ((\x -> if x `elem` axes then sin a' else 0) <$> [X,Y,Z])

где отсутствует коэффициент нормализации.

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

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

Во-первых, нам нужна некоторая кватернионная инфраструктура:

{-#  LANGUAGE  ScopedTypeVariables  #-}
{-#  LANGUAGE  ExplicitForAll       #-}

type GLfloat   = Float
type GLfloatV3 = (GLfloat, GLfloat, GLfloat)
type QuatFloat = [GLfloat]

data Axes =  X | Y | Z  deriving  Eq

qmul :: QuatFloat -> QuatFloat -> QuatFloat
qmul  [qa0, qa1, qa2, qa3]  [qb0, qb1, qb2, qb3] =
    [
       qa0*qb0 - qa1*qb1 - qa2*qb2 - qa3*qb3 ,
       qa0*qb1 + qa1*qb0 + qa2*qb3 - qa3*qb2 ,
       qa0*qb2 - qa1*qb3 + qa2*qb0 + qa3*qb1 ,
       qa0*qb3 + qa1*qb2 - qa2*qb1 + qa3*qb0
    ]
qmul _ _  =  error "Quaternion length differs from 4"

qconj :: QuatFloat -> QuatFloat
qconj q = (head q) : (map negate (tail q)) -- q-conjugation

rotate :: [GLfloat] -> [GLfloat] -> GLfloatV3 -> GLfloatV3
rotate q q' (x,y,z) = let  p             = [0, x,y,z]
                           [q0,q1,q2,q3] = qmul (qmul q p) q'
                      in  (q1, q2, q3)

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

Затем, собственно код вращения. Функция rotQuat0 - ваш исходный алгоритм, который воспроизводит числовые результаты, указанные в конце вашего вопроса. Функция rotQuat1 - это модифицированная версия, дающая 1-нормализованный кватернион.

-- original code:
rotQuat0 :: [Axes] -> GLfloat -> QuatFloat
rotQuat0 axes angle = let  fn x = if (x `elem` axes) then (sin angle) else 0
                      in   (cos angle) : (map fn [X,Y,Z])

-- modified code:
rotQuat1 :: [Axes] -> GLfloat -> QuatFloat
rotQuat1 axes angle = let  corr = 1.0 / sqrt (fromIntegral (length axes))
                           fn x = if (x `elem` axes) then corr*(sin angle) else 0
                      in   (cos angle) : (map fn [X,Y,Z])

Код с использованием rotQuat1:

rotate3f :: GLfloatV3 -> [Axes] -> GLfloat -> [GLfloatV3] -> [GLfloatV3]
rotate3f _ _ _ [] = []
rotate3f _ [] _ _ = []
rotate3f org axes degθ pts =
    let   -- convert to radians and divide by 2, as all q components take θ/2
          a' = (degθ * (pi / 180)) / 2
          u :: [GLfloatV3]
          u = [org, (0,0,0)]
          -- translation if I don't want to rotate it by the origin
          p' = translate3f pts u
          -- if the axis is set, then its related component is
          -- equal to sin θ/2, otherwise it will be zero
          ---- q = cos a' : ((\x -> if x `elem` axes then sin a' else 0) <$> [X,Y,Z])
          q = rotQuat1 axes a'  -- modified version
          q' = qconj q
         -- rotate and translate again to put the object where it belongs
    in   translate3f ((rotate q q') <$> p') [(0,0,0), org] 

             
translate3f :: [GLfloatV3] -> [GLfloatV3] -> [GLfloatV3]
translate3f  pts  [(ax,ay,az), (bx,by,bz)]  =
    let   dx = bx - ax
          dy = by - ay
          dz = bz - az
    in   map  (\(x,y,z) -> (x + dx, y + dy, z + dz))  pts

Код тестирования:

sqNorm3 :: GLfloatV3 -> GLfloat
sqNorm3 (x,y,z) = x*x + y*y +z*z

printAsLines :: Show α => [α] -> IO ()
printAsLines xs = mapM_  (putStrLn . show)  xs

main = do
    let  pt  = (1,2,3) :: GLfloatV3
         pt1 = rotate3f (0,0,0) [X]     90 [pt]
         pt2 = rotate3f (0,0,0) [X,Y]   90 [pt]
         pt3 = rotate3f (0,0,0) [X,Y,Z] 90 [pt]
         pts = map head [pt1, pt2, pt3]
         ptN = map sqNorm3 pts
    printAsLines pts
    putStrLn " "
    printAsLines ptN

Давайте проверим, что с помощью функции rotQuat1 квадрат нормы вашего исходного (1,2,3) входного вектора (то есть 1 + 4 + 9 = 13) остается неизменным, как и положено правильному вращению:

$ ghc opengl00.hs -o ./opengl00.x && ./opengl00.x
[1 of 1] Compiling Main             ( opengl00.hs, opengl00.o )
Linking ./opengl00.x ...

(0.99999994,-3.0,2.0)
(3.6213198,-0.62132025,0.70710695)
(2.5773501,0.84529924,2.5773501)

14.0
13.999995
13.999998
$ 

К сожалению, у меня нет времени, чтобы установить инфраструктуру OpenGL и воспроизвести анимацию. Сообщите нам, полностью ли это исправит.

...