haskell - 四元数旋转有一个奇怪的行为(Haskell OpenGL)

Question

我一直在关注 Haskell OpenGL tutorial .3D 空间中的旋转引起了我的兴趣，所以我开始学习欧拉角，最后学习四元数。

我想使用四元数实现我自己的函数来执行旋转(在立方体上)，我基于这两篇论文:mostly this one和 this one .

当我只在一个轴上执行旋转时，我的函数工作正常，但是当我在 X 和 Y 上执行旋转时，例如，立方体开始随机前进并在旋转时被“阻挡”。

Video of the cube performing rotation on XY .

当我设置三个轴(X、Y、Z)时，它会放大得更多(但没有那个奇怪的阻挡物):video .

这是我的程序代码:

这里是主文件，它创建一个窗口，设置空闲函数并在屏幕上输出旋转角度 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

这是定义旋转函数的核心文件(还有一些其他函数)。我添加了一些评论，因为它可能是一些低质量的 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)]

最后，如果它可以帮助人们查看我的算法是否存在问题，这里有一些使用我的函数的旋转示例:

X 轴上的点 (1, 2, 3) 围绕点 (0, 0, 0)(原点)旋转 90° 得到:(0.99999994,-3.0,2.0)

相同的旋转但在 X 和 Y 轴上给出:(5.4999995,-0.99999994,-0.49999988)

同样的旋转，但在 X、Y 和 Z 轴上给出:(5.9999995,1.9999999,3.9999995)

Answer 1

您指向的关于四元数旋转的第二篇论文有这样一句话:

“(x̂, ŷ, ẑ) 是定义旋转轴的单位向量。”。

所以四元数必须归一化，分量的平方和等于1。

例如，如果您涉及所有 3 个轴，则它必须是 (cos θ/2, r3sin θ/2, r3sin θ/2, r3*sin θ/2)其中 r3 是 3 的平方根的倒数。这就是我如何解释您在帖子末尾提到的旋转结果在涉及多个轴时无法保持矢量长度的原因。

因此，关键部分是函数 rotate3f 中的这一行:

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

缺少归一化因子。

您的代码提供了许多提高可读性的机会。您可以考虑使用 CodeReview了解更多详情。

一个主要问题是源代码行太宽。如果读者必须使用水平 slider ，那么理解代码和发现错误就困难得多。下面，我将尽量避免超过 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)

请注意，定义特殊类型的想法不仅可以减少代码宽度，还可以提供额外的灵 active 。如果有一天您决定用比普通列表更有效的其他数据结构来表示四元数，则可以在保持客户端代码不变的情况下完成。

接下来，正确的旋转代码。函数 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 基础设施和重现动画。请让我们知道这是否解决了整个问题。