프로그램 사용/octave

octave로 공간 fft 돌려보기

구차니 2026. 5. 19. 16:19

귀찮으니(!) gpt를 통해서 딸깍~ 시전

% 이미지 읽기
img = imread("test.jpg");

% RGB -> grayscale 직접 변환
if ndims(img) == 3

    r = double(img(:,:,1));
    g = double(img(:,:,2));
    b = double(img(:,:,3));

    gray = 0.299*r + 0.587*g + 0.114*b;
else
    gray = double(img);
end

% 블럭 크기
blockSize = 32;

[h, w] = size(gray);

% 결과 저장
scoreMap = zeros(floor(h/blockSize), floor(w/blockSize));

% 블럭별 FFT 처리
for by = 1:blockSize:(h - blockSize + 1)
    for bx = 1:blockSize:(w - blockSize + 1)

        % 블럭 추출
        block = gray(by:by+blockSize-1, bx:bx+blockSize-1);

        % 평균 제거
        block = block - mean(block(:));

        % FFT
        F = fft2(block);

        % magnitude
        mag = abs(F);

        % 중앙 이동
        mag = fftshift(mag);

        center = floor(blockSize/2) + 1;

        radius = 4;

        [X, Y] = meshgrid(1:blockSize, 1:blockSize);

        dist = sqrt((X-center).^2 + (Y-center).^2);

        % 고주파만 사용
        mask = dist > radius;

        % 고주파 에너지
        energy = sum(mag(mask).^2);

        % diversity score
        score = log(1 + energy);

        iy = (by-1)/blockSize + 1;
        ix = (bx-1)/blockSize + 1;

        scoreMap(iy, ix) = score;
    end
end

% 정규화
scoreMap = scoreMap - min(scoreMap(:));
scoreMap = scoreMap / max(scoreMap(:));

% 결과 표시
figure;
imshow(uint8(gray));
title("Original");

figure;
imagesc(scoreMap);
axis image;
colorbar;
title("FFT Diversity Map");

 

원본

그레이스케일

 

원인은 모르겠으나 blocksize 5 이하로는 검은색으로 나온다.

 

blocksize = 6 blocksize = 8 blocksize = 10
blocksize = 12 blocksize = 16 blocksize = 32

 

원본

 

그레이스케일

 

blocksize = 6 blocksize = 10 blocksize = 12

blocksize = 16 blocksize = 18 blocksize = 20
blocksize = 24 blocksize = 28 blocksize = 32

 

이미지 폭이 큰 차이가 없어서 그런진 모르겠지만..

blocksize 10~12 정도가 가장 무난한 것 같다.

 

rgb에 대해서 각각 처리 후 합산인데 눈에 보이는건 크게 달라 보이진 않았는데

img = imread("test.jpg");

img = double(img);

blockSize = 32;

[h, w, c] = size(img);

scoreMap = zeros(floor(h/blockSize), floor(w/blockSize));

for by = 1:blockSize:(h - blockSize + 1)
    for bx = 1:blockSize:(w - blockSize + 1)

        totalEnergy = 0;

        % RGB 각각 처리
        for ch = 1:3

            block = img(by:by+blockSize-1, ...
                        bx:bx+blockSize-1, ch);

            block = block - mean(block(:));

            F = fft2(block);

            mag = abs(fftshift(F));

            center = floor(blockSize/2) + 1;

            radius = 4;

            [X, Y] = meshgrid(1:blockSize, 1:blockSize);

            dist = sqrt((X-center).^2 + (Y-center).^2);

            mask = dist > radius;

            energy = sum(mag(mask).^2);

            totalEnergy = totalEnergy + energy;
        end

        score = log(1 + totalEnergy);

        iy = (by-1)/blockSize + 1;
        ix = (bx-1)/blockSize + 1;

        scoreMap(iy, ix) = score;
    end
end

scoreMap = scoreMap - min(scoreMap(:));
scoreMap = scoreMap / max(scoreMap(:));

