언어/C

C언어 영상처리 최종 정리 - filter, padding, upsampling

이게될까 2024. 4. 18. 19:19
728x90
728x90
#define _CRT_SECURE_NO_WARNINGS
#include<stdlib.h>
#include<stdio.h>
#include<windows.h>
#include<math.h>

void histo(int height,int width,double* y,int size, int stride, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);
void CumulativeHisto(int height, int width, double* y, int size, int stride, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);

void GammaOut(double*, double, BITMAPFILEHEADER, BITMAPINFOHEADER);

void applyPaddingWithAverage(double* paddedImg, int height, int width, int padding, double* origImg);
void getPadding(double* py, int height, int width, int padding, int pheight, int pwidth, double* y2);

void calFilter(double* y3, double* py, int fsize, double* filter, int height, int width, int pwidth);
double median(double* arr, int size);
void calMedianFilter(double* y3, double* py, int fsize, int height, int width, int pwidth);

void getQulity(double* y, double* y3, int width, int height, int bitCnt);

void downsampling(int raito, double* org, double* dws, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);
void averageDownsampling(int ratio, double* org, double* dws, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);

void upsampling(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);
void bilinearUpsampling(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo);

void makeOutFile(BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo, char* fileName, double* y);

int main() {
    BITMAPFILEHEADER bmpFile, bmpFile2;
    BITMAPINFOHEADER bmpInfo, bmpInfo2;
    FILE* inputFile = NULL, * inputFile2 = NULL;
    inputFile = fopen("AICenterY.bmp", "rb");

    inputFile2 = fopen("AICenterYBIUPS.bmp", "rb");

    fread(&bmpFile, sizeof(BITMAPFILEHEADER), 1, inputFile);
    fread(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, inputFile);
    fread(&bmpFile2, sizeof(BITMAPFILEHEADER), 1, inputFile2);
    fread(&bmpInfo2, sizeof(BITMAPINFOHEADER), 1, inputFile2);

    int width = bmpInfo.biWidth, height = bmpInfo.biHeight, size = bmpInfo.biSizeImage, bitCnt = bmpInfo.biBitCount, stride = (((bitCnt / 8) * width) + 3) / 4 * 4,max;
    int width2 = bmpInfo2.biWidth, height2 = bmpInfo2.biHeight, size2 = bmpInfo2.biSizeImage, bitCnt2 = bmpInfo2.biBitCount, stride2 = (((bitCnt2 / 8) * width2) + 3) / 4 * 4;

    double *y,*y2,mse=0,psnr=0;
    y = (double*)calloc(width * height, sizeof(double));
    y2 = (double*)calloc(width * height, sizeof(double));

    printf("width,height,size,bitCnt,stride\n");
    printf("%d %d %d %d %d\n",width,height,size,bitCnt,stride);
    printf("%d %d %d %d %d\n", width2, height2, size2, bitCnt2, stride2);

    unsigned char* inputImg = NULL, * inputImg2 = NULL;
    inputImg = (unsigned char*)calloc(size, sizeof(unsigned char));
    inputImg2 = (unsigned char*)calloc(size, sizeof(unsigned char));

    fread(inputImg, sizeof(unsigned char), size, inputFile);
    fread(inputImg2, sizeof(unsigned char), size, inputFile2);
    // input 
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            y[j * width + i] = inputImg[j * stride + 3 * i ];
            //y[j * width + i] = 0.299 * inputImg[j * stride + 3 * i + 2] + 0.587 * inputImg[j * stride + 3 * i + 1] + 0.114 * inputImg[j * stride + 3 * i + 0];
            y2[j * width + i] = inputImg2[j * stride + 3 * i ];
        }
    }


    // Algorithm

    // sampling
    int raito = 2;
    double* dws;
    dws = (double*)calloc(width * height/ raito/ raito, sizeof(double));
    //downsampling(raito, y, dws, bmpFile, bmpInfo);
    averageDownsampling(raito, y, dws, bmpFile, bmpInfo);
    double* ups;
    ups = (double*)calloc(width * height * raito * raito, sizeof(double));
    upsampling(raito, y, ups, bmpFile, bmpInfo);
    //bilinearUpsampling(raito, y, ups, bmpFile, bmpInfo);

    //padding
    int padding = 1, pwidth = width + padding * 2, pheight = height + padding * 2, psize = pwidth * pheight;
    double* py;
    py = (double*)calloc(pwidth * pheight, sizeof(double));

    getPadding(py,height,width, padding, pheight, pwidth, y2);
    //applyPadding(py, height, width, padding, y2);
    //applyPaddingWithAverage(py, height, width, padding, y2);


    // filter

    int fsize = 3;

    // mean filter
    double* filter;
    filter = (double*)calloc(fsize * fsize, sizeof(double));
    for (int i = 0; i < fsize * fsize; i++) filter[i] = 1.0 / (fsize * fsize);

    // Gussian Filter
    double G3filter[] = {1.0/16,2.0/16,1.0/16,2.0/16,4.0/16,2.0/16,1.0/16,2.0/16,1.0/16};
    double G5filter[] = {1/273.0, 4/ 273.0, 7/ 273.0, 4/ 273.0, 1/ 273.0, 4/ 273.0, 16/ 273.0, 26/ 273.0, 16/ 273.0, 4/ 273.0, 7/ 273.0, 26/ 273.0, 41/ 273.0, 26/ 273.0, 7/ 273.0, 4 / 273.0, 16 / 273.0, 26 / 273.0, 16 / 273.0, 4 / 273.0,1 / 273.0, 4 / 273.0, 7 / 273.0, 4 / 273.0, 1 / 273.0 };
    double* y3;
    y3 = (double*)calloc(width * height, sizeof(double));
    calFilter(y3,py, fsize,G3filter,height,width,pwidth);
    //calMedianFilter( y3,py,  fsize,  height,  width,  pwidth);

    //get Quality
    //getQulity(y,y3,width,height, bitCnt);

    //histo

    //histo(height, width, y2, size2, stride2, bmpFile2, bmpInfo2);
    //CumulativeHisto(height, width, y, size, stride, bmpFile, bmpInfo);

    // output

    makeOutFile(bmpFile,  bmpInfo, "Output.bmp",  y);
    GammaOut( y, 2.5, bmpFile,bmpInfo);

    free(inputImg);
    free(inputImg2);

    free(y);
    free(y2);
    free(py);
    free(y3);
    fclose(inputFile);
    fclose(inputFile2);
    return 0; 
}



void getPadding(double* py, int height, int width, int padding, int pheight, int pwidth, double* y2){
// 일단 원본 값 집어넣기
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            py[(j + padding) * pwidth + (i + padding)] = y2[j * width + i];
        }
    }
    // 패딩 값 집어넣기
    //하
    for (int p = 0; p < padding; p++) {
        for (int i = 0; i < width; i++) {
            py[(p)*pwidth + (i + padding)] = y2[0 * width + i];
        }
    }
    //좌
    for (int p = 0; p < padding; p++) {
        for (int j = 0; j < height; j++) {
            py[(j + padding) * pwidth + p] = y2[j * width];
        }
    }
    //우
    for (int p = 0; p < padding; p++) {
        for (int j = 0; j < height; j++) {
            py[(j + padding) * pwidth + (width + padding) + p] = y2[(j + 1) * width - 1];
        }
    }

    //상
    for (int p = 0; p < padding; p++) {
        for (int i = 0; i < width; i++) {
            py[(p + width + padding) * pwidth + (i + padding)] = y2[(height - 1) * width + i];
        }
    }

    // 좌하 ( 0,0) ~ (padding-1,padding-1)
    for (int j = 0; j < padding; j++) {
        for (int i = 0; i < padding; i++) {
            py[(j)*pwidth + (i)] = y2[0 * width + 0];
        }
    }
    // 우하 (0,pwidth-1 -padding) ~ (padding-1,pwidth-1)
    for (int j = 0; j < padding; j++) {
        for (int i = 0; i < padding; i++) {
            py[(j)*pwidth + (width + padding + i)] = y2[1 * width - 1];
        }
    }
    // 좌상
    for (int j = 0; j < padding; j++) {
        for (int i = 0; i < padding; i++) {
            py[(height + padding + j) * pwidth + (i)] = y2[(height - 1) * width];
        }
    }
    //우상
    for (int j = 0; j < padding; j++) {
        for (int i = 0; i < padding; i++) {
            py[(height + padding + j) * pwidth + (width + padding + i)] = y2[(height)*width - 1];
        }
    }
}

void calFilter(double* y3, double* py, int fsize, double* filter, int height, int width, int pwidth) {
    double sumf = 0;
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            sumf = 0;
            for (int k = 0; k < fsize; k++) {
                for (int l = 0; l < fsize; l++) {
                    sumf += filter[k * fsize + l] * py[(j + k) * pwidth + (i + l)];
                }
            }
            y3[j * width + i] = sumf;
        }
    }
}

void histo(int height, int width, double* y, int size,int stride, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    double* histogram;
    histogram = (double*)calloc(256, sizeof(double));

    // get Histogram
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))]++;
        }
    }
    // out Histogram
    double histoMax = 0;
    for (int i = 0; i < 256; i++) {
        if (histogram[i] > histoMax) histoMax = histogram[i];
    }
    for (int i = 0; i < 256; i++)histogram[i] = histogram[i] / histoMax * 450;
    unsigned char* outHisto;
    outHisto = (unsigned char*)calloc(size, sizeof(unsigned char));

    for (int j = 0; j < 512; j++) {
        for (int i = 0; i < 512; i++) {
            if (histogram[i / 2] > j) continue;
            outHisto[j * stride + 3 * i + 0] = 255;
            outHisto[j * stride + 3 * i + 1] = 255;
            outHisto[j * stride + 3 * i + 2] = 255;
        }
    }
    FILE* outputHisto = fopen("OutHisto.bmp", "wb");

    bmpInfo.biWidth = 512;
    bmpInfo.biHeight = 512;
    bmpInfo.biSizeImage = 512*3*512;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 512 * 3 * 512;

    fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputHisto);
    fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputHisto);
    fwrite(outHisto, sizeof(unsigned char), size, outputHisto);
    free(histogram);
    fclose(outputHisto);
    free(outHisto);
}


