REMOVED: the DrawModel function from classifierInterface (they all do exactly the...
[mldemos:baraks-mldemos.git] / Core / drawUtils.cpp
1 /*********************************************************************\r
2 MLDemos: A User-Friendly visualization toolkit for machine learning\r
3 Copyright (C) 2010  Basilio Noris\r
4 Contact: mldemos@b4silio.com\r
5 \r
6 This library is free software; you can redistribute it and/or\r
7 modify it under the terms of the GNU Lesser General Public\r
8 License as published by the Free Software Foundation; either\r
9 version 2.1 of the License, or (at your option) any later version.\r
10 \r
11 This library is distributed in the hope that it will be useful,\r
12 but WITHOUT ANY WARRANTY; without even the implied warranty of\r
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\r
14 Library General Public License for more details.\r
15 \r
16 You should have received a copy of the GNU Lesser General Public\r
17 License along with this library; if not, write to the Free\r
18 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.\r
19 *********************************************************************/\r
20 #include "drawUtils.h"\r
21 #include "glUtils.h"\r
22 #include "basicMath.h"\r
23 #include <QBitmap>\r
24 #include <QDebug>\r
25 #include "qcontour.h"\r
26 #include "kmeans.h"\r
27 #include <jacgrid/jacgrid.h>\r
28 using namespace std;\r
29 \r
30 #define FOUR(a) {a;a;a;a;}\r
31 #define RES 256\r
32 void DrawEllipse(float *mean, float *sigma, float rad, QPainter *painter, QSize size)\r
33 {\r
34     if(mean[0] != mean[0] || mean[1] != mean[1]) return; // nan\r
35     float a = sigma[0], b = sigma[1], c = sigma[2];\r
36     float L[4];\r
37     L[0] = a; L[1] = 0; L[2] = b; L[3] = sqrtf(c*a-b*b);\r
38     if(L[3] != L[3]) L[3] = 0;\r
39     FOR(i,4) L[i] /= sqrtf(a);\r
40 \r
41     const int segments = 64;\r
42     float oldX = FLT_MAX, oldY = FLT_MAX;\r
43     for (float theta=0; theta <= PIf*2.f; theta += (PIf*2.f)/segments)\r
44     {\r
45         float x = cosf(theta)*rad;\r
46         float y = sinf(theta)*rad;\r
47         float nx = L[0]*x;\r
48         float ny = L[2]*x + L[3]*y;\r
49         nx += mean[0];\r
50         ny += mean[1];\r
51         if(oldX != FLT_MAX)\r
52         {\r
53             painter->drawLine(\r
54                         QPointF(nx*size.width(),ny*size.height()),\r
55                         QPointF(oldX*size.width(),oldY*size.height())\r
56                         );\r
57         }\r
58         oldX = nx;\r
59         oldY = ny;\r
60     }\r
61 }\r
62 \r
63 void DrawEllipse(float *mean, float *sigma, float rad, QPainter *painter, Canvas *canvas)\r
64 {\r
65     if(mean[0] != mean[0] || mean[1] != mean[1]) return; // nan\r
66     float a = sigma[0], b = sigma[1], c = sigma[2];\r
67     float L[4];\r
68     L[0] = a; L[1] = 0; L[2] = b; L[3] = sqrtf(c*a-b*b);\r
69     if(L[3] != L[3]) L[3] = 0;\r
70     FOR(i,4) L[i] /= sqrtf(a);\r
71 \r
72     const int segments = 64;\r
73     float oldX = FLT_MAX, oldY = FLT_MAX;\r
74     for (float theta=0; theta <= PIf*2.f; theta += (PIf*2.f)/segments)\r
75     {\r
76         float x = cosf(theta)*rad;\r
77         float y = sinf(theta)*rad;\r
78         float nx = L[0]*x;\r
79         float ny = L[2]*x + L[3]*y;\r
80         nx += mean[0];\r
81         ny += mean[1];\r
82         if(oldX != FLT_MAX)\r
83         {\r
84             painter->drawLine(canvas->toCanvasCoords(nx,ny), canvas->toCanvasCoords(oldX, oldY));\r
85         }\r
86         oldX = nx;\r
87         oldY = ny;\r
88     }\r
89 }\r
90 \r
91 void DrawArrow( const QPointF &ppt, const QPointF &pt, double sze, QPainter &painter)\r
92 {\r
93     QPointF pd, pa, pb;\r
94     double tangent;\r
95 \r
96     pd = ppt - pt;\r
97     if (pd.x() == 0 && pd.y() == 0)\r
98         return;\r
99     tangent = atan2 ((double) pd.y(), (double) pd.x());\r
100     pa.setX(sze * cos (tangent + PIf / 7.f) + pt.x());\r
101     pa.setY(sze * sin (tangent + PIf / 7.f) + pt.y());\r
102     pb.setX(sze * cos (tangent - PIf / 7.f) + pt.x());\r
103     pb.setY(sze * sin (tangent - PIf / 7.f) + pt.y());\r
104     //-- connect the dots...\r
105     painter.drawLine(pt, ppt);\r
106     painter.drawLine(pt, pa);\r
107     painter.drawLine(pt, pb);\r
108 }\r
109 \r
110 QColor ColorFromVector(fvec a)\r
111 {\r
112     // angle is between 0 and 1;\r
113     float angle = atan2(a[0], a[1]) / (2*PIf) + 0.5f;\r
114     vector<fvec> colors;\r
115 #define Col2Col(r,g,b) {fvec c;c.resize(3); c[0] = r;c[1] = g;c[2] = b; colors.push_back(c);}\r
116 \r
117     Col2Col(0,0,255);\r
118     Col2Col(255,0,255);\r
119     Col2Col(255,0,0);\r
120     Col2Col(255,255,0);\r
121     Col2Col(0,255,0);\r
122     Col2Col(0,255,255);\r
123 \r
124     // find where the angle fits in the color list\r
125     int index = (int)(angle*(colors.size())) % colors.size();\r
126     fvec c1 = colors[index];\r
127     fvec c2 = colors[(index+1)%colors.size()];\r
128 \r
129     // compute the ratio between c1 and c2\r
130     float remainder = angle*(colors.size()) - (int)(angle*(colors.size()));\r
131     fvec c3 = c1*(1-remainder) + c2*remainder;\r
132     return QColor(c3[0],c3[1],c3[2]);\r
133 }\r
134 \r
135 QPixmap RocImage(std::vector< std::vector<f32pair> > rocdata, std::vector<const char *> roclabels, QSize size)\r
136 {\r
137     QPixmap pixmap(size);\r
138     pixmap.fill(Qt::white);\r
139     QPainter painter(&pixmap);\r
140     painter.setRenderHint(QPainter::Antialiasing);\r
141 \r
142     int w = pixmap.width(), h = pixmap.height();\r
143     int PAD = 16;\r
144 \r
145     QFont font = painter.font();\r
146     font.setPointSize(12);\r
147     font.setBold(true);\r
148     painter.setFont(font);\r
149 \r
150     FOR(d, rocdata.size())\r
151     {\r
152         int minCol = 128;\r
153         int color = (rocdata.size() == 1) ? 255 : (255 - minCol)*(rocdata.size() - d -1)/(rocdata.size()-1) + minCol;\r
154         color = 255 - color;\r
155 \r
156         std::vector<f32pair> data = rocdata[d];\r
157         if(!data.size()) continue;\r
158         std::sort(data.begin(), data.end());\r
159 \r
160         std::vector<fvec> allData;\r
161         FOR(i, data.size())\r
162         {\r
163             float thresh = data[i].first;\r
164             u32 tp = 0, fp = 0;\r
165             u32 fn = 0, tn = 0;\r
166 \r
167             FOR(j, data.size())\r
168             {\r
169                 if(data[j].second == 1)\r
170                 {\r
171                     if(data[j].first >= thresh) tp++;\r
172                     else fn++;\r
173                 }\r
174                 else\r
175                 {\r
176                     if(data[j].first >= thresh) fp++;\r
177                     else tn++;\r
178                 }\r
179             }\r
180             fVec val;\r
181             float fmeasure = 0;\r
182             if((fp+tn)>0 && (tp+fn)>0 && (tp+fp)>0)\r
183             {\r
184                 val=fVec(fp/float(fp+tn), 1 - tp/float(tp+fn));\r
185                 float precision = tp / float(tp+fp);\r
186                 float recall = tp /float(tp+fn);\r
187                 fmeasure = tp == 0 ? 0 : 2 * (precision * recall) / (precision + recall);\r
188             }\r
189 \r
190             fvec dat;\r
191             dat.push_back(val.x);\r
192             dat.push_back(val.y);\r
193             dat.push_back(data[i].first);\r
194             dat.push_back(fmeasure);\r
195             allData.push_back(dat);\r
196         }\r
197 \r
198         painter.setPen(QPen(QColor(color,color,color), 1.f));\r
199 \r
200         fVec pt1, pt2;\r
201         FOR(i, allData.size()-1)\r
202         {\r
203             pt1 = fVec(allData[i][0]*size.width(), allData[i][1]*size.height());\r
204             pt2 = fVec(allData[i+1][0]*size.width(), allData[i+1][1]*size.height());\r
205             painter.drawLine(QPointF(pt1.x+PAD, pt1.y+PAD),QPointF(pt2.x+PAD, pt2.y+PAD));\r
206         }\r
207         pt1 = fVec(0,size.width());\r
208         painter.drawLine(QPointF(pt1.x+PAD, pt1.y+PAD),QPointF(pt2.x+PAD, pt2.y+PAD));\r
209 \r
210         if(d < roclabels.size())\r
211         {\r
212             QPointF pos(3*size.width()/4,size.height() - (d+1)*16);\r
213             painter.drawText(pos,QString(roclabels[d]));\r
214         }\r
215     }\r
216 \r
217     font = painter.font();\r
218     font.setPointSize(10);\r
219     font.setBold(false);\r
220     font.setCapitalization(QFont::SmallCaps);\r
221     painter.setFont(font);\r
222     painter.setPen(Qt::black);\r
223     painter.drawText(0, 0, size.width(), 16, Qt::AlignCenter, "False Positives");\r
224     painter.translate(0, size.height());\r
225     painter.rotate(-90);\r
226     painter.drawText(0,0, size.height(), 16, Qt::AlignCenter, "True Positives");\r
227 \r
228     return pixmap;\r
229 }\r
230 \r
231 QPixmap BoxPlot(std::vector<fvec> allData, QSize size, float maxVal, float minVal)\r
232 {\r
233     QPixmap boxplot(size);\r
234     if(!allData.size()) return boxplot;\r
235     QBitmap bitmap;\r
236     bitmap.clear();\r
237     boxplot.setMask(bitmap);\r
238     boxplot.fill(Qt::transparent);\r
239     QPainter painter(&boxplot);\r
240 \r
241     //  painter.setRenderHint(QPainter::Antialiasing);\r
242 \r
243     FOR(d,allData.size())\r
244     {\r
245         fvec data = allData[d];\r
246         if(!data.size()) continue;\r
247         FOR(i, data.size()) maxVal = max(maxVal, data[i]);\r
248         FOR(i, data.size()) minVal = min(minVal, data[i]);\r
249     }\r
250     if(minVal == maxVal)\r
251     {\r
252         minVal = minVal/2;\r
253         minVal = minVal*3/2;\r
254     }\r
255 \r
256     FOR(d,allData.size())\r
257     {\r
258         int minCol = 70;\r
259         int color = (allData.size() == 1) ? minCol : (255-minCol) * d / allData.size() + minCol;\r
260 \r
261         fvec data = allData[d];\r
262         if(!data.size()) continue;\r
263         int hpad = 15 + (d*size.width()/(allData.size()));\r
264         int pad = -16;\r
265         int res = size.height()+2*pad;\r
266         int nanCount = 0;\r
267         FOR(i, data.size()) if(data[i] != data[i]) nanCount++;\r
268 \r
269         float mean = 0;\r
270         float sigma = 0;\r
271         FOR(i, data.size()) if(data[i]==data[i]) mean += data[i] / (data.size()-nanCount);\r
272         FOR(i, data.size()) if(data[i]==data[i]) sigma += powf(data[i]-mean,2);\r
273         sigma = sqrtf(sigma/(data.size()-nanCount));\r
274 \r
275         float edge = minVal;\r
276         float delta = maxVal - minVal;\r
277 \r
278         float top, bottom, median, quartLow, quartHi;\r
279         vector<float> outliers;\r
280         vector<float> sorted;\r
281 \r
282         if(data.size() > 1)\r
283         {\r
284             if(sigma==0)\r
285             {\r
286                 sorted = data;\r
287             }\r
288             else\r
289             {\r
290                 // we look for outliers using the 3*sigma rule\r
291                 FOR(i, data.size())\r
292                 {\r
293                     if(data[i]!=data[i]) continue;\r
294                     if (data[i] - mean < 3*sigma)\r
295                         sorted.push_back(data[i]);\r
296                     else outliers.push_back(data[i]);\r
297                 }\r
298             }\r
299             if(!sorted.size()) return boxplot;\r
300             sort(sorted.begin(), sorted.end());\r
301             int count = sorted.size();\r
302             int half = count/2;\r
303             bottom = sorted[0];\r
304             top = sorted[sorted.size()-1];\r
305 \r
306             median = count%2 ? sorted[half] : (sorted[half] + sorted[half - 1])/2;\r
307             if(count < 4)\r
308             {\r
309                 quartLow = bottom;\r
310                 quartHi = top;\r
311             }\r
312             else\r
313             {\r
314                 quartLow = half%2 ? sorted[half/2] : (sorted[half/2] + sorted[half/2 - 1])/2;\r
315                 quartHi = half%2 ? sorted[half*3/2] : (sorted[half*3/2] + sorted[half*3/2 - 1])/2;\r
316             }\r
317         }\r
318         else\r
319         {\r
320             top = bottom = median = quartLow = quartHi = data[0];\r
321         }\r
322 \r
323         QPointF bottomPoint = QPointF(0, size.height() - (int)((bottom-edge)/delta*res) + pad);\r
324         QPointF topPoint = QPointF(0, size.height() - (int)((top-edge)/delta*res) + pad);\r
325         QPointF medPoint = QPointF(0, size.height() - (int)((median-edge)/delta*res) + pad);\r
326         QPointF lowPoint = QPointF(0, size.height() - (int)((quartLow-edge)/delta*res) + pad);\r
327         QPointF highPoint = QPointF(0, size.height() - (int)((quartHi-edge)/delta*res) + pad);\r
328 \r
329         painter.setPen(Qt::black);\r
330         painter.drawLine(QPointF(hpad+35, bottomPoint.y()),     QPointF(hpad+65, bottomPoint.y()));\r
331 \r
332         painter.setPen(Qt::black);\r
333         painter.drawLine(QPointF(hpad+35, topPoint.y()), QPointF(hpad+65, topPoint.y()));\r
334 \r
335         painter.setPen(Qt::black);\r
336         painter.drawLine(QPointF(hpad+50, bottomPoint.y()),     QPointF(hpad+50, topPoint.y()));\r
337 \r
338         painter.setBrush(QColor(color,color,color));\r
339         painter.drawRect(hpad+30, lowPoint.y(), 40, highPoint.y() - lowPoint.y());\r
340 \r
341         painter.setPen(Qt::black);\r
342         painter.drawLine(QPointF(hpad+30, medPoint.y()),        QPointF(hpad+70, medPoint.y()));\r
343 \r
344         const char *longFormat = "%.3f";\r
345         const char *shortFormat = "%.0f";\r
346         const char *format = (maxVal - minVal) > 100 ? shortFormat : longFormat;\r
347         painter.setPen(Qt::black);\r
348         char text[255];\r
349         sprintf(text, format, median);\r
350         painter.drawText(QPointF(hpad-8,medPoint.y()+6), QString(text));\r
351         sprintf(text, format, top);\r
352         painter.drawText(QPointF(hpad+36,topPoint.y()-6), QString(text));\r
353         sprintf(text, format, bottom);\r
354         painter.drawText(QPointF(hpad+36,bottomPoint.y()+12), QString(text));\r
355     }\r
356     return boxplot;\r
357 }\r
358 \r
359 QPixmap Histogram(std::vector<fvec> allData, QSize size, float maxVal, float minVal)\r
360 {\r
361     QPixmap histogram(size);\r
362     if(!allData.size()) return histogram;\r
363     QBitmap bitmap;\r
364     bitmap.clear();\r
365     histogram.setMask(bitmap);\r
366     histogram.fill(Qt::transparent);\r
367     QPainter painter(&histogram);\r
368 \r
369     //  painter.setRenderHint(QPainter::Antialiasing);\r
370 \r
371     FOR(d,allData.size())\r
372     {\r
373         fvec data = allData[d];\r
374         if(!data.size()) continue;\r
375         FOR(i, data.size()) if(data[i]==data[i]) maxVal = max(maxVal, data[i]);\r
376         FOR(i, data.size()) if(data[i]==data[i]) minVal = min(minVal, data[i]);\r
377     }\r
378     if(minVal == maxVal)\r
379     {\r
380         minVal = minVal/2;\r
381         minVal = minVal*3/2;\r
382     }\r
383 \r
384     FOR(d,allData.size())\r
385     {\r
386         int minCol = 70;\r
387         int color = (allData.size() == 1) ? minCol : (255-minCol) * d / allData.size() + minCol;\r
388 \r
389         fvec data = allData[d];\r
390         if(!data.size()) continue;\r
391         int hpad = 15 + (d*size.width()/(allData.size()));\r
392         int pad = -16;\r
393         int res = size.height()+2*pad;\r
394 \r
395         int nanCount = 0;\r
396         FOR(i, data.size()) if(data[i] != data[i]) nanCount++;\r
397         float mean = 0;\r
398         float sigma = 0;\r
399         FOR(i, data.size()) if(data[i]==data[i]) mean += data[i] / (data.size()-nanCount);\r
400         FOR(i, data.size()) if(data[i]==data[i]) sigma += powf(data[i]-mean,2);\r
401         sigma = sqrtf(sigma/(data.size()-nanCount));\r
402 \r
403         float edge = minVal;\r
404         float delta = maxVal - minVal;\r
405         float bottom = 0;\r
406 \r
407         QPointF bottomPoint = QPointF(0, size.height() - (int)((bottom-edge)/delta*res) + pad);\r
408         QPointF topPoint = QPointF(0, size.height() - (int)((mean-edge)/delta*res) + pad);\r
409         QPointF plusPoint = QPointF(0, size.height() - (int)((mean+sigma-edge)/delta*res) + pad);\r
410         QPointF minusPoint = QPointF(0, size.height() - (int)((mean-sigma-edge)/delta*res) + pad);\r
411 \r
412         painter.setPen(Qt::black);\r
413         painter.drawLine(QPointF(hpad+35, bottomPoint.y()),     QPointF(hpad+65, bottomPoint.y()));\r
414 \r
415         painter.setPen(Qt::black);\r
416         painter.drawLine(QPointF(hpad+35, topPoint.y()), QPointF(hpad+65, topPoint.y()));\r
417 \r
418         painter.setBrush(QColor(color,color,color));\r
419         painter.drawRect(hpad+30, topPoint.y(), 40, bottomPoint.y()-topPoint.y());\r
420 \r
421         painter.setPen(Qt::black);\r
422         painter.drawLine(QPointF(hpad+50, plusPoint.y()), QPointF(hpad+50, minusPoint.y()));\r
423 \r
424         painter.setPen(Qt::black);\r
425         painter.drawLine(QPointF(hpad+40, plusPoint.y()),       QPointF(hpad+60, plusPoint.y()));\r
426 \r
427         painter.setPen(Qt::black);\r
428         painter.drawLine(QPointF(hpad+40, minusPoint.y()),      QPointF(hpad+60, minusPoint.y()));\r
429 \r
430         const char *longFormat = "%.3f";\r
431         const char *shortFormat = "%.0f";\r
432         const char *format = (maxVal - minVal) > 10 ? shortFormat : longFormat;\r
433         painter.setPen(Qt::black);\r
434         char text[255];\r
435         sprintf(text, format, mean);\r
436         painter.drawText(QPointF(hpad-8,topPoint.y()+6), QString(text));\r
437         sprintf(text, format, mean+sigma);\r
438         painter.drawText(QPointF(hpad+36,plusPoint.y()-6), QString(text));\r
439         sprintf(text, format, mean-sigma);\r
440         painter.drawText(QPointF(hpad+36,minusPoint.y()+12), QString(text));\r
441     }\r
442     return histogram;\r
443 }\r
444 \r
445 QPixmap RawData(std::vector<fvec> allData, QSize size, float maxVal, float minVal)\r
446 {\r
447     QPixmap rawData(size);\r
448     if(!allData.size()) return rawData;\r
449     QBitmap bitmap;\r
450     bitmap.clear();\r
451     rawData.setMask(bitmap);\r
452     rawData.fill(Qt::transparent);\r
453     QPainter painter(&rawData);\r
454 \r
455     painter.setRenderHint(QPainter::Antialiasing);\r
456 \r
457     FOR(d,allData.size())\r
458     {\r
459         fvec data = allData[d];\r
460         if(!data.size()) continue;\r
461         FOR(i, data.size()) if(data[i]==data[i]) maxVal = max(maxVal, data[i]);\r
462         FOR(i, data.size()) if(data[i]==data[i]) minVal = min(minVal, data[i]);\r
463     }\r
464     if(minVal == maxVal)\r
465     {\r
466         minVal = minVal/2;\r
467         minVal = minVal*3/2;\r
468     }\r
469 \r
470     FOR(d,allData.size())\r
471     {\r
472         int minCol = 70;\r
473         int color = (allData.size() == 1) ? minCol : (255-minCol) * d / allData.size() + minCol;\r
474 \r
475         fvec data = allData[d];\r
476         if(!data.size()) continue;\r
477         int hpad = 15 + (d*size.width()/(allData.size()));\r
478         int hsize = (size.width()/allData.size() - 15);\r
479         int pad = -16;\r
480         int res = size.height()+2*pad;\r
481         int nanCount = 0;\r
482         FOR(i, data.size()) if(data[i] != data[i]) nanCount++;\r
483 \r
484         float mean = 0;\r
485         float sigma = 0;\r
486         FOR(i, data.size()) if(data[i]==data[i]) mean += data[i] / (data.size()-nanCount);\r
487         FOR(i, data.size()) if(data[i]==data[i]) sigma += powf(data[i]-mean,2);\r
488         sigma = sqrtf(sigma/(data.size()-nanCount));\r
489 \r
490         float edge = minVal;\r
491         float delta = maxVal - minVal;\r
492 \r
493         QPointF topPoint = QPointF(0, size.height() - (int)((mean-edge)/delta*res) + pad);\r
494         QPointF plusPoint = QPointF(0, size.height() - (int)((mean+sigma-edge)/delta*res) + pad);\r
495         QPointF minusPoint = QPointF(0, size.height() - (int)((mean-sigma-edge)/delta*res) + pad);\r
496 \r
497         FOR(i, data.size())\r
498         {\r
499             QPointF point = QPointF(hpad + (drand48() - 0.5)*hsize/2 + hsize/2, size.height() - (int)((data[i]-edge)/delta*res) + pad);\r
500             painter.setPen(QPen(Qt::black, 0.5));\r
501             painter.setBrush(QColor(color,color,color));\r
502             painter.drawEllipse(point, 5, 5);\r
503         }\r
504         const char *longFormat = "%.3f";\r
505         const char *shortFormat = "%.0f";\r
506         const char *format = (maxVal - minVal) > 10 ? shortFormat : longFormat;\r
507         painter.setPen(Qt::black);\r
508         char text[255];\r
509         sprintf(text, format, mean);\r
510         painter.drawText(QPointF(hpad-8,topPoint.y()+6), QString(text));\r
511         sprintf(text, format, mean+sigma);\r
512         painter.drawText(QPointF(hpad-8,plusPoint.y()-6), QString(text));\r
513         sprintf(text, format, mean-sigma);\r
514         painter.drawText(QPointF(hpad-8,minusPoint.y()+12), QString(text));\r
515     }\r
516     return rawData;\r
517 }\r
518 \r
519 void Draw3DRegressor(GLWidget *glw, Regressor *regressor)\r
520 {\r
521     // we get the boundaries of our axes\r
522     vector<fvec> samples = glw->canvas->data->GetSamples();\r
523     int dim = glw->canvas->data->GetDimCount();\r
524     fvec mins(dim, FLT_MAX), maxes(dim, -FLT_MAX);\r
525     FOR(i, samples.size())\r
526     {\r
527         FOR(d, dim)\r
528         {\r
529             mins[d] = min(mins[d], samples[i][d]);\r
530             maxes[d] = max(maxes[d], samples[i][d]);\r
531         }\r
532     }\r
533     fvec center = (maxes + mins)*0.5f;\r
534     fvec dists = (maxes - mins)*0.5f;\r
535     float maxDist = dists[0];\r
536     FOR(d, dim) maxDist = max(dists[d], maxDist);\r
537     dists = fvec(dim, maxDist);\r
538     // we double them just to be safe\r
539     mins = center - dists;\r
540     maxes = center + dists;\r
541 \r
542     // and now we draw a nice grid\r
543     int xInd = 0, yInd = 1, zInd = regressor->outputDim;\r
544     if(regressor->outputDim == yInd) yInd = 2;\r
545     else if(regressor->outputDim == xInd) xInd = 2;\r
546     int xSteps = 128, ySteps = 128;\r
547     fvec point(dim,0);\r
548     fvec gridPoints(xSteps*ySteps);\r
549     qDebug() << "Generating regression surface";\r
550     FOR(y, ySteps)\r
551     {\r
552         point[yInd] = y/(float)ySteps*(maxes[yInd]-mins[yInd]) + mins[yInd];\r
553         FOR(x, xSteps)\r
554         {\r
555             point[xInd] = x/(float)xSteps*(maxes[xInd]-mins[xInd]) + mins[xInd];\r
556             // we get the value of the regressor at the coordinates in the meshgrid\r
557             fvec res = regressor->Test(point);\r
558             gridPoints[x+y*xSteps] = res[0];\r
559         }\r
560     }\r
561     qDebug() << "Creating GLObject structure";\r
562     GLObject o = GenerateMeshGrid(gridPoints, xSteps, mins, maxes, xInd, yInd, zInd);\r
563     qDebug() << "Done.";\r
564     o.style = "smooth,transparent";\r
565     o.style += QString(",isolines:%1").arg(zInd);\r
566     o.style += ",blurry:3,color:1.0:1.0:1.0:0.4";\r
567     glw->objects.push_back(o);\r
568 }\r
569 \r
570 void Draw3DClassifier(GLWidget *glw, Classifier *classifier)\r
571 {\r
572     // we get the boundaries of our axes\r
573     vector<fvec> samples = glw->canvas->data->GetSamples();\r
574     int dim = glw->canvas->data->GetDimCount();\r
575     fvec mins(dim, FLT_MAX), maxes(dim, -FLT_MAX);\r
576     int xIndex = glw->canvas->xIndex;\r
577     int yIndex = glw->canvas->yIndex;\r
578     int zIndex = glw->canvas->zIndex;\r
579     FOR(i, samples.size())\r
580     {\r
581         FOR(d, dim)\r
582         {\r
583             mins[d] = min(mins[d], samples[i][d]);\r
584             maxes[d] = max(maxes[d], samples[i][d]);\r
585         }\r
586     }\r
587     fvec center = (maxes + mins)*0.5f;\r
588     fvec dists = (maxes - mins)*0.5f;\r
589     float maxDist = dists[0];\r
590     FOR(d, dim) maxDist = max(dists[d], maxDist);\r
591     dists = fvec(dim, maxDist);\r
592     // we double them just to be safe\r
593     mins = center - dists*2;\r
594     maxes = center + dists*2;\r
595 \r
596     bool bMultiClass = false;\r
597     // and now we draw a volume\r
598     int steps = 64;\r
599     fvec sample(dim);\r
600     float *minVals = new float[steps];\r
601     float *maxVals = new float[steps];\r
602     double *values = new double[steps*steps*steps];\r
603     printf("Generating volumetric data: ");\r
604     FOR(y, steps)\r
605     {\r
606         //        values[y] = new double[steps*steps];\r
607         minVals[y] = FLT_MAX;\r
608         maxVals[y] = -FLT_MAX;\r
609         sample[yIndex] = y/(float)(steps)*(maxes[yIndex]-mins[yIndex]) + mins[yIndex];\r
610         FOR(z, steps)\r
611         {\r
612             sample[zIndex] = z/(float)steps*(maxes[zIndex]-mins[zIndex]) + mins[zIndex];\r
613             FOR(x, steps)\r
614             {\r
615                 sample[xIndex] = x/(float)steps*(maxes[xIndex]-mins[xIndex]) + mins[xIndex];\r
616                 if(classifier->IsMultiClass())\r
617                 {\r
618                     fvec res = classifier->TestMulti(sample);\r
619                     if(res.size() == 1)\r
620                     {\r
621                         values[x + (y + z*steps)*steps] = res[0];\r
622                     }\r
623                     else\r
624                     {\r
625                         // we mostly want to know if the sample is close to the boundary\r
626                         float maxVal = res[0];\r
627                         int maxInd = 0;\r
628                         FOR(d, res.size())\r
629                         {\r
630                             if(maxVal < res[d])\r
631                             {\r
632                                 maxInd = d;\r
633                                 maxVal = res[d];\r
634                             }\r
635                         }\r
636                         // we keep the class with the highest score\r
637                         values[x + (y + z*steps)*steps] = maxInd;\r
638                         bMultiClass = true;\r
639                     }\r
640                 }\r
641                 else\r
642                 {\r
643                     float res = classifier->Test(sample);\r
644                     values[x + (y + z*steps)*steps] = res;\r
645                 }\r
646             }\r
647         }\r
648     }\r
649     printf("done.\n");\r
650     fflush(stdout);\r
651 \r
652     if(bMultiClass)\r
653     {\r
654         int classCount = DatasetManager::GetClassCount(glw->canvas->data->GetLabels());\r
655         FOR(c,classCount-1)\r
656         {\r
657             gridT valueGrid(0.f, steps, steps, steps);\r
658             /*\r
659             FOR(i, steps*steps*steps)\r
660             {\r
661                 if(values[i] == c) valueGrid[i] = 1.f;\r
662                 else\r
663                 {\r
664                     valueGrid[i] = 2.f;\r
665                     for(int c2=c+1; c2<classCount; c2++)\r
666                         if(values[i] == c2) valueGrid[i] = -1.f;\r
667                 }\r
668             }\r
669             */\r
670             FOR(i, steps*steps*steps) valueGrid[i] = values[i] == c ? 1.f : -1.f;\r
671             FOR(d, 3) valueGrid.unit[d] = (maxes[d] - mins[d])/steps;   /* length of a single edge in each dimension*/\r
672             FOR(d, 3) valueGrid.size[d] = maxes[d] - mins[d];           /* length of entire grid in each dimension */\r
673             FOR(d, 3) valueGrid.org[d] = mins[d];                       /* the origin of the grid i.e. coords of (0,0,0) */\r
674             FOR(d, 3) valueGrid.center[d] = (maxes[d] + mins[d])/2;     /* coords of center of grid */\r
675 \r
676             surfaceT surf;\r
677             float surfThreshold = 0.f;\r
678             unsigned int surfIType = JACSurfaceTypes::SURF_CONTOUR;\r
679             printf("Generating isosurfaces: ");\r
680             fflush(stdout);\r
681             JACMakeSurface(surf, surfIType, valueGrid, surfThreshold);\r
682             JACSmoothSurface(surf);\r
683             JACSmoothSurface(surf);\r
684             JACSmoothSurface(surf);\r
685             //surf.Reduce(0.05);\r
686             printf("done.\n");\r
687             fflush(stdout);\r
688 \r
689             printf("Generating mesh: ");\r
690 \r
691             GLObject o;\r
692 \r
693             std::vector<float> &vertices = surf.vertices;\r
694             //std::vector<float> &normals = surf.normals;\r
695             for (int i=0; i<surf.nconn; i += 3)\r
696             {\r
697                 int index = surf.triangles[i];\r
698                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
699                 index = surf.triangles[i+1];\r
700                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
701                 index = surf.triangles[i+2];\r
702                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
703             }\r
704 \r
705             QColor color = SampleColor[(c+1)%SampleColorCnt];\r
706             o.objectType = "Surfaces";\r
707             o.style = "smooth,transparent,blurry";\r
708             o.style += QString("color:%1:%2:%3:0.4").arg(color.redF()).arg(color.greenF()).arg(color.blueF());\r
709             o.style += QString(",offset:%1").arg(c*0.5f,0,'f',2);\r
710             glw->objects.push_back(o);\r
711             printf("done.\n");\r
712             fflush(stdout);\r
713         }\r
714     }\r
715     else\r
716     {\r
717         gridT valueGrid(0.f, steps, steps, steps);\r
718         FOR(i, steps*steps*steps) valueGrid[i] = values[i];\r
719         FOR(d, 3) valueGrid.unit[d] = (maxes[d] - mins[d])/steps;   /* length of a single edge in each dimension*/\r
720         FOR(d, 3) valueGrid.size[d] = maxes[d] - mins[d];           /* length of entire grid in each dimension */\r
721         FOR(d, 3) valueGrid.org[d] = mins[d];                       /* the origin of the grid i.e. coords of (0,0,0) */\r
722         FOR(d, 3) valueGrid.center[d] = (maxes[d] + mins[d])/2;     /* coords of center of grid */\r
723 \r
724         surfaceT surf;\r
725         float surfThreshold = 0.f;\r
726         unsigned int surfIType = JACSurfaceTypes::SURF_CONTOUR;\r
727         printf("Generating isosurfaces: ");\r
728         fflush(stdout);\r
729         JACMakeSurface(surf, surfIType, valueGrid, surfThreshold);\r
730         JACSmoothSurface(surf);\r
731         JACSmoothSurface(surf);\r
732         JACSmoothSurface(surf);\r
733         //surf.Reduce(0.05);\r
734         printf("done.\n");\r
735         fflush(stdout);\r
736 \r
737         printf("Generating mesh: ");\r
738 \r
739         GLObject o;\r
740 \r
741         std::vector<float> &vertices = surf.vertices;\r
742         //std::vector<float> &normals = surf.normals;\r
743         for (int i=0; i<surf.nconn; i += 3)\r
744         {\r
745             int index = surf.triangles[i];\r
746             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
747             index = surf.triangles[i+1];\r
748             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
749             index = surf.triangles[i+2];\r
750             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
751         }\r
752 \r
753         o.objectType = "Surfaces";\r
754         o.style = "smooth,transparent,blurry:1,color:0:0:0:0.3";\r
755         glw->objects.push_back(o);\r
756         printf("done.\n");\r
757         fflush(stdout);\r
758     }\r
759 }\r
760 \r
761 void Draw3DClusterer(GLWidget *glw, Clusterer *clusterer)\r
762 {\r
763     // we get the boundaries of our axes\r
764     vector<fvec> samples = glw->canvas->data->GetSamples();\r
765     int dim = glw->canvas->data->GetDimCount();\r
766     fvec mins(dim, FLT_MAX), maxes(dim, -FLT_MAX);\r
767     int xIndex = glw->canvas->xIndex;\r
768     int yIndex = glw->canvas->yIndex;\r
769     int zIndex = glw->canvas->zIndex;\r
770     FOR(i, samples.size())\r
771     {\r
772         FOR(d, dim)\r
773         {\r
774             mins[d] = min(mins[d], samples[i][d]);\r
775             maxes[d] = max(maxes[d], samples[i][d]);\r
776         }\r
777     }\r
778     fvec center = (maxes + mins)*0.5f;\r
779     fvec dists = (maxes - mins)*0.5f;\r
780     float maxDist = dists[0];\r
781     FOR(d, dim) maxDist = max(dists[d], maxDist);\r
782     dists = fvec(dim, maxDist);\r
783     // we double them just to be safe\r
784     mins = center - dists*2;\r
785     maxes = center + dists*2;\r
786 \r
787     bool bMultiClass = false;\r
788     // and now we draw a volume\r
789     int steps = 64;\r
790     fvec sample(dim);\r
791     float *minVals = new float[steps];\r
792     float *maxVals = new float[steps];\r
793     double *values = new double[steps*steps*steps];\r
794     printf("Generating volumetric data: ");\r
795     bool bOneClass = true;\r
796     int clusterCount = clusterer->NbClusters();\r
797     FOR(y, steps)\r
798     {\r
799         //        values[y] = new double[steps*steps];\r
800         minVals[y] = FLT_MAX;\r
801         maxVals[y] = -FLT_MAX;\r
802         sample[yIndex] = y/(float)(steps)*(maxes[yIndex]-mins[yIndex]) + mins[yIndex];\r
803         FOR(z, steps)\r
804         {\r
805             sample[zIndex] = z/(float)steps*(maxes[zIndex]-mins[zIndex]) + mins[zIndex];\r
806             FOR(x, steps)\r
807             {\r
808                 sample[xIndex] = x/(float)steps*(maxes[xIndex]-mins[xIndex]) + mins[xIndex];\r
809                 fvec res = clusterer->Test(sample);\r
810                 if(res.size() == 1) values[x + (y + z*steps)*steps] = res[0];\r
811                 else\r
812                 {\r
813 \r
814                     float maxVal = res[0];\r
815                     int maxInd = 0;\r
816                     FOR(d, res.size())\r
817                     {\r
818                         if(maxVal < res[d])\r
819                         {\r
820                             maxInd = d;\r
821                             maxVal = res[d];\r
822                         }\r
823                     }\r
824                     // we keep the class with the highest score\r
825                     values[x + (y + z*steps)*steps] = maxInd;\r
826                     bOneClass = false;\r
827                 }\r
828             }\r
829         }\r
830     }\r
831     printf("done.\n");\r
832     fflush(stdout);\r
833 \r
834     if(!bOneClass)\r
835     {\r
836         FOR(c,clusterCount-1)\r
837         {\r
838             gridT valueGrid(0.f, steps, steps, steps);\r
839             FOR(i, steps*steps*steps)\r
840             {\r
841                 if(values[i] == c) valueGrid[i] = 1.f;\r
842                 else\r
843                 {\r
844                     valueGrid[i] = 2.f;\r
845                     for(int c2=c+1; c2<clusterCount; c2++)\r
846                         if(values[i] == c2) valueGrid[i] = -1.f;\r
847                 }\r
848             }\r
849             FOR(d, 3) valueGrid.unit[d] = (maxes[d] - mins[d])/steps;   /* length of a single edge in each dimension*/\r
850             FOR(d, 3) valueGrid.size[d] = maxes[d] - mins[d];           /* length of entire grid in each dimension */\r
851             FOR(d, 3) valueGrid.org[d] = mins[d];                       /* the origin of the grid i.e. coords of (0,0,0) */\r
852             FOR(d, 3) valueGrid.center[d] = (maxes[d] + mins[d])/2;     /* coords of center of grid */\r
853 \r
854             surfaceT surf;\r
855             float surfThreshold = 0.f;\r
856             unsigned int surfIType = JACSurfaceTypes::SURF_CONTOUR;\r
857             printf("Generating isosurfaces: ");\r
858             fflush(stdout);\r
859             JACMakeSurface(surf, surfIType, valueGrid, surfThreshold);\r
860             JACSmoothSurface(surf);\r
861             JACSmoothSurface(surf);\r
862             JACSmoothSurface(surf);\r
863             //surf.Reduce(0.05);\r
864             printf("done.\n");\r
865             fflush(stdout);\r
866 \r
867             printf("Generating mesh: ");\r
868 \r
869             GLObject o;\r
870 \r
871             std::vector<float> &vertices = surf.vertices;\r
872             //std::vector<float> &normals = surf.normals;\r
873             for (int i=0; i<surf.nconn; i += 3)\r
874             {\r
875                 int index = surf.triangles[i];\r
876                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
877                 index = surf.triangles[i+1];\r
878                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
879                 index = surf.triangles[i+2];\r
880                 o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
881             }\r
882 \r
883             QColor color = SampleColor[(c+1)%SampleColorCnt];\r
884             o.objectType = "Surfaces";\r
885             o.style = "smooth,transparent,blurry:1";\r
886             o.style += QString("color:%1:%2:%3:0.3").arg(color.redF()).arg(color.greenF()).arg(color.blueF());\r
887             o.style += QString(",offset:%1").arg((float)c,0,'f',2);\r
888             glw->objects.push_back(o);\r
889             printf("done.\n");\r
890             fflush(stdout);\r
891         }\r
892     }\r
893     else\r
894     {\r
895         gridT valueGrid(0.f, steps, steps, steps);\r
896         FOR(i, steps*steps*steps) valueGrid[i] = values[i];\r
897         FOR(d, 3) valueGrid.unit[d] = (maxes[d] - mins[d])/steps;   /* length of a single edge in each dimension*/\r
898         FOR(d, 3) valueGrid.size[d] = maxes[d] - mins[d];           /* length of entire grid in each dimension */\r
899         FOR(d, 3) valueGrid.org[d] = mins[d];                       /* the origin of the grid i.e. coords of (0,0,0) */\r
900         FOR(d, 3) valueGrid.center[d] = (maxes[d] + mins[d])/2;     /* coords of center of grid */\r
901 \r
902         surfaceT surf;\r
903         float surfThreshold = 0.f;\r
904         unsigned int surfIType = JACSurfaceTypes::SURF_CONTOUR;\r
905         printf("Generating isosurfaces: ");\r
906         fflush(stdout);\r
907         JACMakeSurface(surf, surfIType, valueGrid, surfThreshold);\r
908         JACSmoothSurface(surf);\r
909         JACSmoothSurface(surf);\r
910         JACSmoothSurface(surf);\r
911         //surf.Reduce(0.05);\r
912         printf("done.\n");\r
913         fflush(stdout);\r
914 \r
915         printf("Generating mesh: ");\r
916 \r
917         GLObject o;\r
918 \r
919         std::vector<float> &vertices = surf.vertices;\r
920         //std::vector<float> &normals = surf.normals;\r
921         for (int i=0; i<surf.nconn; i += 3)\r
922         {\r
923             int index = surf.triangles[i];\r
924             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
925             index = surf.triangles[i+1];\r
926             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
927             index = surf.triangles[i+2];\r
928             o.vertices.append(QVector3D(vertices[index*3],vertices[index*3+1],vertices[index*3+2]));\r
929         }\r
930 \r
931         o.objectType = "Surfaces";\r
932         o.style = "smooth,transparent,blurry:1,color:1:1:1:0.4";\r
933         glw->objects.push_back(o);\r
934         printf("done.\n");\r
935         fflush(stdout);\r
936     }\r
937 }\r
938 \r
939 struct Streamline\r
940 {\r
941     std::vector<fvec> trajectory;\r
942     int cluster;\r
943     int length;\r
944 public:\r
945     Streamline() : length(0), cluster(0) {}\r
946     Streamline(std::vector<fvec> trajectory)\r
947         : trajectory(trajectory), length(trajectory.size()), cluster(0) {}\r
948     fvec &operator[](int i) {return trajectory[i];}\r
949     fvec &operator()(int i) {return trajectory[i];}\r
950     void push_back(const fvec point) {trajectory.push_back(point); length++;}\r
951     void clear() {trajectory.clear(); length = 0; cluster = 0;}\r
952     int size() {return length;}\r
953     fvec &back() {return trajectory.back();}\r
954     fvec &front() {return trajectory.front();}\r
955 };\r
956 \r
957 inline float StreamDistance(Streamline &s, vector<fvec> &mean)\r
958 {\r
959     float res = 0;\r
960     int dim = mean[0].size();\r
961     int length = min((int)mean.size(), s.length);\r
962     FOR(i, length)\r
963     {\r
964         FOR(d, dim) res += (mean[i][d]-s[i][d])*(mean[i][d]-s[i][d]);\r
965     }\r
966     return res;\r
967 }\r
968 \r
969 ivec ClusterStreams(std::vector<Streamline> streams, int nbClusters, int maxIterations=5)\r
970 {\r
971     if(!streams.size()) return ivec();\r
972     int count = streams.size();\r
973     nbClusters = min(count, nbClusters);\r
974     vector< vector<fvec> > means(nbClusters); // the clusters 'mean' trajectories\r
975     ivec clusters(count); // the responsiblity (which cluster each stream belongs to\r
976 \r
977     // first thing we do is to invert all streams so that we can follow them backwards\r
978     FOR(i, count)\r
979     {\r
980         std::reverse(streams[i].trajectory.begin(), streams[i].trajectory.end());\r
981     }\r
982 \r
983     // we initialize by picking actual trajectories\r
984     FOR(i, nbClusters)\r
985     {\r
986         int index = i * count / nbClusters;\r
987         means[i] = streams[index].trajectory;\r
988     }\r
989 \r
990     FOR(it, maxIterations)\r
991     {\r
992         // we recompute the responsiblity for each cluster\r
993         qDebug() << "E Step" << it;\r
994         ivec newClusters(count,0);\r
995         FOR(i, count)\r
996         {\r
997             float minDist = FLT_MAX;\r
998             int minInd = 0;\r
999             FOR(j, nbClusters)\r
1000             {\r
1001                 float dist = StreamDistance(streams[i], means[j]);\r
1002                 if(dist < minDist)\r
1003                 {\r
1004                     minInd = j;\r
1005                     minDist = dist;\r
1006                 }\r
1007             }\r
1008             newClusters[i] = minInd;\r
1009         }\r
1010         if(clusters == newClusters) break; // we've converged!\r
1011         clusters = newClusters;\r
1012 \r
1013         // and now we compute the new means\r
1014         qDebug() << "M Step" << it;\r
1015         ivec counts(nbClusters,0);\r
1016         vector<ivec> meanCounts(nbClusters);\r
1017         FOR(i, count)\r
1018         {\r
1019             int c = clusters[i];\r
1020             if(!counts[c])\r
1021             {\r
1022                 means[c] = streams[i].trajectory;\r
1023                 meanCounts[c].resize(streams[i].length,1);\r
1024             }\r
1025             else\r
1026             {\r
1027                 FOR(j, streams[i].size())\r
1028                 {\r
1029                     if(j < means[c].size())\r
1030                     {\r
1031                         means[c][j] += streams[i][j];\r
1032                         meanCounts[c][j]++;\r
1033                     }\r
1034                     else\r
1035                     {\r
1036                         means[c].push_back(streams[i][j]);\r
1037                         meanCounts[c].push_back(1);\r
1038                     }\r
1039                 }\r
1040             }\r
1041             counts[c]++;\r
1042         }\r
1043         FOR(c, nbClusters)\r
1044         {\r
1045             FOR(j, means[c].size())\r
1046             {\r
1047                 means[c][j] /= meanCounts[c][j];\r
1048             }\r
1049         }\r
1050     }\r
1051     return clusters;\r
1052 }\r
1053 \r
1054 unsigned int tessIndices[8][3] = {{0,1,2},{0,2,3},{0,3,4},{0,4,1},{5,2,1},{5,3,2},{5,4,3},{5,1,4}};\r
1055 float tessVerts[6][3] = {{0,0,-1},{1,0,0},{0,-1,0},{-1,0,0},{0,1,0},{0,0,1}};\r
1056 void normalize_vert(float *a)\r
1057 {\r
1058     float d=sqrtf(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);\r
1059     d = 1.f/d;\r
1060     a[0]*=d; a[1]*=d; a[2]*=d;\r
1061 }\r
1062 \r
1063 void draw_recursive_tri(float *a,float *b,float *c,unsigned int div, vector<fvec> &vertices)\r
1064 {\r
1065     if (div==0)\r
1066     {\r
1067         fvec v(3);\r
1068         v[0] = (a[0]+b[0]+c[0])/3.f;\r
1069         v[1] = (a[1]+b[1]+c[1])/3.f;\r
1070         v[2] = (a[2]+b[2]+c[2])/3.f;\r
1071         vertices.push_back(v);\r
1072     }\r
1073     else\r
1074     {\r
1075         register unsigned int i;\r
1076         float ab[3],ac[3],bc[3];\r
1077         for (i=0; i<3; i++)\r
1078         {\r
1079             ab[i]=(a[i]+b[i])/2.0f;\r
1080             ac[i]=(a[i]+c[i])/2.0f;\r
1081             bc[i]=(b[i]+c[i])/2.0f;\r
1082         }\r
1083         normalize_vert(ab);\r
1084         normalize_vert(ac);\r
1085         normalize_vert(bc);\r
1086         draw_recursive_tri(a,ab,ac,div-1,vertices);\r
1087         draw_recursive_tri(b,bc,ab,div-1,vertices);\r
1088         draw_recursive_tri(c,ac,bc,div-1,vertices);\r
1089         draw_recursive_tri(ab,bc,ac,div-1,vertices);\r
1090     }\r
1091 }\r
1092 \r
1093 float **tessellatedSphere(unsigned int detail)\r
1094 {\r
1095     vector<fvec> vertices;\r
1096     FOR(i, 8)\r
1097         draw_recursive_tri\r
1098                 (\r
1099                     tessVerts[tessIndices[i][0]],\r
1100                     tessVerts[tessIndices[i][1]],\r
1101                     tessVerts[tessIndices[i][2]],\r
1102                     detail,vertices\r
1103                     );\r
1104     float **tess = new float*[vertices.size()];\r
1105     FOR(i, vertices.size())\r
1106     {\r
1107         tess[i] = new float[3];\r
1108         tess[i][0] = vertices[i][0];\r
1109         tess[i][1] = vertices[i][1];\r
1110         tess[i][2] = vertices[i][2];\r
1111     }\r
1112     return tess;\r
1113 }\r
1114 \r
1115 // level 1: 32, level 2: 128\r
1116 float **tesssphere = 0;\r
1117 int tesssize = 32;\r
1118 inline int binFromVector(float *v)\r
1119 {\r
1120     if(!tesssphere) tesssphere = tessellatedSphere(1);\r
1121     int bin = 0;\r
1122     float minDist = FLT_MAX;\r
1123     FOR(i, tesssize)\r
1124     {\r
1125         float *t = tesssphere[i];\r
1126         float dist = (t[0]-v[0])*(t[0]-v[0]) + (t[1]-v[1])*(t[1]-v[1]) + (t[2]-v[2])*(t[2]-v[2]);\r
1127         if(dist < minDist)\r
1128         {\r
1129             minDist = dist;\r
1130             bin = i;\r
1131         }\r
1132     }\r
1133     return bin;\r
1134 }\r
1135 \r
1136 fvec ComputeDynamicalEntropy(Dynamical *dynamical, fvec mins, fvec maxes, int gridSteps = 64, int hSteps = 16)\r
1137 {\r
1138     // we begin by generating a dense grid of values\r
1139     qDebug() << "dumping vectors to memory";\r
1140     vector<fvec> grid(gridSteps*gridSteps*gridSteps);\r
1141     fvec sample(3);\r
1142     FOR(z, gridSteps)\r
1143     {\r
1144         sample[2] = z / (float)gridSteps * (maxes[2]-mins[2]) + mins[2];\r
1145         FOR(y, gridSteps)\r
1146         {\r
1147             sample[1] = y / (float)gridSteps * (maxes[1]-mins[1]) + mins[1];\r
1148             FOR(x, gridSteps)\r
1149             {\r
1150                 sample[0] = x / (float)gridSteps * (maxes[0]-mins[0]) + mins[0];\r
1151                 fvec res = dynamical->Test(sample);\r
1152 //                res = res/sqrtf(res*res); // we normalize it\r
1153                 grid[x + (y + z*gridSteps)*gridSteps] = res;\r
1154             }\r
1155         }\r
1156     }\r
1157 \r
1158     if(!tesssphere) tesssphere = tessellatedSphere(1);\r
1159 \r
1160     // bins: binCount*binCount*binCount\r
1161     int binCount = tesssize;\r
1162     int ratio = gridSteps/hSteps;\r
1163     fvec H(hSteps*hSteps*hSteps);\r
1164 \r
1165     FOR(i, hSteps)\r
1166     {\r
1167         FOR(j, hSteps)\r
1168         {\r
1169             FOR(k, hSteps)\r
1170             {\r
1171                 int bins[32];\r
1172                 FOR(d, 32) bins[d] = 0;\r
1173                 // we get the histogram for the current subcube\r
1174                 FOR(z,ratio)\r
1175                 {\r
1176                     FOR(y,ratio)\r
1177                     {\r
1178                         FOR(x,ratio)\r
1179                         {\r
1180                             float *val = &grid[(x + k*ratio) + (y+j*ratio + (z+i*ratio)*gridSteps)*gridSteps][0];\r
1181                             int bin = binFromVector(val);\r
1182                             bins[bin] += 1;\r
1183                         }\r
1184                     }\r
1185                 }\r
1186                 float sum = ratio*ratio*ratio;\r
1187                 float entropy = 0;\r
1188                 FOR(d, binCount)\r
1189                 {\r
1190                     if(bins[d] == 0) continue;\r
1191                     float p = bins[d]/sum;\r
1192                     float e = p*log2(p);\r
1193                     entropy -= e;\r
1194                 }\r
1195                 H[k + (j + i*hSteps)*hSteps] = entropy;\r
1196             }\r
1197         }\r
1198     }\r
1199     return H;\r
1200 }\r
1201 \r
1202 GLObject DrawEntropyField(fvec H, float minv, float maxv, int hSteps)\r
1203 {\r
1204     qDebug() << "drawing entropy field";\r
1205     GLObject o;\r
1206     o.objectType = "Dynamize,Surfaces,quads";\r
1207     o.style = "smooth";\r
1208     float minH = FLT_MAX, maxH = -FLT_MAX;\r
1209     FOR(i, hSteps*hSteps*hSteps)\r
1210     {\r
1211         minH = min(minH, H[i]);\r
1212         maxH = max(maxH, H[i]);\r
1213     }\r
1214 \r
1215     FOR(i, hSteps)\r
1216     {\r
1217         float z = i/(float)hSteps*(maxv-minv) + minv;\r
1218         FOR(j, hSteps)\r
1219         {\r
1220             float y = j/(float)hSteps*(maxv-minv) + minv;\r
1221             FOR(k, hSteps)\r
1222             {\r
1223                 float v = H[k + (j + i*hSteps)*hSteps];\r
1224                 v = (v-minH)/(maxH-minH);\r
1225                 if(v < 0.01) continue;\r
1226                 QColor color(Canvas::GetColorMapValue(v,2));\r
1227                 float x = k/(float)hSteps*(maxv-minv) + minv;\r
1228                 float r = 0.02;\r
1229                 o.vertices.push_back(QVector3D(x-r,y-r,z-r));\r
1230                 o.vertices.push_back(QVector3D(x+r,y-r,z-r));\r
1231                 o.vertices.push_back(QVector3D(x+r,y+r,z-r));\r
1232                 o.vertices.push_back(QVector3D(x-r,y+r,z-r));\r
1233                 FOUR(o.normals.push_back(QVector3D(0,0,1)));\r
1234                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1235                 o.vertices.push_back(QVector3D(x-r,y-r,z+r));\r
1236                 o.vertices.push_back(QVector3D(x+r,y-r,z+r));\r
1237                 o.vertices.push_back(QVector3D(x+r,y+r,z+r));\r
1238                 o.vertices.push_back(QVector3D(x-r,y+r,z+r));\r
1239                 FOUR(o.normals.push_back(QVector3D(0,0,-1)));\r
1240                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1241                 o.vertices.push_back(QVector3D(x-r,y-r,z-r));\r
1242                 o.vertices.push_back(QVector3D(x-r,y-r,z+r));\r
1243                 o.vertices.push_back(QVector3D(x-r,y+r,z+r));\r
1244                 o.vertices.push_back(QVector3D(x-r,y+r,z-r));\r
1245                 FOUR(o.normals.push_back(QVector3D(1,0,0)));\r
1246                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1247                 o.vertices.push_back(QVector3D(x+r,y-r,z-r));\r
1248                 o.vertices.push_back(QVector3D(x+r,y-r,z+r));\r
1249                 o.vertices.push_back(QVector3D(x+r,y+r,z+r));\r
1250                 o.vertices.push_back(QVector3D(x+r,y+r,z-r));\r
1251                 FOUR(o.normals.push_back(QVector3D(-1,0,0)));\r
1252                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1253                 o.vertices.push_back(QVector3D(x-r,y-r,z-r));\r
1254                 o.vertices.push_back(QVector3D(x-r,y-r,z+r));\r
1255                 o.vertices.push_back(QVector3D(x+r,y-r,z+r));\r
1256                 o.vertices.push_back(QVector3D(x+r,y-r,z-r));\r
1257                 FOUR(o.normals.push_back(QVector3D(0,1,0)));\r
1258                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1259                 o.vertices.push_back(QVector3D(x-r,y+r,z-r));\r
1260                 o.vertices.push_back(QVector3D(x-r,y+r,z+r));\r
1261                 o.vertices.push_back(QVector3D(x+r,y+r,z+r));\r
1262                 o.vertices.push_back(QVector3D(x+r,y+r,z-r));\r
1263                 FOUR(o.normals.push_back(QVector3D(0,-1,0)));\r
1264                 FOUR(o.colors.push_back(QVector4D(color.redF(), color.greenF(), color.blueF(), 1.f)));\r
1265             }\r
1266         }\r
1267     }\r
1268     return o;\r
1269 }\r
1270 \r
1271 GLObject DrawStreamTubes(vector<Streamline> streams, float diff, float maxSpeed, int xInd, int yInd, int zInd)\r
1272 {\r
1273     GLObject o;\r
1274     o.objectType = "Dynamize,Surfaces,quads";\r
1275     o.style = "smooth";\r
1276     FOR(i, streams.size())\r
1277     {\r
1278         if(streams[i].length < 2) continue;\r
1279         float radius = diff*0.001;\r
1280         QVector3D p1, p2;\r
1281         vector<QVector3D> oldCircle;\r
1282         vector<QVector3D> oldNormals;\r
1283         QVector3D oldP, pos;\r
1284         float oldSpeed = 0.f;\r
1285         FOR(j, streams[i].length-1)\r
1286         {\r
1287             QVector3D d;\r
1288             p1 = QVector3D(streams[i][j][xInd],streams[i][j][yInd],streams[i][j][zInd]);\r
1289             p2 = QVector3D(streams[i][j+1][xInd],streams[i][j+1][yInd],streams[i][j+1][zInd]);\r
1290             QVector3D dn = (p2-p1);\r
1291             d = dn.normalized();\r
1292             float speed = dn.length();\r
1293             if(j==0)\r
1294             {\r
1295                 if(d.z() != 0)\r
1296                 {\r
1297                     float x=drand48(), y=drand48();\r
1298                     float z=-(x*d.x() + y*d.y())/d.z();\r
1299                     pos = QVector3D(x,y,z).normalized();\r
1300                 }\r
1301                 else\r
1302                 {\r
1303                     float z=drand48(), y=drand48();\r
1304                     float x=-(y*d.y() + z*d.z())/d.x();\r
1305                     pos = QVector3D(x,y,z).normalized();\r
1306                 }\r
1307             }\r
1308             else\r
1309             {\r
1310                 float dDot = QVector3D::dotProduct(d, pos);\r
1311                 pos = (pos - d*dDot).normalized();\r
1312             }\r
1313 \r
1314             // we want to generate a circle orthogonal to the segment direction\r
1315             int segs=10;\r
1316             QQuaternion quat(M_PI, d);\r
1317             //QQuaternion quat(M_PI/(float)segs, d);\r
1318             vector<QVector3D> circle(segs);\r
1319             vector<QVector3D> normals(segs);\r
1320             FOR(s, segs)\r
1321             {\r
1322                 circle[s] = pos*radius*(speed*100.f);\r
1323                 normals[s] = pos;\r
1324                 pos = quat.rotatedVector(pos).normalized();\r
1325             }\r
1326             if(j > 0)\r
1327             {\r
1328                 // we look for the closest point to each of the circle\r
1329 \r
1330                 int offset = 0;\r
1331 \r
1332                 float minDist=FLT_MAX;\r
1333                 int minInd = 0;\r
1334                 FOR(s, segs)\r
1335                 {\r
1336                     QVector3D v = circle[0]-oldCircle[s];\r
1337                     float dist = v.lengthSquared();\r
1338                     if(dist < minDist)\r
1339                     {\r
1340                         minInd = s;\r
1341                         minDist = dist;\r
1342                     }\r
1343                 }\r
1344                 offset = minInd;\r
1345                 QColor cVal(Canvas::GetColorMapValue((speed+oldSpeed)*0.5/maxSpeed,2)); // jet color\r
1346                 QVector4D color(cVal.redF(),cVal.greenF(),cVal.blueF(),1);\r
1347                 //QVector4D color(d.x(),d.y(),d.z(),1);\r
1348 \r
1349                 FOR(s, segs)\r
1350                 {\r
1351                     int i1 = s%segs;\r
1352                     int i2 = (s+1)%segs;\r
1353                     int oi1 = (s+offset)%segs;\r
1354                     int oi2 = (s+1+offset)%segs;\r
1355                     o.vertices.append(oldCircle[oi1] + p1);\r
1356                     o.vertices.append(circle[i1] + p2);\r
1357                     o.vertices.append(circle[i2] + p2);\r
1358                     o.vertices.append(oldCircle[oi2] + p1);\r
1359                     o.normals.append(oldNormals[i1]);\r
1360                     o.normals.append(normals[i1]);\r
1361                     o.normals.append(normals[i2]);\r
1362                     o.normals.append(oldNormals[i2]);\r
1363                     o.colors.append(color);\r
1364                     o.colors.append(color);\r
1365                     o.colors.append(color);\r
1366                     o.colors.append(color);\r
1367                 }\r
1368             }\r
1369             oldP = p1;\r
1370             oldCircle = circle;\r
1371             oldNormals = normals;\r
1372             oldSpeed = speed;\r
1373         }\r
1374     }\r
1375     return o;\r
1376 }\r
1377 \r
1378 GLObject DrawStreamRibbon(vector<Streamline> streams, Dynamical *dynamical, vector<Obstacle> obstacles, float diff, float maxSpeed, int xInd, int yInd, int zInd)\r
1379 {\r
1380     GLObject o;\r
1381     o.objectType = "Dynamize,Surfaces,quads";\r
1382     o.style = "smooth";\r
1383     FOR(i, streams.size())\r
1384     {\r
1385         if(streams[i].length < 2) continue;\r
1386         float radius = diff*0.05;\r
1387         QVector3D p1, p2, q1, q2;\r
1388         QVector3D pos;\r
1389 \r
1390         // we generate the twin trajectory\r
1391         vector<fvec> twin;\r
1392         fvec sample(3);\r
1393         p1 = QVector3D(streams[i][0][xInd],streams[i][0][yInd],streams[i][0][zInd]);\r
1394         p2 = QVector3D(streams[i][1][xInd],streams[i][1][yInd],streams[i][1][zInd]);\r
1395         QVector3D d = (p2-p1).normalized();\r
1396         if(d.z() != 0)\r
1397         {\r
1398             float x=drand48(), y=drand48();\r
1399             float z=-(x*d.x() + y*d.y())/d.z();\r
1400             pos = QVector3D(x,y,z).normalized();\r
1401         }\r
1402         else\r
1403         {\r
1404             float z=drand48(), y=drand48();\r
1405             float x=-(y*d.y() + z*d.z())/d.x();\r
1406             pos = QVector3D(x,y,z).normalized();\r
1407         }\r
1408         q1 = pos + p1;\r
1409         sample[xInd] = pos.x()*radius + p1.x();\r
1410         sample[yInd] = pos.y()*radius + p1.y();\r
1411         sample[zInd] = pos.z()*radius + p1.z();\r
1412         twin.push_back(sample);\r
1413         FOR(j, streams[i].length-1)\r
1414         {\r
1415             fvec res = dynamical->Test(sample);\r
1416             if(dynamical->avoid)\r
1417             {\r
1418                 dynamical->avoid->SetObstacles(obstacles);\r
1419                 fvec newRes = dynamical->avoid->Avoid(sample, res);\r
1420                 res = newRes;\r
1421             }\r
1422             sample += res*dynamical->dT;\r
1423             fvec diff = sample - streams[i][j];\r
1424             float norm = sqrtf(diff*diff);\r
1425             if(norm > radius)\r
1426             {\r
1427                 sample = streams[i][j] + diff/norm*radius;\r
1428             }\r
1429             twin.push_back(sample);\r
1430         }\r
1431 \r
1432         float oldSpeed = 0;\r
1433         FOR(j, streams[i].length-1)\r
1434         {\r
1435             p1 = QVector3D(streams[i][j][xInd],streams[i][j][yInd],streams[i][j][zInd]);\r
1436             p2 = QVector3D(streams[i][j+1][xInd],streams[i][j+1][yInd],streams[i][j+1][zInd]);\r
1437             q1 = QVector3D(twin[j][xInd],twin[j][yInd],twin[j][zInd]);\r
1438             q2 = QVector3D(twin[j+1][xInd],twin[j+1][yInd],twin[j+1][zInd]);\r
1439             float speed = (p2-p1).length();\r
1440             speed += (q2-q1).length();\r
1441             speed /= 2;\r
1442 \r
1443             QColor cVal(Canvas::GetColorMapValue((speed+oldSpeed)*0.5/maxSpeed,2)); // jet color\r
1444             QVector4D color(cVal.redF(),cVal.greenF(),cVal.blueF(),1);\r
1445 \r
1446             //QVector4D color(d1.x(),d1.y(),d1.z(),1);\r
1447 \r
1448             o.vertices.append(p1);\r
1449             o.vertices.append(p2);\r
1450             o.vertices.append(q2);\r
1451             o.vertices.append(q1);\r
1452             o.colors.append(color);\r
1453             o.colors.append(color);\r
1454             o.colors.append(color);\r
1455             o.colors.append(color);\r
1456         }\r
1457     }\r
1458     return o;\r
1459 }\r
1460 \r
1461 GLObject DrawStreamLines(vector<Streamline> streams, int xInd, int yInd, int zInd)\r
1462 {\r
1463     GLObject o;\r
1464     o.objectType = "Dynamize,Lines";\r
1465     o.style = "";\r
1466 //    o.style = QString("fading:%1").arg(steps);\r
1467     FOR(i, streams.size())\r
1468     {\r
1469         QColor c = SampleColor[streams[i].cluster%(SampleColorCnt-1)+1];\r
1470         FOR(j, streams[i].length-1)\r
1471         {\r
1472             o.vertices.append(QVector3D(streams[i][j][xInd],streams[i][j][yInd],streams[i][j][zInd]));\r
1473             o.vertices.append(QVector3D(streams[i][j+1][xInd],streams[i][j+1][yInd],streams[i][j+1][zInd]));\r
1474             o.colors.append(QVector4D(c.redF(), c.greenF(), c.blueF(),1));\r
1475             o.colors.append(QVector4D(c.redF(), c.greenF(), c.blueF(),1));\r
1476         }\r
1477     }\r
1478     return o;\r
1479 }\r
1480 \r
1481 void Draw3DDynamical(GLWidget *glw, Dynamical *dynamical, int displayStyle)\r
1482 {\r
1483     if(!dynamical) return;\r
1484     int dim = glw->canvas->data->GetDimCount();\r
1485     if(dim != 3) return;\r
1486     float dT = dynamical->dT*2; // in 3d we want longer 'trails'\r
1487     int xInd = glw->canvas->xIndex;\r
1488     int yInd = glw->canvas->yIndex;\r
1489     int zInd = glw->canvas->zIndex;\r
1490     vector<fvec> samples = glw->canvas->data->GetSamples();\r
1491     vector< vector<fvec> > trajectories = glw->canvas->data->GetTrajectories(glw->canvas->trajectoryResampleType, glw->canvas->trajectoryResampleCount, glw->canvas->trajectoryCenterType, dT, true);\r
1492     vector<Obstacle> obstacles = glw->canvas->data->GetObstacles();\r
1493 \r
1494     fvec sample(dim,0);\r
1495     float minv=FLT_MAX, maxv=-FLT_MAX;\r
1496     FOR(i, samples.size())\r
1497     {\r
1498         FOR(d, dim)\r
1499         {\r
1500             minv = min(minv, samples[i][d]);\r
1501             maxv = max(maxv, samples[i][d]);\r
1502         }\r
1503     }\r
1504     float diff = maxv-minv;\r
1505     minv = minv - diff*0.5;\r
1506     maxv = maxv + diff*0.5;\r
1507     diff = maxv - minv;\r
1508 \r
1509     qDebug() << "computing field entropy";\r
1510     int gridSteps = 80;\r
1511     int hSteps = 20;\r
1512     int hRatio = gridSteps/hSteps;\r
1513     fvec H = ComputeDynamicalEntropy(dynamical, fvec(dim,minv), fvec(dim,maxv), gridSteps, hSteps);\r
1514     vector< pair<float,int> > hList(H.size());\r
1515     float hSum = 0;\r
1516     float hmin = FLT_MAX, hmax = -FLT_MAX;\r
1517     FOR(i, H.size())\r
1518     {\r
1519         hmin = min(hmin, H[i]);\r
1520         hmax = max(hmax, H[i]);\r
1521     }\r
1522 \r
1523     FOR(i, H.size())\r
1524     {\r
1525         hList[i] = make_pair((H[i]-hmin)/(hmax-hmin), i);\r
1526         hSum += hList[i].first;\r
1527     }\r
1528     //sort(hList.begin(), hList.end(), std::greater< pair<float,int> >());\r
1529     sort(hList.begin(), hList.end()); // smallest to largest\r
1530     reverse(hList.begin(), hList.end());\r
1531 \r
1532     /*\r
1533     FOR(i, hList.size())\r
1534     {\r
1535         qDebug() << "h values ["<< i <<"]:" << hList[i].first;\r
1536     }\r
1537     */\r
1538     // we generate seed points proportionally to the field entropy\r
1539     int seedlings = 80;\r
1540     vector<fvec> seeds(seedlings);\r
1541     float voxelSize = (maxv-minv)/(float)hSteps;\r
1542     FOR(s, seedlings)\r
1543     {\r
1544         float r = drand48()*hSum;\r
1545         float sum = 0;\r
1546         FOR(i, hList.size())\r
1547         {\r
1548             sum += hList[i].first;\r
1549             if(r < sum || i==hList.size()-1)\r
1550             {\r
1551                 int index = hList[i].second;\r
1552                 int x = index%hSteps;\r
1553                 int y = (index/hSteps)%hSteps;\r
1554                 int z = index/(hSteps*hSteps);\r
1555                 // we generate a random sample inside the selected voxel\r
1556                 sample[0] = drand48()*voxelSize + (x/(float)hSteps*(maxv-minv) + minv);\r
1557                 sample[1] = drand48()*voxelSize + (y/(float)hSteps*(maxv-minv) + minv);\r
1558                 sample[2] = drand48()*voxelSize + (z/(float)hSteps*(maxv-minv) + minv);\r
1559                 seeds[s] = sample;\r
1560                 break;\r
1561             }\r
1562         }\r
1563     }\r
1564 \r
1565     fvec mins(3,minv), maxes(3,maxv), origin=(maxes+mins)*0.5f;\r
1566     // we generate the trajectories\r
1567     int steps = 400;\r
1568     float maxSpeed = -FLT_MAX;\r
1569     vector<Streamline> streams(seeds.size());\r
1570     FOR(i, seeds.size())\r
1571     {\r
1572         sample = seeds[i];\r
1573         Streamline s;\r
1574         s.push_back(sample);\r
1575         FOR(j, steps)\r
1576         {\r
1577             fvec res = dynamical->Test(sample);\r
1578             if(dynamical->avoid)\r
1579             {\r
1580                 dynamical->avoid->SetObstacles(obstacles);\r
1581                 fvec newRes = dynamical->avoid->Avoid(sample, res);\r
1582                 res = newRes;\r
1583             }\r
1584             float speed = sqrtf((res*dT)*(res*dT));\r
1585             if(speed > maxSpeed) maxSpeed = speed;\r
1586             sample += res*dT;\r
1587             if(speed < 1e-4) break;\r
1588             if(sqrtf((sample - origin)*(sample - origin)) > diff*0.7) break;\r
1589             s.push_back(sample);\r
1590         }\r
1591         s.cluster = i;\r
1592         streams[i] = s;\r
1593     }\r
1594 \r
1595     GLObject o;\r
1596     switch(displayStyle)\r
1597     {\r
1598     case 0: // entropy field\r
1599         o = DrawEntropyField(H, minv, maxv, hSteps);\r
1600         break;\r
1601     case 1: // tubes\r
1602         o = DrawStreamTubes(streams, diff, maxSpeed, xInd, yInd, zInd);\r
1603         break;\r
1604     case 2: // ribbons\r
1605         o = DrawStreamRibbon(streams, dynamical, obstacles, diff, maxSpeed, xInd, yInd, zInd);\r
1606         break;\r
1607     case 3: // Animation\r
1608         o = DrawStreamLines(streams, xInd, yInd, zInd);\r
1609         break;\r
1610     }\r
1611 \r
1612     // we replace the old vector field (if there is one) with the new one\r
1613     glw->mutex->lock();\r
1614     int oInd = -1;\r
1615     FOR(i, glw->objects.size())\r
1616     {\r
1617         if(glw->objects[i].objectType.contains("Dynamize"))\r
1618         {\r
1619             oInd = i;\r
1620             break;\r
1621         }\r
1622     }\r
1623     if(oInd != -1) glw->objects[oInd] = o;\r
1624     else glw->objects.push_back(o);\r
1625     glw->mutex->unlock();\r
1626 \r
1627     qDebug() << "done.";\r
1628     return;\r
1629 \r
1630 /*\r
1631     int steps = 400;\r
1632     int gridSize = 10;\r
1633     int count = gridSize*gridSize*gridSize;\r
1634 //    int count = 1024;\r
1635 \r
1636     // we generate a bunch of trajectories\r
1637     qDebug() << "generating streamlines";\r
1638     vector<Streamline> streams(count);\r
1639     FOR(i, count)\r
1640     {\r
1641         int index = i;\r
1642         int x = index % gridSize;\r
1643         int y = (index / gridSize) % gridSize;\r
1644         int z = index / (gridSize*gridSize);\r
1645         sample[xInd] = (x/(float)gridSize)*diff + minv;\r
1646         sample[yInd] = (y/(float)gridSize)*diff + minv;\r
1647         sample[zInd] = (z/(float)gridSize)*diff + minv;\r
1648         //FOR(d, dim) sample[d] = drand48()*diff+ minv;\r
1649 \r
1650         Streamline s;\r
1651         s.push_back(sample);\r
1652         FOR(j, steps)\r
1653         {\r
1654             fvec res = dynamical->Test(sample);\r
1655             if(dynamical->avoid)\r
1656             {\r
1657                 dynamical->avoid->SetObstacles(obstacles);\r
1658                 fvec newRes = dynamical->avoid->Avoid(sample, res);\r
1659                 res = newRes;\r
1660             }\r
1661             sample += res*dT;\r
1662             if(res*res < 1e-5) break;\r
1663             s.push_back(sample);\r
1664         }\r
1665         streams[i] = s;\r
1666     }\r
1667 \r
1668     qDebug() << "clustering streamlines";\r
1669 \r
1670     // now we "cluster" them together\r
1671     int nbClusters = 12;\r
1672     int maxIterations = 40;\r
1673     ivec clusters = ClusterStreams(streams, nbClusters, maxIterations);\r
1674     ivec counts(nbClusters,0);\r
1675     FOR(i, nbClusters) seeds[i].resize(dim);\r
1676     FOR(i, streams.size())\r
1677     {\r
1678         streams[i].cluster = clusters[i];\r
1679         seeds[clusters[i]] += streams[i].front();\r
1680         counts[clusters[i]]++;\r
1681     }\r
1682     FOR(i, nbClusters) seeds[i] /= counts[i];\r
1683 \r
1684 \r
1685     // we add explicitly the training trajectories\r
1686     FOR(i, trajectories.size())\r
1687     {\r
1688         seeds.push_back(trajectories[i][0]);\r
1689     }\r
1690 \r
1691     // and now we regenerate our streamlines from the seeds\r
1692     qDebug() << "generating clustered seeds";\r
1693     vector<Streamline> clusteredStreams(seeds.size());\r
1694 \r
1695     qDebug() << "generating GLObjects";\r
1696 \r
1697     int streamStyle = 1;\r
1698     //GLObject o;\r
1699 \r
1700     switch(streamStyle)\r
1701     {\r
1702     case 0: // vector field\r
1703     {\r
1704         o.objectType = "Dynamize,Lines";\r
1705         o.style = "";\r
1706     //    o.style = QString("fading:%1").arg(steps);\r
1707         FOR(i, streams.size())\r
1708         {\r
1709             QColor c = SampleColor[streams[i].cluster%(SampleColorCnt-1)+1];\r
1710             FOR(j, streams[i].length-1)\r
1711             {\r
1712                 o.vertices.append(QVector3D(streams[i][j][xInd],streams[i][j][yInd],streams[i][j][zInd]));\r
1713                 o.vertices.append(QVector3D(streams[i][j+1][xInd],streams[i][j+1][yInd],streams[i][j+1][zInd]));\r
1714                 o.colors.append(QVector4D(c.redF(), c.greenF(), c.blueF(),1));\r
1715                 o.colors.append(QVector4D(c.redF(), c.greenF(), c.blueF(),1));\r
1716             }\r
1717         }\r
1718     }\r
1719         break;\r
1720     case 1: // tubes\r
1721     {\r
1722         o = DrawStreamTubes(streams);\r
1723     }\r
1724         break;\r
1725     case 2: // stripes\r
1726     {\r
1727         o = DrawStreamRibbons(streams);\r
1728     }\r
1729         break;\r
1730     }\r
1731 \r
1732     // we replace the old vector field (if there is one) with the new one\r
1733     glw->mutex->lock();\r
1734     int oIndex = -1;\r
1735     FOR(i, glw->objects.size())\r
1736     {\r
1737         if(glw->objects[i].objectType.contains("Dynamize"))\r
1738         {\r
1739             oIndex = i;\r
1740             break;\r
1741         }\r
1742     }\r
1743     if(oIndex != -1) glw->objects[oIndex] = o;\r
1744     else glw->objects.push_back(o);\r
1745     glw->mutex->unlock();\r
1746     */\r
1747 }\r
1748 \r
1749 void Draw3DMaximizer(GLWidget *glw, Maximizer *maximizer){}\r
1750 void Draw3DProjector(GLWidget *glw, Projector *projector){}\r
1751 void Draw3DReinforcement(GLWidget *glw, Reinforcement *reinforcement){}\r
1752 \r
1753 QLabel *label = 0;\r
1754 void Draw2DDynamical(Canvas *canvas, Dynamical *dynamical)\r
1755 {\r
1756     int w = canvas->width();\r
1757     int h = canvas->height();\r
1758     // we start by generating random noise\r
1759     QImage pixels(w,h,QImage::Format_RGB32);\r
1760     FOR(i, w*h)\r
1761     {\r
1762         int x = i%w;\r
1763         int y = i/w;\r
1764         float value = max(0.f,min(1.f,RandN(0.5f,0.5f)));\r
1765         pixels.setPixel(x,y,qRgb(value*255,value*255,value*255));\r
1766     }\r
1767 \r
1768     // now we "process" the noise by iteratively applying the DS to it\r
1769     QPainter painter(&pixels);\r
1770     painter.setRenderHint(QPainter::Antialiasing);\r
1771     int iterations = 4;\r
1772     float dT = 0.004;\r
1773     vector<Obstacle> obstacles = canvas->data->GetObstacles();\r
1774     qDebug() << "processing noise";\r
1775     FOR(i, iterations)\r
1776     {\r
1777         qDebug() << "iteration" << i;\r
1778         FOR(y, h)\r
1779         {\r
1780             FOR(x, w)\r
1781             {\r
1782                 QPoint point(x,y);\r
1783                 QRgb val = pixels.pixel(point);\r
1784                 fvec sample = canvas->fromCanvas(x,y);\r
1785                 fvec res = dynamical->Test(sample);\r
1786                 if(dynamical->avoid)\r
1787                 {\r
1788                     dynamical->avoid->SetObstacles(obstacles);\r
1789                     fvec newRes = dynamical->avoid->Avoid(sample, res);\r
1790                     res = newRes;\r
1791                 }\r
1792                 sample += res*dT;\r
1793                 QPointF point2 = canvas->toCanvasCoords(sample);\r
1794                 painter.setPen(QColor(val));\r
1795                 painter.drawLine(point, point2);\r
1796                 //if(point.x() < 0 || point.x() >= w || point.y() < 0 || point.y() >= h) continue;\r
1797                 //pixels.setPixel(point.x(), point.y(), val);\r
1798             }\r
1799         }\r
1800     }\r
1801 \r
1802     QPixmap pixmap = QPixmap::fromImage(pixels);\r
1803 //    QPixmap pixmap = QPixmap::fromImage(pixels).scaled(w*2, h*2,Qt::IgnoreAspectRatio, Qt::SmoothTransformation);\r
1804 \r
1805     if(!label)\r
1806     {\r
1807         label = new QLabel();\r
1808         label->setScaledContents(true);\r
1809     }\r
1810     label->setPixmap(pixmap);\r
1811     label->show();\r
1812 }\r