00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: VolumeTexture.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.29 $ $Date: 2020年10月22日 03:43:24 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * Class for managing volumetric texture maps for use by the various 00019 * DrawMolItem representation methods. 00020 ***************************************************************************/ 00021 00022 #include "VolumeTexture.h" 00023 #include "VolumetricData.h" 00024 #include "AtomColor.h" 00025 #include "Scene.h" 00026 #include "VMDApp.h" 00027 #include "Inform.h" 00028 #include "WKFUtils.h" 00029 00030 #include <math.h> 00031 #include <stdio.h> 00032 00033 VolumeTexture::VolumeTexture() 00034 : v(NULL), texmap(NULL), texid(0) { 00035 size[0] = size[1] = size[2] = 0; 00036 } 00037 00038 00039 VolumeTexture::~VolumeTexture() { 00040 if (texmap) vmd_dealloc(texmap); 00041 } 00042 00043 00044 void VolumeTexture::setGridData(VolumetricData *voldata) { 00045 if (v == voldata) return; 00046 v = voldata; 00047 if (texmap) { 00048 vmd_dealloc(texmap); 00049 texmap = NULL; 00050 } 00051 size[0] = size[1] = size[2] = 0; 00052 texid = 0; 00053 } 00054 00055 00056 int VolumeTexture::allocateTextureMap(ptrdiff_t ntexels) { 00057 texid = 0; 00058 if (texmap) { 00059 vmd_dealloc(texmap); 00060 texmap = NULL; 00061 } 00062 ptrdiff_t sz = ntexels*3L*sizeof(unsigned char); 00063 texmap = (unsigned char *) vmd_alloc(sz); 00064 if (texmap == NULL) { 00065 msgErr << "Texture map allocation failed, out of memory?" << sendmsg; 00066 msgErr << "Texture map texel count: " << ntexels << sendmsg; 00067 msgErr << "Failed allocation of size: " << sz / (1024L*1024L) 00068 << "MB" <<sendmsg; 00069 return FALSE; 00070 } 00071 memset(texmap, 0, ntexels*3L*sizeof(unsigned char)); 00072 texid = VMDApp::get_texserialnum(); 00073 return TRUE; 00074 } 00075 00076 00077 void VolumeTexture::generatePosTexture() { 00078 // nice small texture that will work everywhere 00079 size[0] = size[1] = size[2] = 32; 00080 int x, y, z; 00081 ptrdiff_t addr, addr2; 00082 ptrdiff_t num = num_texels(); 00083 00084 if (!allocateTextureMap(num)) return; 00085 00086 ptrdiff_t planesz = size[0] * size[1]; 00087 for (z=0; z<size[2]; z++) { 00088 for (y=0; y<size[1]; y++) { 00089 addr = z * planesz + y * size[0]; 00090 for (x=0; x<size[0]; x++) { 00091 addr2 = (addr + x) * 3L; 00092 texmap[addr2 ] = (unsigned char) (((float) x / (float) size[0]) * 255.0f); 00093 texmap[addr2 + 1] = (unsigned char) (((float) y / (float) size[1]) * 255.0f); 00094 texmap[addr2 + 2] = (unsigned char) (((float) z / (float) size[2]) * 255.0f); 00095 } 00096 } 00097 } 00098 } 00099 00100 00101 // convert Hue/Saturation/Value to RGB colors 00102 static void HSItoRGB(float h, float s, float i, 00103 unsigned char *r, unsigned char *g, unsigned char *b) { 00104 float rv, gv, bv, t; 00105 00106 t=float(VMD_TWOPI)*h; 00107 rv=(float) (1 + s*sinf(t - float(VMD_TWOPI)/3.0f)); 00108 gv=(float) (1 + s*sinf(t)); 00109 bv=(float) (1 + s*sinf(t + float(VMD_TWOPI)/3.0f)); 00110 00111 t=(float) (254.9999 * i / 2); 00112 00113 *r=(unsigned char)(rv*t); 00114 *g=(unsigned char)(gv*t); 00115 *b=(unsigned char)(bv*t); 00116 } 00117 00118 00119 // return the smallest power of two greater than size, up to 2^16. 00120 static int nextpower2(int size) { 00121 int i; 00122 int power; 00123 00124 if (size == 1) 00125 return 1; 00126 00127 power=1; 00128 for (i=0; i<16; i++) { 00129 power <<= 1; 00130 if (power >= size) 00131 return power; 00132 } 00133 00134 return power; 00135 } 00136 00137 00138 void VolumeTexture::generateIndexTexture() { 00139 // nice small texture that will work everywhere 00140 size[0] = size[1] = size[2] = 32; 00141 int x, y, z; 00142 ptrdiff_t addr, addr2, addr3, index; 00143 ptrdiff_t num = num_texels(); 00144 unsigned char coltable[3 * 4096]; 00145 00146 if (!allocateTextureMap(num)) return; 00147 00148 // build a fast color lookup table 00149 for (index=0; index<4096; index++) { 00150 addr = index * 3L; 00151 HSItoRGB(8.0f * index / 4096.0f, 0.75, 1.0, 00152 coltable+addr, coltable+addr+1, coltable+addr+2); 00153 } 00154 00155 ptrdiff_t planesz = size[0] * size[1]; 00156 for (z=0; z<size[2]; z++) { 00157 for (y=0; y<size[1]; y++) { 00158 addr = z * planesz + y * size[0]; 00159 for (x=0; x<size[0]; x++) { 00160 index = addr + x; 00161 addr2 = index * 3L; 00162 addr3 = ((int) ((index / (float) num) * 4095)) * 3L; 00163 texmap[addr2 ] = coltable[addr3 ]; 00164 texmap[addr2 + 1] = coltable[addr3 + 1]; 00165 texmap[addr2 + 2] = coltable[addr3 + 2]; 00166 } 00167 } 00168 } 00169 } 00170 00171 00172 void VolumeTexture::generateChargeTexture(float vmin, float vmax) { 00173 // need a volumetric dataset for this 00174 if (!v) return; 00175 00176 int x, y, z; 00177 ptrdiff_t addr, addr2; 00178 ptrdiff_t daddr; 00179 float vscale, vrange; 00180 00181 size[0] = v->xsize; 00182 size[1] = v->ysize; 00183 size[2] = v->zsize; 00184 for (int i=0; i<3; i++) { 00185 size[i] = nextpower2(size[i]); 00186 } 00187 ptrdiff_t num = num_texels(); 00188 if (!allocateTextureMap(num)) return; 00189 00190 vrange = vmax - vmin; 00191 if (fabs(vrange) < 0.00001) 00192 vscale = 0.0f; 00193 else 00194 vscale = 1.00001f / vrange; 00195 00196 // map volume data scalars to colors 00197 ptrdiff_t planesz = size[0] * size[1]; 00198 ptrdiff_t vplanesz = v->xsize * v->ysize; 00199 for (z=0; z<v->zsize; z++) { 00200 for (y=0; y<v->ysize; y++) { 00201 addr = z * planesz + y * size[0]; 00202 daddr = z * vplanesz + y * v->xsize; 00203 for (x=0; x<v->xsize; x++) { 00204 addr2 = (addr + x) * 3L; 00205 float level, r, g, b; 00206 00207 // map data to range 0->1 00208 level = (v->data[daddr + x] - vmin) * vscale; 00209 level = level < 0 ? 0 : 00210 level > 1 ? 1 : level; 00211 00212 // low values are mapped to red, high to blue 00213 r = (1.0f - level) * 255.0f; 00214 b = level * 255.0f; 00215 if (level < 0.5f) { 00216 g = level * 2.0f * 128.0f; 00217 } else { 00218 g = (0.5f - (level - 0.5f)) * 2.0f * 128.0f; 00219 } 00220 00221 texmap[addr2 ] = (unsigned char) r; 00222 texmap[addr2 + 1] = (unsigned char) g; 00223 texmap[addr2 + 2] = (unsigned char) b; 00224 } 00225 } 00226 } 00227 } 00228 00229 00230 void VolumeTexture::generateHSVTexture(float vmin, float vmax) { 00231 int x, y, z; 00232 ptrdiff_t index, addr, addr2, addr3; 00233 ptrdiff_t daddr; 00234 float vscale, vrange; 00235 unsigned char coltable[3 * 4096]; 00236 00237 size[0] = v->xsize; 00238 size[1] = v->ysize; 00239 size[2] = v->zsize; 00240 for (int i=0; i<3; i++) { 00241 size[i] = nextpower2(size[i]); 00242 } 00243 ptrdiff_t num = num_texels(); 00244 if (!allocateTextureMap(num)) return; 00245 00246 // build a fast color lookup table 00247 for (index=0; index<4096; index++) { 00248 addr = index * 3; 00249 HSItoRGB(4.0f * index / 4096.0f, 0.75, 1.0, 00250 coltable+addr, coltable+addr+1, coltable+addr+2); 00251 } 00252 00253 // calculate scaling factors 00254 vrange = vmax - vmin; 00255 if (fabs(vrange) < 0.00001) 00256 vscale = 0.0f; 00257 else 00258 vscale = 1.00001f / vrange; 00259 00260 // map volume data scalars to colors 00261 ptrdiff_t planesz = size[0] * size[1]; 00262 ptrdiff_t vplanesz = v->xsize * v->ysize; 00263 for (z=0; z<v->zsize; z++) { 00264 for (y=0; y<v->ysize; y++) { 00265 addr = z * planesz + y * size[0]; 00266 daddr = z * vplanesz + y * v->xsize; 00267 for (x=0; x<v->xsize; x++) { 00268 addr2 = (addr + x) * 3L; 00269 float level; 00270 00271 // map data to range 0->1 00272 level = (v->data[daddr + x] - vmin) * vscale; 00273 00274 // Conditional range test written in terms of inclusion 00275 // within the 0:1 range so that cases that encounter a level 00276 // value of NaN will be clamped to 0 rather than remaining NaN 00277 // and subsequently calculating an out-of-bounds color map address. 00278 // Writing the range test this way works because IEEE FP is defined 00279 // such that all comparisons vs. NaN return false. 00280 level = (level >= 0) ? ((level <= 1) ? level : 1) : 0; 00281 00282 // map values to an HSV color map 00283 addr3 = ((int) (level * 4095)) * 3L; 00284 texmap[addr2 ] = coltable[addr3 ]; 00285 texmap[addr2 + 1] = coltable[addr3 + 1]; 00286 texmap[addr2 + 2] = coltable[addr3 + 2]; 00287 } 00288 } 00289 } 00290 } 00291 00292 00293 void VolumeTexture::generateColorScaleTexture(float vmin, float vmax, const Scene *scene) { 00294 int i, x, y, z; 00295 ptrdiff_t addr; 00296 ptrdiff_t daddr; 00297 float vscale, vrange; 00298 00299 wkf_timerhandle timer=wkf_timer_create(); 00300 wkf_timer_start(timer); 00301 00302 size[0] = v->xsize; 00303 size[1] = v->ysize; 00304 size[2] = v->zsize; 00305 for (i=0; i<3; i++) { 00306 size[i] = nextpower2(size[i]); 00307 } 00308 ptrdiff_t num = num_texels(); 00309 if (!allocateTextureMap(num)) return; 00310 00311 vrange = vmax - vmin; 00312 if (fabs(vrange) < 0.00001) 00313 vscale = 0.0f; 00314 else 00315 vscale = 1.00001f / vrange; 00316 00317 // build an unsigned byte RGB color table to speed up conversion 00318 // by eliminating floating point to unsigned char conversions in the 00319 // innermost loops over voxels 00320 unsigned char * rgb3ub = (unsigned char *) calloc(1, MAPCLRS * 3 * sizeof(unsigned char)); 00321 for (i=0; i<MAPCLRS; i++) { 00322 const float *rgb = scene->color_value(MAPCOLOR(i)); 00323 int idx = i*3; 00324 rgb3ub[idx + 0] = (unsigned char)(rgb[0]*255.0f); 00325 rgb3ub[idx + 1] = (unsigned char)(rgb[1]*255.0f); 00326 rgb3ub[idx + 2] = (unsigned char)(rgb[2]*255.0f); 00327 } 00328 00329 // map volume data scalars to colors 00330 ptrdiff_t planesz = size[0] * size[1]; 00331 ptrdiff_t vplanesz = v->xsize * v->ysize; 00332 const int maxcolorindex = MAPCLRS-1; 00333 const float colscale = vscale * maxcolorindex; 00334 for (z=0; z<v->zsize; z++) { 00335 for (y=0; y<v->ysize; y++) { 00336 addr = (z * planesz + y * size[0]) * 3L; 00337 daddr = z * vplanesz + y * v->xsize; 00338 const float *data = &v->data[daddr]; // simplify inner loop ptr arith 00339 unsigned char *tmaprow = &texmap[addr]; // simplify inner loop ptr arith 00340 for (x=0; x<v->xsize; x++,tmaprow+=3) { 00341 // map data min/max to range 0->1 00342 // values must be clamped before use, since user-specified 00343 // min/max can cause out-of-range color indices to be generated 00344 int colindex = (int) ((data[x] - vmin) * colscale); 00345 00346 // This code isn't vulnerable to the effects of NaN inputs 00347 // because the value clamping logic is performed in the integer 00348 // domain, so regardless what comes out of the colindex calculation, 00349 // it will be clamped to the legal color index range. 00350 00351 // clamp integer color index to 0..maxcolorindex 00352 #if 1 00353 colindex = (colindex < 0) ? 0 : colindex; 00354 colindex = (colindex > maxcolorindex) ? maxcolorindex : colindex; 00355 #elif 1 00356 int mask = 0 - (colindex > maxcolorindex); 00357 colindex = (maxcolorindex & mask) | (colindex & ~mask); 00358 #else 00359 if (colindex < 0) 00360 colindex = 0; 00361 else if (colindex > maxcolorindex) 00362 colindex = maxcolorindex; 00363 #endif 00364 00365 colindex *= 3; // compute offset in color table 00366 tmaprow[0] = rgb3ub[colindex ]; 00367 tmaprow[1] = rgb3ub[colindex + 1]; 00368 tmaprow[2] = rgb3ub[colindex + 2]; 00369 } 00370 } 00371 } 00372 00373 free(rgb3ub); 00374 00375 wkf_timer_stop(timer); 00376 double t = wkf_timer_time(timer); 00377 if (t > 5.0) 00378 msgInfo << "Texture generation time: " << t << " sec." << sendmsg; 00379 wkf_timer_destroy(timer); 00380 } 00381 00382 00383 void VolumeTexture::generateContourLineTexture(float densityperline, float linewidth) { 00384 int x, y, z; 00385 ptrdiff_t addr, addr2; 00386 float xp, yp, zp; 00387 00388 float datamin, datamax; 00389 v->datarange(datamin, datamax); 00390 printf("Contour lines...\n"); 00391 printf("range / densityperline: %f\n", log(datamax - datamin) / densityperline); 00392 00393 size[0] = nextpower2(v->xsize*2); 00394 size[1] = nextpower2(v->ysize*2); 00395 size[2] = nextpower2(v->zsize*2); 00396 ptrdiff_t num = num_texels(); 00397 if (!allocateTextureMap(num)) return; 00398 00399 // map volume data scalars to contour line colors 00400 ptrdiff_t planesz = size[0] * size[1]; 00401 for (z=0; z<size[2]; z++) { 00402 zp = ((float) z / size[2]) * v->zsize; 00403 for (y=0; y<size[1]; y++) { 00404 addr = z * planesz + y * size[0]; 00405 yp = ((float) y / size[1]) * v->ysize; 00406 for (x=0; x<size[0]; x++) { 00407 addr2 = (addr + x) * 3L; 00408 xp = ((float) x / size[0]) * v->xsize; 00409 float level; 00410 00411 level = float(fmod(log(v->voxel_value_interpolate(xp,yp,zp)), densityperline) / densityperline); 00412 00413 if (level < linewidth) { 00414 texmap[addr2 ] = 0; 00415 texmap[addr2 + 1] = 0; 00416 texmap[addr2 + 2] = 0; 00417 } else { 00418 texmap[addr2 ] = 255; 00419 texmap[addr2 + 1] = 255; 00420 texmap[addr2 + 2] = 255; 00421 } 00422 } 00423 } 00424 } 00425 } 00426 00427 00428 void VolumeTexture::calculateTexgenPlanes(float v0[3], float v1[3], float v2[3], float v3[3]) const { 00429 00430 int i; 00431 if (!texmap || !v) { 00432 // do something sensible 00433 vec_zero(v0); 00434 vec_zero(v1); 00435 vec_zero(v2); 00436 vec_zero(v3); 00437 v1[0] = v2[1] = v3[2] = 1; 00438 return; 00439 } 00440 00441 // rescale texture coordinates by the portion of the 00442 // entire texture volume they reference 00443 // XXX added an additional scale factor to keep "nearest" texture modes 00444 // rounding into the populated area rather than catching black 00445 // texels in the empty part of the texture volume 00446 float tscale[3]; 00447 tscale[0] = (v->xsize / (float)size[0]) * 0.99999f; 00448 tscale[1] = (v->ysize / (float)size[1]) * 0.99999f; 00449 tscale[2] = (v->zsize / (float)size[2]) * 0.99999f; 00450 00451 // calculate length squared of volume axes 00452 float lensq[3]; 00453 vec_zero(lensq); 00454 for (i=0; i<3; i++) { 00455 lensq[0] += float(v->xaxis[i] * v->xaxis[i]); 00456 lensq[1] += float(v->yaxis[i] * v->yaxis[i]); 00457 lensq[2] += float(v->zaxis[i] * v->zaxis[i]); 00458 } 00459 00460 // Calculate reciprocal space lattice vectors, which are used 00461 // in the OpenGL texgen eye space plane equations in order to transform 00462 // incoming world coordinates to the correct texture coordinates. 00463 // This code should work for both orthogonal and non-orthogonal volumes. 00464 // The last step adds in the NPOT texture scaling where necessary. 00465 // Reference: Introductory Solid State Physics, H.P.Myers, page 43 00466 float xaxdir[3], yaxdir[3], zaxdir[3]; 00467 float nxaxdir[3], nyaxdir[3], nzaxdir[3]; 00468 float bxc[3], cxa[3], axb[3]; 00469 float tmp; 00470 00471 // copy axis direction vectors 00472 for (i=0; i<3; i++) { 00473 xaxdir[i] = float(v->xaxis[i]); 00474 yaxdir[i] = float(v->yaxis[i]); 00475 zaxdir[i] = float(v->zaxis[i]); 00476 } 00477 00478 // calculate reciprocal lattice vector for X texture coordiante 00479 cross_prod(bxc, yaxdir, zaxdir); 00480 tmp = dot_prod(xaxdir, bxc); 00481 for (i=0; i<3; i++) { 00482 nxaxdir[i] = bxc[i] / tmp; 00483 } 00484 00485 // calculate reciprocal lattice vector for Y texture coordiante 00486 cross_prod(cxa, zaxdir, xaxdir); 00487 tmp = dot_prod(yaxdir, cxa); 00488 for (i=0; i<3; i++) { 00489 nyaxdir[i] = cxa[i] / tmp; 00490 } 00491 00492 // calculate reciprocal lattice vector for Z texture coordiante 00493 cross_prod(axb, xaxdir, yaxdir); 00494 tmp = dot_prod(zaxdir, axb); 00495 for (i=0; i<3; i++) { 00496 nzaxdir[i] = axb[i] / tmp; 00497 } 00498 00499 // negate and transform the volume origin to reciprocal space 00500 // for use in the OpenGL texgen plane equation 00501 float norigin[3]; 00502 for (i=0; i<3; i++) 00503 norigin[i] = float(v->origin[i]); 00504 00505 v0[0] = -dot_prod(norigin, nxaxdir) * tscale[0]; 00506 v0[1] = -dot_prod(norigin, nyaxdir) * tscale[1]; 00507 v0[2] = -dot_prod(norigin, nzaxdir) * tscale[2]; 00508 00509 // scale the volume axes for the OpenGL texgen plane equation 00510 for (i=0; i<3; i++) { 00511 v1[i] = nxaxdir[i] * tscale[0]; 00512 v2[i] = nyaxdir[i] * tscale[1]; 00513 v3[i] = nzaxdir[i] * tscale[2]; 00514 } 00515 } 00516 00517 00518