contrib/haskell/src/Hkl/Projects/Sixs.hs

   1 {-# LANGUAGE GADTs             #-}
   2 {-# LANGUAGE OverloadedStrings #-}
   3 {-
   4     Copyright  : Copyright (C) 2014-2019 Synchrotron SOLEIL
   5                                          L'Orme des Merisiers Saint-Aubin
   6                                          BP 48 91192 GIF-sur-YVETTE CEDEX
   7     License    : GPL3+
   8
   9     Maintainer : Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  10     Stability  : Experimental
  11     Portability: GHC only (not tested)
  12 -}
  13 module Hkl.Projects.Sixs
  14     ( main_sixs )
  15         where
  16 import           Control.Concurrent.Async          (mapConcurrently)
  17 import           Control.Monad                     (forM_, forever)
  18 import           Control.Monad.IO.Class            (MonadIO (liftIO))
  19 import           Data.Array.Repa                   (Array, Shape, extent,
  20                                                     listOfShape, size)
  21 import           Data.Array.Repa.Index             (DIM1, DIM2, DIM3, Z)
  22 import           Data.Array.Repa.Repr.ForeignPtr   (F, toForeignPtr)
  23 import           Data.Vector.Storable              (concat, head)
  24 import           Data.Word                         (Word16)
  25 import           Foreign.C.Types                   (CInt (..))
  26 import           Foreign.ForeignPtr                (ForeignPtr, withForeignPtr)
  27 import           Foreign.Marshal.Array             (withArrayLen)
  28 import           Foreign.Ptr                       (Ptr)
  29 import           Foreign.Storable                  (peek)
  30 import           Hkl
  31 import           Numeric.LinearAlgebra             (Matrix)
  32 import           Numeric.Units.Dimensional.NonSI   (angstrom)
  33 import           Numeric.Units.Dimensional.Prelude (Angle, Length, degree,
  34                                                     meter, (*~), (/~))
  35 import           Pipes                             (Pipe, await, each, yield,
  36                                                     (>->))
  37 import           Pipes.Prelude                     (toListM)
  38 import           Pipes.Safe                        (SafeT, runSafeT)
  39 import           Text.Printf                       (printf)
  40
  41 data DataFrame
  42     = DataFrame
  43       Int -- n
  44       Geometry -- geometry
  45       (Matrix Double) -- ub
  46       (Array F DIM2 Word16) -- image
  47
  48 instance Show DataFrame where
  49   show (DataFrame i g m _) = unwords [show i, show g, show m]
  50
  51 class FramesP a where
  52   framesP :: a -> Detector b DIM2 -> Pipe FilePath DataFrame (SafeT IO) ()
  53
  54 data DataFrameHklH5Path
  55   = DataFrameHklH5Path
  56     (Hdf5Path DIM3 Word16) -- Image
  57     (Hdf5Path DIM1 Double) -- Mu
  58     (Hdf5Path DIM1 Double) -- Omega
  59     (Hdf5Path DIM1 Double) -- Delta
  60     (Hdf5Path DIM1 Double) -- Gamma
  61     (Hdf5Path DIM2 Double) -- UB
  62     (Hdf5Path Z Double) -- Wavelength
  63     (Hdf5Path DIM1 Char) -- DiffractometerType
  64
  65 instance FramesP DataFrameHklH5Path where
  66   framesP (DataFrameHklH5Path i m o d g u w t) det = forever $ do
  67     fp <- await
  68     withFileP (openH5 fp) $ \f ->
  69       withHdf5PathP f i $ \i' ->
  70       withHdf5PathP f m $ \m' ->
  71       withHdf5PathP f o $ \o' ->
  72       withHdf5PathP f d $ \d' ->
  73       withHdf5PathP f g $ \g' ->
  74       withHdf5PathP f u $ \u' ->
  75       withHdf5PathP f w $ \w' ->
  76       withHdf5PathP f t $ \_t' -> do
  77       (Just n) <- liftIO $ lenH5Dataspace m'
  78       forM_ [0..n-1] (\j -> yield =<< liftIO
  79                        (do
  80                            mu <- get_position m' j
  81                            omega <- get_position o' j
  82                            delta <- get_position d' j
  83                            gamma <- get_position g' j
  84                            wavelength <- get_position w' 0
  85                            image <- get_image det i' j
  86                            ub <- get_ub u'
  87                            let positions = Data.Vector.Storable.concat [mu, omega, delta, gamma]
  88                                source = Source (Data.Vector.Storable.head wavelength *~ angstrom)
  89                            pure $ DataFrame j (Geometry Uhv source positions Nothing) ub image))
  90
  91 -- | DataFrameSpace
  92
  93 data DataFrameSpace sh = DataFrameSpace DataFrame (Space sh)
  94   deriving Show
  95
  96 type Resolutions = [Double]
  97
  98 space :: Detector a DIM2 -> Array F DIM3 Double -> Resolutions -> DataFrame -> IO (DataFrameSpace DIM3)
  99 space detector pixels rs df@(DataFrame _ g@(Geometry _ (Source w) _ _) _ub img) = do
 100   let k = 2 * pi / (w /~ angstrom)
 101   let nPixels = size . shape $ detector
 102   let pixelsDims = map toEnum $ listOfShape . extent $ pixels :: [CInt]
 103   withGeometry g $ \geometry ->
 104     withForeignPtr (toForeignPtr pixels) $ \pix ->
 105     withArrayLen rs $ \nr r ->
 106     withArrayLen pixelsDims $ \ndim dims ->
 107     withForeignPtr (toForeignPtr img) $ \i -> do
 108       p <- {-# SCC "hkl_binoculars_space_q" #-} hkl_binoculars_space_q geometry k i (toEnum nPixels) pix (toEnum ndim) dims r (toEnum nr)
 109       s <- peek p
 110       return (DataFrameSpace df s)
 111
 112 -- | Create the Cube
 113
 114 withForeignPtrs :: [ForeignPtr a] -> ([Ptr a] -> IO r) -> IO r
 115 withForeignPtrs []       f = f []
 116 withForeignPtrs (fp:fps) f =
 117   withForeignPtr fp $ \p ->
 118   withForeignPtrs fps $ \ps -> f (p:ps)
 119
 120 mkCube :: Shape sh => [DataFrameSpace sh] -> IO (Cube sh)
 121 mkCube dfs = do
 122   let spaces = [spaceHklPointer s | (DataFrameSpace _ s) <- dfs]
 123   let images = [toForeignPtr img | (DataFrameSpace (DataFrame _ _ _ img) _) <- dfs]
 124   let (DataFrameSpace (DataFrame _ _ _ img) _ )= Prelude.head dfs
 125   let nPixels = size . extent $ img
 126   withForeignPtrs spaces $ \pspaces ->
 127     withForeignPtrs images $ \pimages ->
 128     withArrayLen pspaces $ \nSpaces' spaces' ->
 129     withArrayLen pimages $ \_ images' -> do
 130     peek =<< {-# SCC "hkl_binoculars_cube_new" #-} hkl_binoculars_cube_new (toEnum nSpaces') spaces' (toEnum nPixels) images'
 131
 132 type Template = String
 133
 134 data InputFn = InputFn FilePath
 135              | InputRange Template Int Int
 136
 137
 138 toList :: InputFn -> [FilePath]
 139 toList (InputFn f)           = [f]
 140 toList (InputRange tmpl f t) = [printf tmpl i | i <- [f..t]]
 141
 142 data Input = Input { filename     :: InputFn
 143                    , h5path       :: DataFrameHklH5Path
 144                    , output       :: FilePath
 145                    , resolutions  :: [Double]
 146                    , centralPixel :: (Int, Int)  -- x, y
 147                    , sdd          :: (Length Float)  -- sample to detector distance
 148                    , detrot       :: Angle Float
 149                    }
 150
 151 _manip1 :: Input
 152 _manip1 = Input { filename = InputFn "/nfs/ruche-sixs/sixs-soleil/com-sixs/2015/Shutdown4-5/XpadAu111/align_FLY2_omega_00045.nxs"
 153                 , h5path = DataFrameHklH5Path
 154                            (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "xpad_image")
 155                            (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "UHV_MU")
 156                            (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "UHV_OMEGA")
 157                            (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "UHV_DELTA")
 158                            (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "UHV_GAMMA")
 159                            (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "I14-C-CX2__EX__DIFF-UHV__#1/" $ datasetp "UB")
 160                            (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "Monochromator" $ datasetp "wavelength")
 161                            (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "I14-C-CX2__EX__DIFF-UHV__#1" $ datasetp "type")
 162                 , output = "test1.hdf5"
 163                 , resolutions = [0.0002, 0.0002, 0.0002]
 164                 , centralPixel = (0, 0)
 165                 , sdd = 1 *~ meter
 166                 , detrot = 90 *~ degree
 167                 }
 168
 169 manip2 :: Input
 170 manip2 = Input { filename = InputRange "/nfs/ruche-sixs/sixs-soleil/com-sixs/2019/Run3/FeSCO_Cu111/sample2_ascan_omega_%05d.nxs" 77 93
 171                , h5path = DataFrameHklH5Path
 172                           (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "xpad_image")
 173                           (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "mu")
 174                           (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "omega")
 175                           (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "delta")
 176                           (hdf5p $ grouppat 0 $ groupp "scan_data" $ datasetp "gamma")
 177                           (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "i14-c-cx2-ex-diff-uhv" $ datasetp "UB")
 178                           (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "i14-c-c02-op-mono" $ datasetp "lambda")
 179                           (hdf5p $ grouppat 0 $ groupp "SIXS" $ groupp "i14-c-cx2-ex-diff-uhv" $ datasetp "type")
 180                , output = "test2.hdf5"
 181                , resolutions = [0.003, 0.01, 0.003]
 182                , centralPixel = (352, 112)
 183                , sdd = 1.162 *~ meter
 184                , detrot = 90 *~ degree
 185                }
 186
 187 main_sixs :: IO ()
 188 main_sixs = do
 189   let input = manip2
 190
 191   pixels <- getPixelsCoordinates ImXpadS140 (centralPixel input) (sdd input) (detrot input)
 192   r <- runSafeT $ toListM $
 193       each (toList $ filename input)
 194       >-> framesP (h5path input) ImXpadS140
 195   r' <- mapConcurrently (space ImXpadS140 pixels (resolutions input)) r
 196   c <- mkCube r'
 197   saveHdf5 (output input) c
 198   print c
 199
 200   return ()