Basic as it is, however, in my experience, it is very important to have a detailed appreciation of how the method works in order to really get into the 3D world. Here is a worked example using Octave/Matlab of my simplified, geometry-based implementation of how to calculate torsion angle, including how to determine its sign. No theory or (complicated) mathematical formula, just a step-by-step illustration of how I solve this problem.

- Coordinates of four points in variable abcd:
abcd = [ 21.350 31.325 22.681

22.409 31.286 21.483

22.840 29.751 21.498

23.543 29.175 22.594 ]; - Two auxiliary functions: norm_vec() to normalize a vector; get_orth_norm_vec() to get the orthogonal component (normalized) of a vector with reference to another vector, which should have been normalized.
function ovec = norm_vec(vec)

ovec = vec / norm(vec);

endfunction

function ovec = get_orth_norm_vec(vec, vref)

temp = vec - vref * dot(vec, vref);

ovec = norm_vec(temp);

endfunction - Get three vectors: b_c is the normalized vector b→c; b_a_orth is the orthogonal component (normalized) of vector b→a with reference to b→c; c_d_orth is similarly defined, as the orthogonal component (normalized) of vector c→d with reference to b→c.
b_c = norm_vec(abcd(3, :) - abcd(2, :))

% [0.2703158 -0.9627257 0.0094077]

b_a_orth = get_orth_norm_vec(abcd(1, :) - abcd(2, :), b_c)

% [-0.62126 -0.16696 0.76561]

c_d_orth = get_orth_norm_vec(abcd(4, :) - abcd(3, :), b_c)

% [0.41330 0.12486 0.90199] - Now the torsion angle is the angle between the two vectors, b_a_orth and c_d_orth, and can be easily calculated by their dot product. The sign of the torsion angle is determined by the relative orientation of the cross product of the same two vectors with reference to the middle vector b→c. Here they are in opposite direction, thus the torsion angle is negative.
angle_deg = acos(dot(b_a_orth, c_d_orth)) * 180 / pi

% 65.609

sign = dot(cross(b_a_orth, c_d_orth), b_c)

% -0.91075

if (sign < 0)

ang_deg = -angle_deg % -65.609

endif

Moreover, the method to calculate twist angles of helical nucleic acid structures in SCHNAaP and 3DNA is essentially the same.