Eddy_Em

/*
 * Медианная фильтрация
 * однопоточная:
 *      127секунд для 3Кх3К с фильтром 32х32
 *      33секунды для 3Кх3К с фильтром 16x16
 *      4.4секунды 3Кх3К с фильтром 5х5
 * четырехпоточная:
 *      1974секунды для 3Кх3К с фильтром 256x256
 *      40секунд для 3Кх3К с фильтром 32х32
 *      10секунд для 3Кх3К с фильтром 16x16
 *      1.4секунды 3Кх3К с фильтром 5х5
 * вход:
 *      ima - изображение (free выполнять в вызвавшей программе)
 *      f - фильтр
 *      sizex, sizey - размеры изображения
 * выход:
 *      result - отфильтрованное изображение, память выделяется в этой процедуре
 *              края изображения, не попавшие в фильтр (+- 1/2 размера фильтра),
 *              заполняются нулями
 */
int MedFilter(float *ima __attribute__((unused)),
                    float **result  __attribute__((unused)),
                    Filter *f __attribute__((unused)),
                    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
    *result = calloc(ssize, sizeof(float));
    if(!result) return FALSE;
    /*
     * выбор элементов изображения по квадрату wxh в массив arr
     * x0,y0 - координаты центра квадрата
     * cntr - текущая заменяемая строка, либо -1, если необходимо заполнять arr целиком
     * при cntr != -1 очередная строка заменяет строку cntr
     */
    float selarr(int x0, int y0, float *arr, int *cntr){
        int x,y,tx,ty;
        float *tmp, *ptr;
        tx=x0+xhigh; ty=y0+yhigh;
        if(*cntr==-1){ // начало очередного столбца
            tmp = arr;
            for(y=y0-ylow; y<ty; y++){ // полностью заполняем массив
                ptr = &ima[y*sizex+x0-xlow]; // ЗАПОЛНЯЕМ ПО СТРОКАМ!!!
                for(x=x0-xlow; x<tx; x++)
                    *tmp++ = *ptr++;
            }
            *cntr=0; // в следующий раз заменяем первый столбец
        }else{
            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 quick_select(arr, fsz);
    }
    #pragma omp parallel
    {
        float *arr = calloc(fsz, sizeof(float)); // массив для хранения выборки
        if(arr){
        #pragma omp for
        for(x = xlow; x < xm; x++){
            int y, cntr = -1;
            for(y = ylow; y < ym; y++){
                (*result)[y*sizex + x] = selarr(x,y, arr, &cntr);
            }
        }
        free(arr);
        }
    }
    DBG("time=%f\n", dtime()-t0);
    return TRUE;
}