void CumulativeHisto(int height, int width, double* y, int size, int stride, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    double* histogram;
    histogram = (double*)calloc(256, sizeof(double));
    int width2 = 512, height2 = 512, stride2 = (((bmpInfo.biBitCount / 8) * width2) + 3) / 4 * 4, size2 = stride2 * height2;
    printf("%d %d %d %d\n", width2, height2, size2, stride2);

    // get Histogram
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))]++;
        }
    }
    for (int i = 0; i < 255; i++) {
        histogram[i + 1] += histogram[i];
    }

    // out Histogram
    double histoMax = 0;
    for (int i = 0; i < 256; i++) {
        if (histogram[i] > histoMax) histoMax = histogram[i];
    }
    unsigned char* outHisto, * HE;
    outHisto = (unsigned char*)calloc(size2, sizeof(unsigned char));
    HE = (unsigned char*)calloc(size, sizeof(unsigned char));


    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            HE[j * stride + 3 * i + 0] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
            HE[j * stride + 3 * i + 1] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
            HE[j * stride + 3 * i + 2] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
        }
    }
    for (int i = 0; i < 256; i++)histogram[i] = histogram[i] / histoMax * 450;

    for (int j = 0; j < height2; j++) {
        for (int i = 0; i < width2; i++) {
            if (histogram[i / 2] > j) continue;
            outHisto[j * stride2 + 3 * i + 0] = 255;
            outHisto[j * stride2 + 3 * i + 1] = 255;
            outHisto[j * stride2 + 3 * i + 2] = 255;
        }
    }


    FILE* outputHisto = fopen("OutHisto.bmp", "wb");
    bmpInfo.biWidth = 512;
    bmpInfo.biHeight = 512;
    bmpInfo.biSizeImage = 512 * 3 * 512;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 512 * 3 * 512;


    fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputHisto);
    fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputHisto);
    fwrite(outHisto, sizeof(unsigned char), size2, outputHisto);

    FILE* outputHE = fopen("OutHE.bmp", "wb");

    bmpInfo.biWidth = width;
    bmpInfo.biHeight = height;
    bmpInfo.biSizeImage = size;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size;

    fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputHE);
    fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputHE);
    fwrite(HE, sizeof(unsigned char), size, outputHE);
    fclose(outputHE);
    free(HE);



    free(histogram);
    fclose(outputHisto);
    free(outHisto);
}

void getQulity(double* y, double* y3,int width,int height,int bitCnt) {
    double mse=0,max,psnr;
    unsigned char org, out;
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            org = (unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]));
            out = (unsigned char)(y3[j * width + i] > 255 ? 255 : (y3[j * width + i] < 0 ? 0 : y3[j * width + i]));
            mse += (double)((org - out) * (org - out));

        }
    }
    mse /= (height * width);
    max = pow(2, bitCnt / 3) - 1;
    psnr = mse != 0.0 ? 10.0 * log10(max * max / mse) : 99.99;
    printf("mse = %.2lf\npsnr = %.2lf\n", mse, psnr);
}

void applyPaddingWithAverage(double* paddedImg, int height, int width, int padding, double* origImg) {
    int pheight = height + 2 * padding;
    int pwidth = width + 2 * padding;

    // 원본 이미지 데이터를 중앙에 배치
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            paddedImg[(j + padding) * pwidth + (i + padding)] = origImg[j * width + i];
        }
    }

    // 상하좌우 패딩에 평균 값 적용
    for (int j = 0; j < pheight; j++) {
        for (int i = 0; i < pwidth; i++) {
            // 좌측 패딩 처리
            if (i < padding) {
                double sum = 0;
                int count = 0;
                for (int k = -2; k <= 2; k++) {
                    int idx = j + k; // 현재 행에서 상하로 2개씩
                    if (idx >= padding && idx < height + padding) { // 범위 체크
                        sum += paddedImg[idx * pwidth + padding];
                        count++;
                    }
                }
                paddedImg[j * pwidth + i] = count > 0 ? sum / count : 0; // 평균 값 적용
            }

            // 우측 패딩 처리
            if (i >= width + padding) {
                double sum = 0;
                int count = 0;
                for (int k = -2; k <= 2; k++) {
                    int idx = j + k; // 현재 행에서 상하로 2개씩
                    if (idx >= padding && idx < height + padding) { // 범위 체크
                        sum += paddedImg[idx * pwidth + (width + padding - 1)];
                        count++;
                    }
                }
                paddedImg[j * pwidth + i] = count > 0 ? sum / count : 0; // 평균 값 적용
            }
            // 상단 패딩 처리
            for (int p = 0; p < padding; p++) {
                double sum = 0;
                int count = 0;
                for (int k = -2; k <= 2; k++) {
                    int idx = i + k; // 현재 컬럼에서 좌우로 2개씩
                    if (idx >= padding && idx < width + padding) { // 범위 체크
                        sum += paddedImg[padding * pwidth + idx];
                        count++;
                    }
                }
                paddedImg[p * pwidth + i] = count > 0 ? sum / count : 0; // 평균 값 적용
            }

            // 하단 패딩 처리
            for (int p = 0; p < padding; p++) {
                double sum = 0;
                int count = 0;
                for (int k = -2; k <= 2; k++) {
                    int idx = i + k; // 현재 컬럼에서 좌우로 2개씩
                    if (idx >= padding && idx < width + padding) { // 범위 체크
                        sum += paddedImg[(height + padding - 1) * pwidth + idx];
                        count++;
                    }
                }
                paddedImg[(height + padding + p) * pwidth + i] = count > 0 ? sum / count : 0; // 평균 값 적용
            }
        }
    }

}



void downsampling(int raito, double* org, double* dws, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int width2 = bmpInfo.biWidth/raito, height2 = bmpInfo.biHeight / raito, stride2 = (((bmpInfo.biBitCount/8)*width2)+3)/4*4,size2 = stride2*height2;
    int width = bmpInfo.biWidth, height = bmpInfo.biHeight;

    for (int j = 0; j < height2; j++) {
        for (int i = 0; i < width2; i++) {
            dws[j * width2 + i] = org[(j*raito)*width + i*raito];
        }
    }
    bmpInfo.biWidth = width2;
    bmpInfo.biHeight = height2;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutDWS.bmp", dws);

}

void averageDownsampling(int ratio, double* org, double* dws, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int width2 = bmpInfo.biWidth / ratio, height2 = bmpInfo.biHeight / ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * width2) + 3) / 4 * 4;
    int size2 = stride2 * height2;
    int width = bmpInfo.biWidth, height = bmpInfo.biHeight;

    for (int j = 0; j < height2; j++) {
        for (int i = 0; i < width2; i++) {
            double sum = 0;
            int count = 0;
            // Compute the average of pixels in the block
            for (int y = 0; y < ratio; y++) {
                for (int x = 0; x < ratio; x++) {
                    int ny = j * ratio + y;
                    int nx = i * ratio + x;
                    if (ny < height && nx < width) {
                        sum += org[ny * width + nx];
                        count++;
                    }
                }
            }
            dws[j * width2 + i] = sum / count;  // Store the average in the downsampled image
        }
    }
    bmpInfo.biWidth = width2;
    bmpInfo.biHeight = height2;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutDWS.bmp", dws);
}


void upsampling(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int width = bmpInfo.biWidth, height = bmpInfo.biHeight;
    int width2 = bmpInfo.biWidth * ratio, height2 = bmpInfo.biHeight * ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * width2) + 3) / 4 * 4, size2 = stride2 * height2;
    int nearestJ = 0, nearestI=0;
    for (int j = 0; j < height2; j++) {
        for (int i = 0; i < width2; i++) {
            nearestJ = j / ratio;
            nearestI = i / ratio;
            ups[j * width2 + i] = org[nearestJ * width + nearestI];
        }
    }


    bmpInfo.biWidth = width2;
    bmpInfo.biHeight = height2;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);


}

void bilinearUpsampling(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int originalWidth = bmpInfo.biWidth, originalHeight = bmpInfo.biHeight;
    int upscaledWidth = bmpInfo.biWidth * ratio, upscaledHeight = bmpInfo.biHeight * ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * upscaledWidth) + 3) / 4 * 4, size2 = stride2 * upscaledHeight;

    for (int y = 0; y < upscaledHeight; y++) {
        for (int x = 0; x < upscaledWidth; x++) {
            float gx = ((float)x / (float)upscaledWidth) * (float)(originalWidth - 1);
            float gy = ((float)y / (float)upscaledHeight) * (float)(originalHeight - 1);
            int gxi = (int)gx;
            int gyi = (int)gy;

            double top = (gx - gxi) * org[gyi * originalWidth + (gxi + 1)] + (gxi + 1 - gx) * org[gyi * originalWidth + gxi];
            double bottom = (gx - gxi) * org[(gyi + 1) * originalWidth + (gxi + 1)] + (gxi + 1 - gx) * org[(gyi + 1) * originalWidth + gxi];

            ups[y * upscaledWidth + x] = (gy - gyi) * bottom + (gyi + 1 - gy) * top;
        }
    }

    bmpInfo.biWidth = upscaledWidth;
    bmpInfo.biHeight = upscaledHeight;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);

}

