EMAN2
glutil.cpp
Go to the documentation of this file.
1/*
2 * Author: David Woolford, 11/06/2007 (woolford@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#ifdef USE_OPENGL
33
34#ifdef _WIN32
35 #include <windows.h>
36#endif
37
38#ifndef GL_GLEXT_PROTOTYPES
39 #define GL_GLEXT_PROTOTYPES
40#endif //GL_GLEXT_PROTOTYPES
41
42#include <utility>
43
44#include "glutil.h"
45#include "marchingcubes.h"
46
47#ifdef __APPLE__
48 #include "OpenGL/gl.h"
49 #include "OpenGL/glu.h"
50 #include "OpenGL/glext.h"
51#else // WIN32, LINUX
52 #include "GL/gl.h"
53 #include "GL/glu.h"
54 #include "GL/glext.h"
55#endif //__APPLE__
56
57using namespace EMAN;
58
59// By depfault we need to first bind data to the GPU
60GLuint GLUtil::buffer[2] = {0, 0};
61
62unsigned int GLUtil::gen_glu_mipmaps(const EMData* const emdata)
63{
64 if (emdata->get_data() == 0) {
65 throw NullPointerException("Error, attempt to create an OpenGL mipmap "
66 "without internally stored data");
67 }
68
69 ENTERFUNC;
70
71 unsigned int tex_name;
72 glGenTextures(1, &tex_name);
73
74 if (emdata->ny == 1 && emdata->nz == 1) {
75 glBindTexture(GL_TEXTURE_1D, tex_name);
76 gluBuild1DMipmaps(GL_TEXTURE_1D, GL_LUMINANCE, emdata->nx, GL_LUMINANCE,
77 GL_FLOAT, (void*)(emdata->get_data()));
78 } else if (emdata->nz == 1) {
79 glBindTexture(GL_TEXTURE_2D, tex_name);
80 gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, emdata->nx, emdata->ny,
81 GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
82 }
83 else {
84#ifdef _WIN32
85 // There is no gluBuild3DMipmaps() function in glu.h for VS2003 and VS2005
86
87 printf("3D OpenGL mipmaps are not available on this platform.\n");
88#else
89 glBindTexture(GL_TEXTURE_3D, tex_name);
90 gluBuild3DMipmaps(GL_TEXTURE_3D, GL_LUMINANCE, emdata->nx, emdata->ny,
91 emdata->nz, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
92#endif //_WIN32
93 }
94
95 EXITFUNC;
96
97 return tex_name;
98}
99
100unsigned int GLUtil::gen_gl_texture(const EMData* const emdata, GLenum format)
101{
102 if (emdata->get_data() == 0) {
103 throw NullPointerException("Error, attempt to create an OpenGL texture"
104 " without internally stored data");
105 }
106
107 ENTERFUNC;
108
109 unsigned int tex_name;
110 glGenTextures(1, &tex_name);
111
112 if (emdata->ny == 1 && emdata->nz == 1) {
113 glBindTexture(GL_TEXTURE_1D, tex_name);
114 glTexImage1D(GL_TEXTURE_1D, 0, format, emdata->nx, 0, format, GL_FLOAT,
115 (void*)(emdata->get_data()));
116 } else if (emdata->nz == 1) {
117 glBindTexture(GL_TEXTURE_2D, tex_name);
118 glTexImage2D(GL_TEXTURE_2D, 0, format, emdata->nx, emdata->ny, 0, format,
119 GL_FLOAT, (void*)(emdata->get_data()));
120 }
121 else {
122 glBindTexture(GL_TEXTURE_3D, tex_name);
123#ifdef _WIN32
124 PFNGLTEXIMAGE3DPROC glTexImage3D;
125#endif
126 glTexImage3D(GL_TEXTURE_3D, 0, format,
127 emdata->nx, emdata->ny, emdata->nz, 0, format, GL_FLOAT,
128 (void*)(emdata->get_data()));
129 }
130
131 EXITFUNC;
132
133 return tex_name;
134}
135
136unsigned int GLUtil::render_amp8_gl_texture(EMData* emdata,
137 int x0, int y0, int ixsize, int iysize, int bpl, float scale,
138 int mingray, int maxgray, float render_min, float render_max,
139 float gamma, int flags)
140{
141 if (emdata==NULL) return 9999999;
142 string pixels = render_amp8(emdata, x0, y0, ixsize,iysize, bpl, scale,
143 mingray, maxgray, render_min, render_max, gamma, flags);
144
145 unsigned int tex_name;
146 glGenTextures(1, &tex_name);
147
148 glBindTexture(GL_TEXTURE_2D, tex_name);
149 glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, ixsize, iysize , 0,
150 GL_LUMINANCE, GL_UNSIGNED_BYTE, pixels.c_str());
151
152 return tex_name;
153}
154
155// undef GL_GLEXT_PROTOTYPES
156
157#ifdef GL_GLEXT_PROTOTYPES
158#undef GL_GLEXT_PROTOTYPES
159#endif
160
161int GLUtil::nearest_projected_points(const vector<float>& model_matrix,
162 const vector<float>& proj_matrix, const vector<int>& view_matrix,
163 const vector<Vec3f>& points, const float mouse_x, const float mouse_y,
164 const float& nearness)
165{
166 double proj[16];
167 double model[16];
168 int view[4];
169
170 copy(proj_matrix.begin(), proj_matrix.end(), proj);
171 copy(model_matrix.begin(), model_matrix.end(), model);
172 copy(view_matrix.begin(), view_matrix.end(), view);
173
174 vector<Vec3f> unproj_points;
175 double x,y,z;
176 double r,s,t;
177
178 for (vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
179 r = (double) (*it)[0];
180 s = (double) (*it)[1];
181 t = (double) (*it)[2];
182
183 gluProject(r,s,t, model, proj, view, &x, &y, &z);
184 unproj_points.push_back(Vec3f(x, y, z));
185 }
186
187 vector<int> intersections;
188
189 float n_squared = nearness * nearness;
190
191 for(unsigned int i = 0; i < unproj_points.size(); ++i) {
192 Vec3f& v = unproj_points[i];
193 float dx = v[0] - mouse_x;
194 dx *= dx;
195 float dy = v[1] - mouse_y;
196 dy *= dy;
197
198 if ((dx+dy) <= n_squared) intersections.push_back((int)i);
199 }
200
201 int intersection = -1;
202 float near_z = 0;
203
204 for(vector<int>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
205 if (intersection == -1 || unproj_points[*it][2] < near_z) {
206 intersection = *it;
207 near_z = unproj_points[*it][2];
208 }
209 }
210
211 return intersection;
212}
213
214void GLUtil::colored_rectangle(const vector<float>& data,
215 const float& alpha, const bool center_point)
216{
217 glBegin(GL_LINE_LOOP);
218 glColor4f(data[0], data[1], data[2], alpha);
219 glVertex2i(data[3], data[4]);
220 glVertex2i(data[5], data[4]);
221 glVertex2i(data[5], data[6]);
222 glVertex2i(data[3], data[6]);
223 glEnd();
224
225 if (center_point) {
226 glBegin(GL_POINTS);
227 glVertex2f( (data[3]+data[5])/2, (data[4]+data[6])/2);
228 glEnd();
229 }
230}
231
232void GLUtil::mx_bbox(const vector<float>& data,
233 const vector<float>& text_color, const vector<float>& bg_color)
234{
235 glDisable(GL_TEXTURE_2D);
236 glBegin(GL_QUADS);
237
238 if (bg_color.size() == 4) {
239 glColor4f(bg_color[0], bg_color[1], bg_color[2], bg_color[3]);
240 }
241 else {
242 glColor3f(bg_color[0], bg_color[1], bg_color[2]);
243 }
244
245 glVertex3f(data[0]-1, data[1]-1, -.1f);
246 glVertex3f(data[3]+1, data[1]-1, -.1f);
247 glVertex3f(data[3]+1, data[4]+1, -.1f);
248 glVertex3f(data[0]-1, data[4]+1, -.1f);
249 glEnd();
250
251// glEnable(GL_TEXTURE_2D);
252
253 if (text_color.size() == 4) {
254 glColor4f(text_color[0], text_color[1], text_color[2], text_color[3]);
255 }
256 else {
257 glColor3f(text_color[0], text_color[1], text_color[2]);
258 }
259}
260
261EMBytes GLUtil::render_amp8(EMData* emdata, int x0, int y0, int ixsize,
262 int iysize, int bpl, float scale, int min_gray, int max_gray,
263 float render_min, float render_max, float gamma, int flags)
264{
265 ENTERFUNC;
266
267// printf("%f\t%f\t",(float)emdata->get_attr("sigma"),(float)emdata->get_attr("mean"));
268// printf("%d %d %d %d %d %f %d %d %f %f %f %d\n",x0,y0,ixsize,iysize,bpl,
269// scale,min_gray,max_gray,render_min,render_max,gamma,flags);
270
271 if (emdata==NULL) return EMBytes();
272 bool invert = (min_gray > max_gray);
273 int mingray, maxgray;
274
275 if (invert) {
276 mingray = max_gray;
277 maxgray = min_gray;
278 }
279 else {
280 mingray = min_gray;
281 maxgray = max_gray;
282 }
283
284 int asrgb;
285 int hist = (flags & 2)/2;
286 int invy = (flags & 4)?1:0;
287
288 int nx = emdata->nx;
289 int ny = emdata->ny;
290// int nz = emdata->nz;
291 int nxy = emdata->nx * emdata->ny;
292
293 if (emdata->get_ndim() > 2) {
294 throw ImageDimensionException("1D/2D only");
295 }
296
297 if (emdata->is_complex()) {
298 emdata->ri2ap();
299 }
300
301 if (render_max <= render_min) {
302 render_max = render_min + 0.01f;
303 }
304
305 if (gamma <= 0) gamma = 1.0;
306
307 // Calculating a full floating point gamma for
308 // each piGLUtil::xel in the image slows rendering unacceptably
309 // however, applying a gamma-mapping to an 8 bit colorspace
310 // has unacceptable coarse accuracy. So, we oversample the 8 bit colorspace
311 // as a 12 bit colorspace and apply the gamma mapping to that
312 // This should produce good accuracy for gamma values
313 // larger than 0.5 (and a high upper limit)
314 static int smg0 = 0, smg1 = 0; // while this destroys threadsafety in the rendering process
315 static float sgam=0; // it is necessary for speed when rendering large numbers of small images
316 static unsigned char gammamap[4096];
317
318 if (gamma != 1.0 && (smg0 != mingray || smg1 != maxgray || sgam != gamma)) {
319 for (int i=0; i<4096; i++) {
320 if (mingray<maxgray) {
321 gammamap[i] = (unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
322 }
323 else {
324 gammamap[4095-i] = (unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
325 }
326 }
327 }
328 smg0 = mingray; // so we don't recompute the map unless something changes
329 smg1 = maxgray;
330 sgam = gamma;
331
332 if (flags & 8) asrgb = 4;
333 else if (flags & 1) asrgb = 3;
334 else asrgb = 1;
335
336 EMBytes ret=EMBytes();
337// ret.resize(iysize*bpl);
338
339 ret.assign(iysize*bpl + hist*1024, char(invert ? maxgray : mingray));
340
341 unsigned char *data = (unsigned char *)ret.data();
342 unsigned int *histd = (unsigned int *)(data + iysize*bpl);
343
344 if (hist) {
345 for (int i=0; i<256; i++) histd[i]=0;
346 }
347
348 float rm = render_min;
349 float inv_scale = 1.0f / scale;
350 int ysize = iysize;
351 int xsize = ixsize;
352
353 int ymin = 0;
354
355 if (iysize * inv_scale > ny) {
356 ymin = (int) (iysize - ny / inv_scale);
357 }
358
359 float gs = (maxgray - mingray) / (render_max - render_min);
360 float gs2 = 4095.999f / (render_max - render_min);
361// float gs2 = 1.0 / (render_max - render_min);
362
363 if (render_max < render_min) {
364 gs = 0;
365 rm = FLT_MAX;
366 }
367
368 int dsx = -1;
369 int dsy = 0;
370 int remx = 0;
371 int remy = 0;
372 const int scale_n = 100000;
373
374 int addi = 0;
375 int addr = 0;
376
377 if (inv_scale == floor(inv_scale)) {
378 dsx = (int) inv_scale;
379 dsy = (int) (inv_scale * nx);
380 }
381 else {
382 addi = (int) floor(inv_scale);
383 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
384 }
385
386 int xmin = 0;
387
388 if (x0 < 0) {
389 xmin = (int) (-x0 / inv_scale);
390 xsize -= (int) floor(x0 / inv_scale);
391 x0 = 0;
392 }
393
394 if ((xsize - xmin) * inv_scale > (nx - x0)) {
395 xsize = (int) ((nx - x0) / inv_scale + xmin);
396 }
397
398 int ymax = ysize - 1;
399
400 if (y0 < 0) {
401 ymax = (int) (ysize + y0 / inv_scale - 1);
402 ymin += (int) floor(y0 / inv_scale);
403 y0 = 0;
404 }
405
406 if (xmin < 0) xmin = 0;
407 if (ymin < 0) ymin = 0;
408 if (xsize > ixsize) xsize = ixsize;
409 if (ymax > iysize) ymax = iysize;
410
411 int lmax = nx * ny - 1;
412
413 int mid=nx*ny/2;
414 float * image_data = emdata->get_data();
415
417
418 const int rangemax = 4082;
419 int gaussianwide = 5000;
420 unsigned int* grayhe = NULL;
421 float* gaussianpdf = NULL;
422 float* gaussiancdf = NULL;
423 int* gaussianlookup = NULL;
424
425 unsigned int* graypdftwo = NULL;
426
427 bool binflag = 0;
428 int binlocation = -1;
429
430 if (flags & 32) {
431 int graypdf[rangemax] = {0}; // 256
432 int graycdf[rangemax-2] = {0}; // 254
433
434 // unsigned int grayhe[rangemax-2]={0};//#number=254
435
436 graypdftwo = new unsigned int[maxgray-mingray];
437 grayhe = new unsigned int[rangemax-2];
438 // unsigned char graylookup[(int)(render_max-render_min)];//render_max-render_min
439
440 for (int i=0; i<(int)(rangemax-2); i++) {
441 grayhe[i] = 0; // Initialize all elements to zero.
442 }
443
444 for (int i=0; i<(maxgray-mingray); i++) { // 0~253
445 graypdftwo[i] = 0;
446 }
447
448 if (dsx != -1) {
449 int l = x0 + y0 * nx;
450
451 for (int j = ymax; j >= ymin; j--) {
452 int br = l;
453 for (int i = xmin; i < xsize; i++) {
454 if (l > lmax) {
455 break;
456 }
457
458 int k = 0;
459 int p;
460 float t;
461
462 if (dsx == 1) t=image_data[l];
463 else { // This block does local pixel averaging for nicer reduced views
464 t=0;
465
466 if ((l+dsx+dsy) > lmax) {
467 break;
468 }
469
470 for (int iii=0; iii<dsx; iii++) {
471 for (int jjj=0; jjj<dsy; jjj+=nx) {
472 t += image_data[l+iii+jjj];
473 }
474 }
475
476 t /= dsx*(dsy/nx);
477 }
478
479 if (t <= rm) graypdf[0]++;
480 else if (t >= render_max) graypdf[rangemax-1]++;
481 else {
482 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
483 graypdftwo[(unsigned char) (gs * (t - render_min))]++;
484 }
485 l += dsx;
486 }
487
488 l = br + dsy;
489 }
490 }
491 else {
492 remy = 10;
493 int l = x0 + y0 * nx;
494
495 for (int j = ymax; j >= ymin; j--) {
496 int addj = addi;
497
498 // There seems to be some overflow issue happening
499 // where the statement if (l > lmax) break (below) doesn't work
500 // because the loop that iterates jjj can inadvertantly go out of bounds
501 if ((l + addi*nx) >= nxy) {
502 addj = (nxy-l)/nx;
503
504 if (addj <= 0) continue;
505 }
506
507 int br = l;
508 remx = 10;
509
510 for (int i = xmin; i < xsize; i++) {
511 if (l > lmax) break;
512 int p = 0;
513 float t;
514 if (addi <= 1) t = image_data[l];
515 else { // This block does local pixel averaging for nicer reduced views
516 t = 0;
517
518 for (int jjj=0; jjj<addj; jjj++) {
519 for (int iii=0; iii<addi; iii++) {
520 t += image_data[l+iii+jjj*nx];
521 }
522 }
523
524 t /= addi*addi;
525 }
526
528 if (t <= rm) graypdf[0]++;
529 else if (t >= render_max) graypdf[rangemax-1]++;
530 else {
531 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
532 }
533
534// data[i * asrgb + j * bpl] = p;
535// if (hist) histd[p]++;
536 l += addi;
537 remx += addr;
538
539 if (remx > scale_n) {
540 remx -= scale_n;
541 l++;
542 }
543 }
544
545 l = br + addi * nx;
546 remy += addr;
547
548 if (remy > scale_n) {
549 remy -= scale_n;
550 l += nx;
551 }
552 }
553 }
554
555 for (int i=0; i<(rangemax-2); i++) { // 0~253
556 for (int j=0;j<(i+1);j++) {
557 graycdf[i]=graycdf[i]+graypdf[j+1];
558 }
559 }
560
561 // graypdftwo binflag
562 binflag = 0;
563
564 for (int i=1; i<(maxgray-mingray); i++) { // 0~253
565 if (((float)graypdftwo[i]/graycdf[rangemax-3])>0.2) {
566 binflag = 1;
567 binlocation = i;
568
569 break;
570 }
571 }
572
573 if (binflag == 1) {
574 for (int i=(binlocation*16+1); i<((binlocation+1)*16); i++) {
575 graypdf[i] = 0;
576 // graypdf[i]=(graycdf[rangemax-3]-graypdftwo[binlocation])/(maxgray-mingray);
577 }
578
579 for (int i=0; i<(rangemax-2); i++) { // 0~253
580 graycdf[i] = 0;
581 }
582
583 for (int i=0; i<(rangemax-2); i++) { // 0~253
584 for (int j=0;j<(i+1);j++) {
585 graycdf[i] = graycdf[i]+graypdf[j+1];
586 }
587 }
588 }
589
590 // start gaussian matching
591 float mean = abs(rangemax-2)/2;
592 float standdv = abs(mean)/3;
593
594 gaussianpdf = new float[rangemax-2];
595 gaussiancdf = new float[rangemax-2];
596 gaussianlookup = new int[rangemax-2];
597
598 for (int i=0; i<(rangemax-2); i++) {
599 gaussianpdf[i]=exp(-(i-mean)*(i-mean)/(2*standdv*standdv))/sqrt(standdv * standdv * 2 * M_PI);
600
601 if (i != 0) {
602 gaussiancdf[i] = gaussiancdf[i-1]+gaussianpdf[i];
603 }
604 else {
605 gaussiancdf[i] = gaussianpdf[i];
606 }
607 }
608
609 for (int i=0; i<(rangemax-2); i++) {
610 gaussiancdf[i] = graycdf[rangemax-3]*gaussiancdf[i]/gaussiancdf[rangemax-3];
611 }
612
613 for (int i=0; i<(rangemax-2); i++) {
614 for (int j=0; j<(rangemax-2); j++) {
615 if (graycdf[i] <= gaussiancdf[j]) {
616 gaussianlookup[i] = j;
617 break;
618 }
619 }
620 }
621
622 for (int i=0; i<(rangemax-2); i++) {
623 grayhe[i] = floor(0.5+(((double)(rangemax-3)*graycdf[i])/graycdf[rangemax-3]));
624 }
625 }
627
628 if (emdata->is_complex()) {
629 if (dsx != -1) {
630 int l = y0 * nx;
631
632 for (int j = ymax; j >= ymin; j--) {
633 int ll = x0;
634
635 for (int i = xmin; i < xsize; i++) {
636 if (l + ll > lmax || ll >= nx - 2) break;
637
638 int k = 0;
639 unsigned char p;
640
641 if (ll >= nx / 2) {
642 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
643 else k = 2 * (ll - nx / 2) + l + 2 + nx;
644 }
645 else k = nx * ny - (l + 2 * ll) - 2;
646
647 if (k >= mid) k -= mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
648 else k += mid;
649
650 float t = image_data[k];
651 int ph;
652
653 // in color mode
654 if (flags & 16 && asrgb>2) {
655// if (l >= (ny - inv_scale) * nx) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767;
656 if (ll >= nx / 2) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767;
657 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767;
658 }
659
660 if (t <= rm) p = mingray;
661 else if (t >= render_max) p = maxgray;
662 else if (gamma != 1.0) {
663 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
664 p = gammamap[k]; // apply gamma using precomputed gamma map
665// p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
666// p += mingray;
667// k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
668// k += mingray;
669 }
670 else {
671 p = (unsigned char) (gs * (t - render_min));
672 p += mingray;
673 }
674
675 // color rendering
676 if (flags & 16 && asrgb>2) {
677 if (ph<256) {
678 data[i * asrgb + j * bpl] = p*(255-ph)/256;
679 data[i * asrgb + j * bpl+1] = p*ph/256;
680 data[i * asrgb + j * bpl+2] = 0;
681 }
682 else if (ph<512) {
683 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
684 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
685 data[i * asrgb + j * bpl] = 0;
686 }
687 else {
688 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
689 data[i * asrgb + j * bpl] = p*(ph-512)/256;
690 data[i * asrgb + j * bpl+1] = 0;
691 }
692 }
693 else data[i * asrgb + j * bpl] = p;
694 if (hist) histd[p]++;
695 ll += dsx;
696 }
697
698 l += dsy;
699 }
700 }
701 else {
702 remy = 10;
703 int l = y0 * nx;
704
705 for (int j = ymax; j >= ymin; j--) {
706 int br = l;
707 remx = 10;
708 int ll = x0;
709
710 for (int i = xmin; i < xsize - 1; i++) {
711 if (l + ll > lmax || ll >= nx - 2) {
712 break;
713 }
714
715 int k = 0;
716 unsigned char p;
717
718 if (ll >= nx / 2) {
719 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
720 else k = 2 * (ll - nx / 2) + l + 2 + nx;
721 }
722 else k = nx * ny - (l + 2 * ll) - 2;
723
724 if (k >= mid) k -= mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
725 else k += mid;
726
727 float t = image_data[k];
728 // in color mode
729 int ph;
730
731 if (flags & 16 && asrgb>2) {
732 if (l >= (ny * nx - nx)) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767;
733 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767;
734 }
735
736 if (t <= rm)
737 p = mingray;
738 else if (t >= render_max) {
739 p = maxgray;
740 }
741 else if (gamma != 1.0) {
742 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
743 p = gammamap[k]; // apply gamma using precomputed gamma map
744
745// p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
746// p += mingray;
747// k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
748// k += mingray;
749 }
750 else {
751 p = (unsigned char) (gs * (t - render_min));
752 p += mingray;
753 }
754
755 if (flags & 16 && asrgb>2) {
756 if (ph<256) {
757 data[i * asrgb + j * bpl] = p*(255-ph)/256;
758 data[i * asrgb + j * bpl+1] = p*ph/256;
759 data[i * asrgb + j * bpl+2] = 0;
760 }
761 else if (ph<512) {
762 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
763 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
764 data[i * asrgb + j * bpl] = 0;
765 }
766 else {
767 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
768 data[i * asrgb + j * bpl] = p*(ph-512)/256;
769 data[i * asrgb + j * bpl+1] = 0;
770 }
771 }
772 else data[i * asrgb + j * bpl] = p;
773
774 if (hist) histd[p]++;
775
776 ll += addi;
777 remx += addr;
778
779 if (remx > scale_n) {
780 remx -= scale_n;
781 ll++;
782 }
783 }
784
785 l = br + addi * nx;
786 remy += addr;
787
788 if (remy > scale_n) {
789 remy -= scale_n;
790 l += nx;
791 }
792 }
793 }
794 }
795 else {
796 if (dsx != -1) {
797 int l = x0 + y0 * nx;
798
799 for (int j = ymax; j >= ymin; j--) {
800 int br = l;
801
802 for (int i = xmin; i < xsize; i++) {
803 if (l > lmax) {
804 break;
805 }
806
807 int k = 0;
808 unsigned char p;
809 float t;
810
811 if (dsx == 1) t=image_data[l];
812 else { // This block does local pixel averaging for nicer reduced views
813 t=0;
814
815 if ((l+dsx+dsy) > lmax) {
816 break;
817 }
818
819 for (int iii=0; iii<dsx; iii++) {
820 for (int jjj=0; jjj<dsy; jjj+=nx) {
821 t += image_data[l+iii+jjj];
822 }
823 }
824
825 t /= dsx*(dsy/nx);
826 }
827
828 if (t <= rm) p = mingray;
829 else if (t >= render_max) p = maxgray;
830 else if (gamma != 1.0) {
831 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
832 p = gammamap[k];// apply gamma using precomputed gamma map
833// k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
834// k += mingray;
835// k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
836// k += mingray;
837 }
838 else {
839 if (flags & 32) {
840 //p = graylookup[(int)(t - render_min)];
841 //graylookup[i]=(unsigned char)((maxgray-mingray-2)*grayhe[(int)(i*(rangemax-3)/(render_max-render_min))]/(rangemax-3)+1);
842 if (flags & 64)
843 {
844 p=(unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
845 }
846 else {
847 p=(unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);}
848 //p=(unsigned char)gaussianlookup[(int)(ceil)((t - render_min)*(rangemax-3)/(render_max-render_min))]+1;
849 }
850 else
851 {
852 p=(unsigned char) (gs * (t - render_min));
853 }
854
855 p += mingray;
856 }
857
858 if (invert) {
859 p = mingray + maxgray - p;
860 }
861
862 data[i * asrgb + j * bpl] = p;
863
864 if (hist) histd[p]++;
865
866 l += dsx;
867 }
868
869 l = br + dsy;
870 }
871 }
872 else {
873 remy = 10;
874 int l = x0 + y0 * nx;
875
876 for (int j = ymax; j >= ymin; j--) {
877 int addj = addi;
878
879 // There seems to be some overflow issue happening
880 // where the statement if (l > lmax) break (below) doesn't work
881 // because the loop that iterates jjj can inadvertantly go out of bounds
882 if ((l + addi*nx) >= nxy) {
883 addj = (nxy-l)/nx;
884
885 if (addj <= 0) continue;
886 }
887
888 int br = l;
889 remx = 10;
890
891 for (int i = xmin; i < xsize; i++) {
892 if (l > lmax) break;
893 int k = 0;
894 unsigned char p;
895 float t;
896
897 if (addi <= 1) t = image_data[l];
898 else { // This block does local pixel averaging for nicer reduced views
899 t=0;
900 for (int jjj=0; jjj<addj; jjj++) {
901 for (int iii=0; iii<addi; iii++) {
902 t += image_data[l+iii+jjj*nx];
903 }
904 }
905
906 t /= addi*addi;
907 }
908
909 if (t <= rm) p = mingray;
910 else if (t >= render_max) p = maxgray;
911 else if (gamma != 1.0) {
912 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
913 p = gammamap[k]; // apply gamma using precomputed gamma map
914// k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
915// k += mingray;
916// k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
917// k += mingray;
918 }
919 else {
920 if (flags & 32) {
921 // p = graylookup[(int)(t - render_min)];
922
923 if (flags & 64) {
924 p = (unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
925 }
926 else {
927 p = (unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);
928 }
929
930 // p = (unsigned char)gaussianlookup[(int)(ceil)((t - render_min)*(rangemax-3)/(render_max-render_min))]+1;
931 }
932 else {
933 p = (unsigned char) (gs * (t - render_min));
934 }
935
936 p += mingray;
937 }
938
939 if (invert) {
940 p = mingray + maxgray - p;
941 }
942
943 data[i * asrgb + j * bpl] = p;
944
945 if (hist) histd[p]++;
946
947 l += addi;
948 remx += addr;
949
950 if (remx > scale_n) {
951 remx -= scale_n;
952 l++;
953 }
954 }
955
956 l = br + addi * nx;
957 remy += addr;
958
959 if (remy > scale_n) {
960 remy -= scale_n;
961 l += nx;
962 }
963 }
964 }
965 }
966
967 // this replicates r -> g,b
968 if (asrgb == 3 && !(flags & 16)) {
969 for (int j=ymin*bpl; j <= ymax*bpl; j+=bpl) {
970 for (int i=xmin; i<xsize*3; i+=3) {
971 data[i+j+1] = data[i+j+2] = data[i+j];
972 }
973 }
974 }
975
976 if (asrgb == 4 && !(flags & 16)) {
977 for (int j=ymin*bpl; j <= ymax*bpl; j+=bpl) {
978 for (int i=xmin; i<xsize*4; i+=4) {
979 data[i+j+1] = data[i+j+2] = data[i+j+3] = data[i+j];
980 data[i+j+3] = 255;
981 }
982 }
983 }
984
985 EXITFUNC;
986
987 // ok, ok, not the most efficient place to do this, but it works
988 if (invy) {
989 for (int y=0; y<iysize/2; y++)
990 for (int x=0; x<ixsize; x++)
991 std::swap(ret[y*bpl+x], ret[(iysize-y-1)*bpl+x]);
992 }
993
994 // return PyString_FromStringAndSize((const char*) data,iysize*bpl);
995
996 if (flags & 16) {
997 glDrawPixels(ixsize, iysize, GL_LUMINANCE, GL_UNSIGNED_BYTE,
998 (const GLvoid *)ret.data());
999 }
1000
1001 delete [] grayhe;
1002 delete [] gaussianpdf;
1003 delete [] gaussiancdf;
1004 delete [] gaussianlookup;
1005 delete [] graypdftwo;
1006
1007 return ret;
1008}
1009
1010// DEPRECATED
1011
1012unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc,
1013 unsigned int tex_id, bool surface_face_z, bool recontour)
1014{
1015 if (mc->_isodl != 0) glDeleteLists(mc->_isodl,1);
1016 if (recontour) mc->calculate_surface();
1017 if (surface_face_z) mc->surface_face_z();
1018 if (mc->getRGBmode()) mc->color_vertices();
1019
1020#if MARCHING_CUBES_DEBUG
1021 cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() <<
1022 " points and " << nn.elem() << " normals to render in generate dl" << endl;
1023#endif
1024
1025 int maxf;
1026//#ifdef _WIN32
1027// glGetIntegerv(GL_MAX_ELEMENTS_INDICES_WIN, & maxf);
1028//#else
1029 glGetIntegerv(GL_MAX_ELEMENTS_INDICES, & maxf);
1030//#endif //_WIN32
1031
1032#if MARCHING_CUBES_DEBUG
1033 int maxv;
1034 glGetIntegerv(GL_MAX_ELEMENTS_VERTICES, & maxv);
1035
1036 cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
1037 cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
1038#endif
1039
1040 if (recontour) {
1041 for (unsigned int i = 0; i < mc->ff.elem(); ++i) mc->ff[i] /= 3;
1042 }
1043
1044 if (maxf % 3 != 0) {
1045 maxf = maxf - (maxf%3);
1046 }
1047
1048 if (tex_id != 0) {
1049 // Normalize the coordinates to be on the interval 0,1
1050 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(),
1051 1.0f/(float)mc->_emdata->get_ysize(),
1052 1.0f/(float)mc->_emdata->get_zsize());
1053 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
1054 glDisableClientState(GL_NORMAL_ARRAY);
1055 glTexCoordPointer(3, GL_FLOAT, 0, mc->pp.get_data());
1056 }
1057 else {
1058 glEnableClientState(GL_NORMAL_ARRAY);
1059 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
1060 glNormalPointer(GL_FLOAT, 0, mc->nn.get_data());
1061 }
1062
1063 glEnableClientState(GL_VERTEX_ARRAY);
1064 glVertexPointer(3, GL_FLOAT, 0, mc->pp.get_data());
1065
1066 if (mc->getRGBmode()) {
1067 glEnableClientState(GL_COLOR_ARRAY);
1068 glColorPointer(3, GL_FLOAT, 0, mc->cc.get_data());
1069 }
1070
1071 mc->_isodl = glGenLists(1);
1072
1073#if MARCHING_CUBES_DEBUG
1074 int time0 = clock();
1075#endif
1076
1077 glNewList(mc->_isodl,GL_COMPILE);
1078
1079 if (tex_id != 0) {
1080 glEnable(GL_TEXTURE_3D);
1081 glBindTexture(GL_TEXTURE_3D, tex_id);
1082 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
1083 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
1084 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
1085 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
1086 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1087 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
1088 }
1089
1090 // Drawing range elements based on the output of
1091 // glGetIntegerv(GL_MAX_ELEMENTS_INDICES, & maxf);
1092 // Saved about 60% of the time... drawRange should probably always be true
1093 bool drawRange = true;
1094
1095 if (drawRange == false) {
1096 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT, mc->ff.get_data());
1097 } else {
1098 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
1099 {
1100 if ((i+maxf) > mc->ff.elem())
1101 glDrawElements(GL_TRIANGLES, mc->ff.elem()-i, GL_UNSIGNED_INT, &(mc->ff[i]));
1102 else
1103 glDrawElements(GL_TRIANGLES, maxf, GL_UNSIGNED_INT, &(mc->ff[i]));
1104
1105 // glDrawRangeElements is part of the extensions,
1106 // we might want to experiment with its performance at some stage,
1107 // so please leave this code here, commented out.
1108 // This is an either or situation, so if glDrawRangeElements is used,
1109 // glDrawElements above would have to be commented out.
1110 // glDrawRangeElements(GL_TRIANGLES, 0, 0, maxf, GL_UNSIGNED_INT, &ff[i]);
1111 }
1112 }
1113
1114 if (tex_id != 0) glDisable(GL_TEXTURE_3D);
1115
1116 glEndList();
1117
1118#if MARCHING_CUBES_DEBUG
1119 int time1 = clock();
1120
1121 cout << "It took " << (time1-time0) << " " <<
1122 (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
1123#endif
1124
1125 return mc->_isodl;
1126}
1127
1128void GLUtil::contour_isosurface(MarchingCubes* mc)
1129{
1130 mc->calculate_surface();
1131
1132 // What does this do???
1133 for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
1134
1135 // Need to rebind data (to the GPU)
1136 mc->needtobind = true;
1137}
1138
1139void GLUtil::render_using_VBOs(MarchingCubes* mc, unsigned int tex_id,
1140 bool surface_face_z)
1141{
1142 // In current version Texture is not supported b/c it is not used... EVER
1143
1144#ifdef _WIN32
1145 typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
1146 PFNGLGENBUFFERSPROC glGenBuffers;
1147 glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
1148
1149 typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
1150 PFNGLISBUFFERPROC glIsBuffer;
1151 glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
1152#endif //_WIN32
1153
1154 if (surface_face_z) mc->surface_face_z();
1155
1156 // Bug in OpenGL, sometimes glGenBuffers doesn't work.
1157 // Try again, if we still fail then return...
1158
1159// printf("%d %d %d %d ",mc->buffer[0],mc->buffer[1],mc->buffer[2],mc->buffer[3]);
1160 if (!glIsBuffer(mc->buffer[0])) glGenBuffers(4, mc->buffer);
1161// printf("%d %d %d %d\n",mc->buffer[0],mc->buffer[1],mc->buffer[2],mc->buffer[3]);
1162
1163 // whenever something changes, like color mode or color scale (or threshold),
1164 // we need to recolor
1165 if (mc->getRGBmode() && (mc->rgbgenerator.getNeedToRecolor() ||
1166 mc->needtobind)) {
1167 mc->color_vertices();
1168 mc->needtobind = true;
1169 }
1170
1171 int maxf;
1172 glGetIntegerv(GL_MAX_ELEMENTS_INDICES, & maxf);
1173
1174 if (maxf % 3 != 0) {
1175 maxf = maxf - (maxf%3);
1176 }
1177
1178#ifdef _WIN32
1179 typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
1180 PFNGLBINDBUFFERPROC glBindBuffer;
1181 glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
1182
1183 typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size,
1184 const GLvoid *data, GLenum usage);
1185 PFNGLBUFFERDATAPROC glBufferData;
1186 glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
1187#endif //_WIN32
1188
1189 // Normal vectors
1190 glEnableClientState(GL_NORMAL_ARRAY);
1191 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
1192 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[2]);
1193
1194 if (mc->needtobind) {
1195 glBufferData(GL_ARRAY_BUFFER, 4*mc->nn.elem(), mc->nn.get_data(),
1196 GL_STATIC_DRAW);
1197 }
1198
1199 glNormalPointer(GL_FLOAT,0,0);
1200
1201 // Vertex vectors
1202 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[0]);
1203
1204 if (mc->needtobind) {
1205 glBufferData(GL_ARRAY_BUFFER, 4*mc->pp.elem(), mc->pp.get_data(),
1206 GL_STATIC_DRAW);
1207 }
1208
1209 glEnableClientState(GL_VERTEX_ARRAY);
1210 glVertexPointer(3,GL_FLOAT,0,0);
1211
1212 // Color vectors
1213 if (mc->getRGBmode()) {
1214 glEnableClientState(GL_COLOR_ARRAY);
1215 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[3]);
1216
1217 if (mc->needtobind) {
1218 glBufferData(GL_ARRAY_BUFFER, 4*mc->cc.elem(), mc->cc.get_data(),
1219 GL_STATIC_DRAW);
1220 }
1221
1222 glColorPointer(3,GL_FLOAT,0, 0);
1223 }
1224
1225 // Indices
1226 glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mc->buffer[1]);
1227
1228 if (mc->needtobind) {
1229 glBufferData(GL_ELEMENT_ARRAY_BUFFER, 4*mc->ff.elem(), mc->ff.get_data(),
1230 GL_STATIC_DRAW);
1231 }
1232
1233 // This lets us know if buffers an not implmemted
1234 if (!glIsBuffer(mc->buffer[0])) {
1235 cout << "Can't Generate GL Vertex Buffers. glGenBuffer error" << endl;
1236
1237 return;
1238 }
1239
1240 // finally draw the elemenets
1241 glDrawElements(GL_TRIANGLES, mc->ff.elem(), GL_UNSIGNED_INT, 0);
1242
1243 // No longer need to bind data to the GPU
1244
1245 mc->needtobind = false;
1246}
1247
1248static void transpose_4_by_4_matrix (vector <float> & v)
1249{
1250 std::swap (v [ 1], v [ 4]);
1251 std::swap (v [ 2], v [ 8]);
1252 std::swap (v [ 3], v [12]);
1253 std::swap (v [ 6], v [ 9]);
1254 std::swap (v [ 7], v [13]);
1255 std::swap (v [11], v [14]);
1256}
1257
1258void GLUtil::glLoadMatrix (const Transform & xform)
1259{
1260 vector <float> xformlist = xform.get_matrix_4x4 ();
1261
1262 transpose_4_by_4_matrix (xformlist);
1263
1264 glLoadMatrixf (reinterpret_cast <GLfloat *> (& xformlist [0]));
1265}
1266
1267void GLUtil::glMultMatrix (const Transform & xform)
1268{
1269 vector <float> xformlist = xform.get_matrix_4x4 ();
1270
1271 transpose_4_by_4_matrix (xformlist);
1272
1273 glMultMatrixf (reinterpret_cast <GLfloat *> (& xformlist [0]));
1274}
1275
1276// This draws a bounding box, or any box
1277void GLUtil::glDrawBoundingBox(float width, float height, float depth)
1278{
1279 float x = width/2.0f;
1280 float y = height/2.0f;
1281 float z = depth/2.0f;
1282
1283 float vertices[72] = {-x,-y,-z, x,-y,-z, x,-y,-z, x,y,-z, x,y,-z, -x,y,-z, -x,y,-z, -x,-y,-z,
1284 -x,-y,-z, -x,-y,z, x,-y,-z, x,-y,z, x,y,-z, x,y,z, -x,y,-z, -x,y,z,
1285 -x,-y,z, x,-y,z, x,-y,z, x,y,z, x,y,z, -x,y,z, -x,y,z, -x,-y,z };
1286
1287 glBegin(GL_LINES);
1288 for (int i=0; i<72; i+=3) glVertex3f(vertices[i],vertices[i+1],vertices[i+2]);
1289 glEnd();
1290
1291// float vertices[24] = {-w2, h2, d2, w2, h2, d2, w2, -h2, d2, -w2,
1292// -h2, d2, -w2, h2, -d2, w2, h2, -d2, w2, -h2, -d2, -w2, -h2, -d2};
1293// int indices[24] = {0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6,
1294// 6, 7, 7, 4, 0, 4, 3, 7, 1, 5, 2, 6};
1295//
1296// #ifdef _WIN32
1297// typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
1298// PFNGLGENBUFFERSPROC glGenBuffers;
1299// glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
1300//
1301// typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
1302// PFNGLISBUFFERPROC glIsBuffer;
1303// glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
1304//
1305// typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
1306// PFNGLBINDBUFFERPROC glBindBuffer;
1307// glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
1308//
1309// typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size,
1310// const GLvoid *data, GLenum usage);
1311// PFNGLBUFFERDATAPROC glBufferData;
1312// glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
1313// #endif //_WIN32
1314//
1315// if (glIsBuffer(GLUtil::buffer[0]) == GL_FALSE) {
1316// glGenBuffers(2, GLUtil::buffer);
1317// }
1318//
1319// // Could use dirty bit here but not worth my time to implment
1320//
1321// glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
1322// glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
1323// glEnableClientState(GL_VERTEX_ARRAY);
1324// glVertexPointer(3,GL_FLOAT,0,0);
1325//
1326// glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, GLUtil::buffer[1]);
1327// glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices,
1328// GL_STATIC_DRAW);
1329// glDrawElements(GL_LINES,24,GL_UNSIGNED_INT,0);
1330}
1331
1332void GLUtil::glDrawDisk(float radius, int spokes)
1333{
1334#ifdef _WIN32
1335 typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
1336 PFNGLGENBUFFERSPROC glGenBuffers;
1337 glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
1338
1339 typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
1340 PFNGLISBUFFERPROC glIsBuffer;
1341 glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
1342
1343 typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
1344 PFNGLBINDBUFFERPROC glBindBuffer;
1345 glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
1346
1347 typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size,
1348 const GLvoid *data, GLenum usage);
1349 PFNGLBUFFERDATAPROC glBufferData;
1350 glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
1351#endif //_WIN32
1352
1353 //This code is experimental
1354 int arraysize = 4*pow((double)2, (double)spokes);
1355 int sideofarray = pow((double)2, (double)spokes);
1356
1357// float vertices[3*arraysize + 3];
1358 vector<float> vertices(3*arraysize + 3);
1359
1360 // last vertex is center
1361 vertices[3*arraysize] = 0.0;
1362 vertices[3*arraysize+1] = 0.0;
1363 vertices[3*arraysize+2] = 0.0;
1364
1365 vertices[0] = 0.0;
1366 vertices[1] = radius;
1367 vertices[2] = 0.0;
1368
1369 vertices[sideofarray*3] = radius;
1370 vertices[sideofarray*3 + 1] = 0.0;
1371 vertices[sideofarray*3 + 2] = 0.0;
1372
1373 vertices[sideofarray*6] = 0.0;
1374 vertices[sideofarray*6 + 1] = -radius;
1375 vertices[sideofarray*6 + 2] = 0.0;
1376
1377 vertices[sideofarray*9] = -radius;
1378 vertices[sideofarray*9 + 1] = 0.0;
1379 vertices[sideofarray*9 + 2] = 0.0;
1380
1381 // This could also be implemented recursively
1382 for (int step = 0; step < spokes; step++) {
1383 // starting location
1384 int x = sideofarray/pow(2.0,(double)(step+1));
1385
1386 for (int i = 1; i <= 4*pow(2.0,(double)step); i++) {
1387 // take the necessary steps
1388 int index = x + 2*x*(i-1);
1389 int index_f = (index + x) % arraysize;
1390 int index_i = index - x;
1391
1392 cout << index << " " << index_f << " " << index_i << endl;
1393
1394 // need to resclae length to that of radius
1395 vertices[index_f*3] = (vertices[index_f*3] - vertices[index_i*3])/2.0f;
1396 vertices[index_f*3 + 1] = (vertices[index_f*3 + 1] - vertices[index_i*3 + 1])/2.0f;
1397 vertices[index_f*3 + 2] = (vertices[index_f*3 + 2] - vertices[index_i*3 + 2])/2.0f;
1398 }
1399 }
1400
1401 // GL stuff
1402 if (glIsBuffer(GLUtil::buffer[0]) == GL_FALSE) {
1403 glGenBuffers(2, GLUtil::buffer);
1404 }
1405
1406 // Could use dirty bit here but not worth my time to implement
1407 glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
1408
1409 // glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
1410 glBufferData(GL_ARRAY_BUFFER, vertices.size(), &(vertices[0]),
1411 GL_STATIC_DRAW);
1412 glEnableClientState(GL_VERTEX_ARRAY);
1413 glVertexPointer(3, GL_FLOAT, 0, 0);
1414
1415 // Code to select indices
1416}
1417#endif // USE_OPENGL
EMAN::ColorRGBGenerator::getNeedToRecolor
bool getNeedToRecolor()
Lets us know if we need to recalculate colors.
Definition: marchingcubes.h:243
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::nx
int nx
image size
Definition: emdata.h:848
EMAN::GLUtil::mx_bbox
static void mx_bbox(const vector< float > &data, const vector< float > &text_color, const vector< float > &bg_color)
EMAN::GLUtil::render_using_VBOs
static void render_using_VBOs(MarchingCubes *mc, unsigned int tex_id=0, bool surface_face_z=false)
Render a isosurface using buffer objects, this uses non-deprecated methods and improves performance.
EMAN::GLUtil::colored_rectangle
static void colored_rectangle(const vector< float > &data, const float &alpha, const bool center_point=false)
EMAN::GLUtil::render_amp8
static EMBytes render_amp8(EMData *emdata, int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, float gamma, int flags)
Render the image into an 8-bit image.
EMAN::GLUtil::glDrawBoundingBox
static void glDrawBoundingBox(float width, float height, float depth)
Draw a bounding box.
EMAN::GLUtil::glMultMatrix
static void glMultMatrix(const Transform &xform)
Mult a EMAN style transform to open GL w/o having to go through python Calls glMultTransposeMatrix ra...
EMAN::GLUtil::nearest_projected_points
static int nearest_projected_points(const vector< float > &model_matrix, const vector< float > &proj_matrix, const vector< int > &view_matrix, const vector< Vec3f > &points, const float mouse_x, const float mouse_y, const float &nearnes)
Determine the intersection of ....
EMAN::GLUtil::render_amp8_gl_texture
static unsigned int render_amp8_gl_texture(EMData *emdata, int x0, int y0, int ixsize, int iysize, int bpl, float scale, int mingray, int maxgray, float render_min, float render_max, float gamma, int flags)
create an OpenGL texture using render_amp8
EMAN::GLUtil::contour_isosurface
static void contour_isosurface(MarchingCubes *mc)
Recountour isosurface, for use with VBOs.
EMAN::GLUtil::glLoadMatrix
static void glLoadMatrix(const Transform &xform)
Load a EMAN style transform to open GL w/o having to go through python Calls glLoadTransposeMatrix ra...
EMAN::GLUtil::glDrawDisk
static void glDrawDisk(float radius, int spokes)
EMAN::GLUtil::gen_gl_texture
static unsigned int gen_gl_texture(const EMData *const emdata, GLenum format=GL_LUMINANCE)
create an OpenGL texture
EMAN::GLUtil::buffer
static GLuint buffer[2]
Definition: glutil.h:133
EMAN::GLUtil::get_isosurface_dl
static unsigned long get_isosurface_dl(MarchingCubes *mc, unsigned int tex_id=0, bool surface_face_z=false, bool recontour=true)
Get an isosurface display list Traverses the tree, marches the cubes, and renders a display list usin...
EMAN::Isosurface::_emdata
EMData * _emdata
Definition: isosurface.h:84
EMAN::MarchingCubes::calculate_surface
void calculate_surface()
Calculate and generate the entire set of vertices and normals using current states Calls draw_cube(0,...
EMAN::MarchingCubes::_isodl
unsigned long _isodl
Definition: marchingcubes.h:370
EMAN::MarchingCubes::rgbgenerator
ColorRGBGenerator rgbgenerator
Color by radius generator.
Definition: marchingcubes.h:453
EMAN::MarchingCubes::cc
CustomVector< float > cc
Definition: marchingcubes.h:447
EMAN::MarchingCubes::nn
CustomVector< float > nn
Definition: marchingcubes.h:449
EMAN::MarchingCubes::getRGBmode
int getRGBmode()
Return RGB mode.
Definition: marchingcubes.h:351
EMAN::MarchingCubes::ff
CustomVector< unsigned int > ff
Definition: marchingcubes.h:450
EMAN::MarchingCubes::color_vertices
void color_vertices()
Color the vertices.
EMAN::MarchingCubes::pp
CustomVector< float > pp
.Custom vectors for storing points, normals and faces
Definition: marchingcubes.h:446
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::get_matrix_4x4
vector< float > get_matrix_4x4() const
Get the 4x4 transformation matrix using a vector.
Definition: transform.cpp:192
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
EMAN
E2Exception class.
Definition: aligner.h:40
EMAN::Vec3f
Vec3< float > Vec3f
Definition: vec3.h:693
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
Generated by doxygen 1.9.2