Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

DrawMolItemVolume.C

Go to the documentation of this file.
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: DrawMolItemVolume.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.183 $      $Date: 2021/10/28 21:12:15 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   A continuation of rendering types from DrawMolItem
00019  *
00020  *   This file only contains representations for visualizing volumetric data
00021  ***************************************************************************/
00022 
00023 #include <math.h>
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 
00027 #include "DrawMolItem.h"
00028 #include "DrawMolecule.h"
00029 #include "BaseMolecule.h" // need volume data definitions
00030 #include "Inform.h"
00031 #include "Isosurface.h"
00032 #include "MoleculeList.h"
00033 #include "Scene.h"
00034 #include "VolumetricData.h"
00035 #include "VMDApp.h"
00036 #include "WKFUtils.h"
00037 
00038 #define MYSGN(a) (((a) > 0) ? 1 : -1)
00039 
00040 int DrawMolItem::draw_volume_get_colorid(void) {
00041   int colorid = 0;
00042   int catindex;
00043 
00044   switch (atomColor->method()) {
00045     case AtomColor::MOLECULE:
00046     case AtomColor::COLORID:
00047       colorid = atomColor->color_index();
00048       break;
00049 
00050     default:
00051       catindex = scene->category_index("Display");
00052       colorid = scene->category_item_value(catindex, "Foreground"); 
00053       break;
00054   }
00055 
00056   return colorid;
00057 }
00058 
00059 
00060 void DrawMolItem::draw_volume_box_solid(VolumetricData * v) {
00061   float v0[3], v1[3], v2[3];
00062   float vorigin[3], vxaxis[3], vyaxis[3], vzaxis[3];
00063   int usecolor;
00064 
00065   int i;
00066   for (i=0; i<3; i++) {
00067     vorigin[i] = float(v->origin[i]);
00068     vxaxis[i] = float(v->xaxis[i]);
00069     vyaxis[i] = float(v->yaxis[i]);
00070     vzaxis[i] = float(v->zaxis[i]);
00071   }
00072 
00073   append(DMATERIALON); // enable lighting and shading
00074   usecolor=draw_volume_get_colorid();
00075   cmdColorIndex.putdata(usecolor, cmdList);
00076 
00077   // Draw XY plane
00078   v0[0] = vorigin[0];
00079   v0[1] = vorigin[1];
00080   v0[2] = vorigin[2];
00081 
00082   v1[0] = v0[0] + vxaxis[0]; 
00083   v1[1] = v0[1] + vxaxis[1]; 
00084   v1[2] = v0[2] + vxaxis[2]; 
00085 
00086   v2[0] = v1[0] + vyaxis[0];
00087   v2[1] = v1[1] + vyaxis[1];
00088   v2[2] = v1[2] + vyaxis[2];
00089   cmdSquare.putdata(v2, v1, v0, cmdList);
00090 
00091   // Draw XZ plane
00092   v2[0] = v1[0] + vzaxis[0];
00093   v2[1] = v1[1] + vzaxis[1];
00094   v2[2] = v1[2] + vzaxis[2];
00095   cmdSquare.putdata(v0, v1, v2, cmdList);
00096 
00097   // Draw YZ plane
00098   v1[0] = v0[0] + vyaxis[0]; 
00099   v1[1] = v0[1] + vyaxis[1]; 
00100   v1[2] = v0[2] + vyaxis[2]; 
00101 
00102   v2[0] = v1[0] + vzaxis[0];
00103   v2[1] = v1[1] + vzaxis[1];
00104   v2[2] = v1[2] + vzaxis[2];
00105   cmdSquare.putdata(v2, v1, v0, cmdList);
00106 
00107   // Draw XY +Z plane
00108   v0[0] = vorigin[0] + vzaxis[0];
00109   v0[1] = vorigin[1] + vzaxis[1];
00110   v0[2] = vorigin[2] + vzaxis[2];
00111 
00112   v1[0] = v0[0] + vxaxis[0]; 
00113   v1[1] = v0[1] + vxaxis[1]; 
00114   v1[2] = v0[2] + vxaxis[2]; 
00115 
00116   v2[0] = v1[0] + vyaxis[0];
00117   v2[1] = v1[1] + vyaxis[1];
00118   v2[2] = v1[2] + vyaxis[2];
00119   cmdSquare.putdata(v0, v1, v2, cmdList);
00120 
00121   // Draw XZ +Y plane
00122   v0[0] = vorigin[0] + vyaxis[0];
00123   v0[1] = vorigin[1] + vyaxis[1];
00124   v0[2] = vorigin[2] + vyaxis[2];
00125 
00126   v1[0] = v0[0] + vxaxis[0]; 
00127   v1[1] = v0[1] + vxaxis[1]; 
00128   v1[2] = v0[2] + vxaxis[2]; 
00129 
00130   v2[0] = v1[0] + vzaxis[0];
00131   v2[1] = v1[1] + vzaxis[1];
00132   v2[2] = v1[2] + vzaxis[2];
00133   cmdSquare.putdata(v2, v1, v0, cmdList);
00134 
00135   // Draw YZ +X plane
00136   v0[0] = vorigin[0] + vxaxis[0];
00137   v0[1] = vorigin[1] + vxaxis[1];
00138   v0[2] = vorigin[2] + vxaxis[2];
00139 
00140   v1[0] = v0[0] + vyaxis[0]; 
00141   v1[1] = v0[1] + vyaxis[1]; 
00142   v1[2] = v0[2] + vyaxis[2]; 
00143 
00144   v2[0] = v1[0] + vzaxis[0];
00145   v2[1] = v1[1] + vzaxis[1];
00146   v2[2] = v1[2] + vzaxis[2];
00147   cmdSquare.putdata(v0, v1, v2, cmdList);
00148 
00149   draw_volume_box_lines(v);
00150 }
00151 
00152 void DrawMolItem::draw_volume_box_lines(VolumetricData * v) {
00153   float start[3], end[3];
00154   int xaxiscolor, yaxiscolor, zaxiscolor, axiscolor;
00155   int catindex;
00156   float vorigin[3], vxaxis[3], vyaxis[3], vzaxis[3];
00157 
00158   int i;
00159   for (i=0; i<3; i++) {
00160     vorigin[i] = float(v->origin[i]);
00161     vxaxis[i] = float(v->xaxis[i]);
00162     vyaxis[i] = float(v->yaxis[i]);
00163     vzaxis[i] = float(v->zaxis[i]);
00164   }
00165 
00166   // get colors from the scene pointer
00167   catindex = scene->category_index("Axes");
00168   xaxiscolor = scene->category_item_value(catindex, "X"); 
00169   yaxiscolor = scene->category_item_value(catindex, "Y"); 
00170   zaxiscolor = scene->category_item_value(catindex, "Z"); 
00171 
00172   catindex = scene->category_index("Display");
00173   axiscolor = scene->category_item_value(catindex, "Foreground"); 
00174  
00175   append(DMATERIALOFF);
00176   cmdLineType.putdata(SOLIDLINE, cmdList);
00177   cmdLineWidth.putdata(3, cmdList);
00178 
00179   start[0] = vorigin[0];
00180   start[1] = vorigin[1];
00181   start[2] = vorigin[2];
00182 
00183   // Draw X axis of volume box
00184   cmdColorIndex.putdata(xaxiscolor, cmdList);
00185   end[0] = start[0] + vxaxis[0];
00186   end[1] = start[1] + vxaxis[1];
00187   end[2] = start[2] + vxaxis[2];
00188   cmdLine.putdata(start, end, cmdList);
00189 
00190   // Draw Y axis of volume box
00191   cmdColorIndex.putdata(yaxiscolor, cmdList);
00192   end[0] = start[0] + vyaxis[0];
00193   end[1] = start[1] + vyaxis[1];
00194   end[2] = start[2] + vyaxis[2];
00195   cmdLine.putdata(start, end, cmdList);
00196 
00197   // Draw Z axis of volume box
00198   cmdColorIndex.putdata(zaxiscolor, cmdList);
00199   end[0] = start[0] + vzaxis[0];
00200   end[1] = start[1] + vzaxis[1];
00201   end[2] = start[2] + vzaxis[2];
00202   cmdLine.putdata(start, end, cmdList);
00203 
00204   // Draw remaining outline of volume box
00205   cmdLineWidth.putdata(1, cmdList);
00206   cmdColorIndex.putdata(axiscolor, cmdList);
00207  
00208   start[0] = vorigin[0] + vxaxis[0];    
00209   start[1] = vorigin[1] + vxaxis[1];    
00210   start[2] = vorigin[2] + vxaxis[2];    
00211 
00212   end[0] = start[0] + vyaxis[0];
00213   end[1] = start[1] + vyaxis[1];
00214   end[2] = start[2] + vyaxis[2];
00215   cmdLine.putdata(start, end, cmdList);
00216 
00217   end[0] = start[0] + vzaxis[0];
00218   end[1] = start[1] + vzaxis[1];
00219   end[2] = start[2] + vzaxis[2];
00220   cmdLine.putdata(start, end, cmdList);
00221 
00222   start[0] = vorigin[0] + vyaxis[0];    
00223   start[1] = vorigin[1] + vyaxis[1];    
00224   start[2] = vorigin[2] + vyaxis[2];    
00225 
00226   end[0] = start[0] + vxaxis[0];
00227   end[1] = start[1] + vxaxis[1];
00228   end[2] = start[2] + vxaxis[2];
00229   cmdLine.putdata(start, end, cmdList);
00230 
00231   end[0] = start[0] + vzaxis[0];
00232   end[1] = start[1] + vzaxis[1];
00233   end[2] = start[2] + vzaxis[2];
00234   cmdLine.putdata(start, end, cmdList);
00235 
00236 
00237   start[0] = vorigin[0] + vzaxis[0];    
00238   start[1] = vorigin[1] + vzaxis[1];    
00239   start[2] = vorigin[2] + vzaxis[2];    
00240 
00241   end[0] = start[0] + vxaxis[0];
00242   end[1] = start[1] + vxaxis[1];
00243   end[2] = start[2] + vxaxis[2];
00244   cmdLine.putdata(start, end, cmdList);
00245 
00246   end[0] = start[0] + vyaxis[0];
00247   end[1] = start[1] + vyaxis[1];
00248   end[2] = start[2] + vyaxis[2];
00249   cmdLine.putdata(start, end, cmdList);
00250  
00251 
00252   start[0] = vorigin[0] + vxaxis[0] + vyaxis[0] + vzaxis[0];
00253   start[1] = vorigin[1] + vxaxis[1] + vyaxis[1] + vzaxis[1];
00254   start[2] = vorigin[2] + vxaxis[2] + vyaxis[2] + vzaxis[2];
00255  
00256   end[0] = start[0] - vxaxis[0];
00257   end[1] = start[1] - vxaxis[1];
00258   end[2] = start[2] - vxaxis[2];
00259   cmdLine.putdata(start, end, cmdList);
00260 
00261   end[0] = start[0] - vyaxis[0];
00262   end[1] = start[1] - vyaxis[1];
00263   end[2] = start[2] - vyaxis[2];
00264   cmdLine.putdata(start, end, cmdList);
00265 
00266   end[0] = start[0] - vzaxis[0];
00267   end[1] = start[1] - vzaxis[1];
00268   end[2] = start[2] - vzaxis[2];
00269   cmdLine.putdata(start, end, cmdList);
00270 }
00271 
00272 
00273 void DrawMolItem::draw_volume_isosurface_points(const VolumetricData * v,
00274                                 float isovalue, int stepsize, int thickness) {
00275   int x,y,z;
00276   float xax[3], yax[3], zax[3];
00277   int usecolor;
00278   ResizeArray<float> centers;
00279   ResizeArray<float> colors;
00280 
00281   int i;
00282   float vorigin[3];
00283   for (i=0; i<3; i++) {
00284     vorigin[i] = float(v->origin[i]);
00285   }
00286 
00287   // calculate cell axes
00288   v->cell_axes(xax, yax, zax);
00289   ptrdiff_t zplanesz = v->xsize * v->ysize;
00290  
00291   append(DMATERIALOFF);
00292   usecolor = draw_volume_get_colorid();
00293   const float *cp = scene->color_value(usecolor);
00294   cmdColorIndex.putdata(usecolor, cmdList);
00295 
00296   int pointcount = 0;
00297   for (z=0; z<v->zsize; z+=stepsize) {
00298     // Pre-extend our output buffers so we don't have to call append().
00299     // As long as the output buffers are large enough to hold a full 
00300     // additional plane of points, we can avoid individually appending points.
00301     // Since this resizing scheme bypasses the normal bookkeeping done by
00302     // append() methods, we are responsible for setting the number of 
00303     // occupied array elements via set_size(), such that the results 
00304     // returned by the num() method are correct.
00305     ptrdiff_t newsz = (pointcount + zplanesz) * 3L;
00306     if (newsz >= centers.num()) {
00307       colors.extend(newsz);
00308       centers.extend(newsz);
00309     }
00310 
00311     // set starting position of output pointers
00312     int idx = pointcount * 3;
00313     float *pos=&centers[0];
00314     float *cols=&colors[0];
00315     for (y=0; y<v->ysize; y+=stepsize) {
00316       ptrdiff_t rowidx = z*zplanesz + y*v->xsize;
00317       const float *rowaddr = &(v->data[rowidx]);  
00318 
00319       // compute starting coordinate of row...
00320       float rowx = vorigin[0] + y * yax[0] + z * zax[0];
00321       float rowy = vorigin[1] + y * yax[1] + z * zax[1];
00322       float rowz = vorigin[2] + y * yax[2] + z * zax[2];
00323 
00324       // loop through xsize - 1 rather than the full range
00325       int lastgrid = v->xsize - 1;
00326       for (x=0; x<lastgrid; x+=stepsize) {
00327         // draw a point if the isovalue falls between neighboring X samples
00328         // minimal test criterion: 
00329         //   both terms have to be true, or both false, indicating that
00330         //   the isovalue lies in between the two nearest voxel values
00331         if (!((rowaddr[x] <= isovalue) ^ (isovalue < rowaddr[x+1]))) {
00332           // draw a point there
00333           pos[idx  ] = rowx + x * xax[0];
00334           pos[idx+1] = rowy + x * xax[1];
00335           pos[idx+2] = rowz + x * xax[2];
00336 
00337           // XXX per-vertex colors are redundant here and we could get a 2x 
00338           //     draw perf gain here by eliminating them...
00339           cols[idx  ] = cp[0];
00340           cols[idx+1] = cp[1];
00341           cols[idx+2] = cp[2];
00342 
00343           idx+=3;
00344           pointcount++;
00345         }
00346       } 
00347     } 
00348 
00349     // sync up ResizeArray occupied array size
00350     centers.set_size(idx);
00351     colors.set_size(idx);
00352 
00353     // Avoid integer overflow by submitting geometry batches after 
00354     // every 10M primitives, well before the point of int32 wraparound.
00355     if (pointcount > 10000000) {
00356       cmdPointArray.putdata((float *) &centers[0],
00357                             (float *) &colors[0],
00358                             (float) thickness,
00359                             pointcount,
00360                             cmdList);
00361 
00362       // reset tmp buffers
00363       centers.clear();
00364       colors.clear(); // XXX get rid of per-vertex colors
00365       pointcount=0;
00366     }
00367   }
00368 
00369   if (pointcount > 0) {
00370     cmdPointArray.putdata((float *) &centers[0],
00371                           (float *) &colors[0],
00372                           (float) thickness,
00373                           pointcount,
00374                           cmdList);
00375   }
00376 }
00377 
00378 
00379 void DrawMolItem::draw_volume_isosurface_lit_points(VolumetricData * v, 
00380                                  float isovalue, int stepsize, int thickness) {
00381   int x,y,z;
00382   float xax[3], yax[3], zax[3];
00383   ResizeArray<float> centers;
00384   ResizeArray<float> normals;
00385   ResizeArray<float> colors;
00386 
00387   int i;
00388   float vorigin[3];
00389   for (i=0; i<3; i++) {
00390     vorigin[i] = float(v->origin[i]);
00391   }
00392 
00393   // calculate cell axes
00394   v->cell_axes(xax, yax, zax);
00395   ptrdiff_t zplanesz = v->xsize * v->ysize;
00396 
00397   // get direct access to gradient data for speed, and force
00398   // generation of the gradients if they don't already exist
00399   const float *volgradient = v->access_volume_gradient();
00400 
00401   append(DMATERIALON);
00402   int usecolor = draw_volume_get_colorid();
00403   const float *cp = scene->color_value(usecolor);
00404   cmdColorIndex.putdata(usecolor, cmdList);
00405 
00406   int pointcount = 0;
00407   for (z=0; z<v->zsize; z+=stepsize) {
00408     // Pre-extend our output buffers so we don't have to call append().
00409     // As long as the output buffers are large enough to hold a full
00410     // additional plane of points, we can avoid individually appending points.
00411     // Since this resizing scheme bypasses the normal bookkeeping done by
00412     // append() methods, we are responsible for setting the number of 
00413     // occupied array elements via set_size(), such that the results 
00414     // returned by the num() method are correct.
00415     ptrdiff_t newsz = (pointcount + zplanesz) * 3L;
00416     if ((pointcount + zplanesz)*3L >= centers.num()) {
00417       colors.extend(newsz);
00418       normals.extend(newsz);
00419       centers.extend(newsz);
00420     }
00421 
00422     // set starting position of output pointers
00423     int idx = pointcount * 3;
00424     float *pos=&centers[0];
00425     float *cols=&colors[0];
00426     float *norms=&normals[0]; 
00427     for (y=0; y<v->ysize; y+=stepsize) {
00428       ptrdiff_t rowidx = z*zplanesz + y*v->xsize;
00429       float *rowaddr = &(v->data[rowidx]);  
00430 
00431       // compute starting coordinate of row...
00432       float rowx = vorigin[0] + y * yax[0] + z * zax[0];
00433       float rowy = vorigin[1] + y * yax[1] + z * zax[1];
00434       float rowz = vorigin[2] + y * yax[2] + z * zax[2];
00435 
00436       // loop through xsize - 1 rather than the full range
00437       int lastgrid = v->xsize - 1;
00438       for (x=0; x<lastgrid; x+=stepsize) {
00439         // draw a point if the isovalue falls between neighboring X samples
00440         // minimal test criterion:
00441         //   both terms have to be true, or both false, indicating that
00442         //   the isovalue lies in between the two nearest voxel values
00443         if (!((rowaddr[x] <= isovalue) ^ (isovalue < rowaddr[x+1]))) {
00444           // draw a point there
00445           pos[idx  ] = rowx + x * xax[0];
00446           pos[idx+1] = rowy + x * xax[1];
00447           pos[idx+2] = rowz + x * xax[2];
00448 
00449           ptrdiff_t volgradidx = (rowidx + x) * 3L;
00450           norms[idx  ] = volgradient[volgradidx    ];
00451           norms[idx+1] = volgradient[volgradidx + 1];
00452           norms[idx+2] = volgradient[volgradidx + 2];
00453  
00454           // XXX per-vertex colors are redundant here and we could get a 1.3x 
00455           //     draw perf gain here by eliminating them...
00456           cols[idx  ] = cp[0];
00457           cols[idx+1] = cp[1];
00458           cols[idx+2] = cp[2];
00459 
00460           idx+=3;
00461           pointcount++;
00462         }
00463       } 
00464     } 
00465 
00466     // sync up ResizeArray occupied array size
00467     centers.set_size(idx);
00468     normals.set_size(idx);
00469     colors.set_size(idx);
00470 
00471     // Avoid integer overflow by submitting geometry batches after 
00472     // every 10M primitives, well before the point of int32 wraparound.
00473     if (pointcount > 10000000) {
00474       DispCmdLitPointArray cmdLitPointArray;
00475       cmdLitPointArray.putdata((float *) &centers[0],
00476                                (float *) &normals[0],
00477                                (float *) &colors[0],
00478                                (float) thickness,
00479                                pointcount,
00480                                cmdList);
00481 
00482       // reset tmp buffers
00483       centers.clear();
00484       normals.clear();
00485       colors.clear(); // XXX get rid of per-vertex colors
00486       pointcount=0;
00487     }
00488   }
00489 
00490   if (pointcount > 0) {
00491     DispCmdLitPointArray cmdLitPointArray;
00492     cmdLitPointArray.putdata((float *) &centers[0],
00493                              (float *) &normals[0],
00494                              (float *) &colors[0],
00495                              (float) thickness,
00496                              pointcount,
00497                              cmdList);
00498   }
00499 }
00500 
00501 
00502 void DrawMolItem::draw_volume_isosurface_lines(VolumetricData * v, 
00503                                  float isovalue, int stepsize, int thickness) {
00504   IsoSurface *s = new IsoSurface;
00505 
00506   // XXX We must handle multi-billion triangle isosurface meshes -- needs work
00507   // Now that VMD can load multi-gigavoxel volumes it isn't a safe assumption
00508   // that we can do a single-pass MC isosurface extraction, as this could
00509   // result in billions of triangles, and that creates problems for triangle
00510   // mesh indexing on GPUs that can't do 64-bit indices.  So this code must 
00511   // now decompose the MC isosurface calculation into multiple passes that
00512   // produce less than 2 billion triangles per pass.  Since these are 
00513   // inherently data-parallel, there's also an opportunity to pick up speed
00514   // for large volumes with the use of a persistent CPU thread pool and 
00515   // dicing the problem up into multiple segments that won't break indexing
00516   // and allow us to exploit parallelism.
00517   s->clear();
00518   s->compute(v, isovalue, stepsize);
00519 
00520   if (s->numtriangles > 0) {
00521     append(DMATERIALOFF); // enable lighting and shading
00522     cmdLineType.putdata(SOLIDLINE, cmdList);
00523     cmdLineWidth.putdata(thickness, cmdList);
00524 
00525     int usecolor = draw_volume_get_colorid();
00526     cmdColorIndex.putdata(usecolor, cmdList);
00527 
00528     // draw triangles
00529     ptrdiff_t i;
00530     for (i=0; i<s->numtriangles; i++) {
00531       float *addr;
00532       addr = &(s->v[i * 9]); 
00533       cmdLine.putdata(&addr[0], &addr[3], cmdList);
00534       cmdLine.putdata(&addr[3], &addr[6], cmdList);
00535       cmdLine.putdata(&addr[6], &addr[0], cmdList);
00536     }
00537   }
00538 
00539   delete s; // we don't need this stuff after this point 
00540 }
00541 
00542 
00543 
00544 void DrawMolItem::draw_volume_isosurface_trimesh(VolumetricData * v, 
00545                                                  float isovalue, int stepsize,
00546                                                  const float *voltex) {
00547   IsoSurface s;
00548 
00549   // XXX We must handle multi-billion triangle isosurface meshes -- needs work
00550   // Now that VMD can load multi-gigavoxel volumes it isn't a safe assumption
00551   // that we can do a single-pass MC isosurface extraction, as this could
00552   // result in billions of triangles, and that creates problems for triangle
00553   // mesh indexing on GPUs that can't do 64-bit indices.  So this code must 
00554   // now decompose the MC isosurface calculation into multiple passes that
00555   // produce less than 2 billion triangles per pass.  Since these are 
00556   // inherently data-parallel, there's also an opportunity to pick up speed
00557   // for large volumes with the use of a persistent CPU thread pool and 
00558   // dicing the problem up into multiple segments that won't break indexing
00559   // and allow us to exploit parallelism.
00560   s.clear();                 // initialize isosurface data
00561   s.compute(v, isovalue, stepsize); // compute the isosurface
00562   s.vertexfusion(36, 36);    // identify and eliminate duplicated vertices
00563   s.normalize();             // normalize interpolated gradient/surface normals
00564 
00565 #if 1
00566   if (s.numtriangles > 0) {
00567     append(DMATERIALON); // enable lighting and shading
00568     int usecolor = draw_volume_get_colorid();
00569     cmdColorIndex.putdata(usecolor, cmdList);
00570 
00571     if (voltex != NULL) {
00572       // assign per-vertex colors by a 3-D texture map
00573       s.set_color_voltex_rgb3fv(voltex); 
00574     } else {
00575       // use a single color for the entire mesh
00576       s.set_color_rgb3fv(scene->color_value(usecolor)); 
00577     }
00578 
00579     // Create a triangle mesh
00580     // XXX don't try to stripify it since this triggers a crash in ACTC for
00581     //     unknown reasons
00582     // XXX ACTC is undoubtably unable to deal with billion-triangle meshes
00583     // XXX The DispCmd code won't yet scale past 2 billion triangle meshes
00584     cmdTriMesh.putdata(&s.v[0], &s.n[0], &s.c[0], s.v.num() / 3,
00585                        &s.f[0], s.numtriangles, 0, cmdList);
00586   }
00587 #else
00588   if (s.numtriangles > 0) {
00589     append(DMATERIALON); // enable lighting and shading
00590     int usecolor = draw_volume_get_colorid();
00591     cmdColorIndex.putdata(usecolor, cmdList);
00592 
00593     // draw surface with per-vertex normals using a vertex array
00594     float *c = new float[s.numtriangles * 9L];
00595     const float *fp = scene->color_value(usecolor);
00596     ptrdiff_t i;
00597     for (i=0; i<s.numtriangles; i++) { 
00598       int ind = i * 9;
00599 
00600       c[ind    ] = fp[0]; // Red
00601       c[ind + 1] = fp[1]; // Green
00602       c[ind + 2] = fp[2]; // Blue
00603 
00604       ind+=3;
00605       c[ind    ] = fp[0]; // Red
00606       c[ind + 1] = fp[1]; // Green
00607       c[ind + 2] = fp[2]; // Blue
00608 
00609       ind+=3;
00610       c[ind    ] = fp[0]; // Red
00611       c[ind + 1] = fp[1]; // Green
00612       c[ind + 2] = fp[2]; // Blue
00613     }
00614 
00615     // Create a triangle mesh
00616     // XXX don't try to stripify it since this triggers a crash in ACTC for
00617     //     unknown reasons
00618     // XXX ACTC is undoubtably unable to deal with billion-triangle meshes
00619     // XXX The DispCmd code won't yet scale past 2 billion triangle meshes
00620     cmdTriMesh.putdata(&s.v[0], &s.n[0], c, s.v.num() / 3,
00621                        &s.f[0], s.numtriangles, 0, cmdList);
00622 
00623     delete [] c;
00624   }
00625 #endif
00626 }
00627 
00628 
00629 // recompute sliceTextureCoord, sliceVertexCoord, and sliceNormal
00630 // using the given volumetric data set, axis, and offset
00631 // sliceAxis and sliceOffset.  
00632 static void prepare_texture_coordinates(const VolumetricData *v, 
00633     float loc, int axis, float *sliceNormal, float *sliceTextureCoords,
00634     float *sliceVertexes) {
00635 
00636   float t0[3], t1[3], t2[3], t3[3];
00637   float v0[3], v1[3], v2[3], v3[3];
00638   float normal[3];
00639 
00640   float vorigin[3], vxaxis[3], vyaxis[3], vzaxis[3];
00641   int i;
00642   for (i=0; i<3; i++) {
00643     vorigin[i] = float(v->origin[i]);
00644     vxaxis[i] = float(v->xaxis[i]);
00645     vyaxis[i] = float(v->yaxis[i]);
00646     vzaxis[i] = float(v->zaxis[i]);
00647   }
00648 
00649 
00650   if (loc < 0.0f)
00651       loc = 0.0f;
00652 
00653   if (loc > 1.0f)
00654       loc = 1.0f;
00655 
00656   switch (axis) {
00657     // X-Axis
00658     case 0:
00659     default:
00660       t0[0] = loc;
00661       t0[1] = 0.0f;
00662       t0[2] = 0.0f;
00663 
00664       t1[0] = loc;
00665       t1[1] = 0.0f;
00666       t1[2] = 1.0f;
00667 
00668       t2[0] = loc;
00669       t2[1] = 1.0f;
00670       t2[2] = 1.0f;
00671 
00672       t3[0] = loc;
00673       t3[1] = 1.0f;
00674       t3[2] = 0.0f;
00675 
00676       v0[0] = vorigin[0] + (vxaxis[0] * loc);
00677       v0[1] = vorigin[1] + (vxaxis[1] * loc);
00678       v0[2] = vorigin[2] + (vxaxis[2] * loc);
00679 
00680       v1[0] = v0[0] + vzaxis[0];
00681       v1[1] = v0[1] + vzaxis[1];
00682       v1[2] = v0[2] + vzaxis[2];
00683 
00684       v2[0] = v0[0] + vzaxis[0] + vyaxis[0];
00685       v2[1] = v0[1] + vzaxis[1] + vyaxis[1];
00686       v2[2] = v0[2] + vzaxis[2] + vyaxis[2];
00687 
00688       v3[0] = v0[0]             + vyaxis[0];
00689       v3[1] = v0[1]             + vyaxis[1];
00690       v3[2] = v0[2]             + vyaxis[2];
00691 
00692       normal[0] = vxaxis[0];
00693       normal[1] = vxaxis[1];
00694       normal[2] = vxaxis[2];
00695       vec_normalize(&normal[0]);
00696       break;
00697 
00698     // Y-Axis
00699     case 1:
00700       t0[0] = 0.0f;
00701       t0[1] = loc;
00702       t0[2] = 0.0f;
00703 
00704       t1[0] = 1.0f;
00705       t1[1] = loc;
00706       t1[2] = 0.0f;
00707 
00708       t2[0] = 1.0f;
00709       t2[1] = loc;
00710       t2[2] = 1.0f;
00711 
00712       t3[0] = 0.0f;
00713       t3[1] = loc;
00714       t3[2] = 1.0f;
00715 
00716       v0[0] = vorigin[0] + (vyaxis[0] * loc);
00717       v0[1] = vorigin[1] + (vyaxis[1] * loc);
00718       v0[2] = vorigin[2] + (vyaxis[2] * loc);
00719 
00720       v1[0] = v0[0] + vxaxis[0];
00721       v1[1] = v0[1] + vxaxis[1];
00722       v1[2] = v0[2] + vxaxis[2];
00723 
00724       v2[0] = v0[0] + vxaxis[0] + vzaxis[0];
00725       v2[1] = v0[1] + vxaxis[1] + vzaxis[1];
00726       v2[2] = v0[2] + vxaxis[2] + vzaxis[2];
00727 
00728       v3[0] = v0[0]             + vzaxis[0];
00729       v3[1] = v0[1]             + vzaxis[1];
00730       v3[2] = v0[2]             + vzaxis[2];
00731 
00732       normal[0] = vyaxis[0];
00733       normal[1] = vyaxis[1];
00734       normal[2] = vyaxis[2];
00735       vec_normalize(&normal[0]);
00736       break;
00737 
00738     // Z-Axis
00739     case 2:
00740       t0[0] = 0.0f;
00741       t0[1] = 0.0f;
00742       t0[2] = loc;
00743 
00744       t1[0] = 1.0f;
00745       t1[1] = 0.0f;
00746       t1[2] = loc;
00747 
00748       t2[0] = 1.0f;
00749       t2[1] = 1.0f;
00750       t2[2] = loc;
00751 
00752       t3[0] = 0.0f;
00753       t3[1] = 1.0f;
00754       t3[2] = loc;
00755 
00756       v0[0] = vorigin[0] + (vzaxis[0] * loc);
00757       v0[1] = vorigin[1] + (vzaxis[1] * loc);
00758       v0[2] = vorigin[2] + (vzaxis[2] * loc);
00759 
00760       v1[0] = v0[0] + vxaxis[0];
00761       v1[1] = v0[1] + vxaxis[1];
00762       v1[2] = v0[2] + vxaxis[2];
00763 
00764       v2[0] = v0[0] + vxaxis[0] + vyaxis[0];
00765       v2[1] = v0[1] + vxaxis[1] + vyaxis[1];
00766       v2[2] = v0[2] + vxaxis[2] + vyaxis[2];
00767 
00768       v3[0] = v0[0]             + vyaxis[0];
00769       v3[1] = v0[1]             + vyaxis[1];
00770       v3[2] = v0[2]             + vyaxis[2];
00771 
00772       normal[0] = vzaxis[0];
00773       normal[1] = vzaxis[1];
00774       normal[2] = vzaxis[2];
00775       vec_normalize(&normal[0]);
00776       break;
00777   }
00778 
00779   vec_copy(sliceTextureCoords  , t0);
00780   vec_copy(sliceTextureCoords+3, t1);
00781   vec_copy(sliceTextureCoords+6, t2);
00782   vec_copy(sliceTextureCoords+9, t3);
00783   vec_copy(sliceVertexes  , v0);
00784   vec_copy(sliceVertexes+3, v1);
00785   vec_copy(sliceVertexes+6, v2);
00786   vec_copy(sliceVertexes+9, v3);
00787   vec_copy(sliceNormal, normal);
00788 }
00789 
00790 void DrawMolItem::updateVolumeTexture() {
00791   float vmin, vmax;
00792   int volid = atomColor->volume_index();
00793   if (atomRep->method()==AtomRep::VOLSLICE) {
00794     // the volslice rep has its own volid specification
00795     volid = (int)atomRep->get_data(AtomRep::SPHERERES);
00796   }
00797 
00798   VolumetricData *v = mol->modify_volume_data(volid);
00799   if (v == NULL) {
00800     msgInfo << "No volume data loaded at index " << volid << sendmsg;
00801     return;
00802   } 
00803 
00804   // determine if a new 3D texture needs to be prepared and uploaded.  This
00805   // occurs if:
00806   // (1) the coloring method has changed
00807   // (2) the choice of volumetric data set has changed.
00808   atomColor->get_colorscale_minmax(&vmin, &vmax);
00809   if (!vmin && !vmax) {
00810     v->datarange(vmin, vmax);
00811   }
00812 
00813   if (volumeTexture.getTextureMap() && !(needRegenerate & COL_REGEN) &&
00814       volid == voltexVolid && atomColor->method() == voltexColorMethod &&
00815       voltexDataMin == vmin && voltexDataMax == vmax) {
00816     // nothing to do
00817     return;
00818   }
00819 
00820   voltexColorMethod = atomColor->method();
00821   voltexVolid = volid;
00822   voltexDataMin = vmin; 
00823   voltexDataMax = vmax;
00824 
00825   volumeTexture.setGridData(v);
00826 
00827   // update the volumeTexture instance
00828   switch (atomColor->method()) {
00829     case AtomColor::POS:
00830     case AtomColor::POSX:
00831     case AtomColor::POSY:
00832     case AtomColor::POSZ:
00833       volumeTexture.generatePosTexture();
00834       break;
00835 
00836     case AtomColor::INDEX:
00837       volumeTexture.generateIndexTexture();
00838       break;
00839 
00840     case AtomColor::CHARGE:
00841       volumeTexture.generateChargeTexture(vmin, vmax);
00842       break;
00843 
00844     case AtomColor::VOLUME:
00845       volumeTexture.generateColorScaleTexture(vmin, vmax, scene);
00846       break;
00847 
00848 #if 0
00849     case AtomColor::STRUCTURE:
00850       // 3-D Contour lines
00851       volumeTexture.generateContourLineTexture(0.5, 0.5);
00852       break;
00853 #endif
00854       
00855     case AtomColor::NAME:
00856     default:
00857       // HSV color ramp
00858       volumeTexture.generateHSVTexture(vmin, vmax);
00859       break;
00860   }
00861 }
00862 
00863 
00864 void DrawMolItem::draw_volslice(int volid, float slice, int axis, int texmode) {
00865   const VolumetricData *v = mol->modify_volume_data(volid);
00866   float sliceNormal[3];
00867   float sliceTextureCoords[12];
00868   float sliceVertexes[12];
00869   if (!volumeTexture.getTextureMap()) {
00870     msgErr << "draw_volslice: no texture map has been generated!" << sendmsg;
00871     return;
00872   }
00873 
00874   sprintf(commentBuffer, "Mol[%d] Rep[%d] VolSlice", mol->id(), repNumber);
00875   cmdCommentX.putdata(commentBuffer, cmdList);
00876 
00877   prepare_texture_coordinates(v, slice, axis, sliceNormal, sliceTextureCoords,
00878                               sliceVertexes);
00879   
00880   // Rescale the texture coordinates so that they include just the 
00881   // part of the texture map to which we mapped our volume data.
00882   // Clamp range to just below 1.0 to prevent
00883   // textures from being clamped to the
00884   // border color on some video cards.
00885   float tscale[3] = { float(v->xsize), float(v->ysize), float(v->zsize) };
00886   const int *size = volumeTexture.getTextureSize();
00887   for (int i=0; i<3; i++) {
00888     float rescale = (tscale[i] / (float)size[i]) * 0.99999f;
00889     sliceTextureCoords[i  ] *= rescale;
00890     sliceTextureCoords[i+3] *= rescale;
00891     sliceTextureCoords[i+6] *= rescale;
00892     sliceTextureCoords[i+9] *= rescale;
00893   }
00894   // add command to draw the slice itself.
00895   append(DMATERIALON); // enable lighting and shading
00896 
00897   // Pass instructions for the slice itself.
00898   cmdVolSlice.putdata(texmode, sliceNormal, sliceVertexes, sliceTextureCoords, cmdList);
00899 }
00900 
00901 
00902 void DrawMolItem::draw_isosurface(int volid, float isovalue, int drawbox, int style, int stepsize, int thickness) {
00903   VolumetricData * v = NULL;
00904 
00905   v = mol->modify_volume_data(volid);
00906   if (v == NULL) {
00907     msgInfo << "No volume data loaded at index " << volid << sendmsg;
00908     return;
00909   } 
00910 
00911   sprintf(commentBuffer, "Mol[%d] Rep[%d] Isosurface", mol->id(), repNumber);
00912   cmdCommentX.putdata(commentBuffer, cmdList);
00913 
00914   // Safety checks to prevent stepsize from cratering surface extraction code
00915   if (stepsize >= v->xsize)
00916     stepsize = (v->xsize - 1);
00917   if (stepsize >= v->ysize)
00918     stepsize = (v->ysize - 1);
00919   if (stepsize >= v->zsize)
00920     stepsize = (v->zsize - 1);
00921   if (stepsize < 2)
00922     stepsize = 1;
00923 
00924   if (drawbox > 0) {
00925     // don't texture the box if color by volume is active
00926     if (atomColor->method() == AtomColor::VOLUME) {
00927       append(DVOLTEXOFF);
00928     }
00929     // wireframe only?  or solid?
00930     if (style > 0 || drawbox == 2) {
00931       draw_volume_box_lines(v);
00932     } else {
00933       draw_volume_box_solid(v);
00934     }
00935     if (atomColor->method() == AtomColor::VOLUME) {
00936       append(DVOLTEXON);
00937     }
00938   } 
00939 
00940   if ((drawbox == 2) || (drawbox == 0)) {
00941     switch (style) {
00942       case 3:
00943         // shaded points isosurface looping over X-axis, 1 point max per voxel
00944         draw_volume_isosurface_lit_points(v, isovalue, stepsize, thickness);
00945         break;
00946 
00947       case 2:
00948         // points isosurface looping over X-axis, 1 point max per voxel
00949         draw_volume_isosurface_points(v, isovalue, stepsize, thickness);
00950         break;
00951 
00952       case 1:
00953         // lines implementation, max of 18 line per voxel (3-per triangle)
00954         draw_volume_isosurface_lines(v, isovalue, stepsize, thickness);
00955         break;
00956 
00957       case 0:
00958       default:
00959         // trimesh polygonalized surface, max of 6 triangles per voxel
00960         // XXX to enable per-vertex colors rather than 3-D texture mapping,
00961         //     we need to expose a setting for this purpose, for use when
00962         //     people need to export to packages incapable of 3-D texturing.
00963         draw_volume_isosurface_trimesh(v, isovalue, stepsize);
00964         break;
00965     }
00966   }
00967 }
00968 
00969 
00970 #if 0
00971 // draw a 2-D contour line in a slice plane
00972 void DrawMolItem::draw_volslice_contour_lines(int volid, float slice, int axis) {
00973   const VolumetricData * v = NULL;
00974 
00975   v = mol->get_volume_data(volid);
00976 
00977   if (v == NULL) {
00978     msgInfo << "No volume data loaded at index " << volid << sendmsg;
00979     return;
00980   }
00981 
00982   msgInfo << "Placeholder for contour slice not implemented yet..."
00983 }
00984 #endif
00985 
00986 
00987 // calculate seed voxels for field lines based on gradient magnitude
00988 int DrawMolItem::calcseeds_grid(VolumetricData * v, ResizeArray<float> *seeds, int maxseedcount) {
00989   int i;
00990   float vorigin[3];
00991   for (i=0; i<3; i++) {
00992     vorigin[i] = float(v->origin[i]);
00993   }
00994 
00995   // calculate cell axes
00996   float xax[3], yax[3], zax[3];
00997   v->cell_axes(xax, yax, zax);
00998 
00999   int seedcount = maxseedcount+1; // force loop to run once
01000   int stepsize = 1;
01001 
01002   // get direct access to gradient data for speed, and force
01003   // generation of the gradients if they don't already exist
01004   const float *volgradient = v->access_volume_gradient();
01005 
01006   // iterate if we generate more seeds than we can really use
01007   while (seedcount > maxseedcount) {
01008     seedcount=0;
01009     seeds->clear();
01010 
01011     int x,y,z;
01012     float pos[3];
01013     ptrdiff_t zplanesz = v->xsize*v->ysize;
01014     for (z=0; z<v->zsize && seedcount < maxseedcount; z+=stepsize) {
01015       for (y=0; y<v->ysize; y+=stepsize) {
01016         ptrdiff_t rowidx = z*zplanesz + y*v->xsize;
01017 
01018         // compute starting coordinate of row...
01019         float rowx = vorigin[0] + y * yax[0] + z * zax[0];
01020         float rowy = vorigin[1] + y * yax[1] + z * zax[1];
01021         float rowz = vorigin[2] + y * yax[2] + z * zax[2];
01022 
01023         for (x=0; x<v->xsize; x+=stepsize) {
01024           float grad[3];
01025           vec_copy(grad, &volgradient[(rowidx + x) * 3L]);
01026 
01027           pos[0] = rowx + x * xax[0];
01028           pos[1] = rowy + x * xax[1];
01029           pos[2] = rowz + x * xax[2];
01030   
01031           seedcount++;
01032           seeds->append3(&pos[0]);
01033         }
01034       }
01035     }
01036     stepsize++; // increase stepsize if we have to try again
01037   }
01038   
01039   return seedcount;
01040 }
01041 
01042 
01043 // calculate seed voxels for field lines based on gradient magnitude
01044 int DrawMolItem::calcseeds_gradient_magnitude(VolumetricData * v, ResizeArray<float> *seeds, float seedmin, float seedmax, int maxseedcount) {
01045   float seedmin2 = seedmin*seedmin;
01046   float seedmax2 = seedmax*seedmax;
01047 
01048   int i;
01049   float vorigin[3];
01050   for (i=0; i<3; i++) {
01051     vorigin[i] = float(v->origin[i]);
01052   }
01053 
01054   // calculate cell axes
01055   float xax[3], yax[3], zax[3];
01056   v->cell_axes(xax, yax, zax);
01057 
01058   int seedcount = maxseedcount+1; // force loop to run once
01059   int stepsize = 1;
01060 
01061   // get direct access to gradient data for speed, and force
01062   // generation of the gradients if they don't already exist
01063   const float *volgradient = v->access_volume_gradient();
01064 
01065   // iterate if we generate more seeds than we can really use
01066   while (seedcount > maxseedcount) {
01067     seedcount=0;
01068     seeds->clear();
01069 
01070     int x,y,z;
01071     float pos[3];
01072     ptrdiff_t zplanesz = v->xsize*v->ysize;
01073     for (z=0; z<v->zsize && seedcount < maxseedcount; z+=stepsize) {
01074       for (y=0; y<v->ysize; y+=stepsize) {
01075         ptrdiff_t rowidx = z*zplanesz + y*v->xsize;
01076 
01077         // compute starting coordinate of row...
01078         float rowx = vorigin[0] + y * yax[0] + z * zax[0];
01079         float rowy = vorigin[1] + y * yax[1] + z * zax[1];
01080         float rowz = vorigin[2] + y * yax[2] + z * zax[2];
01081 
01082         for (x=0; x<v->xsize; x+=stepsize) {
01083           float grad[3];
01084           float gradmag2;
01085           vec_copy(grad, &volgradient[(rowidx + x) * 3L]);
01086           gradmag2 = dot_prod(grad, grad);
01087 
01088           if ((gradmag2 <= seedmax2) &&
01089               (gradmag2 >= seedmin2)) {
01090             pos[0] = rowx + x * xax[0];
01091             pos[1] = rowy + x * xax[1];
01092             pos[2] = rowz + x * xax[2];
01093   
01094             seedcount++;
01095             seeds->append3(&pos[0]);
01096           }
01097         }
01098       }
01099     }
01100     stepsize++; // increase stepsize if we have to try again
01101   }
01102   
01103   return seedcount;
01104 }
01105 
01106 
01107 // draw a 3-D field lines that follow the volume gradient
01108 void DrawMolItem::draw_volume_field_lines(int volid, int seedusegrid, int maxseeds, 
01109                                           float seedval, float deltacell, 
01110                                           float minlen, float maxlen, 
01111                                           int drawstyle, int tuberes, float thickness) {
01112   VolumetricData * v = NULL;
01113   v = mol->modify_volume_data(volid);
01114   int printdonemesg=0;
01115 
01116   if (v == NULL) {
01117     msgInfo << "No volume data loaded at index " << volid << sendmsg;
01118     return;
01119   }
01120 
01121   int seedcount = 0;
01122   int pointcount = 0;
01123   int totalpointcount = 0;
01124   int usecolor;
01125   ResizeArray<float> seeds;
01126 
01127   sprintf(commentBuffer,"Mol[%d] Rep[%d] FieldLines", mol->id(), repNumber);
01128   cmdCommentX.putdata(commentBuffer, cmdList);
01129 
01130 
01131   DispCmdSphereRes cmdSphereRes;
01132   if (drawstyle != 0) {
01133     cmdSphereRes.putdata(tuberes, cmdList);
01134     append(DMATERIALON); // enable lighting and shading
01135     thickness *= 0.05f;  // XXX hack until we have a better GUI
01136   } else {
01137     append(DMATERIALOFF);
01138   }
01139 
01140   usecolor = draw_volume_get_colorid();
01141   cmdColorIndex.putdata(usecolor, cmdList);
01142 
01143   if (seedusegrid)
01144     seedcount = calcseeds_grid(v, &seeds, maxseeds);
01145   else
01146     seedcount = calcseeds_gradient_magnitude(v, &seeds, seedval*0.5f, seedval*1.5f, maxseeds);
01147 
01148   // Integrate field lines starting with each of the seeds to simulate
01149   // particle advection.
01150   // Uses Euler's approximation for solving the initial value problem.
01151   // We could get a more accurate solution using a fourth order Runge-Kutta
01152   // method, but with more math per iteration.  We may want to implement 
01153   // the integrator as a user selected option.
01154 
01155   // The choice of integration step size is currently arbitrary,
01156   // but will become a user-defined parameter, since it affects speed
01157   // and accuracy.  A good default might be 0.25 times the smallest
01158   // grid cell spacing axis.
01159   float lx, ly, lz;
01160   v->cell_lengths(&lx, &ly, &lz);
01161   float mincelllen=lx;
01162   mincelllen = (mincelllen < ly) ? mincelllen : ly;
01163   mincelllen = (mincelllen < lz) ? mincelllen : lz;
01164   float delta=mincelllen * deltacell; // delta per step (compensates gradient magnitude)
01165 
01166   // minimum gradient magnitude, before we consider that we've found
01167   // a critical point in the dataset.
01168   float mingmag =  0.0001f;
01169 
01170   // max gradient magnitude, before we consider it a source/sink
01171   float maxgmag = 5;
01172 
01173   ResizeArray<float> points;
01174   
01175   // ensure that the volume gradient has been computed prior to
01176   // sampling it for field line construction... (discard pointer)
01177   v->access_volume_gradient();
01178 
01179   // For each seed point, integrate in both positive and
01180   // negative directions for a field line length up to
01181   // the maxlen criterion.
01182   wkfmsgtimer *msgt = wkf_msg_timer_create(1);
01183   int seed;
01184   for (seed=0; seed < seedcount; seed++) {
01185     // emit UI messages as integrator runs, for long calculations...
01186     if (!(seed & 7) && wkf_msg_timer_timeout(msgt)) {
01187       char tmpbuf[128] = { 0 };
01188       sprintf(tmpbuf, "%6.2f %% complete", (100.0f * seed) / (float) seedcount);
01189       msgInfo << "integrating " << seedcount << " field lines: " << tmpbuf << sendmsg;
01190       printdonemesg=1;
01191     }
01192  
01193     int direction;
01194     for (direction=-1; direction != 1; direction=1) {
01195       float pos[3], comsum[3];
01196       vec_copy(pos, &seeds[seed*3]); // integration starting point is the seed
01197 
01198       // init the arrays
01199       points.clear();
01200 
01201       // main integration loop
01202       pointcount=0;
01203       totalpointcount++;
01204       float len=0;
01205       int iterations=0;
01206       float dir = (float) direction;
01207 
01208       vec_zero(comsum); // clear center of mass accumulator
01209 
01210       while ((len<maxlen) && (totalpointcount < 100000000)) {
01211         float grad[3];
01212 
01213         // sample gradient at the current position
01214         v->voxel_gradient_interpolate_from_coord(pos, grad);
01215 
01216         // Early-exit if we run out of bounds (gradient returned will
01217         // be a vector of NANs), run into a critical point (zero gradient)
01218         // or a huge gradient at a source/sink point in the dataset.
01219         // Since IEEE FP defines that all tests against NaN should return
01220         // false, we invert the the bound tests with so that a comparison
01221         // with a gmag of NaN will early-exit the advection loop.
01222         float gmag = norm(grad);
01223         if (!(gmag >= mingmag && gmag <= maxgmag))
01224            break;
01225 
01226         // Draw the current point only after the gradient value
01227         // has been checked, so we don't end up with out-of-bounds
01228         // vertices.
01229         // Only emit a fraction of integration points for display since
01230         // the integrator stepsize needs to be small for more numerical
01231         // accuracy, but the field lines themselves can be well 
01232         // represented with fewer sample points.
01233         if (!(iterations & 1)) {
01234           // Add a vertex for this field line
01235           points.append3(&pos[0]);
01236 
01237           vec_incr(comsum, pos);
01238  
01239           pointcount++;
01240           totalpointcount++;
01241         }
01242 
01243         // adjust integration stepsize so we never move more than 
01244         // the distance specified by delta at each step, to compensate
01245         // for varying gradient magnitude
01246         vec_scaled_add(pos, dir * delta / gmag, grad); // integrate position
01247         len += delta; // accumulate distance
01248 
01249         iterations++;
01250       }
01251 
01252       int drawfieldline = 1;
01253 
01254       // only draw the field line for this seed if we have enough points.
01255       // If we haven't reached the minimum field line length, we'll
01256       // drop the whole field line.
01257       if (pointcount < 2 || len < minlen)
01258         drawfieldline = 0;
01259 
01260       // only draw if bounding sphere diameter exceeds minlen
01261       if (drawfieldline) {       
01262         float com[3];
01263         vec_scale(com, 1.0f / (float) pointcount, comsum);
01264         float minlen2 = minlen*minlen;
01265 
01266         drawfieldline = 0;
01267         int p;
01268         for (p=0; p<pointcount; p++) {
01269           if ((2.0f * distance2(com, &points[p*3])) > minlen2) {
01270             drawfieldline = 1;
01271             break;
01272           }
01273         }
01274       }
01275 
01276       // if drawing style is tubes or spheres we enable the alternate
01277       // rendering path
01278       if (drawstyle != 0) {
01279         // only draw the field line if it met all selection criteria
01280         if (drawfieldline) {
01281           cmdColorIndex.putdata(usecolor, cmdList);
01282           DispCmdCylinder cmdCyl;
01283           int maxcylidx = (pointcount - 1) * 3;
01284           int p;
01285           if (drawstyle == 1) {
01286             for (p=0; p<maxcylidx; p+=3) {
01287               cmdCyl.putdata(&points[p], &points[p+3], thickness,
01288                              tuberes, 0, cmdList);
01289             }
01290             maxcylidx++;
01291           }
01292 
01293           for (p=0; p<maxcylidx; p+=3) {
01294             cmdSphere.putdata(&points[p], thickness, cmdList);
01295           }
01296         }
01297       } else if (drawfieldline) {
01298         // only draw the field line if it met all selection criteria
01299         cmdLineType.putdata(SOLIDLINE, cmdList);
01300         cmdLineWidth.putdata((int) thickness, cmdList);
01301         cmdColorIndex.putdata(usecolor, cmdList);
01302         DispCmdPolyLineArray cmdPolyLineArray;
01303         cmdPolyLineArray.putdata(&points[0], pointcount, cmdList);
01304       } 
01305     }
01306   }
01307   wkf_msg_timer_destroy(msgt);
01308 
01309   if (printdonemesg) 
01310     msgInfo << "field line integration complete." << sendmsg;
01311 }
01312 
01313 

Generated on Wed Dec 4 02:43:38 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002