Let ξ = { x j } j = 1 n be a set of n sample points in the d -cube I d ≔ [ 0 , 1 ] d , and Φ = { φ j } j = 1 n a family of n functions on I d . We define the linear sampling algorithm L n ( Φ , ξ , ⋅ ) for an approximate recovery of a continuous function f on I d from the sampled values f ( x 1 ) , … , f ( x n ) , by L n ( Φ , ξ , f ) ≔ ∑ j = 1 n f ( x j ) φ j . For the Besov class B p , θ α of mixed smoothness α , to study optimality of L n ( Φ , ξ , ⋅ ) in L q ( I d ) we use the quantity r n ( B p , θ α ) q ≔ inf ξ , Φ sup f ∈ B p , θ α ‖ f − L n ( Φ , ξ , f ) ‖ q , where the infimum is taken over all sets of n sample points ξ = { x j } j = 1 n and all families Φ = { φ j } j = 1 n in L q ( I d ) . We explicitly constructed linear sampling algorithms L n ( Φ , ξ , ⋅ ) on the set of sample points ξ = G d ( m ) ≔ { ( 2 − k 1 s 1 , … , 2 − k d s d ) ∈ I d : k 1 + ⋯ + k d ≤ m } , with Φ a family of linear combinations of mixed B-splines which are mixed tensor products of either integer or half integer translated dilations of the centered B-spline of order r . For various 0 < p , q , θ ≤ ∞ and 1 / p < α < r , we proved upper bounds for the worst case error sup f ∈ B p , θ α ‖ f − L n ( Φ , ξ , f ) ‖ q which coincide with the asymptotic order of r n ( B p , θ α ) q in some cases. A key role in constructing these linear sampling algorithms, plays a quasi-interpolant representation of functions f ∈ B p , θ α by mixed B-spline series with the coefficient functionals which are explicitly constructed as linear combinations of an absolute constant number of values of functions. Moreover, we proved that the quasi-norm of the Besov space M B p , θ α is equivalent to a discrete quasi-norm in terms of the coefficient functionals.