2016-11-02 262 views
1

爲了在圖像處理中完成我的任務,我使用了R Studio。我正在使用'EBImage,fftw,...'庫。我有一個關於傅立葉分析和功率譜的問題。R - 圖像的FFT傅里葉頻譜

1)我有2D矩陣,這是一個圖像。此圖像由黑色和白色的水平線組成。我想進行傅里葉變換並通過繪圖顯示其大小。由於圖像中的黑線是水平的,所以功率譜將具有垂直線。我找不到實現它的方法。

我現在用的這個圖像進行測試: Square

當它是真實的,這個形象對水平和垂直軸一些週期性的頻率。然後。 FFT頻譜應該看起來像一個「+」。但是我發現的東西有點不同。

這裏是我的代碼:

setwd(".../Project/R/Workspace/Task1") 
library("EBImage" , lib.loc="~/R/win-library/3.2") 
library("fftwtools", lib.loc="~/R/win-library/3.2") 
library("fftw", lib.loc="~/R/win-library/3.2") 

# Image Acquisition 
img <- readImage(".../Project/Beispielbilder/drmcircle.jpg") 
display(img, title='Image') 

# Grayscaled 
img_gray<-channel(img,"gray") 

# FFT 
img_ff <- fft(img_gray) #fftw2d 

magntd <- sqrt(Re(img_ff)^2+Im(img_ff)^2) 
phase <- atan(Im(img_ff)/Re(img_ff)) 

plot(log(magntd),main="FFT") 

結果,這是我所: FFT Spectrum

這裏是我的questoions:

1)我怎樣才能得到一個正確的頻譜圖片? 2)我怎樣才能將它看作是一幅圖像,而不是情節? (請參閱上面的示例鏈接。)

感謝您提前給予的幫助。

+0

你能爲我們提供一個小的可重複的例子(如與您的代碼樣本二維矩陣一起)的http://stackoverflow.com/questions/5963269/ how-to-make-a-r-reproducible-example – Eugen

+0

你好@Eugen,我已經更正了我的問題以獲得清晰的定義。這裏是示例鏈接,我想要達到的。 [http://matlabgeeks.com/tips-tutorials/how-to-do-a-2-d-fourier-transform-in-matlab/] – EMI

+0

我解決了我的問題,在這裏我給代碼。 – EMI

回答

2

我解決了它。我在這裏添加代碼。請看到的是,從拍攝fftshift功能:How to write fftshift and ifftshift in R?

setwd(".../Project/R/Workspace/Task1") 
library("EBImage" , lib.loc="~/R/win-library/3.2") 

# Image 
img <- readImage(".../Project/Beispielbilder/drmtri.jpg") 
display(img, title='Image') 

# Grayscaled 
img_gray<-channel(img,"gray") 

# FFT 
img_ff <- fft(img_gray) #fftw2d 


################################################### 
################################################### 
# FFT SHIFT 
fftshift <- function(img_ff, dim = -1) { 

    rows <- dim(img_ff)[1]  
    cols <- dim(img_ff)[2]  

    swap_up_down <- function(img_ff) { 
    rows_half <- ceiling(rows/2) 
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)])) 
    } 

    swap_left_right <- function(img_ff) { 
    cols_half <- ceiling(cols/2) 
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half])) 
    } 

    if (dim == -1) { 
    img_ff <- swap_up_down(img_ff) 
    return(swap_left_right(img_ff)) 
    } 
    else if (dim == 1) { 
    return(swap_up_down(img_ff)) 
    } 
    else if (dim == 2) { 
    return(swap_left_right(img_ff)) 
    } 
    else { 
    stop("Invalid dimension parameter") 
    } 
} 

ifftshift <- function(img_ff, dim = -1) { 

    rows <- dim(img_ff)[1]  
    cols <- dim(img_ff)[2]  

    swap_up_down <- function(img_ff) { 
    rows_half <- floor(rows/2) 
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)])) 
    } 

    swap_left_right <- function(img_ff) { 
    cols_half <- floor(cols/2) 
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half])) 
    } 

    if (dim == -1) { 
    img_ff <- swap_left_right(img_ff) 
    return(swap_up_down(img_ff)) 
    } 
    else if (dim == 1) { 
    return(swap_up_down(img_ff)) 
    } 
    else if (dim == 2) { 
    return(swap_left_right(img_ff)) 
    } 
    else { 
    stop("Invalid dimension parameter") 
    } 
} 
################################################### 
################################################### 
# FFT SHIFT 


# Magnitude and Phase 
magntd <- sqrt(Re(img_ff)^2+Im(img_ff)^2) 
phase <- atan(Im(img_ff)/Re(img_ff)) 

img_fftsh <- fftshift(magntd) 

display(log(img_fftsh),title="FFT") 
+0

謝謝你這個非常有用的答案。 – SeldomSeenSlim