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

Segmentation.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: Segmentation.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.29 $        $Date: 2020/10/15 16:07:31 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Scale-space variant of Watershed image segmentation intended for use  
00019  *   on 3-D cryo-EM density maps.
00020  ***************************************************************************/
00021 
00022 #include <limits.h>
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <cmath>
00026 #include <string.h>
00027 #include "Segmentation.h"
00028 #include "Watershed.h"
00029 #include "CUDASegmentation.h"
00030 #include "CUDAGaussianBlur.h"
00031 #include "Inform.h"
00032 #include "ProfileHooks.h"
00033 
00034 #ifdef TIMER
00035 #include "WKFUtils.h"
00036 #endif
00037 
00038 #define INST_GET_RESULTS(G_T) \
00039   template void Segmentation::get_results<G_T>(G_T* results);
00040 
00041 INST_GET_RESULTS(unsigned short)
00042 INST_GET_RESULTS(short)
00043 INST_GET_RESULTS(unsigned int)
00044 INST_GET_RESULTS(int)
00045 INST_GET_RESULTS(unsigned long)
00046 INST_GET_RESULTS(long)
00047 INST_GET_RESULTS(float)
00048 
00049 
00050 Segmentation::Segmentation(const VolumetricData *map, bool cuda) {
00051   PROFILE_PUSH_RANGE("Segmentation::Segmentation(MAP)", 6);
00052 
00053   init(map->xsize, map->ysize, map->zsize, cuda);
00054 
00055   //XXX in future, check VolData to see which type we should use here
00056   image_type = DT_FLOAT;
00057   gaussian_f = new GaussianBlur<float>(map->data, width, height, depth, cuda);
00058 
00059   /* Example
00060    * image_type = DT_UCHAR;
00061    * gaussian_uc = new GaussianBlur<unsigned char>(map->data_c, width, height, depth, cuda);
00062    * */
00063 
00064   PROFILE_POP_RANGE();
00065 }
00066 
00067 
00068 void Segmentation::init(unsigned long w,
00069                         unsigned long h,
00070                         unsigned long d,
00071                         bool cuda) {
00072   width     = w;
00073   height    = h;
00074   depth     = d;
00075   nVoxels   = h * w * d;
00076   nGroups   = nVoxels;
00077   use_cuda  = cuda;
00078   group_type = get_new_group_type();
00079   image_type = DT_ERROR;
00080   segments_us = NULL;
00081   segments_ui = NULL;
00082   segments_ul = NULL;
00083   gaussian_f = NULL;
00084   gaussian_us = NULL;
00085   gaussian_uc = NULL;
00086 
00087   switch (group_type) {
00088     case DT_ULONG:
00089       segments_ul  = allocate_array<unsigned long>(nVoxels);
00090       break;
00091     case DT_UINT:
00092       segments_ui  = allocate_array<unsigned int>(nVoxels);
00093       break;
00094     case DT_USHORT:
00095       segments_us  = allocate_array<unsigned short>(nVoxels);
00096       break;
00097     default:
00098       fprintf(stderr, "Error: invalid segment type.\n");
00099       return;
00100   }
00101 }
00102 
00103 
00104 
00105 Segmentation::Segmentation(float* data,
00106                            molfile_volumetric_t* metadata,
00107                            bool cuda) {
00108   PROFILE_PUSH_RANGE("Segmentation::Segmentation(float *data)", 6);
00109 
00110   init(metadata->xsize, metadata->ysize, metadata->zsize, cuda);
00111   image_type = DT_FLOAT;
00112 
00113   gaussian_f = new GaussianBlur<float>(data, width, height, depth, use_cuda);
00114 
00115   PROFILE_POP_RANGE();
00116 }
00117 
00118 
00119 Segmentation::Segmentation(unsigned short* data,
00120                            molfile_volumetric_t* metadata,
00121                            bool cuda) {
00122   PROFILE_PUSH_RANGE("Segmentation::Segmentation(float *data)", 6);
00123 
00124   init(metadata->xsize, metadata->ysize, metadata->zsize, cuda);
00125   image_type = DT_USHORT;
00126 
00127   gaussian_us = new GaussianBlur<unsigned short>(data, width, height, depth, use_cuda);
00128 
00129   PROFILE_POP_RANGE();
00130 }
00131 
00132 Segmentation::Segmentation(unsigned char* data,
00133                            molfile_volumetric_t* metadata,
00134                            bool cuda) {
00135   PROFILE_PUSH_RANGE("Segmentation::Segmentation(float *data)", 6);
00136 
00137   init(metadata->xsize, metadata->ysize, metadata->zsize, cuda);
00138   image_type = DT_UCHAR;
00139 
00140   gaussian_uc = new GaussianBlur<unsigned char>(data, width, height, depth, use_cuda);
00141 
00142   PROFILE_POP_RANGE();
00143 }
00144 
00145 
00146 Segmentation::~Segmentation() {
00147   if (gaussian_f != NULL) {
00148     delete gaussian_f;
00149   }
00150   if (gaussian_us != NULL) {
00151     delete gaussian_us;
00152   }
00153   if (gaussian_uc != NULL) {
00154     delete gaussian_uc;
00155   }
00156   if (segments_ul != NULL) {
00157     free_array(segments_ul);
00158   }
00159   if (segments_ui != NULL) {
00160     free_array(segments_ui);
00161   }
00162   if (segments_us != NULL) {
00163     free_array(segments_us);
00164   }
00165 }
00166 
00167 
00168 
00169 double Segmentation::segment(int   num_final_groups,
00170                              float watershed_blur_sigma,
00171                              float blur_initial_sigma,
00172                              float blur_multiple,
00173                              MERGE_POLICY policy,
00174                              const bool isverbose) {
00175   PROFILE_PUSH_RANGE("Segmentation::segment()", 6);
00176 
00177   double time_sec = 0.0;
00178 #ifdef TIMER
00179   wkf_timerhandle timer = wkf_timer_create();    
00180   wkf_timer_start(timer);
00181 #endif
00182   
00183   verbose = isverbose;
00184 
00185   // Perform watershed segmentation
00186   switch (group_type) {
00187     case DT_ULONG:
00188       if (verbose)
00189         printf("Using unsigned long representation.\n");
00190       get_watershed_segmentation<unsigned long>(watershed_blur_sigma, segments_ul);
00191       break;
00192     case DT_UINT:
00193       if (verbose)
00194         printf("Using unsigned int representation.\n");
00195       get_watershed_segmentation<unsigned int>(watershed_blur_sigma, segments_ui);
00196       break;
00197     case DT_USHORT:
00198       if (verbose)
00199         printf("Using unsigned short representation.\n");
00200       get_watershed_segmentation<unsigned short>(watershed_blur_sigma, segments_us);
00201       break;
00202     default:
00203       fprintf(stderr, "Error: invalid segmentation group type.\n");
00204       PROFILE_POP_RANGE(); // error return point
00205       return 0;
00206   }
00207 
00208 
00209   // this would switch to smaller datatype for scale-space filtering
00210   // currently not implemented
00211   update_type();
00212 
00213   // merge groups until we have <= num_final_groups
00214   switch (group_type) {
00215     case DT_ULONG:
00216       merge_groups<unsigned long>(num_final_groups, blur_initial_sigma, blur_multiple, policy, segments_ul);
00217       break;
00218     case DT_UINT:
00219       merge_groups<unsigned int>(num_final_groups, blur_initial_sigma, blur_multiple, policy, segments_ui);
00220       break;
00221     case DT_USHORT:
00222       merge_groups<unsigned short>(num_final_groups, blur_initial_sigma, blur_multiple, policy, segments_us);
00223       break;
00224     default:
00225       fprintf(stderr, "Error: invalid segment type.\n");
00226       PROFILE_POP_RANGE(); // error return point
00227       return 0;
00228   }
00229 
00230 #ifdef TIMER
00231   wkf_timer_stop(timer);
00232   time_sec = wkf_timer_time(timer);
00233   wkf_timer_destroy(timer);
00234   if (verbose) {
00235     char msgbuf[2048];
00236     sprintf(msgbuf, "Total segmentation time: %.3f seconds.", time_sec);
00237     msgInfo << msgbuf << sendmsg;
00238   }
00239 #endif
00240 
00241   PROFILE_POP_RANGE(); // final return point
00242 
00243   return time_sec;
00244 }
00245 
00246 
00247 
00248 template <typename GROUP_T>
00249 void Segmentation::merge_groups(unsigned long num_final_groups,
00250                                 float blur_initial_sigma,
00251                                 float blur_multiple,
00252                                 MERGE_POLICY policy,
00253                                 GROUP_T* segments) {
00254   PROFILE_PUSH_RANGE("Segmentation::merge_groups()", 0);
00255 
00256   switch (image_type) {
00257     case DT_FLOAT:
00258       merge_groups_helper<GROUP_T, float>(num_final_groups, blur_initial_sigma,
00259                                           blur_multiple, segments, policy, gaussian_f);
00260       break;
00261     case DT_USHORT:
00262       merge_groups_helper<GROUP_T, unsigned short>(num_final_groups, blur_initial_sigma,
00263                                                    blur_multiple, segments, policy, gaussian_us);
00264       break;
00265     case DT_UCHAR:
00266       merge_groups_helper<GROUP_T, unsigned char>(num_final_groups, blur_initial_sigma,
00267                                                   blur_multiple, segments, policy, gaussian_uc);
00268       break;
00269     default:
00270       fprintf(stderr, "Error: invalid image type.\n");
00271       PROFILE_POP_RANGE(); // error return point
00272       break;
00273   }
00274 
00275   PROFILE_POP_RANGE();
00276 }
00277 
00278 template <typename GROUP_T, typename IMAGE_T>
00279 void Segmentation::merge_groups_helper(unsigned long num_final_groups,
00280                                        float blur_initial_sigma,
00281                                        float blur_multiple,
00282                                        GROUP_T* segments,
00283                                        MERGE_POLICY policy,
00284                                        GaussianBlur<IMAGE_T>* g) {
00285   PROFILE_PUSH_RANGE("Segmentation::merge_groups_helper()", 0);
00286 
00287   if (verbose)
00288     printf("Num groups = %ld\n", nGroups);
00289 
00290   ScaleSpaceFilter<GROUP_T, IMAGE_T> filter(width, height, depth, 
00291                                    nGroups, blur_initial_sigma,
00292                                    blur_multiple, use_cuda);
00293 
00294   // XXX right now we have one policy for the entire merge.
00295   // It might give a better/faster segmentation to start out with MERGE_HILL_CLIMB
00296   // and transition to a MERGE_WATERSHED* policy for integer
00297   // image types once we stop being able to merge groups.
00298   // This might be more work than it is worth though.
00299   while (nGroups > num_final_groups) {
00300     nGroups = filter.merge(segments, g, policy);
00301     if (verbose)
00302       printf("Num groups = %ld\n", nGroups);
00303   }
00304 
00305   PROFILE_POP_RANGE();
00306 }
00307 
00308 
00309 void Segmentation::switch_to_cpu() {
00310   switch(group_type) {
00311     case DT_ULONG:
00312       {
00313       unsigned long* new_seg_array = new unsigned long[nVoxels];
00314       copy_and_convert_type<unsigned long>(segments_ul, new_seg_array, nVoxels, true);
00315       free_array(segments_ul);
00316       segments_ul = new_seg_array;
00317       break;
00318       }
00319     case DT_UINT:
00320       {
00321       unsigned int* new_seg_array = new unsigned int[nVoxels];
00322       copy_and_convert_type<unsigned int>(segments_ui, new_seg_array, nVoxels, true);
00323       free_array(segments_ui);
00324       segments_ui = new_seg_array;
00325       break;
00326       }
00327     case DT_USHORT:
00328       {
00329       unsigned short* new_seg_array = new unsigned short[nVoxels];
00330       copy_and_convert_type<unsigned short>(segments_us, new_seg_array, nVoxels, true);
00331       free_array(segments_us);
00332       segments_us = new_seg_array;
00333       break;
00334       }
00335 
00336     // this should really never happen but xlC complains if we don't
00337     // include all enums in the switch block
00338     case DT_ERROR:
00339     default:
00340       printf("Segmentation::switch_to_cpu(): Invalid voxel type was set\n");
00341       break;
00342   }
00343 
00344   //TODO switch gaussian_fblur to cpu
00345   use_cuda = false;
00346 }
00347 
00348 
00349 
00350 
00351 template <typename GROUP_T>
00352 void Segmentation::get_watershed_segmentation(float watershed_blur_sigma,
00353                                               GROUP_T* segments) {
00354   switch (image_type) {
00355     case DT_FLOAT:
00356       {
00357         gaussian_f->blur(watershed_blur_sigma);
00358         int imageongpu=1;
00359         float *imageptr = gaussian_f->get_image_d();
00360         if (!imageptr) {
00361           imageptr = gaussian_f->get_image();
00362           imageongpu=0;
00363         }
00364         Watershed<GROUP_T, float> watershed(height, width, depth, use_cuda);
00365         watershed.watershed(imageptr, imageongpu, segments, verbose);
00366         break;
00367       }
00368     case DT_USHORT:
00369       {
00370         gaussian_us->blur(watershed_blur_sigma);
00371         int imageongpu=1;
00372         unsigned short *imageptr = gaussian_us->get_image_d();
00373         if (!imageptr) {
00374           imageptr = gaussian_us->get_image();
00375           imageongpu=0;
00376         }
00377         Watershed<GROUP_T, unsigned short> watershed(height, width, depth, use_cuda);
00378         watershed.watershed(imageptr, imageongpu, segments, verbose);
00379         break;
00380       }
00381     case DT_UCHAR:
00382       {
00383         gaussian_uc->blur(watershed_blur_sigma);
00384         int imageongpu=1;
00385         unsigned char *imageptr = gaussian_uc->get_image_d();
00386         if (!imageptr) {
00387           imageptr = gaussian_uc->get_image();
00388           imageongpu=0;
00389         }
00390         Watershed<GROUP_T, unsigned char> watershed(height, width, depth, use_cuda);
00391         watershed.watershed(imageptr, imageongpu, segments, verbose);
00392         break;
00393       }
00394     default:
00395       fprintf(stderr, "Invalid image type.\n");
00396       break;
00397   }
00398 
00399   GROUP_T* group_map = allocate_array<GROUP_T>(nGroups);
00400   sequentialize_groups<GROUP_T>(segments, group_map);
00401   free_array(group_map);
00402 }
00403 
00404 
00405 
00406 unsigned long Segmentation::get_num_groups() {
00407   return nGroups;
00408 }
00409 
00410 
00411 
00412 template <typename OUT_T>
00413 void Segmentation::get_results(OUT_T* results) {
00414   PROFILE_PUSH_RANGE("Segmentation::get_results()", 1);
00415 
00416   switch (group_type) {
00417     case DT_ULONG:
00418       copy_and_convert_type<unsigned long, OUT_T>(segments_ul, results, nVoxels, use_cuda);
00419       break;
00420     case DT_UINT:
00421       copy_and_convert_type<unsigned int, OUT_T>(segments_ui, results, nVoxels, use_cuda);
00422       break;
00423     case DT_USHORT:
00424       copy_and_convert_type<unsigned short, OUT_T>(segments_us, results, nVoxels, use_cuda);
00425       break;
00426     default:
00427       fprintf(stderr, "Error: invalid segment type.\n");
00428       break;
00429   }
00430 
00431   PROFILE_POP_RANGE();
00432 }
00433 
00434 
00435 
00436 template <typename IN_T, typename OUT_T>
00437 void Segmentation::copy_and_convert_type(IN_T* in, OUT_T* out, unsigned long num_elems, bool copy_from_gpu) {
00438 #if defined(VMDCUDA)
00439   if (copy_from_gpu) {
00440     IN_T* local_arr = new IN_T[num_elems];
00441     copy_array_from_gpu(local_arr, in, num_elems * sizeof(IN_T));
00442 
00443 #ifdef WATERSHED_COMPATABILITY_TEST
00444     // dirty hack because gpu and cpu sequentialization don't give same ordering
00445     IN_T* group_map = new IN_T[nGroups];
00446     sequentialize_groups_cpu<IN_T>(local_arr, group_map);
00447     delete [] group_map;
00448 #endif
00449 
00450     for (unsigned long i = 0; i < num_elems; i++) {
00451       out[i] = (OUT_T)local_arr[i];
00452     }
00453     delete [] local_arr;
00454   } else 
00455 #endif
00456   {
00457 #if defined(VMDCUDA)
00458     if (use_cuda) {
00459       copy_and_convert_type_cuda<IN_T, OUT_T>(in, out, num_elems);
00460     } else 
00461 #endif
00462     {
00463       for (unsigned long i = 0; i < num_elems; i++) {
00464         out[i] = (OUT_T)in[i];
00465       }
00466     }
00467   }
00468 }
00469 
00470 
00471 
00472 Segmentation::Datatype Segmentation::get_new_group_type() {
00473   Datatype type;
00474   if (nGroups < USHRT_MAX) {
00475     type = DT_USHORT;
00476   } else if (nGroups < UINT_MAX) {
00477     type = DT_UINT;
00478   } else if (nGroups < ULONG_MAX) {
00479     type = DT_ULONG;
00480   } else {
00481     fprintf(stderr, "Error: images this large are not currently supported.\n");
00482     fprintf(stderr, "Warning: Undefined behavior if you continue to use this object.\n");
00483     type = DT_ERROR;
00484   }
00485   return type;
00486 }
00487 
00488 // these next two functions are only needed once after
00489 // we get the intitial segmentation from watershed.
00490 // They are duplicated in the ScaleSpaceFilter class
00491 // so would be nice to find a way to remove them.
00492 // We need them before we create a ScaleSpaceFilter obj.
00493 // so we know know the max group # and can use decalare
00494 // the obj with appropriate datatype
00495 
00496 
00497 template <typename GROUP_T>
00498 void Segmentation::sequentialize_groups(GROUP_T* segments, GROUP_T* group_map) {
00499 #if defined(VMDCUDA)
00500   if (use_cuda) {
00501     nGroups = sequentialize_groups_cuda<GROUP_T>(segments, group_map, nVoxels, nGroups);
00502   } else 
00503 #endif
00504   {
00505     nGroups = sequentialize_groups_cpu<GROUP_T>(segments, group_map);
00506   }
00507 }
00508 
00509 
00510 
00511 template <typename GROUP_T>
00512 unsigned long Segmentation::sequentialize_groups_cpu(GROUP_T* segments, GROUP_T* group_map) {
00513   unsigned long numGroups = 0; //TODO changed from 0, update CUDA code
00514   GROUP_T max_val = GROUP_T(0)-1;
00515   // GROUP_T is only unsigned integer types so this ^ is defined behavior of C++ standard
00516   memset(group_map, max_val, nGroups * sizeof(GROUP_T));
00517 
00518   for (unsigned long i = 0; i < nVoxels; ++i) {
00519     const GROUP_T curr_group = segments[i];
00520     if (group_map[curr_group] == max_val) {
00521       group_map[curr_group] = GROUP_T(numGroups++);
00522     }
00523     segments[i] = group_map[curr_group];
00524   }
00525   return numGroups;
00526 }
00527 
00528 
00529 // TODO check if we can use a smaller type,
00530 // and if so transition segments to the new type
00531 
00532 void Segmentation::update_type() {
00533   Datatype new_type = get_new_group_type();
00534   if (new_type != group_type) {
00535     switch (group_type) {
00536       case DT_ULONG:
00537         if (new_type == DT_UINT) {
00538           if (verbose) {
00539             printf("Converting unsigned ulong to uint.\n");
00540           }
00541           segments_ui = allocate_array<unsigned int>(nVoxels);
00542           copy_and_convert_type<unsigned long, unsigned int>(segments_ul, segments_ui, nVoxels);
00543         } else {
00544           if (verbose) {
00545             printf("Converting unsigned ulong to ushort.\n");
00546           }
00547           segments_us = allocate_array<unsigned short>(nVoxels);
00548           copy_and_convert_type<unsigned long, unsigned short>(segments_ul, segments_us, nVoxels);
00549         }
00550         free_array(segments_ul);
00551         break;
00552       case DT_UINT:
00553         if (verbose) {
00554           printf("Converting uint to ushort.\n");
00555         }
00556         segments_us = allocate_array<unsigned short>(nVoxels);
00557         copy_and_convert_type<unsigned int, unsigned short>(segments_ui, segments_us, nVoxels);
00558         free_array(segments_ui);
00559         break;
00560       default:
00561         printf("Error in type conversion.\n");
00562         break;
00563     }
00564   }
00565   group_type = new_type;
00566 }
00567 
00568 
00569 
00570 template <typename GROUP_T>
00571 GROUP_T* Segmentation::allocate_array(unsigned long num_elements) {
00572   GROUP_T* arr;
00573 #if defined(VMDCUDA)
00574   if (use_cuda) {
00575     arr = (GROUP_T*) alloc_cuda_array(num_elements * sizeof(GROUP_T));
00576   } else 
00577 #endif
00578   {
00579     arr = new GROUP_T[num_elements];
00580   }
00581   return arr;
00582 }
00583 
00584 
00585 
00586 template <typename GROUP_T>
00587 void Segmentation::free_array(GROUP_T*& arr) {
00588 #if defined(VMDCUDA)
00589   if (use_cuda) {
00590     free_cuda_array(arr);
00591   } else 
00592 #endif
00593   {
00594     delete [] arr;
00595   }
00596   arr = NULL;
00597 }
00598 

Generated on Sat Nov 9 02:45:30 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002