void bilinearUpsample(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int originalWidth = bmpInfo.biWidth, originalHeight = bmpInfo.biHeight;
    int upscaledWidth = originalWidth * ratio, upscaledHeight = originalHeight * ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * upscaledWidth) + 3) / 4 * 4, size2 = stride2 * upscaledHeight;

    for (int j = 0; j < upscaledHeight; j++) {
        for (int i = 0; i < upscaledWidth; i++) {
            float x = (float)i / ratio;
            float y = (float)j / ratio;

            int x1 = (int)x;
            int y1 = (int)y;
            int x2 = x1 + 1 < originalWidth ? x1 + 1 : x1;
            int y2 = y1 + 1 < originalHeight ? y1 + 1 : y1;

            double Q11 = org[y1 * originalWidth + x1];
            double Q12 = org[y2 * originalWidth + x1];
            double Q21 = org[y1 * originalWidth + x2];
            double Q22 = org[y2 * originalWidth + x2];

            double R1 = (x2 - x) / (x2 - x1) * Q11 + (x - x1) / (x2 - x1) * Q21;
            double R2 = (x2 - x) / (x2 - x1) * Q12 + (x - x1) / (x2 - x1) * Q22;

            ups[j * upscaledWidth + i] = (y2 - y) / (y2 - y1) * R1 + (y - y1) / (y2 - y1) * R2;
        }
    }

    bmpInfo.biWidth = upscaledWidth;
    bmpInfo.biHeight = upscaledHeight;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);

}

void nTapInterpolation(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int originalWidth = bmpInfo.biWidth, originalHeight = bmpInfo.biHeight;
    int upscaledWidth = originalWidth * ratio, upscaledHeight = originalHeight * ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * upscaledWidth) + 3) / 4 * 4, size2 = stride2 * upscaledHeight;

    // Example of 3-TAP filter coefficients
    double coefficients[3] = { 0.25, 0.5, 0.25 };  // Simple averaging coefficients

    for (int j = 0; j < upscaledHeight; j++) {
        for (int i = 0; i < upscaledWidth; i++) {
            float x = (float)i / ratio;
            float y = (float)j / ratio;

            int x1 = fmax(0, (int)x - 1);
            int x2 = (int)x;
            int x3 = fmin(originalWidth - 1, (int)x + 1);

            int y1 = fmax(0, (int)y - 1);
            int y2 = (int)y;
            int y3 = fmin(originalHeight - 1, (int)y + 1);

            double sum = 0;
            sum += org[y1 * originalWidth + x1] * coefficients[0];
            sum += org[y2 * originalWidth + x2] * coefficients[1];
            sum += org[y3 * originalWidth + x3] * coefficients[2];

            ups[j * upscaledWidth + i] = sum;
        }
    }

    bmpInfo.biWidth = upscaledWidth;
    bmpInfo.biHeight = upscaledHeight;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);

}



void fourTapInterpolation(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int originalWidth = bmpInfo.biWidth, originalHeight = bmpInfo.biHeight;
    int upscaledWidth = originalWidth * ratio, upscaledHeight = originalHeight * ratio;
    int stride2 = (((bmpInfo.biBitCount / 8) * upscaledWidth) + 3) / 4 * 4, size2 = stride2 * upscaledHeight;

    // Coefficients for 4-Tap interpolation (simple averaging example)
    double coefficients[4] = { 0.2, 0.3, 0.3, 0.2 }; // These are example coefficients

    for (int j = 0; j < upscaledHeight; j++) {
        for (int i = 0; i < upscaledWidth; i++) {
            float x = (float)i / ratio;
            float y = (float)j / ratio;

            int x1 = fmax(0, (int)x - 1);
            int x2 = (int)x;
            int x3 = fmin(originalWidth - 1, (int)x + 1);
            int x4 = fmin(originalWidth - 1, (int)x + 2);

            int y1 = fmax(0, (int)y - 1);
            int y2 = (int)y;
            int y3 = fmin(originalHeight - 1, (int)y + 1);
            int y4 = fmin(originalHeight - 1, (int)y + 2);

            double sum = 0;
            sum += org[y1 * originalWidth + x1] * coefficients[0];
            sum += org[y2 * originalWidth + x2] * coefficients[1];
            sum += org[y3 * originalWidth + x3] * coefficients[2];
            sum += org[y4 * originalWidth + x4] * coefficients[3];

            ups[j * upscaledWidth + i] = sum;
        }
    }

    bmpInfo.biWidth = upscaledWidth;
    bmpInfo.biHeight = upscaledHeight;
    bmpInfo.biSizeImage = size2;
    bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
    makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);

}


// Median filter 함수
double median(double* arr, int size) {
    // 배열을 정렬
    for (int i = 0; i < size - 1; i++) {
        for (int j = i + 1; j < size; j++) {
            if (arr[i] > arr[j]) {
                double tmp = arr[i];
                arr[i] = arr[j];
                arr[j] = tmp;
            }
        }
    }
    // 중간값 반환
    if (size % 2 == 0) {
        return (arr[size / 2 - 1] + arr[size / 2]) / 2.0;
    }
    else {
        return arr[size / 2];
    }
}

// 미디언 필터 적용 함수
void calMedianFilter(double* y3, double* py, int fsize, int height, int width, int pwidth) {
    int filterSize = fsize * fsize;
    double* window = (double*)malloc(filterSize * sizeof(double));

    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            int count = 0;
            // 각 픽셀 주변의 값을 window 배열에 저장
            for (int k = 0; k < fsize; k++) {
                for (int l = 0; l < fsize; l++) {
                    int xIdx = i + l - fsize / 2;
                    int yIdx = j + k - fsize / 2;
                    // 이미지 경계 처리
                    if (xIdx >= 0 && xIdx < width && yIdx >= 0 && yIdx < height) {
                        window[count++] = py[yIdx * pwidth + xIdx];
                    }
                }
            }
            // 미디언 값을 결과 배열에 저장
            y3[j * width + i] = median(window, count);
        }
    }

    free(window);
}

void makeOutFile(BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo, char* fileName, double* y) {
    int width = bmpInfo.biWidth, height = bmpInfo.biHeight, size = bmpInfo.biSizeImage, bitCnt = bmpInfo.biBitCount, stride = (((bitCnt / 8) * width) + 3) / 4 * 4;
    unsigned char* outputImg = NULL;
    outputImg = (unsigned char*)calloc(size, sizeof(unsigned char));
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            for (int k = 0; k < 3; k++) outputImg[j * stride + 3 * i + k] = (unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]));

            //outputImg[j * stride + 3 * i + 1] = (unsigned char)(y3[j * width + i] > 255 ? 255 : (y3[j * width + i] < 0 ? 0 : y3[j * width + i]));
            //outputImg[j * stride + 3 * i + 2] = (unsigned char)(y3[j * width + i] > 255 ? 255 : (y3[j * width + i] < 0 ? 0 : y3[j * width + i]));
        }
    }
    FILE* outputFile = fopen(fileName, "wb");
    fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputFile);
    fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputFile);
    fwrite(outputImg, sizeof(unsigned char), size, outputFile);
    fclose(outputFile);
    free(outputImg);
}

void GammaOut(double* y, double g, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
    int width = bmpInfo.biWidth, height = bmpInfo.biHeight, size = bmpInfo.biSizeImage;
    double* gy;
    gy = (double*)calloc(width * height, sizeof(double));

    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            gy[j * width + i] = 255 * pow(y[j * width + i] / 255.0, 1/g);
        }
    }
    makeOutFile(bmpFile, bmpInfo, "gammaOut.bmp", gy);
    free(gy);

}

노이즈 제거를 위해 사용할 수 있는 필터는 다양하며, 각각 특정 유형의 노이즈에 대해 더 적합하거나 특별한 기능을 수행할 수 있습니다. 다음은 평균 필터와 가우시안 필터 외에 추천할 수 있는 10가지 필터입니다:

  1. 미디언 필터 (Median Filter):
    • 솔트 앤 페퍼 노이즈 같은 극단적인 값에 강합니다. 각 픽셀 주변의 중간 값을 선택합니다.

// Median filter 함수
double median(double* arr, int size) {
    // 배열을 정렬
    for (int i = 0; i < size - 1; i++) {
        for (int j = i + 1; j < size; j++) {
            if (arr[i] > arr[j]) {
                double tmp = arr[i];
                arr[i] = arr[j];
                arr[j] = tmp;
            }
        }
    }
    // 중간값 반환
    if (size % 2 == 0) {
        return (arr[size / 2 - 1] + arr[size / 2]) / 2.0;
    }
    else {
        return arr[size / 2];
    }
}

// 미디언 필터 적용 함수
void calMedianFilter(double* y3, double* py, int fsize, int height, int width, int pwidth) {
    int filterSize = fsize * fsize;
    double* window = (double*)malloc(filterSize * sizeof(double));

    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            int count = 0;
            // 각 픽셀 주변의 값을 window 배열에 저장
            for (int k = 0; k < fsize; k++) {
                for (int l = 0; l < fsize; l++) {
                    int xIdx = i + l - fsize / 2;
                    int yIdx = j + k - fsize / 2;
                    // 이미지 경계 처리
                    if (xIdx >= 0 && xIdx < width && yIdx >= 0 && yIdx < height) {
                        window[count++] = py[yIdx * pwidth + xIdx];
                    }
                }
            }
            // 미디언 값을 결과 배열에 저장
            y3[j * width + i] = median(window, count);
        }
    }

    free(window);
}
  1. 양방향 필터 (Bilateral Filter):
    • 에지를 보존하면서 노이즈를 제거합니다. 공간적 근접성과 픽셀 간의 강도 차이 모두를 고려합니다.
double gaussian(double x, double sigma) {
    return exp(-(x * x) / (2 * sigma * sigma));
}

