chiark / gitweb /
working on rotation - this is not the decomposition I wanted
[moebius2.git] / project.c
1 /*
2  * Displays a conformation
3  */
4
5 #include <X11/Xlib.h>
6 #include <X11/Xutil.h>
7
8 #include "mgraph.h"
9
10 #define MAXTRIS (N*2)
11
12 typedef struct { double vertex[3][D3]; } Triangle;
13
14 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
15 static int ntris;
16 static Vertices conformation;
17
18 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
19 static GSL_MATRIX(transform);
20
21 const char *input_filename;
22
23 static void read_input(void) {
24   FILE *f;
25   int r;
26   
27   f= fopen(input_filename, "rb");  if (!f) diee("input file");
28   errno= 0;
29   r= fread(&conformation,sizeof(conformation),1,f);  if (r!=1) diee("fread");
30   fclose(f);
31 }
32
33 static void transform_coordinates(void) {
34   double result[D3];
35   GSL_VECTOR(result);
36   gsl_vector input_gsl= { D3,1 };
37
38   int v, k;
39   
40   FOR_VERTEX(v) {
41     input_gsl.data= &conformation[v][0];
42     GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
43                        0.0, &result_gsl) );
44     K conformation[v][k]= result[k];
45   }
46 }
47
48 static void addtriangle(int va, int vb, int vc) {
49   Triangle *t= &trisbuffer[ntris];
50   int k;
51   
52   assert(ntris < MAXTRIS);
53   K {
54     t->vertex[0][k]= conformation[va][k];
55     t->vertex[1][k]= conformation[vb][k];
56     t->vertex[2][k]= conformation[vc][k];
57   }
58   displaylist[ntris++]= t;
59 }
60
61 static void generate_display_list(void) {
62   int vb, ve[3], e;
63
64   ntris= 0;
65   FOR_VERTEX(vb) {
66     /* We use the two triangles in the parallelogram vb, vb+e1, vb+e0, vb+e2.
67      * We go round each triangle clockwise (although our surface is non-
68      * orientable so it shouldn't matter).
69      */
70     for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
71     if (ve[0]>=0) {
72       if (ve[1]>=0) addtriangle(vb,ve[0],ve[1]);
73       if (ve[2]>=0) addtriangle(vb,ve[2],ve[0]);
74     }
75   }
76 }    
77
78 static int dl_compare(const void *tav, const void *tbv) {
79   const Triangle *const *tap= tav, *ta= *tap;
80   const Triangle *const *tbp= tbp, *tb= *tbp;
81   double za= ta->vertex[0][2];
82   double zb= tb->vertex[0][2];
83   return za > zb ? -1 :
84          za < zb ? +1 : 0;
85 }
86
87 static void sort_display_list(void) {
88   qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
89 }
90
91 /*---------- X stuff ----------*/
92
93 #define WSZ 400
94
95 static Display *display;
96 static Pixmap pixmap, doublebuffers[2];
97 static Window window;
98 static GC linegc, fillgc;
99 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
100 static int ncut, currentbuffer, x11depth, x11screen;
101
102 static double scale= 0.3;
103 static double eye_z= -10, eye_x= 0;
104 static double cut_z= -9;
105
106 static void drawtriangle(const Triangle *t) {
107   XPoint points[4];
108   int i;
109
110   for (i=0; i<3; i++) {
111     double *v= t->vertex[i];
112     double x= v[0];
113     double y= v[1];
114     double z= v[2];
115
116     if (z < cut_z) { ncut++; return; }
117     
118     double zezezp= eye_z / (eye_z - z);
119     points[i].x= scale * wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
120     points[i].y= scale * wmindim * (zezezp *  y                 ) + wheight/2;
121   }
122   points[3]= points[0];
123
124   XA( XFillPolygon(display,pixmap,fillgc, points,3,Convex,CoordModeOrigin) );
125   XA( XDrawLines(display,pixmap,linegc,points, 4,CoordModeOrigin) );
126 }
127
128 static const unsigned long core_event_mask=
129   ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
130
131 static void mkpixmaps(void) {
132   for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
133     XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
134     doublebuffers[currentbuffer]= pixmap;
135   }
136   currentbuffer= 0;
137 }
138
139 static void display_prepare(void) {
140   XGCValues gcv;
141   XSetWindowAttributes wa;
142   XVisualInfo vinfo;
143   XSizeHints hints;
144   
145   XA( display= XOpenDisplay(0) );
146   x11screen= DefaultScreen(display);
147   x11depth= DefaultDepth(display,x11screen);
148   XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&vinfo) );
149   
150   wa.event_mask= core_event_mask;
151   XA( window= XCreateWindow(display, DefaultRootWindow(display),
152                             0,0, wwidth,wheight, 0,x11depth,
153                             InputOutput, vinfo.visual,
154                             CWEventMask, &wa) );
155
156   hints.flags= USPosition;
157   hints.x= 10;
158   hints.y= 10;
159   XSetWMNormalHints(display,window,&hints);
160
161   mkpixmaps();
162
163   gcv.function= GXcopy;
164   gcv.plane_mask= ~0UL;
165   gcv.foreground= WhitePixel(display,x11screen);
166   linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
167   
168   gcv.function= GXclear;
169   gcv.plane_mask= ~0UL;
170   fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
171 }
172
173 static void display_conformation(void) {
174   Triangle *const *t;
175   int i;
176   
177   pixmap= doublebuffers[currentbuffer];
178   XA( XFillRectangle(display,pixmap,fillgc,0,0,wwidth,wheight) );
179   for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
180     drawtriangle(*t);
181   printf("shown, %d/%d triangles cut\n", ncut, ntris);
182   
183   XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
184   XA( XClearWindow(display,window) );
185   currentbuffer= !currentbuffer;
186 }
187
188 static void show(void) {
189   read_input();
190   transform_coordinates();
191   generate_display_list();
192   sort_display_list();
193   display_conformation();
194 }
195
196 typedef struct {
197   const char *name;
198   void (*start)(void);
199   void (*delta)(double dx, double dy);
200   void (*conclude)(void);
201   void (*abandon)(void);
202 } Drag;
203
204 #define DRAG(x)                                 \
205   static const Drag drag_##x= {                 \
206     #x, drag_##x##_start, drag_##x##_delta,     \
207     drag_##x##_conclude, drag_##x##_abandon     \
208   }
209
210 #define DRAG_SAVING(x, thing)                           \
211   static typeof(thing) original_##thing;                \
212   static void drag_##x##_start(void) {                  \
213     memcpy(&original_##thing, &thing, sizeof(thing));   \
214   }                                                     \
215   static void drag_##x##_conclude(void) { }             \
216   static void drag_##x##_abandon(void) {                \
217     memcpy(&thing, &original_##thing, sizeof(thing));   \
218     show();                                             \
219   }                                                     \
220   DRAG(x)
221
222 static void drag_none_start(void) { }
223 static void drag_none_delta(double dx, double dy) { }
224 static void drag_none_conclude(void) { }
225 static void drag_none_abandon(void) { }
226 DRAG(none);
227
228 static void pmatrix(const char *n, double m[D3][D3]) {
229   int j,k;
230   for (j=0; j<D3; j++) {
231     printf("%10s [ ",n);
232     K printf("%# f ",m[j][k]);
233     printf("]\n");
234     n="";
235   }
236 }
237
238 #define PMATRIX(x) pmatrix(#x,x);
239
240 static void drag_rotate_delta(double dx, double dy) {
241   /* We multiple our transformation matrix by this matrix:
242    *     [  1   0   0 ]
243    *     [  0   1   0 ]
244    *     [ dx  dy   1 ]
245    * and then renormalise.
246    */
247
248   static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
249   static GSL_MATRIX(rotateby);
250
251   double qr[D3][D3];
252   gsl_matrix qr_gsl= { D3,D3,D3,&qr[0][0] };
253
254   double tau[D3];
255   GSL_VECTOR(tau);
256   
257   rotateby[2][0]= dx;
258   rotateby[2][1]= dy;
259
260   PMATRIX(rotateby);
261   
262   GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
263                      &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
264
265   PMATRIX(rotateby);
266
267   GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
268   GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
269                            &transform_gsl, &rotateby_gsl /*dummy*/) );
270
271   PMATRIX(transform);
272   PMATRIX(rotateby);
273
274   printf("drag_rotate_delta...\n");
275   show();
276 }
277 DRAG_SAVING(rotate, transform);
278
279 static void drag_scale_delta(double dx, double dy) {
280   scale *= pow(3.0, dy);
281   show();
282 }
283 DRAG_SAVING(scale, scale);
284
285 static const Drag *drag= &drag_none;
286
287 static int drag_last_x, drag_last_y;
288
289 static void drag_position(int x, int y) {
290   drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
291               (y - drag_last_y) * 1.0 / wmaxdim);
292   drag_last_x= x;
293   drag_last_y= y;
294 }
295
296 static void event_button(XButtonEvent *e) {
297   if (e->window != window || !e->same_screen) return;
298   if (e->type == ButtonPress) {
299     if (e->state || drag != &drag_none) {
300       printf("drag=%s press state=0x%lx abandon\n",
301              drag->name, (unsigned long)e->state);
302       drag->abandon();
303       drag= &drag_none;
304       return;
305     }
306     switch (e->button) {
307     case Button1: drag= &drag_rotate;  break;
308     case Button2: drag= &drag_scale;   break;
309     default: printf("unknown drag start %d\n", e->button);
310     }
311     printf("drag=%s press button=%lu start %d,%d\n",
312            drag->name, (unsigned long)e->button, e->x, e->y);
313     drag_last_x= e->x;
314     drag_last_y= e->y;
315     drag->start();
316   }
317   if (e->type == ButtonRelease) {
318     printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
319     drag_position(e->x, e->y);
320     drag->conclude();
321     drag= &drag_none;
322   }
323 }
324
325 static void event_motion(int x, int y) {
326   printf("drag=%s motion %d,%d\n", drag->name, x, y);
327   drag_position(x,y);
328 }
329
330 static void event_config(XConfigureEvent *e) {
331   if (e->width == wwidth && e->height == wheight)
332     return;
333
334   wwidth= e->width;  wheight= e->height;
335   wmaxdim= wwidth > wheight ? wwidth : wheight;
336   wmindim= wwidth < wheight ? wwidth : wheight;
337   
338   XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
339   for (currentbuffer=0; currentbuffer<2; currentbuffer++)
340     XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
341
342   mkpixmaps();
343   show();
344 }
345
346 int main(int argc, const char *const *argv) {
347   XEvent event;
348   int k;
349   int motion_deferred=0, motion_x=-1, motion_y=-1;
350   
351   if (argc != 2 || argv[1][0]=='-') {
352     fputs("need filename\n",stderr); exit(8);
353   }
354   input_filename= argv[1];
355
356   read_input();
357   K transform[k][k]= 1.0;
358   display_prepare();
359   show();
360
361   XMapWindow(display,window);
362   for (;;) {
363     if (motion_deferred) {
364       int r= XCheckMaskEvent(display,~0UL,&event);
365       if (!r) {
366         event_motion(motion_x, motion_y);
367         motion_deferred=0;
368         continue;
369       }
370     } else {
371       XNextEvent(display,&event);
372     }
373     switch (event.type) {
374
375     case ButtonPress:
376     case ButtonRelease:     event_button(&event.xbutton);     break;
377       
378     case ConfigureNotify:   event_config(&event.xconfigure);  break;
379
380     case MotionNotify:
381 #if 0
382       motion_x= event.xmotion.x;
383       motion_y= event.xmotion.y;
384       motion_deferred= 1;
385 #endif
386       continue;
387       
388     default:
389       printf("unknown event type %u 0x%x\n", event.type,event.type);
390     }
391   }
392 }
393