顯示具有 ImageJ 標籤的文章。 顯示所有文章
顯示具有 ImageJ 標籤的文章。 顯示所有文章

2021年1月11日

用imagej高斯模糊萃取輪廓

最近被諮詢了兩個imagej的分析案例,雖然影像和目的需求都不一樣,但分析的手法卻是雷同的,都牽涉到高斯模糊。

什麼時候會用到高斯模糊呢?「瞇著眼睛看,好像就會看到答案」的影像,應該都會用到高斯模糊。這個對有近視但沒戴眼鏡的人最有感覺了,當你看到一個人,但是你並不想注意她的眼睛鼻子長什麼樣子,你只想關注他的外型輪廓、體型、身高等特徵,那麼你把眼鏡拔下就有高斯模糊的效果了。


案例一,左圖是一個SEM的圖片,提問者想要分析平滑區域和有孔洞的區域兩者面積的差異。這個案例需要先把兩個區域切割清楚,但如果直接用灰階值做閾值的劃分,一定會有誤。

仔細看看這張圖,如果你瞇著眼睛看,不要把細節看得那麼仔細,就會變成右圖。你就會發現那就是分析的目標了。

高斯

方法是先將影像做高斯模糊,可以消除細節,呈現輪廓

 Process › Filters › Gaussian Blur...

再做中值濾波

Process › Filters › Median...



第二個案例是研究珊瑚的年層線,左圖可以看見珊瑚切片上一圈一圈的年層線,如同上例,也是瞇著眼睛消除細節後,就比較能明顯區分出來。

Montage

方法一樣是先做高斯模糊,然後再做尋找邊緣等分析程序。

run("Gaussian Blur...", "sigma=5");
run("Find Edges");
setMinAndMax(0, 10);
setOption("BlackBackground", false);
run("Convert to Mask");
run("Dilate");
run("Erode");
run("Erode");
run("Despeckle");  

2020年6月4日

Imagej入門課程影片

上週四應生物學科中心之邀,講了一堂Imagej的入門課程。

目的希望用短短一小時內的時間,至少知道以下問題如何分析
1.測量圖像的長度、面積等資訊
2.測量圖像內的顏色資訊
3.手動或自動計算顆粒數量

為此,我難得作了入門講義給與會者,因為我覺得我一步一步講完之後,初學者應該會不知道剛剛作了哪些步驟。

講義的電子檔(ppt、pdf、png)在這裡



imagej入門講義



演講完之後,想說講義都作了,不然再錄個教學影片吧!





片長約38分鐘,不是一鏡到底,是分成了幾段個別錄製,最後再剪接成一部影片的。以往我會分開成好幾部影片,不過Youtube最近更改了上傳介面,每上傳一部影片就要個別設定細節,所以我乾脆剪接成一部吧。

剪接的時候,本來是用Openshot,不過剪接轉檔的時間太久了,後來改用 Avidemux ,在不更改影音編碼的前提下,剪接輸出的速度只要一兩秒就完成了,對我真是相當適用。

2020年5月18日

用imagej分析水蚤心跳-kymograph

之前有做過 用Tracker分析水蚤心跳,但其實imagej 就可以做了,忘記當時為什麼只用 Tracker來做?

首先要把分析的影片匯入到imagej當中,這要先安裝ffmpeg才可以。這在之前的imagej教學文章已經提過好多次了。

接下來是分析的方式,基本上就是看圖片的明暗變化,來判斷心臟跳動與否。使用到的功能是stacks裡的plot z-axis profile

Plot of Dynamic Reslice of 水蚤心跳 

 另外還可用 kymograph的方法來視覺化心臟跳動的過程,先在圖像上標示分析的線段後,用image/stacks/reslice,它就會把stacks的每一個frame都切下線段的部份重組成一張圖,觀察圖上的明暗變化,也可以看出心臟跳動的規律性。
  Reslice of 水蚤心跳


橫軸的部份是像素變化,縱軸則是時間(影格)

如果需要即時改變選擇區域,也可以改用 image/stacks/dynamic reslice

這種reslice做出的 kymograph,也可以同時觀察到不同位置的時間變化。

使用說明請見影片

2019年12月30日

用imagej做曲率分析的原理與實作

因為有朋友的研究需進行物體的曲率測量,因此周末研究了一下到底怎麼做。
(急著想知道怎麼做的,直接到頁面最底下看。前面一大段都是原理和原理實作)

要描述一個曲線有多彎,會用到一個密切圓的概念,曲線上取三個點,三點可成一個圓(三點都在圓弧上)。三點非常靠近時形成的圓,就是曲線上該點的密切圓。

當曲線越「彎」,密切圓就越小,而曲線越「直」,密切圓就越大。一個超大的圓,大到極大,圓弧看起來就變成直線。就像地球超大,你看海面就是幾乎是水平的。


以下兩張圖,分別可以看出在這個函數的曲線上,不同點的密切圓大小不同。







那麼曲率要怎麼描述呢?可以用1/r來表示,r是密切圓的半徑。下圖的黑線處就是這個sin函數在各點的曲率變化。



上述的過程,可以看這部影片的演示




其實要分析物體的曲率變化,其實就是要用函數圖形去貼合要測量的部分,然後去分析函數的曲率變化。接下來就是函數怎麼找出來。

方法是弄出一些控制點,用這些控制點擬合出一個多項式函數。下面的影片演示的是用ggb的FitPloy來擬合數個點來產生函數,然後再對函數求出一階微分和二階微分,再用那些微分方程式算出曲率。看起來好難,但其實意外的簡單。

這個連結裡有求曲率的公式,也詳細說明了公式怎麼出來的。
如何通俗地理解曲率?


我也繼續以ggb為工具,實際計算出雞蛋的表面曲率



用數個點去找出曲線,不只一種方法。可以用貝茲曲線或是B-spine(B樣條曲線)來做。
貝茲曲線在向量繪圖軟體,如inkscape會用到。我自己以前只是會用,但是實際原理不甚清楚,這次趁著機會也下手實作了解一下。

以下影片是用ggb來實作貝茲曲線


本來也想實作B-spine(B樣條曲線),但是發現有一點複雜,可能以後用python來實作看看吧。樣條曲線的原理說明,這篇應該是最清楚的
簡單粗暴 B樣條曲線入門