void bilateralFilter(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingSize) {
    int filterRadius = 3;
    double sigmaSpace = 10.0;
    double sigmaColor = 25.0;

    // 이중 반복문을 통해 이미지의 각 픽셀을 처리
    for (int i = 0; i < imgHeight; i++) {
        for (int j = 0; j < imgWidth; j++) {
            double weightedSum = 0.0;
            double normalizationFactor = 0.0;

            // 중심 픽셀의 밝기값
            double centerPixel = paddingImg[i * paddingSize + j];

            // 필터의 윈도우 내의 픽셀에 대해 반복
            for (int di = -filterRadius; di <= filterRadius; di++) {
                for (int dj = -filterRadius; dj <= filterRadius; dj++) {
                    int ni = i + di;
                    int nj = j + dj;

                    // 경계 검사
                    if (ni >= 0 && ni < imgHeight && nj >= 0 && nj < imgWidth) {
                        double neighborPixel = paddingImg[ni * paddingSize + nj];
                        double spaceWeight = gaussian(sqrt(di * di + dj * dj), sigmaSpace);
                        double colorWeight = gaussian(centerPixel - neighborPixel, sigmaColor);

                        double weight = spaceWeight * colorWeight;
                        weightedSum += neighborPixel * weight;
                        normalizationFactor += weight;
                    }
                }
            }

            // 가중치에 따른 평균을 계산하여 출력 이미지 설정
            outputImg[i * imgWidth + j] = weightedSum / normalizationFactor;
        }
    }
}
  1. 와이너 필터 (Wiener Filter):
    • 통계적 접근을 사용하여 전체 이미지의 국소 평균과 분산을 기반으로 노이즈를 제거합니다.

void wienerFilter(double* outImg, double* paddedImg, int width, int height, int pheight, int pwidth) {
// 임시 버퍼 생성
double* tempImg = (double*)malloc(pwidth * pheight * sizeof(double));
int i, j, m, n;

// 필터 크기 설정
int filterSize = 3; // 3x3 필터를 예로 들겠습니다.
int halfSize = filterSize / 2;

// 노이즈 추정 (이 예제에서는 간단히 평균 노이즈를 사용합니다)
double noise_variance = 0.02; // 임의의 노이즈 분산값

// 이미지 전체를 반복
for (i = 0; i < height; i++) {
    for (j = 0; j < width; j++) {
        double sum = 0.0;
        double sum_weights = 0.0;

        // 로컬 패치에 대해 와이너 필터 계산
        for (m = -halfSize; m <= halfSize; m++) {
            for (n = -halfSize; n <= halfSize; n++) {
                int px = j + n;
                int py = i + m;

                // 경계 검사
                if (px >= 0 && px < pwidth && py >= 0 && py < pheight) {
                    double pixel = paddedImg[py * pwidth + px];
                    double weight = 1 / (noise_variance + fabs(pixel));
                    sum += pixel * weight;
                    sum_weights += weight;
                }
            }
        }

        // 정규화된 가중치로 출력 이미지 계산
        outImg[i * width + j] = sum / sum_weights;
    }
}

free(tempImg);

}

4. **안시로트로픽 디퓨전 필터 (Anisotropic Diffusion Filter)**:
   - 이미지의 에지를 유지하면서 노이즈를 제거합니다. 각 픽셀의 이웃과의 차이를 기반으로 확산 과정을 조절합니다.

double gradient(double x) {
return exp(-x * x / (2 * 0.1 * 0.1)); // Gaussian function, K=0.1
}

void anisotropicDiffusionFilter(double* outImg, double* paddedImg, int width, int height, int pheight, int pwidth) {
int i, j, iter;
int num_iterations = 15; // Number of diffusion iterations
double lambda = 0.25; // Diffusion coefficient, usually in the range 0.1-0.25
double* newImg = (double*)malloc(pwidth * pheight * sizeof(double));

// Initialize newImg from paddedImg
for (i = 0; i < pheight; i++) {
    for (j = 0; j < pwidth; j++) {
        newImg[i * pwidth + j] = paddedImg[i * pwidth + j];
    }
}

for (iter = 0; iter < num_iterations; iter++) {
    for (i = 1; i < height + 1; i++) {
        for (j = 1; j < width + 1; j++) {
            // Compute gradients
            double north = newImg[(i - 1) * pwidth + j] - newImg[i * pwidth + j];
            double south = newImg[(i + 1) * pwidth + j] - newImg[i * pwidth + j];
            double east = newImg[i * pwidth + (j + 1)] - newImg[i * pwidth + j];
            double west = newImg[i * pwidth + (j - 1)] - newImg[i * pwidth + j];

            // Update the image using anisotropic diffusion
            newImg[i * pwidth + j] += lambda * (
                gradient(north) * north + 
                gradient(south) * south +
                gradient(east) * east +
                gradient(west) * west
            );
        }
    }

    // Copy newImg back to paddedImg for the next iteration
    for (i = 1; i < height + 1; i++) {
        for (j = 1; j < width + 1; j++) {
            paddedImg[i * pwidth + j] = newImg[i * pwidth + j];
        }
    }
}

// Copy the internal image (without padding) to outImg
for (i = 0; i < height; i++) {
    for (j = 0; j < width; j++) {
        outImg[i * width + j] = newImg[(i + 1) * pwidth + (j + 1)];
    }
}

free(newImg);

}


설명된 주요 구성 요소:
그래디언트 함수: 이 함수는 강도의 그래디언트를 기반으로 확산 계수를 계산합니다. 종종 강도 그래디언트가 작은 영역을 매끄럽게 하면서 가장자리를 보존할 수 있는 지수 붕괴 함수로 정의됩니다.
람다(λ): 이것은 확산 계수입니다. 이것은 반복당 확산 과정이 얼마나 큰 영향을 미치는지를 제어합니다. 값이 작으면 확산 속도가 느려지고 그 반대의 경우도 있습니다.
반복: 이방성 확산은 반복되는 과정입니다. 반복 횟수에 따라 이미지가 얼마나 부드러워지는지가 결정됩니다. 반복 횟수가 많으면 이미지가 더 부드러워지지만 지나치게 반복되면 중요한 세부 사항이 흐려질 수 있습니다.
경계 조건: 이 구현은 에지 픽셀이 처리되지 않는 간단한 경계 조건을 가정하여 이미지 경계를 보존합니다.
이 기능은 기존 C 프로젝트에 통합되어 이미지를 처리하는 데 사용할 수 있으며 에지를 보존하면서 노이즈를 줄여 이미지의 중요한 구조적 세부 사항을 유지할 수 있습니다.


5. **논-로컬 민스 필터 (Non-Local Means Filter)**:
   - 이미지 내의 유사한 패치를 찾아 평균을 취함으로써 노이즈를 제거합니다. 이는 특히 질감이 복잡한 이미지에서 효과적입니다.

#define FILTER_SIZE 3
#define SEARCH_SIZE 21
#define SIGMA 10.0

double euclideanDistance(double* a, double* b, int size) {
double sum = 0.0;
for (int i = 0; i < size; i++) {
sum += (a[i] - b[i]) * (a[i] - b[i]);
}
return sqrt(sum);
}

double weight(double distance, double sigma) {
return exp(-(distance * distance) / (sigma * sigma));
}

void nonLocalMeansFilter(double* outImg, double* paddedImg, int width, int height, int pheight, int pwidth) {
int dx, dy, nx, ny, kx, ky;
double filterWindow[FILTER_SIZE * FILTER_SIZE];
double neighborWindow[FILTER_SIZE * FILTER_SIZE];
int halfFilterSize = FILTER_SIZE / 2;
int halfSearchSize = SEARCH_SIZE / 2;
double h = 2.0 * SIGMA * SIGMA;

for (int i = 0; i < height; i++) {
    for (int j = 0; j < width; j++) {
        double weightsSum = 0.0;
        double weightedSum = 0.0;
        double normFactor = 0.0;

        // Extract the local neighborhood around (i, j) in the padded image
        for (dy = -halfFilterSize; dy <= halfFilterSize; dy++) {
            for (dx = -halfFilterSize; dx <= halfFilterSize; dx++) {
                int imageIndex = (i + halfSearchSize + dy) * pwidth + (j + halfSearchSize + dx);
                filterWindow[(dy + halfFilterSize) * FILTER_SIZE + (dx + halfFilterSize)] = paddedImg[imageIndex];
            }
        }

        // Compare to all neighborhoods in the search window
        for (ny = -halfSearchSize; ny <= halfSearchSize; ny++) {
            for (nx = -halfSearchSize; nx <= halfSearchSize; nx++) {
                // Extract the neighbor window centered at (i + ny, j + nx)
                for (ky = -halfFilterSize; ky <= halfFilterSize; ky++) {
                    for (kx = -halfFilterSize; kx <= halfFilterSize; kx++) {
                        int neighborIndex = (i + halfSearchSize + ny + ky) * pwidth + (j + halfSearchSize + nx + kx);
                        neighborWindow[(ky + halfFilterSize) * FILTER_SIZE + (kx + halfFilterSize)] = paddedImg[neighborIndex];
                    }
                }

                double dist = euclideanDistance(filterWindow, neighborWindow, FILTER_SIZE * FILTER_SIZE);
                double w = weight(dist, h);
                int centralPixelIndex = (i + halfSearchSize + ny) * pwidth + (j + halfSearchSize + nx);
                weightedSum += w * paddedImg[centralPixelIndex];
                weightsSum += w;
            }
        }

        normFactor = weightsSum > 0 ? 1.0 / weightsSum : 0;
        outImg[i * width + j] = weightedSum * normFactor;
    }
}

}


