00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifdef EMAN2_USING_OPENGL
00037
00038 #include "glutil.h"
00039 #include "emdata.h"
00040 #include "marchingcubes.h"
00041
00042 #ifndef _WIN32
00043 #ifndef GL_GLEXT_PROTOTYPES
00044 #define GL_GLEXT_PROTOTYPES
00045 #endif //GL_GLEXT_PROTOTYPES
00046 #endif //_WIN32
00047
00048 #ifdef __APPLE__
00049 #include "OpenGL/gl.h"
00050 #include "OpenGL/glu.h"
00051 #include "OpenGL/glext.h"
00052 #else // WIN32, LINUX
00053 #include "GL/gl.h"
00054 #include "GL/glu.h"
00055 #include "GL/glext.h"
00056 #endif //__APPLE__
00057
00058 using namespace EMAN;
00059
00060 unsigned int GLUtil::gen_glu_mipmaps(const EMData* const emdata)
00061 {
00062 if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL mipmap without internally stored data");
00063 ENTERFUNC;
00064
00065 unsigned int tex_name;
00066 glGenTextures(1, &tex_name);
00067
00068 if ( emdata->ny == 1 && emdata->nz == 1 )
00069 {
00070 glBindTexture(GL_TEXTURE_1D, tex_name);
00071 gluBuild1DMipmaps(GL_TEXTURE_1D, GL_LUMINANCE, emdata->nx, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00072 } else if (emdata->nz == 1) {
00073 glBindTexture(GL_TEXTURE_2D, tex_name);
00074 gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, emdata->nx, emdata->ny, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00075 }
00076 else {
00077 #ifdef _WIN32
00078
00079 printf("3D OpenGL mipmaps are not available on this platform.\n");
00080 #else
00081 glBindTexture(GL_TEXTURE_3D, tex_name);
00082 gluBuild3DMipmaps(GL_TEXTURE_3D, GL_LUMINANCE, emdata->nx, emdata->ny, emdata->nz, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00083 #endif //_WIN32
00084 }
00085
00086 EXITFUNC;
00087 return tex_name;
00088 }
00089
00090 unsigned int GLUtil::gen_gl_texture(const EMData* const emdata)
00091 {
00092 if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL texture without internally stored data");
00093 ENTERFUNC;
00094
00095 unsigned int tex_name;
00096 glGenTextures(1, &tex_name);
00097
00098 if ( emdata->ny == 1 && emdata->nz == 1 )
00099 {
00100 glBindTexture(GL_TEXTURE_1D, tex_name);
00101 glTexImage1D(GL_TEXTURE_1D,0,GL_LUMINANCE,emdata->nx,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00102 } else if (emdata->nz == 1) {
00103 glBindTexture(GL_TEXTURE_2D, tex_name);
00104 glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,emdata->nx,emdata->ny,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00105 }
00106 else {
00107 glBindTexture(GL_TEXTURE_3D, tex_name);
00108 #ifdef _WIN32
00109 PFNGLTEXIMAGE3DPROC glTexImage3D;
00110 #endif
00111 glTexImage3D(GL_TEXTURE_3D,0, GL_LUMINANCE,emdata->nx,emdata->ny,emdata->nz,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00112 }
00113
00114 EXITFUNC;
00115 return tex_name;
00116 }
00117
00118 unsigned int GLUtil::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) {
00119
00120 string pixels = render_amp8(emdata, x0, y0, ixsize,iysize, bpl, scale, mingray, maxgray, render_min, render_max, gamma,flags);
00121
00122 unsigned int tex_name;
00123 glGenTextures(1, &tex_name);
00124
00125 glBindTexture(GL_TEXTURE_2D,tex_name);
00126 glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,ixsize,iysize,0,GL_LUMINANCE, GL_UNSIGNED_BYTE, pixels.c_str() );
00127
00128 return tex_name;
00129 }
00130
00131
00132 #ifdef GL_GLEXT_PROTOTYPES
00133 #undef GL_GLEXT_PROTOTYPES
00134 #endif
00135
00136 int GLUtil::nearest_projected_points(const vector<float>& model_matrix, const vector<float>& proj_matrix, const vector<int>& view_matrix,
00137 const vector<Vec3f>& points, const float mouse_x, const float mouse_y,const float& nearness)
00138 {
00139 double proj[16];
00140 double model[16];
00141 int view[4];
00142
00143 copy(proj_matrix.begin(), proj_matrix.end(), proj);
00144 copy(model_matrix.begin(), model_matrix.end(), model);
00145 copy(view_matrix.begin(), view_matrix.end(), view);
00146
00147 vector<Vec3f> unproj_points;
00148 double x,y,z;
00149 double r,s,t;
00150 for(vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
00151 r = (double) (*it)[0];
00152 s = (double) (*it)[1];
00153 t = (double) (*it)[2];
00154 gluProject(r,s,t, model, proj, view, &x,&y,&z);
00155 unproj_points.push_back(Vec3f(x,y,z));
00156 }
00157
00158 vector<int> intersections;
00159
00160 float n_squared = nearness*nearness;
00161 for(unsigned int i = 0; i < unproj_points.size(); ++i){
00162 Vec3f& v = unproj_points[i];
00163 float dx = v[0] - mouse_x;
00164 dx *= dx;
00165 float dy = v[1] - mouse_y;
00166 dy *= dy;
00167 if ((dx+dy) <= n_squared) intersections.push_back((int)i);
00168 }
00169
00170 int intersection = -1;
00171 float near_z = 0;
00172 for(vector<int>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
00173 if (intersection == -1 || unproj_points[*it][2] < near_z) {
00174 intersection = *it;
00175 near_z = unproj_points[*it][2];
00176 }
00177 }
00178
00179 return intersection;
00180 }
00181
00182 void GLUtil::colored_rectangle(const vector<float>& data,const float& alpha,const bool center_point){
00183
00184 glBegin(GL_LINE_LOOP);
00185 glColor4f(data[0],data[1],data[2],alpha);
00186 glVertex2i(data[3],data[4]);
00187 glVertex2i(data[5],data[4]);
00188 glVertex2i(data[5],data[6]);
00189 glVertex2i(data[3],data[6]);
00190 glEnd();
00191
00192 if (center_point) {
00193 glBegin(GL_POINTS);
00194 glVertex2f( (data[3]+data[5])/2, (data[4]+data[6])/2);
00195 glEnd();
00196 }
00197 }
00198
00199 void GLUtil::mx_bbox(const vector<float>& data, const vector<float>& text_color, const vector<float>& bg_color) {
00200 glColor4f(bg_color[0],bg_color[1],bg_color[2],bg_color[3]);
00201 glDisable(GL_TEXTURE_2D);
00202 glBegin(GL_QUADS);
00203
00204 glVertex3f(data[0]-1,data[1]-1,-.1);
00205 glVertex3f(data[3]+1,data[1]-1,-.1);
00206 glVertex3f(data[3]+1,data[4]+1,-.1);
00207 glVertex3f(data[0]-1,data[4]+1,-.1);
00208 glEnd();
00209 glEnable(GL_TEXTURE_2D);
00210 glColor4f(text_color[0],text_color[1],text_color[2],text_color[3]);
00211 }
00212
00213 std::string GLUtil::render_amp8(EMData* emdata, int x0, int y0, int ixsize, int iysize,
00214 int bpl, float scale, int mingray, int maxgray,
00215 float render_min, float render_max,float gamma,int flags)
00216 {
00217 ENTERFUNC;
00218
00219 int asrgb;
00220 int hist=(flags&2)/2;
00221 int invy=(flags&4)?1:0;
00222
00223 int nx = emdata->nx;
00224 int ny = emdata->ny;
00225
00226 int nxy = emdata->nx * emdata->ny;
00227
00228 if (emdata->get_ndim() > 2) {
00229 throw ImageDimensionException("1D/2D only");
00230 }
00231
00232 if (emdata->is_complex()) {
00233 emdata->ri2ap();
00234 }
00235
00236 if (render_max <= render_min) {
00237 render_max = render_min + 0.01f;
00238 }
00239
00240 if (gamma<=0) gamma=1.0;
00241
00242
00243
00244
00245
00246
00247
00248
00249 static int smg0=0,smg1=0;
00250 static float sgam=0;
00251 static unsigned char gammamap[4096];
00252 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00253 for (int i=0; i<4096; i++) {
00254 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00255 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00256 }
00257 }
00258 smg0=mingray;
00259 smg1=maxgray;
00260 sgam=gamma;
00261
00262 if (flags&8) asrgb=4;
00263 else if (flags&1) asrgb=3;
00264 else asrgb=1;
00265
00266 std::string ret=std::string();
00267
00268 ret.assign(iysize*bpl+hist*1024,char(mingray));
00269 unsigned char *data=(unsigned char *)ret.data();
00270 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00271 if (hist) {
00272 for (int i=0; i<256; i++) histd[i]=0;
00273 }
00274
00275 float rm = render_min;
00276 float inv_scale = 1.0f / scale;
00277 int ysize = iysize;
00278 int xsize = ixsize;
00279
00280 int ymin = 0;
00281 if (iysize * inv_scale > ny) {
00282 ymin = (int) (iysize - ny / inv_scale);
00283 }
00284
00285 float gs = (maxgray - mingray) / (render_max - render_min);
00286 float gs2 = 4095.999f / (render_max - render_min);
00287
00288 if (render_max < render_min) {
00289 gs = 0;
00290 rm = FLT_MAX;
00291 }
00292
00293 int dsx = -1;
00294 int dsy = 0;
00295 int remx = 0;
00296 int remy = 0;
00297 const int scale_n = 100000;
00298
00299 int addi = 0;
00300 int addr = 0;
00301 if (inv_scale == floor(inv_scale)) {
00302 dsx = (int) inv_scale;
00303 dsy = (int) (inv_scale * nx);
00304 }
00305 else {
00306 addi = (int) floor(inv_scale);
00307 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00308 }
00309
00310 int xmin = 0;
00311 if (x0 < 0) {
00312 xmin = (int) (-x0 / inv_scale);
00313 xsize -= (int) floor(x0 / inv_scale);
00314 x0 = 0;
00315 }
00316
00317 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00318 xsize = (int) ((nx - x0) / inv_scale + xmin);
00319 }
00320 int ymax = ysize - 1;
00321 if (y0 < 0) {
00322 ymax = (int) (ysize + y0 / inv_scale - 1);
00323 ymin += (int) floor(y0 / inv_scale);
00324 y0 = 0;
00325 }
00326
00327 if (xmin < 0) xmin = 0;
00328 if (ymin < 0) ymin = 0;
00329 if (xsize > ixsize) xsize = ixsize;
00330 if (ymax > iysize) ymax = iysize;
00331
00332 int lmax = nx * ny - 1;
00333
00334 int mid=nx*ny/2;
00335 float * image_data = emdata->get_data();
00336 if (emdata->is_complex()) {
00337 if (dsx != -1) {
00338 int l = y0 * nx;
00339 for (int j = ymax; j >= ymin; j--) {
00340 int ll = x0;
00341 for (int i = xmin; i < xsize; i++) {
00342 if (l + ll > lmax || ll >= nx - 2) break;
00343
00344 int k = 0;
00345 unsigned char p;
00346 if (ll >= nx / 2) {
00347 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00348 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00349 }
00350 else k = nx * ny - (l + 2 * ll) - 2;
00351 if (k>=mid) k-=mid;
00352 else k+=mid;
00353 float t = image_data[k];
00354 if (t <= rm) p = mingray;
00355 else if (t >= render_max) p = maxgray;
00356 else if (gamma!=1.0) {
00357 k=(int)(gs2 * (t-render_min));
00358 p = gammamap[k];
00359
00360
00361
00362
00363 }
00364 else {
00365 p = (unsigned char) (gs * (t - render_min));
00366 p += mingray;
00367 }
00368 data[i * asrgb + j * bpl] = p;
00369 if (hist) histd[p]++;
00370 ll += dsx;
00371 }
00372 l += dsy;
00373 }
00374 }
00375 else {
00376 remy = 10;
00377 int l = y0 * nx;
00378 for (int j = ymax; j >= ymin; j--) {
00379 int br = l;
00380 remx = 10;
00381 int ll = x0;
00382 for (int i = xmin; i < xsize - 1; i++) {
00383 if (l + ll > lmax || ll >= nx - 2) {
00384 break;
00385 }
00386 int k = 0;
00387 unsigned char p;
00388 if (ll >= nx / 2) {
00389 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00390 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00391 }
00392 else k = nx * ny - (l + 2 * ll) - 2;
00393 if (k>=mid) k-=mid;
00394 else k+=mid;
00395
00396 float t = image_data[k];
00397 if (t <= rm)
00398 p = mingray;
00399 else if (t >= render_max) {
00400 p = maxgray;
00401 }
00402 else if (gamma!=1.0) {
00403 k=(int)(gs2 * (t-render_min));
00404 p = gammamap[k];
00405
00406
00407
00408
00409 }
00410 else {
00411 p = (unsigned char) (gs * (t - render_min));
00412 p += mingray;
00413 }
00414 data[i * asrgb + j * bpl] = p;
00415 if (hist) histd[p]++;
00416 ll += addi;
00417 remx += addr;
00418 if (remx > scale_n) {
00419 remx -= scale_n;
00420 ll++;
00421 }
00422 }
00423 l = br + addi * nx;
00424 remy += addr;
00425 if (remy > scale_n) {
00426 remy -= scale_n;
00427 l += nx;
00428 }
00429 }
00430 }
00431 }
00432 else {
00433 if (dsx != -1) {
00434 int l = x0 + y0 * nx;
00435 for (int j = ymax; j >= ymin; j--) {
00436 int br = l;
00437 for (int i = xmin; i < xsize; i++) {
00438 if (l > lmax) {
00439 break;
00440 }
00441 int k = 0;
00442 unsigned char p;
00443 float t;
00444 if (dsx==1) t=image_data[l];
00445 else {
00446 t=0;
00447 for (int iii=0; iii<dsx; iii++) {
00448 for (int jjj=0; jjj<dsy; jjj+=nx) {
00449 t+=image_data[l+iii+jjj];
00450 }
00451 }
00452 t/=dsx*(dsy/nx);
00453 }
00454
00455 if (t <= rm) p = mingray;
00456 else if (t >= render_max) p = maxgray;
00457 else if (gamma!=1.0) {
00458 k=(int)(gs2 * (t-render_min));
00459 p = gammamap[k];
00460
00461
00462
00463
00464 }
00465 else {
00466 p = (unsigned char) (gs * (t - render_min));
00467 p += mingray;
00468 }
00469 data[i * asrgb + j * bpl] = p;
00470 if (hist) histd[p]++;
00471 l += dsx;
00472 }
00473 l = br + dsy;
00474 }
00475 }
00476 else {
00477 remy = 10;
00478 int l = x0 + y0 * nx;
00479 for (int j = ymax; j >= ymin; j--) {
00480 int addj = addi;
00481
00482
00483
00484 if (( l + addi*nx ) >= nxy ) {
00485 addj = (nxy-l)/nx;
00486 if (addj <= 0) continue;
00487 }
00488 int br = l;
00489 remx = 10;
00490 for (int i = xmin; i < xsize; i++) {
00491 if (l > lmax) break;
00492 int k = 0;
00493 unsigned char p;
00494 float t;
00495 if (addi<=1) t = image_data[l];
00496 else {
00497 t=0;
00498 for (int jjj=0; jjj<addj; jjj++) {
00499 for (int iii=0; iii<addi; iii++) {
00500 t+=image_data[l+iii+jjj*nx];
00501 }
00502 }
00503 t/=addi*addi;
00504 }
00505 if (t <= rm) p = mingray;
00506 else if (t >= render_max) p = maxgray;
00507 else if (gamma!=1.0) {
00508 k=(int)(gs2 * (t-render_min));
00509 p = gammamap[k];
00510
00511
00512
00513
00514 }
00515 else {
00516 p = (unsigned char) (gs * (t - render_min));
00517 p += mingray;
00518 }
00519 data[i * asrgb + j * bpl] = p;
00520 if (hist) histd[p]++;
00521 l += addi;
00522 remx += addr;
00523 if (remx > scale_n) {
00524 remx -= scale_n;
00525 l++;
00526 }
00527 }
00528 l = br + addi * nx;
00529 remy += addr;
00530 if (remy > scale_n) {
00531 remy -= scale_n;
00532 l += nx;
00533 }
00534 }
00535 }
00536 }
00537
00538
00539 if (asrgb==3) {
00540 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00541 for (int i=xmin; i<xsize*3; i+=3) {
00542 data[i+j+1]=data[i+j+2]=data[i+j];
00543 }
00544 }
00545 }
00546 if (asrgb==4) {
00547 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00548 for (int i=xmin; i<xsize*4; i+=4) {
00549 data[i+j+1]=data[i+j+2]=data[i+j+3]=data[i+j];
00550 data[i+j+3]=255;
00551 }
00552 }
00553 }
00554
00555 EXITFUNC;
00556
00557
00558 if (invy) {
00559 int x,y;
00560 char swp;
00561 for (y=0; y<iysize/2; y++) {
00562 for (x=0; x<ixsize; x++) {
00563 swp=ret[y*bpl+x];
00564 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00565 ret[(iysize-y-1)*bpl+x]=swp;
00566 }
00567 }
00568 }
00569
00570
00571 if (flags&16) {
00572 glDrawPixels(ixsize,iysize,GL_LUMINANCE,GL_UNSIGNED_BYTE,(const GLvoid *)ret.data());
00573 }
00574
00575 return ret;
00576 }
00577
00578
00579 unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z)
00580 {
00581 if ( mc->_isodl != 0 ) glDeleteLists(mc->_isodl,1);
00582
00583
00584 mc->calculate_surface();
00585 if ( surface_face_z ) mc->surface_face_z();
00586 #if MARCHING_CUBES_DEBUG
00587 cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() << " points and " << nn.elem() << " normals to render in generate dl" << endl;
00588 #endif
00589 int maxf;
00590
00591
00592
00593 glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00594
00595 #if MARCHING_CUBES_DEBUG
00596 int maxv;
00597 glGetIntegerv(GL_MAX_ELEMENTS_VERTICES,&maxv);
00598 cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
00599 cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
00600 #endif
00601
00602 for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00603
00604 if ( maxf % 3 != 0 )
00605 {
00606 maxf = maxf - (maxf%3);
00607 }
00608
00609 if ( tex_id != 0 ) {
00610
00611 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(), 1.0f/(float)mc->_emdata->get_ysize(), 1.0f/(float)mc->_emdata->get_zsize());
00612 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
00613 glDisableClientState(GL_NORMAL_ARRAY);
00614 glTexCoordPointer(3,GL_FLOAT,0,mc->pp.get_data());
00615 }
00616 else {
00617 glEnableClientState(GL_NORMAL_ARRAY);
00618 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00619 glNormalPointer(GL_FLOAT,0,mc->nn.get_data());
00620 }
00621
00622 glEnableClientState(GL_VERTEX_ARRAY);
00623 glVertexPointer(3,GL_FLOAT,0,mc->pp.get_data());
00624
00625 mc->_isodl = glGenLists(1);
00626
00627 #if MARCHING_CUBES_DEBUG
00628 int time0 = clock();
00629 #endif
00630 glNewList(mc->_isodl,GL_COMPILE);
00631
00632 if ( tex_id != 0 ) {
00633 glEnable(GL_TEXTURE_3D);
00634 glBindTexture(GL_TEXTURE_3D, tex_id);
00635 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
00636 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
00637 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
00638 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
00639 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
00640 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
00641 }
00642
00643
00644 bool drawRange = true;
00645 if ( drawRange == false ) {
00646 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT,mc->ff.get_data());
00647 } else {
00648 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
00649 {
00650 if ( (i+maxf) > mc->ff.elem())
00651 glDrawElements(GL_TRIANGLES,mc->ff.elem()-i,GL_UNSIGNED_INT,&(mc->ff[i]));
00652 else
00653 glDrawElements(GL_TRIANGLES,maxf,GL_UNSIGNED_INT,&(mc->ff[i]));
00654
00655
00656
00657
00658 }
00659 }
00660
00661 if ( tex_id != 0 ) glDisable(GL_TEXTURE_3D);
00662
00663 glEndList();
00664 #if MARCHING_CUBES_DEBUG
00665 int time1 = clock();
00666 cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
00667 #endif
00668 return mc->_isodl;
00669 }
00670
00671
00672 #endif // EMAN2_USING_OPENGL