寫了一大堆之後,終於要進入用imagej測量曲率的部分了。
這裡需要使用的plugin叫做Kappa,已經有人寫得非常清楚了
ImageJ实用技巧——曲率计算与拟合(插件篇)
kappa 的github https://github.com/brouhardlab/Kappa/tree/master/docs
這裡就放一段實作的教學影片吧





2019年12月22日

用imagej做熱像影像的衍射映射變形合成

直接先來看結果吧!
以下這兩張圖,一個是手機相機拍攝的影像,另一張是熱像儀(Seek thermal)拍攝的影像。兩張影像來自不同相機、不同角度、不同焦距。影像大小也不一樣。




熱像影像有物體的溫度分佈,而相機影像有物體的紋理和邊緣,如果想要同時擁有兩者的優點,就像下圖的影像,該怎麼做呢?




我們使用imagej就可以做到了。

以下是影片教學。我也同時以文字紀錄在下



首先將兩張待處理的影像用imagej開啟,在imagej上用multi-point tool在定位點上做標記。我是在拍攝前就先以錢幣定位,分析時就是用錢幣位置來定位。做影像定位的時候,要注意每個點的順序是有意義的,兩張圖的定位點要能互相對應。

定位好了,就選擇Plugin/Transform/Landmark Correspondences,注意Transformation class要選擇affine。

如此就可以將圖片變形一致了,如下圖


下一步,由於希望將兩圖合成時,由手機影像提供物體的邊緣和紋理,所以選擇手機影像檢測邊緣。方法是Process/Find edge
圖片就會成為這樣



接下來關閉其他圖片,只留下變形後的熱像影像和邊緣化的手機影像,然後將兩張圖片疊成一張stack。
方法  Image/stack/Images to stack
選項就用預設值就好

接著就可以做疊合
方法是選擇這個stack之後,選擇Image/stack/Z project
Projection type選擇 Sum slices
這樣就可以做出精確疊合的影像了




2019年11月23日

用Imagej的Macro製作花磚圖片

上個月的時候,看到一則新聞《這些花磚圖案超美 竟全靠飛蛾貢獻》-暨南大學陳彥錚教授將蛾類圖片利用程式製作成花磚。我看到的時候,覺得這好有創意,很想也跟著做看看。本來想用python來寫,不過想到還有GUI要處理,乾脆試著用Imagej加上Macro來處理看看好了。


程式放在GitHub
主要就是兩個程式
  • makePolygon.ijm
  • afterPolygon.ijm

