使用Haskell修復陣列實現相位解纏算法

Question

我正在嘗試使用修復陣列實現Haskell中三相結構光掃描的相位解纏算法。我想實現從點（寬度/ 2，高度/ 2）向外遞歸的基於洪泛填充的解包算法。不幸的是，使用這種遞歸方法，我得到了一個內存不足的例外。我是Haskell和Repa庫的新手，所以我想知道它是否看起來像我做任何明顯錯誤的事情。任何幫助，將不勝感激！

更新（@leventov）：

我現在正在考慮實施以下使用可變數組在Yarr算法如下路徑。 （出版物：K. Chen，J. Xi，Y. Yu & JF Chicharo，「用於三維條紋圖形輪廓術的快速質量引導填充相位解纏算法」，在工業應用的光學計量和檢驗中， 2010， 1-9。）

{-# OPTIONS_GHC -Odph -rtsopts -fno-liberate-case -fllvm -optlo-O3 -XTypeOperators -XNoMonomorphismRestriction #-} 

module Scanner where 

import Data.Word 
import Data.Fixed 
import Data.Array.Repa.Eval 
import qualified Data.Array.Repa as R 
import qualified Data.Array.Repa.Repr.Unboxed as U 
import qualified Data.Array.Repa.Repr.ForeignPtr as P 
import Codec.BMP 
import Data.Array.Repa.IO.BMP 
import Control.Monad.Identity (runIdentity) 
import System.Environment(getArgs) 

type ImRead = Either Error Image 
type Avg = P.Array R.U R.DIM2 (ImageT, ImageT, ImageT) 
type ImageT = (Word8, Word8, Word8) 
type PhaseT = (Float, Float, Float) 
type WrapT = (Float, Int) 
type Image = P.Array R.U R.DIM2 (Word8, Word8, Word8) 
type Phase = P.Array R.U R.DIM2 (Float, Float, Float) 
type Wrap = P.Array R.U R.DIM2 (Float, Int) 
type UWrapT = (Float, Int, [(Int, Int)], String) 
type DepthT = (Float, Int, String) 

{-# INLINE noise #-} 
{-# INLINE zskew #-} 
{-# INLINE zscale #-} 
{-# INLINE compute #-} 
{-# INLINE main #-} 
{-# INLINE doMain #-} 
{-# INLINE zipImg #-} 
{-# INLINE mapWrap #-} 
{-# INLINE avgPhase #-} 
{-# INLINE doAvg #-} 
{-# INLINE doWrap #-} 
{-# INLINE doPhase #-} 
{-# INLINE isPhase #-} 
{-# INLINE diffPhase #-} 
{-# INLINE shape #-} 
{-# INLINE countM #-} 
{-# INLINE inArr #-} 
{-# INLINE idx #-} 
{-# INLINE getElem #-} 
{-# INLINE start #-} 
{-# INLINE unwrap #-} 
{-# INLINE doUnwrap #-} 
{-# INLINE doDepth #-} 
{-# INLINE write #-} 

noise :: Float 
noise = 0.1 

zskew :: Float 
zskew = 24 

zscale :: Float 
zscale = 130 

compute :: (R.Shape sh, U.Unbox e) => P.Array R.D sh e -> P.Array R.U sh e 
compute a = runIdentity (R.computeP a) 

main :: IO() 
main = do 
     commandArguments <- getArgs 
     case commandArguments of 
      (file1 : file2 : file3 : _) -> do 
       image1 <- readImageFromBMP file1 
       image2 <- readImageFromBMP file2 
       image3 <- readImageFromBMP file3 
       doMain image1 image2 image3            
      _ -> putStrLn "Not enough arguments" 

doMain :: ImRead -> ImRead -> ImRead -> IO() 
doMain (Right i1) (Right i2) (Right i3) = write 
     where 
      write = writeFile "out.txt" str 
      (p, m, d, str) = start $ mapWrap i1 i2 i3 
doMain _ _ _ = putStrLn "Error loading image" 

zipImg :: Image -> Image -> Image -> Avg 
zipImg i1 i2 i3 = U.zip3 i1 i2 i3 

mapWrap :: Image -> Image -> Image -> Wrap 
mapWrap i1 i2 i3 = compute $ R.map wrap avg 
     where 
      wrap = (doWrap . avgPhase) 
      avg = zipImg i1 i2 i3 

avgPhase :: (ImageT, ImageT, ImageT) -> PhaseT 
avgPhase (i1, i2, i3) = (doAvg i1, doAvg i2, doAvg i3) 

doAvg :: ImageT -> Float 
doAvg (r, g, b) = (r1 + g1 + b1)/d1 
     where 
      r1 = fromIntegral r 
      g1 = fromIntegral g 
      b1 = fromIntegral b 
      d1 = fromIntegral 765 

doWrap :: PhaseT -> WrapT 
doWrap (p1, p2, p3) = (wrap, mask) 
     where 
      wrap = isPhase $ doPhase (p1, p2, p3) 
      mask = isNoise $ diffPhase [p1, p2, p3] 

doPhase :: PhaseT -> (Float, Float) 
doPhase (p1, p2, p3) = (x1, x2) 
     where 
      x1 = sqrt 3 * (p1 - p3) 
      x2 = 2 * p2 - p1 - p3 

isPhase :: (Float, Float) -> Float 
isPhase (x1, x2) = atan2 x1 x2/(2 * pi) 

diffPhase :: [Float] -> Float 
diffPhase phases = maximum phases - minimum phases 

isNoise :: Float -> Int 
isNoise phase = fromEnum $ phase <= noise 

shape :: Wrap -> [Int] 
shape wrap = R.listOfShape $ R.extent wrap 

countM :: Wrap -> (Float, Int) 
countM wrap = R.foldAllS count (0,0) wrap 
     where count = (\(x, y) (i, j) -> (x, y)) 

start :: Wrap -> UWrapT 
start wrap = unwrap wrap (x, y) (ph, m, [], "") 
     where 
      [x0, y0] = shape wrap 
      x  = quot x0 2 
      y  = quot y0 2 
      (ph, m) = getElem wrap (x0, y0) 

inArr :: Wrap -> (Int, Int) -> Bool 
inArr wrap (x,y) = x >= 0 && y >= 0 && x < x0 && y < y0 
     where 
      [x0, y0] = shape wrap 

idx :: (Int, Int) -> (R.Z R.:. Int R.:. Int) 
idx (x, y) = (R.Z R.:. x R.:. y) 

getElem :: Wrap -> (Int, Int) -> WrapT 
getElem wrap (x, y) = wrap R.! idx (x, y) 

unwrap :: Wrap -> (Int, Int) -> UWrapT -> UWrapT 
unwrap wrap (x, y) (ph, m, done, str) = 
     if 
      not $ inArr wrap (x, y) || 
      (x, y) `elem` done || 
      toEnum m::Bool 
     then 
      (ph, m, done, str) 
     else 
      up 
      where 
       unwrap' = doUnwrap wrap (x, y) (ph, m, done, str) 
       right = unwrap wrap (x+1, y) unwrap' 
       left = unwrap wrap (x-1, y) right 
       down = unwrap wrap (x, y+1) left 
       up  = unwrap wrap (x, y-1) down 

doUnwrap :: Wrap -> (Int, Int) -> UWrapT -> UWrapT 
doUnwrap wrap (x, y) (ph, m, done, str) = unwrapped 
     where 
      unwrapped = (nph, m, (x, y):done, out) 
      (phase, mask) = getElem wrap (x, y) 
      rph = fromIntegral $ round ph 
      off = phase - (ph - rph) 
      nph = ph + (mod' (off + 0.5) 1) - 0.5 
      out = doDepth wrap (x, y) (nph, m, str) 

doDepth :: Wrap -> (Int, Int) -> DepthT -> String 
doDepth wrap (x, y) (ph, m, str) = write (x, ys, d, str) 
     where 
      [x0, y0] = shape wrap 
      ys  = y0 - y 
      ydiff = fromIntegral (y - (quot y0 2)) 
      plane = 0.5 - ydiff/zskew 
      d  = (ph - plane) * zscale 


write :: (Int, Int, Float, String) -> String 
write (x, y, depth, str) = str ++ vertex 
     where 
      vertex = xstr ++ ystr ++ zstr 
      xstr = show x ++ " " 
      ystr = show y ++ " " 
      zstr = show depth ++ "\n" 

Answer 1

對不起，我的第一個建議誤導一些浪費你的時間。

，因爲後者需要線性時間應該使用像素狀態的另一個2維陣列（已經訪問過或不）代替

(x, y) `elem` done 

。

解決幾乎相同的任務的示例：對於repa and vector和yarr。

也許，你有內存不足的例外，因爲通過追加到最後（write函數）構建一個字符串 - 最糟糕的解決方案，線性時間和內存消耗。您最好使用cons（:）來彙總結果，並以相反的順序將其寫入輸出文件。甚至更好 - 將結果寫入另一個未裝箱的(Int, Int, Float)元素（分配width*height大小的矢量 - 作爲可能大小的上限）。