TclSegmentation.C Source 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: 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