簡單筆記一下今天工作的成果:將 Matlab 的矩陣噴出成 Amira 的 .am
檔案,也就是 Amira 3D mesh。
這個成果是根據我手上有的 .am
檔案還有 Matlab Central 上的 LoadData_Amira.m
逆推的,特此致謝。
雖然從 LoadData_Amira.m
看起來 AM 檔案有支援 uint8 之外的其他數值格式,但是我剛才弄半天弄不太出來,所以為了簡便,就先強制輸出為 uint8 吧,科科!
剛才突然想到,之前一直有個很差的習慣,就是遇到2D的矩陣,就用很彆扭的方式去解決,例如要解2D矩陣的 max
,就寫
A = rand(6, 3);
max_num = max(max(A));
照理來說min
和mean
這類 會從 column 開始做起然後才做更高維度 的函數,我也都是用一樣的方法來解決。
可是俗話說的好, 夜路走多了會遇到鬼,如果遇到 std
我就沒轍了,如果你這樣寫你就準備出糗了
有在用 Matlab 的人可能有用過 polar
這個畫極座標的繪圖函數,但是很令人崩擴的是 polar
這個函數本身功能非常有限,如果要多別的事情,都要靠 handler 還有 set
搭配來完成,不像 plot
早就已經做到萬能的地步,幾乎想要做什麼,都可以直接在 argument 裡面完成。
現在的問題是,假設我想要把兩個函數都畫在同一張 polar 上面,聰明的各位可能馬上就會想到可以透過 hold on,
來實作,這是沒錯的,但是如果實際去做,會發現一個很糟糕的事情:polar的limit會被第一個畫的函數所限制住,例如假設我想要畫以下的函數:
最簡單的寫法大概是這樣:
close all;
N = 101;
theta = linspace(0, 2*pi, N);
y1 = abs(cos(theta));
y2 = 1.5 * abs(sin(theta));
figure(1),
polar(theta, y1, 'k');
hold on,
polar(theta, y2, 'r');
hold off;
但是畫出來卻會變這樣:
囧!因為 limit 已經被第一個polar決定了,所以如果後面再hold on
上去的圖案更大,就只能悲劇。
Amira 是一個許多做顯微術或者生物影像的人都會用到的軟體,可以用它來具象化許多3D的data。
除了原生的AM檔案之外,有時候我們也會希望藉由 Amira來顯現來自 Matlab我們加工後或者計算後的檔案,因此就會需要由 Amira來打開 .mat
檔案;照理來說 Amira 是支援的,但是我自己實際測試,卻會噴出這樣的錯誤訊息:
Reading minivirus_A.mat ...
opening connection to Matlab engine...
Couldn't open libeng.dll:
File does not exist.
Couldn't open libmat.dll:
File does not exist.
Couldn't open libmx.dll:
File does not exist.
Failed to connect to Matlab
Be sure Matlab is correctly installed, and the path is valid
Error opening file G:/minivirus_A.mat
網路上查半天,都找不到問題所在或者解法,但是自己嘗試去解決,卻發現解法超乎想像地簡單。
首先先試著去找 libeng.dll
、libmat.dll
或libmx.dll
這三個鬼東西在哪。答案很簡單,就在
C:\Program Files\MATLAB\R2012b\bin\win64
的下面。我想 Amira 抓不到的原因非常簡單,是因為 Amira 不知道 64-bit 的 Matlab 已經改變存放這些lib的位置了。
解決方法
解法非常簡單,就是改系統變數中的 PATH
值。
首先對著 我的電腦 按右鍵,選 內容 ,接著選 進階系統設定 ,然後按 進階 這個頁籤,最後按最下面的按鈕 環境變數 。接著在下面的小框框中找到PATH
,點選 編輯,然後在最後面加上
;C:\Program Files\MATLAB\R2012b\bin\win64
就大功告成,此時重開 Amira 後就可以開啟 .mat
檔案了。
之前就有想過要用解析解求一個均勻的圓球的質量投影,但是推半天沒推出來,只好先用MATLAB做出一個3D的volume,再用 sum(V,3)
來求投影質量。我知道這樣超白癡的,而且要是球稍微大一點,電腦就記憶體不足當掉給你看了。
剛才突然... 就夢到正確的解析解做法了......
原理
已知圓球球面的公式是\
其中 為一常數,是球面的半徑。
現在我們將其推廣,把球面推廣成球體,則參數 要滿足
其中
所以要快速求出一個球體的質量投影,只要移項還有乘以2就好了,得到:
不過記得要把超出 support 的地方手動設為零,不然可是會出現 complex 的。請直接參考下面的實作程式碼。
MATLAB 實作
function Z = sphere_proj(DIA)
RAD = DIA/2; % radius
x = linspace(-RAD, RAD, DIA);
[X, Y] = meshgrid(x, x);
R = sqrt( X.^2 + Y.^2);
Z = 2 * sqrt( RAD^2 - R.^2 );
Z( R > RAD) = 0;
簡介
訓練目標
以1D為切入點,了解 EEMD 之基本原理,並以 Matlab 語言實作之。
Matlab
完成這份訓練講義之後,預期學習者將會應用到或者學到如何
- 創造並使用
function
- 找出資料向量中的極值
- 使用
interp1
內插出平滑的資料曲線 - 打造 EMD 引擎
- 將 white noise 加入資料之中
- 完成 EEMD 函數
習題
尋找極值函數 findExtrema
在這題當中,我們要學會寫出一個傳入資料向量 data
,可以回傳資料中極值的 值 以及 位置 的函數。以下將分為四個小題進行。
學會如何自定函數
試找尋資料,學習如何創造一個自定函數,或查詢 Matlab 內建教學
doc function
或參考以下範例:
function y = numAdd(a, b)
y = a + b;
尋找最大值
試寫一個函數 findMax.m
,輸入任意長度向量 data
可傳回兩組輸出:
- 最大值所在的位置
max_ind
- 最大值的值
max_val
提示: 函數的第一行可能長成這個樣子:
function [max_ind, max_val] = findMax(data)
% coding here
請注意,一串資料可以有多個 local maxima。此題要求的是找出所有的 local maxima,而非僅是資料中的最大值,否則用內建函數 max()
即可達成。
尋找最小值
承上題,找出「所有」最小值 min_val
及其分別之所在位置 min_ind
。
提示: 試想,是不是將上一題的傳入資料由 data
改為 -data
就可以達成目標了?
尋找所有極值
綜合前面三題所學,請完成一個可以找出所有極值,並分別回傳最大最小值的值與位置的函數 findExtrema.m
function [max_ind, max_val, min_ind, min_val] = findExtrema(data)
% coding here
注意: data(1)
與 data(end)
是否也須考慮為「極值」?如果是,其原因為何?
內插曲線
下面這兩題,我們將會利用前一大題所得到的資料,以及 Matlab 內建的 interp1
函數,得到 EMD 所需要的「趨勢」曲線。
使用離散資料內插曲線
利用 interp1
,將 max_ind
和 max_val
內插成一條長度和取樣率皆與 data
相同之「連續」曲線 max_trend
。
利用 interp1
,將 min_ind
和 min_val
內插成一條長度和取樣率皆與 data
相同之「連續」曲線 min_trend
。
注意: max_trend
或min_trend
的頭尾,都要和 end
一樣(可以想成是 boundary conditions),否則會有些區域沒有定義。
求得「平均」曲線
將上小題的兩條曲線做點對點平均,得到 avg_trend
。
打造 EMD 引擎
參考 paper 上的演算法,以及上面所練習的技巧與成品,應該已經可以完成 EMD 的引擎。其界面如下:
function modes = emd(data, nmode)
% coding here
其中 nmode
是所要求的模態數量。請注意,如果原本的 data
長度是 的話,則得到的輸出資料 modes
的大小會是 。
打造 EEMD 引擎
如果前面的 emd.m
已經驗證完成,那麼接下來就是讓它進化成 EEMD 的時候了。比起 EMD , EEMD 多了下面兩項特色
- 外加了額外的 white noise 進去
- 先跑出 N 組各別使用不同的 white noise 所解出來的結果,再一起取平均。
下面我們只拆成兩小題練習,不足之處請直接參閱 paper。
加入 white noise
假定使用者指定一個特定的 white noise 強度 wn_std
,我們有兩種方式可以將等效強度的 white noise Wn
加到 data
裡面去,第一種是:
其中 指 data
, 是求標準差的函數;或者第二種方法
試問:兩種方法何種較優?較劣者,有什麼致命的缺點?
請利用 Matlab 內建的 randn()
完成此題。
完成 EEMD 架構
參考 paper 上的演算法,以及上面所練習的技巧與成品,終於我們可以完成 EEMD 的全面架構。完成品的輸入資料除了 data
以及nmode
之外,還多了指定 ensemble 數量的 nens
與 white noise 的標準化強度 wn_std
。其界面如下:
function modes = eemd(data, nmode, nens, wn_std)
% coding here