Source

papl / src / PaplTransform.ml

Full commit
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
(*
  Copyright (c) 2012 Anders Lau Olsen.
  See LICENSE file for terms and conditions.
*)

(* ---------------------------------------------------------------------- *)

type rng_t = PaplRandom.rng_t

let pi = BatFloat.pi
let pi_x_2 = 2. *. pi

module SO2 = struct
  module V2D = PaplVector.V2D

  type t = float
  type vec = V2D.t

  let rotate a = a

  let step a b = mod_float (b -. a) pi_x_2

  let flip_snd (a, b) =
    let d = a -. b in
    let r = d /. pi_x_2 in
    let n = floor (0.5 +. r) in
    let b = b +. n *. pi_x_2 in
      if abs_float (b -. a) > pi then
        failwith "Internal error: Angle distance greater than pi."
      else b

  let in_range x (a, b) =
    ceil ((a -. x) /. pi_x_2) <= floor ((b -. x) /. pi_x_2)

  let interpolate a b =
    let b = flip_snd (a, b) in
      PaplInterpolate.Float.interpolate a b

  (* Same as
       abs_float (flip_snd (a, b) -. a)
     but slightly faster.
  *)
  let dist a b =
    let d = a -. b in
    let abs_d = abs_float d in
      if abs_d < pi then abs_d (* Specialization for speed. *)
      else
        let r = d /. pi_x_2 in
        let n = floor (0.5 +. r) in
          abs_float (n *. pi_x_2 -. d)

  let id = 0.

  let mul a b = a +. b

  let inverse a = -. a

  let apply a (x, y) =
    let ca = cos a in
    let sa = sin a in
      (ca *. x -. sa *. y, sa *. x +. ca *. y)

  let angle a = a

  let complex a = { Complex.re = cos a; im = sin a }

  module Sampler = struct
    let uniform ?rng range =
      BatEnum.map rotate (PaplSampler.uniform ?rng range)

    let uniform_all ?rng () = uniform ?rng (-. pi, pi)
  end
end

(* ---------------------------------------------------------------------- *)