figure;
imshow(uint8(img));

figure;
imagesc(scoreMap);
axis image;
colorbar;
title("RGB FFT Diversity");

 

blocksize = 10

 

두개를 비교해보면 하단에 파란색 마트 정도는 좀 더 잘 보이는 차이가 있는데

처리 시간대비로는 그리 큰 효과를 얻긴 힘들지도 모르겠다.

grayscale rgb 각각 처리후 합산

 

 LAB 으로 전환해서 색상에 대해서 처리, blocksize = 10

% ==========================================
% LAB 기반 Fabric Diversity / Anomaly Map
% ==========================================

img = imread("test.jpg");

img = double(img) / 255.0;

%------------------------------------------
% RGB -> LAB 변환
%------------------------------------------

% sRGB gamma correction 제거
mask = img <= 0.04045;

img_linear = zeros(size(img));

img_linear(mask) = img(mask) / 12.92;
img_linear(~mask) = ((img(~mask)+0.055)/1.055).^2.4;

% RGB -> XYZ
M = [ ...
    0.4124564 0.3575761 0.1804375;
    0.2126729 0.7151522 0.0721750;
    0.0193339 0.1191920 0.9503041];

[h, w, ~] = size(img);

rgb_reshaped = reshape(img_linear, [], 3);

xyz = rgb_reshaped * M';

X = xyz(:,1) / 0.95047;
Y = xyz(:,2);
Z = xyz(:,3) / 1.08883;

% XYZ -> LAB helper
f = @(t) ((t > 0.008856).*t.^(1/3) + ...
         (t <= 0.008856).*(7.787*t + 16/116));

fx = f(X);
fy = f(Y);
fz = f(Z);

L = 116*fy - 16;
A = 500*(fx - fy);
B = 200*(fy - fz);

lab = zeros(h, w, 3);

lab(:,:,1) = reshape(L, h, w);
lab(:,:,2) = reshape(A, h, w);
lab(:,:,3) = reshape(B, h, w);

%------------------------------------------
% FFT Block Analysis
%------------------------------------------

blockSize = 32;
step = 16;   % overlap sliding

outH = floor((h - blockSize)/step) + 1;
outW = floor((w - blockSize)/step) + 1;

scoreMap = zeros(outH, outW);

% FFT mask 생성
center = floor(blockSize/2) + 1;

[Xm, Ym] = meshgrid(1:blockSize, 1:blockSize);

dist = sqrt((Xm-center).^2 + (Ym-center).^2);

radius = 4;

highMask = dist > radius;

%------------------------------------------
% Sliding Window
%------------------------------------------

oy = 1;

for by = 1:step:(h - blockSize + 1)

    ox = 1;

    for bx = 1:step:(w - blockSize + 1)

        totalEnergy = 0;

        % LAB 각 채널 FFT
        for ch = 1:3

            block = lab( ...
                by:by+blockSize-1, ...
                bx:bx+blockSize-1, ...
                ch);

            % 평균 제거
            block = block - mean(block(:));

            % FFT
            F = fft2(block);

            mag = abs(fftshift(F));

            % 고주파 에너지
            energy = sum(mag(highMask).^2);

            totalEnergy = totalEnergy + energy;
        end

        % 로그 스케일
        score = log(1 + totalEnergy);

        scoreMap(oy, ox) = score;

        ox = ox + 1;
    end

    oy = oy + 1;
end

%------------------------------------------
% Normalize
%------------------------------------------

scoreMap = scoreMap - min(scoreMap(:));
scoreMap = scoreMap / max(scoreMap(:));

%------------------------------------------
% Display
%------------------------------------------

figure;
imagesc(img);
axis image;
title("Original");

figure;
imagesc(scoreMap);
axis image;
colorbar;
title("LAB FFT Diversity / Anomaly Map");

 

block 10, step 16

 

block 10, step 10

 

rgb 각각 처리후 합산 RGB->LAB, 색상만 처리