6. **Total Variation Denoising**:
   - 이미지의 전체 변동을 최소화하면서 노이즈를 제거합니다. 이 기법은 이미지의 주요 구조를 유지하면서 스무딩합니다.

void totalVariationDenoising(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
double lambda = 0.1; // 정규화 파라미터
int iterations = 100; // 반복 횟수
double tau = 0.125; // 시간 간격 (stable choice for tau is 1/8)

// 출력 이미지에 초기 입력 이미지를 복사
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth; j++) {
        outputImg[i * imgWidth + j] = paddingImg[i * paddingWidth + j];
    }
}

// 임시 버퍼
double* newImg = (double*)malloc(imgHeight * imgWidth * sizeof(double));

// TV denoising 알고리즘 실행
for (int iter = 0; iter < iterations; iter++) {
    for (int i = 1; i < imgHeight - 1; i++) {
        for (int j = 1; j < imgWidth - 1; j++) {
            // 근접 픽셀에서의 차이 계산
            double dx = outputImg[i * imgWidth + (j + 1)] - outputImg[i * imgWidth + j];
            double dy = outputImg[(i + 1) * imgWidth + j] - outputImg[i * imgWidth + j];
            double dx_prev = outputImg[i * imgWidth + j] - outputImg[i * imgWidth + (j - 1)];
            double dy_prev = outputImg[i * imgWidth + j] - outputImg[(i - 1) * imgWidth + j];

            // Div(P) 계산
            double divP = (dx - dx_prev) + (dy - dy_prev);

            // 이미지 업데이트
            newImg[i * imgWidth + j] = outputImg[i * imgWidth + j] + tau * lambda * divP;
        }
    }

    // 버퍼 스왑
    double* temp = outputImg;
    outputImg = newImg;
    newImg = temp;
}

// 결과 복사
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth; j++) {
        outputImg[i * imgWidth + j] = newImg[i * imgWidth + j];
    }
}

// 메모리 해제
free(newImg);

}


7. **Wavelet Transform**:
   - 이미지를 다양한 주파수 성분으로 분해하고, 노이즈가 생길 가능성이 높은 성분만을 제거합니다.

void waveletTransform(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
// Haar Wavelet 변환에 사용되는 계수
const double sqrt2 = sqrt(2.0);

// 중간 계산을 위한 임시 버퍼
double* tempRow = (double*)malloc(sizeof(double) * paddingWidth);
double* tempCol = (double*)malloc(sizeof(double) * imgHeight);

// 각 행에 대한 웨이블릿 변환
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth / 2; j++) {
        int index1 = i * paddingWidth + 2 * j;
        int index2 = i * paddingWidth + 2 * j + 1;
        tempRow[2 * j] = (paddingImg[index1] + paddingImg[index2]) / sqrt2;     // 평균 (근사)
        tempRow[2 * j + 1] = (paddingImg[index1] - paddingImg[index2]) / sqrt2; // 차이 (세부)
    }
    memcpy(&outputImg[i * imgWidth], tempRow, imgWidth * sizeof(double));
}

// 각 열에 대한 웨이블릿 변환
for (int j = 0; j < imgWidth; j++) {
    for (int i = 0; i < imgHeight / 2; i++) {
        int index1 = 2 * i * imgWidth + j;
        int index2 = (2 * i + 1) * imgWidth + j;
        tempCol[2 * i] = (outputImg[index1] + outputImg[index2]) / sqrt2;        // 평균 (근사)
        tempCol[2 * i + 1] = (outputImg[index1] - outputImg[index2]) / sqrt2;   // 차이 (세부)
    }
    for (int i = 0; i < imgHeight; i++) {
        outputImg[i * imgWidth + j] = tempCol[i];
    }
}

free(tempRow);
free(tempCol);

}

void haarWaveletTransform(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingSize) {
const double sqrt2 = sqrt(2.0); // 분모에 sqrt(2)가 자주 사용되므로 미리 계산

// 각 행에 대한 웨이블릿 변환
for (int i = 0; i < imgHeight; i++) {
    int halfWidth = imgWidth / 2;
    for (int j = 0; j < halfWidth; j++) {
        int indexA = i * paddingSize + 2 * j; // 현재 위치
        int indexB = indexA + 1;             // 바로 다음 위치
        outputImg[i * imgWidth + j] = (paddingImg[indexA] + paddingImg[indexB]) / sqrt2;       // 평균 (근사)
        outputImg[i * imgWidth + halfWidth + j] = (paddingImg[indexA] - paddingImg[indexB]) / sqrt2; // 차이 (세부)
    }
}

// 임시 배열을 사용하여 각 열에 대해 변환 수행
double* tempCol = (double*)malloc(sizeof(double) * imgHeight);

for (int j = 0; j < imgWidth; j++) {
    int halfHeight = imgHeight / 2;
    for (int i = 0; i < halfHeight; i++) {
        int indexA = 2 * i * imgWidth + j;
        int indexB = (2 * i + 1) * imgWidth + j;
        tempCol[i] = (outputImg[indexA] + outputImg[indexB]) / sqrt2;  // 평균 (근사)
        tempCol[halfHeight + i] = (outputImg[indexA] - outputImg[indexB]) / sqrt2;  // 차이 (세부)
    }
    for (int i = 0; i < imgHeight; i++) {
        outputImg[i * imgWidth + j] = tempCol[i];
    }
}

free(tempCol);

}

void waveletTransform(double* output, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
// 가정: output은 이미 메모리가 할당되어 있고, 충분한 크기를 가지고 있음.
// 이 함수에서는 Haar 웨이블릿을 사용하여 간단한 2D 변환을 수행합니다.
int halfWidth = imgWidth / 2;
int halfHeight = imgHeight / 2;

// 임시 버퍼 생성
double* tempRow = (double*)malloc(sizeof(double) * paddingWidth);
double* tempCol = (double*)malloc(sizeof(double) * imgHeight);

// 각 행에 대해 1D 웨이블릿 변환 수행
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < halfWidth; j++) {
        // 저주파수 (평균)
        tempRow[2 * j] = (paddingImg[i * paddingWidth + 2 * j] + paddingImg[i * paddingWidth + 2 * j + 1]) / sqrt(2);
        // 고주파수 (차이)
        tempRow[2 * j + 1] = (paddingImg[i * paddingWidth + 2 * j] - paddingImg[i * paddingWidth + 2 * j + 1]) / sqrt(2);
    }
    memcpy(paddingImg + i * paddingWidth, tempRow, sizeof(double) * paddingWidth);
}

// 각 열에 대해 1D 웨이블릿 변환 수행
for (int j = 0; j < imgWidth; j++) {
    for (int i = 0; i < halfHeight; i++) {
        // 저주파수 (평균)
        tempCol[2 * i] = (paddingImg[2 * i * paddingWidth + j] + paddingImg[(2 * i + 1) * paddingWidth + j]) / sqrt(2);
        // 고주파수 (차이)
        tempCol[2 * i + 1] = (paddingImg[2 * i * paddingWidth + j] - paddingImg[(2 * i + 1) * paddingWidth + j]) / sqrt(2);
    }
    for (int i = 0; i < imgHeight; i++) {
        output[i * paddingWidth + j] = tempCol[i];
    }
}

free(tempRow);
free(tempCol);

}


8. **Guided Filter**:
   - 가이드 이미지를 사용하여 주 이미지의 필터링을 수행합니다. 이는 에지 보존에 탁월하며, 디테일과 에지를 유지하면서 스무딩을 수행합니다.

void guidedFilter(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
int r = 4; // 지역 윈도우 반경
double eps = 0.0001; // 정규화 파라미터

// 지역 평균을 계산하기 위한 임시 배열
double* meanI = (double*)malloc(sizeof(double) * imgHeight * imgWidth);
double* meanP = (double*)malloc(sizeof(double) * imgHeight * imgWidth);
double* corrI = (double*)malloc(sizeof(double) * imgHeight * imgWidth);
double* corrIP = (double*)malloc(sizeof(double) * imgHeight * imgWidth);

double* varI = (double*)malloc(sizeof(double) * imgHeight * imgWidth);  // 입력의 분산
double* covIP = (double*)malloc(sizeof(double) * imgHeight * imgWidth); // 입력과 출력의 공분산

double* a = (double*)malloc(sizeof(double) * imgHeight * imgWidth);     // 선형 계수
double* b = (double*)malloc(sizeof(double) * imgHeight * imgWidth);     // 선형 계수

// 임시 배열 초기화
memset(meanI, 0, sizeof(double) * imgHeight * imgWidth);
memset(meanP, 0, sizeof(double) * imgHeight * imgWidth);
memset(corrI, 0, sizeof(double) * imgHeight * imgWidth);
memset(corrIP, 0, sizeof(double) * imgHeight * imgWidth);

// 원본 이미지의 평균 및 자기상관을 계산
for (int i = r; i < imgHeight - r; i++) {
    for (int j = r; j < imgWidth - r; j++) {
        double sumI = 0, sumP = 0, sumII = 0, sumIP = 0;
        for (int k = -r; k <= r; k++) {
            for (int l = -r; l <= r; l++) {
                double I = paddingImg[(i + k) * paddingWidth + (j + l)];
                double P = paddingImg[(i + k) * paddingWidth + (j + l)]; // P와 I가 같은 경우
                sumI += I;
                sumP += P;
                sumII += I * I;
                sumIP += I * P;
            }
        }
        int size = (2 * r + 1) * (2 * r + 1);
        meanI[i * imgWidth + j] = sumI / size;
        meanP[i * imgWidth + j] = sumP / size;
        corrI[i * imgWidth + j] = sumII / size;
        corrIP[i * imgWidth + j] = sumIP / size;
    }
}

// 입력의 분산과 입력-출력의 공분산을 계산
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth; j++) {
        varI[i * imgWidth + j] = corrI[i * imgWidth + j] - meanI[i * imgWidth + j] * meanI[i * imgWidth + j];
        covIP[i * imgWidth + j] = corrIP[i * imgWidth + j] - meanI[i * imgWidth + j] * meanP[i * imgWidth + j];
        a[i * imgWidth + j] = covIP[i * imgWidth + j] / (varI[i * imgWidth + j] + eps);
        b[i * imgWidth + j] = meanP[i * imgWidth + j] - a[i * imgWidth + j] * meanI[i * imgWidth + j];
    }
}

