pub struct OpenGL {
pub points: Vec<f32>,
geom: bool,
pub size: usize,
current: Option<(Vector, Vector, Vector)>,
pub sphere_n: usize,
pub backbone_n: usize,
pub pair_n: usize,
}
impl OpenGL {
pub fn new(
use_geometry_shader: bool,
sphere_n: usize,
backbone_n: usize,
pair_n: usize,
) -> Self {
OpenGL {
size: if use_geometry_shader { 5 } else { 8 },
geom: use_geometry_shader,
points: Vec::new(),
current: None,
sphere_n,
backbone_n,
pair_n,
}
}
pub fn reset<A, B>(&mut self, design: &crate::Design<A, B>) {
self.points.clear();
self.current = None;
for i in 0..design.strands.len() {
self.add_strand(design, self.geom, i);
}
let mut nucl = HashMap::new();
for i in 0..design.strands.len() {
self.add_pairs(&design, self.geom, &mut nucl, i);
}
}
fn push_p64(&mut self, x: [f64; 3]) {
self.points.push(x[0] as f32);
self.points.push(x[1] as f32);
self.points.push(x[2] as f32);
}
fn push_int(&mut self, x: u32) {
self.points.push(unsafe { std::mem::transmute(x) });
}
fn add_cylinder(&mut self, i: usize, point: [f64; 3], color: u32, r: f64, n: usize) {
assert!(!self.geom);
let p = Vector([point[0], point[1], point[2]]);
if let Some((p0, vx0, vy0)) = self.current {
let vec = (p - p0).normalize();
let vx = vec.normal();
let vy = vec.vec_product(vx);
let _vx0 = if vx0.is_null() { vx } else { vx0 };
let _vy0 = if vy0.is_null() { vy } else { vy0 };
let mut t = Vector::NULL;
let mut normal = Vector::NULL;
let mut is_first = true;
for rot in 0..=2 * n + 1 {
let angle = -((rot / 2) as f64) * 2. * PI / (n as f64);
normal = r * vx * angle.cos() + r * vy * angle.sin();
t = if rot % 2 == 0 {
p0 + normal
} else {
p + normal
};
if is_first {
self.push_p64(t.0);
self.push_p64(normal.normalize().0);
self.push_int(i as u32);
self.push_int(color);
is_first = false;
}
self.push_p64(t.0);
self.push_p64(normal.normalize().0);
self.push_int(i as u32);
self.push_int(color);
}
self.push_p64(t.0);
self.push_p64(normal.normalize().0);
self.push_int(i as u32);
self.push_int(color);
self.current = Some((p, vx, vy));
} else {
self.current = Some((p, Vector::NULL, Vector::NULL));
}
}
fn add_sphere(&mut self, k: usize, point: [f64; 3], color: u32, r: f64, n: usize) {
assert!(!self.geom);
let p = Vector([point[0], point[1], point[2]]);
let mut last = p;
let mut normal = Vector::NULL;
for i in 0..n {
let alpha = (i as f64) * PI / (n as f64);
let alpha_ = ((i + 1) as f64) * PI / (n as f64);
let r0 = r * alpha.sin();
let r1 = r * alpha_.sin();
for j in 0..=2 * n + 1 {
let beta = -2. * ((j / 2) as f64) * PI / (n as f64);
normal = if j & 1 == 0 {
Vector([r0 * beta.sin(), -r * alpha.cos(), r0 * beta.cos()])
} else {
Vector([r1 * beta.sin(), -r * alpha_.cos(), r1 * beta.cos()])
};
last = p + normal;
if j == 0 {
self.push_p64(last.0);
self.push_p64(normal.normalize().0);
self.push_int(k as u32);
self.push_int(color);
}
self.push_p64(last.0);
self.push_p64(normal.normalize().0);
self.push_int(k as u32);
self.push_int(color)
}
self.push_p64(last.0);
self.push_p64(normal.normalize().0);
self.push_int(k as u32);
self.push_int(color)
}
}
fn add_strand<A, B>(&mut self, design: &crate::Design<A, B>, geom: bool, i: usize) {
let ref s = design.strands[i];
self.current = None;
let color = if let Some(ref color) = s.color {
color.as_int()
} else {
0x888888
};
for d in s.domains.iter() {
if d.helix < 0 {
continue;
}
let h = d.helix as usize;
if d.forward {
for x in d.start..d.end {
let point = design.helices[h].space_pos(
design
.parameters
.as_ref()
.unwrap_or(&crate::Parameters::DEFAULT),
x,
true,
);
if geom {
self.push_p64(point);
self.push_int(i as u32);
self.push_int(color);
} else {
self.add_cylinder(i, point, color, 0.1, self.backbone_n);
self.add_sphere(i, point, color, 0.2, self.sphere_n);
}
}
} else {
for x in (d.start..d.end).rev() {
let point = design.helices[h].space_pos(
design
.parameters
.as_ref()
.unwrap_or(&crate::Parameters::DEFAULT),
x,
false,
);
if geom {
self.push_p64(point);
self.push_int(i as u32);
self.push_int(color);
} else {
self.add_cylinder(i, point, color, 0.1, self.backbone_n);
self.add_sphere(i, point, color, 0.2, self.sphere_n);
}
}
}
}
if geom {
self.push_p64([1. / 0., 0., 0.]);
self.push_int(i as u32);
self.push_int(color);
}
}
fn add_pairs<A, B>(
&mut self,
design: &crate::Design<A, B>,
geom: bool,
nucl: &mut HashMap<(isize, isize, bool), usize>,
i: usize,
) {
let ref s = design.strands[i];
for dom in s.domains.iter() {
if dom.helix < 0 {
continue;
}
for x in dom.start..dom.end {
nucl.insert((dom.helix, x, dom.forward), i);
if let Some(&j) = nucl.get(&(dom.helix, x, !dom.forward)) {
let a = design.helices[dom.helix as usize].space_pos(
design
.parameters
.as_ref()
.unwrap_or(&crate::Parameters::DEFAULT),
x,
dom.forward,
);
let b = design.helices[dom.helix as usize].space_pos(
design
.parameters
.as_ref()
.unwrap_or(&crate::Parameters::DEFAULT),
x,
!dom.forward,
);
let mut a = Vector([a[0], a[1], a[2]]);
let mut b = Vector([b[0], b[1], b[2]]);
let v = (b - a).normalize();
a = a + v * 0.19;
b = b - v * 0.19;
if geom {
let k = (i << 16) | j;
self.push_p64(a.0);
self.push_int(k as u32);
self.push_int(0xaaaaaa);
self.push_p64(b.0);
self.push_int(k as u32);
self.push_int(0xaaaaaa);
self.push_p64([1. / 0., 0., 0.]);
self.push_int(k as u32);
self.push_int(0);
} else {
self.current = Some((a, Vector::NULL, Vector::NULL));
self.add_cylinder((i << 16) | j, b.0, 0xaaaaaa, 0.05, self.pair_n);
}
}
}
}
}
}
use std::f64::consts::PI;
use std::collections::HashMap;
#[derive(Debug, Clone, Copy, PartialEq)]
struct Vector(pub [f64; 3]);
impl Vector {
fn normal(self) -> Self {
let mut arg = 0;
if self.0[1].abs() < self.0[0].abs() {
arg = 1;
}
if self.0[2].abs() < self.0[arg] {
arg = 2;
}
(Vector(if arg == 0 {
[0., self.0[2], -self.0[1]]
} else if arg == 1 {
[self.0[2], 0., -self.0[0]]
} else {
[self.0[1], -self.0[0], 0.]
}))
.normalize()
}
fn normalize(&self) -> Self {
let norm = (self.0[0] * self.0[0] + self.0[1] * self.0[1] + self.0[2] * self.0[2]).sqrt();
Vector([self.0[0] / norm, self.0[1] / norm, self.0[2] / norm])
}
fn vec_product(&self, b: Self) -> Self {
Vector([
self.0[1] * b.0[2] - self.0[2] * b.0[1],
self.0[2] * b.0[0] - self.0[0] * b.0[2],
self.0[0] * b.0[1] - self.0[1] * b.0[0],
])
}
const NULL: Vector = Vector([0., 0., 0.]);
fn is_null(&self) -> bool {
self.0 == Self::NULL.0
}
}
impl std::ops::Sub<Vector> for Vector {
type Output = Vector;
fn sub(self, b: Vector) -> Self::Output {
Vector([self.0[0] - b.0[0], self.0[1] - b.0[1], self.0[2] - b.0[2]])
}
}
impl std::ops::Add<Vector> for Vector {
type Output = Vector;
fn add(self, b: Vector) -> Self::Output {
Vector([self.0[0] + b.0[0], self.0[1] + b.0[1], self.0[2] + b.0[2]])
}
}
impl std::ops::Mul<f64> for Vector {
type Output = Vector;
fn mul(self, b: f64) -> Self::Output {
Vector([self.0[0] * b, self.0[1] * b, self.0[2] * b])
}
}
impl std::ops::Mul<Vector> for f64 {
type Output = Vector;
fn mul(self, b: Vector) -> Self::Output {
Vector([b.0[0] * self, b.0[1] * self, b.0[2] * self])
}
}