1 // OGL_RIBBON.CPP
2
3 // Copyright (C) 1998 Tommi Hassinen, Jarno Huuskonen.
4
5 // This package is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9
10 // This package is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with this package; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 /*################################################################################################*/
20
21 #include "ogl_ribbon.h" // config.h is here -> we get ENABLE-macros here...
22
23 #include <sstream>
24 using namespace std;
25
26 #include <ghemical/v3d.h>
27 #include <ghemical/eng1_sf.h>
28
29 #include "spline.h"
30
31 /*################################################################################################*/
32
33 const i32s ogl_ribbon::resol = 10;
34
35 const fGL ogl_ribbon::width = 0.15;
36 const fGL ogl_ribbon::helix = 0.10;
37
ogl_ribbon(project * p1,color_mode * p9,i32s p2,i32s order)38 ogl_ribbon::ogl_ribbon(project * p1, color_mode * p9, i32s p2, i32s order) :
39 ogl_smart_object()
40 {
41 prj = p1; cmode = p9;
42 chn = p2; extra_points = order - 2;
43
44 vector<chn_info> & ci_vector = (* prj->GetCI());
45
46 if (ci_vector[chn].GetType() != chn_info::amino_acid)
47 {
48 assertion_failed(__FILE__, __LINE__, "chn type != AA");
49 }
50
51 if (ci_vector[chn].GetLength() < 3)
52 {
53 assertion_failed(__FILE__, __LINE__, "chn length < 3");
54 }
55
56 const char * tmp_sss = ci_vector[chn].GetSecStrStates();
57
58 char * state = new char[ci_vector[chn].GetLength() + 1];
59 state[ci_vector[chn].GetLength()] = 0;
60
61 // the default state is always loop. strands have the same places as in K&S.
62 // helices are shifted: the smallest helical peptide is "4...." -> "LHHHL".
63
64 setup * su = prj->GetCurrentSetup();
65 setup1_sf * susf = dynamic_cast<setup1_sf *>(su);
66 if (susf != NULL) // optionally, if we use SF, get the constraints...
67 {
68 bool error = false;
69 if (susf->chn_vector.size() != ci_vector.size()) error = true;
70 if (susf->chn_vector[chn].res_vector.size() != (i32u) ci_vector[chn].GetLength()) error = true;
71
72 if (error)
73 {
74 assertion_failed(__FILE__, __LINE__, "susf error.");
75 }
76
77 for (i32u n1 = 0;n1 < susf->chn_vector[chn].res_vector.size();n1++)
78 {
79 state[n1] = 'L';
80 if (susf->chn_vector[chn].res_vector[n1].GetState() == SF_STATE_STRAND) state[n1] = 'S';
81 }
82
83 for (i32s n1 = 0;n1 < ((i32s) susf->chn_vector[chn].res_vector.size()) - 3;n1++)
84 {
85 if (susf->chn_vector[chn].res_vector[n1].GetState() == SF_STATE_HELIX4)
86 {
87 state[n1 + 1] = 'H';
88 state[n1 + 2] = 'H';
89 state[n1 + 3] = 'H';
90 }
91 }
92 }
93 else // normally, use the DSSP results...
94 {
95 for (i32s n1 = 0;n1 < ci_vector[chn].GetLength();n1++)
96 {
97 state[n1] = 'L';
98 if (tmp_sss[n1] == 'S') state[n1] = 'S';
99 }
100
101 for (i32s n1 = 0;n1 < ((i32s) ci_vector[chn].GetLength()) - 3;n1++)
102 {
103 if (tmp_sss[n1] == '4')
104 {
105 state[n1 + 1] = 'H';
106 state[n1 + 2] = 'H';
107 state[n1 + 3] = 'H';
108 }
109 }
110 }
111
112 iter_al range1[2];
113 prj->GetRange(1, chn, range1);
114
115 length = ci_vector[chn].GetLength() - 1;
116
117 cp1 = new fGL[length * 3];
118 cp2 = new fGL[length * 3];
119
120 data3 = new fGL[length * 3];
121
122 v3d<fGL> old_tmpv7; bool flag = true;
123 for (i32s n1 = 0;n1 < length;n1++)
124 {
125 fGL color1[4]; color1[3] = 1.0;
126 fGL color2[4]; color2[3] = 1.0;
127
128 if (susf != NULL) // optional constraint colours...
129 {
130 color1[0] = 0.5; color1[1] = 0.0; color1[2] = 0.8;
131 if (state[n1 + 0] == 'H') { color1[0] = 0.8; color1[1] = 0.5; color1[2] = 0.0; }
132 if (state[n1 + 0] == 'S') { color1[0] = 0.0; color1[1] = 0.8; color1[2] = 0.5; }
133
134 color2[0] = 0.5; color2[1] = 0.0; color2[2] = 0.8;
135 if (state[n1 + 1] == 'H') { color2[0] = 0.8; color2[1] = 0.5; color2[2] = 0.0; }
136 if (state[n1 + 1] == 'S') { color2[0] = 0.0; color2[1] = 0.8; color2[2] = 0.5; }
137 }
138 else // the normal colours...
139 {
140 color1[0] = 0.0; color1[1] = 0.0; color1[2] = 1.0;
141 if (state[n1 + 0] == 'H') { color1[0] = 1.0; color1[1] = 0.0; color1[2] = 0.0; }
142 if (state[n1 + 0] == 'S') { color1[0] = 0.0; color1[1] = 1.0; color1[2] = 0.0; }
143
144 color2[0] = 0.0; color2[1] = 0.0; color2[2] = 1.0;
145 if (state[n1 + 1] == 'H') { color2[0] = 1.0; color2[1] = 0.0; color2[2] = 0.0; }
146 if (state[n1 + 1] == 'S') { color2[0] = 0.0; color2[1] = 1.0; color2[2] = 0.0; }
147 }
148
149 iter_al c_alpha_1;
150 iter_al c_alpha_2;
151
152 v3d<fGL> ca_to_ca;
153 v3d<fGL> normal(1.0, 0.0, 0.0);
154
155 iter_al range2a[2]; prj->GetRange(2, range1, n1 + 0, range2a);
156 iter_al range2b[2]; prj->GetRange(2, range1, n1 + 1, range2b);
157
158 c_alpha_2 = range2b[0]; while (c_alpha_2 != range2b[1] && ((* c_alpha_2).builder_res_id & 0xFF) != 0x01) c_alpha_2++;
159 if (c_alpha_2 == range2b[1])
160 {
161 assertion_failed(__FILE__, __LINE__, "C_alpha_2 not found.");
162 }
163
164 c_alpha_1 = range2a[0]; while (c_alpha_1 != range2a[1] && ((* c_alpha_1).builder_res_id & 0xFF) != 0x01) c_alpha_1++;
165 if (c_alpha_1 == range2a[1])
166 {
167 assertion_failed(__FILE__, __LINE__, "C_alpha_1 not found.");
168 }
169
170 iter_al c_carboxyl = range2a[0]; while (c_carboxyl != range2a[1] && ((* c_carboxyl).builder_res_id & 0xFF) != 0x02) c_carboxyl++;
171 if (c_carboxyl == range2a[1])
172 {
173 assertion_failed(__FILE__, __LINE__, "C_carboxyl not found.");
174 }
175
176 ca_to_ca = v3d<fGL>((* c_alpha_1).GetCRD(0), (* c_alpha_2).GetCRD(0));
177 v3d<fGL> tmpv1 = v3d<fGL>((* c_alpha_1).GetCRD(0), (* c_carboxyl).GetCRD(0));
178 normal = ca_to_ca.vpr(tmpv1); normal = normal / normal.len();
179
180 v3d<fGL> tmpv7 = normal.vpr(ca_to_ca); tmpv7 = tmpv7 / tmpv7.len();
181
182 v3d<fGL> midpoint = v3d<fGL>((* c_alpha_1).GetCRD(0));
183 midpoint = midpoint + (ca_to_ca / 2.0);
184
185 i32s helix_count = 0;
186 if (state[n1 + 0] == 'H') helix_count++;
187 if (state[n1 + 1] == 'H') helix_count++;
188 midpoint = midpoint + (normal * ((fGL) helix_count * helix));
189
190 v3d<fGL> tmpv8 = midpoint + (tmpv7 * width);
191 v3d<fGL> tmpv9 = midpoint - (tmpv7 * width);
192
193 if (n1 != 0 && tmpv7.ang(old_tmpv7) > M_PI / 2.0) flag = !flag;
194 old_tmpv7 = tmpv7;
195
196 for (i32s n2 = 0;n2 < 3;n2++)
197 {
198 if (flag)
199 {
200 cp1[n1 * 3 + n2] = tmpv8[n2];
201 cp2[n1 * 3 + n2] = tmpv9[n2];
202 }
203 else
204 {
205 cp1[n1 * 3 + n2] = tmpv9[n2];
206 cp2[n1 * 3 + n2] = tmpv8[n2];
207 }
208
209 data3[n1 * 3 + n2] = (color1[n2] + color2[n2]) / 2.0;
210 }
211 }
212
213 ref1 = new spline(order, length + extra_points * 2);
214
215 head_points1 = new fGL[extra_points * 3];
216 tail_points1 = new fGL[extra_points * 3];
217
218 head_refs1[0] = & cp1[1 * 3]; head_refs1[1] = & cp1[0 * 3];
219 tail_refs1[0] = & cp1[(length - 2) * 3]; tail_refs1[1] = & cp1[(length - 1) * 3];
220
221 ref2 = new spline(order, length + extra_points * 2);
222
223 head_points2 = new fGL[extra_points * 3];
224 tail_points2 = new fGL[extra_points * 3];
225
226 head_refs2[0] = & cp2[1 * 3]; head_refs2[1] = & cp2[0 * 3];
227 tail_refs2[0] = & cp2[(length - 2) * 3]; tail_refs2[1] = & cp2[(length - 1) * 3];
228
229 UpdateExtraPoints(head_refs1, head_points1, tail_refs1, tail_points1);
230 UpdateExtraPoints(head_refs2, head_points2, tail_refs2, tail_points2);
231
232 for (i32s n1 = 0;n1 < extra_points;n1++)
233 {
234 ref1->SetPoint(n1, & head_points1[(extra_points - (n1 + 1)) * 3]);
235 ref2->SetPoint(n1, & head_points2[(extra_points - (n1 + 1)) * 3]);
236 }
237
238 for (i32s n1 = 0;n1 < extra_points;n1++)
239 {
240 ref1->SetPoint(extra_points + length + n1, & tail_points1[n1 * 3]);
241 ref2->SetPoint(extra_points + length + n1, & tail_points2[n1 * 3]);
242 }
243
244 for (i32s n1 = 0;n1 < length;n1++)
245 {
246 ref1->SetPoint(extra_points + n1, & cp1[n1 * 3]);
247 ref2->SetPoint(extra_points + n1, & cp2[n1 * 3]);
248 }
249
250 for (i32s n1 = 0;n1 < length + extra_points * 2 + order;n1++)
251 {
252 ref1->SetKnot(n1, (fGL) (n1 - extra_points) - (fGL) order / 2.0);
253 ref2->SetKnot(n1, (fGL) (n1 - extra_points) - (fGL) order / 2.0);
254 }
255
256 #ifdef RIBBON_USE_DISPLSTS
257 list_id = prj->GetDisplayListIDs(1);
258 #endif // RIBBON_USE_DISPLSTS
259
260 i32s np1 = resol * (length - 1);
261 i32s np2 = np1 + 1;
262
263 data1 = new fGL[np2];
264 data2a = new fGL[np2 * 3];
265 data2b = new fGL[np2 * 3];
266
267 for (i32s n1 = 0;n1 < np2;n1++)
268 {
269 data1[n1] = ((fGL) (length - 1) * n1) / (fGL) np1;
270 i32s index = (i32s) floor(data1[n1] + 0.5);
271
272 fGL width1 = 0.5; if (state[index + 0] == 'H' || state[index + 0] == 'S') width1 = 1.0;
273 fGL width2 = 0.5; if (state[index + 1] == 'H' || state[index + 1] == 'S') width2 = 1.0;
274
275 fGL ang = M_PI * ((data1[n1] + 0.5) - (fGL) index) / 2.0;
276 fGL tmp1 = cos(ang); fGL w1 = width1 * tmp1 * tmp1;
277 fGL tmp2 = sin(ang); fGL w2 = width2 * tmp2 * tmp2;
278 fGL ribbon_width = w1 + w2;
279
280 fGL pos1[3]; ref1->Compute(data1[n1], pos1);
281 fGL pos2[3]; ref2->Compute(data1[n1], pos2);
282
283 v3d<fGL> v1 = v3d<fGL>(pos1);
284 v3d<fGL> v2 = v3d<fGL>(pos1, pos2); v2 = v2 / 2.0;
285
286 v3d<fGL> v3 = (v1 + v2) - (v2 * ribbon_width);
287 v3d<fGL> v4 = (v1 + v2) + (v2 * ribbon_width);
288
289 for (i32s n2 = 0;n2 < 3;n2++)
290 {
291 data2a[n1 * 3 + n2] = v3[n2];
292 data2b[n1 * 3 + n2] = v4[n2];
293 }
294 }
295
296 delete[] state;
297 }
298
~ogl_ribbon(void)299 ogl_ribbon::~ogl_ribbon(void)
300 {
301 delete[] cp1;
302 delete[] cp2;
303
304 delete ref1;
305 delete[] head_points1;
306 delete[] tail_points1;
307
308 delete ref2;
309 delete[] head_points2;
310 delete[] tail_points2;
311
312 delete[] data1;
313 delete[] data2a;
314 delete[] data2b;
315 delete[] data3;
316
317 #ifdef RIBBON_USE_DISPLSTS
318 prj->DeleteDisplayLists(list_id, 1);
319 #endif // RIBBON_USE_DISPLSTS
320 }
321
UpdateExtraPoints(fGL ** hr,fGL * hp,fGL ** tr,fGL * tp)322 void ogl_ribbon::UpdateExtraPoints(fGL ** hr, fGL * hp, fGL ** tr, fGL * tp)
323 {
324 v3d<fGL> tmpv1; v3d<fGL> tmpv2;
325
326 tmpv1 = v3d<fGL>(hr[1]);
327 tmpv2 = v3d<fGL>(hr[0], hr[1]);
328 for (i32s n1 = 0;n1 < extra_points;n1++)
329 {
330 v3d<fGL> tmpv3 = tmpv1 + (tmpv2 * ((fGL) n1 + 1));
331 for (i32s n2 = 0;n2 < 3;n2++) hp[n1 * 3 + n2] = tmpv3[n2];
332 }
333
334 tmpv1 = v3d<fGL>(tr[1]);
335 tmpv2 = v3d<fGL>(tr[0], tr[1]);
336 for (i32s n1 = 0;n1 < extra_points;n1++)
337 {
338 v3d<fGL> tmpv3 = tmpv1 + (tmpv2 * ((fGL) n1 + 1));
339 for (i32s n2 = 0;n2 < 3;n2++) tp[n1 * 3 + n2] = tmpv3[n2];
340 }
341 }
342
Render(void)343 void ogl_ribbon::Render(void)
344 {
345 #ifdef RIBBON_USE_DISPLSTS
346 if (glIsList(list_id) == GL_TRUE) glCallList(list_id);
347 else
348 {
349 glNewList(list_id, GL_COMPILE_AND_EXECUTE);
350 #endif // RIBBON_USE_DISPLSTS
351
352 i32s np1 = resol * (length - 1);
353
354 glEnable(GL_LIGHTING);
355 glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, true);
356
357 glBegin(GL_TRIANGLES);
358 for (i32s n1 = 0;n1 < np1;n1++)
359 {
360 v3d<fGL> va = v3d<fGL>(& data2a[(n1 + 0) * 3], & data2b[(n1 + 0) * 3]);
361 v3d<fGL> vb = v3d<fGL>(& data2a[(n1 + 1) * 3], & data2b[(n1 + 1) * 3]);
362
363 v3d<fGL> vc = v3d<fGL>(& data2a[(n1 + 0) * 3], & data2a[(n1 + 1) * 3]);
364 v3d<fGL> vd = v3d<fGL>(& data2b[(n1 + 0) * 3], & data2b[(n1 + 1) * 3]);
365
366 v3d<fGL> vA1 = va.vpr(vc); v3d<fGL> vA2 = vb.vpr(vc);
367 v3d<fGL> vB1 = va.vpr(vd); v3d<fGL> vB2 = vb.vpr(vd);
368
369 if (n1 != 0)
370 {
371 v3d<fGL> ve = v3d<fGL>(& data2a[(n1 + 0) * 3], & data2a[(n1 - 1) * 3]);
372 v3d<fGL> vf = v3d<fGL>(& data2b[(n1 + 0) * 3], & data2b[(n1 - 1) * 3]);
373 vA1 = vA1 + (ve.vpr(va)); vB1 = vB1 + (vf.vpr(va));
374 }
375
376 if (n1 != (np1 - 1))
377 {
378 v3d<fGL> vg = v3d<fGL>(& data2a[(n1 + 1) * 3], & data2a[(n1 + 2) * 3]);
379 v3d<fGL> vh = v3d<fGL>(& data2b[(n1 + 1) * 3], & data2b[(n1 + 2) * 3]);
380 vA2 = vA2 + (vb.vpr(vg)); vB2 = vB2 + (vb.vpr(vh));
381 }
382
383 vA1 = vA1 / vA1.len(); vA2 = vA2 / vA2.len();
384 vB1 = vB1 / vB1.len(); vB2 = vB2 / vB2.len();
385
386 i32s index = (i32s) floor(data1[n1]);
387 fGL color1[3]; fGL mod1 = data1[n1 + 0] - (fGL) index;
388 fGL color2[3]; fGL mod2 = data1[n1 + 1] - (fGL) index;
389
390 for (i32s n2 = 0;n2 < 3;n2++)
391 {
392 color1[n2] = data3[index * 3 + n2] * (1.0 - mod1) + data3[(index + 1) * 3 + n2] * mod1;
393 color2[n2] = data3[index * 3 + n2] * (1.0 - mod2) + data3[(index + 1) * 3 + n2] * mod2;
394 }
395
396 glColor3fv(color2);
397 glNormal3fv(vA2.data); glVertex3fv(& data2a[(n1 + 1) * 3]);
398
399 glColor3fv(color1);
400 glNormal3fv(vA1.data); glVertex3fv(& data2a[(n1 + 0) * 3]);
401 glNormal3fv(vB1.data); glVertex3fv(& data2b[(n1 + 0) * 3]);
402
403 glColor3fv(color1);
404 glNormal3fv(vB1.data); glVertex3fv(& data2b[(n1 + 0) * 3]);
405
406 glColor3fv(color2);
407 glNormal3fv(vB2.data); glVertex3fv(& data2b[(n1 + 1) * 3]);
408 glNormal3fv(vA2.data); glVertex3fv(& data2a[(n1 + 1) * 3]);
409 }
410
411 glEnd(); // GL_TRIANGLES
412
413 glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, false);
414 glDisable(GL_LIGHTING);
415
416 #ifdef RIBBON_USE_DISPLSTS
417 glEndList();
418 }
419 #endif // RIBBON_USE_DISPLSTS
420 }
421
422 /*################################################################################################*/
423
424 // eof
425