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: Voltool.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.11 $ $Date: 2020年10月15日 17:52:06 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * General volumetric data processing routines, particularly supporting MDFF 00019 * 00020 ***************************************************************************/ 00021 00022 #include <stdio.h> 00023 #include <stdlib.h> 00024 #include <string.h> 00025 #include <float.h> // FLT_MAX etc 00026 #include "Inform.h" 00027 #include "utilities.h" 00028 //#include "SymbolTable.h" 00029 00030 #include "AtomSel.h" 00031 #include "VMDApp.h" 00032 #include "MoleculeList.h" 00033 #include "VolumetricData.h" 00034 #include "VolMapCreate.h" // volmap_write_dx_file() 00035 00036 #include "CUDAMDFF.h" 00037 #include "MDFF.h" 00038 #include <math.h> 00039 #include <tcl.h> 00040 #include "TclCommands.h" 00041 #include "Measure.h" 00042 #include "MolFilePlugin.h" 00043 00044 #include <iostream> 00045 #include <string> 00046 #include <sstream> 00047 #include "Voltool.h" 00048 00051 void init_from_intersection(VolumetricData *mapA, const VolumetricData *mapB, VolumetricData *newvol) { 00052 int d; 00053 00054 // Find intersection of A and B 00055 // The following has been verified for orthog. cells 00056 // (Does not work for non-orthog cells) 00057 00058 for (d=0; d<3; d++) { 00059 newvol->origin[d] = MAX(mapA->origin[d], mapB->origin[d]); 00060 newvol->xaxis[d] = MAX(MIN(mapA->origin[d]+mapA->xaxis[d], mapB->origin[d]+mapB->xaxis[d]), newvol->origin[d]); 00061 newvol->yaxis[d] = MAX(MIN(mapA->origin[d]+mapA->yaxis[d], mapB->origin[d]+mapB->yaxis[d]), newvol->origin[d]); 00062 newvol->zaxis[d] = MAX(MIN(mapA->origin[d]+mapA->zaxis[d], mapB->origin[d]+mapB->zaxis[d]), newvol->origin[d]); 00063 } 00064 00065 vec_sub(newvol->xaxis, newvol->xaxis, newvol->origin); 00066 vec_sub(newvol->yaxis, newvol->yaxis, newvol->origin); 00067 vec_sub(newvol->zaxis, newvol->zaxis, newvol->origin); 00068 00069 newvol->xsize = (int) MAX(dot_prod(newvol->xaxis,mapA->xaxis)*mapA->xsize/dot_prod(mapA->xaxis,mapA->xaxis), \ 00070 dot_prod(newvol->xaxis,mapB->xaxis)*mapB->xsize/dot_prod(mapB->xaxis,mapB->xaxis)); 00071 newvol->ysize = (int) MAX(dot_prod(newvol->yaxis,mapA->yaxis)*mapA->ysize/dot_prod(mapA->yaxis,mapA->yaxis), \ 00072 dot_prod(newvol->yaxis,mapB->yaxis)*mapB->ysize/dot_prod(mapB->yaxis,mapB->yaxis)); 00073 newvol->zsize = (int) MAX(dot_prod(newvol->zaxis,mapA->zaxis)*mapA->zsize/dot_prod(mapA->zaxis,mapA->zaxis), \ 00074 dot_prod(newvol->zaxis,mapB->zaxis)*mapB->zsize/dot_prod(mapB->zaxis,mapB->zaxis)); 00075 // Create map... 00076 if (newvol->data) delete[] newvol->data; 00077 newvol->data = new float[newvol->gridsize()]; 00078 } 00079 00082 void init_from_intersection(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol) { 00083 int d; 00084 00085 // Find intersection of A and B 00086 // The following has been verified for orthog. cells 00087 // (Does not work for non-orthog cells) 00088 00089 for (d=0; d<3; d++) { 00090 newvol->origin[d] = MAX(mapA->origin[d], mapB->origin[d]); 00091 newvol->xaxis[d] = MAX(MIN(mapA->origin[d]+mapA->xaxis[d], mapB->origin[d]+mapB->xaxis[d]), newvol->origin[d]); 00092 newvol->yaxis[d] = MAX(MIN(mapA->origin[d]+mapA->yaxis[d], mapB->origin[d]+mapB->yaxis[d]), newvol->origin[d]); 00093 newvol->zaxis[d] = MAX(MIN(mapA->origin[d]+mapA->zaxis[d], mapB->origin[d]+mapB->zaxis[d]), newvol->origin[d]); 00094 } 00095 00096 vec_sub(newvol->xaxis, newvol->xaxis, newvol->origin); 00097 vec_sub(newvol->yaxis, newvol->yaxis, newvol->origin); 00098 vec_sub(newvol->zaxis, newvol->zaxis, newvol->origin); 00099 00100 newvol->xsize = (int) MAX(dot_prod(newvol->xaxis,mapA->xaxis)*mapA->xsize/dot_prod(mapA->xaxis,mapA->xaxis), \ 00101 dot_prod(newvol->xaxis,mapB->xaxis)*mapB->xsize/dot_prod(mapB->xaxis,mapB->xaxis)); 00102 newvol->ysize = (int) MAX(dot_prod(newvol->yaxis,mapA->yaxis)*mapA->ysize/dot_prod(mapA->yaxis,mapA->yaxis), \ 00103 dot_prod(newvol->yaxis,mapB->yaxis)*mapB->ysize/dot_prod(mapB->yaxis,mapB->yaxis)); 00104 newvol->zsize = (int) MAX(dot_prod(newvol->zaxis,mapA->zaxis)*mapA->zsize/dot_prod(mapA->zaxis,mapA->zaxis), \ 00105 dot_prod(newvol->zaxis,mapB->zaxis)*mapB->zsize/dot_prod(mapB->zaxis,mapB->zaxis)); 00106 // Create map... 00107 if (newvol->data) delete[] newvol->data; 00108 newvol->data = new float[newvol->gridsize()]; 00109 } 00110 00111 00114 void init_from_union(VolumetricData *mapA, const VolumetricData *mapB, VolumetricData *newvol) { 00115 // Find union of A and B 00116 // The following has been verified for orthog. cells 00117 // (Does not work for non-orthog cells) 00118 00119 vec_zero(newvol->xaxis); 00120 vec_zero(newvol->yaxis); 00121 vec_zero(newvol->zaxis); 00122 00123 int d; 00124 for (d=0; d<3; d++) { 00125 newvol->origin[d] = MIN(mapA->origin[d], mapB->origin[d]); 00126 } 00127 d=0; 00128 newvol->xaxis[d] = MAX(MAX(mapA->origin[d]+mapA->xaxis[d], mapB->origin[d]+mapB->xaxis[d]), newvol->origin[d]); 00129 d=1; 00130 newvol->yaxis[d] = MAX(MAX(mapA->origin[d]+mapA->yaxis[d], mapB->origin[d]+mapB->yaxis[d]), newvol->origin[d]); 00131 d=2; 00132 newvol->zaxis[d] = MAX(MAX(mapA->origin[d]+mapA->zaxis[d], mapB->origin[d]+mapB->zaxis[d]), newvol->origin[d]); 00133 00134 newvol->xaxis[0] -= newvol->origin[0]; 00135 newvol->yaxis[1] -= newvol->origin[1]; 00136 newvol->zaxis[2] -= newvol->origin[2]; 00137 00138 newvol->xsize = (int) MAX(dot_prod(newvol->xaxis,mapA->xaxis)*mapA->xsize/dot_prod(mapA->xaxis,mapA->xaxis), \ 00139 dot_prod(newvol->xaxis,mapB->xaxis)*mapB->xsize/dot_prod(mapB->xaxis,mapB->xaxis)); 00140 newvol->ysize = (int) MAX(dot_prod(newvol->yaxis,mapA->yaxis)*mapA->ysize/dot_prod(mapA->yaxis,mapA->yaxis), \ 00141 dot_prod(newvol->yaxis,mapB->yaxis)*mapB->ysize/dot_prod(mapB->yaxis,mapB->yaxis)); 00142 newvol->zsize = (int) MAX(dot_prod(newvol->zaxis,mapA->zaxis)*mapA->zsize/dot_prod(mapA->zaxis,mapA->zaxis), \ 00143 dot_prod(newvol->zaxis,mapB->zaxis)*mapB->zsize/dot_prod(mapB->zaxis,mapB->zaxis)); 00144 // Create map... 00145 if (newvol->data) delete[] newvol->data; 00146 newvol->data = new float[newvol->gridsize()]; 00147 } 00148 00149 void init_from_identity(VolumetricData *mapA, VolumetricData *newvol) { 00150 00151 vec_copy(newvol->origin, mapA->origin); 00152 vec_copy(newvol->xaxis, mapA->xaxis); 00153 vec_copy(newvol->yaxis, mapA->yaxis); 00154 vec_copy(newvol->zaxis, mapA->zaxis); 00155 00156 newvol->xsize = mapA->xsize; 00157 newvol->ysize = mapA->ysize; 00158 newvol->zsize = mapA->zsize; 00159 // Create map... 00160 if (newvol->data) delete[] newvol->data; 00161 newvol->data = new float[newvol->gridsize()]; 00162 } 00163 00164 00167 void init_from_union(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol) { 00168 // Find union of A and B 00169 // The following has been verified for orthog. cells 00170 // (Does not work for non-orthog cells) 00171 00172 vec_zero(newvol->xaxis); 00173 vec_zero(newvol->yaxis); 00174 vec_zero(newvol->zaxis); 00175 00176 int d; 00177 for (d=0; d<3; d++) { 00178 newvol->origin[d] = MIN(mapA->origin[d], mapB->origin[d]); 00179 } 00180 d=0; 00181 newvol->xaxis[d] = MAX(MAX(mapA->origin[d]+mapA->xaxis[d], mapB->origin[d]+mapB->xaxis[d]), newvol->origin[d]); 00182 d=1; 00183 newvol->yaxis[d] = MAX(MAX(mapA->origin[d]+mapA->yaxis[d], mapB->origin[d]+mapB->yaxis[d]), newvol->origin[d]); 00184 d=2; 00185 newvol->zaxis[d] = MAX(MAX(mapA->origin[d]+mapA->zaxis[d], mapB->origin[d]+mapB->zaxis[d]), newvol->origin[d]); 00186 00187 newvol->xaxis[0] -= newvol->origin[0]; 00188 newvol->yaxis[1] -= newvol->origin[1]; 00189 newvol->zaxis[2] -= newvol->origin[2]; 00190 00191 newvol->xsize = (int) MAX(dot_prod(newvol->xaxis,mapA->xaxis)*mapA->xsize/dot_prod(mapA->xaxis,mapA->xaxis), \ 00192 dot_prod(newvol->xaxis,mapB->xaxis)*mapB->xsize/dot_prod(mapB->xaxis,mapB->xaxis)); 00193 newvol->ysize = (int) MAX(dot_prod(newvol->yaxis,mapA->yaxis)*mapA->ysize/dot_prod(mapA->yaxis,mapA->yaxis), \ 00194 dot_prod(newvol->yaxis,mapB->yaxis)*mapB->ysize/dot_prod(mapB->yaxis,mapB->yaxis)); 00195 newvol->zsize = (int) MAX(dot_prod(newvol->zaxis,mapA->zaxis)*mapA->zsize/dot_prod(mapA->zaxis,mapA->zaxis), \ 00196 dot_prod(newvol->zaxis,mapB->zaxis)*mapB->zsize/dot_prod(mapB->zaxis,mapB->zaxis)); 00197 // Create map... 00198 if (newvol->data) delete[] newvol->data; 00199 newvol->data = new float[newvol->gridsize()]; 00200 } 00201 00202 VolumetricData * init_new_volume(){ 00203 double origin[3] = {0., 0., 0.}; 00204 double xaxis[3] = {0., 0., 0.}; 00205 double yaxis[3] = {0., 0., 0.}; 00206 double zaxis[3] = {0., 0., 0.}; 00207 int numvoxels [3] = {0, 0, 0}; 00208 float *data = NULL; 00209 VolumetricData *newvol = new VolumetricData("density map", origin, xaxis, yaxis, zaxis, 00210 numvoxels[0], numvoxels[1], numvoxels[2], 00211 data); 00212 return newvol; 00213 } 00214 00215 00216 int init_new_volume_molecule(VMDApp *app, VolumetricData *newvol, const char *name){ 00217 int newvolmolid = app->molecule_new(name,0,1); 00218 00219 app->molecule_add_volumetric(newvolmolid, "density newvol", newvol->origin, 00220 newvol->xaxis, newvol->yaxis, newvol->zaxis, newvol->xsize, newvol->ysize, 00221 newvol->zsize, newvol->data); 00222 app->molecule_set_style("Isosurface"); 00223 app->molecule_addrep(newvolmolid); 00224 00225 return newvolmolid; 00226 } 00227 00228 00229 void vol_com(VolumetricData *vol, float *com){ 00230 float ix,iy,iz; 00231 00232 vec_zero(com); 00233 double mass = 0.0; 00234 00235 for (int i = 0; i < vol->gridsize(); i++) { 00236 float m = vol->data[i]; 00237 vol->voxel_coord(i, ix, iy, iz); 00238 float coord[3] = {ix,iy,iz}; 00239 mass = mass+m; 00240 vec_scale(coord, m, coord); 00241 vec_add(com, com, coord); 00242 00243 } 00244 00245 float scale = 1.0f/float(mass); 00246 vec_scale(com, scale, com); 00247 } 00248 00249 void add(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol, bool interp, bool USE_UNION) { 00250 00251 // adding maps by spatial coords is slower than doing it directly, but 00252 // allows for precisely subtracting unaligned maps, and/or maps of 00253 // different resolutions 00254 00255 if ( USE_UNION) { 00256 // UNION VERSION 00257 init_from_union(mapA, mapB, newvol); 00258 } else { 00259 // INTERSECTION VERSION 00260 init_from_intersection(mapA, mapB, newvol); 00261 } 00262 for (long i=0; i<newvol->gridsize(); i++){ 00263 float x, y, z; 00264 newvol->voxel_coord(i, x, y, z); 00265 00266 if (interp) newvol->data[i] = \ 00267 mapA->voxel_value_interpolate_from_coord_safe(x,y,z) + \ 00268 mapB->voxel_value_interpolate_from_coord_safe(x,y,z); 00269 else newvol->data[i] = \ 00270 mapA->voxel_value_from_coord_safe(x,y,z) + \ 00271 mapB->voxel_value_from_coord_safe(x,y,z); 00272 } 00273 00274 } 00275 00276 void subtract(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol, bool interp, bool USE_UNION) { 00277 00278 // adding maps by spatial coords is slower than doing it directly, but 00279 // allows for precisely subtracting unaligned maps, and/or maps of 00280 // different resolutions 00281 00282 if ( USE_UNION) { 00283 // UNION VERSION 00284 init_from_union(mapA, mapB, newvol); 00285 } else { 00286 // INTERSECTION VERSION 00287 init_from_intersection(mapA, mapB, newvol); 00288 } 00289 for (long i=0; i<newvol->gridsize(); i++){ 00290 float x, y, z; 00291 newvol->voxel_coord(i, x, y, z); 00292 00293 if (interp) newvol->data[i] = \ 00294 mapA->voxel_value_interpolate_from_coord_safe(x,y,z) - \ 00295 mapB->voxel_value_interpolate_from_coord_safe(x,y,z); 00296 else newvol->data[i] = \ 00297 mapA->voxel_value_from_coord_safe(x,y,z) - \ 00298 mapB->voxel_value_from_coord_safe(x,y,z); 00299 } 00300 00301 } 00302 00303 void multiply(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol, bool interp, bool USE_UNION) { 00304 00305 // adding maps by spatial coords is slower than doing it directly, but 00306 // allows for precisely subtracting unaligned maps, and/or maps of 00307 // different resolutions 00308 00309 if ( USE_UNION) { 00310 // UNION VERSION 00311 init_from_union(mapA, mapB, newvol); 00312 for (long i=0; i<newvol->gridsize(); i++){ 00313 float x, y, z; 00314 float voxelA, voxelB; 00315 newvol->voxel_coord(i, x, y, z); 00316 if (interp) { 00317 voxelA = mapA->voxel_value_interpolate_from_coord(x,y,z); 00318 voxelB = mapB->voxel_value_interpolate_from_coord(x,y,z); 00319 } else { 00320 voxelA = mapA->voxel_value_from_coord(x,y,z); 00321 voxelB = mapB->voxel_value_from_coord(x,y,z); 00322 } 00323 00324 if (!myisnan(voxelA) && !myisnan(voxelB)) 00325 newvol->data[i] = voxelA * voxelB; 00326 else if (!myisnan(voxelA) && myisnan(voxelB)) 00327 newvol->data[i] = voxelA; 00328 else if (myisnan(voxelA) && !myisnan(voxelB)) 00329 newvol->data[i] = voxelB; 00330 else 00331 newvol->data[i] = 0.; 00332 } 00333 } else { 00334 // INTERSECTION VERSION 00335 init_from_intersection(mapA, mapB, newvol); 00336 00337 for (long i=0; i<newvol->gridsize(); i++){ 00338 float x, y, z; 00339 newvol->voxel_coord(i, x, y, z); 00340 00341 if (interp) newvol->data[i] = \ 00342 mapA->voxel_value_interpolate_from_coord(x,y,z) * \ 00343 mapB->voxel_value_interpolate_from_coord(x,y,z); 00344 else newvol->data[i] = \ 00345 mapA->voxel_value_from_coord(x,y,z) * \ 00346 mapB->voxel_value_from_coord(x,y,z); 00347 } 00348 00349 } 00350 } 00351 00352 void average(VolumetricData *mapA, VolumetricData *mapB, VolumetricData *newvol, bool interp, bool USE_UNION) { 00353 00354 // adding maps by spatial coords is slower than doing it directly, but 00355 // allows for precisely subtracting unaligned maps, and/or maps of 00356 // different resolutions 00357 00358 if ( USE_UNION) { 00359 // UNION VERSION 00360 init_from_union(mapA, mapB, newvol); 00361 } else { 00362 // INTERSECTION VERSION 00363 init_from_intersection(mapA, mapB, newvol); 00364 } 00365 for (long i=0; i<newvol->gridsize(); i++){ 00366 float x, y, z; 00367 newvol->voxel_coord(i, x, y, z); 00368 00369 if (interp) newvol->data[i] = \ 00370 (mapA->voxel_value_interpolate_from_coord_safe(x,y,z) + \ 00371 mapB->voxel_value_interpolate_from_coord_safe(x,y,z))*0.5f; 00372 else newvol->data[i] = \ 00373 (mapA->voxel_value_from_coord_safe(x,y,z) + \ 00374 mapB->voxel_value_from_coord_safe(x,y,z))*0.5f; 00375 } 00376 00377 } 00378 00379 void vol_moveto(VolumetricData *vol, float *com, float *pos){ 00380 float origin[3] = {0.0, 0.0, 0.0}; 00381 origin[0] = (float)vol->origin[0]; 00382 origin[1] = (float)vol->origin[1]; 00383 origin[2] = (float)vol->origin[2]; 00384 00385 float transvector[3]; 00386 vec_sub(transvector, pos, com); 00387 vec_add(origin, origin, transvector); 00388 vol->origin[0] = origin[0]; 00389 vol->origin[1] = origin[1]; 00390 vol->origin[2] = origin[2]; 00391 } 00392 /* 00393 void vectrans(float *npoint, float *mat, double *vec){ 00394 npoint[0]=vec[0]*mat[0]+vec[1]*mat[4]+vec[2]*mat[8]; 00395 npoint[1]=vec[0]*mat[1]+vec[1]*mat[5]+vec[2]*mat[9]; 00396 npoint[2]=vec[0]*mat[2]+vec[1]*mat[6]+vec[2]*mat[10]; 00397 } 00398 */ 00399 void vol_move(VolumetricData *vol, float *mat){ 00400 float origin[3] = {0.0, 0.0, 0.0}; 00401 origin[0] = (float)vol->origin[0]; 00402 origin[1] = (float)vol->origin[1]; 00403 origin[2] = (float)vol->origin[2]; 00404 00405 float transvector[3] = {mat[12], mat[13], mat[14]}; 00406 00407 float npoint[3]; 00408 00409 //deal with origin transformation 00410 //vectrans 00411 npoint[0]=origin[0]*mat[0]+origin[1]*mat[4]+origin[2]*mat[8]; 00412 npoint[1]=origin[0]*mat[1]+origin[1]*mat[5]+origin[2]*mat[9]; 00413 npoint[2]=origin[0]*mat[2]+origin[1]*mat[6]+origin[2]*mat[10]; 00414 00415 vec_add(origin, npoint, transvector); 00416 vol->origin[0] = origin[0]; 00417 vol->origin[1] = origin[1]; 00418 vol->origin[2] = origin[2]; 00419 00420 //deal with delta transformation 00421 double deltax[3] = {vol->xaxis[0],vol->xaxis[1],vol->xaxis[2]}; 00422 double deltay[3] = {vol->yaxis[0],vol->yaxis[1],vol->yaxis[2]}; 00423 double deltaz[3] = {vol->zaxis[0],vol->zaxis[1],vol->zaxis[2]}; 00424 00425 float npointx[3]; 00426 float npointy[3]; 00427 float npointz[3]; 00428 vectrans(npointx, mat, deltax); 00429 vectrans(npointy, mat, deltay); 00430 vectrans(npointz, mat, deltaz); 00431 00432 for (int i = 0; i<3; i++){ 00433 vol->xaxis[i] = npointx[i]; 00434 vol->yaxis[i] = npointy[i]; 00435 vol->zaxis[i] = npointz[i]; 00436 } 00437 00438 } 00439 00443 void histogram( VolumetricData *vol, int nbins, long *bins, float *midpts) { 00444 //get minmax values of map 00445 float min, max; 00446 vol->datarange(min, max); 00447 // Calculate the width of each bin 00448 double binwidth = (max-min)/nbins; 00449 //precompute inverse 00450 double binwidthinv = 1/binwidth; 00451 // Allocate array that will contain the number of voxels in each bin 00452 //int *bins = (int*) malloc(nbins*sizeof(int)); 00453 memset(bins, 0, nbins*sizeof(long)); 00454 00455 // Calculate histogram 00456 for (long i=0; i<vol->gridsize(); i++) 00457 bins[long((vol->data[i]-min)*binwidthinv)]++; 00458 00459 for (int j = 0; j < nbins; j++) 00460 midpts[j] = float(min + (0.5f*binwidth) + (j*binwidth)); 00461 } 00462 00464 int write_file(VMDApp *app, Molecule *volmol, int volid, const char* filename) { 00465 00466 /* Get the extension */ 00467 // XXX - casting to avoid a compilation error on SOLARIS. 00468 // char* ext = strrchr(filename, '.'); 00469 char* ext = (char *)strrchr(filename, '.'); 00470 // if (!ext) { 00471 // fprintf(stderr, "Couldn't find an extension in the filename '%s'\n", filename); 00472 // return NULL; 00473 // } 00474 00475 ext = &ext[1]; 00476 00477 vmdplugin_t *p = app->get_plugin("mol file reader", ext); 00478 00479 if (!p) { 00480 p = app->get_plugin("mol file converter", ext); 00481 } 00482 00483 MolFilePlugin *plugin = NULL; 00484 00485 if (p) { 00486 plugin = new MolFilePlugin(p); 00487 } 00488 00489 if (plugin != NULL){ 00490 int natom = 0; 00491 int filehandler = plugin->init_write(filename, natom); 00492 if (filehandler == MOLFILE_ERROR) { 00493 return 0; 00494 } 00495 00496 plugin->write_volumetric(volmol, volid); 00497 } 00498 00499 return 1; 00500 00501 }