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 #ifdef _WIN32
00039 #include <windows.h>
00040 #endif
00041
00042 #ifndef GL_GLEXT_PROTOTYPES
00043 #define GL_GLEXT_PROTOTYPES
00044 #endif //GL_GLEXT_PROTOTYPES
00045
00046 #include "glutil.h"
00047 #include "emdata.h"
00048 #include "marchingcubes.h"
00049
00050 #ifdef __APPLE__
00051 #include "OpenGL/gl.h"
00052 #include "OpenGL/glu.h"
00053 #include "OpenGL/glext.h"
00054 #else // WIN32, LINUX
00055 #include "GL/gl.h"
00056 #include "GL/glu.h"
00057 #include "GL/glext.h"
00058 #endif //__APPLE__
00059
00060 using namespace EMAN;
00061
00062
00063 GLuint GLUtil::buffer[2] = {0, 0};
00064
00065 unsigned int GLUtil::gen_glu_mipmaps(const EMData* const emdata)
00066 {
00067 if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL mipmap without internally stored data");
00068 ENTERFUNC;
00069
00070 unsigned int tex_name;
00071 glGenTextures(1, &tex_name);
00072
00073 if ( emdata->ny == 1 && emdata->nz == 1 )
00074 {
00075 glBindTexture(GL_TEXTURE_1D, tex_name);
00076 gluBuild1DMipmaps(GL_TEXTURE_1D, GL_LUMINANCE, emdata->nx, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00077 } else if (emdata->nz == 1) {
00078 glBindTexture(GL_TEXTURE_2D, tex_name);
00079 gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, emdata->nx, emdata->ny, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00080 }
00081 else {
00082 #ifdef _WIN32
00083
00084 printf("3D OpenGL mipmaps are not available on this platform.\n");
00085 #else
00086 glBindTexture(GL_TEXTURE_3D, tex_name);
00087 gluBuild3DMipmaps(GL_TEXTURE_3D, GL_LUMINANCE, emdata->nx, emdata->ny, emdata->nz, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00088 #endif //_WIN32
00089 }
00090
00091 EXITFUNC;
00092 return tex_name;
00093 }
00094
00095 unsigned int GLUtil::gen_gl_texture(const EMData* const emdata, GLenum format)
00096 {
00097 if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL texture without internally stored data");
00098 ENTERFUNC;
00099
00100 unsigned int tex_name;
00101 glGenTextures(1, &tex_name);
00102
00103 if ( emdata->ny == 1 && emdata->nz == 1 )
00104 {
00105 glBindTexture(GL_TEXTURE_1D, tex_name);
00106 glTexImage1D(GL_TEXTURE_1D,0,format,emdata->nx,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00107 } else if (emdata->nz == 1) {
00108 glBindTexture(GL_TEXTURE_2D, tex_name);
00109 glTexImage2D(GL_TEXTURE_2D,0,format,emdata->nx,emdata->ny,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00110 }
00111 else {
00112 glBindTexture(GL_TEXTURE_3D, tex_name);
00113 #ifdef _WIN32
00114 PFNGLTEXIMAGE3DPROC glTexImage3D;
00115 #endif
00116 glTexImage3D(GL_TEXTURE_3D,0, format,emdata->nx,emdata->ny,emdata->nz,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00117 }
00118
00119 EXITFUNC;
00120 return tex_name;
00121 }
00122
00123 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) {
00124
00125 string pixels = render_amp8(emdata, x0, y0, ixsize,iysize, bpl, scale, mingray, maxgray, render_min, render_max, gamma,flags);
00126
00127 unsigned int tex_name;
00128 glGenTextures(1, &tex_name);
00129
00130 glBindTexture(GL_TEXTURE_2D,tex_name);
00131 glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,ixsize,iysize,0,GL_LUMINANCE, GL_UNSIGNED_BYTE, pixels.c_str() );
00132
00133 return tex_name;
00134 }
00135
00136
00137 #ifdef GL_GLEXT_PROTOTYPES
00138 #undef GL_GLEXT_PROTOTYPES
00139 #endif
00140
00141 int GLUtil::nearest_projected_points(const vector<float>& model_matrix, const vector<float>& proj_matrix, const vector<int>& view_matrix,
00142 const vector<Vec3f>& points, const float mouse_x, const float mouse_y,const float& nearness)
00143 {
00144 double proj[16];
00145 double model[16];
00146 int view[4];
00147
00148 copy(proj_matrix.begin(), proj_matrix.end(), proj);
00149 copy(model_matrix.begin(), model_matrix.end(), model);
00150 copy(view_matrix.begin(), view_matrix.end(), view);
00151
00152 vector<Vec3f> unproj_points;
00153 double x,y,z;
00154 double r,s,t;
00155 for(vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
00156 r = (double) (*it)[0];
00157 s = (double) (*it)[1];
00158 t = (double) (*it)[2];
00159 gluProject(r,s,t, model, proj, view, &x,&y,&z);
00160 unproj_points.push_back(Vec3f(x,y,z));
00161 }
00162
00163 vector<int> intersections;
00164
00165 float n_squared = nearness*nearness;
00166 for(unsigned int i = 0; i < unproj_points.size(); ++i){
00167 Vec3f& v = unproj_points[i];
00168 float dx = v[0] - mouse_x;
00169 dx *= dx;
00170 float dy = v[1] - mouse_y;
00171 dy *= dy;
00172 if ((dx+dy) <= n_squared) intersections.push_back((int)i);
00173 }
00174
00175 int intersection = -1;
00176 float near_z = 0;
00177 for(vector<int>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
00178 if (intersection == -1 || unproj_points[*it][2] < near_z) {
00179 intersection = *it;
00180 near_z = unproj_points[*it][2];
00181 }
00182 }
00183
00184 return intersection;
00185 }
00186
00187 void GLUtil::colored_rectangle(const vector<float>& data,const float& alpha,const bool center_point){
00188
00189 glBegin(GL_LINE_LOOP);
00190 glColor4f(data[0],data[1],data[2],alpha);
00191 glVertex2i(data[3],data[4]);
00192 glVertex2i(data[5],data[4]);
00193 glVertex2i(data[5],data[6]);
00194 glVertex2i(data[3],data[6]);
00195 glEnd();
00196
00197 if (center_point) {
00198 glBegin(GL_POINTS);
00199 glVertex2f( (data[3]+data[5])/2, (data[4]+data[6])/2);
00200 glEnd();
00201 }
00202 }
00203
00204 void GLUtil::mx_bbox(const vector<float>& data, const vector<float>& text_color, const vector<float>& bg_color) {
00205 glDisable(GL_TEXTURE_2D);
00206 glBegin(GL_QUADS);
00207 if (bg_color.size()==4) glColor4f(bg_color[0],bg_color[1],bg_color[2],bg_color[3]);
00208 else glColor3f(bg_color[0],bg_color[1],bg_color[2]);
00209 glVertex3f(data[0]-1,data[1]-1,-.1f);
00210 glVertex3f(data[3]+1,data[1]-1,-.1f);
00211 glVertex3f(data[3]+1,data[4]+1,-.1f);
00212 glVertex3f(data[0]-1,data[4]+1,-.1f);
00213 glEnd();
00214
00215 if (text_color.size()==4) glColor4f(text_color[0],text_color[1],text_color[2],text_color[3]);
00216 else glColor3f(text_color[0],text_color[1],text_color[2]);
00217
00218 }
00219
00220 std::string GLUtil::render_amp8(EMData* emdata, int x0, int y0, int ixsize, int iysize,
00221 int bpl, float scale, int mingray, int maxgray,
00222 float render_min, float render_max,float gamma,int flags)
00223 {
00224 ENTERFUNC;
00225
00226 int asrgb;
00227 int hist=(flags&2)/2;
00228 int invy=(flags&4)?1:0;
00229
00230 int nx = emdata->nx;
00231 int ny = emdata->ny;
00232
00233 int nxy = emdata->nx * emdata->ny;
00234
00235 if (emdata->get_ndim() > 2) {
00236 throw ImageDimensionException("1D/2D only");
00237 }
00238
00239 if (emdata->is_complex()) {
00240 emdata->ri2ap();
00241 }
00242
00243 if (render_max <= render_min) {
00244 render_max = render_min + 0.01f;
00245 }
00246
00247 if (gamma<=0) gamma=1.0;
00248
00249
00250
00251
00252
00253
00254
00255
00256 static int smg0=0,smg1=0;
00257 static float sgam=0;
00258 static unsigned char gammamap[4096];
00259 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00260 for (int i=0; i<4096; i++) {
00261 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00262 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00263 }
00264 }
00265 smg0=mingray;
00266 smg1=maxgray;
00267 sgam=gamma;
00268
00269 if (flags&8) asrgb=4;
00270 else if (flags&1) asrgb=3;
00271 else asrgb=1;
00272
00273 std::string ret=std::string();
00274
00275 ret.assign(iysize*bpl+hist*1024,char(mingray));
00276 unsigned char *data=(unsigned char *)ret.data();
00277 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00278 if (hist) {
00279 for (int i=0; i<256; i++) histd[i]=0;
00280 }
00281
00282 float rm = render_min;
00283 float inv_scale = 1.0f / scale;
00284 int ysize = iysize;
00285 int xsize = ixsize;
00286
00287 int ymin = 0;
00288 if (iysize * inv_scale > ny) {
00289 ymin = (int) (iysize - ny / inv_scale);
00290 }
00291
00292 float gs = (maxgray - mingray) / (render_max - render_min);
00293 float gs2 = 4095.999f / (render_max - render_min);
00294
00295 if (render_max < render_min) {
00296 gs = 0;
00297 rm = FLT_MAX;
00298 }
00299
00300 int dsx = -1;
00301 int dsy = 0;
00302 int remx = 0;
00303 int remy = 0;
00304 const int scale_n = 100000;
00305
00306 int addi = 0;
00307 int addr = 0;
00308 if (inv_scale == floor(inv_scale)) {
00309 dsx = (int) inv_scale;
00310 dsy = (int) (inv_scale * nx);
00311 }
00312 else {
00313 addi = (int) floor(inv_scale);
00314 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00315 }
00316
00317 int xmin = 0;
00318 if (x0 < 0) {
00319 xmin = (int) (-x0 / inv_scale);
00320 xsize -= (int) floor(x0 / inv_scale);
00321 x0 = 0;
00322 }
00323
00324 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00325 xsize = (int) ((nx - x0) / inv_scale + xmin);
00326 }
00327 int ymax = ysize - 1;
00328 if (y0 < 0) {
00329 ymax = (int) (ysize + y0 / inv_scale - 1);
00330 ymin += (int) floor(y0 / inv_scale);
00331 y0 = 0;
00332 }
00333
00334 if (xmin < 0) xmin = 0;
00335 if (ymin < 0) ymin = 0;
00336 if (xsize > ixsize) xsize = ixsize;
00337 if (ymax > iysize) ymax = iysize;
00338
00339 int lmax = nx * ny - 1;
00340
00341 int mid=nx*ny/2;
00342 float * image_data = emdata->get_data();
00343
00345
00346 const int rangemax=4082;
00347 int gaussianwide=5000;
00348 unsigned int* grayhe=NULL;
00349 float* gaussianpdf=NULL;
00350 float* gaussiancdf=NULL;
00351 int* gaussianlookup=NULL;
00352
00353 unsigned int* graypdftwo=NULL;
00354
00355 bool binflag=0;
00356 int binlocation=-1;
00357
00358 if (flags&32){
00359 int graypdf[rangemax]={0};
00360 int graycdf[rangemax-2]={0};
00361
00362
00363 graypdftwo = new unsigned int[maxgray-mingray];
00364 grayhe = new unsigned int[rangemax-2];
00365
00366
00367 for (int i=0; i<(int)(rangemax-2); i++) {
00368 grayhe[i] = 0;
00369 }
00370
00371 for (int i=0; i<(maxgray-mingray); i++) {
00372 graypdftwo[i]=0;
00373 }
00374
00375 if (dsx != -1) {
00376 int l = x0 + y0 * nx;
00377 for (int j = ymax; j >= ymin; j--) {
00378 int br = l;
00379 for (int i = xmin; i < xsize; i++) {
00380 if (l > lmax) {
00381 break;
00382 }
00383 int k = 0;
00384 int p;
00385 float t;
00386 if (dsx==1) t=image_data[l];
00387 else {
00388 t=0;
00389 for (int iii=0; iii<dsx; iii++) {
00390 for (int jjj=0; jjj<dsy; jjj+=nx) {
00391 t+=image_data[l+iii+jjj];
00392 }
00393 }
00394 t/=dsx*(dsy/nx);
00395 }
00396
00397 if (t <= rm) graypdf[0]++;
00398 else if (t >= render_max) graypdf[rangemax-1]++;
00399 else {
00400 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
00401 graypdftwo[(unsigned char) (gs * (t - render_min))]++;
00402 }
00403 l += dsx;
00404 }
00405 l = br + dsy;
00406 }
00407 }
00408 else {
00409 remy = 10;
00410 int l = x0 + y0 * nx;
00411 for (int j = ymax; j >= ymin; j--) {
00412 int addj = addi;
00413
00414
00415
00416 if (( l + addi*nx ) >= nxy ) {
00417 addj = (nxy-l)/nx;
00418 if (addj <= 0) continue;
00419 }
00420 int br = l;
00421 remx = 10;
00422 for (int i = xmin; i < xsize; i++) {
00423 if (l > lmax) break;
00424 int k = 0;
00425 int p;
00426 float t;
00427 if (addi<=1) t = image_data[l];
00428 else {
00429 t=0;
00430 for (int jjj=0; jjj<addj; jjj++) {
00431 for (int iii=0; iii<addi; iii++) {
00432 t+=image_data[l+iii+jjj*nx];
00433 }
00434 }
00435 t/=addi*addi;
00436 }
00438 if (t <= rm) graypdf[0]++;
00439 else if (t >= render_max) graypdf[rangemax-1]++;
00440 else {
00441 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
00442 }
00443
00444 data[i * asrgb + j * bpl] = p;
00445 if (hist) histd[p]++;
00446 l += addi;
00447 remx += addr;
00448 if (remx > scale_n) {
00449 remx -= scale_n;
00450 l++;
00451 }
00452 }
00453 l = br + addi * nx;
00454 remy += addr;
00455 if (remy > scale_n) {
00456 remy -= scale_n;
00457 l += nx;
00458 }
00459 }
00460 }
00461
00462 for (int i=0; i<(rangemax-2); i++) {
00463 for (int j=0;j<(i+1);j++) {
00464 graycdf[i]=graycdf[i]+graypdf[j+1];
00465 }
00466 }
00467
00468
00469 binflag=0;
00470 for (int i=1; i<(maxgray-mingray); i++) {
00471
00472 if (((float)graypdftwo[i]/graycdf[rangemax-3])>0.2){
00473 binflag=1;
00474 binlocation=i;
00475 break;
00476 }
00477 }
00478
00479 if (binflag==1){
00480 for (int i=(binlocation*16+1); i<((binlocation+1)*16); i++) {
00481 graypdf[i]=0;
00482
00483 }
00484
00485 for (int i=0; i<(rangemax-2); i++) {
00486 graycdf[i]=0;
00487 }
00488
00489 for (int i=0; i<(rangemax-2); i++) {
00490 for (int j=0;j<(i+1);j++) {
00491 graycdf[i]=graycdf[i]+graypdf[j+1];
00492 }
00493 }
00494 }
00495
00496
00497
00498 float mean=abs(rangemax-2)/2;
00499 float standdv=abs(mean)/3;
00500 gaussianpdf=new float[rangemax-2];
00501 gaussiancdf=new float[rangemax-2];
00502 gaussianlookup=new int[rangemax-2];
00503
00504 for (int i=0; i<(rangemax-2); i++){
00505 gaussianpdf[i]=exp(-(i-mean)*(i-mean)/(2*standdv*standdv))/sqrt(standdv * standdv * 2 * M_PI);
00506
00507 if(i!=0){
00508 gaussiancdf[i]=gaussiancdf[i-1]+gaussianpdf[i];
00509 }
00510 else{
00511 gaussiancdf[i]=gaussianpdf[i];
00512 }
00513 }
00514
00515 for (int i=0; i<(rangemax-2); i++) {
00516 gaussiancdf[i]=graycdf[rangemax-3]*gaussiancdf[i]/gaussiancdf[rangemax-3];
00517 }
00518
00519 for (int i=0; i<(rangemax-2); i++){
00520 for (int j=0; j<(rangemax-2); j++){
00521 if (graycdf[i]<=gaussiancdf[j]){
00522 gaussianlookup[i]=j;
00523 break;
00524 }
00525 }
00526 }
00527
00528 for (int i=0; i<(rangemax-2); i++) {
00529 grayhe[i]=floor(0.5+(((double)(rangemax-3)*graycdf[i])/graycdf[rangemax-3]));
00530 }
00531
00532 }
00533
00535
00536
00537 if (emdata->is_complex()) {
00538 if (dsx != -1) {
00539 int l = y0 * nx;
00540 for (int j = ymax; j >= ymin; j--) {
00541 int ll = x0;
00542 for (int i = xmin; i < xsize; i++) {
00543 if (l + ll > lmax || ll >= nx - 2) break;
00544
00545 int k = 0;
00546 unsigned char p;
00547 if (ll >= nx / 2) {
00548 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00549 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00550 }
00551 else k = nx * ny - (l + 2 * ll) - 2;
00552 if (k>=mid) k-=mid;
00553 else k+=mid;
00554 float t = image_data[k];
00555 int ph;
00556
00557 if (flags&16 && asrgb>2) {
00558
00559 if (ll >= nx / 2) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00560 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00561 }
00562 if (t <= rm) p = mingray;
00563 else if (t >= render_max) p = maxgray;
00564 else if (gamma!=1.0) {
00565 k=(int)(gs2 * (t-render_min));
00566 p = gammamap[k];
00567
00568
00569
00570
00571 }
00572 else {
00573 p = (unsigned char) (gs * (t - render_min));
00574 p += mingray;
00575 }
00576
00577 if (flags&16 && asrgb>2) {
00578 if (ph<256) {
00579 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00580 data[i * asrgb + j * bpl+1] = p*ph/256;
00581 data[i * asrgb + j * bpl+2] = 0;
00582 }
00583 else if (ph<512) {
00584 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00585 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00586 data[i * asrgb + j * bpl] = 0;
00587 }
00588 else {
00589 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00590 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00591 data[i * asrgb + j * bpl+1] = 0;
00592 }
00593 }
00594 else data[i * asrgb + j * bpl] = p;
00595 if (hist) histd[p]++;
00596 ll += dsx;
00597 }
00598 l += dsy;
00599 }
00600 }
00601 else {
00602 remy = 10;
00603 int l = y0 * nx;
00604 for (int j = ymax; j >= ymin; j--) {
00605 int br = l;
00606 remx = 10;
00607 int ll = x0;
00608 for (int i = xmin; i < xsize - 1; i++) {
00609 if (l + ll > lmax || ll >= nx - 2) {
00610 break;
00611 }
00612 int k = 0;
00613 unsigned char p;
00614 if (ll >= nx / 2) {
00615 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00616 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00617 }
00618 else k = nx * ny - (l + 2 * ll) - 2;
00619 if (k>=mid) k-=mid;
00620 else k+=mid;
00621
00622 float t = image_data[k];
00623
00624 int ph;
00625 if (flags&16 && asrgb>2) {
00626 if (l >= (ny * nx - nx)) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00627 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00628 }
00629 if (t <= rm)
00630 p = mingray;
00631 else if (t >= render_max) {
00632 p = maxgray;
00633 }
00634 else if (gamma!=1.0) {
00635 k=(int)(gs2 * (t-render_min));
00636 p = gammamap[k];
00637
00638
00639
00640
00641 }
00642 else {
00643 p = (unsigned char) (gs * (t - render_min));
00644 p += mingray;
00645 }
00646 if (flags&16 && asrgb>2) {
00647 if (ph<256) {
00648 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00649 data[i * asrgb + j * bpl+1] = p*ph/256;
00650 data[i * asrgb + j * bpl+2] = 0;
00651 }
00652 else if (ph<512) {
00653 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00654 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00655 data[i * asrgb + j * bpl] = 0;
00656 }
00657 else {
00658 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00659 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00660 data[i * asrgb + j * bpl+1] = 0;
00661 }
00662 }
00663 else data[i * asrgb + j * bpl] = p;
00664 if (hist) histd[p]++;
00665 ll += addi;
00666 remx += addr;
00667 if (remx > scale_n) {
00668 remx -= scale_n;
00669 ll++;
00670 }
00671 }
00672 l = br + addi * nx;
00673 remy += addr;
00674 if (remy > scale_n) {
00675 remy -= scale_n;
00676 l += nx;
00677 }
00678 }
00679 }
00680 }
00681 else {
00682 if (dsx != -1) {
00683 int l = x0 + y0 * nx;
00684 for (int j = ymax; j >= ymin; j--) {
00685 int br = l;
00686 for (int i = xmin; i < xsize; i++) {
00687 if (l > lmax) {
00688 break;
00689 }
00690 int k = 0;
00691 unsigned char p;
00692 float t;
00693 if (dsx==1) t=image_data[l];
00694 else {
00695 t=0;
00696 for (int iii=0; iii<dsx; iii++) {
00697 for (int jjj=0; jjj<dsy; jjj+=nx) {
00698 t+=image_data[l+iii+jjj];
00699 }
00700 }
00701 t/=dsx*(dsy/nx);
00702 }
00703
00704 if (t <= rm) p = mingray;
00705 else if (t >= render_max) p = maxgray;
00706 else if (gamma!=1.0) {
00707 k=(int)(gs2 * (t-render_min));
00708 p = gammamap[k];
00709
00710
00711
00712
00713 }
00714 else {
00715 if (flags&32){
00716
00717
00718 if (flags&64)
00719 {
00720 p=(unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
00721 }
00722 else{
00723 p=(unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);}
00724
00725 }
00726 else
00727 {
00728 p=(unsigned char) (gs * (t - render_min));
00729 }
00730
00731 p += mingray;
00732 }
00733 data[i * asrgb + j * bpl] = p;
00734 if (hist) histd[p]++;
00735 l += dsx;
00736 }
00737 l = br + dsy;
00738 }
00739 }
00740 else {
00741 remy = 10;
00742 int l = x0 + y0 * nx;
00743 for (int j = ymax; j >= ymin; j--) {
00744 int addj = addi;
00745
00746
00747
00748 if (( l + addi*nx ) >= nxy ) {
00749 addj = (nxy-l)/nx;
00750 if (addj <= 0) continue;
00751 }
00752 int br = l;
00753 remx = 10;
00754 for (int i = xmin; i < xsize; i++) {
00755 if (l > lmax) break;
00756 int k = 0;
00757 unsigned char p;
00758 float t;
00759 if (addi<=1) t = image_data[l];
00760 else {
00761 t=0;
00762 for (int jjj=0; jjj<addj; jjj++) {
00763 for (int iii=0; iii<addi; iii++) {
00764 t+=image_data[l+iii+jjj*nx];
00765 }
00766 }
00767 t/=addi*addi;
00768 }
00769 if (t <= rm) p = mingray;
00770 else if (t >= render_max) p = maxgray;
00771 else if (gamma!=1.0) {
00772 k=(int)(gs2 * (t-render_min));
00773 p = gammamap[k];
00774
00775
00776
00777
00778 }
00779 else {
00780 if (flags&32){
00781
00782 if (flags&64){
00783 p=(unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
00784 }
00785 else{
00786 p=(unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);
00787 }
00788
00789 }
00790 else{
00791 p=(unsigned char) (gs * (t - render_min));
00792 }
00793
00794 p += mingray;
00795 }
00796 data[i * asrgb + j * bpl] = p;
00797 if (hist) histd[p]++;
00798 l += addi;
00799 remx += addr;
00800 if (remx > scale_n) {
00801 remx -= scale_n;
00802 l++;
00803 }
00804 }
00805 l = br + addi * nx;
00806 remy += addr;
00807 if (remy > scale_n) {
00808 remy -= scale_n;
00809 l += nx;
00810 }
00811 }
00812 }
00813 }
00814
00815
00816 if (asrgb==3 && !(flags&16)) {
00817 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00818 for (int i=xmin; i<xsize*3; i+=3) {
00819 data[i+j+1]=data[i+j+2]=data[i+j];
00820 }
00821 }
00822 }
00823 if (asrgb==4 && !(flags&16)) {
00824 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00825 for (int i=xmin; i<xsize*4; i+=4) {
00826 data[i+j+1]=data[i+j+2]=data[i+j+3]=data[i+j];
00827 data[i+j+3]=255;
00828 }
00829 }
00830 }
00831
00832 EXITFUNC;
00833
00834
00835 if (invy) {
00836 int x,y;
00837 char swp;
00838 for (y=0; y<iysize/2; y++) {
00839 for (x=0; x<ixsize; x++) {
00840 swp=ret[y*bpl+x];
00841 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00842 ret[(iysize-y-1)*bpl+x]=swp;
00843 }
00844 }
00845 }
00846
00847
00848 if (flags&16) {
00849 glDrawPixels(ixsize,iysize,GL_LUMINANCE,GL_UNSIGNED_BYTE,(const GLvoid *)ret.data());
00850 }
00851
00852
00853 delete[] grayhe;
00854 delete[] gaussianpdf;
00855 delete[] gaussiancdf;
00856 delete[] gaussianlookup;
00857 delete[] graypdftwo;
00858
00859 return ret;
00860 }
00861
00862
00863 unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z, bool recontour)
00864 {
00865 if ( mc->_isodl != 0 ) glDeleteLists(mc->_isodl,1);
00866
00867 if ( recontour ) mc->calculate_surface();
00868 if ( surface_face_z ) mc->surface_face_z();
00869 if( mc->getRGBmode() ) mc->color_vertices();
00870
00871 #if MARCHING_CUBES_DEBUG
00872 cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() << " points and " << nn.elem() << " normals to render in generate dl" << endl;
00873 #endif
00874 int maxf;
00875
00876
00877
00878 glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00879
00880 #if MARCHING_CUBES_DEBUG
00881 int maxv;
00882 glGetIntegerv(GL_MAX_ELEMENTS_VERTICES,&maxv);
00883 cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
00884 cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
00885 #endif
00886
00887 if ( recontour ) for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00888
00889 if ( maxf % 3 != 0 )
00890 {
00891 maxf = maxf - (maxf%3);
00892 }
00893
00894 if ( tex_id != 0 ) {
00895
00896 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(), 1.0f/(float)mc->_emdata->get_ysize(), 1.0f/(float)mc->_emdata->get_zsize());
00897 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
00898 glDisableClientState(GL_NORMAL_ARRAY);
00899 glTexCoordPointer(3,GL_FLOAT,0,mc->pp.get_data());
00900 }
00901 else {
00902 glEnableClientState(GL_NORMAL_ARRAY);
00903 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00904 glNormalPointer(GL_FLOAT,0, mc->nn.get_data());
00905 }
00906
00907 glEnableClientState(GL_VERTEX_ARRAY);
00908 glVertexPointer(3,GL_FLOAT,0, mc->pp.get_data());
00909
00910 if (mc->getRGBmode()){
00911 glEnableClientState(GL_COLOR_ARRAY);
00912 glColorPointer(3,GL_FLOAT,0, mc->cc.get_data());
00913 }
00914
00915 mc->_isodl = glGenLists(1);
00916
00917 #if MARCHING_CUBES_DEBUG
00918 int time0 = clock();
00919 #endif
00920
00921 glNewList(mc->_isodl,GL_COMPILE);
00922
00923 if ( tex_id != 0 ) {
00924 glEnable(GL_TEXTURE_3D);
00925 glBindTexture(GL_TEXTURE_3D, tex_id);
00926 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
00927 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
00928 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
00929 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
00930 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
00931 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
00932 }
00933
00934
00935
00936 bool drawRange = true;
00937 if ( drawRange == false ) {
00938 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT, mc->ff.get_data());
00939 } else {
00940 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
00941 {
00942 if ( (i+maxf) > mc->ff.elem())
00943 glDrawElements(GL_TRIANGLES,mc->ff.elem()-i,GL_UNSIGNED_INT,&(mc->ff[i]));
00944 else
00945 glDrawElements(GL_TRIANGLES,maxf,GL_UNSIGNED_INT,&(mc->ff[i]));
00946
00947
00948
00949
00950 }
00951 }
00952
00953 if ( tex_id != 0 ) glDisable(GL_TEXTURE_3D);
00954
00955 glEndList();
00956
00957 #if MARCHING_CUBES_DEBUG
00958 int time1 = clock();
00959 cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
00960 #endif
00961 return mc->_isodl;
00962 }
00963
00964 void GLUtil::contour_isosurface(MarchingCubes* mc)
00965 {
00966 mc->calculate_surface();
00967
00968 for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00969
00970 mc->needtobind = true;
00971 }
00972
00973 void GLUtil::render_using_VBOs(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z)
00974 {
00975
00976
00977 #ifdef _WIN32
00978 typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
00979 PFNGLGENBUFFERSPROC glGenBuffers;
00980 glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
00981
00982 typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
00983 PFNGLISBUFFERPROC glIsBuffer;
00984 glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
00985 #endif //_WIN32
00986
00987 if ( surface_face_z ) mc->surface_face_z();
00988
00989
00990 if (!glIsBuffer(mc->buffer[0])) glGenBuffers(4, mc->buffer);
00991
00992
00993 if( mc->getRGBmode() && (mc->rgbgenerator.getNeedToRecolor() || mc->needtobind)){
00994 mc->color_vertices();
00995 mc->needtobind = true;
00996 }
00997
00998 int maxf;
00999 glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
01000
01001 if ( maxf % 3 != 0 )
01002 {
01003 maxf = maxf - (maxf%3);
01004 }
01005
01006 #ifdef _WIN32
01007 typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01008 PFNGLBINDBUFFERPROC glBindBuffer;
01009 glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01010
01011 typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01012 PFNGLBUFFERDATAPROC glBufferData;
01013 glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01014 #endif //_WIN32
01015
01016
01017 glEnableClientState(GL_NORMAL_ARRAY);
01018 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
01019 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[2]);
01020 if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->nn.elem(), mc->nn.get_data(), GL_STATIC_DRAW);
01021 glNormalPointer(GL_FLOAT,0,0);
01022
01023
01024 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[0]);
01025 if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->pp.elem(), mc->pp.get_data(), GL_STATIC_DRAW);
01026 glEnableClientState(GL_VERTEX_ARRAY);
01027 glVertexPointer(3,GL_FLOAT,0,0);
01028
01029
01030 if (mc->getRGBmode()){
01031 glEnableClientState(GL_COLOR_ARRAY);
01032 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[3]);
01033 if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->cc.elem(), mc->cc.get_data(), GL_STATIC_DRAW);
01034 glColorPointer(3,GL_FLOAT,0, 0);
01035 }
01036
01037
01038 glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mc->buffer[1]);
01039 if (mc->needtobind) glBufferData(GL_ELEMENT_ARRAY_BUFFER, 4*mc->ff.elem(), mc->ff.get_data(), GL_STATIC_DRAW);
01040
01041
01042 if (!glIsBuffer(mc->buffer[0])){
01043 cout << "Can't Generate GL Vertex Buffers. glGenBuffer error" << endl;
01044 return;
01045 }
01046
01047
01048 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT,0);
01049
01050 mc->needtobind = false;
01051
01052 }
01053
01054
01055 void GLUtil::glLoadMatrix(const Transform& xform)
01056 {
01057 #ifdef _WIN32
01058 typedef void (APIENTRYP PFNGLLOADTRANSPOSEMATRIXFPROC) (const GLfloat *m);
01059 PFNGLLOADTRANSPOSEMATRIXFPROC glLoadTransposeMatrixf = NULL;
01060 glLoadTransposeMatrixf = (PFNGLLOADTRANSPOSEMATRIXFPROC) wglGetProcAddress("glLoadTransposeMatrixf");
01061 #endif //_WIN32
01062
01063 vector<float> xformlist = xform.get_matrix_4x4();
01064 glLoadTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
01065 }
01066
01067 void GLUtil::glMultMatrix(const Transform& xform)
01068 {
01069 #ifdef _WIN32
01070 typedef void (APIENTRYP PFNGLMULTTRANSPOSEMATRIXFPROC) (const GLfloat *m);
01071 PFNGLMULTTRANSPOSEMATRIXFPROC glMultTransposeMatrixf = NULL;
01072 glMultTransposeMatrixf = (PFNGLMULTTRANSPOSEMATRIXFPROC) wglGetProcAddress("glMultTransposeMatrixf");
01073 #endif //_WIN32
01074
01075 vector<float> xformlist = xform.get_matrix_4x4();
01076 glMultTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
01077 }
01078
01079
01080 void GLUtil::glDrawBoundingBox(float width, float height, float depth)
01081 {
01082
01083 float w2 = width/2.0f;
01084 float h2 = height/2.0f;
01085 float d2 = depth/2.0f;
01086
01087 float vertices[24] = {-w2, h2, d2, w2, h2, d2, w2, -h2, d2, -w2, -h2, d2, -w2, h2, -d2, w2, h2, -d2, w2, -h2, -d2, -w2, -h2, -d2};
01088 int indices[24] = {0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6, 6, 7, 7, 4, 0, 4, 3, 7, 1, 5, 2, 6};
01089
01090 #ifdef _WIN32
01091 typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
01092 PFNGLGENBUFFERSPROC glGenBuffers;
01093 glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
01094
01095 typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
01096 PFNGLISBUFFERPROC glIsBuffer;
01097 glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
01098
01099 typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01100 PFNGLBINDBUFFERPROC glBindBuffer;
01101 glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01102
01103 typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01104 PFNGLBUFFERDATAPROC glBufferData;
01105 glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01106 #endif //_WIN32
01107
01108 if ( glIsBuffer(GLUtil::buffer[0]) == GL_FALSE ){
01109 glGenBuffers(2, GLUtil::buffer);
01110 }
01111
01112
01113 glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
01114 glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
01115 glEnableClientState(GL_VERTEX_ARRAY);
01116 glVertexPointer(3,GL_FLOAT,0,0);
01117
01118 glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, GLUtil::buffer[1]);
01119 glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices, GL_STATIC_DRAW);
01120 glDrawElements(GL_LINES,24,GL_UNSIGNED_INT,0);
01121
01122 }
01123
01124 void GLUtil::glDrawDisk(float radius, int spokes)
01125 {
01126 #ifdef _WIN32
01127 typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
01128 PFNGLGENBUFFERSPROC glGenBuffers;
01129 glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
01130
01131 typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
01132 PFNGLISBUFFERPROC glIsBuffer;
01133 glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
01134
01135 typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01136 PFNGLBINDBUFFERPROC glBindBuffer;
01137 glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01138
01139 typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01140 PFNGLBUFFERDATAPROC glBufferData;
01141 glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01142 #endif //_WIN32
01143
01144
01145 int arraysize = 4*pow((double)2,(double)spokes);
01146 int sideofarray = pow((double)2,(double)spokes);
01147
01148
01149 vector<float> vertices(3*arraysize + 3);
01150
01151
01152 vertices[3*arraysize] = 0.0;
01153 vertices[3*arraysize+1] = 0.0;
01154 vertices[3*arraysize+2] = 0.0;
01155
01156 vertices[0] = 0.0;
01157 vertices[1] = radius;
01158 vertices[2] = 0.0;
01159
01160 vertices[sideofarray*3] = radius;
01161 vertices[sideofarray*3 + 1] = 0.0;
01162 vertices[sideofarray*3 + 2] = 0.0;
01163
01164 vertices[sideofarray*6] = 0.0;
01165 vertices[sideofarray*6 + 1] = -radius;
01166 vertices[sideofarray*6 + 2] = 0.0;
01167
01168 vertices[sideofarray*9] = -radius;
01169 vertices[sideofarray*9 + 1] = 0.0;
01170 vertices[sideofarray*9 + 2] = 0.0;
01171
01172
01173 for(int step = 0; step < spokes; step++){
01174
01175 int x = sideofarray/pow(2.0,(double)(step+1));
01176 for(int i = 1; i <= 4*pow(2.0,(double)step); i++ ){
01177
01178
01179 int index = x + 2*x*(i-1);
01180 int index_f = (index + x) % arraysize;
01181 int index_i = index - x;
01182 cout << index << " " << index_f << " " << index_i << endl;
01183
01184
01185 vertices[index_f*3] = (vertices[index_f*3] - vertices[index_i*3])/2.0f;
01186 vertices[index_f*3 + 1] = (vertices[index_f*3 + 1] - vertices[index_i*3 + 1])/2.0f;
01187 vertices[index_f*3 + 2] = (vertices[index_f*3 + 2] - vertices[index_i*3 + 2])/2.0f;
01188 }
01189 }
01190
01191
01192 if ( glIsBuffer(GLUtil::buffer[0]) == GL_FALSE ){
01193 glGenBuffers(2, GLUtil::buffer);
01194 }
01195
01196
01197 glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
01198
01199 glBufferData(GL_ARRAY_BUFFER, vertices.size(), &(vertices[0]), GL_STATIC_DRAW);
01200 glEnableClientState(GL_VERTEX_ARRAY);
01201 glVertexPointer(3,GL_FLOAT,0,0);
01202
01203
01204
01205
01206 }
01207 #endif // EMAN2_USING_OPENGL