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: TclSegmentation.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.15 $ $Date: 2020年10月15日 17:44:27 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * Tcl bindings for density map segmentation functions 00019 * 00020 ***************************************************************************/ 00021 #include <stdio.h> 00022 #include <stdlib.h> 00023 #include <string.h> 00024 #include <float.h> // FLT_MAX etc 00025 00026 #include "Inform.h" 00027 #include "utilities.h" 00028 00029 #include "AtomSel.h" 00030 #include "VMDApp.h" 00031 #include "MoleculeList.h" 00032 #include "VolumetricData.h" 00033 #include "VolMapCreate.h" // volmap_write_dx_file() 00034 00035 #include "Segmentation.h" 00036 #include "MDFF.h" 00037 #include <math.h> 00038 #include <tcl.h> 00039 #include "TclCommands.h" 00040 00041 00042 int return_segment_usage(Tcl_Interp *interp) { 00043 Tcl_SetResult(interp, (char *) "usage: segmentation " 00044 "segment: -mol <mol> -vol <vol> [options]\n" 00045 " options:\n" 00046 " -mol <molid> specifies an already loaded density's molid\n" 00047 " for use as the segmentation volume source\n" 00048 " -vol <volume id> specifies an already loaded density's\n" 00049 " volume id for use as the volume source. Defaults to 0.\n" 00050 " -groups <count> iterate merging groups until no more than the\n" 00051 " target number of segment groups remain\n" 00052 " -watershed_blur <sigma> initial Gaussian blur sigma to use\n" 00053 " prior to the first pass of the Watershed algorithm\n" 00054 " -starting_sigma <sigma> initial Gaussian blur sigma to use\n" 00055 " for the iterative scale-space segmentation algorithm\n" 00056 " -blur_multiple <multfactor> multiplicative constant to apply\n" 00057 " to the scale-space blur sigma at each iteration\n" 00058 " -separate_groups flag causes each segment group to be emitted\n" 00059 " as an additional density map, with all other group voxels\n" 00060 " masked out, maintaining the dimensions of the original map\n" 00061 , 00062 TCL_STATIC); 00063 00064 return 0; 00065 } 00066 00067 00068 int segment_volume(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) { 00069 int verbose = 0; 00070 if (argc < 3) { 00071 return_segment_usage(interp); 00072 return TCL_ERROR; 00073 } 00074 00075 int seg_def_groups = 128; 00076 float def_blur_sigma = 1.0f; 00077 float def_initial_sigma = 1.0f; 00078 float def_blur_multiple = 1.5f; 00079 int def_separate_groups = 0; 00080 00081 int num_groups_target = seg_def_groups; 00082 double watershed_blur_sigma = def_blur_sigma; 00083 double blur_starting_sigma = def_initial_sigma; 00084 double blur_multiple = def_blur_multiple; 00085 int separate_groups = def_separate_groups; 00086 00087 int molid=-1; 00088 int volid=0; 00089 for (int i=0; i < argc; i++) { 00090 char *opt = Tcl_GetStringFromObj(objv[i], NULL); 00091 00092 if (!strcmp(opt, "-mol")) { 00093 if (i == argc-1) { 00094 Tcl_AppendResult(interp, "No molid specified",NULL); 00095 return TCL_ERROR; 00096 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) { 00097 Tcl_AppendResult(interp, "\n molid incorrectly specified",NULL); 00098 return TCL_ERROR; 00099 } 00100 } 00101 00102 if (!strcmp(opt, "-vol")) { 00103 if (i == argc-1){ 00104 Tcl_AppendResult(interp, "No volume id specified",NULL); 00105 return TCL_ERROR; 00106 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &volid) != TCL_OK) { 00107 Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL); 00108 return TCL_ERROR; 00109 } 00110 } 00111 00112 if (!strcmp(opt, "-groups")) { 00113 if (i == argc-1){ 00114 Tcl_AppendResult(interp, "No target group count specified",NULL); 00115 return TCL_ERROR; 00116 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &num_groups_target) != TCL_OK) { 00117 Tcl_AppendResult(interp, "\n target group count incorrectly specified",NULL); 00118 return TCL_ERROR; 00119 } 00120 } 00121 00122 if (!strcmp(opt, "-watershed_blur")) { 00123 if (i == argc-1){ 00124 Tcl_AppendResult(interp, "No target blur sigma specified",NULL); 00125 return TCL_ERROR; 00126 } else if ( Tcl_GetDoubleFromObj(interp, objv[i+1], &watershed_blur_sigma) != TCL_OK) { 00127 Tcl_AppendResult(interp, "\n watershed blur sigma incorrectly specified",NULL); 00128 return TCL_ERROR; 00129 } 00130 } 00131 00132 if (!strcmp(opt, "-starting_sigma")) { 00133 if (i == argc-1){ 00134 Tcl_AppendResult(interp, "No starging sigma specified",NULL); 00135 return TCL_ERROR; 00136 } else if ( Tcl_GetDoubleFromObj(interp, objv[i+1], &watershed_blur_sigma) != TCL_OK) { 00137 Tcl_AppendResult(interp, "\n starting sigma incorrectly specified",NULL); 00138 return TCL_ERROR; 00139 } 00140 } 00141 00142 if (!strcmp(opt, "-blur_multiple")) { 00143 if (i == argc-1){ 00144 Tcl_AppendResult(interp, "No blur multiple specified",NULL); 00145 return TCL_ERROR; 00146 } else if ( Tcl_GetDoubleFromObj(interp, objv[i+1], &blur_multiple) != TCL_OK) { 00147 Tcl_AppendResult(interp, "\n blur multiple incorrectly specified",NULL); 00148 return TCL_ERROR; 00149 } 00150 } 00151 00152 if (!strcmp(opt, "-separate_groups")) { 00153 separate_groups = 1; 00154 } 00155 00156 if (!strcmp(opt, "-verbose") || (getenv("VMDSEGMENTATIONVERBOSE") != NULL)) { 00157 verbose = 1; 00158 } 00159 } 00160 00161 if (verbose) 00162 msgInfo << "Verbose segmentation diagnostic output enabled." << sendmsg; 00163 00164 MoleculeList *mlist = app->moleculeList; 00165 const VolumetricData *map = NULL; 00166 00167 if (molid > -1) { 00168 Molecule *volmol = mlist->mol_from_id(molid); 00169 if (volmol != NULL) 00170 map = volmol->get_volume_data(volid); 00171 } else { 00172 Tcl_AppendResult(interp, "\n no target volume specified",NULL); 00173 return TCL_ERROR; 00174 } 00175 00176 if (map == NULL) { 00177 Tcl_AppendResult(interp, "\n no target volume correctly specified",NULL); 00178 return TCL_ERROR; 00179 } 00180 00181 // perform the segmentation 00182 msgInfo << "Performing scale-space image segmentation on mol[" 00183 << molid << "], volume[" << volid << "]..." 00184 << sendmsg; 00185 00186 PROFILE_PUSH_RANGE("Segmentation (All Steps)", 6); 00187 00188 // initialize a segmentation object, perform memory allocations, etc. 00189 Segmentation seg(map, 1); 00190 00191 // compute the segmentation 00192 seg.segment(num_groups_target, float(watershed_blur_sigma), 00193 float(blur_starting_sigma), float(blur_multiple), 00194 MERGE_HILL_CLIMB); 00195 00196 // get the group results back 00197 00198 // allocate output array for volumetric map of segment group indices, and 00199 // perform any required type conversions and changes in storage locality 00200 float *fp32seggrpids = new float[map->gridsize()]; 00201 seg.get_results<float>(fp32seggrpids); 00202 00203 PROFILE_POP_RANGE(); 00204 00205 // propagate results to persistent storage 00206 app->molecule_add_volumetric(molid, "Segmentation", map->origin, 00207 map->xaxis, map->yaxis, map->zaxis, 00208 map->xsize, map->ysize, map->zsize, 00209 fp32seggrpids); 00210 00211 // Use segmentation group information to mask and split the 00212 // segmented density map into multiple multiple density maps 00213 if (separate_groups) { 00214 PROFILE_PUSH_RANGE("Segmentation split groups", 1); 00215 msgInfo << "Splitting density map into individual group density maps" << sendmsg; 00216 // allocate output array for volumetric map of segment group indices, and 00217 // perform any required type conversions and changes in storage locality 00218 long *l64seggrpids = new long[map->gridsize()]; 00219 seg.get_results<long>(l64seggrpids); 00220 00221 char title[4096]; 00222 int nGroups = seg.get_num_groups(); 00223 long nVoxels = map->gridsize(); 00224 for (int i=0; i<nGroups; i++) { 00225 // each new density map will be owned/managed by VMD after addition 00226 // to the target molecule 00227 float *output = new float[nVoxels]; 00228 for (long v=0; v<nVoxels; v++) { 00229 output[v] = l64seggrpids[v] == i ? map->data[v] : 0.0f; 00230 } 00231 printf("Masking density map by individual groups: %.2f%% \r", 100 * i / (double)nGroups); 00232 snprintf(title, sizeof(title), "Segment Group %d", i); 00233 00234 // propagate results to persistent storage 00235 app->molecule_add_volumetric(molid, "Segmentation", map->origin, 00236 map->xaxis, map->yaxis, map->zaxis, 00237 map->xsize, map->ysize, map->zsize, 00238 output); 00239 } 00240 00241 msgInfo << "Finished processing segmentation groups." << sendmsg; 00242 PROFILE_POP_RANGE(); 00243 } 00244 00245 msgInfo << "Segmentation tasks complete." << sendmsg; 00246 return TCL_OK; 00247 } 00248 00249 00250 int obj_segmentation(ClientData cd, Tcl_Interp *interp, int argc, Tcl_Obj * const objv[]){ 00251 if (argc < 2) { 00252 return_segment_usage(interp); 00253 return TCL_ERROR; 00254 } 00255 char *argv1 = Tcl_GetStringFromObj(objv[1],NULL); 00256 00257 VMDApp *app = (VMDApp *)cd; 00258 if (!strupncmp(argv1, "segment", CMDLEN)) 00259 return segment_volume(app, argc-1, objv+1, interp); 00260 00261 return_segment_usage(interp); 00262 return TCL_OK; 00263 } 00264 00265