以下處理的範例-銀目天蠶蛾的圖片來自飛蛾資訊分享站(http://twmoth.tesri.gov.tw/peo/FBMothInfo/23805)

先來分析看看怎麼處理,這是最後的成果,捨去背景不看,其實內部就是同一張圖片旋轉不同角度拼接起來的。方塊的方向不同(蛾頭方向朝內或朝外),做出來的效果也不同



再看原圖一下,也就是說要在原圖的適當位置上找出一個合適角度的正方形

 

但是要用imagej畫這種正方形很難。第一、它的邊長不是平行x軸和y軸,第二、角度又要符合自己需求,第三、大小又要可調整。

後來我想到一個方法,既然不能直接畫出正方形,那我就指定兩點,再用這兩點來畫出正方形吧。哪兩點呢?就用直線工具畫出中心點到正方形的其中一個頂點,剩下再用程式去找出來另外三個頂點。

其實原理在高中數學就有,把頂點一的座標減去中心點的座標,就可以得到一個中心點到頂點的向量,接下來再用旋轉矩陣的公式,找到旋轉90、180、270度之後的新向量,再把這些向量加回原來的中心點座標,那就可以得到另外三個頂點的座標了。


這就旋轉矩陣




但是找出四個座標之後,還不能直接畫出方形,因為四個點的座標可能會超出畫面,所以必須調整畫布大小。所以下一步就是分別計算四個點的x和y是不是有超過圖片的寬和高,如果有,就調整畫布大小。最後才可以把選取區畫出來了。


以上這些步驟就是在makePolygon.ijm這個程式所執行的

第二個程式makePolygon.ijm就只是把找到的正方形貼到一張新的大圖上,拼出一個花磚
看起來很簡單,但在imagej上要實現,花了一些功夫(我想應該有更快的作法,但是目前還不會就是了)

以下第一部影片介紹怎麼操作



第二部影片就是說明程式製作的原理是什麼,但是第二支程式的細節太多,多到我忘了,反正就是不斷的位移和旋轉....

2019年9月22日

用imagej的macro製作的無腦光譜分析程式


想要產生出如下圖這樣各種光譜分析圖,或是利用CFL燈泡(省電燈泡、日光燈等)來校正一整個資料夾的光譜照片,然後也是產生像下面這樣的圖。你覺得需要多少步驟呢?
只要兩個按鍵就好了喔,連輸入光譜值找對應的X座標都不用喔。





如何完成光譜分析呢?這部影片告訴你如何按兩個鈕就完成


程式碼

為什麼可以按兩個鈕就完成,連輸入波長都不用呢?因為我已經在程式中去做了那些工作。讓程式自動用hue值辨識光譜顏色、尋找波峰,再用預設的光譜值套入來進行校正。這樣的程式到底是怎麼製作呢,我也將程式說明錄了下來(其實這是避免自己將來萬一忘了)







如果是想要用imagej自己手動作這些工作,我在幾年前的文章曾經寫過,可以參考一下
用imagej作分光光度計與光譜儀的量測分析

2019年8月4日

用Imagej製作灰階影像加彩色線條的視覺效果


前幾天看到網路上流傳了一個在灰階影像上加上彩色線條之後,看起來就好像是彩色圖片的視覺效果,我看了覺得很有趣,很想自己來動手玩一玩。於是就用imagej玩看看,順便寫了一個code script放在文章底下,有興趣的可以複製貼上玩玩看。

可以做什麼效果呢?可以針對個別區域製作橫紋效果。下面組圖是用同樣粗細的線寬,但是不同的線距


也可以做出斜紋效果



或是直紋效果



這些線寬都不一樣,但是線寬和線距比例都是1:3



同樣的code改一改,也可以做出網紋的效果



以下就是教學影片了,前兩部是從原理製作開始講,如果要直接使用code來做的人,就從第三部來看就可以了。





以下就是imagej的macro囉
===========================
lineW = 2;
space = 8;
angle = 0;


run("Select None");
rename("color");
run("Duplicate...", "title=grey");

selectWindow("grey");
run("8-bit");
run("RGB Color");
run("Set Measurements...", "area bounding redirect=None decimal=3");
roiManager("Select", 0);
run("Rotate...", "angle="+ angle +"");
roiManager("Add");
run("To Bounding Box");
run("Measure");

BX = getResult("BX", nResults-1);
BY = getResult("BY", nResults-1);
Width = getResult("Width", nResults-1);
Height = getResult("Height", nResults-1);

lineN = floor(Height/(lineW+space));
print(lineN);
for (i = 0; i < lineN; i++) {
//建立選取區
run("Specify...", "width="+ Width +" height="+ lineW +" x="+ BX +" y="+ BY+(space+lineW)*i +"");
//run("Rotate...", "angle="+ angle +"");
roiManager("Add");
//交集運算
roiManager("Select", newArray(1,roiManager("count")-1));
roiManager("AND");
//產生要繪圖的選取區ROI
roiManager("Add");
roiManager("deselect")
roiManager("Select", roiManager("count")-2);
roiManager("Delete");

//roiManager("Select", roiManager("count")-1);
//roiManager("Delete");
}

a1 = newArray(roiManager("count")-2);
for (i=0; i<a1.length; i++)
  a1[i] = i+2;
roiManager("Select", a1);
roiManager("Combine");
roiManager("Add");

a1 = newArray(roiManager("count")-2);
for (i=0; i<a1.length; i++)
  a1[i] = i+1;
roiManager("Select", a1);
roiManager("Delete");


roiManager("Select",roiManager("count")-1);
run("Rotate...", "angle="+ -angle +"");
roiManager("Add");

//====複製原圖,貼到黑白圖====
selectWindow("color");
roiManager("Select", roiManager("count")-1);
run("Copy");

selectWindow("grey");
roiManager("Select", roiManager("count")-1);
run("Paste");
//=========================

roiManager("Select", newArray(roiManager("count")-1,roiManager("count")-2));
roiManager("Delete");

selectWindow("color");
run("Select None");

selectWindow("grey");
run("Select None");
rename(lineW+":"+space);

text = "w"+lineW+":space"+space;
setFont("SansSerif", 14, " antialiased");
makeText(text, 20, getHeight()-20);
run("Add Selection...", "stroke=yellow fill=#660000ff new");
run("Select None");

2019年8月1日

用Imagej做顯微追焦

顯微追焦是什麼?其實我也沒想過要怎麼去命名這種影片。我用以下的影片來說明好了,它就是把在顯微影片下,本來四處跑來跑去的生物,用影像處理讓它固定在同一個地方,甚至身體方向也固定下來。這跟攝影的追焦有點像啦,所以就叫顯微追焦囉。
處理成這樣的影片有什麼好處呢?我實作的結果,發現這樣更能夠更仔細看到這些生物的細節,或是運動方式。如果是原先到處跑的影片,你又要追它的位置,又要看清楚,那就實在太累了。



處理的流程是先標記生物的位置,然後用程式碼位移每個影格,讓標記的位置都放在同一個地方,然後你就可以裁切獲得追焦後的影像。
這系列的教學影片共有四集,會用到兩種簡單的程式碼,分別是「用點標記」和「用線標記」的兩種。程式我就貼在文章最底下
要使用這些程式,就只要在Imagej中開啟File/New/Text Window
把以下的程式按照需求貼在那window就好了(要確認功能表上的Language,看一下勾選IJM喔)。











=====用點標記的======
origin_x= getWidth/2;
origin_y= getHeight/2;
for (i = 0; i <  nResults ; i++) {
slice_x = getResult("X", i);
slice_y = getResult("Y", i);

translate_x= origin_x-slice_x;
translate_y= origin_y-slice_y;

setSlice(getResult("Slice", i));
run("Translate...", "x="+translate_x+" y="+translate_y+" interpolation=None slice");
}

====用線標記的=========
origin_x= getWidth/2;
origin_y= getHeight/2;
for (i = 0; i <  nResults ; i++) {
BX =  getResult("BX", i);
BY =  getResult("BY", i);
W =  getResult("Width", i);
H =  getResult("Height", i);
A = getResult("Angle", i);

if(A>-180 && A<=-90){
X=BX+W;
Y=BY;
}
if(A>-90 && A<=0){
X=BX;
Y=BY;
}
if(A>0 && A<=90){
X=BX;
Y=BY+H;
}
if(A>90 && A<=180){
X=BX+W;
Y=BY+H;
}


translate_x= origin_x-X;
translate_y= origin_y-Y;
rotate_A = A + 90;
setSlice(getResult("Slice", i));
run("Translate...", "x="+translate_x+" y="+translate_y+" interpolation=None slice");
run("Rotate... ", "angle="+ rotate_A +" grid=1 interpolation=Bilinear");
}
==============

2019年7月8日

imagej+colab機器學習作細胞分割

這是延續前篇的續作。

前一篇用了imagej的macro做了細胞的分割,不過實際上我認為誤差還是蠻大的,心裡想著應該還有更好的方法。想著想著突然想到這個問題,其實用機器學習的監督式學習可以解決。

只要把每個細胞團的參數拿出來餵給機器去學習,告訴機器「這種細胞團參數叫做一個細胞,那個細胞團參數叫做兩個細胞....」,當機器學會之後,就可以讓機器去預測那些細胞團是由幾個細胞組成的就可以。不過要做這樣的人工智慧,前提得先有工人智慧,我需要先人工辨識標記好每個細胞團到底是幾個細胞組成。

標記好或是機器學習預測好之後,我還需要有檢核的程式,讓我可以肉眼比對到底分類的結果如何。

整個機器學習的檔案架構是這樣的

├── dataSet_predict   用來預測的文件
│   ├── Cells
│   ├── Results
│   └── Roi
├── dataSet_train   用來訓練的文件
│   ├── Cells
│   ├── Results
│   └── Roi
├── image     圖像的資料夾
│   ├── predict
│   └── train
├── ImageCombine
└── macro
    ├── cellSegment.ijm
    ├── checkModel.ijm
    ├── makeModel.ijm
    └── MakePredic.ijm


看起來很複雜,不過使用者只要記得把細胞的圖像分成兩部份,一部分放在image/predict,一部分放在image/train。當然train量一定是比較少的啊,我們是要用train裡的圖像來預測predict裡頭的圖像啊。

程式的部份我做了四支程式,三支是在imagej下執行,一支是python,讓它在colab來執行。先介紹三支imagej的程式
makeModel.ijm  執行之後選擇圖片所在的資料夾,自動會產生相關的文件,放在train的資料夾裡
checkModel.ijm  用來檢查那些人工標記或是電腦辨識之後的分析結果是否正確
MakePredic.ijm 用來產生預測用的檔案

決定把機器學習的檔案放在colab上執行,是想說這樣使用者就省去自己部屬python環境的時間,反正只要上傳訓練檔案,等個十幾秒後自動會下載分析後的檔案



整個完整的分析、人工標記、上傳的流程,我照慣例錄製成影片說明。有需要的可以看看



整個程式檔和範例都放在雲端
https://drive.google.com/drive/folders/1m4MxOdCAziTr3tdNQhKP8H6y6nqfCNoz?usp=sharing


2019年6月27日

以imagej的macro實作GUI,進行重疊細胞的分割

上週五收到一個生技公司人員的來信,信中問到要作細胞的分割來計數,但是作分水嶺分割時出現問題。我寫imagej作細胞分割的文章其實寫了好幾篇,本來想說請他看看之前文章就好。不過因為給的圖片很漂亮,剛好議題又是我感興趣的,所以就試試看來解決了。

來看一下圖片。這是利用不同劑量的藥劑處理非癌化細胞後,再以螢光染劑進行染色。研究者的目的是要計算這有多少細胞。


由於有些細胞是重疊的,一般想法就是直接作分水嶺切割,然後再作analyze particles。不過這樣的作法會有大問題。


這是我們期待看到的,兩個細胞被分水嶺算法剛好切開。



但實際上會發生的是,單獨一個細胞往往會被切成好幾塊



或是兩個細胞被切成三塊



如何解決這樣的問題呢?
解決流程是這樣的
(1)先把單獨未重疊的單顆細胞抓出,並計數。
(2)只針對重疊的細胞作分水嶺切割,然後作計數。

現在問題來了,怎麼知道誰是單顆細胞,這裡會運用一個形態學的運算方式來幫助判斷。
以下兩張圖是先將圖片做threshold處理後,轉二值化影像,再用wand tool選取後,然後用Edit/Selection/Convex hull做出的新選取區。

convex hull是個重點!它就像是用條橡皮筋套住你的選取區,產生的新選取區就是convex hull。你看第一張圖片,是兩個細胞重疊,用convex圈起之後,會留下很多沒圈到的地方。第二張圖片是單獨的細胞,圈起來之後留下的空白很少。

所以啊,我們可以用convex hull和細胞的面積做運算,看看細胞在hull裡頭佔面積的多寡來推論那是不是單獨的細胞。








要計算這個東西,其實不用自己手動一個一個圈細胞。在measure中的測量值solidity就是這個啦。你可以勾選Analyze/Set Measurements/Shape descriptors,就可以呈現出這個數值了。





接下來的分析和製作Macro,我就交給影片說明了


影片中所使用的Macro 就是以下的文字啦

=============================================
Dialog.create("Preference");
Dialog.addNumber("solidity:", 0.88);
Dialog.addNumber("cell area min :", 2000);
Dialog.addCheckbox("Single cell fill color", true);
Dialog.addCheckbox("Watershed cell fill color", true);
Dialog.show();
solidityValue = Dialog.getNumber();
cellAreaMin = Dialog.getNumber();
singleCheck = Dialog.getCheckbox();
watershedCheck = Dialog.getCheckbox();

//關閉所有圖檔
run("Close All");
//製作結果表
Table.create("cellCounts");

//打開資料夾,取得檔案list
dir = getDirectory("Choose a Directory ");
list = getFileList(dir);
filecount = 1;
listFiles(dir);

function listFiles(dir) {

for (i=0; i<list.length; i++) {
if (endsWith(list[i], "/")){}

// listFiles(""+dir+list[i]);
else{
//print((filecount++) + ": " + dir + list[i]);
open(dir + list[i] );
bname=File.nameWithoutExtension;
cellCounts();
}
}
}

function cellCounts(){
//細胞有幾個
cellNum_single=0;     //一個單獨的細胞
cellNum_overlap=0;     //重疊的細胞團
  cellNum_watershed=0;  //分水嶺切割後,單獨的細胞

rename("original");

//製作空白的檔案
run("Duplicate...", "title=cell_single");
run("Select All");
setBackgroundColor(255, 255, 255);
run("Clear", "slice");
run("Duplicate...", "title=cell_overlap");
run("Duplicate...", "title=cell_watershed");

//分析前處理
selectWindow("original");
run("Duplicate...", "title=a");
run("8-bit");
setAutoThreshold("Default dark");
run("Convert to Mask");

//把所有細胞都抓出,設定measure參數
run("Set Measurements...", "area perimeter shape redirect=None decimal=3");
run("Analyze Particles...", "size="+cellAreaMin+"-Infinity circularity=0.10-1.00 display exclude clear include add");

//用Solidity推斷細胞是否有重疊
numROIs = roiManager("count");
for(i=0; i<numROIs;i++) {// loop through ROIs
solidity = getResult("Solidity", i);

if(solidity>solidityValue){//Solidity大於0.88者,應該是單顆細胞
cellNum_single =cellNum_single +1;
//塗色
if (singleCheck){
selectWindow("cell_single");
roiManager("Select", i);
setForegroundColor(255*random, 255*random,255*random); //隨機填入顏色
run("Fill");
}
}

else {//如果是小於0.88,可能是多個細胞重疊,進入迴圈作分水嶺切割
cellNum_overlap =cellNum_overlap +1;
selectWindow("cell_overlap");
roiManager("Select", i);
setForegroundColor(0,0,0); //填入黑色
run("Fill");

}
}

//到cell_overlap的圖片中,執行分水嶺算法
selectWindow("cell_overlap");
run("Make Binary");
run("Watershed");
run("Select All");
run("Analyze Particles...", "size="+cellAreaMin+"-Infinity circularity=0.10-1.00 display exclude clear include add");
numROIs = roiManager("count");
for(i=0; i<numROIs;i++) {// loop through ROIs
cellNum_watershed =cellNum_watershed +1;
if(watershedCheck){
selectWindow("cell_watershed");
roiManager("Select", i);
setForegroundColor(255*random, 255*random,255*random); //隨機填入顏色
run("Fill");
}
}

//輸出細胞的計數結果
rowsize = Table.size("cellCounts");
Table.set("filename", rowsize, bname,"cellCounts");
Table.set("cell_single",rowsize, cellNum_single,"cellCounts");
Table.set("cell_overlap",rowsize, cellNum_overlap,"cellCounts");
Table.set("cell_watershed",rowsize, cellNum_watershed,"cellCounts");
Table.update("cellCounts");

//關閉圖檔
selectWindow("cell_single");
saveAs("Jpeg", dir+"1/" + bname+"_1_single"+".Jpg");
selectWindow("cell_overlap");
saveAs("Jpeg", dir+"1/" + bname+"_2_overlap"+".Jpg");
selectWindow("cell_watershed");
saveAs("Jpeg", dir+"1/" + bname+"_3_watershed"+".Jpg");


close("cell_single");
close("cell_overlap");
close("cell_watershed");
close("a");
close("original");
close("Results");
close("ROI Manager");
}

2019年2月6日

使用imagej Macro客製化imagej的LUT(lookup table)

這兩天在我的imagej教學影片《05 imagej的LUT》上,有人留言詢問了問題:

你好,我想请问一下如何只是显示比如说25到150之间的lut,而小于25和大于150点不显示,而且在添加calibration bar的时候也是显示25-150呢?谢谢啦!
我现在的问题是我把8bit的图片的lut设置成了jet,也就是最低的灰度值为蓝色,最高的255为红色。但是呢我的图片最高的灰度值为150,最低的为25,我想在这个范围内设置为jet的lut可以吗?也就是25的灰度值为蓝色,150的灰度值为红色。且中间的颜色变化是像jet中连续的。我的问题有点儿长,不知道是不是说明白了。谢谢您!


lut是用來視覺化圖片的工具,作用就是把8bit灰階圖片著色,lut提供的就是著色的規則表。例如。例如下圖左下角的LUT規定了「72這個灰階值,要用R0 G159 B255的顏色來著色」


不同的顏色規則就會形成不同的規則表,在imagej裡頭預設就有數十種LUT可以使用。利用Image/Color/Display LUTs 可以顯示這些LUT


其實我並不太在youtube回答留言者的問題,因為光是我Gmail裡的問題我就回答不完了,再者是youtube的留言提問,並沒有辦法附圖去說明。不過因為留言者的問題剛好跟我過去想解決的問題有關,他寫出的問題有我需要的關鍵字,所以我就投入去找資料了。

留言回應者提到的LUT名稱叫做「jet」,我翻了好久倒是沒看到。後來去查才知道是Matlab這個程式裡頭的LUT,它的顏色配置就是本文第一張圖左上角的樣子,就像是彩虹的顏色。


要解決留言回應者的問題前,得先找到方法去做出這個規則。一開始是在一個論壇上看到有imagej Macro的解法,幾行就可以生成jet lut

r=newArray(256);g=newArray(256);b=newArray(256); 
for (i=0;i<256;i++) { 
    i4=4*i/256; 
    r[i]=255*minOf(maxOf(minOf(i4-1.5,-i4+4.5),0),1); 
    g[i]=255*minOf(maxOf(minOf(i4-0.5,-i4+3.5),0),1); 
    b[i]=255*minOf(maxOf(minOf(i4+0.5,-i4+2.5),0),1); 

setLut(r,g,b); 


雖然可以生出jet Colormap,但是中間的i4-1.5或是-i4+4.5,弄不懂怎麼修改,所以暫時放下,繼續找其他文章。後來在stackoverflow的討論串裡找到解法。


這樣的解法就容易修改了,我只要縮窄原有的範圍(從256縮到151),再做位移(從0往後位移到25)就可以了。根據這位提問者的需求,做出來的就會是這樣。



在imagej點選File/New/Text Window,把下面這段貼到視窗裡,然後把Text Window的Language設定為IJ1 Macro,再按下Run/Run就可以修改Lut了。

在前幾行的xrange、xoffset和NanFill的數值都可以依照需求修改。
==============================================

run("8-bit");

xrange = 126; //灰階值範圍,不需特製的話,使用256
xoffset = 25; //最低灰階值,不需特製的話使用0
NaNFill = 255; //沒有灰階的部份填充的顏色,白色為255,黑色為0

r=newArray(256);
g=newArray(256);
b=newArray(256);

for (i=0;i<xrange;i++) {
x= -1 + i/(xrange/2);
r0 = 1.5 - abs(2*x -1);
    g0 = 1.5 - abs(2*x);
    b0 = 1.5 - abs(2*x +1);
   
    //clamp
    r[i+xoffset]=255*minOf(maxOf(r0,0),1);
    g[i+xoffset]=255*minOf(maxOf(g0,0),1);
    b[i+xoffset]=255*minOf(maxOf(b0,0),1);
}

//填充灰階值下界
for (i=0;i<xoffset;i++) {
r[i]=NaNFill;
    g[i]=NaNFill;
    b[i]=NaNFill;
}
//填充灰階值上界
for (i=xrange+xoffset;i<256;i++) {
r[i]=NaNFill;
    g[i]=NaNFill;
    b[i]=NaNFill;
}

setLut(r,g,b);  //設定LUT
run("Show LUT");  //顯示LUT的list



2019年2月5日

使用Imagej做影像分割的流程-使用PlantCV的流程教學

農業研究上,常常需要分析大量植物的表型,如果用傳統的個別測量方式,費時不少。因此目前已經有開發了自動機器人,固定為大量植栽攝影,而攝影所得的影像資料也是透過程式來加速處理影像分割和分析。這次介紹的網站裡的程式-PLANTCV就是一個用在此目的的PYTHON程式。
https://plantcv.readthedocs.io/en/latest/


我閱讀了其中流程的教學文件,覺得很實用,雖然是PYTHON的程式,不過影像分割的流程其實也能應用在imagej。因此本文就是將其Python流程轉換成Imagej流程的紀錄,如果是用imagej來處理大量影像資料,例如幾百張圖,那麼也只要把這些步驟錄製在macro裡再做修改就可以大量處理了。

以下是一些步驟的個別說明:

轉換色彩空間
以往我使用imagej做影像分割時,都是先用原始影像轉8 bits灰階,再用threshold或是原始彩色影像用color threshold來取得目標物。不過其實可以先將原始影像轉成HSV或是Lab的影像。因為有些目標物在RGB的channel之中未必會明顯突出,但是在HSV或是LAB影像中反而會變得明顯,因此轉成這類型影像後,再針對其中比較明顯的channel來做影像分割。
imagej的指令
Image/Type/RGB Stack
Image/Type/HSB Stack
Image/Type/Lab Stack


針對不同部位使用不同的影像空間的channel
由於植物並非色彩均質的物體,有些地方比較綠、有些地方黃甚至紅,因此可以使用不同色彩空間的channel來分割不同部位,最後再把這些組合在一起就可以。例如在Lab的色彩空間中,a channel是偏紅色綠色的,而b channel是偏藍色黃色的,所以就可以從這些部位個別取出要的部位再組合。


設定閾值做二值化影像
這部份就是直接設定Thershold,選擇要的部位來做二值化影像
imagej的指令
Image/Adjust/Threshold
Image/Adjust/Color Threshold




去除雜點、去噪
在二值化的影像中,難免會有一些噪點,可以使用中值濾波器,或是其他濾波器將雜點去除
imagej的指令
Process/Noise/Despeckle...
Process/Filters/Median..




影像邏輯運算
通常使用or運算做聯集,例如前述的過程得到一個以葉子為主的二值化影像,和一個以莖為主的二值化影像,接下來就可以用這兩個影像做聯集,也就是or運算,將兩個影像合併成一組。
imagej的指令
Process/Image Calculator


影像相減與相加
影像相減subtract,常用於減去背景,例如兩張圖片一張為沒有植株的盆栽,另一個是有植株的,兩者影像相減可得到植株的影像(但通常還需要再處理才能使用)
影像相加add,例如產生兩張圖,一張是邊緣強化的(透過sobel運算),另一張是增強對比的(經過laplace濾波),把兩張影像相加,可以得到兩者的優點,邊緣清楚且對比明顯
imagej的指令
Process/Image Calculator
Process/Find edge
Process/Enhance Contrast


從遮罩mask取得原圖目標物的選取區
例如有一個已經經過多次處理得到的二值化影像,在目標物的部份已經有了mask,可以將此mask轉成selection,再於原圖中套用此選取區
imagej的指令
Edit/Selection/Creation Selection(使用於二值化的影像,可得到選取區)
Edit/Selection/Restore selection(將其他影像的選取區套用於此)


以下用PlantCV的範例來做imagej的示範,共有兩個,並且用影片來說明。
---

---

流程:
1.原圖轉HSB --> 取得S channel影像 --> 設定閾值做二值化影像 --> 中值濾波去雜點  -->得到s_blu圖
2.原圖轉Lab --> 取得b channel影像 --> 設定閾值做二值化影像 --> 得到b圖
3.將s_blu圖和b圖做or運算,取得聯集後的影像,得到mask圖
4.將mask圖轉出選取區,套用在原圖上(或加入ROI manager)


第二個示範是紅外線攝影下的植物,有兩個影像,其中一個是有植物的,另一個是沒有植物的

--
--
流程:
1.原圖減去背景 -->  設定閾值做二值化影像  --->得到圖a
2.原圖做對比增強,利用laplace filter ,再用原圖減去這張圖,或是直接用Process/Enhance Contrast  --> 得到圖b
3.原圖做Find Edge尋找邊緣  --> 中值濾波  -->圖片反相 Edit/Invert   --> 得到圖c
4.圖b和圖c做add運算--> 設定閾值做二值化影像  ---> 侵蝕  --> 得到圖d
5.圖a和圖d做 or 運算 -->得到遮罩mask,再轉選取區套用在原圖上




形狀樣式(patten)的量化分析-使用imagej的套件 PAT-GEOM

最近看到一篇論文,是發展了一套Imagej的套件PAT-GEOM,可以用在分析生物身上的斑紋、花紋或條紋(以下這些統稱為樣式,patten)。

論文網址:
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13131

我仔細翻閱實作後,覺得是很有用的工具,以往我們在分析這些樣式的時候,通常也只分析面積有多大、顏色是什麼,再不然就只能說這個和那個「看起來」像,或是「看起來」比較亂...。其實不是我們分析的時候不想量化,而是根本舊不知道怎麼量化這些東西。我很高興能看看到這套工具,而且論文和guide都寫得很清楚,不僅能學習工具怎麼使用,也知道測量的原理和分析之後可以拿來做什麼,同時也提供了分析的例子,讓使用者能弄清楚工具的使用。

這套工具能檢測量化的項目有七項,包括樣式的形狀、方向性、尺寸、對比、分佈、整體樣式的分佈、隨機性等。在作者的網站上可以下載這個套件,也有使用說明和範例可以下載
http://ianzwchan.com/my-research/pat-geom/

安裝順序
1.安裝imagej(或是fiji)
2.將下載的套件解壓縮後放在imagej的plugins資料夾中
3.前往以下網址,下載橢圓傅立葉分析的套件。下載後一樣放在plugins資料夾
http://imagejdocu.tudor.lu/doku.php?id=plugin:analysis:fourier_shape_analysis:start
4.啟動imagej
5.在imagej工具列的Plugins最底下可以找到PAT-GEOM的選項


以下是這套工具能夠分析的項目說明,以及分析的例子。
除了以下文字說明外,我也錄製影片來說明分析方式和例子
---

---




一、形狀(Marking shape)
分析方式:利用橢圓傅立葉分析(Elliptical Fourier analysis)將樣式的形狀做量化。先將要分析的樣式圈出後成封閉的曲線,然後用橢圓傅立葉分析去用大量的函數擬合那個封閉曲線。這些擬合的函數都會有量化的數字,這些數字越接近,代表樣式越相似
例子:
比較杜鵑鳥和其宿主的蛋殼表面花紋
比較兩族群生物的樣式平均形狀,例如長頸鹿
鑑定有特殊色彩斑塊的個體(例如鯨鯊)
比較螃蟹甲殼上斑塊和背景的差異


二、樣式的方向性(Marking shape directionality)
分析方式:各個樣式的長寬比例和方向
例子:
比較食蚜蠅和黃蜂的身體條紋
比較動物身上的條紋和其背景的條紋(例如斑馬)
測量在基因調控下或是選汰壓力下的蝶翅眼點形狀
測量老虎的條紋的差異


三、尺寸(Marking size)
分析方式:除分析樣式的面積外,也會計算形心(幾何中心,Centroid)到樣式邊緣的距離平方總和,再平均之後取平方根的值(averaged centroid size)
例子:
比較人造的獵物和其模擬的對象差異(例如帝王蝶的捕食者實驗)
比較兩個不同族群的生物斑點尺寸的差異(例如七星瓢蟲的斑點)
比較動物和其背景的形狀尺寸


四、對比(Pattern contrast)
分析方式:使用變異係數(Coefficient of variation),用標準差除平均值的值。當對比越明顯,變異係數就越大
例子:
測定比目魚身上的色彩斑塊是否符合隨機取樣的背景基質
比較某一動物的兩個不同區域差異(例如會變色的花枝)


五、分佈(Marking distribution)
分析方式:像素陣列(Pixel matrix)。把生物的整個表面的樣式做成二值化影像(樣式處為1、其餘為0),成為Pixel matrix。把不同個體體表的樣式都做成Pixel matrix之後,可以疊加在一起成為單張圖,即可觀察「這群生物合起來的樣子長什麼樣」
例子:
將一整個族群的生物的「平均樣式」視覺化,例如把一群螃蟹的背甲做成pixel matrix看「這群螃蟹背甲長什麼樣」,觀察不同族群之間的差異
設計擬真的獵物,例如設計海參的警戒色樣式,取得一群海參的平均像素陣列後,觀察海參的平均樣式長什麼樣。


六、整體樣式分佈的方向性(Marking distribution directionality)
分析:分析體表全部的斑塊的角度和排列方向
例子:
在選汰壓力下,某個族群生物的斑塊是否會發展成線性排列(例如某些魚類體表的斑塊排列,或是蝶翅上的斑塊、眼斑)
比較兩種生物是否有相似的身體形狀

七、樣式的隨機性(Pattern randomness)
分析:將樣式存成GIF,用檔案的尺寸來表示隨機性。由於GIF檔案壓縮的特性,越隨機的圖樣,檔案尺寸會越大,因此可以用檔案的尺寸來代表隨機性
比較同種生物的不同型態的樣式(例如彩虹䗉螺Umbonium vestiarium此種生物具有多種不同的斑紋)
測定擬態的品質(例如杜鵑鳥和其宿主的蛋殼、岸蟹和其棲息的背景)

2019年2月3日

imagej中,進行二值化影像處理的一些方法

在imagej處理影像時,這個二值化影像處理的流程,我覺得是蠻實用的,包括形態學影像處理,如侵蝕、膨脹、開運算、閉運算、分水嶺、孔洞填充、邊界提取等。


在分析「影像中感興趣的位置的面積有多大」這樣的問題時,常常會遇到區塊之間彼此相連,無法分割計算,或是影像中有雜點、空洞,這時候我們需要進行影像的預處理來解決問題。
首先須將影像轉換成二值化影像,讓影像只有黑或白兩種顏色。在這之前要先決定「哪些顏色成為黑色,哪些成為白色」,所以需要一個界線數值,超過這個就變成黑(或白)。這個界線的數值稱為閾值,這個步驟在【Process/Binary/Make Binary】中進行

接下來需要一些形態學的處理,包含侵蝕、膨脹、開運算、閉運算。

侵蝕【Process/Binary/Erode】是將區塊的邊界向內收縮,所以小雜點就會消失。

膨脹Process/Binary/Dilate】是使邊界向外膨脹,可以填補區塊內的空洞。

開運算【Process/Binary/Open】,是先侵蝕再膨脹的過程,可以消除小的雜點,或是斷開兩個有細線相連的目標物。

閉運算【Process/Binary/Close】,是先膨脹再侵蝕,可以用來填充目標物內的細小孔洞,連接兩個鄰近的目標物。

提取影像的邊緣【Process/Binary/Outline】,實際的運算其實是把原圖減去侵蝕後的影像,如此可以得到影像的邊緣。

【Process/Binary/Watershed】,用在把兩個兩團相連的東西切開來,比方說細胞分裂時,中間會有牽絲的部份,這個指定可以斷開連結

【Process/Binary/Fill Holes】,顧名思義就是把影像中的洞填起來,將二值化影像中被黑色包圍的白色轉成黑色。

2018年12月24日

使用imagej的plugin TrackMate追蹤細胞

最近一個網友來信問用imagej追蹤細胞的問題,他的材料是一群T cell。在fiji的plugin中,有數個做tracking的plugin,設定最簡單的應該是MTrack 2,不過我一直沒有用這個東西做成功,而且MTrack2的網路資料也好少,很難找到參考文章。

所以後來我用一個設定很複雜的plugin,叫做TrackMate。雖然設定很複雜,但是還好有夠多的文章可以看怎麼做。例如這個 Getting started with TrackMate

做細胞或是物體追蹤的流程,大致上是這樣,要先偵測追蹤的物體在哪,然後把不同frame的同一物體辨識出來成為一個track。

其中還會遇到一些問題,比方說有些物體在某些frame會消失或是快速移動,想像一下拍攝車輛的影片,可能車輛會突然加速或被遮擋。或是某些物體會在途中突然變成兩個,例如細胞分裂的時候,從一個細胞變成兩個細胞。上述這些問題在TrackMate中,都有解決的方式,不過以下的教學影片就沒講到這個東西,反正還沒遇到。


在影片中會提到偵測物體的演算法,有幾個如下。不同的方法偵測的物體,各有所長,也些適合偵測大的,像是The Downsample LoG detector,而DoG detector就適合偵測小的。

  1. DoG detector:對於小尺寸(5 pixels)的追蹤很快
  2. LoG detector:5~20 pixels
  3. The Downsample LoG detector:大於20 pixels
詳細的處理就請看影片

2018年11月27日

用imagej的z project分析小魚尾鰭血液流動影片,...

數週前第一次用imagej的Z project的standard deviation模式做影片疊圖,覺得在處理影片上能夠看到很特別的資訊。(那次的經驗是處理小p的影片,請見此文章

這個standard deviation的特色是讓影片變動的部份變得明顯,而不變的部份則不明顯。

最近讓學生拍小魚尾鰭的血液流動影片時,突然想起可以用imagej來分析看看,也許可以血管抓出來?不過這次分析的影片是用一年前拍攝的孔雀魚尾鰭血液流動,可見微血管


流程就如上次做《用Imagej疊圖看水波干涉的腹線(antinodal line)》那樣,到最後一步選的是standard deviation。

影片長度很長,代表影格也很多,不過並不是每個影格都適合拿來做處理,因為我需要的穩定拍攝的畫面,也就是魚鰭不動,純粹只有血液流動的畫面。

以下分別是取出100張、500張、1000張、2000張影格疊出的畫面。從100張影格疊出的畫面,的確很容易就看到近乎平行的數條血管,而且也能夠看出橫向的微血管。隨著疊加的影格越多,我發現血管之外的好幾個圓圈圈也越來越明顯,意味著在長時間的變化下(約1分鐘),那個圓圈圈是會變動的!


100 frames

500 frames

1000 frames

2000 frames


欸?引發我的好奇心,那到底是什麼?回去影片仔細觀察,原來是魚鰭上的色素細胞。在觀察血液的時候,根本就沒有留意到這個構造是會改變形狀的。當我再仔細觀察影片,才發現它們可真的會變大變小的,而且有件事讓我疑惑,本來我以為所有的色素細胞改變會是一致的,不過仔細去看才發現是有的會變大,有的會變小,看來有些事情可以拿來研究研究了。

2018年11月17日

imagej分離色頻後選取區域的方法

附圖是一個放射醫學的研究生問的問題,想要從這張圖裡把綠色的弧形區域選取出(第二張圖那樣的區域)

這樣的問題分成幾個部份:
1.怎麼圈出綠色的框框?
2.圈出綠色框之後,中間的空洞怎麼填滿?




第一個問題就用threshold就行,不過像這類黑白圖上再套疊彩色的圖,我發現把圖片轉成Lab stake的效果會比較好。(請看影片示範)

而第二個問題,會用到幾個二值化影像處理的指令,像是Dilate膨脹、Erode侵蝕、Fill Hole填洞或是Despeckle去雜點。(一樣請看影片示範)

最後得到的二值化影像(如上面第二張圖,就可以用make selection去轉成選取區)


2018年11月5日

用imagej加手機前鏡頭偵測脈搏






有很多手機app都可以利用手機的前後鏡頭來偵測脈搏,只要手指蓋住鏡頭,App可以自動偵測出脈搏。其原理是計算每幀影格的像素強度變化,可能做些濾波和數值處理。

其實好幾年前就想自己寫個程式來偵測看看,但是當時對於電腦視覺根本不熟,現在想想其實imagej就可以做處理了。

讀入影像
1.使用fiji,先在Help/Update..,更新imagej,然後在Imagej Updater的視窗左下角選擇Manager Update sites,勾選FFMPEG,然後Apply Changes。這會加入一個匯入mp4影像的外掛。安裝後重新啟動fjij



2.開啟fiji,在File/Import/Movie(ffmpeg),把拍攝的影片匯入成stack。
3.由於影片的前後不一定是需要的部份,比方說還沒拍到手指,或是手指已經離開。因此需要剪接stack。不過imagej對於stack內的slices刪除不是很好用(或者是說我還沒找到好方法),所以建議是用複製影格的方式來做,比方說第10影格到第90影格才是需要的部份,那就在Image/Duplicate...把要的影格複製出來成一個stack。

繪製像素強度變化圖(即脈搏圖)
1.利用方框選取工具,然後選擇Image/Stack/Plot z-axis profile,即可以畫出如下圖的脈搏圖




如果不嫌費工夫,也可以把stack在分離出不同色頻來分析,例如這個是分析R Channel的強度變化,雜訊會少很多。方法有點複雜
1.利用Image/Type/RGB Stack,可以把原始的stack分離出不同的色頻,然後再用Image/Stack/Plot z-axis profile來繪製。




同樣也可以分離出HSB,這個是亮度B( brightness )的圖



以下是處理過程的教學影片

2018年10月30日

用Imagej疊圖看水波干涉的腹線(antinodal line)

前幾天看到小p拿出他的兵器-做水波干涉實驗的兵器,錄製了一段影片。突然我想到我在這篇《imagej將時間資料視覺化》寫到imagej可以將時間資料濃縮,可以從中看出特別的時間資訊。

一時手癢就來把小p拍攝的影片用imagej疊圖分析看看。

1.先用youtube-dl把影片用mp4格式下載回來。
2.用fiji(特製的imagej)的file/import/Movie(FFMPEG)將mp4檔讀成stacks
3.用Image/Stacks/Z Project製作疊圖

以下分別是用Min intensity和standard deviation製作的疊圖






Min intensity是把影像中每個隨著時間變化的像素,用最小的數值(相當於最暗的)去疊出來。

而standard deviation的圖我是第一次用這疊,發現還蠻有趣的。疊出來的影像中,黑色的部份代表標準差很小,意思就是數值變動不大,所以在影片中沒什麼變化的點就會呈現黑色。而動得越厲害的部份,就會越彩色。所以左方的水波產生器顏色就很彩色。

把影片做成疊圖來分析,在看水波干涉時產生的腹線(antinodal line),如中央腹線、第一腹線、第二腹線等可是特別明顯呢。