Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /*
- * Функции быстрого расчета медианы украдены из
- *
- * Fast median search: an ANSI C implementation
- * Nicolas Devillard - ndevilla AT free DOT fr
- * July 1998
- */
- #define PIX_SORT(a,b) { if ((**a)>(**b)) ELEM_SWAP((a),(b)); }
- #define ELEM_SWAP(a,b) { register float **t=(a);(a)=(b);(b)=t; }
- float opt_med3(float ***p, int n __attribute__((unused))){
- PIX_SORT(p[0],p[1]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[0],p[1]) ;
- return(**p[1]) ;
- }float opt_med5(float ***p, int n __attribute__((unused))){
- PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ;
- PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ;
- PIX_SORT(p[1],p[2]) ; return(**p[2]) ;
- }float opt_med7(float ***p, int n __attribute__((unused))){
- PIX_SORT(p[0], p[5]) ; PIX_SORT(p[0], p[3]) ; PIX_SORT(p[1], p[6]) ;
- PIX_SORT(p[2], p[4]) ; PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[5]) ;
- PIX_SORT(p[2], p[6]) ; PIX_SORT(p[2], p[3]) ; PIX_SORT(p[3], p[6]) ;
- PIX_SORT(p[4], p[5]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[1], p[3]) ;
- PIX_SORT(p[3], p[4]) ; return (**p[3]) ;
- }float opt_med9(float ***p, int n __attribute__((unused))){
- PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
- PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ;
- PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
- PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ;
- PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ;
- PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ;
- PIX_SORT(p[4], p[2]) ; return(**p[4]) ;
- } float opt_med25(float ***p, int n __attribute__((unused))){
- PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[2], p[4]) ;
- PIX_SORT(p[2], p[3]) ; PIX_SORT(p[6], p[7]) ; PIX_SORT(p[5], p[7]) ;
- PIX_SORT(p[5], p[6]) ; PIX_SORT(p[9], p[10]) ; PIX_SORT(p[8], p[10]) ;
- PIX_SORT(p[8], p[9]) ; PIX_SORT(p[12], p[13]); PIX_SORT(p[11], p[13]) ;
- PIX_SORT(p[11], p[12]); PIX_SORT(p[15], p[16]); PIX_SORT(p[14], p[16]) ;
- PIX_SORT(p[14], p[15]); PIX_SORT(p[18], p[19]); PIX_SORT(p[17], p[19]) ;
- PIX_SORT(p[17], p[18]); PIX_SORT(p[21], p[22]); PIX_SORT(p[20], p[22]) ;
- PIX_SORT(p[20], p[21]); PIX_SORT(p[23], p[24]); PIX_SORT(p[2], p[5]) ;
- PIX_SORT(p[3], p[6]) ; PIX_SORT(p[0], p[6]) ; PIX_SORT(p[0], p[3]) ;
- PIX_SORT(p[4], p[7]) ; PIX_SORT(p[1], p[7]) ; PIX_SORT(p[1], p[4]) ;
- PIX_SORT(p[11], p[14]); PIX_SORT(p[8], p[14]) ; PIX_SORT(p[8], p[11]) ;
- PIX_SORT(p[12], p[15]); PIX_SORT(p[9], p[15]) ; PIX_SORT(p[9], p[12]) ;
- PIX_SORT(p[13], p[16]); PIX_SORT(p[10], p[16]); PIX_SORT(p[10], p[13]) ;
- PIX_SORT(p[20], p[23]); PIX_SORT(p[17], p[23]); PIX_SORT(p[17], p[20]) ;
- PIX_SORT(p[21], p[24]); PIX_SORT(p[18], p[24]); PIX_SORT(p[18], p[21]) ;
- PIX_SORT(p[19], p[22]); PIX_SORT(p[8], p[17]) ; PIX_SORT(p[9], p[18]) ;
- PIX_SORT(p[0], p[18]) ; PIX_SORT(p[0], p[9]) ; PIX_SORT(p[10], p[19]) ;
- PIX_SORT(p[1], p[19]) ; PIX_SORT(p[1], p[10]) ; PIX_SORT(p[11], p[20]) ;
- PIX_SORT(p[2], p[20]) ; PIX_SORT(p[2], p[11]) ; PIX_SORT(p[12], p[21]) ;
- PIX_SORT(p[3], p[21]) ; PIX_SORT(p[3], p[12]) ; PIX_SORT(p[13], p[22]) ;
- PIX_SORT(p[4], p[22]) ; PIX_SORT(p[4], p[13]) ; PIX_SORT(p[14], p[23]) ;
- PIX_SORT(p[5], p[23]) ; PIX_SORT(p[5], p[14]) ; PIX_SORT(p[15], p[24]) ;
- PIX_SORT(p[6], p[24]) ; PIX_SORT(p[6], p[15]) ; PIX_SORT(p[7], p[16]) ;
- PIX_SORT(p[7], p[19]) ; PIX_SORT(p[13], p[21]); PIX_SORT(p[15], p[23]) ;
- PIX_SORT(p[7], p[13]) ; PIX_SORT(p[7], p[15]) ; PIX_SORT(p[1], p[9]) ;
- PIX_SORT(p[3], p[11]) ; PIX_SORT(p[5], p[17]) ; PIX_SORT(p[11], p[17]) ;
- PIX_SORT(p[9], p[17]) ; PIX_SORT(p[4], p[10]) ; PIX_SORT(p[6], p[12]) ;
- PIX_SORT(p[7], p[14]) ; PIX_SORT(p[4], p[6]) ; PIX_SORT(p[4], p[7]) ;
- PIX_SORT(p[12], p[14]); PIX_SORT(p[10], p[14]); PIX_SORT(p[6], p[7]) ;
- PIX_SORT(p[10], p[12]); PIX_SORT(p[6], p[10]) ; PIX_SORT(p[6], p[17]) ;
- PIX_SORT(p[12], p[17]); PIX_SORT(p[7], p[17]) ; PIX_SORT(p[7], p[10]) ;
- PIX_SORT(p[12], p[18]); PIX_SORT(p[7], p[12]) ; PIX_SORT(p[10], p[18]) ;
- PIX_SORT(p[12], p[20]); PIX_SORT(p[10], p[20]); PIX_SORT(p[10], p[12]) ;
- return (**p[12]);
- }
- float quick_select(float ***arr, int n){
- int low, high;
- int median;
- int middle, ll, hh;
- float ret;
- low = 0 ; high = n-1 ; median = (low + high) / 2;
- for(;;){
- if(high <= low) /* One element only */
- break;
- if(high == low + 1){ /* Two elements only */
- PIX_SORT(arr[low], arr[high]) ;
- break;
- }
- /* Find median of low, middle and high items; swap into position low */
- middle = (low + high) / 2;
- PIX_SORT(arr[middle], arr[high]) ;
- PIX_SORT(arr[low], arr[high]) ;
- PIX_SORT(arr[middle], arr[low]) ;
- /* Swap low item (now in position middle) into position (low+1) */
- ELEM_SWAP(arr[middle], arr[low+1]) ;
- /* Nibble from each end towards middle, swapping items when stuck */
- ll = low + 1;
- hh = high;
- for(;;){
- do ll++; while (**arr[low] > **arr[ll]);
- do hh--; while (**arr[hh] > **arr[low]);
- if(hh < ll) break;
- ELEM_SWAP(arr[ll], arr[hh]) ;
- }
- /* Swap middle item (in position low) back into correct position */
- ELEM_SWAP(arr[low], arr[hh]) ;
- /* Re-set active partition */
- if (hh <= median) low = ll;
- if (hh >= median) high = hh - 1;
- }
- ret = **arr[median];
- return ret;
- }
- #undef COPYARR
- #undef PIX_SORT
- #undef ELEM_SWAP
- /*
- * Медианная фильтрация
- * однопоточная:
- * 127секунд для 3Кх3К с фильтром 32х32
- * 33секунды для 3Кх3К с фильтром 16x16
- * 4.4секунды 3Кх3К с фильтром 5х5
- * четырехпоточная:
- * 45секунд для 3Кх3К с фильтром 41x41
- * 27секунд для 3Кх3К с фильтром 32х32
- * 7.3секунд для 3Кх3К с фильтром 16x16
- * 1.5секунд 3Кх3К с фильтром 6x6
- * 1.1секунда 3Кх3К с фильтром 5х5 (med_opt)
- * вход:
- * ima - изображение (free выполнять в вызвавшей программе)
- * f - фильтр
- * sizex, sizey - размеры изображения
- * выход:
- * result - отфильтрованное изображение, память выделяется в этой процедуре
- * края изображения, не попавшие в фильтр (+- 1/2 размера фильтра),
- * заполняются нулями
- */
- int MedFilter(float *ima,
- float **result,
- Filter *f,
- int sizex,
- int sizey
- ){
- int xlow = f->w / 2, xhigh = f->w - xlow;// границы области фильтра [x0-xlow, x0+xhigh)
- int ylow = f->h / 2, yhigh = f->h - ylow; // [y0-ylow, y0+yhigh)
- int x, xm=sizex-xhigh, ym=sizey-yhigh; // xm,ym - верхние границы фильтруемого изображения
- int ssize=sizex*sizey, fsz=f->w*f->h, H = f->h - 1;
- #ifdef EBUG
- double t0 = dtime();
- #endif
- float (*medfn)(float ***p, int n);
- *result = calloc(ssize, sizeof(float));
- if(!result) return FALSE;
- switch(f->w*f->h){
- case 3:
- medfn = opt_med3;
- break;
- case 5:
- medfn = opt_med5;
- break;
- case 7:
- medfn = opt_med7;
- break;
- case 9:
- medfn = opt_med9;
- break;
- case 25:
- medfn = opt_med25;
- break;
- default:
- medfn = quick_select;
- break;
- }
- /*
- * выбор элементов изображения по квадрату wxh в массив arr
- * x0,y0 - координаты центра квадрата
- * cntr - текущая заменяемая строка, либо -1, если необходимо заполнять arr целиком
- * при cntr != -1 очередная строка заменяет строку cntr
- */
- float selarr0(int x0, int y0, float **arr, float ***sel){
- int x,y,tx,ty;
- float **tmp, *ptr;
- tx=x0+xhigh; ty=y0+yhigh;
- tmp = arr;
- for(y=y0-ylow; y<ty; y++){ // полностью заполняем массив
- ptr = &ima[y*sizex+x0-xlow]; // ЗАПОЛНЯЕМ ПО СТРОКАМ!!!
- for(x=x0-xlow; x<tx; x++)
- *tmp++ = ptr++;
- }
- return medfn(sel, fsz);
- }
- float selarr1(int x0, int y0, float **arr, float ***sel, int *cntr){
- int x,tx,ty;
- float **tmp, *ptr;
- tx=x0+xhigh; ty=y0+yhigh;
- tmp = &arr[*cntr*f->w]; // указатель на заменяемую строку
- ptr = &ima[(ty-1)*sizex+x0-xlow];
- for(x=x0-xlow; x<tx; x++) // заменяем данные в столбце
- *tmp++ = ptr++;
- if(++(*cntr) > H) *cntr = 0;
- return medfn(sel, fsz);
- }
- #pragma omp parallel
- {
- float **arr = calloc(fsz, sizeof(float*)); // массив для хранения выборки
- float ***sel = calloc(fsz, sizeof(float*)); // массив для хранения отсортированной выборки
- for(x = 0; x < fsz; x++) sel[x] = arr + x;
- if(arr){
- #pragma omp for
- for(x = xlow; x < xm; x++){
- int y, cntr = 0;
- (*result)[ylow*sizex + x] = selarr0(x,ylow, arr, sel);
- for(y = ylow+1; y < ym; y++){
- (*result)[y*sizex + x] = selarr1(x,y, arr, sel, &cntr);
- }
- }
- free(arr);
- }
- }
- DBG("time=%f\n", dtime()-t0);
- return TRUE;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement