API tools faq
paste
Guest User

Eddy_Em

a guest
Sep 1st, 2011
201
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 9.11 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  * Функции быстрого расчета медианы украдены из
  3.  *
  4.  *   Fast median search: an ANSI C implementation
  5.  *    Nicolas Devillard - ndevilla AT free DOT fr
  6.  *                  July 1998
  7.  */
  8. #define PIX_SORT(a,b) { if ((**a)>(**b)) ELEM_SWAP((a),(b)); }
  9. #define ELEM_SWAP(a,b) { register float **t=(a);(a)=(b);(b)=t; }
  10. float opt_med3(float ***p, int n __attribute__((unused))){
  11.     PIX_SORT(p[0],p[1]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[0],p[1]) ;
  12.     return(**p[1]) ;
  13. }float opt_med5(float ***p, int n __attribute__((unused))){
  14.     PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ;
  15.     PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ;
  16.     PIX_SORT(p[1],p[2]) ; return(**p[2]) ;
  17. }float opt_med7(float ***p, int n __attribute__((unused))){
  18.     PIX_SORT(p[0], p[5]) ; PIX_SORT(p[0], p[3]) ; PIX_SORT(p[1], p[6]) ;
  19.     PIX_SORT(p[2], p[4]) ; PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[5]) ;
  20.     PIX_SORT(p[2], p[6]) ; PIX_SORT(p[2], p[3]) ; PIX_SORT(p[3], p[6]) ;
  21.     PIX_SORT(p[4], p[5]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[1], p[3]) ;
  22.     PIX_SORT(p[3], p[4]) ; return (**p[3]) ;
  23. }float opt_med9(float ***p, int n __attribute__((unused))){
  24.     PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
  25.     PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ;
  26.     PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
  27.     PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ;
  28.     PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ;
  29.     PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ;
  30.     PIX_SORT(p[4], p[2]) ; return(**p[4]) ;
  31. } float opt_med25(float ***p, int n __attribute__((unused))){
  32.     PIX_SORT(p[0], p[1])  ; PIX_SORT(p[3], p[4])  ; PIX_SORT(p[2], p[4]) ;
  33.     PIX_SORT(p[2], p[3])  ; PIX_SORT(p[6], p[7])  ; PIX_SORT(p[5], p[7]) ;
  34.     PIX_SORT(p[5], p[6])  ; PIX_SORT(p[9], p[10]) ; PIX_SORT(p[8], p[10]) ;
  35.     PIX_SORT(p[8], p[9])  ; PIX_SORT(p[12], p[13]); PIX_SORT(p[11], p[13]) ;
  36.     PIX_SORT(p[11], p[12]); PIX_SORT(p[15], p[16]); PIX_SORT(p[14], p[16]) ;
  37.     PIX_SORT(p[14], p[15]); PIX_SORT(p[18], p[19]); PIX_SORT(p[17], p[19]) ;
  38.     PIX_SORT(p[17], p[18]); PIX_SORT(p[21], p[22]); PIX_SORT(p[20], p[22]) ;
  39.     PIX_SORT(p[20], p[21]); PIX_SORT(p[23], p[24]); PIX_SORT(p[2], p[5]) ;
  40.     PIX_SORT(p[3], p[6])  ; PIX_SORT(p[0], p[6])  ; PIX_SORT(p[0], p[3]) ;
  41.     PIX_SORT(p[4], p[7])  ; PIX_SORT(p[1], p[7])  ; PIX_SORT(p[1], p[4]) ;
  42.     PIX_SORT(p[11], p[14]); PIX_SORT(p[8], p[14]) ; PIX_SORT(p[8], p[11]) ;
  43.     PIX_SORT(p[12], p[15]); PIX_SORT(p[9], p[15]) ; PIX_SORT(p[9], p[12]) ;
  44.     PIX_SORT(p[13], p[16]); PIX_SORT(p[10], p[16]); PIX_SORT(p[10], p[13]) ;
  45.     PIX_SORT(p[20], p[23]); PIX_SORT(p[17], p[23]); PIX_SORT(p[17], p[20]) ;
  46.     PIX_SORT(p[21], p[24]); PIX_SORT(p[18], p[24]); PIX_SORT(p[18], p[21]) ;
  47.     PIX_SORT(p[19], p[22]); PIX_SORT(p[8], p[17]) ; PIX_SORT(p[9], p[18]) ;
  48.     PIX_SORT(p[0], p[18]) ; PIX_SORT(p[0], p[9])  ; PIX_SORT(p[10], p[19]) ;
  49.     PIX_SORT(p[1], p[19]) ; PIX_SORT(p[1], p[10]) ; PIX_SORT(p[11], p[20]) ;
  50.     PIX_SORT(p[2], p[20]) ; PIX_SORT(p[2], p[11]) ; PIX_SORT(p[12], p[21]) ;
  51.     PIX_SORT(p[3], p[21]) ; PIX_SORT(p[3], p[12]) ; PIX_SORT(p[13], p[22]) ;
  52.     PIX_SORT(p[4], p[22]) ; PIX_SORT(p[4], p[13]) ; PIX_SORT(p[14], p[23]) ;
  53.     PIX_SORT(p[5], p[23]) ; PIX_SORT(p[5], p[14]) ; PIX_SORT(p[15], p[24]) ;
  54.     PIX_SORT(p[6], p[24]) ; PIX_SORT(p[6], p[15]) ; PIX_SORT(p[7], p[16]) ;
  55.     PIX_SORT(p[7], p[19]) ; PIX_SORT(p[13], p[21]); PIX_SORT(p[15], p[23]) ;
  56.     PIX_SORT(p[7], p[13]) ; PIX_SORT(p[7], p[15]) ; PIX_SORT(p[1], p[9]) ;
  57.     PIX_SORT(p[3], p[11]) ; PIX_SORT(p[5], p[17]) ; PIX_SORT(p[11], p[17]) ;
  58.     PIX_SORT(p[9], p[17]) ; PIX_SORT(p[4], p[10]) ; PIX_SORT(p[6], p[12]) ;
  59.     PIX_SORT(p[7], p[14]) ; PIX_SORT(p[4], p[6])  ; PIX_SORT(p[4], p[7]) ;
  60.     PIX_SORT(p[12], p[14]); PIX_SORT(p[10], p[14]); PIX_SORT(p[6], p[7]) ;
  61.     PIX_SORT(p[10], p[12]); PIX_SORT(p[6], p[10]) ; PIX_SORT(p[6], p[17]) ;
  62.     PIX_SORT(p[12], p[17]); PIX_SORT(p[7], p[17]) ; PIX_SORT(p[7], p[10]) ;
  63.     PIX_SORT(p[12], p[18]); PIX_SORT(p[7], p[12]) ; PIX_SORT(p[10], p[18]) ;
  64.     PIX_SORT(p[12], p[20]); PIX_SORT(p[10], p[20]); PIX_SORT(p[10], p[12]) ;
  65.     return (**p[12]);
  66. }
  67. float quick_select(float ***arr, int n){
  68.     int low, high;
  69.     int median;
  70.     int middle, ll, hh;
  71.     float ret;
  72.     low = 0 ; high = n-1 ; median = (low + high) / 2;
  73.     for(;;){
  74.         if(high <= low) /* One element only */
  75.             break;
  76.         if(high == low + 1){ /* Two elements only */
  77.             PIX_SORT(arr[low], arr[high]) ;
  78.             break;
  79.         }
  80.         /* Find median of low, middle and high items; swap into position low */
  81.         middle = (low + high) / 2;
  82.         PIX_SORT(arr[middle], arr[high]) ;
  83.         PIX_SORT(arr[low], arr[high]) ;
  84.         PIX_SORT(arr[middle], arr[low]) ;
  85.         /* Swap low item (now in position middle) into position (low+1) */
  86.         ELEM_SWAP(arr[middle], arr[low+1]) ;
  87.         /* Nibble from each end towards middle, swapping items when stuck */
  88.         ll = low + 1;
  89.         hh = high;
  90.         for(;;){
  91.             do ll++; while (**arr[low] > **arr[ll]);
  92.             do hh--; while (**arr[hh] > **arr[low]);
  93.             if(hh < ll) break;
  94.             ELEM_SWAP(arr[ll], arr[hh]) ;
  95.         }
  96.         /* Swap middle item (in position low) back into correct position */
  97.         ELEM_SWAP(arr[low], arr[hh]) ;
  98.         /* Re-set active partition */
  99.         if (hh <= median) low = ll;
  100.         if (hh >= median) high = hh - 1;
  101.     }
  102.     ret = **arr[median];
  103.     return ret;
  104. }
  105. #undef COPYARR
  106. #undef PIX_SORT
  107. #undef ELEM_SWAP
  108.  
  109. /*
  110.  * Медианная фильтрация
  111.  * однопоточная:
  112.  *      127секунд для 3Кх3К с фильтром 32х32
  113.  *      33секунды для 3Кх3К с фильтром 16x16
  114.  *      4.4секунды 3Кх3К с фильтром 5х5
  115.  * четырехпоточная:
  116.  *      45секунд для 3Кх3К с фильтром 41x41
  117.  *      27секунд для 3Кх3К с фильтром 32х32
  118.  *      7.3секунд для 3Кх3К с фильтром 16x16
  119.  *      1.5секунд 3Кх3К с фильтром 6x6
  120.  *      1.1секунда 3Кх3К с фильтром 5х5 (med_opt)
  121.  * вход:
  122.  *      ima - изображение (free выполнять в вызвавшей программе)
  123.  *      f - фильтр
  124.  *      sizex, sizey - размеры изображения
  125.  * выход:
  126.  *      result - отфильтрованное изображение, память выделяется в этой процедуре
  127.  *              края изображения, не попавшие в фильтр (+- 1/2 размера фильтра),
  128.  *              заполняются нулями
  129.  */
  130. int MedFilter(float *ima,
  131.                     float **result,
  132.                     Filter *f,
  133.                     int sizex,
  134.                     int sizey
  135.                     ){
  136.     int xlow = f->w / 2, xhigh = f->w - xlow;// границы области фильтра [x0-xlow, x0+xhigh)
  137.     int ylow = f->h / 2, yhigh = f->h - ylow; // [y0-ylow, y0+yhigh)
  138.     int x, xm=sizex-xhigh, ym=sizey-yhigh; // xm,ym - верхние границы фильтруемого изображения
  139.     int ssize=sizex*sizey, fsz=f->w*f->h, H = f->h - 1;
  140.     #ifdef EBUG
  141.     double t0 = dtime();
  142.     #endif
  143.     float (*medfn)(float ***p, int n);
  144.     *result = calloc(ssize, sizeof(float));
  145.     if(!result) return FALSE;
  146.     switch(f->w*f->h){
  147.         case 3:
  148.             medfn = opt_med3;
  149.         break;
  150.         case 5:
  151.             medfn = opt_med5;
  152.         break;
  153.         case 7:
  154.             medfn = opt_med7;
  155.         break;
  156.         case 9:
  157.             medfn = opt_med9;
  158.         break;
  159.         case 25:
  160.             medfn = opt_med25;
  161.         break;
  162.         default:
  163.             medfn = quick_select;
  164.         break;
  165.     }
  166.     /*
  167.      * выбор элементов изображения по квадрату wxh в массив arr
  168.      * x0,y0 - координаты центра квадрата
  169.      * cntr - текущая заменяемая строка, либо -1, если необходимо заполнять arr целиком
  170.      * при cntr != -1 очередная строка заменяет строку cntr
  171.      */
  172.     float selarr0(int x0, int y0, float **arr, float ***sel){
  173.         int x,y,tx,ty;
  174.         float **tmp, *ptr;
  175.         tx=x0+xhigh; ty=y0+yhigh;
  176.         tmp = arr;
  177.         for(y=y0-ylow; y<ty; y++){ // полностью заполняем массив
  178.             ptr = &ima[y*sizex+x0-xlow]; // ЗАПОЛНЯЕМ ПО СТРОКАМ!!!
  179.             for(x=x0-xlow; x<tx; x++)
  180.                 *tmp++ = ptr++;
  181.         }
  182.         return medfn(sel, fsz);
  183.     }
  184.     float selarr1(int x0, int y0, float **arr, float ***sel, int *cntr){
  185.         int x,tx,ty;
  186.         float **tmp, *ptr;
  187.         tx=x0+xhigh; ty=y0+yhigh;
  188.         tmp = &arr[*cntr*f->w]; // указатель на заменяемую строку
  189.         ptr = &ima[(ty-1)*sizex+x0-xlow];
  190.         for(x=x0-xlow; x<tx; x++) // заменяем данные в столбце
  191.             *tmp++ = ptr++;
  192.         if(++(*cntr) > H) *cntr = 0;
  193.         return medfn(sel, fsz);
  194.     }
  195.     #pragma omp parallel
  196.     {
  197.         float **arr = calloc(fsz, sizeof(float*)); // массив для хранения выборки
  198.         float ***sel = calloc(fsz, sizeof(float*)); // массив для хранения отсортированной выборки
  199.         for(x = 0; x < fsz; x++) sel[x] = arr + x;
  200.         if(arr){
  201.         #pragma omp for
  202.         for(x = xlow; x < xm; x++){
  203.             int y, cntr = 0;
  204.             (*result)[ylow*sizex + x] = selarr0(x,ylow, arr, sel);
  205.             for(y = ylow+1; y < ym; y++){
  206.                 (*result)[y*sizex + x] = selarr1(x,y, arr, sel, &cntr);
  207.             }
  208.         }
  209.         free(arr);
  210.         }
  211.     }
  212.     DBG("time=%f\n", dtime()-t0);
  213.     return TRUE;
  214. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!