// 선형 모델을 사용하여 출력을 계산
for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth; j++) {
        outputImg[i * imgWidth + j] = a[i * imgWidth + j] * paddingImg[i * paddingWidth + j] + b[i * imgWidth + j];
    }
}

// 할당된 메모리 해제
free(meanI);
free(meanP);
free(corrI);
free(corrIP);
free(varI);
free(covIP);
free(a);
free(b);

}


9. **Adaptive Median Filter**:
   - 필터의 크기를 동적으로 조절하여, 노이즈의 양에 따라 최적의 결과를 제공합니다.

// 중앙값을 찾는 함수
double median(double* arr, int size) {
for (int i = 0; i < size - 1; i++) {
for (int j = i + 1; j < size; j++) {
if (arr[i] > arr[j]) {
double temp = arr[i];
arr[i] = arr[j];
arr[j] = temp;
}
}
}
if (size % 2 == 0)
return (arr[size / 2 - 1] + arr[size / 2]) / 2;
else
return arr[size / 2];
}

// Adaptive Median Filter 함수
void adaptiveMedianFilter(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
int maxWindowSize = 7; // 최대 윈도우 크기 설정

for (int i = 0; i < imgHeight; i++) {
    for (int j = 0; j < imgWidth; j++) {
        int windowSize = 3;
        double Zmed = 0;
        double Zmin = 0;
        double Zmax = 0;
        double Zxy = paddingImg[i * paddingWidth + j];
        double Smax = 0;

        while (windowSize <= maxWindowSize) {
            int k = windowSize / 2;
            int count = 0;
            double* window = (double*)malloc(sizeof(double) * windowSize * windowSize);

            Zmin = INFINITY;
            Zmax = -INFINITY;
            for (int m = -k; m <= k; m++) {
                for (int n = -k; n <= k; n++) {
                    int idxX = i + m;
                    int idxY = j + n;
                    // 경계 처리
                    if (idxX < 0) idxX = 0;
                    if (idxY < 0) idxY = 0;
                    if (idxX >= imgHeight) idxX = imgHeight - 1;
                    if (idxY >= imgWidth) idxY = imgWidth - 1;

                    double pixelValue = paddingImg[idxX * paddingWidth + idxY];
                    window[count++] = pixelValue;
                    if (pixelValue < Zmin) Zmin = pixelValue;
                    if (pixelValue > Zmax) Zmax = pixelValue;
                }
            }

            Zmed = median(window, count);
            Smax = Zmax - Zmin;

            // A1 조건
            if (Zmed > Zmin && Zmed < Zmax) {
                // A2 조건
                if (Zxy > Zmin && Zxy < Zmax) {
                    outputImg[i * imgWidth + j] = Zxy;
                } else {
                    outputImg[i * imgWidth + j] = Zmed;
                }
                free(window);
                break;
            }
            windowSize += 2; // 윈도우 크기 증가
            free(window);
        }

        // 최대 윈도우 크기에 도달했는데 A1이 충족되지 않는 경우
        if (windowSize > maxWindowSize) {
            outputImg[i * imgWidth + j] = Zmed;
        }
    }
}

}


10. **Recursive Least Squares Filter (RLS Filter)**:
    - 이미지의 노이즈 제거를 위해 재귀적 최소 제곱 추정을 사용합니다. 이 필터는 동적 시스템에서 주로 사용되지만, 이미지에서 노이즈 제거에도 적용할 수 있습니다.

void rlsFilter(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingWidth) {
double lambda = 0.99; // forgetting factor, lambda는 일반적으로 0.9에서 1 사이의 값
double delta = 0.01; // 초기화 파라미터

// 초기화
double** P = (double**)malloc(imgHeight * sizeof(double*));
for (int i = 0; i < imgHeight; i++) {
    P[i] = (double*)malloc(imgWidth * sizeof(double));
    for (int j = 0; j < imgWidth; j++) {
        P[i][j] = 1.0 / delta; // 대각 요소를 delta로 초기화
    }
}

// 필터링을 위한 버퍼
double* y_est = (double*)calloc(imgHeight * imgWidth, sizeof(double)); // 추정값
double* e = (double*)calloc(imgHeight * imgWidth, sizeof(double)); // 오차값

// RLS 필터 알고리즘
for (int i = 1; i < imgHeight - 1; i++) {
    for (int j = 1; j < imgWidth - 1; j++) {
        // 이웃 픽셀에서 입력을 추출
        double y = paddingImg[i * paddingWidth + j]; // 현재 픽셀 값
        double y_n = paddingImg[(i-1) * paddingWidth + j]; // 이전 픽셀 값
        double y_est_n = y_est[(i-1) * imgWidth + j]; // 이전 추정값

        // 이전 오차 업데이트
        e[(i-1) * imgWidth + j] = y_n - y_est_n;

        // 게인 계산
        double k = P[i][j] * y_n / (lambda + y_n * P[i][j] * y_n);

        // 추정값 업데이트
        y_est[i * imgWidth + j] = y_est_n + k * e[(i-1) * imgWidth + j];

        // 오차 계산
        e[i * imgWidth + j] = y - y_est[i * imgWidth + j];

        // P 매트릭스 업데이트
        P[i][j] = (1.0 / lambda) * (P[i][j] - k * y_n * P[i][j]);

        // 출력 이미지에 결과 저장
        outputImg[i * imgWidth + j] = y_est[i * imgWidth + j];
    }
}

// 메모리 해제
for (int i = 0; i < imgHeight; i++) {
    free(P[i]);
}
free(P);
free(y_est);
free(e);

}


여긴 upsampling

영상처리에서 upsampling을 수행하는 다양한 방식이 있습니다. 아래에 10가지 추천드리겠습니다:

Nearest Neighbor (최근접 이웃) 방식: 가장 가까운 픽셀 값을 사용하여 이미지를 확대합니다.

void upsampling(int ratio, double* org, double* ups, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
int width = bmpInfo.biWidth, height = bmpInfo.biHeight;
int width2 = bmpInfo.biWidth * ratio, height2 = bmpInfo.biHeight * ratio;
int stride2 = (((bmpInfo.biBitCount / 8) * width2) + 3) / 4 * 4, size2 = stride2 * height2;
int nearestJ = 0, nearestI=0;
for (int j = 0; j < height2; j++) {
for (int i = 0; i < width2; i++) {
nearestJ = j / ratio;
nearestI = i / ratio;
ups[j * width2 + i] = org[nearestJ * width + nearestI];
}
}
bmpInfo.biWidth = width2;
bmpInfo.biHeight = height2;
bmpInfo.biSizeImage = size2;
bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size2;
makeOutFile(bmpFile, bmpInfo, "OutUPS.bmp", ups);
}

Bilinear Interpolation (양선형 보간): 주변 4개 픽셀의 가중 평균을 사용하여 이미지를 확대합니다.

void bilinearUpsample(double* inputImg, int imgHeight, int imgWidth, int ratio, double** outputImg) {
int outputHeight = imgHeight * ratio;
int outputWidth = imgWidth * ratio;

// 메모리 할당
*outputImg = (double*)malloc(outputHeight * outputWidth * sizeof(double));

// 업샘플링 수행
for (int i = 0; i < outputHeight; i++) {
    for (int j = 0; j < outputWidth; j++) {
        // 원본 이미지에서의 위치 계산
        float x = (float)j / ratio;
        float y = (float)i / ratio;

        int x1 = (int)x;
        int y1 = (int)y;
        int x2 = x1 + 1 < imgWidth ? x1 + 1 : x1;
        int y2 = y1 + 1 < imgHeight ? y1 + 1 : y1;

        // 각 방향으로의 비율 계산
        float x2_ratio = x - x1;
        float x1_ratio = 1 - x2_ratio;
        float y2_ratio = y - y1;
        float y1_ratio = 1 - y2_ratio;

        // Bilinear interpolation
        double result = inputImg[y1 * imgWidth + x1] * x1_ratio * y1_ratio +
                        inputImg[y1 * imgWidth + x2] * x2_ratio * y1_ratio +
                        inputImg[y2 * imgWidth + x1] * x1_ratio * y2_ratio +
                        inputImg[y2 * imgWidth + x2] * x2_ratio * y2_ratio;

        (*outputImg)[i * outputWidth + j] = result;
    }
}

}


Bicubic Interpolation (바이큐빅 보간): 주변 16개 픽셀의 가중 평균을 사용하여 이미지를 확대합니다.

double cubicInterpolate(double p[4], double x) {
return p[1] + 0.5 * x * (p[2] - p[0] + x * (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] + x * (3.0 * (p[1] - p[2]) + p[3] - p[0])));
}

double bicubicInterpolate(double p[4][4], double x, double y) {
double arr[4];
arr[0] = cubicInterpolate(p[0], y);
arr[1] = cubicInterpolate(p[1], y);
arr[2] = cubicInterpolate(p[2], y);
arr[3] = cubicInterpolate(p[3], y);
return cubicInterpolate(arr, x);
}