module SO3 = struct
  module V3D = PaplVector.V3D
  module V4D = PaplVector.V4D

  type t = V4D.t
  type vec = PaplVector.V3D.t

  let is_zero x eps = abs_float x < eps

  let dot = V4D.dot
  let (+:) = V4D.(+:)
  let ( *: ) = V4D.( *: )
  let ( /: ) = V4D.( /: )

  let unit = V4D.unit
  let norm_sqr = V4D.norm2_sqr
  let norm = V4D.norm2
  let neg = V4D.neg
  let map2 = V4D.map2

  let conj (a, b, c, d) = (a, -.b, -.c, -.d)

  let id = 1., 0., 0., 0.

  let ca_sa angle = let a = angle /. 2. in cos a, sin a

  let rotate_x angle = let ca, sa = ca_sa angle in (ca, sa, 0., 0.)
  let rotate_y angle = let ca, sa = ca_sa angle in (ca, 0., sa, 0.)
  let rotate_z angle = let ca, sa = ca_sa angle in (ca, 0., 0., sa)

  let rotate_angle_axis angle ((x, y, z) as axis) =
    let ca, sa = ca_sa angle in
    let len = V3D.norm2 axis in
    let eps = 1e-12 in
      if len < eps then invalid_arg "rotate_axis: axis is zero"
      else
        let s = sa /. len in
          (ca, s *. x, s *. y, s *. z)

  let mul (a0, b0, c0, d0) (a1, b1, c1, d1) =
    a0 *. a1 -. b0 *. b1 -. c0 *. c1 -. d0 *. d1,
    a0 *. b1 +. b0 *. a1 +. c0 *. d1 -. d0 *. c1,
    a0 *. c1 -. b0 *. d1 +. c0 *. a1 +. d0 *. b1,
    a0 *. d1 +. b0 *. c1 -. c0 *. b1 +. d0 *. a1

  let inverse q = conj q /: norm_sqr q

  let clamp x (a, b) eps =
    let die () = invalid_arg
      (Printf.sprintf
         "clamp: x=%.16f is outside interval (%.16f, %.16f) by more than eps=%.16f"
         x a b eps)
    in
      if x < a then
        if a -. x < eps
        then a
        else die ()
      else if b < x then
        if x -. b < eps
        then b
        else die ()
      else
        x

  let get_angle a = 2.0 *. acos (clamp a (0., 1.) 1e-12)

  let angle (a, _, _, _) = get_angle (abs_float a)

  let dist q0 q1 = get_angle (abs_float (dot q0 q1))

  let angle_axis ((a, b, c, d) as q) =
    let len = V3D.norm2 (b, c, d) in
    let eps = 1e-12 in
      if len < eps then
        (* Rotate zero about an arbitrary vector. *)
        (0., (1., 0., 0.))
      else
        let angle = angle q in
        let axis = (b /. len, c /. len, d /. len) in
          (angle, axis)

  let unit_quaternion q = q

  let rotate_unit_quaternion q = unit q

  let apply q (x, y, z) =
    let (_, x', y', z') = mul (mul q (0., x, y, z)) (inverse q)
    in (x', y', z')

  let interpolate q0 q1 =
    let lerp q0 q1 s = unit
      (map2
         (fun a b -> PaplInterpolate.Float.interpolate a b s)
         q0 q1)
    in
    let slerp ca q0 q1 =
      let eps = 1e-7 in
        if is_zero (1.0 -. ca) eps then lerp q0 q1
        else
          let a = acos (clamp ca (0., 1.) 1e-10) in
            fun s ->
              let f1 = sin ((1. -. s) *. a) in
              let f2 = sin (s *. a) in
              let sa = sin a in
              let e1 = (f1 /. sa) *: q0 in
              let e2 = (f2 /. sa) *: q1 in
                e1 +: e2
    in
    let ca = dot q0 q1 in
      if ca < 0. then
        slerp (-. ca) q0 (neg q1)
      else
        slerp ca q0 q1

  module Sampler = struct

    (*
      A random unit quaternion.

      See Section 5.2.2 of Steven LaValle's Planning Algorithms.
      http://planning.cs.uiuc.edu/node198.html
    *)
    let get_uniform_all ?rng =
      let float_fun = PaplRandom.float ?rng in
      let r () = float_fun 1.0 in
        fun () ->
          let open BatFloat in
          let u1, u2, u3 = r (), r (), r () in
          let a2 = pi_x_2 * u2 in
          let a3 = pi_x_2 * u3 in
          let s2 = sqrt (1. - u1) in
          let s3 = sqrt u1 in
            (s2 * sin a2, s2 * cos a2, s3 * sin a3, s3 * cos a3)

    let uniform_all ?rng () =
      BatEnum.from (get_uniform_all ?rng)

    let inv_pi = 1. /. pi

    (* [angle_to_frequency a] is the fraction of random rotations that will have
       a magnitude less than or equal to [a]. *)
    let angle_to_frequency a = inv_pi *. (a -. sin a)

    (* The derivative of the [angle_to_frequency] function with respect to
       [a]. *)
    let angle_to_frequency_deriv a = inv_pi *. (1. -. cos a)

    (* [frequency_to_angle r] finds the angle [a] for which a fraction [r] of
       all random rotations will have a magnitude less than or equal to [a]. *)
    let frequency_to_angle r =
      (* We wish to solve [g r a = 0] for [a]. *)
      let g a = angle_to_frequency a -. r in
      (* The derivative of [g]. *)
      let dg a = angle_to_frequency_deriv a in
      (* A simple but OK approximation to the inverse of [g] that is used to
         compute a seed to Newton's method. *)
      let approx_inv_g r = sqrt (pi *. pi *. r) in
      (* The maximum error allowed. *)
      let eps = 1e-6 in
      (* Iterate Newton's method until the error is below the threshold. *)
      let rec loop a =
        let e = g a in
          if abs_float e <= eps then a
          else loop (a -. e /. dg a)
      in
        (* By using a precomputed table instead of this approximate solution it
           may be possibly to save 1-2 iterations of Newton's method, but I
           don't think those little savings are worth the effort. *)
        loop (approx_inv_g r)

    (* [uniform_angle_within_offset offset] is an angle [a] of rotation in the
       range [0, offset], 0 <= offset <= pi, randomly selected in such a way
       that rotation matrices generated by rotation with this angle will be
       uniformly distributed within the region constrained by the [offset]
       value. *)
    let uniform_angle_within_offset ?rng =
      let small_offset = pi /. 50. in
      let float_fun = PaplRandom.float ?rng in
        fun offset ->
          if offset < small_offset then
            (* Approximate formula for when the surface is approximately a
               disk. See e.g.
               http://mathworld.wolfram.com/DiskPointPicking.html
               or derive for yourself that taking the square root gives the
               suitable scaling for the sampling of a disk.
            *)
            sqrt (float_fun 1.) *. offset
          else
            let r = angle_to_frequency offset in
              frequency_to_angle (float_fun r)

    let get_uniform_offset ?rng =
      let get_angle = uniform_angle_within_offset ?rng in
      let get_axis = PaplVector.V3D.get_uniform_unit ?rng in
        fun offset ->
          let angle = get_angle offset in
          let axis = get_axis () in
            rotate_angle_axis angle axis

    let uniform_offset ?rng =
      let get = get_uniform_offset ?rng in
        fun offset ->
          BatEnum.from (fun () -> get offset)
  end
end

module MakeSE
  (SO : PaplSpatialGroup.S)
  (V : PaplVector.ALL with type t = SO.vec)
  (Misc : sig
     val zero : V.t
     val dist : SO.t PaplMetric.t
   end) =
struct
  type vec = V.t
  type rot = SO.t

  type t = V.t * SO.t

  let (+:) = V.(+:)

  let interpolate =
    PaplInterpolate.Tuple2.interpolate V.interpolate SO.interpolate

  let make t r = (t, r)

  let make_translate t = make t SO.id

  let make_rotate r = make Misc.zero r

  let rotation (t, r) = r

  let translation (t, r) = t

  let get pair = pair

  let id = make Misc.zero SO.id

  let inverse (t, r) =
    let r_inv = SO.inverse r in
      (SO.apply r_inv (V.neg t), r_inv)

  let mul (t1, r1) (t2, r2) = (t1 +: SO.apply r1 t2, SO.mul r1 r2)

  let apply (t, r) pnt =
    let open V in t +: SO.apply r pnt

  let mp wp = PaplMetric.scale wp V.dist2
  let mo wo = PaplMetric.scale wo Misc.dist

  let dist1 wp wo = PaplMetric.Tuple2.dist1 (mp wp) (mo wo)

  let dist2_sqr wp wo = PaplMetric.Tuple2.dist2_sqr (mp wp) (mo wo)

  let dist2 wp wo = PaplMetric.Tuple2.dist2 (mp wp) (mo wo)

  let dist_inf wp wo = (); (* PaplMetric.Tuple2.dist_inf (mp wp) (mo wo) *)
    fun (p0, r0) (p1, r1) ->
      let dp = V.dist2 p0 p1 in
      let dr = Misc.dist r0 r1 in
        Pervasives.max (wp *. dp) (wo *. dr)

  let dist_pos dist = (); fun (x0, _) (x1, _) -> dist x0 x1
  let dist_rot dist = (); fun (_, x0) (_, x1) -> dist x0 x1
end

module type SE_METRIC = sig
  type t
  type vec
  type rot
  val dist1 : float -> float -> t PaplMetric.t
  val dist2_sqr : float -> float -> t PaplMetric.t
  val dist2 : float -> float -> t PaplMetric.t
  val dist_inf : float -> float -> t PaplMetric.t
  val dist_pos : vec PaplMetric.t -> t PaplMetric.t
  val dist_rot : rot PaplMetric.t -> t PaplMetric.t
end

module type SE = sig
  type t
  type vec
  type rot

  include PaplSpatialGroup.S with type t := t and type vec := vec
  include SE_METRIC with type t := t and type vec := vec and type rot := rot

  val make : vec -> rot -> t
  val make_rotate : rot -> t
  val make_translate : vec -> t

  val rotation : t -> rot
  val translation : t -> vec
  val get : t -> vec * rot
end

module SE2 = struct
  module Misc = struct
    let zero = (0., 0.)
    let dist = SO2.dist
  end

  include MakeSE (SO2) (PaplVector.V2D) (Misc)

  let pos_angle (t, r) = (t, SO2.angle r)
  let pos_complex (t, r) = (t, SO2.complex r)

  let make_pos_angle pos angle = make pos (SO2.rotate angle)

  let make_angle angle = make_rotate (SO2.rotate angle)

  module Sampler = struct
    type box_t = vec * vec

    let make sp so = PaplSampler.product2 sp so

    let uniform ?rng box range =
      make
        (PaplVector.V2D.uniform ?rng box)
        (SO2.Sampler.uniform ?rng range)

    let uniform_all ?rng box =
      make
        (PaplVector.V2D.uniform ?rng box)
        (SO2.Sampler.uniform_all ?rng ())
  end
end

module SE3 = struct
  module Misc = struct
    let zero = (0., 0., 0.)
    let dist = SO3.dist
  end

  include MakeSE (SO3) (PaplVector.V3D) (Misc)

  let to_array ((p0, p1, p2), (w, x, y, z)) =
    let (-) = (-.) in
    let ( * ) = ( *. ) in
    let (+) = (+.) in
    let xx = 2. * (x * x) in
    let yy = 2. * (y * y) in
    let zz = 2. * (z * z) in
    let xy = 2. * (x * y) in
    let xz = 2. * (x * z) in
    let yz = 2. * (y * z) in
    let xw = 2. * (x * w) in
    let yw = 2. * (y * w) in
    let zw = 2. * (z * w) in
      [| 1. - yy - zz; xy - zw     ; xz + yw     ; p0;
         xy + zw     ; 1. - xx - zz; yz - xw     ; p1;
         xz - yw     ; yz + xw     ; 1. - xx - yy; p2 |]

  let of_array xs =
    match xs with
        [| r00; r01; r02; px;
           r10; r11; r12; py;
           r20; r21; r22; pz |] ->
          let trace = r00 +. r11 +. r22 in
          let rot =
            if trace > 0. then
              let r = sqrt (trace +. 1.) in
              let s = 0.5 /. r in
                (0.5 *. r,
                 (r21 -. r12) *. s,
                 (r02 -. r20) *. s,
                 (r10 -. r01) *. s)
            else if r00 > r11 && r00 > r22 then
              let r = sqrt (r00 -. r11 -. r22 +. 1.) in
              let s = 0.5 /. r in
                ((r21 -. r12) *. s,
                 0.5 *. r,
                 (r01 +. r10) *. s,
                 (r02 +. r20) *. s)
            else if r11 > r22 then
              let r = sqrt (1. +. r11 -. r00 -. r22) in
              let s = 0.5 /. r in
                ((r02 -. r20) *. s,
                 (r01 +. r10) *. s,
                 0.5 *. r,
                 (r12 +. r21) *. s)
            else
              let r = sqrt (1. +. r22 -. r00 -. r11) in
              let s = 0.5 /. r in
                ((r10 -. r01) *. s,
                 (r02 +. r20) *. s,
                 (r12 +. r21) *. s,
                 0.5 *. r)
          in
            (* Mathematically the scaling to a unit vector isn't needed, but it
               may help guard against small round off errors in the matrix. *)
            make (px, py, pz) (SO3.unit rot)
      | _ -> invalid_arg "PaplTransform.SE3.of_array: 12 elements expected."

  module Sampler = struct
    type box_t = vec * vec

    let make sp so = PaplSampler.product2 sp so

    let uniform_offset ?rng box offset =
      make
        (PaplVector.V3D.uniform ?rng box)
        (SO3.Sampler.uniform_offset ?rng offset)

    let uniform_all ?rng box =
      make
        (PaplVector.V3D.uniform ?rng box)
        (SO3.Sampler.uniform_all ?rng ())
  end
end