binoculars-ng/src/Hkl/Binoculars/Projections/Angles.hs

   1 {-# LANGUAGE DataKinds           #-}
   2 {-# LANGUAGE FlexibleContexts    #-}
   3 {-# LANGUAGE FlexibleInstances   #-}
   4 {-# LANGUAGE OverloadedStrings   #-}
   5 {-# LANGUAGE RecordWildCards     #-}
   6 {-# LANGUAGE ScopedTypeVariables #-}
   7 {-# LANGUAGE TypeFamilies        #-}
   8
   9 {-# OPTIONS_GHC -fno-warn-orphans #-}
  10
  11 {-
  12     Copyright  : Copyright (C) 2014-2024 Synchrotron SOLEIL
  13                                          L'Orme des Merisiers Saint-Aubin
  14                                          BP 48 91192 GIF-sur-YVETTE CEDEX
  15     License    : GPL3+
  16
  17     Maintainer : Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  18     Stability  : Experimental
  19     Portability: GHC only (not tested)
  20 -}
  21
  22 module Hkl.Binoculars.Projections.Angles
  23     ( newAngles
  24     , processAngles
  25     , updateAngles
  26     ) where
  27
  28 import           Control.Concurrent.Async           (mapConcurrently)
  29 import           Control.Monad.Catch                (MonadThrow)
  30 import           Control.Monad.IO.Class             (MonadIO (liftIO), liftIO)
  31 import           Control.Monad.Logger               (MonadLogger, logDebugN,
  32                                                      logInfoN)
  33 import           Control.Monad.Reader               (MonadReader, ask)
  34 import           Data.HashMap.Lazy                  (fromList)
  35 import           Data.Ini                           (Ini (..))
  36 import           Data.Text                          (pack, unpack)
  37 import           Data.Text.IO                       (putStr)
  38 import           Data.Vector.Storable.Mutable       (unsafeWith)
  39 import           Foreign.C.Types                    (CDouble (..))
  40 import           Foreign.ForeignPtr                 (withForeignPtr)
  41 import           Path                               (Abs, Dir, Path)
  42 import           Pipes                              (each, runEffect, (>->))
  43 import           Pipes.Prelude                      (filter, map, tee, toListM)
  44 import           Pipes.Safe                         (runSafeT)
  45 import           Text.Printf                        (printf)
  46
  47 import           Hkl.Binoculars.Common
  48 import           Hkl.Binoculars.Config
  49 import           Hkl.Binoculars.Config.Common
  50 import           Hkl.Binoculars.Pipes
  51 import           Hkl.Binoculars.Projections
  52 import           Hkl.Binoculars.Projections.QCustom
  53 import           Hkl.C.Binoculars
  54 import           Hkl.Detector
  55 import           Hkl.Geometry
  56 import           Hkl.Image
  57 import           Hkl.Repa
  58 import           Hkl.Utils
  59
  60 ------------
  61 -- Config --
  62 ------------
  63
  64 default'BinocularsConfig'Angles :: Config 'AnglesProjection
  65 default'BinocularsConfig'Angles
  66   = BinocularsConfig'Angles
  67     { binocularsConfig'Angles'Common = default'BinocularsConfig'Common
  68     , binocularsConfig'Angles'ProjectionType = AnglesProjection
  69     , binocularsConfig'Angles'ProjectionResolution = Resolutions3 1 1 1
  70     , binocularsConfig'Angles'ProjectionLimits = Nothing
  71     , binocularsConfig'Angles'DataPath = default'DataSourcePath'DataFrameQCustom
  72     , binocularsConfig'Angles'SampleAxis = SampleAxis "omega"
  73     }
  74
  75 instance HasIniConfig 'AnglesProjection where
  76   data instance Config 'AnglesProjection
  77       = BinocularsConfig'Angles
  78         { binocularsConfig'Angles'Common                 :: BinocularsConfig'Common
  79         , binocularsConfig'Angles'ProjectionType         :: ProjectionType
  80         , binocularsConfig'Angles'ProjectionResolution   :: Resolutions DIM3
  81         , binocularsConfig'Angles'ProjectionLimits       :: Maybe (RLimits DIM3)
  82         , binocularsConfig'Angles'SampleAxis             :: SampleAxis
  83         , binocularsConfig'Angles'DataPath               :: DataSourcePath DataFrameQCustom
  84         } deriving (Show)
  85
  86   newtype Args 'AnglesProjection
  87       = Args'AnglesProjection (Maybe ConfigRange)
  88
  89   getConfig content@(ConfigContent cfg) (Args'AnglesProjection mr) capabilities
  90       = do binocularsConfig'Angles'Common <- parse'BinocularsConfig'Common cfg mr capabilities
  91            binocularsConfig'Angles'ProjectionType <- parseFDef cfg "projection" "type" (binocularsConfig'Angles'ProjectionType default'BinocularsConfig'Angles)
  92            binocularsConfig'Angles'ProjectionResolution <- parseFDef cfg "projection" "resolution" (binocularsConfig'Angles'ProjectionResolution default'BinocularsConfig'Angles)
  93            binocularsConfig'Angles'ProjectionLimits <- parseMb cfg "projection" "limits"
  94            binocularsConfig'Angles'SampleAxis <- parseFDef cfg "input" "sample_axis" $ case binocularsConfig'Angles'ProjectionType of
  95                                                                                         AnglesProjection   -> SampleAxis "omega"
  96                                                                                         Angles2Projection  -> SampleAxis "mu"
  97                                                                                         HklProjection      -> undefined
  98                                                                                         QCustomProjection  -> undefined
  99                                                                                         QIndexProjection   -> undefined
 100                                                                                         QparQperProjection -> undefined
 101                                                                                         QxQyQzProjection   -> undefined
 102                                                                                         RealSpaceProjection -> undefined
 103                                                                                         PixelsProjection -> undefined
 104                                                                                         TestProjection  -> undefined
 105            binocularsConfig'Angles'DataPath <- pure $ eitherF (const $ guess'DataSourcePath'DataFrameQCustom binocularsConfig'Angles'Common Nothing content) (parse' cfg "input" "datapath")
 106                                               (\md -> case md of
 107                                                        Nothing -> guess'DataSourcePath'DataFrameQCustom binocularsConfig'Angles'Common Nothing content
 108                                                        Just d  ->  overload'DataSourcePath'DataFrameQCustom binocularsConfig'Angles'Common Nothing d)
 109            pure BinocularsConfig'Angles{..}
 110
 111 instance ToIni (Config 'AnglesProjection) where
 112   toIni c = toIni (binocularsConfig'Angles'Common c)
 113             `mergeIni`
 114             Ini { iniSections = fromList [ ("input", elemFDef' "datapath" binocularsConfig'Angles'DataPath c default'BinocularsConfig'Angles
 115                                                      <> elemFDef "sample_axis" binocularsConfig'Angles'SampleAxis c default'BinocularsConfig'Angles
 116                                                      [ "the name of the sample axis"
 117                                                      , ""
 118                                                      , "default value: `omega`"
 119                                                      ]
 120                                             )
 121                                          , ("projection", elemFDef' "type" binocularsConfig'Angles'ProjectionType c default'BinocularsConfig'Angles
 122                                                           <> elemFDef' "resolution" binocularsConfig'Angles'ProjectionResolution c default'BinocularsConfig'Angles
 123                                                           <> elemFMbDef' "limits" binocularsConfig'Angles'ProjectionLimits c default'BinocularsConfig'Angles
 124                                            )]
 125                 , iniGlobals = []
 126                 }
 127
 128 -------------------------
 129 -- Angles Projection --
 130 -------------------------
 131
 132 {-# INLINE spaceAngles #-}
 133 spaceAngles :: Detector a DIM2 -> Array F DIM3 Double -> Resolutions DIM3 -> Maybe Mask -> Maybe (RLimits DIM3) -> SampleAxis -> Space DIM2 -> DataFrameQCustom -> IO (DataFrameSpace DIM2)
 134 spaceAngles det pixels rs mmask' mlimits sAxis space@(Space fSpace) (DataFrameQCustom att g img _ _) =
 135   withNPixels det $ \nPixels ->
 136   withGeometry g $ \geometry ->
 137   withForeignPtr (toForeignPtr pixels) $ \pix ->
 138   withResolutions rs $ \nr r ->
 139   withPixelsDims pixels $ \ndim dims ->
 140   withMaybeMask mmask' $ \ mask'' ->
 141   withMaybeLimits mlimits rs $ \nlimits limits ->
 142   withSampleAxis sAxis $ \sampleAxis ->
 143   withForeignPtr fSpace $ \pSpace -> do
 144   case img of
 145     (ImageInt32 arr) -> unsafeWith arr $ \i -> do
 146       {-# SCC "hkl_binoculars_space_angles_int32_t" #-} c'hkl_binoculars_space_angles_int32_t pSpace geometry i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) sampleAxis
 147     (ImageWord16 arr) -> unsafeWith arr $ \i -> do
 148       {-# SCC "hkl_binoculars_space_angles_uint16_t" #-} c'hkl_binoculars_space_angles_uint16_t pSpace geometry i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) sampleAxis
 149     (ImageWord32 arr) -> unsafeWith arr $ \i -> do
 150       {-# SCC "hkl_binoculars_space_angles_uint32_t" #-} c'hkl_binoculars_space_angles_uint32_t pSpace geometry i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) sampleAxis
 151
 152   return (DataFrameSpace img space att)
 153
 154 ----------
 155 -- Pipe --
 156 ----------
 157
 158 processAnglesP :: (MonadIO m, MonadLogger m, MonadReader (Config 'AnglesProjection) m, MonadThrow m)
 159                => m ()
 160 processAnglesP = do
 161   conf :: Config 'AnglesProjection <- ask
 162
 163   -- directly from the common config
 164   let common = binocularsConfig'Angles'Common conf
 165
 166   let overwrite = binocularsConfig'Common'Overwrite common
 167   let det = binocularsConfig'Common'Detector common
 168   let (NCores cap) =  binocularsConfig'Common'NCores common
 169   let destination = binocularsConfig'Common'Destination common
 170   let centralPixel' = binocularsConfig'Common'Centralpixel common
 171   let (Meter sampleDetectorDistance) = binocularsConfig'Common'Sdd common
 172   let (Degree detrot) = binocularsConfig'Common'Detrot common
 173   let mImageSumMax = binocularsConfig'Common'ImageSumMax common
 174   let inputRange = binocularsConfig'Common'InputRange common
 175   let nexusDir = binocularsConfig'Common'Nexusdir common
 176   let tmpl = binocularsConfig'Common'Tmpl common
 177   let maskMatrix = binocularsConfig'Common'Maskmatrix common
 178   let mSkipFirstPoints = binocularsConfig'Common'SkipFirstPoints common
 179   let mSkipLastPoints = binocularsConfig'Common'SkipLastPoints common
 180
 181   -- directly from the specific config
 182   let mlimits = binocularsConfig'Angles'ProjectionLimits conf
 183   let res = binocularsConfig'Angles'ProjectionResolution conf
 184   let datapaths = binocularsConfig'Angles'DataPath conf
 185   let projectionType = binocularsConfig'Angles'ProjectionType conf
 186   let sampleAxis = binocularsConfig'Angles'SampleAxis conf
 187
 188   -- built from the config
 189   output' <- liftIO $ destination' projectionType Nothing inputRange mlimits destination overwrite
 190   filenames <- InputFn'List <$> files nexusDir (Just inputRange) tmpl
 191   mask' <- getMask maskMatrix det
 192   pixels <- liftIO $ getPixelsCoordinates det centralPixel' sampleDetectorDistance detrot Normalisation
 193   let fns = concatMap (replicate 1) (toList filenames)
 194   chunks <- liftIO $ runSafeT $ toListM $ each fns >-> chunkP mSkipFirstPoints mSkipLastPoints datapaths
 195   let ntot = sum (Prelude.map clength chunks)
 196   let jobs = chunk (quot ntot cap) chunks
 197
 198   -- log parameters
 199
 200   logDebugNSH filenames
 201   logDebugNSH datapaths
 202   logDebugNSH chunks
 203   logDebugN "start gessing final cube size"
 204
 205   -- guess the final cube dimensions (To optimize, do not create the cube, just extract the shape)
 206
 207   guessed <- liftIO $ withCubeAccumulator EmptyCube $ \c ->
 208     runSafeT $ runEffect $
 209     each chunks
 210     >-> Pipes.Prelude.map (\(Chunk fn f t) -> (fn, [f, quot (f + t) 4, quot (f + t) 4 * 2, quot (f + t) 4 * 3, t]))
 211     >-> framesP datapaths
 212     >-> project det 3 (spaceAngles det pixels res mask' mlimits sampleAxis)
 213     >-> accumulateP c
 214
 215   logDebugN "stop gessing final cube size"
 216
 217   -- do the final projection
 218
 219   logInfoN (pack $ printf "let's do a Angles projection of %d %s image(s) on %d core(s)" ntot (show det) cap)
 220
 221   liftIO $ withProgressBar ntot $ \pb -> do
 222     r' <- mapConcurrently (\job -> withCubeAccumulator guessed $ \c ->
 223                              runSafeT $ runEffect $
 224                              each job
 225                              >-> Pipes.Prelude.map (\(Chunk fn f t) -> (fn, [f..t]))
 226                              >-> framesP datapaths
 227                              >-> Pipes.Prelude.filter (\(DataFrameQCustom _ _ img _ _) -> filterSumImage mImageSumMax img)
 228                              >-> project det 3 (spaceAngles det pixels res mask' mlimits sampleAxis)
 229                              >-> tee (accumulateP c)
 230                              >-> progress pb
 231                          ) jobs
 232     saveCube output' (unpack . serializeConfig $ conf) r'
 233
 234 ---------
 235 -- Cmd --
 236 ---------
 237
 238 processAngles :: (MonadLogger m, MonadThrow m, MonadIO m) => Maybe FilePath -> Maybe ConfigRange -> m ()
 239 processAngles mf mr = cmd processAnglesP mf (Args'AnglesProjection mr)
 240
 241 newAngles :: (MonadIO m, MonadLogger m, MonadThrow m)
 242           => Path Abs Dir -> m ()
 243 newAngles cwd = do
 244   let conf = default'BinocularsConfig'Angles
 245              { binocularsConfig'Angles'Common = default'BinocularsConfig'Common
 246                                                 { binocularsConfig'Common'Nexusdir = Just cwd }
 247              }
 248   liftIO $ Data.Text.IO.putStr $ serializeConfig conf
 249
 250 updateAngles :: (MonadIO m, MonadLogger m, MonadThrow m)
 251              => Maybe FilePath -> Maybe ConfigRange -> m ()
 252 updateAngles mf mr = cmd (pure ()) mf (Args'AnglesProjection mr)