void bicubicUpsample(double* outputImg, double* paddingImg, int imgHeight, int imgWidth, int paddingSize, int scaleFactor) {
int outputHeight = imgHeight * scaleFactor;
int outputWidth = imgWidth * scaleFactor;

for (int i = 0; i < outputHeight; ++i) {
    for (int j = 0; j < outputWidth; ++j) {
        double x = (double)j / scaleFactor;
        double y = (double)i / scaleFactor;
        int fx = (int)floor(x);
        int fy = (int)floor(y);
        double x_diff = x - fx;
        double y_diff = y - fy;

        double mat[4][4];
        for (int m = -1; m <= 2; ++m) {
            for (int n = -1; n <= 2; ++n) {
                int px = fx + n;
                int py = fy + m;
                if (px < 0) px = 0;
                if (py < 0) py = 0;
                if (px >= imgWidth) px = imgWidth - 1;
                if (py >= imgHeight) py = imgHeight - 1;
                mat[m + 1][n + 1] = paddingImg[py * paddingSize + px];
            }
        }

        outputImg[i * outputWidth + j] = bicubicInterpolate(mat, x_diff, y_diff);
    }
}

}

double cubicInterpolate(double p[4], double x) {
return p[1] + 0.5 * x * (p[2] - p[0] + x * (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] + x * (3.0 * (p[1] - p[2]) + p[3] - p[0])));
}

double bicubicInterpolate(double p[4][4], double x, double y) {
double arr[4];
arr[0] = cubicInterpolate(p[0], y);
arr[1] = cubicInterpolate(p[1], y);
arr[2] = cubicInterpolate(p[2], y);
arr[3] = cubicInterpolate(p[3], y);
return cubicInterpolate(arr, x);
}

double* bicubicUpsample(double* inputImg, int imgHeight, int imgWidth, int ratio) {
int outputHeight = imgHeight * ratio;
int outputWidth = imgWidth * ratio;
double* outputImg = (double*)malloc(outputHeight * outputWidth * sizeof(double));

for (int i = 0; i < outputHeight; ++i) {
    for (int j = 0; j < outputWidth; ++j) {
        double x = (double)j / ratio;
        double y = (double)i / ratio;
        int fx = (int)floor(x);
        int fy = (int)floor(y);
        double x_diff = x - fx;
        double y_diff = y - fy;

        double mat[4][4];
        for (int m = -1; m <= 2; ++m) {
            for (int n = -1; n <= 2; ++n) {
                int px = fx + n;
                int py = fy + m;
                if (px < 0) px = 0;
                if (py < 0) py = 0;
                if (px >= imgWidth) px = imgWidth - 1;
                if (py >= imgHeight) py = imgHeight - 1;
                mat[m + 1][n + 1] = inputImg[py * imgWidth + px];
            }
        }

        outputImg[i * outputWidth + j] = bicubicInterpolate(mat, x_diff, y_diff);
    }
}

return outputImg;

}

Wavelet Transform (웨이블릿 변환): 다중 해상도로 분해하여 이미지를 확대합니다.

void inverseHaar1D(double* line, int length) {
double* temp = (double*)malloc(length * sizeof(double));
int halfLength = length / 2;

for (int i = 0; i < halfLength; i++) {
    temp[2*i] = (line[i] + line[i + halfLength]) / sqrt(2.0);
    temp[2*i + 1] = (line[i] - line[i + halfLength]) / sqrt(2.0);
}

for (int i = 0; i < length; i++) {
    line[i] = temp[i];
}

free(temp);

}

void inverseHaar2D(double* img, int width, int height) {
// Process each row
for (int i = 0; i < height; i++) {
inverseHaar1D(img + i * width, width);
}

// Process each column
double* tempCol = (double*)malloc(height * sizeof(double));

for (int j = 0; j < width; j++) {
    for (int i = 0; i < height; i++) {
        tempCol[i] = img[i * width + j];
    }

    inverseHaar1D(tempCol, height);

    for (int i = 0; i < height; i++) {
        img[i * width + j] = tempCol[i];
    }
}

free(tempCol);

}

double* upsamplingWithHaarWavelet(double* img, int width, int height, int ratio) {
int newWidth = width * ratio;
int newHeight = height * ratio;
double* upsampledImg = (double*)calloc(newWidth * newHeight, sizeof(double));

// Initial copy of image to the top-left corner
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        upsampledImg[y * newWidth + x] = img[y * width + x];
    }
}

// Apply inverse Haar Wavelet Transform
inverseHaar2D(upsampledImg, newWidth, newHeight);

return upsampledImg;

}

Edge-Aware Upsampling: 엣지 정보를 고려하여 이미지를 확대합니다.

double weight(double colorDiff, double spatialDist, double sigma) {
return exp(-(colorDiff * colorDiff) / (2 * sigma * sigma) - spatialDist * spatialDist / (2 * sigma * sigma));
}

void edgeAwareUpsample(double* inputImg, double* outputImg, int imgHeight, int imgWidth, int ratio, double sigma) {
int outputHeight = imgHeight * ratio;
int outputWidth = imgWidth * ratio;

for (int i = 0; i < outputHeight; ++i) {
    for (int j = 0; j < outputWidth; ++j) {
        double x = (double)j / ratio;
        double y = (double)i / ratio;
        int fx = (int)floor(x);
        int fy = (int)floor(y);
        double x_diff = x - fx;
        double y_diff = y - fy;

        int x2 = fx + 1 < imgWidth ? fx + 1 : fx;
        int y2 = fy + 1 < imgHeight ? fy + 1 : fy;

        double Q11 = inputImg[fy * imgWidth + fx];
        double Q12 = inputImg[y2 * imgWidth + fx];
        double Q21 = inputImg[fy * imgWidth + x2];
        double Q22 = inputImg[y2 * imgWidth + x2];

        double weights[4];
        weights[0] = weight(0.0, sqrt(x_diff * x_diff + y_diff * y_diff), sigma);
        weights[1] = weight(fabs(Q12 - Q11), sqrt(x_diff * x_diff + (1.0 - y_diff) * (1.0 - y_diff)), sigma);
        weights[2] = weight(fabs(Q21 - Q11), sqrt((1.0 - x_diff) * (1.0 - x_diff) + y_diff * y_diff), sigma);
        weights[3] = weight(fabs(Q22 - Q11), sqrt((1.0 - x_diff) * (1.0 - x_diff) + (1.0 - y_diff) * (1.0 - y_diff)), sigma);

        double totalWeight = weights[0] + weights[1] + weights[2] + weights[3];
        outputImg[i * outputWidth + j] = (Q11 * weights[0] + Q12 * weights[1] + Q21 * weights[2] + Q22 * weights[3]) / totalWeight;
    }
}

}


영상처리에서 사람의 시각적 만족도를 높이기 위해 사용되는 기법은 종종 전통적인 성능 지표(예: PSNR, MSE)와는 독립적으로 작동합니다. 아래는 영상을 보다 시각적으로 개선하기 위한 다양한 기법들을 소개합니다:

1. **히스토그램 평활화 (Histogram Equalization, HE)**:
   - 히스토그램 평활화는 이미지의 콘트라스트를 개선하여 어두운 영역과 밝은 영역의 디테일을 더 잘 볼 수 있도록 합니다.

void CumulativeHisto(int height, int width, double* y, int size, int stride, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
double* histogram;
histogram = (double*)calloc(256, sizeof(double));
int width2 = 512, height2 = 512, stride2 = (((bmpInfo.biBitCount / 8) * width2) + 3) / 4 * 4, size2 = stride2 * height2;
printf("%d %d %d %d\n", width2, height2, size2, stride2);

// get Histogram
for (int j = 0; j < height; j++) {
    for (int i = 0; i < width; i++) {
        histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))]++;
    }
}
for (int i = 0; i < 255; i++) {
    histogram[i + 1] += histogram[i];
}

// out Histogram
double histoMax = 0;
for (int i = 0; i < 256; i++) {
    if (histogram[i] > histoMax) histoMax = histogram[i];
}
unsigned char* outHisto, * HE;
outHisto = (unsigned char*)calloc(size2, sizeof(unsigned char));
HE = (unsigned char*)calloc(size, sizeof(unsigned char));


for (int j = 0; j < height; j++) {
    for (int i = 0; i < width; i++) {
        HE[j * stride + 3 * i + 0] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
        HE[j * stride + 3 * i + 1] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
        HE[j * stride + 3 * i + 2] = (unsigned char)(255 * histogram[(unsigned char)(y[j * width + i] > 255 ? 255 : (y[j * width + i] < 0 ? 0 : y[j * width + i]))] / (width * height));
    }
}
for (int i = 0; i < 256; i++)histogram[i] = histogram[i] / histoMax * 450;

for (int j = 0; j < height2; j++) {
    for (int i = 0; i < width2; i++) {
        if (histogram[i / 2] > j) continue;
        outHisto[j * stride2 + 3 * i + 0] = 255;
        outHisto[j * stride2 + 3 * i + 1] = 255;
        outHisto[j * stride2 + 3 * i + 2] = 255;
    }
}


FILE* outputHisto = fopen("OutHisto.bmp", "wb");
bmpInfo.biWidth = 512;
bmpInfo.biHeight = 512;
bmpInfo.biSizeImage = 512 * 3 * 512;
bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 512 * 3 * 512;


fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputHisto);
fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputHisto);
fwrite(outHisto, sizeof(unsigned char), size2, outputHisto);

FILE* outputHE = fopen("OutHE.bmp", "wb");

bmpInfo.biWidth = width;
bmpInfo.biHeight = height;
bmpInfo.biSizeImage = size;
bmpFile.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + size;

fwrite(&bmpFile, sizeof(BITMAPFILEHEADER), 1, outputHE);
fwrite(&bmpInfo, sizeof(BITMAPINFOHEADER), 1, outputHE);
fwrite(HE, sizeof(unsigned char), size, outputHE);
fclose(outputHE);
free(HE);



