Advertisement
Guest User

Eddy_Em

a guest
Sep 1st, 2011
200
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 9.11 KB | None | 0 0
  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
Advertisement