00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
00056 image_type = DT_FLOAT;
00057 gaussian_f = new GaussianBlur<float>(map->data, width, height, depth, cuda);
00058
00059
00060
00061
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
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();
00205 return 0;
00206 }
00207
00208
00209
00210
00211 update_type();
00212
00213
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();
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();
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();
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
00295
00296
00297
00298
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
00337
00338 case DT_ERROR:
00339 default:
00340 printf("Segmentation::switch_to_cpu(): Invalid voxel type was set\n");
00341 break;
00342 }
00343
00344
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
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
00489
00490
00491
00492
00493
00494
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;
00514 GROUP_T max_val = GROUP_T(0)-1;
00515
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
00530
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