xpict_core/
rings.rs

1//! Practical SSSR and ring-interior normals for depiction.
2//!
3//! Layout backends own coordinates. This module only supplies drawing helpers:
4//! which bonds sit in small rings, and unit normals toward the ring centroid
5//! so Kekulé offsets stay short and inside (not the acyclic extend-onto-singles
6//! path). Mirrors `python/xpict/draw/rings.py`.
7
8use std::collections::{HashMap, HashSet, VecDeque};
9
10use crate::scene::{AtomIn, MoleculeIn};
11
12/// One SSSR ring as an ordered atom cycle (atom indices).
13#[derive(Debug, Clone, PartialEq, Eq)]
14pub struct Ring {
15    pub atoms: Vec<i32>,
16}
17
18impl Ring {
19    pub fn size(&self) -> usize {
20        self.atoms.len()
21    }
22
23    pub fn bonds(&self) -> HashSet<(i32, i32)> {
24        let n = self.atoms.len();
25        let mut out = HashSet::with_capacity(n);
26        for i in 0..n {
27            out.insert(bond_key(self.atoms[i], self.atoms[(i + 1) % n]));
28        }
29        out
30    }
31}
32
33fn bond_key(a: i32, b: i32) -> (i32, i32) {
34    if a < b {
35        (a, b)
36    } else {
37        (b, a)
38    }
39}
40
41fn adjacency(mol: &MoleculeIn) -> HashMap<i32, Vec<i32>> {
42    let mut adj: HashMap<i32, Vec<i32>> = mol
43        .atoms
44        .iter()
45        .map(|a| (a.index, Vec::new()))
46        .collect();
47    for b in &mol.bonds {
48        adj.entry(b.begin).or_default().push(b.end);
49        adj.entry(b.end).or_default().push(b.begin);
50    }
51    adj
52}
53
54/// Shortest cycle containing edge `u`–`v` (remove the edge, BFS `u`→`v`).
55fn shortest_cycle_through_bond(
56    adj: &HashMap<i32, Vec<i32>>,
57    u: i32,
58    v: i32,
59    max_size: usize,
60) -> Option<Vec<i32>> {
61    let mut prev: HashMap<i32, Option<i32>> = HashMap::new();
62    prev.insert(u, None);
63    let mut q = VecDeque::new();
64    q.push_back(u);
65    let mut found = false;
66
67    while let Some(cur) = q.pop_front() {
68        if found {
69            break;
70        }
71        let mut depth = 0usize;
72        let mut t = Some(cur);
73        while let Some(node) = t {
74            if node == u {
75                break;
76            }
77            depth += 1;
78            t = prev.get(&node).copied().flatten();
79        }
80        if depth > max_size.saturating_sub(1) {
81            continue;
82        }
83        let Some(neighbors) = adj.get(&cur) else {
84            continue;
85        };
86        for &nxt in neighbors {
87            if (cur == u && nxt == v) || (cur == v && nxt == u) {
88                continue;
89            }
90            if nxt == v && cur != u {
91                prev.insert(v, Some(cur));
92                found = true;
93                break;
94            }
95            if prev.contains_key(&nxt) {
96                continue;
97            }
98            prev.insert(nxt, Some(cur));
99            q.push_back(nxt);
100        }
101    }
102
103    if !found || !prev.contains_key(&v) {
104        return None;
105    }
106    let mut path = vec![v];
107    let mut cur = Some(v);
108    while cur != Some(u) {
109        cur = prev.get(&cur?)?.as_ref().copied();
110        let node = cur?;
111        path.push(node);
112        if path.len() > max_size {
113            return None;
114        }
115    }
116    path.reverse();
117    if path.len() < 3 || path.len() > max_size {
118        return None;
119    }
120    Some(path)
121}
122
123/// Practical SSSR for depiction: unique shortest cycles through each bond.
124pub fn find_sssr(mol: &MoleculeIn, max_size: usize) -> Vec<Ring> {
125    if mol.bonds.is_empty() {
126        return Vec::new();
127    }
128    let adj = adjacency(mol);
129    let mut by_bonds: HashMap<Vec<(i32, i32)>, Ring> = HashMap::new();
130    for b in &mol.bonds {
131        let Some(cyc) = shortest_cycle_through_bond(&adj, b.begin, b.end, max_size) else {
132            continue;
133        };
134        let ring = Ring { atoms: cyc };
135        let mut keys: Vec<_> = ring.bonds().into_iter().collect();
136        keys.sort_unstable();
137        match by_bonds.get(&keys) {
138            Some(prev) if prev.size() <= ring.size() => {}
139            _ => {
140                by_bonds.insert(keys, ring);
141            }
142        }
143    }
144    let mut rings: Vec<Ring> = by_bonds.into_values().collect();
145    rings.sort_by(|a, b| {
146        a.size()
147            .cmp(&b.size())
148            .then_with(|| a.atoms.cmp(&b.atoms))
149    });
150    // Greedy cover: prefer smaller rings; keep a ring if it adds a new bond.
151    let mut kept = Vec::new();
152    let mut covered: HashSet<(i32, i32)> = HashSet::new();
153    for ring in &rings {
154        let bonds = ring.bonds();
155        if bonds.iter().any(|b| !covered.contains(b)) {
156            covered.extend(bonds);
157            kept.push(ring.clone());
158        }
159    }
160    if kept.is_empty() {
161        rings
162    } else {
163        kept
164    }
165}
166
167fn ring_centroid(ring: &Ring, coords: &HashMap<i32, (f64, f64)>) -> Option<(f64, f64)> {
168    if ring.atoms.is_empty() {
169        return None;
170    }
171    let mut sx = 0.0;
172    let mut sy = 0.0;
173    for &i in &ring.atoms {
174        let (x, y) = *coords.get(&i)?;
175        sx += x;
176        sy += y;
177    }
178    let n = ring.atoms.len() as f64;
179    Some((sx / n, sy / n))
180}
181
182/// Unit normals for ring bonds pointing toward the owning (smallest) ring centroid.
183pub fn bond_interior_normals(
184    rings: &[Ring],
185    coords: &HashMap<i32, (f64, f64)>,
186) -> HashMap<(i32, i32), (f64, f64)> {
187    let mut normals = HashMap::new();
188    let mut sorted: Vec<&Ring> = rings.iter().collect();
189    sorted.sort_by_key(|r| r.size());
190    for ring in sorted {
191        let Some((cx, cy)) = ring_centroid(ring, coords) else {
192            continue;
193        };
194        let n = ring.size();
195        for i in 0..n {
196            let a = ring.atoms[i];
197            let b = ring.atoms[(i + 1) % n];
198            let key = bond_key(a, b);
199            if normals.contains_key(&key) {
200                continue;
201            }
202            let Some(&(x1, y1)) = coords.get(&a) else {
203                continue;
204            };
205            let Some(&(x2, y2)) = coords.get(&b) else {
206                continue;
207            };
208            let dx = x2 - x1;
209            let dy = y2 - y1;
210            let length = (dx * dx + dy * dy).sqrt().max(1e-9);
211            let mut nx = -dy / length;
212            let mut ny = dx / length;
213            let mx = 0.5 * (x1 + x2);
214            let my = 0.5 * (y1 + y2);
215            if (cx - mx) * nx + (cy - my) * ny < 0.0 {
216                nx = -nx;
217                ny = -ny;
218            }
219            normals.insert(key, (nx, ny));
220        }
221    }
222    normals
223}
224
225/// Fill [`BondIn::interior`] for ring bonds from SSSR normals.
226///
227/// Ring membership **overrules** the acyclic centered/extend path: any bond
228/// that sits in a small ring gets an interior normal (caller value replaced).
229pub fn apply_ring_interiors(mol: &mut MoleculeIn) {
230    let rings = find_sssr(mol, 8);
231    if rings.is_empty() {
232        return;
233    }
234    let coords: HashMap<i32, (f64, f64)> = mol
235        .atoms
236        .iter()
237        .map(|a: &AtomIn| (a.index, (a.x, a.y)))
238        .collect();
239    let normals = bond_interior_normals(&rings, &coords);
240    for bond in &mut mol.bonds {
241        let key = bond_key(bond.begin, bond.end);
242        if let Some(n) = normals.get(&key) {
243            bond.interior = Some(*n);
244        }
245    }
246}
247
248#[cfg(test)]
249mod tests {
250    use super::*;
251    use crate::scene::BondIn as SceneBond;
252
253    fn benzene() -> MoleculeIn {
254        // Regular hexagon, bond length 20.
255        let mut atoms = Vec::new();
256        for i in 0..6 {
257            let ang = std::f64::consts::PI / 2.0 + i as f64 * std::f64::consts::TAU / 6.0;
258            atoms.push(AtomIn {
259                index: i,
260                element: Some("C".into()),
261                z: None,
262                x: 20.0 * ang.cos(),
263                y: 20.0 * ang.sin(),
264                label: None,
265                charge: 0,
266            });
267        }
268        let orders = [1.0, 2.0, 1.0, 2.0, 1.0, 2.0];
269        let mut bonds = Vec::new();
270        for i in 0..6 {
271            bonds.push(SceneBond {
272                index: i,
273                begin: i,
274                end: (i + 1) % 6,
275                order: orders[i as usize],
276                stereo: None,
277                interior: None,
278            });
279        }
280        MoleculeIn {
281            id: Some("benzene".into()),
282            atoms,
283            bonds,
284            color: None,
285            atom_shade: None,
286            bond_shade: None,
287            shade_vmin: None,
288            shade_vmax: None,
289            mark_atoms: vec![],
290            mark_bonds: vec![],
291            weight: 1.0,
292            scale: 1.0,
293        }
294    }
295
296    #[test]
297    fn benzene_sssr_is_one_hexagon() {
298        let rings = find_sssr(&benzene(), 8);
299        assert_eq!(rings.len(), 1);
300        assert_eq!(rings[0].size(), 6);
301    }
302
303    #[test]
304    fn benzene_doubles_get_interior() {
305        let mut mol = benzene();
306        apply_ring_interiors(&mut mol);
307        let doubles: Vec<_> = mol.bonds.iter().filter(|b| b.order >= 1.5).collect();
308        assert_eq!(doubles.len(), 3);
309        for b in doubles {
310            assert!(b.interior.is_some(), "bond {} missing interior", b.index);
311        }
312    }
313
314    fn ethanol_chain() -> MoleculeIn {
315        MoleculeIn {
316            id: Some("etoh".into()),
317            atoms: vec![
318                AtomIn {
319                    index: 0,
320                    element: Some("C".into()),
321                    z: None,
322                    x: 0.0,
323                    y: 0.0,
324                    label: None,
325                    charge: 0,
326                },
327                AtomIn {
328                    index: 1,
329                    element: Some("C".into()),
330                    z: None,
331                    x: 20.0,
332                    y: 0.0,
333                    label: None,
334                    charge: 0,
335                },
336                AtomIn {
337                    index: 2,
338                    element: Some("O".into()),
339                    z: None,
340                    x: 30.0,
341                    y: 10.0,
342                    label: Some("OH".into()),
343                    charge: 0,
344                },
345            ],
346            bonds: vec![
347                SceneBond {
348                    index: 0,
349                    begin: 0,
350                    end: 1,
351                    order: 1.0,
352                    stereo: None,
353                    interior: None,
354                },
355                SceneBond {
356                    index: 1,
357                    begin: 1,
358                    end: 2,
359                    order: 1.0,
360                    stereo: None,
361                    interior: None,
362                },
363            ],
364            color: None,
365            atom_shade: None,
366            bond_shade: None,
367            shade_vmin: None,
368            shade_vmax: None,
369            mark_atoms: vec![],
370            mark_bonds: vec![],
371            weight: 1.0,
372            scale: 1.0,
373        }
374    }
375
376    #[test]
377    fn acyclic_find_sssr_empty_and_apply_is_noop() {
378        let mut mol = ethanol_chain();
379        assert!(find_sssr(&mol, 8).is_empty());
380        apply_ring_interiors(&mut mol);
381        assert!(mol.bonds.iter().all(|b| b.interior.is_none()));
382
383        let empty = MoleculeIn {
384            id: None,
385            atoms: vec![],
386            bonds: vec![],
387            color: None,
388            atom_shade: None,
389            bond_shade: None,
390            shade_vmin: None,
391            shade_vmax: None,
392            mark_atoms: vec![],
393            mark_bonds: vec![],
394            weight: 1.0,
395            scale: 1.0,
396        };
397        assert!(find_sssr(&empty, 8).is_empty());
398    }
399
400    #[test]
401    fn naphthalene_sssr_has_two_hexagons() {
402        // Two fused hexagons sharing edge 0–1; atoms 6..=9 are reflections of 2..=5.
403        let mut atoms = Vec::new();
404        for i in 0..6 {
405            let ang = std::f64::consts::PI / 2.0 + i as f64 * std::f64::consts::TAU / 6.0;
406            atoms.push(AtomIn {
407                index: i,
408                element: Some("C".into()),
409                z: None,
410                x: 20.0 * ang.cos(),
411                y: 20.0 * ang.sin(),
412                label: None,
413                charge: 0,
414            });
415        }
416        let mx = 0.5 * (atoms[0].x + atoms[1].x);
417        let my = 0.5 * (atoms[0].y + atoms[1].y);
418        let reflected: Vec<(f64, f64)> = (2..6)
419            .map(|i| {
420                let a = &atoms[i as usize];
421                (mx - (a.x - mx), my - (a.y - my))
422            })
423            .collect();
424        for (k, &(x, y)) in reflected.iter().enumerate() {
425            atoms.push(AtomIn {
426                index: 6 + k as i32,
427                element: Some("C".into()),
428                z: None,
429                x,
430                y,
431                label: None,
432                charge: 0,
433            });
434        }
435        let mut bonds = Vec::new();
436        let mut bi = 0i32;
437        for i in 0..6 {
438            bonds.push(SceneBond {
439                index: bi,
440                begin: i,
441                end: (i + 1) % 6,
442                order: 1.0,
443                stereo: None,
444                interior: None,
445            });
446            bi += 1;
447        }
448        for &(a, b) in &[(1, 6), (6, 7), (7, 8), (8, 9), (9, 0)] {
449            bonds.push(SceneBond {
450                index: bi,
451                begin: a,
452                end: b,
453                order: 1.0,
454                stereo: None,
455                interior: None,
456            });
457            bi += 1;
458        }
459        let mol = MoleculeIn {
460            id: Some("naph".into()),
461            atoms,
462            bonds,
463            color: None,
464            atom_shade: None,
465            bond_shade: None,
466            shade_vmin: None,
467            shade_vmax: None,
468            mark_atoms: vec![],
469            mark_bonds: vec![],
470            weight: 1.0,
471            scale: 1.0,
472        };
473        let rings = find_sssr(&mol, 8);
474        assert_eq!(rings.len(), 2, "expected two hexagons, got {:?}", rings);
475        assert!(rings.iter().all(|r| r.size() == 6));
476    }
477}
xpict_core/rings.rs

xpict_core/
rings.rs
