rdex-client colour histogram calculation and upload
[rdex:client.git] / src / library.c
1 /* =====================================================================
2 rdex -- reaction-diffusion explorer
3 Copyright (C) 2008  Claude Heiland-Allen <claudiusmaximus@goto10.org>
4 ------------------------------------------------------------------------
5 Library Module
6 ===================================================================== */
7
8 #include <assert.h>
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <string.h>
12 #include <wordexp.h>
13 #include <math.h>
14 #include <gsl/gsl_math.h>
15 #include <gsl/gsl_eigen.h>
16 #include "library.h"
17 #include "upload.h"
18
19 //======================================================================
20 // prototypes
21
22 struct library *library_init(struct library *library) {
23   if (!library) { return 0; }
24   if (! borderwindow_init(&library->borderwindow)) { return 0; }
25   library->count = 0;
26   library->arraysize = 256;
27   library->array = malloc(sizeof(struct libelem) * library->arraysize);
28   library->distmatrix = malloc(sizeof(struct libdist) * library->arraysize * library->arraysize);
29   for (int i = 0; i < 4; ++i) { library->vmin[i] =  1.0e6; }
30   for (int i = 0; i < 4; ++i) { library->vmax[i] = -1.0e6; }
31   for (int i = 0; i < 4; ++i) { library->vsub[i] = 0; }
32   for (int i = 0; i < 4; ++i) { library->vmul[i] = 0; }
33   const GLfloat m0[16] = {1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
34   for (int i = 0; i < 16; ++i) { library->m[i] = m0[i]; }
35   library->d = 4;
36   library->X = 0;
37   library->Y = 0;
38   library->Z = 0;
39   library->width = 1;
40   library->height = 1;
41   library->snap = 0;
42   library->image = image_alloc(128, 128, 3); //FIXME
43   tamura_init(&library->tamura, 128, 128); //FIXME
44   upload_init(&library->upload); //FIXME
45   library->buffersize = library_tex_size * library_tex_size * 3;
46   library->buffer = malloc(library->buffersize);
47   library->session = getenv("RDEX_SESSION");
48   if (! library->session) { library->session = "."; }
49   library_load(library);
50   return library;
51 }
52
53 void library_grow(struct library *library) {
54   int newarraysize = library->arraysize * 2;
55   struct libelem *newarray = malloc(sizeof(struct libelem) * newarraysize);
56   memcpy(newarray, library->array, sizeof(struct libelem) * library->count);
57   free(library->array);
58   struct libdist *newmatrix = malloc(sizeof(struct libdist) * newarraysize * newarraysize);
59   for (int i = 0; i < library->arraysize; ++i) {
60     memcpy(&newmatrix[newarraysize * i], &library->distmatrix[library->arraysize * i], sizeof(struct libdist) * library->count);
61   }
62   free(library->distmatrix);
63   library->arraysize = newarraysize;
64   library->array = newarray;
65   library->distmatrix = newmatrix;
66 }
67
68 void library_append(
69   struct library *library,
70   GLfloat ru, GLfloat rv, GLfloat f, GLfloat k, GLuint t,
71   struct tamura_features *feature
72 ) {
73   assert(library->count < library->arraysize);
74   library->vmin[0] = fmin(library->vmin[0], ru);
75   library->vmin[1] = fmin(library->vmin[1], rv);
76   library->vmin[2] = fmin(library->vmin[2], f);
77   library->vmin[3] = fmin(library->vmin[3], k);
78   library->vmax[0] = fmax(library->vmax[0], ru);
79   library->vmax[1] = fmax(library->vmax[1], rv);
80   library->vmax[2] = fmax(library->vmax[2], f);
81   library->vmax[3] = fmax(library->vmax[3], k);
82   library->array[library->count].v[0] = ru;
83   library->array[library->count].v[1] = rv;
84   library->array[library->count].v[2] = f;
85   library->array[library->count].v[3] = k;
86   library->array[library->count].t = t;
87   library->array[library->count].feature[0] = feature[0].coarseness;
88   library->array[library->count].feature[1] = feature[0].contrast;
89   library->array[library->count].feature[2] = feature[0].directionality;
90   library->array[library->count].feature[3] = feature[1].coarseness;
91   library->array[library->count].feature[4] = feature[1].contrast;
92   library->array[library->count].feature[5] = feature[1].directionality;
93   for (int i = 0; i < library->count; ++i) {
94     float dt = 0;
95     for (int j = 0; j < 6; ++j) {
96       float fd = library->array[i].feature[j] - library->array[library->count].feature[j];
97       dt += fd * fd;
98     }
99     dt = sqrtf(dt);
100     library->distmatrix[i * library->arraysize + library->count].texture = dt;
101     library->distmatrix[i + library->arraysize * library->count].texture = dt;
102     library->distmatrix[library->count * library->arraysize + library->count].texture = 0.0;
103   }
104   library->count += 1;
105   if (library->count == library->arraysize) {
106     library_grow(library);
107   }
108 }
109
110 void library_rerange(struct library *library) {
111   for (int i = 0; i < 4; ++i) {
112     library->vmul[i] = library->vmax[i] - library->vmin[i];
113     library->vsub[i] = library->vmin[i] + library->vmul[i]/2;
114     if (library->vmul[i] > 0) { library->vmul[i] = 1.0/library->vmul[i]; }
115   }
116 /*
117   if (library->count > 0) {
118     double sigma = 0;
119     double oldsigma = 0;
120     double dsigma = 0;
121     double epsilon = 0.0001;
122     gsl_matrix *X[2];
123     X[0] = gls_matrix_alloc(library->count, library->count);
124     X[1] = gls_matrix_alloc(library->count, library->count);
125     int t = 0;
126     gsl_matrix *A = matrix_alloc(library->count, library->count);
127     gsl_matrix *B = matrix_alloc(library->count, library->count);
128     do {
129       oldsigma = sigma;
130
131       // Aij_ii = 1
132       // Aij_jj = 1
133       // Aij_ij = -1
134       // Aij_ji = -1
135       // Aij_?? = 0
136       // sij(Y)_xy = delta_ij / d_ij(Y_xy)   or 0 when /0
137       // B(Y ) = sum_ij (wij sij (Y )Aij
138       // X[1-t] = n^(−1) B(X[t]) X[t]
139       for (int i = 0; i < n; ++i) {
140         for (int j = 0; j < n; ++j) {
141           // euclidean output distance
142           double dX_ij = 0;
143           for (int s = 0; s < p; ++s) {
144             double d = gsl_matrix_get(X[t], i, s) - gsl_matrix_get(X[t], j, s);
145             dX_ij += d * d;
146           }
147            gsl_matrix_set(B,i,j, gsl_matrix_get(B,i,j)
148              + delta_ij_ii / d_ij_ii
149              + delta_ij_jj / d_ij_jj
150              - delta_ij_ij / d_ij_ij
151              - delta_ij_ji / d_ij_ji);
152         }
153       }
154       // calculate stress
155       sigma = 0;
156       for (int i = 0; i < n; ++i) {
157         for (int j = 0; j < n; ++j) {
158           // euclidean output distance
159           double dX_ij = 0;
160           for (int s = 0; s < p; ++s) {
161             double d = gsl_matrix_get(X[t], i, s) - gsl_matrix_get(X[t], j, s);
162             dX_ij += d * d;
163           }
164           dX_ij = sqrt(dX_ij);
165           // difference from input dissimilarity
166           double d = gsl_matrix_get(delta, i, j) - dX_ij;
167           sigma += d * d;
168         }
169       }
170
171       // check convergence
172       dsigma = sigma - oldsigma;
173     } while (dsigma > epsilon);
174
175 */
176 /*
177   // multidimensional scaling of dissimilarities { d_rs }
178   // Find matrix A = [ -0.5 d_rs^2 ]
179   gsl_matrix *A = gsl_matrix_alloc(library->count, library->count);
180   for (int i = 0; i < library->count; ++i) {
181     for (int j = 0; j < library->count; ++j) {
182       double dij = library->distmatrix[i * library->arraysize + j].texture;
183       gsl_matrix_set(A,i,j, -0.5 * dij * dij);
184     }
185   }
186   // Find matrix B = [a_rs - a_r. - a_.s + a_..]
187   // where
188   //  a_r. = n^-1  sum(s)    a_rs
189   //  a_.s = n^-1  sum(r)    a_rs
190   //  a_.. = n^-2  sum(r)(s) a_rs
191   gsl_vector *Ai_ = gsl_vector_alloc(library->count);
192   gsl_vector *A_j = gsl_vector_alloc(library->count);
193   double A__ = 0.0;
194   for (int i = 0; i < library->count; ++i) {
195     gsl_vector_set(Ai_,i,0.0);
196     gsl_vector_set(A_j,i,0.0);
197   }
198   for (int i = 0; i < library->count; ++i) {
199     for (int j = 0; j < library->count; ++j) {
200       double aij = gsl_matrix_get(A,i,j);
201       gsl_vector_set(Ai_,j, gsl_vector_get(Ai_,j) + aij);
202       gsl_vector_set(A_j,i, gsl_vector_get(A_j,i) + aij);
203       A__ += aij;
204     }
205   }
206   for (int i = 0; i < library->count; ++i) {
207     gsl_vector_set(Ai_,i, gsl_vector_get(Ai_,i) / library->count);
208     gsl_vector_set(A_j,i, gsl_vector_get(A_j,i) / library->count);
209   }
210   A__ /= library->count * library->count;
211   gsl_matrix *B = gsl_matrix_alloc(library->count, library->count);
212   for (int i = 0; i < library->count; ++i) {
213     for (int j = 0; j < library->count; ++j) {
214       gsl_matrix_set(B,i,j, gsl_matrix_get(A,i,j) - gsl_vector_get(Ai_, i) - gsl_vector_get(A_j, j) + A__);
215     }
216   }
217   // Find eigenvalues: l_1 >= l_2 >= ... >= l_n-1  
218   // and eigenvectors: v_1    v_2    ...    v_n-1
219   gsl_vector *e = gsl_vector_alloc(library->count);
220   gsl_matrix *E = gsl_matrix_alloc(library->count, library->count);
221   gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(library->count);
222   gsl_eigen_symmv(B, e, E, w);
223   gsl_eigen_symmv_free(w);
224   gsl_eigen_symmv_sort(e, E, GSL_EIGEN_SORT_VAL_DESC);
225
226   for (int i = 0; i < library->count; ++i) {
227     fprintf(stderr, " %f ", gsl_vector_get(e, i));
228   }
229   fprintf(stderr, "\n");
230
231   // where eigenvectors normalized so that v_i^T v_i = l_i
232   gsl_vector *l = gsl_vector_alloc(library->count);
233 */
234 /*
235   for (int i = 0; i < library->count; ++i) {
236     double s = 0.0;
237     for (int j = 0; j < library->count; ++j) {
238       double t = gsl_matrix_get(E,i,j);
239 //      fprintf(stderr, " %f ", t);
240       s += t * t;
241     }
242 //    fprintf(stderr, "\n");
243     gsl_vector_set(l, i, sqrt(s));
244   }
245 //  fprintf(stderr, "\n");
246 */
247 /*
248   for (int i = 0; i < library->count; ++i) {
249     double k = sqrt(sqrt(gsl_vector_get(e, i))); // / gsl_vector_get(l, j);
250     for (int j = 0; j < library->count; ++j) {
251       gsl_matrix_set(E,i,j, gsl_matrix_get(E,i,j) * k);
252     }
253   }
254   // Coordinates of n points in p dimensional Euclidean space are:
255   //  x_ri = v_ir  (r = 1...n, i = 1...p)
256   for (int i = 0; i < library->count; ++i) {
257     for (int j = 0; j < 4; ++j) {
258       if (j < library->count) {
259         library->array[i].mds[j] = gsl_matrix_get(E, j, i);
260       } else {
261         library->array[i].mds[j] = 0.0;
262       }
263       fprintf(stderr, " %f ", library->array[i].mds[j]);
264     }
265 //    library->array[i].mds[3] = 0.0;
266     fprintf(stderr, "\n");
267   }
268   gsl_matrix_free(A);
269   gsl_vector_free(Ai_);
270   gsl_vector_free(A_j);
271   gsl_matrix_free(B);
272   gsl_vector_free(e);
273   gsl_matrix_free(E);
274   gsl_vector_free(l);
275 */
276 //  }
277 }
278
279 int library_load(struct library *library) {
280   glEnable(GL_TEXTURE_2D);
281   wordexp_t p;
282   const char *pattern2 = "/ru_*__rv_*__f_*__k_*.ppm";
283   int patternlen = strlen(library->session) + strlen(pattern2) + 1;
284   char *pattern = malloc(patternlen);
285   snprintf(pattern, patternlen, "%s%s", library->session, pattern2);
286   if (0 == wordexp(pattern, &p, WRDE_NOCMD | WRDE_SHOWERR | WRDE_UNDEF)) {
287     char sru[9]; char srv[9]; char sf[9]; char sk[9];
288     for (int i = 0; i < p.we_wordc; ++i) {
289       int l = strlen(p.we_wordv[i]);
290       const char *pattern3 = "ru_%8c__rv_%8c__f_%8c__k_%8c.ppm";
291       if (
292         l >= 52 &&
293         4 == sscanf(p.we_wordv[i] + l - 52, pattern3, sru, srv, sf, sk)
294       ) {
295         sru[8] = '\0'; srv[8] = '\0'; sf[8] = '\0'; sk[8] = '\0';
296         double ru, rv, f, k;
297         if (1 != sscanf(sru, "%lg", &ru)) { break; }
298         if (1 != sscanf(srv, "%lg", &rv)) { break; }
299         if (1 != sscanf(sf,  "%lg", &f))  { break; }
300         if (1 != sscanf(sk,  "%lg", &k))  { break; }
301         FILE *image = fopen(p.we_wordv[i], "rb");
302         if (!image) { break; }
303         fseek(image, -library->buffersize, SEEK_END);
304         fread(library->buffer, library->buffersize, 1, image);
305         fclose(image);
306         GLuint t;
307         glGenTextures(1, &t);
308         glBindTexture(GL_TEXTURE_2D, t);
309         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
310         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
311         glTexImage2D(
312           GL_TEXTURE_2D, 0, GL_RGB,
313           library_tex_size, library_tex_size, 0,
314           GL_RGB, GL_UNSIGNED_BYTE, library->buffer
315         );
316         struct tamura_features feature[2];
317         { // FIXME this is a bit dodgy, assumes filename stuff
318           int sn = strlen(p.we_wordv[i]);
319           char *filename = malloc(sn + 4);
320           strncpy(filename, p.we_wordv[i], sn + 1);
321           strncpy(&filename[sn-3], "rdex", 5);
322           FILE *featfile = fopen(filename, "rb");
323           if (featfile) {
324             int count;
325             count = fscanf(
326               featfile,
327               "RDEX/0\n%f %f %f %f %f %f\n",
328               &feature[0].coarseness,
329               &feature[0].contrast,
330               &feature[0].directionality,
331               &feature[1].coarseness,
332               &feature[1].contrast,
333               &feature[1].directionality
334             );
335             fclose(featfile);
336             if (count != 6) { break; }
337           }
338         }
339         library_append(library, ru, rv, f, k, t, &feature[0]);
340       }
341     }
342     wordfree(&p);
343   }
344   free(pattern);
345   library_rerange(library);
346   glDisable(GL_TEXTURE_2D);
347   return 1;
348 }
349
350 void library_display_save(
351   struct library *library, GLuint fbo, GLuint texture, GLuint rawtexture,
352   double ru, double rv, double f, double k, const char *behaviour
353 ) {
354   if (! library->snap) { return; }
355   library->snap = 0;
356
357   struct tamura_features feature[2];
358   {
359     char filename[1024];
360     snprintf(
361       filename, 1000,
362       "%s/ru_%f__rv_%f__f_%f__k_%f.rdex",  // FIXME ensure same stem
363       library->session, ru, rv, f, k
364     );
365     image_copytexture(library->image, rawtexture);
366     tamura_calculate(&library->tamura, &feature[0], library->image);
367     FILE *featurefile = fopen(filename, "wb");
368     if (featurefile) {
369       fprintf(featurefile,
370         "RDEX/0\n"                         // magic number/format version
371         "%f %f %f %f %f %f\n",             // texture vector
372      // "%f %f %f %f %f %f %f %f %f\n",    // colour vector FIXME
373         feature[0].coarseness,
374         feature[0].contrast,
375         feature[0].directionality,
376         feature[1].coarseness,
377         feature[1].contrast,
378         feature[1].directionality
379       );
380       fclose(featurefile);
381     }
382   }
383
384   GLuint t;
385   glGenTextures(1, &t);
386   glEnable(GL_TEXTURE_2D);
387   glBindTexture(GL_TEXTURE_2D, t);
388   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
389   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
390   glTexImage2D(
391     GL_TEXTURE_2D, 0, GL_RGBA,
392     library_tex_size, library_tex_size, 0,
393     GL_RGBA, GL_UNSIGNED_BYTE, 0
394   );
395   glMatrixMode(GL_PROJECTION);
396   glLoadIdentity();
397   glViewport(0, 0, library_tex_size, library_tex_size);
398   gluOrtho2D(0, 1, 0, 1);
399   glMatrixMode(GL_MODELVIEW);
400   glLoadIdentity();
401   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo);
402   glFramebufferTexture2DEXT(
403     GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D,
404     t, 0
405   );
406   glBindTexture(GL_TEXTURE_2D, texture);
407   glBegin(GL_QUADS); { glColor4f(1,1,1,1);
408     glTexCoord2f(0, 0); glVertex2f(0, 0);
409     glTexCoord2f(1, 0); glVertex2f(1, 0);
410     glTexCoord2f(1, 1); glVertex2f(1, 1);
411     glTexCoord2f(0, 1); glVertex2f(0, 1);
412   } glEnd();
413   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
414   char filename[1024];
415   snprintf(
416     filename, 1000,
417     "%s/ru_%f__rv_%f__f_%f__k_%f.ppm",  // FIXME ensure same stem
418     library->session, ru, rv, f, k
419   );
420   FILE *imgfile = fopen(filename, "wb");
421   char ppm[UPLOAD_DATA_PPMLEN];
422   int ppm_bytes = 0;
423   if (imgfile) {
424     glBindTexture(GL_TEXTURE_2D, t);
425     glGetTexImage(
426       GL_TEXTURE_2D, 0, GL_RGB, GL_UNSIGNED_BYTE, library->buffer
427     );
428     ppm_bytes = snprintf(ppm, UPLOAD_DATA_PPMLEN, "P6\n%d %d 255\n", library_tex_size, library_tex_size);
429     for (int y = library_tex_size - 1; y >= 0; --y) {
430       memcpy(ppm + ppm_bytes, library->buffer + y * library_tex_size * 3, library_tex_size * 3);
431       ppm_bytes += library_tex_size * 3;
432     }
433     fwrite(ppm, ppm_bytes, 1, imgfile);
434     fclose(imgfile);
435   }
436   struct histogram histogram;
437   struct image *image = image_alloc(library_tex_size, library_tex_size, 3);
438   image_copytexture(image, t);
439   histogram_calc(&histogram, image);
440   image_free(image);
441   struct upload_data *udata = upload_prepare(
442     &library->upload,
443     ru, rv, f, k, behaviour,
444     feature[0].coarseness, feature[0].contrast, feature[0].directionality,
445     feature[1].coarseness, feature[1].contrast, feature[1].directionality,
446     &histogram, ppm, ppm_bytes
447   );
448   upload_upload(&library->upload, udata);
449   upload_cleanup(udata);
450   library_append(library, ru, rv, f, k, t, &feature[0]);
451   library_rerange(library);
452   glDisable(GL_TEXTURE_2D);
453 }
454
455 void library_reshape(struct library *library, int w, int h) {
456   library->width = w;
457   library->height = h;
458 }
459
460 void library_project(struct library *library) {
461   for (int i = 0; i < library->count; ++i) {
462     GLfloat x0 = (library->array[i].v[0] - library->vsub[0]) * library->vmul[0];
463     GLfloat x1 = (library->array[i].v[1] - library->vsub[1]) * library->vmul[1];
464     GLfloat x2 = (library->array[i].v[2] - library->vsub[2]) * library->vmul[2];
465     GLfloat x3 = (library->array[i].v[3] - library->vsub[3]) * library->vmul[3];
466 /*
467     GLfloat x0 = library->array[i].mds[0];
468     GLfloat x1 = library->array[i].mds[1];
469     GLfloat x2 = library->array[i].mds[2];
470     GLfloat x3 = library->array[i].mds[3];
471 */
472     GLfloat y0 = library->m[ 0] * x0 + library->m[ 4]*x1 + library->m[ 8]*x2 + library->m[12]*x3;
473     GLfloat y1 = library->m[ 1] * x0 + library->m[ 5]*x1 + library->m[ 9]*x2 + library->m[13]*x3;
474     GLfloat y2 = library->m[ 2] * x0 + library->m[ 6]*x1 + library->m[10]*x2 + library->m[14]*x3;
475     GLfloat y3 = library->m[ 3] * x0 + library->m[ 7]*x1 + library->m[11]*x2 + library->m[15]*x3;
476     GLfloat k = 3 * library->d / (library->d + y3);
477     y0 *= k;
478     y1 *= k;
479     y2 *= k;
480     library->array[i].w[0] = y0;
481     library->array[i].w[1] = y1;
482     library->array[i].w[2] = y2;
483   }
484 }
485 /*
486 int library_cmpz(struct libelem *a, struct libelem *b) {
487   if (a->w[2] < b->w[2]) return -1;
488   if (a->w[2] > b->w[2]) return  1;
489   return 0;
490 }
491
492 void library_sort(struct library *library) {
493   qsort(library->array, library->count, sizeof(struct libelem), library_cmpz);
494 }
495 */
496 void library_display(struct library *library, GLuint texture) {
497   glClearColor(0,0,0,1);
498   glClear(GL_COLOR_BUFFER_BIT);
499   int width, height;
500   glEnable(GL_TEXTURE_2D);
501   glBindTexture(GL_TEXTURE_2D, texture);
502   glMatrixMode(GL_PROJECTION);
503   glLoadIdentity();
504   gluOrtho2D(0, 1, 0, 1);
505   glMatrixMode(GL_MODELVIEW);
506   glLoadIdentity();
507   glViewport(0, 0, library->width, library->height);
508   glGetTexLevelParameteriv(
509     GL_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &width
510   );
511   glGetTexLevelParameteriv(
512     GL_TEXTURE_2D, 0, GL_TEXTURE_HEIGHT, &height
513   );
514   glBegin(GL_QUADS); {
515     float x = 0.5 * library->width/width;
516     float y = 0.5 * library->height/height;
517     glColor4f(1,1,1,1);
518     glTexCoord2f(0.5 - x, 0.5 - y); glVertex2f(0, 0);
519     glTexCoord2f(0.5 + x, 0.5 - y); glVertex2f(1, 0);
520     glTexCoord2f(0.5 + x, 0.5 + y); glVertex2f(1, 1);
521     glTexCoord2f(0.5 - x, 0.5 + y); glVertex2f(0, 1);
522   } glEnd();
523   glMatrixMode(GL_PROJECTION);
524   glLoadIdentity();
525   float v = (float) library->width / (float) library->height;
526   glFrustum(-v, v, -1, 1, 4, 12);
527   glMatrixMode(GL_MODELVIEW);
528   glViewport(0, 0, library->width, library->height);
529   glLoadIdentity();
530   // { 4D rolling ball
531     float R = 15;
532     float r = sqrt(
533       library->X * library->X +
534       library->Y * library->Y +
535       library->Z * library->Z
536     );
537     float D = sqrt(R*R+r*r);
538     float c = R / D;
539     float s = r / D;
540     float x, y, z;
541     if (r > 0.00001) {
542       x = library->X / r;
543       y = library->Y / r;
544       z = library->Z / r;
545     } else {
546       x = 0;
547       y = 0;
548       z = 0;
549     }
550     GLfloat m[16] = {
551       1-x*x*(1-c), -(1-c)*x*y,  -(1-c)*x*z,  s*x,
552       -(1-c)*x*y,  1-y*y*(1-c), -(1-c)*y*z,  s*y,
553       -(1-c)*x*z,  -(1-c)*y*z,  1-z*z*(1-c), s*z,
554       -s*x,        -s*y,        -s*z,        c
555     };
556     glMultMatrixf(m);
557     glMultMatrixf(library->m);
558     glGetFloatv(GL_MODELVIEW_MATRIX , library->m);
559   // } rolling ball
560   glLoadIdentity();
561   glTranslatef(0,0,-8);
562   library_project(library);
563   //library_sort(library);
564   glEnable(GL_DEPTH_TEST);
565   glClear(GL_DEPTH_BUFFER_BIT);
566   borderwindow_use(&library->borderwindow);
567   for (int i = 0; i < library->count; ++i) {
568     glBindTexture(GL_TEXTURE_2D, library->array[i].t);
569     GLfloat y0 = library->array[i].w[0];
570     GLfloat y1 = library->array[i].w[1];
571     GLfloat y2 = library->array[i].w[2];
572     glBegin(GL_QUADS); {
573       glColor4f(1, 1, 1, 1);
574       glTexCoord2f(0, 0); glVertex3f(y0 - 0.1, y1 - 0.1, y2);
575       glTexCoord2f(1, 0); glVertex3f(y0 + 0.1, y1 - 0.1, y2);
576       glTexCoord2f(1, 1); glVertex3f(y0 + 0.1, y1 + 0.1, y2);
577       glTexCoord2f(0, 1); glVertex3f(y0 - 0.1, y1 + 0.1, y2);
578     } glEnd();
579   }
580   borderwindow_use(0);
581   glDisable(GL_DEPTH_TEST);
582   glDisable(GL_TEXTURE_2D);
583 }
584
585 void library_pmotion(struct library *library, int x, int y) {
586   library->X =  (x - library->width/2.0)  / (float) library->width;
587   library->Y = -(y - library->height/2.0) / (float) library->height;
588   library->Z = 0;
589 }
590
591 void library_amotion(struct library *library, int x, int y) {
592   library->X = (x - library->width/2.0)  / (float) library->width;
593   library->Y = 0;
594   library->Z = (y - library->height/2.0) / (float) library->height;
595 }
596
597 void library_idle(struct library *library) {
598 }
599
600 void library_dump(struct library *library) {
601   fprintf(stderr, "\nLIBRARY DISTANCE MATRIX\n");
602   for (int i = 0; i < library->count; ++i) {
603     fprintf(stderr, "{");
604     for (int j = 0; j < library->count; ++j) {
605       fprintf(stderr, " %f", library->distmatrix[i * library->arraysize + j].texture);
606     }
607     fprintf(stderr, " }\n");
608   }
609   fprintf(stderr, "-----------------------\n");
610 }
611
612 // EOF