contrib/haskell/src/Hkl/XRD/Calibration.hs

   1 {-# LANGUAGE CPP #-}
   2
   3 module Hkl.XRD.Calibration
   4        ( NptExt(..)
   5        , XRDCalibrationEntry(..)
   6        , XRDCalibration(..)
   7        , calibrate
   8        ) where
   9
  10 import Control.Monad.IO.Class (liftIO)
  11 import Data.ByteString.Char8 (pack)
  12 import Data.List (foldl')
  13 import Data.Text (Text)
  14 import Data.Vector.Storable
  15   ( Vector
  16   , head
  17   , concat
  18   , fromList
  19   , slice
  20   , toList
  21   )
  22 import Numeric.LinearAlgebra
  23   ( Matrix
  24   , (<>)
  25   , (#>)
  26   , atIndex
  27   , ident
  28   )
  29 import Numeric.GSL.Minimization
  30   ( MinimizeMethod(NMSimplex2)
  31   , minimizeV
  32   )
  33 import Numeric.Units.Dimensional.Prelude (meter, radian, nano, (/~), (*~))
  34 import Pipes.Safe (MonadSafe(..), runSafeT, bracket)
  35 import Prelude hiding (head, concat, lookup, readFile, writeFile, unlines)
  36
  37 import Hkl.C
  38 import Hkl.Detector
  39 import Hkl.H5
  40 import Hkl.PyFAI
  41 import Hkl.MyMatrix
  42 import Hkl.PyFAI.PoniExt
  43 import Hkl.Types
  44 import Hkl.XRD
  45
  46 -- | Calibration
  47
  48 data NptExt a = NptExt { nptExtNpt :: Npt
  49                        , nptExtMyMatrix :: MyMatrix Double
  50                        , nptExtDetector :: Detector a
  51                        }
  52               deriving (Show)
  53
  54 data XRDCalibrationEntry = XRDCalibrationEntry { xrdCalibrationEntryNxs :: Nxs
  55                                                , xrdCalibrationEntryIdx :: Int
  56                                                , xrdCalibrationEntryNptPath :: FilePath
  57                                                }
  58                            deriving (Show)
  59
  60 data XRDCalibration = XRDCalibration { xrdCalibrationName :: Text
  61                                      , xrdCalibrationOutputDir :: FilePath
  62                                      , xrdCalibrationEntries :: [XRDCalibrationEntry]
  63                                      }
  64                       deriving (Show)
  65
  66 withDataItem :: MonadSafe m => File -> DataItem -> (Dataset -> m r) -> m r
  67 withDataItem hid (DataItem name _) = bracket (liftIO acquire') (liftIO . release')
  68     where
  69       acquire' :: IO Dataset
  70       acquire' = openDataset hid (pack name) Nothing
  71
  72       release' :: Dataset -> IO ()
  73       release' = closeDataset
  74
  75 getM :: File -> DataFrameH5Path -> Int -> IO (MyMatrix Double)
  76 getM f p i' = runSafeT $
  77     withDataItem f (h5pGamma p) $ \g' ->
  78     withDataItem f (h5pDelta p) $ \d' ->
  79     withDataItem f (h5pWavelength p) $ \w' -> liftIO $ do
  80       let mu = 0.0
  81       let komega = 0.0
  82       let kappa = 0.0
  83       let kphi = 0.0
  84       gamma <- get_position g' 0
  85       delta <- get_position d' i'
  86       wavelength <- get_position w' 0
  87       let source = Source (head wavelength *~ nano meter)
  88       let positions = concat [mu, komega, kappa, kphi, gamma, delta]
  89       let geometry = Geometry K6c source positions Nothing
  90       let detector = ZeroD
  91       m <- geometryDetectorRotationGet geometry detector
  92       return (MyMatrix HklB m)
  93
  94 readXRDCalibrationEntry :: Detector a -> XRDCalibrationEntry -> IO (NptExt a)
  95 readXRDCalibrationEntry d e =
  96     withH5File f $ \h5file -> do
  97       m <- getM h5file p idx
  98       npt <- nptFromFile (xrdCalibrationEntryNptPath e)
  99       return (NptExt npt m d)
 100     where
 101       idx = xrdCalibrationEntryIdx e
 102       (Nxs f _ p) = xrdCalibrationEntryNxs e
 103
 104 -- | Poni Calibration
 105
 106 -- The minimized function is the quadratic difference of the
 107 -- theoretical tth angle and for each pixel, the computed tth angle.
 108
 109 -- synonyme types use in order to improve the calibration performance
 110
 111 type NptEntry' = (Double, [Vector Double]) -- tth, detector pixels coordinates
 112 type Npt' = (Double, [NptEntry']) -- wavelength, [NptEntry']
 113 type NptExt' a = (Npt', Matrix Double, Detector a)
 114
 115 poniEntryFromList :: PoniEntry -> [Double] -> PoniEntry
 116 poniEntryFromList p [rot1, rot2, rot3, poni1, poni2, d] =
 117   p { poniEntryDistance = d *~ meter
 118     , poniEntryPoni1 = poni1 *~ meter
 119     , poniEntryPoni2 = poni2 *~ meter
 120     , poniEntryRot1 = rot1 *~ radian
 121     , poniEntryRot2 = rot2 *~ radian
 122     , poniEntryRot3 = rot3 *~ radian
 123     }
 124 poniEntryFromList _ _ = error "Can not convert to a PoniEntry"
 125
 126 poniEntryToList :: PoniEntry -> [Double]
 127 poniEntryToList p = [ poniEntryRot1 p /~ radian
 128                     , poniEntryRot2 p /~ radian
 129                     , poniEntryRot3 p /~ radian
 130                     , poniEntryPoni1 p /~ meter
 131                     , poniEntryPoni2 p /~ meter
 132                     , poniEntryDistance p /~ meter
 133                     ]
 134
 135 calibrate :: XRDCalibration -> PoniExt -> Detector a -> IO PoniExt
 136 calibrate c (PoniExt p _) d =  do
 137   let entry = last p
 138   let guess = fromList $ poniEntryToList entry
 139   -- read all the NptExt
 140   npts <- mapM (readXRDCalibrationEntry d) (xrdCalibrationEntries c)
 141   -- in order to improve computation speed, pre-compute the pixel coodinates.
 142
 143   let (solution, _p) = minimizeV NMSimplex2 1E-16 3000 box (f (preCalibrate npts)) guess
 144   -- mplot $ drop 3 (toColumns p)
 145   print _p
 146   return $ PoniExt [poniEntryFromList entry (toList solution)] (MyMatrix HklB (ident 3))
 147     where
 148       preCalibrate''' :: Detector a -> NptEntry -> NptEntry'
 149       preCalibrate''' detector (NptEntry _ tth _ points) = (tth /~ radian, map (coordinates detector) points)
 150
 151       preCalibrate'' :: Npt -> Detector a -> Npt'
 152       preCalibrate'' n detector = (nptWavelength n /~ meter, map (preCalibrate''' detector) (nptEntries n))
 153
 154       preCalibrate' :: NptExt a -> NptExt' a
 155       preCalibrate' (NptExt n m detector) = (preCalibrate'' n detector, m', detector)
 156         where
 157           (MyMatrix _ m') = changeBase m PyFAIB
 158
 159       preCalibrate :: [NptExt a] -> [NptExt' a]
 160       preCalibrate = map preCalibrate'
 161
 162       box :: Vector Double
 163       box = fromList [0.1, 0.1, 0.1, 0.01, 0.01, 0.01]
 164
 165       f :: [NptExt' a] -> Vector Double -> Double
 166       f ns params = foldl' (f' rotation translation) 0 ns
 167         where
 168             rot1 = params `atIndex` 0
 169             rot2 = params `atIndex` 1
 170             rot3 = params `atIndex` 2
 171
 172             rotations = map (uncurry fromAxisAndAngle)
 173                         [ (fromList [0, 0, 1], rot3 *~ radian)
 174                         , (fromList [0, 1, 0], rot2 *~ radian)
 175                         , (fromList [1, 0, 0], rot1 *~ radian)]
 176
 177             rotation = foldl' (<>) (ident 3) rotations
 178
 179             translation :: Vector Double
 180             translation = slice 3 3 params
 181
 182       f' :: Matrix Double -> Vector Double -> Double -> NptExt' a -> Double
 183       f' rotation translation x ((_wavelength, entries), m, _detector) =
 184         foldl' (f'' translation r) x entries
 185           where
 186             r :: Matrix Double
 187             r = m <> rotation
 188
 189       f'' :: Vector Double -> Matrix Double -> Double -> NptEntry'-> Double
 190       {-# INLINE f'' #-}
 191       f'' translation r x (tth, pixels) = foldl' (f''' translation r tth) x pixels
 192
 193       f''' :: Vector Double -> Matrix Double -> Double -> Double -> Vector Double -> Double
 194       {-# INLINE f''' #-}
 195       f''' translation r tth x pixel = x + dtth * dtth
 196           where
 197             kf = r #> (pixel - translation)
 198             x' = kf `atIndex` 0
 199             y' = kf `atIndex` 1
 200             z' = kf `atIndex` 2
 201
 202             dtth = tth - atan2 (sqrt (x'*x' + y'*y')) (-z')