free(histogram);
fclose(outputHisto);
free(outHisto);

}


2. **적응형 히스토그램 평활화 (Adaptive Histogram Equalization, AHE)**:
   - AHE는 이미지를 여러 작은 블록으로 나누고 각 블록에 대해 히스토그램 평활화를 수행하여 국소적인 콘트라스트를 향상시킵니다.

void adaptiveHistogramEqualization(double* inputImg, int height, int width, double* outputImg, int tileSize) {
int numTilesX = (width + tileSize - 1) / tileSize; // X축 타일 수 계산
int numTilesY = (height + tileSize - 1) / tileSize; // Y축 타일 수 계산
int tileArea = tileSize * tileSize;

for (int ty = 0; ty < numTilesY; ty++) {
    for (int tx = 0; tx < numTilesX; tx++) {
        int startX = tx * tileSize;
        int startY = ty * tileSize;
        int endX = startX + tileSize > width ? width : startX + tileSize;
        int endY = startY + tileSize > height ? height : startY + tileSize;

        double* histogram = calloc(256, sizeof(double));
        double* cdf = calloc(256, sizeof(double));

        // Build histogram for the current tile
        for (int y = startY; y < endY; y++) {
            for (int x = startX; x < endX; x++) {
                int pixelValue = (int)inputImg[y * width + x];
                histogram[pixelValue]++;
            }
        }

        // Compute the Cumulative Distribution Function (CDF)
        double cdfMin = 0;
        double scale = 255.0 / (tileArea - cdfMin);
        for (int i = 0; i < 256; i++) {
            cdf[i] = (i == 0) ? histogram[i] : (cdf[i - 1] + histogram[i]);
            if (cdf[i] != 0 && cdfMin == 0) cdfMin = cdf[i];
            cdf[i] = (cdf[i] - cdfMin) * scale;
        }

        // Map the image using the CDF
        for (int y = startY; y < endY; y++) {
            for (int x = startX; x < endX; x++) {
                int pixelValue = (int)inputImg[y * width + x];
                outputImg[y * width + x] = cdf[pixelValue];
            }
        }

        free(histogram);
        free(cdf);
    }
}

}

이 코드는 이미지를 tileSize로 정의된 크기의 타일로 나누고, 각 타일에 대해 독립적인 히스토그램 평활화를 수행합니다. 각 픽셀은 해당 타일의 누적 분포 함수를 사용하여 새로운 값으로 매핑됩니다. 이 접근 방법은 이미지의 다양한 부분에서 지역적인 콘트라스트를 개선하며, 특히 균일하지 않은 조명이나 대비가 있는 영역에서 유용합니다.

tileSize는 처리할 타일의 크기를 결정하며, 이 값에 따라 콘트라스트 향상의 정도와 세밀함이 달라집니다. 타일이 작을수록 더 세밀한 콘트라스트 조정이 가능하지만, 너무 작으면 노이즈가 증가할 수 있습니다.


3. **Contrast Limited Adaptive Histogram Equalization (CLAHE)**:
   - CLAHE는 AHE의 변형으로, 콘트라스트가 지나치게 높아지는 것을 방지하기 위해 히스토그램 평활화 시 콘트라스트 제한을 둡니다.

CLAHE (Contrast Limited Adaptive Histogram Equalization)는 표준 AHE(Adaptive Histogram Equalization)의 일반적인 문제인 콘트라스트 증가를 완화하기 위해 개발되었습니다. CLAHE는 지역적으로 콘트라스트를 개선하지만, 특정 임계값을 넘지 않도록 하여 노이즈의 증가를 제한합니다.

CLAHE를 구현하기 위해서는 각 타일의 히스토그램을 클리핑한 후 적용하는 방식을 사용합니다. 아래는 이러한 CLAHE의 C 코드 구현 예시입니다.

```c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

void clipHistogram(double* histogram, int numPixels, double clipLimit) {
    double excess = 0;
    for (int i = 0; i < 256; i++) {
        if (histogram[i] > clipLimit) {
            excess += histogram[i] - clipLimit;
            histogram[i] = clipLimit;
        }
    }

    // Redistribute excess
    double increment = excess / 256.0;
    for (int i = 0; i < 256; i++) {
        histogram[i] += increment;
    }
}

void clahe(double* inputImg, int height, int width, double* outputImg, int tileSize, double clipLimit) {
    int numTilesX = (width + tileSize - 1) / tileSize;
    int numTilesY = (height + tileSize - 1) / tileSize;
    int numPixelsPerTile = tileSize * tileSize;

    for (int ty = 0; ty < numTilesY; ty++) {
        for (int tx = 0; tx < numTilesX; tx++) {
            int startX = tx * tileSize;
            int startY = ty * tileSize;
            int endX = startX + tileSize > width ? width : startX + tileSize;
            int endY = startY + tileSize > height ? height : startY + tileSize;

            double* histogram = calloc(256, sizeof(double));
            double* cdf = calloc(256, sizeof(double));

            // Build histogram
            for (int y = startY; y < endY; y++) {
                for (int x = startX; x < endX; x++) {
                    int pixelValue = (int)inputImg[y * width + x];
                    histogram[pixelValue]++;
                }
            }

            // Clip histogram and normalize
            clipHistogram(histogram, numPixelsPerTile, clipLimit);

            // Calculate CDF
            double scale = 255.0 / (numPixelsPerTile);
            for (int i = 0; i < 256; i++) {
                cdf[i] = (i == 0) ? histogram[i] : (cdf[i - 1] + histogram[i]);
                cdf[i] = (cdf[i] - cdf[0]) * scale;
            }

            // Apply CLAHE
            for (int y = startY; y < endY; y++) {
                for (int x = startX; x < endX; x++) {
                    int pixelValue = (int)inputImg[y * width + x];
                    outputImg[y * width + x] = cdf[pixelValue];
                }
            }

            free(histogram);
            free(cdf);
        }
    }
}

int main() {
    int width = 512;
    int height = 512;
    double* inputImg = calloc(width * height, sizeof(double));
    double* outputImg = calloc(width * height, sizeof(double));
    int tileSize = 32;
    double clipLimit = 40.0; // 클리핑 임계값

    // 여기에서 inputImg 초기화

    clahe(inputImg, height, width, outputImg, tileSize, clipLimit);

    // 결과 outputImg 사용

    free(inputImg);
    free(outputImg);

    return 0;
}

이 코드는 각 타일의 히스토그램을 계산하고, 클리핑한 후, 균일하게 재분배하여 타일에 대한 히스토그램 평활화를 수행합니다. clipLimit은 타일당 허용할 최대 픽셀 수를 결정하며, 이 값을 조절함으로써 결과 이미지의 콘트라스트를 세밀하게 조정할 수 있습니다. 이 기법은 의료 영상, 위성 영상 등에서 세밀한 텍스처를 강조할 필요가 있을 때 특히 유용합니다.

  1. 감마 보정 (Gamma Correction):
    • 감마 보정은 이미지의 밝기를 비선형적으로 조정하여, 사람의 눈이 밝기를 인식하는 방식을 모방합니다.

void GammaOut(double* y, double g, BITMAPFILEHEADER bmpFile, BITMAPINFOHEADER bmpInfo) {
int width = bmpInfo.biWidth, height = bmpInfo.biHeight, size = bmpInfo.biSizeImage;
double* gy;
gy = (double*)calloc(width * height, sizeof(double));

for (int j = 0; j < height; j++) {
    for (int i = 0; i < width; i++) {
        gy[j * width + i] = 255 * pow(y[j * width + i] / 255.0, 1/g);
    }
}
makeOutFile(bmpFile, bmpInfo, "gammaOut.bmp", gy);
free(gy);

}

```
5. 양방향 필터 (Bilateral Filter):

  • 양방향 필터는 엣지를 보존하면서 노이즈를 감소시키는 기법으로, 이미지의 세밀한 부분을 더욱 선명하게 만들어 줍니다.
  1. Guided Filter:

    • Guided Filter는 edge-aware smoothing을 제공하며, 이미지의 가장자리를 유지하면서 부드럽게 처리합니다.
  2. Retinex Theory:

    • Retinex 이론은 이미지의 색상을 조정하여 인간의 시각 시스템이 빛의 변화에 적응하는 방식을 모방하여 이미지의 색상을 더 자연스럽게 만듭니다.
  3. High Dynamic Range (HDR) 톤 매핑:

    • HDR 톤 매핑은 다양한 노출의 이미지들을 조합하여 넓은 동적 범위를 갖는 이미지를 생성하며, 세부적인 텍스처와 색상을 강화합니다.
  4. 디헤이징 (Dehazing):

    • 디헤이징은 안개나 흐림 효과를 줄이기 위해 이미지에서 불필요한 흐림을 제거합니다.
  5. 샤프닝 (Sharpening):

    • 샤프닝 필터는 이미지의 경계를 더욱 뚜렷하게 만들어 주어, 전체적인 이미지의 선명도를 향상시킵니다.

이러한 기법들은 사진이나 비디오에서 시각적 품질을 향상시키는 데 주로 사용되며, 특히 저조도 또는 고대비 환경에서 시각적 경험을 크게 개선할 수 있습니다. 각 기법은 특정 상황이나 요구에 따라 적절히 선택되어 사용됩니다.

728x90

'언어 > C' 카테고리의 다른 글

C언어 영상처리 Edge 구하기, threshold 구하기  (2) 2024.05.02
C언어 영상처리 - Rotation  (0) 2024.04.30
C언어 영상 처리 - Gamma  (1) 2024.04.18
C언어 영상처리 예정  (0) 2024.04.17
C언어 영상처리 다양한 필터  (0) 2024.04.17