新四季網

圖象拼合方法與設備的製作方法

2023-10-09 11:44:19

專利名稱:圖象拼合方法與設備的製作方法
相關申請的交叉引用本申請要求題為「Multidimensional Image Mosaics」的美國臨時專利申請續列號60/220,025的優先權,該申請於2000年7月21日提交,其內容通過引用包括在這裡。
關於聯邦政府資助研發的聲明本發明部分得到美國政府「國家科學基礎研究獎NO.IIS-00-85864」的支持,因而美國政府對本發明擁有一定權利。
背景技術:
發明領域本發明一般涉及圖象拼合,尤其涉及應用拼圖技術擴展圖象動態範圍和/或測定接收自景像的幅射信號的附加特性的方法與系統。
相關技術的說明成像系統(如攝像機)的主要局限性包括視場、動態範圍、光譜(如色彩)分辨度和景深(即在圖象平面中保持景點正確聚焦的距離範圍)的限制。另外,常規成像系統一般只測量作為光接收方向的函數的入射光的光強,不能測量其他特性,諸如深度(如目標與攝像機的距離)、以及光的偏振態,它們適用於遠程識別材料,形狀與照明條件並適合反射分析,而且,常規攝像機的測量質量較低,如在一般CCD攝像機中,光強解析度只有8位,光譜解析度極差,僅包括三條寬帶通道,通常為紅綠藍通道。
為克服上述諸局限性,即便作了種種努力,得到的系統仍很複雜,而且只解決了一個小問題而不及其餘,如成像分光儀雖在譜維度上有高的解析度,但未擴展傳感器的光強動態範圍。
在不損害空間解析度的情況下,獲得大視場圖象的一般方法是應用「拼圖法」。這種技術將較小的圖象組合成視場更寬的較大圖象,而各小圖象覆蓋一不同的景觀(view of the scene)。該法一直應用於各種科學領域,諸如射電天文學、合成孔徑雷達(SAR)遙感、光學觀測天文學和地球與太陽系其他天體的光學遙感。為適應攝像機隨機移動,近年開發了若干算法,而且這類算法使拼圖應用於電視攝像機。在重疊了較小部分圖象的區域中,原始數據經處理可提高其空間解析度。然而,就各像素的光譜、偏振與亮度而言,常規的拼圖技術不能提高解析度(如動態範圍)。若對一系列圖象引入視差,則可利用拼圖法恢復深度。但與用聚/散焦插入信號估算深度的方法相比,視差法往往不強健,而且更複雜。
非線性檢測器可用於擴展圖象的光學動態範圍,如業已製造的CMOS檢測器可以(1)得到是光強對數的電子輸出信號,或(2)組合兩副不同積分次數的圖象。這類傳感器的光強動態範圍的量級為1∶106,能對大(即高幅照度)信號作不飽和檢測。不過這類器件中的光強信息受壓縮,以便對高光強範圍作採樣(稀疏地),檢測器應用了較低光強專用的量化電平,因而輸出的光強解析度仍只有8~12位。
對較高光強具有較低透射比的非線性透射比硬體,可擴展任一給定檢測器的動態範圍。然而,光強動態範圍仍要按檢測器有限的解析度作量化,即普通CCD的8位解析度經簡單的非線性延伸而覆蓋更高的幅照度範圍,因而非線性壓縮犧牲了較低光強範圍的解析度。
自動增益控制(AGC)在視頻與數碼攝像機中很常見,且類似於靜像攝像機的自動曝光,但其一大缺點在於它的作用是全面性的,結果增益設定對圖象的某些部分可能過高,而對其他部分則過低。例如,亮點在較暗圖象內會飽和,而暗點在較亮圖象內則暗得無法正確檢出。拼圖可由序列構成,當掃描景像時,AGC自適應地改變傳感器增益。然而,儘管該技術擴展了一些動態範圍,但是這類方法仍無法正確地測量大部分暗圖象的亮點和大部分亮圖象的暗點。
在業餘與專業攝像中,普遍在攝像機上安裝空間變動的濾光器,但這種濾光器主要用來改變原始圖象以產生特殊的視覺效果,未曾與解析度提高算法一起用過。
有人曾提出運用一組不同曝光的多幅圖象來擴展圖象各像素的動態範圍。一種這樣的方法是對各像素估算與其多個樣本的數據最合適的值。另一方法是對各像素選擇使局部反差最大的值。然而,這些方法用固定的攝像機拍攝圖像序列,所以並不擴展視場。
另一種方法應用了小濾光器的拼合陣列,覆蓋成像器的檢測器陣列,各濾光器覆蓋一特定的檢測器像素,得出的空間不均一的遮蔽可調製撞擊檢測器的光。為擴展光強動態範圍,可對傳感器陣列覆蓋一中性(即與色彩無關)密度濾光器的空間拼合陣列。但這種結構為了擴展動態範圍,犧牲了空間解析度。對檢測器覆蓋一拼合的濾色器陣列,能獲得光譜信息,但這種結構為獲得某種光譜解析度(即色彩信息),犧牲了空間解析度。此外,還可對檢測器覆蓋一拼合的以各種不同方向定向的線性偏振器,但這種結構為獲得偏振信息,犧牲了空間解析度。
對檢測器陣列覆蓋一種空間變化的濾譜器,即其光譜通帶在檢測器的垂直和/或水平視角上變化的濾譜器,可得到高解析度濾譜作用。在這種系統中,視場內不同的點被不同地過濾。跨過景象掃描攝像機的視場,得到各點的光譜。但將濾譜器直接放在檢波器陣列上,因難以改變濾譜作用的有效特性或測量接收自景象的光的其他特性,而降低了系統的靈活性。
若線性掃描器逐行掃描景象,可不犧牲空間解析度而得到譜信息,如在三線掃描器中,用紅綠藍濾色器連續檢測圖象的各線性部分。遙感工作中常用的Pushbroom攝像機也如此操作;色散元件把每條景線衍射到2D檢測器陣列上,結果在多條譜通道中同時測出每條線。然而,這類掃描器與Pushbroom限於定速的一維(1D)掃描,而且這種系統形成的圖象沒有凹陷,整幅圖形利用同樣的檢測器特性掃描。相應地,為拍攝大的視場,要求多次採集,因為每次採集只能捕獲1像素寬的柱。
圖象一直用不同的聚焦設定值拍攝,然後組合起來生成大景深的圖象。一種應用傾斜傳感器的方法,在擴展視場的同時成功地拍攝了所有聚焦的景點,但是並未擴展圖象動態範圍。
在光學行業,在成像系統前面或內部轉動空間上變化的遮光器與十字線是常用的方法,但是這類系統要求成像器具有在圖象採集期間移動的附加內部或外部部件。

發明內容
因此,本發明的一個目的是提供一種能擴展視場並擴大亮度動態範圍的成像技術。
本發明的另一目的是提供一種既可擴展視場,又能提供有關景象的光譜、偏振和/或深度信息的成像技術。
各種目的由下述本發明諸方面實現。
根據本發明的一個方面,增強解析度數據的生成方法包括一種成像方法,包括用成像器作第一組組測量以生成第一像值的第一步驟,第一組測量包括至少測量一次來自第一景區的第一射線束的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束具有第一強度靈敏度特性,而且成像器對第一主射線射線束的強度具有第一動態範圍;用成像器作第二組測量以生成第二像值的第二步驟,第二組測量包括至少測量一次第一景區發射的第二射線束的強度,第二射線束在成像器參考幀內具有第二主射線,第二主射線與第一主射線不同,成像器對第二主射線射線束具有第二強度靈敏度特性,第二強度靈敏度特性與第一強度靈敏度特性不同,成像器對第二主射線射線束的強度具有第二動態範圍;而且對第一和第二像值作拼合操作,以至少對第一與第二射線束的強度之一生成與成像器的第三動態範圍相關的第三像值,而第三動態範圍至少大於成像器的第一與第二動態範圍之一。
根據本發明的另一個方面,成像方法包括用成像器作第一組測量以生成第一像值的第一步驟,第一組測量包括至少測量一次第一景區發射的第一射線束中至少一個選擇偏振分量的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束具有第一偏振靈敏度特性,而第一偏振靈敏度特性包括對偏振角在第一角範圍之外的信號分量減小的靈敏度,第一射線束至少一個選擇的偏振分量在第一角範圍內有一偏振角;用成像器作第二組測量以生成第二像值的第二步驟,第二組測量包括至少測量一次第一景區發射的第二射線束中至少一個選擇的偏振分量的強度,第二射線束在成像器參考幀內具有第二主射線,第二主射線與第一主射線不同,成像器對第二主射線射線束具有第二偏振靈敏度特性,第二偏振靈敏度特性對偏振角在第二角範圍之外的信號分量減小了靈敏度,第二射線束中至少一個選擇的偏振分量在第二角範圍內有一偏振角,而第二與第一角範圍不同;移動成像器的第三步驟,包括下列之一在第一與第二步驟之間相對第一景區轉動成像器;在第一與第二步驟之間相對第一景區平移成像器;和運用第一與第二像值判定第一與第二射線束之一的偏振態。
附圖簡介通過結合示明本發明諸示例性實施例的附圖所作的詳述,將明白本發明的其他目的、特徵與優點,其中

圖1是本發明拼圖法示例步驟的流程圖;
圖2是本發明示例校正步驟的流程圖;圖3是本發明另一示例校正步驟的流程圖;圖4是本發明示例拼圖步驟的流程圖;圖5是本發明示例拼圖技術的框圖;圖6是本發明示例拼圖系統的框圖;圖7是本發明另一示例拼圖系統的示圖;圖8是本發明又一拼圖系統的示圖;圖9是本發明再一拼圖系統的示圖;圖10是本發明另一示例拼圖系統的示圖;圖11是本發明另一拼圖系統的示圖;圖12是本發明又一示例拼圖系統的示圖;圖13是本發明再一示例拼圖系統的示圖;圖14是本發明又一示例拼圖系統的示圖;圖15是本發明再一示例拼圖系統的示圖;圖16是本發明又一示例拼圖系統的示圖;圖17是本發明另一拼圖系統的示圖;圖18是本發明又一拼圖系統的示圖;圖19是本發明再一示例拼圖系統的示圖;圖20A是本發明示例拼圖步驟的空間與強度範圍的曲線圖;圖20B是本發明另一示例拼圖步驟的空間與強度範圍的曲線圖;圖21是本發明又一示例拼圖步驟的空間與強度範圍的曲線圖;圖22A是本發明示例濾波器的密度分布與相應有效遮蔽特性的曲線圖;圖22B是本發明另一示例濾波器的密度分布與相應有效遮蔽特性的曲線圖;圖22C是本發明又一示例濾波器的密度分布與相應有效遮蔽特性的曲線圖;圖23是本發明示例拼圖系統的示圖;圖24A是本發明示例成像器附件的有效密度分布曲線圖;圖24B是圖24A中分布曲線的對數函數曲線圖;圖25A是本發明示例濾波器的中心波長分布與相應有效遮蔽特性的曲線圖;
圖25B是本發明另一示例濾波器的中心波長分布與相應有效遮蔽特性的曲線圖;圖25C是本發明再一示例濾波器的中心波長分布與相應有效遮蔽特性的曲線圖;圖26A是本發明一示例濾波器的截止波長分布與相應有效遮蔽特性的曲線圖;圖26B是本發明另一示例濾波器的截止波長分布與相應有效遮蔽特性的曲線圖;圖26C是本發明再一示例濾波器的截止波長分布與相應有效遮蔽特性的曲線圖;圖27是本發明一具有高通濾波器陣列的示例成像器的一部分靈敏度特性曲線圖;圖28是本發明一具有窄帶濾波器的示例成像器的靈敏度特性曲線圖;圖29是本發明一具有一組寬帶濾波器的示例成像器的靈敏度特性曲線圖;圖30A是本發明一示例偏振濾波器陣列示圖;圖30B是本發明一示例偏振濾波器製造示圖;圖30C是本發明另一示例偏振濾波器製造圖;圖31是一示例攝像機的焦點特性圖;圖32A是本發明一示例折射元件的折射分布曲線圖;圖32B是本發明另一示例折射元件的折射分布曲線圖;圖32C是本發明又一示例折射元件的折射分布曲線圖;圖33是本發明一組光學元件的示例結構圖;圖34是本發明一示例光學元件圖;圖35A是本發明一示例圖象部分配準步驟圖;圖35B是本發明另一示例圖象部分配準步驟圖;圖36A是本發明一示例遮蔽的衰減分布曲線;圖36B是圖36A中衰減分布的對數特性曲線;圖37是本發明執行拼圖算法的計算機系統圖;和圖38是用於圖37中計算機系統的處理器部分的框圖。
圖中,除非另有說明,一般用同樣的標號與字符標明所示實施例相同的特徵、單元、元件或部分。而且,現在雖然參照諸附圖並結合所示實施例對本發明作詳細描述,但是可對描述的實施例作變化與更改而不違背所附權項所限定的本發明的實際範圍與精神。
發明的詳細描述根據本發明,成像器經配置,可使其一個或多個檢波特性在垂直和/或水平視角而變化;該成像器可以包括電子靜像攝像機或視頻攝像機等運動圖象攝像機。例如,經配置,成像器的視場左側與右側的特性可以不同,或者視場上部與下部的特性可以相異。例如,成像器的不均一(即空間變化)靈敏度特性可以包括對景象亮度的靈敏度、對特定色光的靈敏度、對特定偏振角和/或焦距(即目標聚焦距離)光的靈敏度。在連續的快照或幀之間轉動和/或平移成像器,可將連續的快照或幀組合成視場更寬的較大圖象。這種將不同景觀的組合稱為「拼圖」。此外,若攝像機在快照或幀之間的移動足夠小,則在每次通過攝像機的不同視場部分時,可多次拍攝某些景區。由於成像器每一部分的視場有一不同的靈敏度特性,所以可用各種成像靈敏度特性攝得最終的「多採樣」景像部分(即多次採樣部分),從而得到景像各部分的附加信息。例如,可用寬範圍的強度靈敏度特性多次採樣某一景部,如在用透鏡上裝有空間不均一衰減器的攝像機跨越景象作垂直或水平(或任意角度)隨動拍攝時,攝取多個幀,這樣就有效地攝取各景部,並擴展了動態範圍。此外,在用具有空間不均一色彩靈敏度特性的成像器作多次快照時,通過跨越景象的隨動拍攝可獲得各景部的譜信息。同樣地,應用具有不均一偏振靈敏度特性的成像器(如帶不均一偏振濾波器的攝像機)可獲得偏振信息,應用有不均一焦距的成像器可獲得深度(即距離)信息。
如圖5的框圖所示,本發明的系統一般可看作有兩個部件。該系統包括一硬體部件,含有用於攝像504的成像器(如攝像機)502。圖象504由圖象處理部件執行軟體處理,後者包括一種或多種圖象處理算法506,提供增強的輸出圖象508和/或有關被成像景象特性的信息510。圖6更詳細地展示了此系統。成像器502包括一個或多個濾波器604,濾除進入成像器502的光信號(如光束)632。成像器502包括成像光學元件606和檢測器陣列608,諸如CCD或CMOS圖象檢測陣列。檢測器陣列608產生可由電路610處理的信號。為轉動和/或平移成像器502,可將成像器502裝在機械化的轉動和/或平移機架630上,和/或由飛機之類的運動平臺載運。
模/數轉換器612處理攝像機電路610產生的模擬信號,把模擬信號轉換成可存入幀存儲器614的數位訊號。圖象用處理器618分析處理,它可以是計算機的處理器,執行本發明的各種算法。處理器618還可用來「提供」(即產生)新景象。轉換器616可將圖象轉換成可由外部設備602處理的格式而生成輸出圖象508。外部設備602可以包括例如轉換器620,用於將圖象轉換成視頻信號、印表機信號或其他輸出格式。處理器618還可提取有關被成像景象特性的信號622。下面詳細描述上面諸技術的實例。
圖8示出一例按本發明作拼圖的成像器。成像器712包括有孔徑710、物鏡802和可以是如CCD檢測器陣列的圖象檢測器708的攝像機702,還包括一空間變化的中性(即與波長無關)密度濾波器704,在圖示實例中,靠近上部816的衰減較小,靠近下部820的衰減較大。成像器712分別從景點A、B、接收射線束810、812與814(本例為光束)。在成像器712的參考幀內,若干主射線例如射線804、806與808,限定了成像器712接收的射線束810、812與814或其他輻射信號組的各自的方向。如光學中所周知,一般認為,光束的主射線可以限定束傳播的路徑。在圖示實例中,束810、812與814的主射線分別是成像器712的主射線804、806與808。再者,儘管圖8隻示出三條主射線804、806與808,但是成像器在理論上在其參考幀內擁有無數條主射線。此外,雖然圖8的示例成像器712主要用於圖象光,但是本發明的技術還適用於任何電磁輻射或其他輻射的成像,包括但不限於紅外(IR)輻射、x射線輻射、合成孔徑雷達(SAR)信號、粒子束(如電子顯微鏡的電子束)與聲(如超聲)輻射。
分別包括景點A、B、C發射的信號的信號組810、812與814,被物鏡802分別聚焦到檢測器708上的點A』,B』與C』上。成像器712可以繞旋轉矢量 旋轉,使成像器712沿任一主射線804、806與808或沿成像器712的視場中任一其他主射線接收任一射線束810、812與814。例如,可用如圖8那樣定向的成像器712作第一次快照,此時沿圖示的主射線804接收信號組810。然後,成像器712以逆時針方向旋轉,從而沿主射線806接收射線束810。接著作另一次快照。於是,成像器712再旋轉,從而沿主射線808接收信號組810,並可作第三次快照。相應地,在各三次快照中,射線束810通過濾波器704的三個部分816、818與820接收,因而在檢測器708的三個點A』,B』與C』接收。另應指出,為了跨越景象隨動拍攝,雖然圖8示出了旋轉的成像器712,但是成像器712也可像圖7那樣平移。其實可用任一種運動使成像器712跨越景象706隨動拍攝,讓成像器712沿其視場的任一主射線接收景點A的輻射信號。對景象706作了多次快照後,就可將快照用於本發明的拼圖步驟。例如,如下面參照圖1詳述的那樣,可以用拼圖步驟擴展代表接收自各種景區706的光信號組強度的圖象的動態範圍。
圖1示出本發明示例性拼圖步驟。較佳地,該步驟包括校正步驟100,下面詳細討論其實例。圖7~19、23和31中標為712的成像器(如帶空間變化衰減濾波器的攝像機)對被成像景象作第一組測量(步驟102)。在成像器712的參考幀中,有第一主射線(如圖8~18中的主射線804),對於沿第一主射線接收的信號,如第一主射線光束,成像器712具有第一強度靈敏度特性與第一動態範圍。成像系統具有第一強度靈敏度特性和第一動態範圍。例如,在成像器712視場內,第一主射線對應於特定的觀察方向,並對該部分視場覆蓋一第一衰減量衰減濾波器。成像器對通過衰減濾波器特定部分的光束的靈敏度,取決於濾波器該特定部分的衰減量。成像器對光束的動態範圍一般由光束在其上聚焦的檢測器部分(如CCD單元)的動態範圍決定。
第一組測量的輸出是第一測量值,例如可代表第一射線束或景象內第一區或點發射的其他輻射的強度(即亮度)。第一射線束的主射線與成像器的第一主射線相對應或等同。
移動成像器712(如轉動和/或平移),可拍攝不同的景部(步驟104)。第二景部與第一景部重疊,使第一景點或景區仍在成像器712的視場內,但是來自第一景區的光射線束沿成像器712的參考幀內的第二主射線(如圖8~18的主射線806)接收。成像器作第二組測量步驟106),產生第二測量值,該值代表來自第一景區的第二光線束的強度。成像器712對沿第二主射線接收的輻射信號組具有第二強度靈敏度特性。例如,若第二主射線通過其第二衰減量不同於第一衰減量的衰減器,則成像器712對第二光線束的靈敏度就不同於它對第一光線束的靈敏度。
對第一和第二測量值作數學運算,就產生第三測量值(步驟108)。例如,若成像器712包括一用於拍攝圖像的CCD檢測器陣列708,就用該陣列的第一單元或第一組單元測量第一光線束而產生第一測量值,並用第二單元或第二組單元測量第二光線束而產生第二測量值。根據成像器712對沿第一與第二主射線接收的信號的特性,第一或第二單元接收的信號可能過亮或過暗,使信號位於該單元的準確範圍以外。若第一測量值表明已準確地測量了第一光線束,但第二測量值指示第二光線束未能準確地測量,則可丟棄第二測量值,而把第一測量值用作第三測量值,並作為代表第一景區的像素的輸出值。若成像器712受第一主射線撞擊的區域與受第二主射線撞擊的區域的特性不同,例如若沿第一主射線接收的光比沿第二主射線的光被衰減得更大,則有效地擴展了成像器712的動態範圍,因為高強度光在沿第一主射線接收時可準確地測量,而低強度光在沿第二主射線接收時可準確地測量。由於第一景區發射的光是沿兩條主射線捕獲的,所以光更容易被檢測器兩區域的至少一個區域準確地測量。相應地,可將第三測量值視作具有第三有效動態範圍,而且大於第一與第二測量值的有效動態範圍之一或二者。
除了上述有效地擴展強度測量的動態範圍的步驟以外,圖1的示例拼圖步驟還可以包括測量接收自第一景區的光線束的其它特性的步驟。這類其它特性可以包括例如譜特性、偏振和/或焦點特性,可用於推導攝像機與各種景象細節的距離。例如,第三組測量可測量來自第一景區的第三光線束中至少一個譜分量的強度(步驟110)。第三組測量產生第四測量值。第三光線束可以沿成像器712的第三主射線接收,該成像器配置成對第三主射線光線束具有第一光譜靈敏度特性。為了選擇第三光線束的譜分量,第一光譜靈敏度特性較佳地包括一具有第一波長靈敏度帶的帶通特性。換言之,所選的譜分量具有位於第一波長靈敏度帶內的一個或多個波長,而且可以擁有或可以不擁有產生高於檢測器噪聲的信號所需的能量。
接著,攝像機轉動或平移(步驟112),作第四組測量(步驟114)。產生第五測量值的第四組測量包括至少測量一次來自第一景區的第四光線束某一譜分量的強度。第四光線束沿成像器712的第四主射線接收。成像器712對第四主射線輻射信號具有第二光譜靈敏度特性,它較佳地包括第二波長靈敏度帶的帶通特性,以便選擇波長位於第二波長靈敏度帶內的諸分量。
此外,本發明的拼圖法還能測量第一景區發射的光的偏振。在這種系統中,第三組測量(步驟110)包括至少測量一次第三光線束中選擇的一偏振分量強度。具有上述第三主射線的第三光線束包括來自第一景區的第三輻射信號。成像器712對沿第三主射線接收的射線束具有第一靈敏度特性,該特性包括對偏振角在第一角範圍之外的光線降低的靈敏度。檢測第三光線束中所選的偏振分量,因其偏振角在第一角範圍內。然後轉動或平移成像器712(步驟112),執行第四組測量,產生第五測量值(步驟114)。第四組測量包括至少測量一次來自第一景區第四光線束所選偏振分量的強度。第四光線束沿成像器712的第四主射線接收。成像器712對沿第四主射線接收的射線束具有第二偏振靈敏度特性,該特性包括對偏振角在第二角範圍外的信號分量減低的靈敏度。成像器712檢測第四光線束的所選偏振分量,因為該分量的偏振角在第二角範圍內。
另外,可用焦點特性不均一的成像器確定第一景區離成像器712有多遠。具體而言,可將成像器712配置成對沿第三主射線接收的光具有第一焦點特性,而對沿第四主射線接收的光具有第二焦點特性。作第三組測量以產生第四測量值,其中第三組測量包括至少測量一次來自第一景區的第三光線束的強度(步驟110)。第三射線束沿第三主射線接收。攝像機的第一焦點特性包括目標聚焦的第一焦距。轉動或平移成像器712(步驟112),作第四組測量而產生第五測量值(步驟114),其中第四組測量包括至少測量一次第一景區第四射線束的強度。第四組信號沿第四主射線接收。若沿第四主射線成像,則攝像機的第二焦點特性包括目標聚焦的第二焦距。第二焦距與第一焦距不同。
對於應用強度靈敏度特性不均一的成像器拼圖,通過獲取攝像機的強度靈敏度特性如何跨越視場而變化的估值有利於對成像器做校正。一種校正技術是用成像器捕獲各種不同的景像與景部,再將檢測器各部分如各檢測器單元產生的測量值相加或求均,結果是一組相對的和/或標定的值,代表沿各主射線的成像器特性。圖2示出的一例校正方法100可用於圖1的方法。與圖1所示方法的其他步驟一樣,圖2所示一系列步驟100由在成像器712的參考幀中接收第一與第二主射線的成像器執行。成像器712測量第一組第一主射線射線束(如光線束)的強度,產生第一組校正測量值(步驟202),用於確定成像器712對沿第一主射線接收的信號的第一強度靈敏度特性的第一估值(步驟204)。該第一估值通過計算第一組校正測量值的總和和/或平均值而確定。成像器712還測量第二組第二主射線射線束的強度,產生第二組校正測量值(步驟206),用於確定成像器712對沿第二主射線接收的信號的第二強度靈敏度特性的第二估值(步驟208)。第二估算通過計算第二組校正測量值的總和和/或平均值確定。
此外,當景區跨越視場運行時,通過對它進行跟蹤可以校正成像器的強度靈敏度特性。例如,圖3示出的校正方法300就應用了沿成像器接收的多束主射線從選擇的景部接收的射線束(如光線束)的測量值。圖3的校正方法100可用於圖3的拼圖方法。在圖示校正方法300中,成像器712作第三組測量,產生第四測量值(步驟302),其中包括至少測量一次第二景區第三光線束的強度。第三光線束沿成像器接收的第一主射線接收。然後,成像器轉動或平移,使成像器沿其第二主射線接收來自第二景區的光線(步驟304)。成像器作第四組測量,產生第五測量值(步驟306),其中第四組測量包括至少測量一次第二景區第四光線束的強度。第四光線束沿成像器接收的第二主射線接收。第四與第五測量值用於評估成像器712的第一強度靈敏度特性(即成像器對第一主射線光線束的靈敏度)和該成像器原第二強度靈敏度特性(即成像器相對於沿第二主射線接收的光線的靈敏度)之間的數學關係(步驟308)。該數學關係通過計算第四與第五測量值之差或它們的比率進行評估。
上述校正法可進一步理解如下。考慮一種只沿x軸改變透射比的中性密度遮蔽M(x)。某景點在圖像k中表示為像點Xk;像點xk的線性化強度為 。同一景點在圖像p中表示為xp;像點xp的強度為 。這兩點都遵循下述關係M(xk)g^p-M(xp)g^k=0.---(1)]]>跟蹤若干圖像中大一些景點可得出許多遮蔽應在各像素x遵循的公式,這些公式可用矩陣形式表示為FM=0。一例矩陣F為 幀數標為p。有些點可能沒有可靠的數據,因而採用附加的平滑度公式有利於使解規範化。例如,使D2M2劣化的公式可以定為尋求LM=0。其中 上述公式通常相互矛盾,因而該算法利用最小二乘方法、耐用統計法或其他方法尋找優化求解法以實現優化。對最小二乘方法,優化解為M^=argminM(MAAM),---(4)]]>其中A=FL,---(5)]]>而β是相對於不一致數據的補償對不平滑解的補償加權的參數。用單值分解(SVD)求出非無效解,再將最大 (即 的最大值)設為1。根據上述公式將M的協方差矩陣估算為Cov(M)=(AA)-1M^AAM^(nr+1-l)-1---(6)]]>式中l為M的元數,nr是A的行數。可以將它視作加權的最小二乘方問題以β加權屬於L的行,而屬於F的行通常對更強的像素(即更大的g)更準確。這相當於應用歸一化行,然後以 加權每行r。相應地,算法應用nr=∑r,cA2(r,c),據此加每行平方的權值。
Cov(m)的對角線給出的M的方差導致 的置信區間。注意,該公式不在logM域內,因而並不嚴厲補償極低M處數據之間的相對不一致或在低M處可能相對顯著的擾動。作為最終的後處理步驟,還要平滑 ,主要對強烈光衰減區域的估算有影響。
應指出,實際上遮蔽自校正以同一圖象序列為基礎,該圖象序列最終經處理而得到代表景象的圖象,此時遮蔽在點x的估值在統計上並非與在特定幀中的點x測得的信號無關。但要指出,根據對若干幀中眾多點的跟蹤,M在各點的估值受到上百甚至上萬個公式(A的行)的影響。因此在配準與融合圖象時,為實用起見,可假定每次特定採樣的信號與估算的遮蔽無關。
本發明的自校正方法已作過實驗測試,測試應用的市售線性可變密度濾波器,長3cm,固定於線性CCD攝像機25mm透鏡的前方約30cm。濾波器的最大密度為2(相當於衰減100倍),儘管有效遮蔽特性因系統的附加漸暈效應而具有較寬範圍。在視場中衰減最小的部分,M接近不變。攝像機在幀間轉動,各點越過視場成像14次。利用M的粗估值,圖象按下述配準方法配準。接著,根據隨機對應不飽和點與非暗點,產生5萬多個公式,用於以614個像素的解析度確定遮蔽特性。圖36A示出上述自校正方法產生的估算遮蔽函數,其對數示於圖36B。
對多項式、S形等參數曲線或仿樣曲線的數據,也可把M確定為最佳擬合解。或者,除了應用MSE判據外,可以應用其他優化法,如增強統計學與選代投影法等。由於M為相乘,所以可將logM用作相加參數,結合I在各景點的估值利用線性優化形式得以優化。
此外,若景象I只沿x軸變化(一維信號),而且各「幀」P相對球坐標系平移tp,則方差和為px(g^p(x)-M(x)I(x+tp))2---(40)]]>對 (x)不飽和。
該算法對平移參數tp、遮蔽M與變量I優化了上述誤差函數。作為優化過程部分作調整的其他參數包括(1)M的平滑度(即M參數化所需的係數數量),和(2)I對幀間的特定移動量(即tp與tp+1之差)的平滑度。I與M較佳地限於正數。公式(7)裡的方差和作為M的函數而加權(使問題非線性),或在對數域中計算該誤差,使上述公式對未知量I與M成為線性優化問題。
遮蔽還可以選代校正,例如校正方法以遮蔽各點的初始估值開始,可用初始評估的遮蔽函數估算各像素強度的融值。融合像值利用上述校正法對遮蔽計算改進的估值,而改進的遮蔽估值可對融合圖象導出較佳的估值,並依次類推。
根據本發明,具有空間變化偏振靈敏度特性的成像器可獲得景象發射的輻射(如光或其他電磁輻射)的偏振信息,不管該成像器是否也具有空間變化強度靈敏度特性。圖4示出一例這樣的方法。圖示方法400應用的成像器712,對沿成像器712第一主射線接收的輻射信號組具有第一偏振靈敏度特性,對沿成像器712第二主射線接收的輻射具有第二偏振靈敏度特性。例如,成像器712包括帶不均一偏振濾波器的攝像機702,該濾波器接納在第一部分視場中有第一偏振角的光,並且接納在第二部分視場中有第二偏振角的光。圖30A示出一例這樣的濾波器。圖示濾波器3010包括三個偏振部分3002、3004與3006,分別發射偏振角為3020、3022與3024的光。濾波器3010還包括可通過任何偏振光的部分3008。圖30A的濾波器3010可在攝像機702內部外或外安裝,形成偏振靈敏度特性在視場左右兩側之間變化的成像器712。當成像器跨越景象隨動拍攝而捕獲景象的多幅快照或幀時,就通過空間變化偏振濾波器3010的部分3002、3004、3006與3008中的一個或多個捕獲各景點。
空間變化偏振濾波器可用幾種不同方法製作,如在圖30B中,這種偏振濾波器3012從盤3018中切割,發射光的偏振角位於與盤3018水平的方向。在圖30C中,空間變化偏振濾波器3016從盤3014中切割,發射偏振角與盤成經向的光。
包括上述濾波器之一的成像器712通常具有空間變化偏振靈敏度特性。在任何場合中,不管成像器如何構成,在圖4的方法400中,都用成像器712作第一組測量而產生第一測量值(步驟402),其中包括至少測量一次第一景區第一光線束中至少一個所選偏振分量的強度。第一光線束的主射線對應於成像器接收的第一主射線。成像器對第一主射線光線束的第一偏振靈敏度特性,包括對偏振角在第一角範圍外的光分量降低的靈敏度。選擇的第一信號偏振分量由成像器檢測,因其偏振角在第一角範圍內。為捕獲景象的不同部分,要轉動或平移成像器(步驟401)。用成像器作第二組測量而產生第二測量值(步驟406),其中包括至少測量一次第一景區第二光線束中至少一個所選偏振分量的強度。第二區域的主射線對應於成像器接收的第二主射線。成像器對第二主射線光線束的第二偏振靈敏度特性包括對偏振角在第二角範圍外的信號分量降低的靈敏度,選擇的第二光線束的偏振分量要予以檢測,因其偏振角在第二角範圍內。
捕獲了一組快照和/或視頻幀後,以幾種不同的方法分析得到的圖象數據。希望以儘可能高的精確度測量景象特性,還希望以最有效的可能方法作必要的測量與計算。在採集的數據量與數據捕獲和分析所需的時間與運算能力之間有一折衷方案。通過在單位視角變化內獲取更多的快照或幀,可實現更大的精密度,不過更大數據量的捕獲與分析更加耗時費錢。捕獲大量圖象的回報正在減少,因為在實踐中,被測的景象特性的自由度通常有限,因而景象的過分採樣會導致冗餘數據。再者,外裝濾波器的效果一般有點模糊,因為濾波器散焦了。因而應用快照之間極小的角度變化並不能獲得足夠的附加信息來證明捕獲與處理數據所需的附加時間與費用是值得的。因此,本發明的拼圖法應該較佳地平衡附加信息與附加時間費用之間的折衷。下面針對應用對強度、波長、偏振和深度具有空間變化靈敏度特性的成像器的拼圖法,討論單位視角變化的有效實用幀速率的測定方法。
例如,考慮一由裝在攝像機上的遮蔽形成其空間變化強度靈敏度特性的示例成像器。設遮蔽的透射比為M,檢測器無遮蔽的光強為I,則濾波後落在檢測器上的光為g(x,y)=M(x,y)l(x,y). (8)觀看高反差景象時,一般用幅值或倍頻程定義強度,為此,通常以倍頻程來配置攝像機孔徑「f光圈」,使每次「光圈」增加對應於加倍的被測強度。在數位相機傳感器中,這在測量的二進位表示法嗯對應於移1位。例如,若8位攝像機在圖象某一像素位置把光強測量00011010,則增大一個光圈將在第二圖象中對應於讀數0011010(0),其中新的最低位是新圖象添加的信息。
例如,在一種實現偶數幅值分度的優化遮蔽中,衰減倍頻程跨越視場線性的變化,即衰減呈指數變化,於是log2M(x)正比於x。在這種結構中,固定掃描增量(如成像器觀察方向的一系列等量變化)將產生被測強度幅值的固定變化,同等地採樣所有的強度範圍。在透鏡前方以某一距離將線性可變密度濾波器裝到攝像機,可近似實現這一特性。但要指出,鑑於暈像、透視與透鏡失真,在logM(x)中不能準確地保持濾波器密度的線性度。
設I是濾波器透射比最大時(即M=1)落在檢測器上的光強(即幅照度)。特定的線性可變濾波器與特定的檢測器共同測定可被系統不飽和檢測的景象幅照度的最小與最大界限。設檢測器能檢測的高於該檢測器噪聲(暗度)的最小幅照度(對指定的攝像機指標而言)為 ,這決定了整個系統能檢測的最小幅照度。設檢測器能不飽和測量的最大幅照度為 。於是,檢測器以倍頻程計的光學動態範圍為DRdetector=log2ImaxdetectorImindetector--(9)]]>通常DRdetector=8位。整個系統能不飽和檢測的最大幅照度是檢測器在最強衰減下(即最小遮蔽M)得出的最大輸出;Imaxsystem=Imaxdetector]]>/最小M。
因此,系統的總光學動態範圍為DRsystem=log2ImaxsystemImaxdetector=DRdetector-log2(minM)=DRdetector+log2max(1/M)---(10)]]>關於如何最有效地採樣景象信息,可假定捕獲的圖象的幅照度範圍在 與 之間。改變觀察方向相當於改變通過其看見某一點的孔徑光圈(允許進入檢測器的光功率),每改變一次全光圈(即衰減改變2倍)相當於在檢測器動態範圍DRdetecd內測量值的、二進位表示法中移1位。對某些場合,在掃描中丟失最小信息可能較佳,此時在完成掃描時,不會「遺漏」每個幀增量獲得的位。但對無冗餘的場合而言,應用「有損」掃描就行了,完成掃描時,只要求無飽和點,而且所有點都高於檢測器的閾值(只要IImindetector]]>)。這種方法可實現最有效的掃描,每次增量可最大程度地擴展測量的光學動態範圍。
例如考慮某一成像系統,其中log2[max(1/M)]<DRdetector,DRdetector=8位,minM=1/8,所以DRsystem=11位。設某景點得出的二進位值在M=1時為11111111,在M=1/2時為10111001,則前者測量飽和,而後者測量不飽和。相應地,M=1/8時(右移2位以上),測量得出00101110。這樣,幅照度的11位表示是0010111001(0)。假定景象包含系統整個光學動態範圍的幅照度值,則保證信息(位)丟失最小的最小冗餘掃描是一種以一個倍數程增量對各點作四次採樣的掃描。若使用了線性可變濾波器,則在下列二者之間存在線性關係(1)連續幀之間的位移,和(2)特定景區被成像的衰減量。
對動態範圍要求不高的場合,例如可對各點只作單次不飽和測量,若IImindetector]]>,測量值就高於檢測器的閾值。此時,只用M=1與M=1/8就夠了。這樣可對上述示例像素產生值00101110(000),而對景象的暗點產生例如00000000011的值。在這種系統中,任一點可用遮蔽的一個極值測量,然後用另一極值測量,可以跳過中間遮蔽值。換言之,對各點採樣兩次就夠了。幀間增量可以很大,只要它小於幀長一半,而且M=1與M=1/8之間的過渡區小於該幀間增量。
現在考慮一log2[max(1/M)]≥DRdetector的系統。包括遮蔽在內,系統的光學動態範圍至少是檢測器本身的兩倍。以一個倍數程增量對每點採樣,保證丟失的位最少。對有損但最有效的掃描而言,測量每點時在連續掃描增量間不重疊強度信息最好。例如,若DRdetector=8位而DRsystem=16位,則每點兩次採樣就足以獲取整個值範圍。第一次採樣高清晰度地檢測系統量化器對其產生標定為 的像值的每點,該值小於256,第二次採樣不飽和地捕獲像值大於255而小於65536的各點。
通常,對「無損失」的有效掃描,各點因採樣1+log2[max(1/M)]次,對最有效掃描,各點應採樣ξ=[DRsystem/DRdetector]次,ξ總是被捨入到最接近的整數。實踐中,希望根據經驗應用上述數字,並且應用較為稠密的採樣速率,因為與較稠密的採樣相關的冗餘度可使幅照度估算不嘈雜,且便於穩定地配準圖象序列內的圖象。還要指出,上述景象採樣方法不是景象採樣的唯一方法,還可採用其他周期和/或非周期的採樣間隔。
本發明的拼圖方法還可在強度域中實現「超解析度」,因為平滑變化的遮蔽可起到模擬計算設備的作用,因而在許多場合中能提供優於單次量化圖象測量的量化精度。例如,考慮一種其整數量化器區分不出像值96.8與97.2的成像器,即這兩個像值都將產生一個輸出值97。運用M=1/2將信號衰減一倍,成像器分別輸出值48與49,因而增強了鑑別力。在另一例中,利用M=2-6的檢測結果,即引用7幅圖象,可區分出95.8和96.2。一般而言,若應用更多圖象,則融合估值的強度解析度更高。[具體而言,要得出精確的結果應如何組合多個採樣?]這是對景象比上述速率更稠密地採樣的另一個優點。[何種速率?]對於測量景象光譜特性的拼圖法,所選視角的優化間隔取決於檢測器各像素位置的強度與波長函數相對于波長的變化速度,這類似於多類型信號的採樣問題若採樣少得難以檢測信號擾動,就會出現混淆。因此,若各景點或景區的光譜平滑,如黑體幅射場合,就可稀疏地採樣。另一方面,若光譜含有細譜線,如原子蒸氣的放射與吸收譜線,該光譜應稠密地採樣,相應地,成像器在幀間的轉動或平移量要很小。
然而,即使來自景象的光譜具有窄譜帶,在測量這些譜帶時,通過使外裝濾波器散焦模糊,也可平滑測得的分布曲線。例如,研究一種線性變化的濾譜器。可將模糊的測量值模擬成譜s(λ)與寬度Δλblur的窗函數h(λ)的卷積。設S(uλ)與H(uλ)分別是s(λ)與h(λ)的富裡葉變換,H是截止頻率≈1/Δλblur的低通濾波器。若濾波器的發射特性H在頻率上無限伸展,可將該截止頻率定義為(1)信號能量極低的頻率,(2)H的第一個零點,或(3)正比於函數H的標準偏差的量。在任一情況下,該「截止」頻率將≈1/Δλblur。要求的採樣間隔為λ(x)sample≈Δλblur/2 (11)可用系統的物理量綱表示上述採樣判據。對於長為L的線性可變濾波器而言,(x)=max+min2+BLx,x-L2,L2,---(12)]]>其中λmax與λmin分別是整個濾波器通過的最大與最小波長,B≡λmax-λmin是濾波器總帶寬。若濾波器散焦模糊模型是寬度Δxblur的窗核的卷積,則ΔXblur=BΔXblur/L,因而sampleBxblur2L---(13)]]>例如,研究圖8的成像器712。若成像器712在無限遠聚焦,則外接濾波器704被核模糊,該核的寬度等於通過孔徑710的目標光束(如光束810、812或814)的寬度。因此,對透鏡802與外接濾波器704圖示結構而言,Δxblur=D,這裡的D為孔徑710的直徑。因此sampleBD2L---(14)]]>上述結果可定性說明如下。若帶B在長濾波器上伸展,使模糊核的大小不重要,即D<<L,則原始光譜(可以有銳譜線)的細節得以保留,因而採樣較佳地應稠密一些。定量地說,從公式(14)可看出,當D<<L時,採樣周期Δλsample變小了。
假定成像器712在圖象採集之間以角增量ΔΦsample轉動,設濾波器704與轉軸的距離為A,例如可以是圖8所示系統的投射中心點0。而且,讓濾波器704與光軸垂直。於是,從景象傳播到投射中心的每條光線在採樣之間位移Δxsample≈AΔΦsample,假定視場全形θ很小,即假定sinθ≈θ。參照公式(11),sampleBLxsampleBLsample---(15)]]>比較公式(14)與(15),注意f光圈數f#≡F/D,F是焦距,可發現sampleD2A=F2Af#,---(16)]]>與濾波器的量綱或總帶寬無關。對應於公式(16)的定性說明是較大的孔徑D引起較大的散焦模糊;相應地,要求較少的採樣,因而可增大ΔΦsample。再者,對給定的轉動增量,濾波器與投射中心之間較小的距離A,導致通過投射中心O的光線的波長λ的變化較小,從而減小了得到的有效波長採樣周期。因此對期望的Δλsample而言,ΔΦλsample較佳地隨A減小而正大。
公式(16)示出應用外裝濾波器與把濾波器直接放在檢測器陣列上的優點利用外裝濾波器,用戶只要改變透鏡孔徑大小或外裝濾波器與透鏡的距離,就能選擇期望的採樣速率。若應用對譜內容要求高解析度,可將濾波器置於遠離投射中心,而且孔徑很小。另一方面,若用戶要快掃景象,須對指定的視場攝較少的幀,則ΔΦsample應較大,因而孔徑要大,而且/或者濾波器要靠近透鏡。
另一例研究要求360°全景景觀的應用,並假定使用的成像器包括一長度為L能覆蓋成像器整個視場的濾波器。成像器長標為Ld,則F/A=Ld/L,因此sample12f#LdL---(17)]]>為捕獲360°全景,要求圖像數為 須指出,上述分析雖然假設了單一透鏡的攝像系統、小角度和作為波長函數呈線性的外部帶通濾波器遮蔽,但是上述原理通常可導出其他系統的類似關係。再者,對單透鏡成像器或任何其他系統而言,採樣間隔可以不同,並且實際上可以是非周期性的。
此外,為測定景象發射幅射(如光)的偏振特性,可分析由具有空間變化偏振靈敏度特性的成像器所捕獲的一組圖象。來自景點的光的偏振態有四種自由度,即4個Stokes參數,已為眾所周知。因此,若檢測器對景點發射光作四次或更多次測量,每次測量有不同的偏振濾波,就能評估光的偏振態。在眾多場合中,可略去橢圓與圓偏振態,只剩三種自由度強度、部分偏振和偏振平面(即偏振角)。將線性偏振器調諧成通過與光軸成α角的偏振分量,使入射光偏振平面成θ角,設入射光強度為I,其部分偏振為P,則通過濾波器的光強為g=C+Acos[2(α-θ)], (19)式中C=I/2,A=PC。此時,三次測量就足以評估光的偏振態。例如,若每次測量(對應於幀P)僅是偏振器的αp角變化,則1cos(2ap)sin(2ap)CAcAs=gp,---(20)]]>式中Ac=Acos(2θ),As=Asin(2θ)。因而M=CAcAs=g,---(21)]]>其中M=1cos(2a1)sin(2a1)1cos(2a2)sin(2a2)1cos(2am)sin(2am),---(22)]]>而m是測量景點的幀數,而且g=g1g2gm---(23)]]>略去橢圓偏振分量,用m=3方程將公式(21)變為CAcAs=M-1g---(24)]]>公式(24)可計算光強值I、部分偏振P和偏振平面角θ例如,若成像器應用覆蓋視場ΦFOV的濾波器,而且偏振濾波器跨越視場ΦFOV而變化,則在每幀之間,成像器最好將其觀察方向改變Δφsample≈φFOV/3。
須指出,在某些情況下,對來自各景點的偏振光作兩次測量,足以得到該景象的偏振信息,若有偏振的先驗信息,則更是如此。這類先驗信息包括景象中出現的典型景物的偏振特性信息。
被檢輻射的偏振角不一定是幀間變化的唯一靈敏度特性,例如輻射總衰減(影響I)也能變化。另外,成像器使用的偏振器不一定完全偏振,可以只偏重一個偏振分量。此時,分析相似,但矩陣M每行必須歸一化,以補償該行相關的衰減。在一些常規系統中,測量偏振時用偏振分束器分裂光束,並將分離的分量送給分離的傳感器。其他系統應用了旋轉線性偏振器或液晶單元一類的電控偏振設備。但根據本發明,偏振測量可用具有空間變化偏振靈敏度特性的成像器如配備空間變化偏振濾波器的攝像機執行。這種成像器作偏振測量時,不用附加傳感器,成像器的任一內部部件在圖象採集時無須移動,也不用電子濾波器。免用附加移動的內部部件,減小了能耗,提高了可靠性。
在實施本發明的拼圖法,應用的成像器712最好包括一種圖7和圖8所示攝像機702外裝的空間不均一濾波器704,因為這種結構便於成像器712通過更改或置換外裝濾波器704而修改。帶外裝濾波器的其他成像器實例示於圖9和10。圖9的成像器712包括彎曲的外裝空間變化濾波器902,其每一部分與孔徑710的投影中心O的距離幾乎一樣。在圖9系統中,濾波器902投射到孔徑710上的任意兩個同尺寸部分,如部分904與906,幾乎有同一尺寸,因而成像器712的靈敏度特性主要取決於形成濾波器902的材料的特性和/或厚度。
圖10系統包括經定向的濾波器1002,濾波器1002與孔徑710投射中心O的距離隨視角而變。尤其是濾波器1002的第一部分1004比第二部分1006更接近投射中心O。如圖10所示,部分1006定向成與通過部分1006和投射中心O的主射線808成一角度1010,部分1004定向成與通過部分1004和投射中心O的主射線804成一不同角度1008。角1010小於角1008,因而儘管部分1006的面積大於部分1004,但是這兩個部分1004與1006在孔徑710的投射中心O上具有同樣的投射面積。
在衰減濾減器的情況下,濾波器特定部分的有效衰減量隨光通過濾波器的角的銳度而增大。在幹擾濾波器等波選濾波器的情況中,濾波器發射的波長取決於入射角,因而外裝濾波器尤其有利,因為若濾波器是柔性或可拆卸安裝的,就容易修正其角定向,所以變於修正成像器的特性。
雖然以上強調了應用帶外裝濾波器的攝像機,但是空間變化濾波器不一定外裝於攝像機,如圖11示出的成像器712將其濾波器1102直接裝在檢測器708上,其優點在於,由於檢測器708處於透鏡802的聚焦平面,而且濾波器1102極靠近檢測器708,因而濾波器1102的特性模糊極少。結果,成像器712的特性由濾波器1102更精確地限定,從而可應用濾波特性變化很微細的濾波器。然而,這一優點的代價是更改濾波器的靈活性喪失了,一般與帶內部濾波器而不用外裝濾波器的攝像機有關。
圖12示出的攝像機中,濾波器1202位於攝像機內部,但不與傳感器708直接接觸。圖12的結構可使用彎曲或傾斜的濾波器。
在圖13所示帶複合透鏡的攝像系統中,空間變化濾波器1302可裝在複合透鏡的元件1306與1304之間。在圖13實例中,物鏡1306形成目標圖象。另一透鏡1304將圖象透射到檢測器708上,濾波器1302較佳地靠近透鏡1304。所有攝像機透鏡幾乎都是用於校正像差等的複合型,通常有兩個以上元件。
如圖14所示,帶複合透鏡的成像器還可包括位於物鏡1404聚焦平面的擴散器1406。具有空間變化特性的濾波器1402儘量靠近擴散器1406。物鏡1404在擴散器1406上形成由空間變化濾波器1402修正的聚焦圖象。另一透鏡1408將聚焦圖象從擴散器1406投射到檢測器708。
本發明的成像器可以配置成利用有空間變化反射比的反射器接收景象的輻射(如光)線。在這種系統中,圖15的示例用具有空間變化反射比的反射器1502將景象的光信號810、812與814(即從點A、B、C)反射入攝像機孔徑710。若空間變化反射比反射器1502部分透射(即發射某種光),它可以置於攝像機孔徑710與參照目標1506之間,而參照目標1506可校正系統或吸收攝像機702的雜散反射,例如已知的校正技術可利用已知光源校正攝像機中檢測器的特性。圖15的結構在功能上相當於虛擬攝像機1504,後者與實際攝像機702的定向角度不一,能直接有效地接收光線束810、812與814,不靠反射。
在成像器712中使用外裝反射器1502的另一優點是反射器1502可任意彎曲,便於控制成像器712的視場和/或放大倍數。
圖16示出的成像器712,有一位於攝像機702外部的濾波器1602和內部的另一濾波器1604。
圖17的成像器712具有在攝像機702外部的濾波器1702與參照目標1780,另一濾波器1704在攝像機702的光學元件1710內,還有一個濾波器1706在攝像機702內部靠近檢測器708。
本發明的成像器,如圖17的成像器712可以平移而掃描景象706,如圖18所示。隨著成像器712沿圖示方向平移,景象706的任一特定點,如圖18的A點,可在成像器712的參考幀裡沿多條主射線804、806與808成像。另如圖19所示,成像器712不一定直線移動,而是能沿路徑1902運行,而且還可以同時平移和轉動。此外,儘管圖7~18的濾波器主要圖示成從上到下變化,但是本發明的濾波器704能沿任何方向變化,例如濾波器704的特性能水平變化,此時拼圖法包括使成像器712繞其光軸轉動以捕獲圖象作處理。再者,雖然濾波器704能以不彎曲和/或不移動的方式裝到攝像機上,但是也可使用相對攝像機702彎曲或移動的濾波器。在這種系統中,濾波器可在一系統圖象的各幀之間彎曲或移動。
如圖20A與20B所示,本發明的拼圖技術能擴展圖象測量的空間範圍(即視場的寬度或高度)和總強度範圍(和相應的動態範圍)。圖20A示出單幀2002與拼圖2004的空間範圍的差別,拼圖2004由多個幀構成。圖20B示出拼圖2004擴展的空間範圍與強度範圍。實際上,由於拼圖2004邊緣附近的像素被採樣較少次數,所以其總強度範圍比拼圖2004中心像素的強度範圍小。該效應可從圖21看出,在靠近拼圖2004空間範圍邊緣的區域2102中,總強度範圍縮小了。但在遠離邊緣的部分2104內,強度範圍均一地更大,這一效應類似於人肉眼的凹型,視場中心附近的解析度更大,這有利於景象中心部分比邊緣位置更受關注的場合,因而成像時比邊緣部分具有更高的解析度。
如上所述,由於外裝濾波器相對於攝像機檢測器容易偏離焦點,所以它對攝像機特性的作用容易變模糊,如在圖22A所示的濾波器密度函數2202中,密度D在第一位置2214與第二位置2216之間跨過濾波器表面一維呈線性變化,但在該範圍外不變。透射比T為10-D,因而所得成像器相應的強度靈敏度特性函數在整個濾波器上呈指數變化。由於濾波器密度函數2202的模糊作用,濾波器的有效密度函數2204有圓角。同樣如圖22B所示,有階躍狀密度函數2206的濾波器具有帶圓角的有效密度函數2208。如圖22C所示攝像機外裝濾波器可以具有任一種任選的密度函數2210,但若函數2210有銳角,則有效函數2212就有圓角。
對於空間變化衰減遮蔽函數f,一般可將有效遮蔽函數M模擬為M(x,y)=f(x,y)*h(x,y) (25)式中h(x,y)是系統對遠景聚焦時,攝像機對濾波器那樣近的目標的散焦模糊點擴展函數(PSF)。對於圓對稱PSF,遮蔽有效地是x的一維函數,f與 (即h的Abbel變換)卷積。例如,若把核模擬為高斯函數,則h是帶標準偏差r的高斯函數,M(x)=erf(-x/r)。若從適當位置通過濾波器觀看,由於M(x)取0與1之間的任一值,則任一景點在理論上都可成像而不飽和,而與該點亮度無關。因此,該簡型系統在理論上有無限大的動態範圍。
以圖23的簡型遮斷器2302為例,可進一步說明外裝附件的模糊作用。雖然遮斷器2302本身具有明顯的急劇空間變化特性,但它對攝像機702特性的作用在攝像機702視場內是漸變的,如圖24A與24B所示。圖24A示出簡型遮斷器2302更圓滑的的有效密度函數2402,圖24B示出簡型遮斷器2302的有效密度函數的對數2406。此外,圖24A還示出了對線性變化密度濾波器的有效密度函數2404的模糊作用,圖24B示出線性變化密度濾波器有效密度函數2404的對數2408,對數函數2408是直線。
上述易於在帶外裝光學元件的成像器中出現的模糊作用,對成像器的可調性很有利,因為只要沿成像器的光軸改變該光學元件的位置,就可更改成像器的光學特性。例如,簡型遮斷器2302有效密度函數2402中間附近的斜率(見圖24A與23),通過將遮斷器2302移得更接近攝像機702的焦距而增大,而通過移得更遠離來使該焦距減小。
儘管以上強調了導致具有空間變化、波長中性強度靈敏度特性的成像器的空間變化、波長中性密度濾波器,但是也可應用空間變化濾色器或其他基于波長的濾波器。例如,這種空間變化波長濾波器可以是一種空間變化中心波長為~的空間變化帶通濾波器(如幹擾濾波器)。圖25A示出濾波器上中心波長λ0與位置x的函數2502,該濾波器的中心波長λ0在視場內線性變化。濾波器函數2502產生一種具有波長靈敏度特性函數2504的成像器,函數2504在視場內也近似線性地變化,但有圓角。圖25B的一例濾波器,其中心波長特性函數2506在視場內呈階躍變化。這種濾波器導致的成像器,其波長靈敏度特性函數2508基本上對應於濾波器的階躍狀函數2506,不過因模糊作用而有圓角,因為該濾波器對攝像機有點散焦。如圖25C所示,可以使用具有任一隨選中心波長函數2510的帶通濾波器,但若濾波器的函數2510有銳角,得出的成像器將具有帶更多圓角的波長靈敏度特性函數2512。
本發明的成像系統可應用空間變化低通或高通濾波器。例如,圖26A~26C示出了低通或高通濾波器上的截止波長λcut與位置x的示例函數2602、2606與2610。與圖22A~22C和25A~25C所示的濾波器相似,模糊作用使得用各自濾波器構成的諸成像器各自的波長靈敏度特性函數2604、2608與2612圓化了。
圖27進一步示出了模糊作用,圖示為檢測器靈敏度與光線束波長λ的曲線圖,該光線束的主射線靠近高通濾波器陣列兩段的邊界。一半束光通過截止波長為 的段,另一半束光通過截止波長為 的段,檢測器的有效截止波長為 ,介於 與 。相對於兩種獨立段各自的過渡帶,邊界上有主射線的光束的過渡帶更寬。
與上述空間變化衰減濾波器一樣,改變濾波器沿成像器光軸的位置,也可調整波長濾波器的有效特性。具體而言,若濾波器更接近攝像機的焦距,其特性變得更銳、不圓和更斜,而隨著濾波器移離攝像機焦距,其特性就變得更圓滑而不很斜。實際上,特性2506與2606分別示於圖25B與26B的階躍狀波長濾波器陣列,若充分偏離焦點,可用於近似為特性2502與2602分別示於圖25A與26A的線性變化濾波器。
本發明的空間變化帶通濾波器還可用來提高彩色成像器對鑑別接收自景象的信號波長的精度(即動態範圍)。如圖28所示,典型的彩色成像器(如彩色攝像機)有三條相當寬的彩色通道2802、2804和2806,可分別接收藍綠紅光。對成像器接一個帶較窄通帶2808的濾波器,可提高其鑑別波長的能力。例如,若窄帶濾波器的通帶2808在攝像機較寬的藍帶2802內,則最後被攝像機檢出的藍光2810要比綠紅光2812與2814強得多。相應地,由於窄帶濾波器特性2808不像原始攝像機特性2802、2804與2806那麼寬,測定的色彩比僅將信號定為「藍」更仔細。窄帶濾波器通帶2802的中心波長λ在視場上變化,因而成像器跨越景象作隨動拍攝並捕獲多幅圖象,各景點在多條色帶內成像,各有一不同的中心波長λ0,因此能以更大的色解析度測定各景點發射光的譜特性。
圖29示出被分成三大部分的外接濾色器的用法,藍色濾波器部分具有藍色帶通特性2904,綠色濾波器部分具有綠色帶通特性2906,紅色濾波器部分具有紅色帶通特性2908。三色濾波器陣列外裝到彩色攝像機,攝像機的一組三色通道分別由通到攝像機中檢測器的藍綠紅光的帶通特性2902、2912與2910限定。當成像器跨越景象隨動拍攝時,各景部通過外裝的藍綠紅色濾波器部分觀察並成像,因而各濾波器部分的各藍綠紅特性2904、2906和2908就同攝像機本身的各藍綠紅特性2902、2912和2910一起應用,所以該成像器實際上有9條色帶,代表了傳感器通道與外裝濾波器元件特性所有可能的組合。相應地,攝像機的色解析度提高了。一般而言,若攝像機本身有Ns條帶,而附加的遮蔽有Nm條帶,則成像器整體能測量NsNm條帶。
根據本發明的另一個方面,具有空間變化聚焦特性的成像器可以跨越景象隨動拍攝,而且可根據各景點在成像器712隨動拍攝時捕獲的成組圖象中每幅圖象內的聚焦情況,來確定諸景點與成像器712的距離。圖31示出成像器712對分別由景象聚焦部3104與散焦部3106發射的光線束3108與3110的響應特性。物鏡3102與攝像機702的焦平面3118的距離為 。聚焦景點3104作為單點3114投射到焦平面3118上,但散焦景點3106以直徑為d的圓盤3116投射到焦平面3118上。
一般而言,與攝像機702的孔徑710的距離為 的目標處於聚焦狀態,相應地,處於聚焦的任一目標的距離均為 。與之相對照,任何散焦的目標都與孔徑710有不同距離(如u)。
成像器經配置,可使目標聚焦的距離 在成像器視場內變化。若聚焦特性變化足夠大,且可作足夠多的快照,則在至少這些快照之一中,所有或大多數景點將處於或接近聚焦狀態。測定成像器沿聚焦的目標接收的特定主射線並掌握了成像器的空間變化聚焦特性後,可按下法測定各景點的距離。
任何光學系統都有一個目標聚焦或接近聚焦的距離範圍,目標移離該範圍就會變模糊。該範圍的大小通常稱為系統的「景深」。本發明的成像系統把這一距離/深度範圍轉換成成像器內沿其光軸的位置範圍。該位置範圍的大小稱為「焦深」。在許多場合中,景物所處的深度範圍寬於景深,因而某些目標會散焦。參照圖31的成像系統,目標點在距離u的模糊直徑d為d=D|uF-v~u+Fv|~Fu,---(26)]]>式中F為焦距,D為透鏡孔徑。能以不同的聚焦設定捕獲多幅目標圖象(如採集期間移動系統內部部件,諸如透鏡或檢測器)。選擇在對應於目標點的像點使圖像銳度最大的圖象,通常可測出該目標點的聚焦狀態。再者,可用「融合」法對整個景象構成完全清晰的圖象。這種方法對每個景區或景點在景區或景點最清晰的圖象中選擇相應的點,再把整組選出的點組合成所有或大部分景點的清晰圖象。
還可用其他技術測量景點的距離/深度。例如,若改變 而改變聚焦設定,並用圖象拉普拉斯算子測出銳度,就可將最佳聚焦地v算為v=argmax|x,y2Iv~(x,y)|---(27)]]>式中 指圖象平面中的拉普拉斯算子,而各圖象I(x,y)用 參數化,即定義為 的函數。要注意,一旦確定了最佳聚焦平面,就可以用它來推導目標點的深度u。另外,該深度不僅可通過求出最佳聚焦平面推導出來,還可通過基於兩個或更多幀運用「散焦深度」法估算模糊直徑d而求出。2001年5月8日頒發給Nayar等人的題為「Apparatus and Methods for Determining theThree-Dimension shape of an Object Using Active Illumination andRelative Blurring in Two-Images Due to Defocus」的美國專利N0.6,229,913,已詳述過該「散焦深度」法,其內容通過引用包括在這裡。
圖象中某一景點由尺寸對應於不確定位置U(即模糊量)的區域代表,而U由檢測器陣列的相鄰元件間的像差、衍射模糊或有限距離ΔXgrid造成。若模糊直徑d不大於U,可認為該景點在系統的景深內,因為認為聚焦了。由於模糊直逕取決於與最佳聚焦態的距離,所以d中的U變化對應於特定的深度(即距離)變化。例如研究這樣一種情況,即為在特定的景深間隔內捕獲一組圖象,要將成像系統某個元件移動Δv。不用物理地改變Δv,可用一塊折射率為n的透明介質改變光程長度。為作調整,要改變介質的厚度t或折射率n,實際上這兩種方法可一起應用。根據視場內的厚度或折射率變化速率,系統計算幀間橫越增量,在選擇的景深間隔內產生期望的有效軸向變化。
將整個光學系統移近或移離景物,可避免移動部件,但這種方法只有在深度範圍極小諸如顯微鏡應用時才實用。另一方面,若把一塊玻璃或聚碳酸酯等透明介質放在檢測器上,檢測器與透鏡間的光程長度就變長。對摺射率為n的材料,沿光線的每毫米材料相當於通過n毫米自由空間的傳播,這有效地將一部分 拉長了n倍(若折射元件在攝像機內部),或拉長了一部分 (若折射元件在攝像機外部),使聚焦的目標點變散焦。反之,原來在檢測器平面後面某一平面上聚焦的某些目標點,現在變成在檢測器上聚焦了。
本發明一例空間變化折射元件是一塊透明材料,厚度在視場內變化。或者,除了空間變化厚度以外,濾波器具有空間變化折射率,這類濾波器的幾例特性示於圖32A~32C。圖32A示出折射率線性變化的稜鏡或折射元件的折射特性。若將這種元件置於光學系統之外,諸如圖8所示的成像器712之外,將基本上只偏轉視場,由此改變了系統的投射中心,很少有助於延長成像器的景深。因此,最好使該元件位於或靠近檢測器,如圖11所示的成像器712中。若在圖14所示成像器712中的擴散器1406上形成中間圖象,則最好將元件1402直接裝在擴散器前面。
若折射元件的折射特性和/或厚度作為視場內的位置的函數而變化,則得到的成像器將具有空間變化聚焦特性。在圖32B的實例中,折射元件的厚度或折射率以階躍狀變化。實際上,折射元件的厚度或折射率可按圖32C所示的任一種函數變化。
以上討論了靈敏度或聚焦特性在視場內變化的成像器的用法。如上所述,這類空間變化靈敏度或聚焦特性可以包括空間變化強度靈敏度特性(如由空間變化衰減器形成)、空間變化波長靈敏度特性(如由空間變化濾色器形成)、空間變化偏振靈敏度特性(如由空間變化偏振濾波器形成)和/或空間變化聚焦特性(如由空間變化折射元件形成)。根據本發明另一個方面,可將二種或更多的上述特性組合在單個成像器內。如圖33所示,攝像機可同時將空間變化衰減器3302、空間變化濾譜器3304、空間變化偏振濾波器3306和空間變化折射元件3308全部或部分內接或外接起來,使形成的成像器具有多種在成像器視場內變化的特性。
在有些場合,光的每種特性與其他特性無關,如在體積目標中,各像點接收來自不同深度的多個(或連續)點的光。任一點發射的光可以包括源自不同景點或景區的光(如不同景點或景區的散射或放射),因而各譜分量具有不同的偏振、亮度與聚焦態。
但在大多數情況下,某些特性會被簡併或高度耦合。如在特定偏振分量與譜帶內,來自特定深度的某種輻射可能強得無法不飽和地檢測,因而最好用寬亮度動態範圍的濾波器檢測強度幅值。此外,若目標不是三維的,則深度測量應與光譜和偏振測量分開。尤其是,為避免檢測器飽和引起的深度估算誤差,在具有可變的厚度或折射率的濾波器中應用可變密度濾波處理,有利於擴展亮度動態範圍。此時最好應用所有可能的帶寬靈敏度、偏振角靈敏度和強度靈敏度組合,在各像素位置作多次獨立的測量。此類測量可用配置成使用複合濾波器的成像器執行,該複合濾波器有各種區域,每個區域為特定靈敏度特性的空間變化專用。如圖34所示,可以使用濾波器3410,其第一段3402用於空間改變成像器的第一靈敏度特性,第二段3404用於改變第二靈敏度特性,第三段3406用於改變第三靈敏度特性,第四段3408用於改變第四靈敏度特性。此外,各段3402、3404、3406與3408不一定限為具有任何特定的形狀或位置,或沿任一特定方向具有各自空間變化的濾波器特性,而可以是任何形狀和/或位置並在濾波器3410上沿任何方向空間變化的濾波器特性的區域。
須指出,各濾波器特性的變化可配置成對濾波器有一不同的空間頻率,例如深度隨cosx而變,中性密度隨cos2x而變,中心波長隨cos4x而變。而且,各濾波器特性的空間變化不一定是正弦形,也可採用方波形、鋸齒形和其他圖形。
另一方面,有些濾波器特性能不要求專用的獨立濾波區域,如來自某景點或景區的光的偏振往往將比較強烈地依賴於光波長。相應地,某個空間變化偏振濾波器重疊一空間變化波長濾波器,很少或不會丟失信息。再者,來自景象所有部分的光不一定按其特性或以同一解析度來測量。此外,對於速度遠比測量完整度重要的場合,可用圖33的濾波器3302、3304、3306與3308等成組重疊濾波器作快速掃描,不使用圖34的複合濾波器3410;複合濾波器3410一般要求較慢的掃描來得到給定的解析度。
用本發明的掃描法獲取圖象後,在相應地像素位置即代表同一景點的像素位置獲得的數據經處理,可提取被處理圖象各像素的高動態範圍強度、偏振態、深度與色彩的值。為了融合各景點得到的所有原始數據,要配準一系列圖象。
圖35A示出一系列快照3502、3504、3506、3508和3510,其中在各幅快照的像素位置3522、3524、3526、3528和3530檢測某景點發出的光。測定隨便哪個代表各快照該景點的像素的位置,就能以數字方法對準圖35B所示的快照3502、3504、3506、3508與3510,並處理相應於各像素位置3522、3524、3526、3528與3530的數據,可產生增強動態範圍、譜信息、偏振信息和/或深度信息的優質像素,如上面針對本發明各種拼圖技術所討論的那樣。配準方法包括下列方法的一種或多種1.可以知道幀間的移動。將成像系統裝在電機上並按連續圖象間觀察方向的偏移校正其速度,就是這種情況。從位置受監視的衛星、飛機或其他運載工具攝取圖象也是如此。
2.人工配準圖象,如使圖象的控制斑點匹配。這些斑點一般包含明顯的圖象特徵,其形態耐受來自這類圖象特徵的光的特性變化。換言之,在不同的波長帶、偏振態和寬範圍的光強分布中,控制斑點的形狀相似。
3.跟蹤第二方法所述的控制斑點,自動配準圖象。這類配準法已為眾所周知,有些方法應用較高程度的特徵描述,但大多數方法使斑點間的相關性最大,或使斑點差最小。若圖象在幀間的變化很小,這類方法效果更佳,否則無法作相似性測量。例如,景象在λ=630nm的內容可同在λ=400nm的內容完全不相關。然而,若採用差別不大的波長捕獲圖象,而且捕獲圖象的光的波長帶足夠寬,則連續圖象將一般具有使它們匹配的充分的相似性。
4.求出能提供最佳整體匹配的圖象間的變形參數而自動配準圖象。這類方法已為眾所周知,常用匹配指標包括MSE判據的質量變化、耐統計誤差指標和相關性。這類方法與方法3相似,即更適合圖象在幀間變化小的場合。這些方法對圖象可作有效的粗細比較與配準。
5.當圖象在幀間變化很大時,匹配問題類似於不同傳感器獲取配準圖象的情況。雖然景象在幀間通過成像器視場移動,但是遮蔽在特徵上控制著圖象,由此會使這種算法配準小於正確平移的圖象平移,或甚至配準不平移。眾所周知,解決該問題的一個辦法是結合應用方法3和4。在這種組合方法中,為增強對在圖象中不容易變化的明顯特徵的依賴性,各圖象要作高通濾波。應用粗細法配準,並優化像間整體變形參數。
6.即使變化很大,如在多傳感器系統中,配準圖象的另一方法是用交互信息,即圖象間特徵的相似性作為匹配標準。該法一般要求圖象在統計上並非完全獨立,這是一種有效假設。
有時為了便於匹配,可以部分補償濾波作用,例如假定成像器應用一種變化的密度濾波器。這種濾波器起遮蔽作用,將視場各部分衰減M倍,M在視場內變化。由於遮蔽此時為相乘,有利於運算各圖象的對數,使遮蔽函數M成為相加分量,於是對得到的對數圖象作高濾波可減小遮蔽作用,因為遮蔽變化通常較慢。若在配準操作開始之前已得到M的估值,則在操作前,可用相應的1/M值放大各圖象的每一像素值,該步驟可減輕短暫恆定遮蔽對運動估值的偏置作用。因強衰減而變得更暗的測量結果,經傳感器量化後將變得較嘈雜。在比較與配準圖象時,最好計入這種與衰減相關的不確定性。
還要注意,若重疊區比各幅圖象的像素格更密集地採樣,能以子像素精度作配準,因此除了得出景點的多維特性以外,還可同時獲得空間超解析度。
一示例算法以遮蔽M的粗估值為基礎,不確定性為ΔM。估值1/M用於使圖象場展平。如上所述,由於映射變換,衰減成極暗的像素(尤其是量化為零的像素)變嘈雜了,因而在後面的所有算法中,測量結果在作比較、區分與融合時,最好考慮到不確定性。該法的原理如下1、粗平整場反應1/M。設落在檢測器上的光強在M=1時為I,測量單位為 ,無量綱。落在檢測器上的的光在光學濾波後為g(x,y)=M(x,y)I(x,y)。對每一測量結果g±Δg,估算出景象幅照度I±ΔI。運用一次近似法傳播不確定性。
2、由於應用了高動態範圍圖象,所以把暗景點誤差為像亮景點誤差一樣重要,例如應將亮點的5%誤差像暗點的同等誤差那樣給予補償。為調節各像素的亮度,要算出各被測景象幅照度的對數,即S=logI,得出S±ΔS。
3、粗/細樣式。眾所周知,可靠有效的配準方法是用圖象四面體(多解析度)表示法。在以粗空間分辨度對圖象配準和彎曲後,該粗移動估值就用作以更細的標度改善轉換估值等的初始條件。但與忽略了數據不確定性的常規方法相比,本發明方法對S圖象建立一個最大似然四面體;在四面體的每一層面,代表某一鄰域的值是在研究了參與像素的不確定性後最可能的值。
4、配準似然性的最大化。估計測量值S具有高斯分布(雖然方法不限於高斯分布)。最佳配準是將圖象間的Mahalanobis距離減至最小的一種配準。
5、拼合配準。在球坐標系中,只配準成對圖象會導致估算圖象位置誤差累加,因此最好對最新更新的拼合配準系列中每幅新圖象,再與更新拼圖融合。
該算法詳述如下。檢測器在其動態範圍內有一響應R(可以為非線性),像素的像值 為g~=R(g)---(28)]]>使響應線性化,檢測器的估算強度 為g^=R-1(g^),---(29)]]>因而g^=|dR-1dg~|g~---(30)]]>
假定估算的遮蔽與測得的信號無關,(I)2=(g^/M)2+(g^M/M2)2,---(31)]]>其中假設g~=0.5]]>,因為攝像機輸出為整數形式(8位檢測器為0~255)。為避免暗點(I=0)對數運算潛在的不穩定性,後續運算可以除去諸暗點。因此,計算中最好應用S=log(1+I),不用I。在任何情況下,s=|dsdI|I---(32)]]>注意,若檢測器響應R為對數,不必應用公式(29)~(32),S可置成 。任何被認為飽和的像素(如對8位檢測器而言, 接近255),被作為極不確定對待。因此,將其相應的ΔS置成極大值(即偏離系統的動態範圍)。
假定S測量結果為高斯型,因而對兩個獨立的測量結果S1±ΔS1,和S2±ΔS2而言,Se值的對數似然特性如下logLikelihood~-E2(se-s1s1)2-(se-s2s2)2---(33)]]>Se的最大似然解使Mahalanobis距離E最小se=s2^(s12s1+s22s2),---(34)]]>其中s2^=(0.5*d2E2dse2)-1=(1s12+1s22)---(35)]]>於是,對應該景點的圖象測量的距離值為E^2=(s^-s1s1)2+(s^-s2s2)2---(36)]]>假定圖象的所有像素均獨立,則整個幀或像素任何其他子組的距離量度 為E^total2=allpixelsE^2(eachpixel).---(37)]]>根據上述目標函數,兩幀間(或新幀與原有拼圖間)的最佳配準是使 最小的配準。這裡每對測量S1與S2對應於相應像素的像值。
注意 一般隨像素數而增大,使配準偏向將公式(37)的和中的像素總數減至最少,減少了圖象間的重疊。為對付這一作用,可用重疊區裡的像素數,即 或 等使公式(37)歸一化,其中stotal2(allpixelss^eachpixel-2)-1---(38)]]>若不同測量或不同像素間的統計相依性不能忽略,則 與 的公式可以統一使用測量結果的協方差矩陣,而不單是它們的方差。
為使配準更加可靠而有效,要從粗到細解析度分層地配準。圖象的特定四面體等級的粗表示,可通過對圖象以某個核作低通濾波後再欠採樣而實現,核的寬度取決於該等級(等級更高/更粗,對原始圖象操作的核越寬)。在任何情況下,這種表示法的像素值是被測像素的加權歸一化之和s=kkskkk---(39)]]>式中ωk是像素值Sk的權重。
須指出,這裡的討論指四面體等級的原始全解析度圖象的結構,諸像素是獨立的,以簡化推導。但四面體通常是迭代構成的,中間等級的像素在統計上不獨立。若在迭代過程中尋求附加精度,則加權不僅取決於中間等級的像素方差,還依賴於其鄰居的全協方差矩陣。因此,該矩陣也要擴展到四面體。
在常規四面體中,權重ωk等於高斯核值ak。但在本發明方法中,可將公式(39)視作公式(34)的統一化。指定給像素的權重ωk隨其高斯權重ak的減小和其不確定性ΔSk的增大而線性地減小,因而k=ak/sk2]]>,而且s=(kk)-1---(40)]]>在上述最大似然四面體中,各分辨級的表示不僅包括各像素的加權平均值,還包括各像素值的不確定性。在該四面體中,由於不確定性較小的點接收更大的加權,因而不確定區(如飽和區或極暗點)在較粗等級時減小了重要性,其表示受到相鄰穩定點更多的影響。因此由一個值表示的區域更可靠。相應地,以每個標度代表S±ΔS,可通過使匹配似然性最大而有效地配準圖象。
須指出,上述方法並非是配準與融合圖象系列的唯一方法,例如可用其他加權法測定各點的值,或用每點合格的單次測量結果代表該點。若運動場被參數化,中性密度遮蔽不隨時間變化,而且景象靜止,則可將遮蔽、運動和未知景象的所有參數當作單一的優化問題。
成組圖象被捕獲配準而且知道了成像器特性(如通過校正)後,本發明的算法產生至少一個數據(如像素)代表景象的每一點或每個小區域。這種融合法用來自每幅原始圖象的該特定像素產生代表性像素,而該特定像素對應於被代表的景點。融合法的選擇一般取決於要成像的景象特性,如亮度、色彩、偏振或深度。
例如,為融合以空間變化衰減法獲取的強度(即亮度)圖像,最好應用下列原則1、該點所選的值要容易將到數據點的Mahalanobis距離減至最小,因而公式(34)和(35)使用sk=Ik,其中Ik=g^k/M]]>。
2、對原始圖象的飽和點(在8位攝像機中, 接近255)指定極大的不確定性,如上所述。
3、用接縫最小化方法避免邊界可能出現的無審美性接縫。拼合圖象邊界上,存在用不同數據源評估的點之間的過渡區。除去拼合接縫有多種方法。一種方法是尋找「優化」接縫線,即對觀眾最不顯眼的接縫線。優化接縫線選擇法已在本領域廣泛採用。另一方法是應用羽化操作法,對圖象按相對於圖象中心或邊界的像素位置來加權。這種加權法容易配合方法1所述的加權平均法,且類似於上述在四面體結構中使用的加權。例如,可將不確定性乘上某個因子,該因子在圖象邊界平滑地增大到∞。
4、在不確定性有突變而 變化通常小的飽和區邊界,也會出現接縫。這種接縫可用方法3討論的方法消除。最好用下述方法平滑飽和區的清晰度。規定一「低飽和值」a(如8位攝像機為230),再規定一「高飽和值」b(如250),若 ,認為該像素飽和了;若 ,而且該點鄰近某飽和點,則認為該點「與飽和相關」;若 ,而且該點鄰近「與飽和相關」點,也被認為「與飽和相關」。這樣,「與飽和相關」的點總被包括在與飽和點連接的組內,因此該方法不影響與飽和點不相關的強度為 的點 。找出了被稱為「飽和相關」的所有點之後,將其不確定性乘上 或另一函數,隨著到達飽和值,逐漸從規則不確定性(乘上1)過渡到極大不確定性。
須指出,本發明的數據融合併不限於上述示例方法,還可運用減少其他效應的持久統計法或運用迭代投射法融合數據。還能用多解析度方法,即通過以若干不同解析度分析圖象來消除接縫,這在本領域是公知的。
還要指出,上述原則3和4對高動態範圍拼圖的產生並非關鍵,只是用於使得出的圖象更具審美性。其他方法包括(1)只選各點單次「合格」的測量結果作為該點的代表值;和(2)選擇產生最大反差的值(如以各種標度[說明])。
對偏振數據,圖象點分析如下。假定可以忽略橢圓偏振分量,並對各點應用三次測量,就可用公式(24)導出各像素的偏振態分量。若應用多於最小要求的像素,則可應用一組超定的方程。通過將方差減至最小求解該組方程,得出CAcAs=(MM)-1Mg~---(41)]]>然後,該方法設I=2C,A=Ac2+As2,P=A/C]]>,而θ=Sin-1(As/A)=Cos-1(Ac/A)。
應當指出,上述方法並非是能用數據提取偏振信息的唯一方法,比如可用持久統計法估算偏振。P、θ與I可以直接估算,即略去中間變量Ac、As與C;可以應用投射迭代法(如迫使I和P為正);或者不必提取偏振信息,只是通過比較諸圖象的拉普拉斯算子四面體表示法並選擇具有最大值的表示法,以產生最大反差的方式來融合諸圖象,這是本領域公知的。
為導出深度信息,要對相應的點實施由焦點得深度的方法,選擇給出最佳焦點的圖象,像常規的由焦點得深度技術一樣。接著融合選出的圖象,只讓「清晰的」像素作拼合。例如,若拉普拉斯算子測量圖象清晰度,則I^(x,y)=Ip(x,y)(x,y)]]>p(x,y)=argmax|2x,yIp~(x,y)|,---(42)]]>式中(x,y)是靜止(外部世界)坐標系裡的坐標, 為幀索引。然而,融合不一定只在像素一級,能以多個解析度標度實施。再者,該技術並不限於聚焦的拉普拉斯算子指標,也可用其他聚焦指標。
為導出光譜,對不同譜帶的同一景點連續作相應的測量。若濾波器不是帶通濾波器,而是高通或低通濾波器,可利用相鄰測量的微分導出窄帶譜信息。光譜可用變寬帶測量,甚至包括多條帶。考慮到帶結構、譜的正性和其他約束條件,求解一組原始數據提供的方程,可提取高解析度光譜。還可利用已知源與材料光譜的先驗信息,確定哪些類型的已知源與照射分量符合測量結果,包括相對強度、溫度和/或測量的空間位置。
應用這裡描述的本發明方法,可以擴展任何視頻攝像機(或靜像攝像機)的動態範圍,顯著提高攝像機的譜解析度,而且/或者將濾波器接在攝像機透鏡前面,或者甚至拓展攝像機原有的漸暈特性,可直接求出進入信號的偏振態與聚焦特性。
景象譜內容知識可以導出拍攝時景象照明光源的信息,也可導出目標反射比,因此可用不同模擬的照明源與特性拍出景象。例如,可在白熾燈下以模擬方式拍出被攝目標,就像在螢光燈或夕陽下拍攝時一樣。
還要指出,大多數檢測器的光譜響應特性,即使在檢測器中使用了RGB濾波器,也與人肉眼的譜特性極不一樣。膠捲與彩印也有同樣的限制。掌握了景象實際的譜內容,就能以更接近肉眼直接觀察景象的方式拍攝圖象。還可根據膠捲、屏幕或印表機的特性提供圖象,而不按照肉眼的響應特性。
本發明的技術在拼圖的同時擴展了成像器的光學動態範圍,動態範圍對各景點擴展,與周圍諸點的亮度無關。成像器不必內裝任何移動部件,成像器的移動與建立普通拼圖的要求一樣。本發明的技術還可結合其他擴展動態範圍的方法,如有源CMOS檢測器或AGC,實現進一步擴展。注意,若密度濾波器有較寬的密度範圍,光學動態範圍就更大。例如,若M∈[10-4,1],則檢測器的動態範圍將擴展得超過其原範圍約13位(相當於提高80dB)。由於濾波器密度是相加特性,所以串接若干這種濾波器可加強空間變化衰減作用。
增大焦深可在寬距離範圍內產生清晰的圖象。採集偏振數據可消除或增強反射與半反射視覺效果。
本領域的技術人員知道,圖1~6、35A與35B所示的方法可在附錄中程序說明的合適軟體控制下運行的各種標準計算機平臺上實施。在有些場合,專用計算機硬體,如留駐在標準個人計算機或工作站的總線上的外設卡,都可提高上述諸方法的操作效率。
圖37和38是適合實施本發明的示例計算機硬體。參照圖37,計算機系統包括處理部3710、顯示器3720、鍵盤3730和滑鼠,還可包括其他輸入裝置,如掃描圖象媒體3700的光學掃描器3750和攝像機3780,此外還可包括印表機3760。該計算機系統一般包括一臺或多臺碟片驅動器3770,可對磁性媒體(如盒帶)或光學媒體(即CD-ROM)等計算機可讀媒體進行讀寫,存貯數據與應用軟體。
圖38是處理器部3710的功能框圖。處理器部3710一般包括處理單元3810、控制邏輯3820與存儲單元3830以及硬碟驅動器與接口,還包括定時器3850(即時鐘電路)與輸入/輸出口3840。根據處理單元使用的微處理器,處理部3710還可包括協同處理器3860。控制邏輯3820與處理單元3810一起提供存儲單元3830與輸入/輸出口3840通信所需的控制。定時器3850對處理單元3810和控制邏輯3820提供時序參考信號。協同處理器3860可增強實時複數運算的能力。現代化計算機比協同處理器具有更多的處理單元。
存儲單元3830包括不同類型的存儲器,如易失與非易失存儲器和只讀與可編程存儲器。如圖38所示,存儲單元3830包括只讀存儲器(ROM)3831、電可擦可編程只讀存儲器(EEPROM)3832和隨機存取存儲器(RAM)3833。可用不同的計算機處理器、存儲器結構、數據結構等實施本發明,而且本發明並不限於特定的平臺。例如,雖然圖37把處理器部3710示為計算機系統的一部分,但是也可將處理部3710和/或其圖示元件包括在靜像攝像機或運動圖象攝像機(如視頻攝像機)之類的成像器中。
本領域技術人員知道,附錄原始碼表示出的軟體可用各種程式語言編制。本發明的示例軟體算法都用MatlabTM語言編寫。附錄提供了若干此類示例算法的計算機原始碼。
雖然已結合諸特定示例實施例描述了本發明,但是應該明白,在不違背所附權項提出的本發明的精神與範圍的情況下,可對揭示的諸實施例作出各種變化、替代與更改。
附錄用於拼合、平整場、圖象融合、高動態範圍、多解析度配準、最大似然四面體、衰減遮蔽自校正、接縫與飽和區羽化等的MATLAB原始碼。形成高動態範圍拼合的原始碼<![CDATA[function m=calchorizprof2(basefilename,postfix,netframes,left,right,bottom,roof);imagehight=roof-bottom+1;basephotoname=′photo′;rawname=′basicim′;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALCULATING THE MEAN PROFILEystep=-1;%frame_inds=36ystep1;frame_inds=length(netframes)ystep1;imagewidth=right-left+1;M=zeros(length(frame_inds),imagewidth);for indframe=1length(frame_inds), yframe_ind=frame_inds(indframe);   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 postfix];   rawim=imread(shem1);   basicim=double(flipud(rawim))-1;   X=basicim(bottomroof,leftright); %LX=log2(X+eps); %saturated=(X>sun); %unsaturated=(~saturated); %unsatinds=find(unsaturated); %XN=NaN*ones(size(X)); %XN(unsatinds)=X(unsatinds); %LXN=NaN*ones(size(X)); %LXN(unsatinds)=LX(unsatinds); OrdM(indframe,)=mean(X);% nonsatM(indframe,)=nanmean(XN); % LM(indframe,)=mean(LX); % nosatLM(indframe,)=nanmean(LXN);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tm=nanmean(OrdM);%The 0.05 is out of O(200) so it is very insignificant in the bright   %points while being still small in the dark O(1)points.mx=max(tm);%tm=0.005+tm/mx;tm=0.001+tm/mx;mx=max(tm);m=tm/mx;%ylm=log2(m);%figure(1)%hold off%plot(m,′k′);%set(gca,′xlim′,[1,imagewidth]);%figure(2)%hold off%plot(ylm,′k′)%set(gca,′ylim′,[-12,0])%set(gca,′xlim′,[1,imagewidth]);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function clipped=clipfrombig(oldcell,Lnew,Rnew,Dnew,Unew)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSphoto=oldcell{1,1}; % the 1st (old) imagecenterx=oldcell{2,1};%coordinates of center of old imagecentery=oldcell{3,1}; [oldrows,oldcols]=size(photo); Lold=centerx-(oldcols-1)/2; %the Left side of the old image   Rold=centerx+(oldcols-1)/2; % the Right side   Dold=centery-(oldrows-1)/2; % the Down of the old image   Uold=centery+(oldrows-1)/2; % the Upper side of the old image   greatL=Lold; %the Left side of the combined-greater image   greatR=Rold;   greatD=Dold;   greatU=Uold;   spanx=greatR-greatL+1;   spany=greatU-greatD+1;   greatX=ones(spany,1)*(greatLgreatR);%the coordinate grid of the combined-greaterimage   greatY=(greatDgreatU)′*ones(1,spanx); %Lnew; %Rnew; %Dnew; %Unew;[newLDy,newLDx]=find((greatX==Lnew) (greatY==Dnew));  [newLUy,newLUx]=find((greatX==Lnew) (greatY==Unew));  [newRDy,newRDx]=find((greatX==Rnew) (greatY=Dnew));  [newRUy,newRUx]=find((greatX=Rnew) (greatY==Unew));  [newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);  [newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);  [newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);  [newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);  clipped=photo(newLD_row-1newLU_row,newLD_colnewRD_col);   clipped=flipud(clipped);clearpackfuzzmarge=45;transit=80;%load ../decayfunm=exp(log(valint)-max(log(valint)));delta_m=(1/max(valint))*ones(size(m));load iterdecaylevels_x=2;levels_y=0;minlevel=0;Lshulaim=26;%Lshulaim=0;left=min(xint)+Lshulaim;right=max(xint);bottom=170;roof=355;imagehight=roof-bottom+1;sun=250;sunerror=1/eps;SE=ones(3,5);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cols=right-left+1;prof=profil2(cols,transit);inv_prof2=(1./prof).^2;FETHER2=ones(imagehight,1)*inv_prof2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%oldlength=length(m);m=m(1+Lshulaimoldlength);delta_m=delta_m(1+Lshulaimoldlength);sig=1;invm=1./m;invm_square=invm.^2;invm_four=invm.^4;a=(sig.^2).*invm_square;b=(delta_m.^2).*invm_four;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%invM=ones(imagehight,1)*invm;A=ones(imagehight,1)*a;B=ones(imagehight,1)*b;load ../coordm18.mat %%%% LOADING THE COORDINATESbasefilename=′../gradfilt1_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[1∶9,11∶37];centerx1=0;%this is alwayscentery1=0;%this is always%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%--%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MERGING IMAGES IN A BIG FRAME%ystep=-1;%frame_inds=36ystep30;ystep=-1;frame_inds=36ystep2;k=1;for indframe=1length(frame_inds),tic yframe_ind=frame_inds(indframe); if k==1   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 postfix];   rawim=imread(shem1);   basicim=double(flipud(rawim))-1;   X=basicim(bottomroof,leftright);   Ya=X.*invM;   saturated=(X>sun);   dilsaturated=double(dilate(saturated,SE));   DELTAYa2=(A+((X.^2).*B)).*(~dilsaturated) + sunerror.*dilsaturated;   DELTAYa2=DELTAYa2.*FETHER2;% [rows,cols]=size(Ya);% CoordYa=ones(rows,1)*(1cols);% fuzzline=round(fuzzmarge*rand(rows,1));% FuzzlineL=fuzzline*ones(1,cols);% fuzzline=cols-round(fuzzmarge*rand(rows,1));% FuzzlineR=fuzzline*ones(1,cols);% lobeseder=(CoordYa<FuzzlineL) | (CoordYa>FuzzlineR);% Ilobeseder=find(lobeseder);% clear FuzzlineR FuzzlineL lobeseder% tmp=Ya;% tmp(Ilobeseder)=NaN;% Ya=tmp;   Acell{1,1}=Ya;   Acell{2,1}=DELTAYa2;   Acell{3,1}=centerx2(yframe_ind);   Acell{4,1}=centery2(yframe_ind);  else   Acell=fusedcell;  end yframe2=netframes(yframe_ind + ystep); sframe2=sprintf(′%d′,yframe2);shem2=[basefilename sframe2 postfix]; rawim=imread(shem2); basicim=double(flipud(rawim))-1; X=basicim(bottomroof,leftright); Yb=X.*invM; saturated=(X>sun); dilsaturated=double(dilate(saturated,SE)); DELTAYb2=(A+((X.^2).*B)).*(~dilsaturated) + sunerror.*dilsaturated; DELTAYb2=DELTAYb2.*FETHER2; clear rawim basicim X% fuzzline=round(fuzzmarge*rand(rows,1));% FuzzlineL=fuzzline*ones(1,cols);% fuzzline=cols-round(fuzzmarge*rand(rows,1));% FuzzlineR=fuzzline*ones(1,cols);% lobeseder=(CoordYa<FuzzlineL) | (CoordYa>FuzzlineR);% Ilobeseder=find(lobeseder);% clear FuzzlineR FuzzlineL lobeseder% tmp=Yb;% tmp(Ilobeseder)=NaN;% Yb=tmp; Bcell(1,1}=Yb; Bcell{2,1}=DELTAYb2; Bcell{3,1}=centerx2(yframe_ind+ystep); Bcell{4,1}=centery2(yframe_ind+ystep); fusedcell=fuseml11(Acell,Bcell);%fused=double(fusedcell{1,1});%deltaz2=fusedcell{2,1};%figure(1)%imagesc(fused);axis xy;axis image%figure(2)%imagesc(deltaz2);axis xy;axis image%pause  k=k+1;  clear Acell Bcelltocendfused=double(fusedcell(1,1});deltaz2=fusedcell{2,1};greatcenterx=fusedcell{3,1};greatcentery=fusedcell{4,1};%save fethermosaic0.mat fused deltaz2 greatcenterx greatcenteryclear fusedcellfigure(3);image(fused); yshow24;set(gcf,′menubar′,′figure′)figure(4)image(150-40*log(deltaz2));yshow24set(gcf,′menubar′,′figure′)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% logfused=log(l+fused);mx=max(logfused);imx8=(2^8)/mx;logfused8=uint8(round(logfused*imx8));figure(5);image(logfused8);yshow24;set(gcf,′menubar′,′figure′)clear saturated logfused valint tikunx tikuny invM invm invm_squareclear invm_square dilsaturated flatcurve flatfun centery2 centerx2clear a b Ya Yb DELTAYb2 DELTAYa2 A B function fusedcell=fusem122(oldcell,newcell)%% oldcell and newcell contain images to put in a single% coordinate system,and the vector-coordinates of the% image centers in (x,y).% fusedcell is the cell of the updated mosaic image% k is the number of images already fused%% merging=′average′- The overlap area is an average of% all the images that exist in it,and the program% makes all weights be the same.% merging=′avdecay′-The overlap area is the average of the% old image or mosaic,withthe new image.So,older% images in the old mosaic decay exponentially in time.% merging=′split′- the overlap is cut by half - each side% comes from a different image.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DECODING THE CELLS INTO THEIR CONSTITUENTSZold =oldcell{1,1}; % the 1st (old) imageDELTAZ2old=oldcell{2,1);oldcenterx=oldcell{3,1};%coordinates of center of old imageoldcentery=oldcell{4,1};Znew=newcell{1,1}; % the 2nd (new) imageDELTAZ2new=newcell{2,1};newcenterx=newcell{3,1};newcentery=newcell{4,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; %the Left side of the new (2nd) imageRnew=newcenterx+(newcols-1)/2;Dnew=newcentery-(newrows-1)/2;Unew=newcentery+(newrows-l)/2;greatL=min([Lold,Lnew]); %the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR);%the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size{greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;DELTAZoldgreat=NaN*zeros(size(greatX));DELTAZoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=DELTAZ2old;%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack((newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;DELTAZnewgreat=NaN*zeros(size(greatX));DELTAZnewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=DELTAZ2new;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TAKING CARE OF THE OVERLAP AREAnanold=uint8(isnan(Zoldgreat));nannew=uint8(isnan(Znewgreat));fused=NaN*zeros(size(greatX));%fused(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;%fused(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;besdernew=find(-nannew);fused(besdernew)=Znewgreat(besdernew);besderold=find(-nanold);fused(besderold)=Zoldgreat(besderold);DELTA2fused=NaN*zeros(size(greatX));%DELTA2fused(newLD_row-1newLU_row,newLD_colnewRD_col)=DELTAZ2new;%DELTA2fused(oldLD_row-loldLU_row,oldLD_cololdRD_col)=DELTAZ2old;DELTA2fused(besdernew)=DELTAZnewgreat(besdernew);DELTA2fused(besderold)=DELTAZoldgreat(besderold);%Sigmaoldgreat=sqrt(DELTAZoldgreat);%Sigmanewgreat=sqrt(DELTAZnewgreat);Sigma2oldgreat=DELTAZoldgreat;Sigma2newgreat=DELTAZnewgreat;nangreat=(nanold | nannew);E=~nangreat;F=find(E);if ~isempty(F)   oldneto=NaN*ones(size(Zoldgreat));   newneto=NaN*ones(size(Znewgreat));   oldsigneto=oldneto;   newsigneto=newneto;   oldneto(F)=Zoldgreat(F);   newneto(F)=Znewgreat(F);   oldsig2neto(F)=Sigma2oldgreat(F);   newsig2neto(F)=Sigma2newgreat(F);   %size(oldneto(F))   %size(newneto(F))   %size(oldsig2neto(F)′)   %size(newsig2neto(F)′)   nominat=(oldneto(F)./oldsig2neto(F)′) + (newneto(F)./newsig2neto(F)′);   denomin=(1./oldsig2neto(F)′) + (1./newsig2neto(E)′);   sig2seam=1./denomin;   seamvalue=sig2seam.*nominat;   fused(F)=seamvalue;  DELTA2fused(F)=sigseam.^2;  DELTA2fused(F)=sig2seam;endfused=flipud(fused);DELTA2fused=flipud(DELTA2fused);%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT=COMBINED IMAGEgreatcenterx=(greatR+greatL)/2;greatcentery=(greatU+greatD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERfusedcell{1,1}=fused;fusedcell{2,1}=DELTA2fused;fusedcell{3,1}=greatcenterx;fusedcell{4,1}=greatcentery;function diffml=generallikely31(oldcell,newcell)Zold=oldcell(1,1}; %THE 1ST SUB-IMAGEZnew=newcell{1,1}; %THE 2ND SUB-IMAGEDELTAZ2old=oldcell{2,1};DELTAZ2new=newcell{2,1};nanold=uint8(isnan(Zold));nannew=uint8(isnan(Znew));nangreat=(nanold | nannew);E=~nangreat;F=find(E);Zoldneto=NaN*ones(size(Zold));Znewneto=NaN*ones(size(Znew));DELTAZ2oldneto=NaN*ones(size(Zold));DELTAZ2newneto=NaN*ones(size(Znew));bestvalue=NaN*ones(size(Zold));if isempty(F) diffml=Inf;else Zoldneto(F)=Zold(F);   Znewneto(F)=Znew(F); DELTAZ2oldneto(F)=DELTAZ2old(F); DELTAZ2newneto(F)=DELTAZ2new(F); Sigma2old=DELTAZ2oldneto; Sigma2new=DELTAZ2newneto; nominat=(Zoldneto(F)./Sigma2old(F)) + (Znewneto(F)./Sigma2new(F)); denomin=(1./Sigma2old(F)) + (1./Sigma2new(F)); sig2seam=1./denomin; bestvalue(F)=sig2seam.*nominat; Eold=((bestvalue(F) -Zoldneto(F)).^2)./Sigma2old(F); Enew=((bestvalue(F) -Znewneto(F)).^2) ./Sigma2new(F); E=Eold+Enew; SSD=sum(E); generalized_denomin=sum(denomin); diffm1=SSD/generalized_denomin;% diffm1=mean(E);endclearload decayfunlevels_x=3;levels_y=0;minlevel=0;left=min(xint);right=max(xint);bottom=170;roof=355;imagehight=roof-bottom+1;sun=250;sunerror=1/eps;SE=ones(3,5);%yepsilon=0.01;yepsilon=eps;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%basefilename=′gradfilt1_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[1∶9,11∶37];m=calchorizprof4(basefilename,postfix,netframes,left,right,bottom,roof);clear xint valint%m=valint/max(valint);%delta_m=1/max(valint);%sig=1;delta_m=0.01*ones(size(m));sig=0.5;invm=1./m;invm_square=invm.^2;c=(delta_m.^2).*invm_square;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%M =ones(imagehight,1)*m;epsM=yepsilon*M;invM=ones(imagehight,1)*invm;C =ones(imagehight,1)*c;SIG2=(sig.^2);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%centerx2=zeros(1,length(netframes));centery2=zeros(1,length(netframes));tikunx=zeros(1,length(netframes));tikuny=zeros(1,length(netframes));ystep=1;centerx2guess=47*ystep;%initializing the disparitycentery2guess=0;%frame_inds=11ystep18;frame_inds=1ystep35;k=1;for indframe=1length(frame_inds),  %%%%%%%%%%%%%%%%%%%%%%%%%%%% READING AND PREPARING THE RAW IMAGES  yframe_ind=frame_inds(indframe);tic  if k==1   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 postfix];   rawim=imread(shem1);   basicim=double(flipud(rawim))-1;   X=basicim(bottomroof,leftright);   Za=log(yepsilon+X.*invM);   saturated=(X>sun);   dilsaturated=double(dilate(saturated,SE));   DELTAZa2=( (1./((epsM+X).^2)) .* (SIG2+(C.*(X-^2))) ).*(~dilsaturated) +sunerror.*dilsaturated;   centerx1=0;centery1=0;   holdcell{1,1}=Za;   holdcell{2,1}=DELTAZa2;   holdcell{3,1}=centerx1;   holdcell{4,1}=centery1;   fusedcell=holdcell;  else   holdcell=fusedcell;  end  yframe2=netframes(yframe_ind+ystep)  sframe2=sprintf(′%d′,yframe2);  shem2=[basefilename sframe2 postfix];  rawim=imread(shem2);  basicim=double(flipud(rawim))-1;  X=basicim(bottomroof,leftright);  Zb=log(yepsilon+X.*invM);  saturated= (X>sun);  dilsaturated=double(dilate(saturated,SE));  DELTAZb2=( (1./((epsM+X).^2)) .* (SIG2+(C.*(X.^2))) ).*(~dilsaturated) +sunerror.*dilsaturated;  hnewcell{1,1}=Zb;  hnewcell{2,1}=DELTAZb2;  if k==1   hnewcell{3,1}=centerx2guess;   hnewcell{4,1}=centery2guess;  else   hnewcell{3,1}=centerx2(old_yframe_ind)+centerx2guess;   hnewcell{4,1}=centery2(old_yframe_ind)+centery2guess;  end  newcell=hnewcell;  %if yframe2>11,   %figure(1); image(42*fusedcell{1,1});colormap(gray(256));axis xy;   %figure(3); imagesc(-log(fusedcell{2,1}));colormap(gray(256));axis xy;   %figure(2); image(42*Zb);colormap(gray(256));axis xy;   %figure(4); imagesc(-log((DELTAZb2)));colormap(gray(256));axis xy;   %drawnow   %1111111   %pause  %endclear rawim basicim %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DERIVING TRANSLATION BETWEEN IMAGES[tikunx(yframe_ind+ystep),tikuny(yframe_ind+ystep)]=gmultireml31(hldcell,hnewcell,levels_x,levels_y,minlevel);% centerx2(yframe_ind+ystep)=newcell{3,1} + tikunx(yframe_ind+ystep); centery2(yframe_ind+ystep)=newcell{4,1} + tikuny(yframe_ind+ystep); [tikunx(yframe_ind+ystep) tikuny(yframe_ind+ystep)] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUSING THE IMAGES newcell{3,1}=centerx2(yframe_ind+ystep); newcell{4,1}=centery2(yframe_ind+ystep); %%%%%%%%% NEED TO CORRECT THIS WITH THE NEW KIND OF FUSION fusedcell=fusem122(fusedcell,newcell); old_yframe_ind=yframe_ind+ystep; k=k+1;toc save yytemp.mat centerx2 centery2 tikunx tikunyend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear newcell oldcellsave coordml40.mat centerx2 centery2 tikunx tikunyfused=fusedcell{1,1};deltaz2=fusedcell{2,1};save fused40.mat fused deltaz2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear fusedcell figure(2); imagesc(fused); yshow24; set(gcf,′menubar′,′figure′) figure(3); imagesc(-log((DELTAZb2))); yshow24; set(gcf,′menubar′, ′figure′) %darklevel=3;%sunlevel=249;% clippedl=validimag(bottomroof,leftright).*flatfun;% hclippedl=log2(clippedl+1);%LOG to deal with the high dynamic rangefunction [xmove,ymove]=gmlscalemov31(Acell,OBcell)%% [xmove,ymove]=mlscalemov4(Acell,OBcell)% Acell and OBcell are cells containing images at a specific% scale,and each of their initial center-coordinates.% The routine searches for the best GLOBAL CORRELATION between% the images,in a 5×5 neighborhood around the initial% centers.% THE IMAGES NEED NOT BE OF THE SAME SIZE.% [xmove,ymove]give the shifting in (x,y) coordinates of% Bcell relative to Acell - that is the correction of the% center vector of Bcell.%% See also scalprod2 putinbig7Bcell=OBcell;%movementx=[-212];%movementy=[-212];movementx=[-111];movementy=[-111];centerx1=Acell{3,1};centery1=Acell{4,1};centerx2=Bcell{3,1};centery2=Bcell{4,1};for xind=1length(movementx), centerx2_updated =centerx2+movementx(xind); Bcell{3,1}=centerx2_updated;  for yind=1length(movementy),   centery2_updated=centery2+movementy(yind);   Bcell{4,1}=centery2_updated;   [oldgreatcell,newgreatcell]=putbigml3(Acell,Bcell);   %photoA=oldgreatcell{1,1}; %THE 1ST SUB-IMAGE   %photoB=newgreatcell{1,1}; %THE 2ND SUB-IMAGE   %delta2A=oldgreatcell{2,1}; %THE 1ST SUB-IMAGE   %delta2B=newgreatcell{2,1}; %THE 2ND SUB-IMAGE   %figure(5); image(42*photoA);colormap(gray(256));axis xy;   %figure(6); image(42*photoB);colormap(gray(256));axis xy;   %figure(7); imagesc(-log((delta2A)));colormap(gray(256));axis xy;   %figure(8); imagesc(-log((delta2B)));colormap(gray(256));axis xy;   %22222   %pause   diffm1(xind,yind)=generallikely31(oldgreatcell,newgreatcell);  endendmindiff=min(diffml);%diffm1[xmovementind,ymovementind]= find(diffm1==mindiff);xmove=round(mean(movementx(xmovementind)));ymove=round(mean(movementy(ymovementind)));function[trans1atex,translatey]=gmultireml31(oldcell,newcell,levels_x,levels_y,minlevel)%% [translatex,translatey]=multiml5(oldcell,newcell,levels_x,levels_y,minlevel)% oldcell and newcell have two input images at raw resolution% who will be registerred.% THE REGISTRATION IS DONE IN A COARSE-TO-FINE MANNER.The output vector% is the (x,y) shift of the image in newcell relative to the one in% oldcell.%% See also scaledmove7 mizoor3 putinbig6shape=′valid′;[oldgreatcell,newgreatcell]=puthullml2(oldcell,newcell);za=oldgreatcell{1,1};DELTAZ2a=oldgreatcell{2,1};centerx1=0;centery1=0;Zb=newgreatcell{1,1};DELTAZ2b=newgreatcell{2,1};centerx2(1)=0;centery2(1)=0;clear oldgreatcell newgreatcellif ~(levels_x==levels_y)  [a,maindirection]=max([levels_x,levels_y]);  if maindirection==1 %%% x-direction is the main movementunceratinty   levelseq=[levels_x-1levels_y+1];   for levelind=1length(levelseq),   level=levelseq(levelind);   if ~(levelind==1)   leveldifference=levelseq(levelind-1)-level;   centerx2(levelind)=centerx2temp*(2^leveldifference);   centery2(levelind)=centery2temp;   end   [photoA,delta2A]=mizoorml21_xy(Za,DELTAZ2a,level,levels_y,shape);   [photoB,delta2B]=mizoorml21_xy(Zb,DELTAZ2b,level,levels_y,shape);   Acell{1,1}=photoA;   Acell{2,1}=delta2A;   Acell{3,1}=centerx1;   Acell{4,1}=centery1; % figure(5); image(42*photoA); colormap(gray(256)); axis xy; % figure(6); image(42*photoB); colormap(gray(256)); axis xy; % figure(7); imagesc(-log((delta2A))); colormap(gray(256)); axis xy; % figure(8); imagesc(-log((delta2B))); colormap(gray(256)); axis xy; % 22222 % pause   Bcell{1,1}=photoB;   Bcell{2,1}=delta2B;   Bcell{3,1}=centerx2(levelind);   Bcell{4,1}=centery2(levelind);[xmove,ymove]=gmlscalemov31(Acell,Bcell);   centerx2temp=centerx2(levelind)+xmove;   centery2temp=centery2(levelind)+ymove;   end else %%% y-direction is the main movement unceratinty   levelseq=[levels_y-1levels_x+1];   for levelind=1length(levelseq),   level=levelseq(levelind);   if -(levelind==1)   leveldifference=levelseq(levelind-1)-level;   centerx2(levelind)=centerx2temp;   centery2(levelind)=centery2temp*(2^leveldifference);   end   [photoA,delta2A]=mizoorml21_xy(Za,DELTAZ2a,levels_x,level,shape);   [photoB,delta2B]=mizoorml21_xy(Zb,DELTAZ2b,levels_x,level,shape);   Acell{1,1}=photoA;   Acell{2,1}=delta2A;   Acell{3,1}=centerx1;   Acell{4,1}=centery1;   Bcell{1,1}=photoB;   Bcell{2,1}=delta2B;   Bcell{3,1}=centerx2(levelind);   Bcell{4,1}=centery2(levelind);   [xmove,ymove]=gmlscalemov31(Acell,Bcell);   centerx2temp=centerx2(levelind)+xmove;   centery2temp=centery2(levelind)+ymove;   end endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if ~(levels_x==levels_y) [a,maindirection]=max([levels_x,levels_y]); if maindirection==1   levelseq=[levels_y-1minlevel];   centerx2(1)=centerx2temp*(2^ieveldifference);   centery2(1)=centery2temp; else   levelseq=[levels_x-1minlevel];   centerx2(1)=centerx2temp;   centery2(1)=centery2temp*(2^leveldifference); endelse levelseq=[levels_x-1minlevel]; centerx2(1)=0; centery2(1)=0;endfor levelind=1length(levelseq), level=levelseq(levelind); if ~(levelind==1)   leveldifference=levelseq(levelind-1)-level;   centerx2(levelind)=centerx2temp*(2^leveldifference);   centery2(levelind)=centery2temp*(2^leveldifference); end[photoA,delta2A]=mizoorml21(Za,DELTAZ2a,level,shape);  [photoB,delta2B]=mizoorml21(Zb,DELTAZ2b,level,shape);  Acell{1,1}=photoA  Acell{2,1}=delta2A;  Acell{3,1}=centerx1;  Acell{4,1}=centery1;  Bcell{1,1}=photoB;  Bcell{2,1}=delta2B;  Bcell{3,1}=centerx2(levelind);  Bcell{4,1}=centery2(levelind);  [xmove,ymove]=gmlscalemov31(Acell,Bcell);  centerx2temp=centerx2(levelind)+xmove;  centery2temp=centery2(levelind)+ymove; end translatex=centerx2temp*(2^minlevel); translatey=centery2temp*(2^minlevel); function[enhanced,mnval,mxval]=histwin4(X,Y,photo,logphoto)XS=sort(X);YS=sort(Y);L=round(XS(2));R=round(XS(3));D=round(YS(2));U=round(YS(3));%LD LU RU RD LDXp= [L L R R L];Yp= [D U U D D];XO= (L+R)/2;YO= (D+U)/2;spanx=size(photo,2);spany=size(photo,1);[LD_row,LD_col]=xy2ij_brack([L,D],spany);[RD_row,RD_col]=xy2ij_brack([R,D],spany);[LU_row,LU_col]=xy2ij_brack([L,U],spany);[RU_row,RU_col]=xy2ij_brack([R,U],spany);stam=flipud(photo);subarea=stam(LD_row-1LU_row,LD_colRD_col);mn=min(subarea);mx=max(subarea);lemata=subarea-mn;mx=max(lemata);imx16=(2^16)/mx;substretch=uint16(round(lemata*imx16));subareahist=histeq(substretch,256);mx=max(subareahist);imx8=(2^8)/double(mx);subenhanced=uint8(round(double(subareahist)*imx8));stam=flipud(logphoto);stam(LD_row-1LU_row,LD_colRD_col)=subenhanced;enhanced=flipud(stam);mnval=min(subarea);mxval=max(subarea);clearclose allpackkamut=1240;%lambda=3.5e+O4;%lambdalog=3.5e4;%lambda=2.5e+O4;lambda=le+O4;lambdalog=6e4;load ../decayfunclear m delta_mload lsatmosaic0.matLshulaim=0;oldLshulaim=0;left=min(xint)+Lshulaim;right=max(xint);bottom=170;roof=355;imagehight=roof-bottom+1;imagewidth=right-1eft+1;%sun=250;%sun=100;sun=190;dark=2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%logfused=log(1+fused);mx=max(logfused);imx8=(2^8)/mx;logfused8=uint8(round(logfused*imx8));[rows,cols]=size(logfused8);L=l;R=cols;D=1;U=rows;figure(1);hold offimage(logfused8);colormap(gray(256));axis(′image′,′off′,′xy′);%Xwide=cols;Ywide=round(Xwide*rows/cols);Xwide=cols/1.8;Ywide=round(Xwide*rows/cols);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*rows/cols;set(gcf,′position′,[1 500 Xwide Ywide]);spanx=R-L+1; spany=U-D+1; X=uint16(ones(spany,1)*(LR)); %the coordinate grid of the combined-greater image Y=uint16((DU) ′*ones(1,spanx)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COORDINATE GRID OF THE MOSAIC [bigrows,bigcols]=size(fused); fusedcell{1,1}=flipud(fused); % the 1st (old) image fusedcell{2,1}=greatcenterx+oldLshulaim/2;%coordinates of center of old image fusedcell{3,1}=greatcentery; mxx=max(fused); ybase=sqrt(2); %ybase=2; denomin=log2(ybase); fracoctav=(log2(mxx/sun)/denomin)+1; intoctav=ceil(fracoctav); goodpts=[]; for octavnum=lintoctav;  if octavnum==l   minlevel=10*dark;   maxlevel=sun;  else   minlevel=sun.*(ybase.^(octavnum-2));   maxlevel=sun.*(ybase.^(octavnum-1));  end  inrange=((fused>minlevel) (~(fused>maxlevel)));  inrange_ind=find(inrange);  lenrange=length(inrange_ind);  if lenrange<kamut   goodpts=[goodpts;inrange_ind];  else   linerand=1+round((lenrange-1)*(rand(kamut,1)));   picked_ind=inrange_ind(linerand);   goodpts=[goodpts;picked_ind]; % x=X(picked_ind); % y=Y(picked_ind); endend   x=X(goodpts);   y=Y(goodpts);   length(x)   figure(1);   hold off   image(logfused8);colormap(gray(256));axis(′image′,′off′,′xy,);   hold on% plot(x,y,′r.′);axis xy   plot(x,y, ′r*′);axis xy  set(gca,′xtick′,0102400);set(gca,′ytick′,610250);%grid on; axis on;  drawnow  pause(10)[rows,cols]=size(fused);greatL=(greatcenterx-(cols-1)/2); %the Left side of the imagegreatR=(greatcenterx+(cols-1)/2); % the Right sidegreatD=(greatcentery-(rows-1)/2); % the Down of the imagegreatU=(greatcentery+(rows-1)/2); % the Upper side of the imageclear deltaz2 fused logfused8 X Y x y fusedcell logfused8 logfusedclear inrange inrange_ind intoctav lenrange linerand maxlevel minlevelclear picked_ind imx8xpix=1imagewidth;ystep=1;frame_inds=1ystep36;load ../coordml8.mat%%%% LOADING THE COORDINATES%basefilename=′../gradfilt1 ′;basefilename=′/proj/cavel/users/yoav/MOSAIC/MAY26/gradfilt1_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[1∶9 ,11∶37];ILUTZIM=sparse([]);tic%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MERGING IMAGES IN A BIG FRAMEfor indframe=1length{frame_inds), yframe_ind=frame_inds(indframe); yframe=netframes(yframe_ind); sframe=sprintf(′%d′,yframe); shem2=[basefilename sframe postfix]; rawim=imread(shem2); basicim=double(flipud(rawim))-1; photo=uint8(basicim(bottomroof,leftright)); centerx=centerx2(yframe_ind); centery=centery2(yframe_ind);   oldcell{1,1}=photo; % the 1st (old) image   oldcell{2,1}=centerx; %coordinates of center of old image oldcell{3,1}=centery;[beseder,xbeseder,ybeseder,valbeseder]=subpts21(oldcell,greatL,greatR,greatD,greatU,goodpts); beseder_ind=find(beseder); xsmall=xbeseder(beseder_ind); ysmall=ybeseder(beseder_ind); valsmall=valbeseder(beseder_ind); [rows,cols]=size(photo); if ~(indframe==1)   bothbeseder=(oldbeseder beseder);   beseder_ind=find(bothbeseder);   OLDxsmall=xbesederold(beseder_ind);   OLDysmall=ybesederold(beseder_ind);   OLDval= valold(beseder_ind);   NEWxsmall=xbeseder(beseder_ind);   NEWysmall=ybeseder(beseder_ind);NEWval =valbeseder(beseder ind);   XPIX=ones(length(OLDxsmall),1)*xpix;   OLDXSMALL=OLDxsmall*ones(1,imagewidth);   IOLD=find(~[abs(XPIX-OLDXSMALL)′));   chunk=sparse((zeros(length(OLDxsmall),imagewidth))′);   chunk(IOLD)=NEWval;% positive sign   NEWXSMALL=NEWxsmall*ones(1,imagewidth);   INEW=find(~(abs(XPIX-NEWXSMALL)′));   chunk(INEW)=-OLDval; %minus sign   chunk=sparse(chunk′);   ILUTZIM=[ILUTZIM;chunk];   %figure(2); hold off; image(oldphoto);   %colormap(gray(256)); axis(′image′,′off′,′xy′);   %set(gcf,′position′,[59 292 612 108]);   %set(gca,′position′,
);hold on   %plot(OLDxsmall,OLDysmall,′r.′);axis xy   %axison;grid on;set(gca,′xtick′,010700);set(gca,′ytick′,710250); %figure(3); hold off; image(photo);   %colormap(gray(256)); axis(′image′,′off′,′xy′);   %Xwide=cols;Ywide=round(Xwide*rows/cols);   %set(gcf,′position′,[59 100 612 108]);   %set(gca,′position′,
); hold on;   %plot(NEWxsmall,NEWysmall,′r.′); axis xy  %axis on;grid on;set(gca,′xtick′,010700);set(gca,′ytick′,710250);  end  oldbeseder=beseder;  xbesederold=xbeseder;  ybesederold=ybeseder;  valold=valbeseder;endtocILUTZ=sparse(ILUTZIM);whos ILUTZclear ILUTZIMclear A oldphoto D U L R INEW IOLD Lshulaim NEWXSMALL NEWval NEWxsmallclear NEWysmall OLDXSMALL OLDval OLDxsmall OLDysmall SE XPIX xpixclear Xwide Ywide ans beseder beseder_ind bigcols bigrows bothbesederclear bottom chunk roof cols rows fracoctav frame_inds fusedcellclear mx mxx netframes octavnum oldLshulaim oldbeseder oldcellclear photo oldcell clipped fusedclipped saturated dilsaturatedclear hetclipped beseder netfused LI LG centerx centery Dnew Doldclear Unew Uold Lnew Lold Rnew Rold centerx2 centery2 oldcols oldrowsclear bigrows bigcols greatcenterx greatcentery ystep fusedcellclear rawim basicim X big2D basefilename basephotoname SE picked_indclear centerx centery imagehight indframe postfix right roof spanx spanyclear rawname sframe shem2 tikunx tikuny yframe yframe_ind sun sunerrorclear tikunx tikuny tzamtzam valbeseder valint valold valsmall vecclear xbeseder xbesederold xsmall ybeseder ybesederold ysmall ystepclose allpack%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ynegat=(ILUTZ<0); yposit=(ILUTZ>0); ynozer=(ynegat | yposit); ymeasured=sum(ynozer); meanmeasured=mean(ymeasured); %saf=min([1+(meanmeasured/10),5]); saf=l; reliable=find(~(ymeasured<saf)); %mnn=max([min(reliable)-1,1]); mnn= min(reliable)-1; %mxx=min([max(reliable)+1,imagewidth]); mxx= max(reliable)+1; %I=find(ymeasured); %endmesures=min ([max(I)+1,imagewidth]); mnnozer=ynozer(,1mnn); J=[]; if mnn>1 J=find ( (sum(mnnozer′) ) ′); elseif mnn==1   J=find(mnnozer); end NEWILUTZ=ILUTZ; NEWILUTZ(J,)=0; mxnozer=ynozer(,mxximagewidth); J=[]; if (mxx<imagewidth) J=find((sum(mxnozer′)) ′); elseif (imagewidth==mxx)   J=find(mxnozer); end NEWILUTZ(J, )=0; ynegat=(NEWILUTZ<0); yposit=(NEWILUTZ>0); ynozer=(ynegat | yposit); ymeasured=sum(ynozer); yesmeasured=find(ymeasured>0); mnn=min(yesmeasured); mxx=max(yesmeasured); oldimagewidth=imagewidth; imagewidth=mxx-mnn+1; NETILUTZ=NEWILUTZ(,mnnmxx); clear I J NEWILUTZ mnnozer mxnozer reliable yesmeasured clear ynegat ynozer yposit meanmeasured %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %imagewidth=oldimagewidth; tic amain=ones(imagewidth,1);Amain=diag(amain);aoff1=-0.5*ones(imagewidth-1,1);aoff1(1)=-1;Aoff1=diag(aoff1,1);aoff_1=-0.5*ones(imaqewidth-1,1);aoff_1(imagewidth-1)=-1;Aoff_1=diag(aoff_1,-1);SM=Amain+Aoff1+Aoff_1;SM(1,)=0;SM(imagewidth,)=0;clear Amain Aoff1 Aoff_1 aoff_1 aoff1 amainsqlambda=sqrt(lambda);HALAKUT=sqlambda*sparse(SM);LOCENERG=NETILUTZ.^2;%locenerg=sqrt(sum(LOCENERG′))′;locenerg=(sum(LOCENERG′))′;maxenerg=sqrt(max(locenerg));NORMILUTZ=NETILUTZ/235;ALLILUTZ=[NORMILUTZ;HALAKUT];clear HALAKUT SM LOCENERG NORMILUTZ NETILUTZpack[U,S,V]=svd(full(ALLILUTZ),0);s=diag(S);[Y,I]=min(s)clear S UM=V(,I);[Mxx,I]=max (abs(M));Imn=min (I); mpeak=M(Imn);m=M/mpeak;m_amud=m;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%invC=(NORMILUTZ′*NORMILUTZ) + 0.5*eye(imagewidth);invC=(ALLILUTZ′*ALLILUTZ) + 0.5*eye(imagewidth);C=inv(invC);%NORMILUTZ=NETILUTZ/235;%Sr=m_amud′*((NORMILUTZ′*NORMILUTZ)*m_amud);Sr=m_amud′*(invC*m_amud);%LOCENERG=NORMILUTZ.^2;LOCENERG=ALLILUTZ.^2;%locenerg=sqrt(sum(LOCENERG′))′;locenerg=(sum(LOCENERG′))′;n_keilu=sum(locenerg);dof=length(m_amud)-1;syx=sqrt(Sr/(n_keilu-dof));syx2=syx^2;C_norm=C* syx2;sig_m=sqrt(diag(C_norm));figure(5)plot(log10(sig_m(1imagewidth)./m_amud(mnn mxx)))grid on%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tocclear Vmlog=log(abs(m));% P=mlog ′;P=zeros(1,oldimagewidth);P(mnn mxx)=mlog;leftmarg=1mnn-1;P(fliplr(leftmarg))=2*mlog(1)-mlog(leftmarg+1);rightmarg=mxx+1oldimagewidth;oldrightmarg=rightmarg-oldimagewidth+imagewidth;P(fliplr(rightmarg))=2*mlog(imagewidth)-mlog(oldrightmarg-1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%L=length(P);Imain=ones(L,1);if mnn>1 Imain(1 mnn)=1;endif mxx<L Imain(mxx+1L)=1;endImat=sparse(diag(Imain));amain=ones(L,1);Amain=diag(amain);aoff1=-0.5*ones(L-1,1);aoff1(1)=-1;Aoff1=diag(aoff1,1);aoff_1=-0.5*ones(L-1,1);aoff_1(L-1)=-1;Aoff_1=diag(aoff_1, -1);S=sparse(Amain+Aoff1+Aoff_1);S(1,)=0;S(L,)=0;clear Amain Aoff1 Aoff_1 aoff_1 aoff1 amainStS=S ′* S;clear S;R=I mat + lambdalog*StS; iR=inv(R); Y=iR*P′; Ynorm=Y-max(Y); smoothm=exp(Ynorm′); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mlogdelta=(sig_m(1imagewidth)./m_amud(mnnmxx)); %P=mlog′; dP=zeros(1,oldimagewidth); dP(mnnmxx)=mlogdelta; leftmarg=1mnn-1; dP(fliplr(leftmarg))=2*mlogdelta(1)-mlogdelta(leftmarg+1); rightmarg=mxx+1oldimagewidth; oldrightmarg=rightmarg-oldimagewidth+imagewidth; dP(fliplr(rightmarg))=2*mlogdelta(imagewidth)-mlogdelta(oldrightmarg-1); dY=iR*dP′; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% longm=exp(P); longdm=smoothm.*dY′; m_lemata=smoothm-longdm; m_lemala=smoothm+longdm; figure(1) hold off plot(longm,′b′) hold on plot(smoothm,′r′) plot(m_lemata,′m′) plot(m_lemala,′m′) set(gca,′xlim′,[1,oldimagewidth]); set(gca,′ylim′,
); dlog2m=dY′/log(2); log2m_lemata=log2(smoothm)-dlog2m; log2m_lemala=log2(smoothm)+dlog2m; figure(2) hold off plot(log2(abs(longm)),′b′) hold on plot(log2(smoothm),′r′) plot(log2m_lemata,′m′) plot(log2m_lemala,′m′) set(gca,′xlim′,[1,oldimagewidth]); m=smoothm; delta_m=longdm; save randecay37.mat xint m delta_m Ynorm smoothm C ILUTZ sig_m mnn mxxload ../decayfunorigm=valint/max(valint);figure(1)hold onplot(origm,′g′)set(gca,′xlim′,[1,oldimagewidth]);figure(2)hold onplot(log2(origm),′g′)set(gca,′xlim′,[1,oldimagewidth]);function[small,delta2small]=mizoorml21(big,delta2big,levels,shape);newA=big;newdelta2A=delta2big;for iteration=1levels, [reducedA,reducedelta2A]=reduceml21(newA,newdelta2A,shape); newA=reducedA; newdelta2A=reducedelta2A;endsmall=newA;delta2small=newdelta2A;function{small,delta2small]=mizoorml21_xy(big,delta2big,levels_x,levels_y,shape);newA=big;newdelta2A=delta2big;if ~(levels_x==levels_y) for iteration=1min([levels_x,levels_y]),   [reducedA_x,reducedelta2A_x]=reduceml21_x(newA,newdelta2A,shape);   [reducedA,reducedelta2A]=reduceml21_y(reducedA_x,reducedelta2A_x,shape);   newA=reducedA;   newdelta2A=reducedelta2A; end [a,maindirection]=max([levels_x,levels_y]); if maindirection==1 %%% x-direction is the main movement unceratinty   for iteration=levels_y+1 levels_x,   [reducedA,reducedelta2A]=reduceml21_x(newA,newdelta2A,shape);   newA=reducedA;   newdelta2A=reducedelta2A;  end  else%%% y-direction is the main movementunceratinty   for iteration=levels_x+1 levels_y,   [reducedA,reducedelta2A]=reduceml21_y(newA,newdelta2A,shape);   newA=reducedA;   newdelta2A=reducedelta2A;   end endelse levels=levels_x; for iteration=1levels,   [reducedA_x,reducedelta2A_x]=reduceml21_x(newA,newdelta2A,shape);[reducedA,reducedelta2A]=reduceml21_y(reducedA_x,reducedelta2A_x,shape);   newA=reducedA;  newdelta2A=reducedelta2A; endendsmall=newA;delta2small=newdelta2A;function diffm1=mostlikely2(oldcell,newcell)%% corrcoef=scalprod2(Zoldgreat,Znewgreat)% Calculates the correlation coefficient (″normalized″)% between partly overlapping images Zoldgreat and Znewgreat.% The result iscorrcoefZold=oldcell{1,1}; %THE 1ST SUB-IMAGEZnew=newcell{i,1}; %THE 2ND SUB-IMAGEDELTAZ2old=oldcell{2,1};DELTAZ2new=newcell{2,1};nanold=uint8(isnan(Zold));nannew=uint8(isnan(Znew));nangreat= (nanold ) nannew);E=~nangreat;F=find(E);Zoldneto=NaN*ones(size(Zold));Znewneto=NaN*ones(size(Znew));DELTAZ2oldneto=NaN*ones(size(Zold));DELTAZ2newneto=NaN*ones(size(Znew));bestvalue=NaN*ones(size(Zold));if isempty(F)   diffml=Inf;else  Zoldneto(F)=Zold(F);  Znewneto(F)=Znew(F);   DELTAZ2oldneto(F)=DELTAZ2old(F);  DELTAZ2newneto(F)=DELTAZ2new(F);% Sigmaold=sqrt(DELTAZ2oldneto);% Sigmanew=sqrt(DELTAZ2newneto);   Sigma2old=DELTAZ2oldneto;   Sigma2new=DELTAZ2newneto;  nominat=(Zoldneto(F)./Sigma2old(F)) + (Znewneto(F)./Sigma2new(F));   denomin=(1./Sigma2old(F)) + (1./Sigma2new(F));   sig2seam=1./denomin;   bestvalue(F)=sig2seam.*nominat;% Eold=( (bestvalue(F) -Zoldneto(F))./Sigmaold(F) ).^2;% Enew=( (bestvalue(F) -Znewneto(F)) ./Sigmanew(F) ).^2;   Eold=((bestvalue(F) -Zoldneto(F)) .^2) ./Sigma2old(F);   Enew=((bestvalue(F) -Znewneto(F)) .^2) ./Sigma2new(F);   E=Eold+Enew;   diffm1=mean(E);endfunction [translatex,translatey]=gmultireml31(oldcell,newcell,levels_x,levels_y,minlevel)%% [translatex,translatey]=multiml5(oldcell,newcell,levels_x,levels_y,minlevel)% oldcell and newcell have two input images at raw resolution% who will be registerred.% THE REGISTRATION IS DONE IN A COARSE-TO-FINE MANNER.The output vector% is the (x,y) shift of the image in newcell relative to the one in% oldcell.%% See also scaledmove7 mizoor3 putinbig6shape=′valid′;[oldgreatcell,newgreatcell]=puthullml2(oldcell,newcell);Za=oldgreatcell{1,1};DELTAZ2a=oldgreatcell{2,1};centerx1=0;centery1=0;Zb=newgreatcell{1,1};DELTAZ2b=newgreatcell{2,1};centerx2(1)=0;centery2(1)=0;clear oldgreatcell newgreatcellif -(levels_x==ievels_y)   [a,maindirection]=max([levels_x,levels_y]);   if maindirection==1%%% x-direction is the main movementunceratinty   levelseq=[levels_x-1levels_y+1];   for levelind=1length(levelseq),   level=levelseq(levelind);   if ~(levelind==1)   leveldifference=levelseq(levelind-1)-level;   centerx2(levelind)=centerx2temp*(2^leveldifference);   centery2(levelind)=centery2temp;   end   [photoA,delta2A]=mizoorml21_xy(Za,DELTAZ2a,level,levels_y,shape);   [photoB,delta2B]=mizoorml21_xy(Zb,DELTAZ2b,level,levels_y,shape);   Acell{1,1}=photoA;   Acell{2,1}=delta2A;   Acell{3,1}=centerxl;   Acell{4,1}=centery1; % figure(5); image(42*photoA); colormap(gray(256)); axis xy; % figure(6); image(42*photoB); colormap(gray(256)); axis xy; % figure(7); imagesc(-log((delta2A))); colormap(gray(256)); axis xy; % figure(8); imagesc(-log((delta2B))); colormap(gray(256)); axis xy; % 22222 % pause   Bcell{1,1}=photoB;   Bcell{2,1}=delta2B;Bcell{3,1}=centerx2(levelind);   Bcell{4,1}=centery2(levelind);   [xmove,ymove]=mlscalemov21(Acell,Bcell);   centerx2temp=centerx2(levelind)+xmove;   centery2temp=centery2(levelind)+ymove;   end else %%% y-direction is the main movement unceratinty   levelseq=[levels_y-1levels_x+l];   for levelind=1length(levelseq),   level=levelseq(levelind);   if ~(levelind==1)   leveldifference=levelseq(levelind-1)-level;   centerx2(levelind)=centerx2temp;   centery2(levelind)=centery2temp*(2^leveldifference);   end   [photoA,delta2A]=mizoorml21_xy(Za,DELTAZ2a,levels_x,level,shape);   [photoB,delta2B]=mizoorml21_xy(Zb,DELTAZ2b,levels_x,level,shape);   ′Acell{1,1}=photoA;   Acell{2,1}=delta2A;   Acell{3,1}=centerxl;   Acell{4,1}=centeryl;   Bcell{1,1}=photoB;   Bcell{2,1}=delta2B;   Bcell{3,1}=centerx2(levelind);   Bcell{4,1}=centery2(levelind);   [xmove,ymove]=mlscalemov21(Acell,Bcell);   centerx2temp=centerx2(levelind)+xmove;   centery2temp=centery2(levelind)+ymove;   end endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if ~(levels_x==levels_y)  [a,maindirection]=max([levels_x,levels_y]);  if maindirection==1   levelseq=[levels_y-1minlevel];   centerx2(1)=centerx2temp*(2^leveldifference);   centery2(1)=centery2temp;  else   levelseq=[levels_x-1minlevel];   centerx2(1)=centerx2temp;   centery2(1)=centery2temp*(2^leveldifference); endelse levelseq=[levels_x-1minlevel]; centerx2(1)=0; centery2(1)=0;endfor levelind=1length(levelseq), level=levelseq(levelind); if ~(levelind==1)   leveldifference=levelseq(levelind-1)-level;centerx2(levelind)=centerx2temp*(2^leveldifference);   centery2(levelind)=centery2temp*(2^leveldifference); end [photoA,delta2A]=mizoorml21(Za,DELTAZ2a,level,shape); [photoB,delta2B]=mizoorml21(Zb,DELTAZ2b,level,shape); Acell{1,1}=photoA; Acell{2,1}=delta2A; Acell{3,1}=centerxl; Acell{4,1}=centeryl; Bcell{1,1}=photoB; Bcell{2,1}=delta2B; Bcell(3,1}=centerx2(levelind); Bcell{4,1}=centery2(levelind); [xmove,ymove]=mlscalemov21(Acell,Bcell); centerx2temp=centerx2(levelind)+xmove; centery2temp=centery2(levelind)+ymove;endtranslatex=centerx2temp*(2^minlevel);translatey=centery2temp*(2^minlevel);clearclose allpackkamut=1240;%lambda=3.5e+04;%lambdalog=3.5e4;%lambda=2.5e+04;lambda=le+04;lambdalog=6e4;load ../decayfunclear m delta_mload lsatmosaic0.matLshulaim=0;oldLshulaim=0;left=min(xint)+Lshulaim;right=max(xint);bottom=170; roof=355;imagehight=roof-bottom+1;imagewidth=right-left+1;%sun=250;%sun=100;sun=190;dark=2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%logfused=log(l+fused);mx=max(logfused);imx8=(2^8)/mx;logfused8=uint8(round(logfused*imx8));[rows,cols]=size(logfused8);L=i;R=cols;D=1;U=rows;figure(1);hold offimage(logfused8);colormap(gray(256));axis(′image′,′off′,′xy′);%Xwide=cols;Ywide=round(Xwide*rows/cols);Xwide=cols/1.8;Ywide=round(Xwide*rows/cols);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*rows/cols;set(gcf,′position′,[1 500 Xwide Ywide]);spanx=R-L+1;spany=U-D+1;X=uint16(ones(spany,1)*(LR)); %the coordinate grid of the combined-greater imageY=uint16((DU)′*ones(1,spanx));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COORDINATE GRID OF THE MOSAIC[bigrows,bigcols]=size(fused);fusedcell{1,1}=flipud(fused); % the 1st (old) imagefusedcell{2,1}=greatcenterx+oldLshulaim/2; %coordinates of center of old imagefusedcell{3,1}=greatcentery;mxx=max(fused);%ybase=sqrt(2);%denomin=log2(ybase);%fracoctav=(log2(mxx/sun)/denomin) +1;%intoctav=ceil(fracoctav);[rows,cols]=size(fused);greatL=(greatcenterx-(cols-1)/2); %the Left side of the imagegreatR=(greatcenterx+(cols-1)/2); % the Right sidegreatD=(greatcentery-(rows-1)/2); % the Down of the imagegreatU={greatcentery+(rows-1)/2); % the Upper side of the imageclear deltaz2 fused logfused8 X Y x y fusedcell logfused8 logfusedclear inrange inrange_ind intoctav lenrange linerand maxlevel minlevelclear picked_ind imx8xpix=1imagewidth;ystep=1;frame_inds=1ystep36;.load ../coordm18.mat%%%% LOADING THE COORDINATESbasefilename=′/proj/cavel/users/yoav/MOSAIC/MAY26/gradfilt1_′;postfix=′.tif′; basephotoname=′photo′; rawname=′basicim′;netframes=[19 ,1137];ILUTZIM=sparse([]);tic%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MERGING IMAGES IN A BIG FRAMEfor indframe=1length(frame_inds)-1, yframe_ind=frame_inds(indframe);yframe=netframes(yfrsme_ind);sframe=sprintf(′%d′,yframe);shem2=[basefilename sframe postfix];rawim=imread(shem2);basicim=double(flipud(rawim))-1;photo=uint8(basicim(bottomroof,leftright));centerx=centerx2(yframe_ind);centery=centery2(yframe_ind);oldcell{1,1}=photo;% the 1st (old) imageoldcell{2,1}=centerx; %coordinates of center of old imageoldcell{3,1}=centery;oldphoto=photo;yframe_ind=frame_inds(indframe+1);yframe=netframes(yframe_ind);sframe=sprintf(′%d′,yframe);shem2=[basefilename sframe postfix];rawim=imread(shem2);basicim=double(flipud(rawim))-1;photo=uint8(basicim(bottomroof,leftright));centerx=centerx2(yframe_ind);centery=centery2(yframe_ind);newcell{1,1}=photo; % the 1st (old) imagenewcell{2,1}=centerx; %coordinates of center of old imagenewcell{3,1}=centery;[oldselect,newselect]=ptspick7(oldcell,newcell,greatL,greatR,greatD,greatU); [rows,cols]=size(photo);oldxsmall=oldselect{1,1};oldysmall=oldselect{2,1};oldvals =oldselect{3,1};newxsmall=newselect{1,1};newysmall=newselect{2,1};newvals =newselect{3,1};chunkleng=length(oldvals);chunk=sparse( (zeros(chunkleng,imagewidth) )′ );XPIX=ones(chunkleng,1)*xpix;OLDXSMALL=oldxsmall*ones(1,imagewidth);IOLD=find(~(abs(XPIX-OLDXSMALL)′));chunk(IOLD)=newvals;% positive signNEWXSMALL=newxsmall*ones(1,imagewidth);INEW=find(~(abs(XPIX-NEWXSMALL)′));chunk(INEW)=-oldvals; %minus signchunk=sparse(chunk′); ILUTZIM=[ILUTZIM;chunk];   figure(2); hold off; image(oldphoto);   colormap(gray(256));axis(′image′,′off′,′xy′);   set(gcf,′position′,[59 292 612 108]);   set(gca,′position′,
); hold on   plot(oldxsmall,oldysmall,′r.′);axis xy   axis on;grid on; set(gca,′xtick′,010700);set(gca,′ytick′,710250);   figure(3); hold off; image(photo);   colormap(gray(256));axis(′image′,′off′,′xy′);   set(gcf,′position′,[59 100 612 108]);   set(gca,′position′,
); hold on;   plot(newxsmall,newysmall,′r.′);axis xy   axis on;grid on; set(gca,′xtick′;010700);set(gca,′ytick′,710250);pauseendtoc5555555pauseILUTZ=sparse(ILUTZIM);whos ILUTZclear ILUTZIMclear A oldphoto D U L R INEW IOLD Lshulaim NEWXSMALL NEWval NEWxsmallclear NEWysmall OLDXSMALL OLDval OLDxsmall OLDysmall SE XPIX xpixclear Xwide Ywide ans beseder beseder_ind bigcols bigrows bothbesederclear bottom chunk roof cols rows fracoctav frame_inds fusedcellclear mx mxx netframes octavnum oldLshulaim oldbeseder oldcellclear photo oldcell clipped fusedclipped saturated dilsaturatedclear netclipped beseder netfused LI LG centerx centery Dnew Doldclear Unew Uold Lnew Lold Rnew Rold centerx2 centery2 oldcols oldrowsclear bigrows bigcols greatcenterx greatcentery ystep fusedcellclear rawim basicim X big2D basefilename basephotoname SE picked_indclear centerx centery imagehight indframe postfix right roof spanx spanyclear rawname sframe shem2 tikunx tikuny yframe yframe_ind sun sunerrorclear tikunx tikuny tzamtzam valbeseder valint valold valsmall vecclear xbeseder xbesederold xsmall ybeseder ybesederold ysmall ystepclose allpack%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ynegat=(ILUTZ<0);yposit=(ILUTZ>0);ynozer=(ynegat | yposit);ymeasured=sum(ynozer);meanmeasured=mean(ymeasured);%saf= min([ 1+(meanmeasured/10), 5]);saf=1;reliable=find(~(ymeasured<saf));%mnn=max([min(reliable)-1, 1]);mnn==min(reliable)-1;%mxx=min ([max(reliable)+1,imagewidth]);mxx= max(reliable)+1;%I=find(ymeasured);%endmesures=min([max(I)+1,imagewidth]);mnnozer=ynozer(,1mnn);J=[ ];if mnn>1 J=find((sum(mnnozer′)) ′);elseif mnn==1 J=find(mnnozer);endNEWILUTZ=ILUTZ;NEWILUTZ(J,)=0;mxnozer=ynozer(,mxximagewidth);J=[ ];if(mxx<imagewidth) J=find((sum(mxnozer′)) ′);elseif (imagewidth==mxx) J=find(mxnozer);endNEWILUTZ(J,)=0;ynegat=(NEWILUTZ<0);yposit=(NEWILUTZ>0);ynozer=(ynegat | yposit);ymeasured=sum(ynozer);yesmeasured=find(ymeasured>0);mnn=min(yesmeasured);mxx=max(yesmeasured);oldimagewidth=imagewidth;imagewidth=mxx-mnn+1;NETILUTZ=NEWILUTZ(,mnnmxx);clear I J NEWILUTZ mnnozer mxnozer reliable yesmeasuredclear ynegat ynozer yposit meanmeasured%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%imagewidth=oldimagewidth;ticamain=ones(imagewidth,1);Amain=diag(amain);aoff1=-0.5*ones(imagewidth-1,1);aoff1(1)=-1;Aoff1=diag(aoff1,1);aoff_1=-0.5*ones(imagewidth-1,1);aoff_1(imagewidth-1)=-1;Aoff_1=diag(aoff_1,-1);SM=Amain+Aoff1+Aoff_1;SM(i,)=0;SM(imagewidth,)=0;clear Amain Aoff1 Aoff_1 aoff_1 aoff1 amainsqlambda=sqrt(lambda);HALAKUT=sqlambda*sparse(SM);LOCENERG=NETILUTZ.^2;%locenerg=sqrt(sum(LOCENERG′))′;locenerg=(sum(LOCENERG′))′;maxenerg=sqrt(max(locenerg));NORMILUTZ=NETILUTZ/235;ALLILUTZ=[NORMILUTZ;HALAKUT];clear HALAKUT SM LOCENERG NORMILUTZ NETILUTZpack[U,S,V]=svd(full(ALLILUTZ),0);s=diag(S);[Y,I]=min(s)clear S UM=V(,I);[Mxx,I]=max(abs(M));Imn=min(I);mpeak=M(Imn);m=M/mpeak;m_amud=m;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%invC=(NORMILUTZ′*NORMILUTZ) + 0.5*eye(imagewidth);invC=(ALLILUTZ′*ALLILUTZ) + 0.5*eye(imagewidth);C=inv(invC);%NORMILUTZ=NETILUTZ/235;%Sr=m_amud′*((NORMILUTZ′*NORMILUTZ)*m_amud);Sr=m_amud′*(invC*m_amud);%LOCENERG=NORMILUTZ.^2;LOCENERG=ALLILUTZ.^2;%locenerg=sqrt(sum(LOCENERG′))′;locenerg=(sum(LOCENERG′))′;n_keilu=sum(locenerg);dof=length(m_amud)-1;syx=sqrt(Sr/(n_keilu-dof));syx2=syx^2;C_norm=C*syx2;sig_m=sqrt(diag(C_norm));figure(5)plot(log10(sig_m(1imagewidth)./m_amud(mnnmxx)))grid on%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tocclear Vmlog=log(abs(m));%P=mlog′;P=zeros(1,oldimagewidth);P(mnnmxx)=mlog;leftmarg=1mnn-1;P(fliplr(leftmarg))=2*mlog(1)-mlog(leftmarg+1);rightmarg=mxx+1oldimagewidth;oldrightmarg=rightmarg-oldimagewidth+imagewidth;P(fliplr(rightmarg))=2*mlog(imagewidth)-mlog(oldrightmarg-1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%L=length(P);Imain=ones(L,1);if mnn>1 Imain(1mnn)=1;endif mxx<L Imain(mxx+1L)=1;endImat=sparse(diag(Imain));amain=ones(L,1);Amain=diag(amain);aoff1=-0.5*ones(L-1,1);aoff1(1)=-1;Aoff1=diag(aoff1,1);aoff_1=-0.5*ones(L-1,1);aoff_1(L-1)=-1;Aoff_1=diag(aoff_1,-1);S=sparse(Amain+Aoff1+Aoff_1);S(1,)=0;S(L,)=0;clear Amain Aoff1 Aoff_1 aoff_1 aoff1 amainStS=S′*S;clear S;R=Imat + lambdalog*StS;iR=inv(R);Y=iR*P′;Ynorm=Y-max(Y);smoothm=exp(Ynorm′);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mlogdelta=(sig_m(1imagewidth)./m_amud(mnnmxx));%P=mlog′;dP=zeros(1,oldimagewidth);dP(mnnmxx)=mlogdelta;leftmarg=1mnn-1;dP(fliplr(leftmarg))=2*mlogdelta(1)-mlogdelta(leftmarg+1);rightmarg=mxx+1oldimagewidth;oldrightmarg=rightmarg-oldimagewidth+imagewidth;dP(fliplr(rightmarg))=2*mlogdelta(imagewidth)-mlogdelta(oldrightmarg-1);dY=iR*dP′;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%longm=exp(P);longdm=smoothm.*dY′; m_lemata=smoothm-longdm;m_lemala=smoothm+longdm;figure(1)hold offplot(longm,′b′)hold onplot(smoothm,′r′)plot(m_lemata,′m′)plot(m_lemala,′m′)set(gca,′xlim′,[1,oldimagewidth]);set(gca,′ylim′,
);dlog2m=dY′/log(2);log2m_lemata=log2(smoothm)-dlog2m;log2m_lemala=log2(smoothm)+dlog2m;figure(2)hold offplot(log2(abs(longm)),′b′)hold onplot(log2(smoothm),′r′)plot(log2m_lemata,′m′)plot(log2m_lemala,′m′)set(gca,′xlim′,[1,oldimagewidth]);m=smoothm;delta_m=longdm;save randecay37.mat xint m delta_m Ynorm smoothm C ILUTZ sig_m mnn mxxload ../decayfunorigm=valint/max(valint);figure(1)hold onplot(origm,′g′)set(gca,′xlim′,[1,oldimagewidth]);figure(2)hold onplot(log2(origm),′g′)set(gca,′xlim′,[1,oldimagewidth]);clearclose allpackwidth=101;half_width=(width-1)/2;load nfuseupdated%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%logfused=log(1+fused);mx=max(logfused);imx16=(2^16)/mx;logfused16=uint16(round(logfused*imx16));histlog=double(histeq(logfused16,256));mn=min(histlog);histlog=histlog-mn;mx=max(histlog);imx8=(2^8)/double(mx);histlog8=uint8(round(double(histlog)*imx8));[rows,cols]=size(logfused16);L=1;R=cols;D=1;U=rows;figure(2);hold offimagesc(histlog8);colormap(gray(256));axis(′image′,′off′,′xy′);Xwide=cols/1.8;Ywide=round(Xwide*rows/cols);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*rows/cols;set(gcf,′position′,[1 500 Xwide Ywide]);clear deltaz2 logfusedbakasha=[′\n PRESS ANY KEY BETWEEN ORDERS′...  ′\n Press the letter ″n″ for normalized histogram′ ...   ′\n Press the letter ″e″ for normalized histogram′...   ′\n Press the digits ″s″,″m″,″l″ for smal//medium/large/windows \n\n′];button=0;sprintf(bakasha)pauseplotedimg=histlog8;figure(2)koteretl=sprintf(′Click LEFT mouse to mark area. Window width is %d′,width);set(gcf,′name′,koteretl)  while ~(button==113),%′q′= botton=113   figure(2)   [x,y,button]=ginput(1);   hold on   if(button==115)%′s′= botton=115  width=41;   koteretl=sprintf(′Click LEFT mouse to mark area. Window width is %d′, width);   set(gcf,′name′,koteretl)   elseif (button==109) %′m′= botton=109   width=75;  koteretl=sprintf(′Click LEFT mouse to mark area. Window width is %d′,width);  set(gcf,′name′,koteretl)   elseif (button==108) %′1′= botton=108   width=101;   koteretl=sprintf(′Click LEFT mouse to mark area. Window width is %d′,width);   set(gcf,′name′,koteretl)   elseif (button==1) | (button==3)   half_width=(width-1)/2;   XL=round(x-half_width);XR=round(x+half_width);   YD=round(y-half_width);YU=round(y+half_width);   XLL=max(XL,L); XRR=min(XR,R);YDD=max(YD,D);YUU=min(YU,U);   X=[xLL XLL XRR XRR]; Y=[YDD YUU YUU YDD];   plotedimg=plotcontrs5(X,Y,plotedimg,histlog8);   figure(2);   hold off   imagesc(plotedimg);   colormap(gray(256)); axis(′image′,′off′,′xy′); %Xwide=cols/1.8; Ywide=round(Xwide*rows/cols); %set(gcf,′position′,[55 109 Xwide Ywide]); %set(gca,′position′,
); %vec=get(gcf,′position′); %Ywide=vec(3)*rows/cols; %set(gcf,′position′,[1 500 Xwide Ywide]);   axis(′image′,′off′,′xy′); set(gca,′position′,
); % ribuax=[X, X(1)]; ribuay=[Y, Y(1)]; % hold on;plot(ribuax,ribuay,′m-′)   [enhanced,subenhanced,mnval,mxval]=stretchwin5(X,Y,fused,histlog8);   smunval=sprintf(′%d′,round(mnval))   smxval=sprintf(′%d′,round(mxval))   figure   imagesc(subenhanced);   colormap(gray(256))   axis(′image′,′off′,′xy′);set(gca,′position′,
);   set(gcf,′position′,[100 100 2*width 2*width]);   end   pause(0.5)   end %tofileimg=flipud(plotedimg); %to_mx=max(tofileimg); %tofileimg=double(tofileimg)/double(to_mx);%imwrite(tofileimg,′plotedmos1.jpg′,′jpeg′,′Quality′,100);function plotedimg=plotcontrs5(X,Y,oldplotimg,logphoto)XS=sort (X);YS=sort (Y);L=round (XS (2));R=round (XS (3));D=round (YS (2));U=round (YS (3));% LD LU RU RD LDXp=[L L R R L];Yp=[D U U D D];X0=(L+R)/2;Y0=(D+U)/2;spanx=size(oldplotimg,2);spany=size(oldplotimg,1);[LD_row,LD_col]=xy2ij_brack([L,D],spany);[RD_row,RD_col]=xy2ij_brack([R,D],spany);[LU_row,LU_col]=xy2ij_brack([L,U],spany);[RU_row,RU_col]=xy2ij_brack([R,U],spany);stam=flipud(logphoto);smol=stam(LD_row-1LU_row ,LD_col);newsmol=uint8(255*double(smol<128));yamin=stam(RD_row-1RU_row,RD_col);newyamin=uint8(255*double(yamin<128));lemata=stam(RD_row ,LD_colRD_col);newlemata=uint8(255*double(lemata<128));lemala=stam(RU_row ,LU_colRU_col);newlemala=uint8(255*double(lemala<128));%whos newsmol newyamin newlemata newlemala%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%stam=flipud(oldplotimg);stam( LD_row-1LU_row ,LD_col)=newsmol;stam( RD_row-1RU_row ,RD_col)=newyamin;stam( RD_row ,LD_colRD_col)=newlemata;stam( RU_row ,LU_colRU_col)=newlemala;plotedimg=flipud(stam);function prof=profil2(cols,transit)mn=eps^7;tmp=ones(1,cols);L=ltransit;R=cols-L+1;%tmp(L)=mn+ (1-mn)*(L-1)/(transit-1);%tmp(R)=mn+ (1-mn)*(L-1)/(transit-1);tmp(L)=mn+ (1-mn) * ( 1 - cos(pi*(L-1)/(transit-1)) )/2;tmp(R)=mn+ (1-mn) * ( 1 - cos(pi*(L-1)/(transit-1)) )/2;prof=tmp;function {oldselect,newselect]=ptspick6(oldcell,newcell,greatL,greatR,greatD,greatU)sun=240;%sun=200;%sun=235;dark=2;kamut=50;segments=20;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold=oldcell{1,1}; % the 1st (old) imageoldcenterx=oldcell{2,1};%coordinates of center of old imageoldcentery=oldcell{3,1};[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old imageZnew=newcell{1,1}; % the 2nd (new) imagenewcenterx=newcell{2,1};%coordinates of center of new imagenewcentery=newcell{3,1};[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; % the Left side of the new imageRnew=newcenterx+(newcols-1)/2; % the Right sideDnew=newcentery-(newrows-1)/2; % the Down of the new imageUnew=newcentery+(newrows-1)/2; % the Upper side of the new imagespanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy, oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy, oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy, oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy, oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);{oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;Zoldgreat=flipud(Zoldgreat);%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find ((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find ((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find ((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat( newLD_row-1newLU_row,newLD_colnewRD_col )=Znew;Znewgreat=flipud(Znewgreat);%%%%%%%%%%%%%%% DIVIDING TO PARTSoldLOnan= (~isnan(Zoldgreat));newLOnan= (~isnan(Znewgreat));LOnan= (oldLOnan newLOnan);oversegm=max(LOnan);I=find(oversegm);shortL=min(I);shortR=max(I);deltax=(shortR-shortL)/segments;%%%%%%%%%%%%%%%% SELECTING GOOD POINTSminlevel=dark;maxlevel=sun;oldinrange=( (Zoldgreat>minlevel) (~(Zoldgreat>maxlevel)) );newinrange=( (Znewgreat>minlevel) (~(Znewgreat>maxlevel)) );inrange= ( (oldinrange newinrange) LOnan);[rows,cols]=size(inrange);greatCOLS=ones(rows,1)*(1cols);%figure(1);%hold off%imagesc(inrange);%colormap(gray(256)); axis(′image′,′off′,′xy′);%Xwide=cols/1.8;Ywide=round(Xwide*rows/cols);%set(gcf,′position′,[1 500 Xwide Ywide]);%set(gca,′position′,
);%vec=get(gcf,′position′);%Ywide=vec(3)*rows/cols;%set(gcf,′position′, [1 500 Xwide Ywide]);%pausegoodpts=[ ];for ypart=1segments,   Lsegm= shortL + round(deltax*(ypart-1));   Rsegm= shortL + round(deltax*(ypart));%insegment= ( (~(greatX < Lsegm)) (~(greatX > Rsegm)) );   insegment= ( (~(greatCOLS < Lsegm)) (~(greatCOLS > Rsegm)) );   beseder= (inrange insegment);%figure(1);%hold off%imagesc(beseder);%colormap(gray(256));axis(′image′,′off′,′xy′);%[rows,cols]=size(inrange);%Xwide=cols/1.8; Ywide=round(Xwide*rows/cols);%set(gcf,′position′,[1 500 Xwide Ywide]);%set(gca,′position′,
);%vec=get(gcf,′position′);%Ywide=vec(3)*rows/cols;%set(gcf,′position′, [1 500 Xwide Ywide]);%pause inrange_ind=find(beseder); lenrange=length(inrange_ind); if lenrange<kamut   goodpts=[goodpts;inrange_ind]; else   linerand=1+round((lenrange-1)*(rand(kamut,1)));   picked_ind=inrange_ind(linerand);   goodpts=[goodpts;picked_ind]; endend%%%%%%%%% THE POINT COORDS VALSx=greatX(goodpts);y=greatY(goodpts);oldxsmall=x-Lold+1;oldysmall=y-Dold+1;newxsmall=x-Lnew+1;newysmall=y-Dnew+1;oldvals=Zoldgreat(goodpts);newvals=Znewgreat(goodpts);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%oldselect{1,1}=oldxsmall;oldselect{2,1}=oldysmall;oldselect{3,1′}=oldvals;newselect{1,1}=newxsmall;newselect{2,1}=newysmall;newselect{3,1}=newvals;function [oldgreatcell,newgreatcell]=putbigml3(oldcell,newcell)%% oldcell and newcell contain images to put in a single% coordinate system, and the vector-coordinates of the% image centers in (x,y).% oldgreatcell and newgreatcell are the cells of the new large% frame with a common coorinate system,in each of which% the single images are embedded. These cells also contain the% coordinates of the center of the big frames.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold=oldcell{1,1}; % the 1st (old) imageDELTAZ2old=oldcell{2,1};oldcenterx=oldcell{3,1}; %coordinates of center of old imageoldcentery=oldcell{4,1};Znew=newcell{1,1};% the 2nd (new) imageDELTAZ2new=newcell{2,1};newcenterx=newcell{3,1};newcentery=newcell{4,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; % the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; %the Left side of the new (2nd) imageRnew=newcenterx+(newcols-1)/2;Dnew=newcentery-(newrows-1)/2;Unew=newcentery+(newrows-1)/2;greatL=min([Lold,Lnew]); %the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat( oldLD_row-1oldLU_row,oldLD_cololdRD_col )=Zold;Zoldgreat=flipud(Zoldgreat);DELTAZoldgreat=NaN*zeros(size(greatX));DELTAZoldgreat( oldLD_row-1oldLU_row,oldLD_cololdRD_col )=DELTAZ2old;DELTAZoldgreat=flipud(DELTAZoldgreat);%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find ((grsatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find ((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find ((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat( newLD_row-1newLU_row,newLD_colnewRD_col )=Znew;Znewgreat=flipud(Znewgreat);DELTAZnewgreat=NaN*zeros(size(greatX));DELTAZnewgreat( newLD_row-1newLU_row,newLD_colnewRD_col )=DELTAZ2new;DELTAZnewgreat=flipud(DELTAZnewgreat);%%%%%%%%%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT-COMBINED IMAGEgreatcenterx=(greatR+greatL)/2;greatcentery=(greatU+greatD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERoldgreatcell{1,1}=Zoldgreat;oldgreatcell{2,1}=DELTAZoldgreat;oldgreatcell{3,1}=greatcenterx;oldgreatcell{4,1}=greatcentery;newgreatcell{1,1}=Znewgreat;newgreatcell{2,1}=DELTAZnewgreat;newgreatcell{3,1}=greatcenterx;newgreatcell{4,1}=greatcentery;function [oldgreatcell,newgreatcell]=puthullmll(oldcell,newcell)%% oldcell and newcell contain images to put in a single% coordinate system, and the vector-coordinates of the% image centers in (x,y).% oldgreatcell and newgreatcell are the cells of the new large% frame with a common coorinate system,in each of which% the single images are embedded. These cells also contain the% coordinates of the center of the big frames.ymargin=100;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold =oldcell{1,1}; % the 1st (old) imageDELTAZ2old=oldcell{2,1};oldcenterx=oldcell{3,1};%coordinates of center of old imageoldcentery=oldcell{4,1};Znew=newcell{1,1}; % the 2nd (new) imageDELTAZ2new=newcell{2,1};newcenterx=newcell{3,1};newcentery=newcell{4,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=ceil(oldcenterx-(oldcols-1)/2)-0.5; % the Left side of the old imageRold=ceil(oldcenterx+(oldcols-1)/2)-0.5; % the Right sideDold=ceil(oldcentery-(oldrows-1)/2)-0.5; % the Down of the old imageUold=ceil(oldcentery+(oldrows-1)/2)-0.5; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=ceil(newcenterx-(newcols-1)/2)-0.5; % the Left side of the new (2nd) imageRnew=ceil(newcenterx+(newcols-1)/2)-0.5; Dnew=ceil(newcentery-(newrows-1)/2)-0.5;Unew=ceil(newcentery+(newrows-1)/2)-0.5;greatL=min([Lold,Lnew]); %the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;DELTAZoldgreat=NaN*zeros(size(greatX));DELTAZoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=DELTAZ2old;%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find ((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find ((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find ((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[nawRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;DELTAZnewgreat=NaN*zeros(size(greatX));DELTAZnewgreat(newLD_row-lnewLU_row,newLD_colnewRD_col)=DELTAZ2new;%%%%%%%%%%%%%%%%%%%% %%%%%%%% DEFINING HULLLnew_margin=Lnew-ymargin;Rnew_margin=Rnew+ymargin;Dnew_margin=Dnew-ymargin;Unew_margin=Unew+ymargin;hullL=min([ max([Lold,Lnew_margin],Lnew]);%the Left side of the combined-greaterimagehullR=max([ min([Rold,Rnew_margin],Rnew]);hullD=min([ max([Dold,Dnew_margin],Dnew]);hullU=max([ min([Uold,Unew_margin],Unew]);[hullLDy,hullLDx]=find ((greatX==hullL) (greatY==hullD));[hullLUy,hullLUx]=find ((greatX==hullL) (greatY==hullU));[hullRDy,hullRDx]=find ((greatX==hullR) (greatY==hullD));[hullRUy,hullRUx]=find ((greatX==hullR) (greatY==hullU));[hullLD_row, hullLD_col]=xy2ij_brack([hullLDx,hullLDy],spany);[hul1RD_row,hullRD_col]=xy2ij_brack([hullRDx,hullRDy],spany);[hullLU_row,hullLU_col]=xy2ij_brack([hullLUx,hullLUy],spany);[hullRU_row,hullRU_col]=xy2ij_brack([hullRUx,hullRUy],spany);Zoldhull=Zoldgreat( hullLD_row-1hullLU_row,hullLD_colhullRD_col);Znewhull=Znewgreat( hullLD_row-1hullLU_row,hullLD_colhullRD_col);DELTAZ2oldhull=DELTAZoldgreat( hullLD_row-1hullLU_row ,hullLD_colhullRD_col);DELTAZ2newhull=DELTAZnewgreat( hullLD_row-1hullLU_row ,hullLD_colhullRD_col);%%%%%%%%%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT-COMBINED IMAGE%figure(1)%imagesc(Zoldhull);axis xy%figure(2)%imagesc(Znewhull);axis xy%Zoldhull=flipud(Zoldhull);%Znewhull=flipud(Znewhull);%figure(3)%imagesc(Zoldhull);axis xy;colormap(gray(256))%figure(4)%imagesc(Znewhull);axis xy;colormap(gray(256))%size(Zoldhull)hullcenterx=(hullR+hullL)/2;hullcentery=(hullU+hullD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERoldgreatcell{1,1}=Zoldhull;oldgreatcell{2,1}=DELTAZ2oldhull;oldgreatcell{3,1}=hullcenterx;oldgreatcell{4,1}=hullcentery;newgreatcell{1,1}=Znewhull;newgreatcell{2,1}=DELTAZ2newhull;newgreatcell{3,1}=hullcenterx;newgreatcell{4,1}=hullcentery;function [oldgreatcell,newgreatcell]=putinbig6(oldcell,newcell)%% oldcell and newcell contain images to put in a single% coordinate system,and the vector-coordinates of the% image centers in (x,y).% oldgreatcell and newgreatcell are the cells of the new large% frame with a common coorinate system,in each of which% the single images are embedded.These cells also contain the% coordinates of the center of the big frames.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold=oldcell{1,1}; % the 1st (old) imageoldcenterx=oldcell{2,1};%coordinates of center of old imageoldcentery=oldcell{3,1};Znew=newcell[1,1}; % the 2nd (new) imagenewcenterx=newcell{2,1};newcentery=newcell{3,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; %the Left side of the new (2nd) imageRnew=newcenterx+(newcols-1)/2;Dnew=newcentery-(newrows-1)/2;Unew=newcentery+(newrows-1)/2;greatL=min([Lold,Lnew]); %the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR);%the coordinate grid of the combined-greater imagegreatY=(greatDgreatU) ′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;Zoldgreat=flipud(Zoldgreat);%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find ((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find ((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find ((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;Znewgreat=flipud(Znewgreat);%%%%%%%%%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT-COMBINED IMAGEgreatcenterx=(greatR+greatL)/2;greatcentery=(greatU+greatD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERoldgreatcell{1,1}=Zoldgreat;oldgreatcell{2,1}=greatcenterx;oldgreatcell{3,1}=greatcentery;newgreatcell{1,1}=Znewgreat;newgreatcell{2,1}=greatcenterx;newgreatcell{3,1}=greatcentery;function [small,smalldelta2]=reduceml21(big,bigdelta2,shape);gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25;gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=1./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(conv2(weightedbig,gs,shape),gs′,shape); %convolution in one lineblurdenomin=conv2(conv2(invsig2,gs,shape),gs′,shape);samples_row=1 2 size(blurnominat,1); %sample grid in one linesamples_col=1 2 size(blurnominat,2); %sample grid in one linesmallnominat=blurnominat(samples_row,samples_col);%sampling in one linesmalldenomin=blurdenomin(samples_row,samples_col);%sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig2;function [small,smalldelta2]=reduceml21_x(big,bigdelta2,shape)gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25; gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=l./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(weightedbig,gs,shape); %convolution in one lineblurdenomin=conv2(invsig2,gs,shape);samples_col=1 2 size(blurnominat,2); %sample grid in one linesmallnominat=blurnominat(,samples_col); %sampling in one linesmalldenomin=blurdenomin(,samples_col); %sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig2=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig2;function [small,smalldelta2]=reduceml21_x(big,bigdelta2,shape)gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25;gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=1./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(weightedbig,gs′,shape); %convolution in one lineblurdenomin=conv2(invsig2,gs′,shape);samples_row=1 2 size(blurnominat, 1); %sample grid in one linesmallnominat=blurnominat(samples_row,); %sampling in one linesmalldenomin=blurdenomin(samples_row,); %sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig2=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig;function Fethsatur=satursmooth7(X,uppersun,lowersun)SE=ones(3);mn=eps;saturated= uint8(X>uppersun);maybesatur=uint8(X>lowersun);s=max (saturated );if s==0 Fethsatur=ones(size (X));else shulaim=uint8(1); oldregion=uint8(saturated); %k=0; while ~(shulaim==0)   %k=k+1   dilsaturated=(dilate(oldregion,SE));   newregion=(dilsaturated maybesatur);   added= newregion (~(oldregion));   shulaim=max(added);   oldregion=newregion; endtmp=ones (size (X)); fact=(1-mn)/(uppersun- (lowersun-1)); for level=lowersunuppersun,   badpixels=((X==level) newregion);   badvalue=mn+ (uppersun-level) * fact;   lobeseder=find (badpixels);   if ~isempty (lobeseder)   tmp (lobeseder) =badvalue;   end  end  lobeseder=find(saturated);  tmp(lobeseder)=mn;  Fethsatur=tmp;endfunction [enhanced,subenhanced,mnval,mxval]=stretchwin5(X,Y,photo,logphoto)XS=sort(X);YS=sort(Y);L=round(XS(2));R=round(XS(3));D=round(YS(2));U=round(YS(3));% LD LU RU RD LDXp=[L L R R L];Yp=[D U U D D];X0=(L+R)/2;Y0=(D+U)/2;spanx=size(photo,2);spany=size(photo,1);[LD_row,LD_col]=xy2ij_brack([L,D],spany);[RD_row,RD_col]=xy2ij_brack([R,D],spany);[LU_row,LU_col]=xy2ij_brack([L,U],spany);[RU_row,RU_col]=xy2ij_brack([R,U],spany);stam=flipud(photo);subarea=stam(LD_row-1LU_row,LD_colRD_col);mn=min(subarea);lemata=subarea-mn;mx=max(lemata);imx8=(2^8)/mx;subenhanced=uint8(round(lemata*imx8));mn=min(logphoto);lemata=double(logphoto)-double(mn);mx=max(lemata);imx8=(2^8)/mx;stam=flipud(lemata*imx8);stam(LD_row-1LU_row,LD_colRD_col)=subenhanced;enhanced=flipud(stam);mnval=min(subarea);mxval=max(subarea);function[beseder,xbeseder,ybeseder,valbeseder]=subpts2l(oldcell,greatL,greatR,greatD,greatu,goodpts)%sun=240;%sun=200;sun=235;dark=2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold=oldcell{1,1}; % the 1st (old) imageoldcenterx=oldcell{2,1};%coordinates of center of old imageoldcentery=oldcell{3,1};[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old imagespanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;Zoldgreat=flipud(Zoldgreat); %%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRIDvals=Zoldgreat(goodpts);LOnan=(-isnan(vals));LObright =(vals<sun);LOdark =(~(vals<dark));beseder= ((LOnan LObright) LOdark);%beseder= (LOnan);beseder_ind=find(beseder);besederpts=goodpts(beseder_ind);x=greatX(besederpts);y=greatY(besederpts);xsmall=x-Lold+1;ysmall=y-Dold+1;xbeseder=NaN*ones(size(goodpts));xbeseder(beseder_ind)=xsmall;ybeseder=NaN*ones(size(goodpts));ybeseder(beseder_ind)=ysmall;valbeseder=NaN*ones(size(goodpts));valbeseder(beseder_ind)=Zoldgreat(besederpts);%spanx=Rold-Lold+1;%spany=Uold-Dold+1;%greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the comoined-greater image%greatY=(greatDgreatU) ′*ones(1,spanx);clearclose allpack%load ../decayfun%clear m delta_mload decay53337pairsclear ILUTZ C Ynorm sig_m smoothm mnn mxx%xint m delta_m Ynorm smoothm C ILUTZ sig_m mnn mxx%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fuzzmarge=45;transit=80;uppersun=250;lowersun=205;levels_x=2;levels_y=0;minlevel=0;Lshulaim=0;left=min(xint)+Lshulaim;right=max(xint);bottom=170;roof=355;imagehight=roof-bottom+1;SE=ones(3,5);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cols=right-left+1;prof=profil2(cols,transit);inv_prof2=(1./prof).^2;FETHER2=ones(imagehight,1)*inv_prof2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%oldlength=length(m);m=m(1+Lshulaimoldlength);delta m=delta_m(1+Lshulaimoldlength);sig=1;invm=1./m;invm_square=invm.^2;invm_four =invm.^4;a=(sig.^2).*invm_square;b=(delta_m.^2).*invm_four;invM=ones(imagehight,1)*invm;A= ones(imagehight,1)*a;B= ones(imagehight,1)*b;load ../coordml10.mat %%%% LOADING THE COORDINATESbasefilename=′../gradfilt1_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[1∶9 , 11∶37];centerx1=0; %this is alwayscentery1=0; %this is always%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MERGING IMAGES IN A BIG FRAMEstam=round(rand(1));ystep=(2*stam)-1;if ystep==-1 frame_inds=36ystep2;elseif ystep==1 frame_inds=1ystep35;end%ystep=-1;%frame_inds=36ystep2;k=1;for indframe=1length(frame_inds),%tic yframe_ind=frame_inds(indframe); if k==1   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 postfix];rawim=imread(shem1);   basicim=double(flipud(rawim))-1;   X=basicim(bottomroof,leftright);   Ya=X.*invM;   Fethsatur=satursmooth7(X,uppersun,lowersun);   IFETHSATUR2=(1./Fethsatur).^2;   DELTAYa2=(A+ ((X.^2).*B)).*IFETHSATUR2;   DELTAYa2=DELTAYa2.*FETHER2;   Acell[1,1}=Ya;   Acell{2,1}=DELTAYa2;   Acell{3,1}=centerx2(yframe_ind);   Acell{4,1}=centery2(yframe_ind);  else   Acell=fusedcell;  end yframe2=netframes(yframe_ind + ystep); sframe2=sprintf(′%d′,yframe2); shem2=[basefilename sframe2 postfix]; rawim=imread(shem2); basicim=double(flipud(rawim))-1; X=basicim(bottomroof,leftright); Yb=X.*invM; Fethsatur=satursmooth7(X,uppersun,lowersun); IFETHSATUR2=(1./Fethsatur).^2; DELTAYb2=(A+ ((X.^2).*B)).*IFETHSATUR2; DELTAYb2=DELTAYb2.*FETHER2; clear rawim basicim X Fethsatur IFETHSATUR2 Bcell{1,1}=Yb; Bcell{2,1}=DELTAYb2; Bcell{3,1}=centerx2(yframe_ind+ystep); Bcell{4,1}=centery2(yframe_ind+ystep); fusedcell=fuseml22(Acell,Bcell);  k=k+1;  clear Acel1 Bcell%tocendfused=double(fusedcell{1,1});deltaz2=fusedcell{2,1};greatcenterx=fusedcell{3,1};greatcentery=fusedcell{4,1};save nfuseupdated fused deltaz2 greatcenterx greatcenteryclear fusedcellclear saturated logfused valint tikunx tikuny invM invm invm_squareclear invm_square dilsaturated flatcurve flatfun centery2 centerx2clear a b Ya Yb DELTAYb2 DELTAYa2 A Bmax(fused)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%load nfuseupdatedlogfused=log(1+fused);mx=max(logfused);%imx8=(2^8)/mx;imx16=(2^16)/mx;logfused16=uint16(round(logfused*imx16));histlog=histeq(logfused16,256);mx=max(histlog);%imx8=(2^8)/double(mx);imx16=(2^16)/double(mx);histlog16=uint16(round(double(histlog)*imx16));[rows,cols]=size(logfused16);figure(2);hold offimagesc(histlog16);colormap(gray(256));axis(′image′,′off′,′xy′);Xwide=cols/1.8;Ywide=round(Xwide*rows/cols);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*rows/cols;set(gcf,′position′,[1 300 Xwide Ywide]);% [I,J]=xy2ij_brack(XYcoord,N);% Converts x-y coordinates to ij coordinates% I is the row-number of the element,% J is the column-number ofthe element% N is the number of rows in the matrix.function [I,J]=xy2ij_brack(XYcoord,N);X=XYcoord(,1);Y=XYcoord(,2);%index in row-numberI=N-Y+1;%index in colsJ=X;%IJcoord=[I,J];形成偏振拼合的原始碼clearclose allpackload fusecomponl[ROWS,COLS]=size(Lfused);I0 =reshape(Lfused,ROWS*COLS,1);I45=reshape(Mfused,ROWS*COLS,1);I90=reshape(Rfused,ROWS*COLS,1);F=[I0,I45,I90]′;alfa=
′*pi/180;M=[ones(3,1),cos(2*alfa),sin(2*alfa)];IM=inv(M);Param=IM*F%%%%%%%%%%%%%%%%%%%%%%%%%%%% finding the net size%%%%%%%%%%%%%bigC= Param(1,);bigC=reshape(bigC,ROWS,COLS);stamC=sum(double(uint8(bigC)));a=find(stamC>0);l=min(a);r=max(a);stamC=sum(double(uint8(bigC)),2);a=find (stamC>0);d=min (a);u=max (a);lfused=Lfused(du,lr);mfused=Mfused(du,lr);rfused=Rfused(du,lr);[rows,cols]=size(lfused);clear bigC Lfused Mfused Rfused a stamC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%I0 =reshape(lfused,rows*cols,1);I45=reshape(mfused,rows*cols,1);I90=reshape(rfused,rows*cols,1);F=[I0,I45,I90]′;alfa=
′*pi/180;M=[ones(3,1),cos(2*alfa),sin(2*alfa)];IM=inv(M);Param=IM*F;C= Param(1,);Ac=Param(2,);As=Param(3,);clear Param%%%%%%%%%%%%%%%%%%%%%%%%% A positive,ambiguous theta%%%%%%%%%%%%%%%%%%%Twotheta=atan(As./Ac).*(~(Ac==0))-pi/2*((Ac==0).*(As<0))+pi/2*((Ac==0).*(~(As<0)));twotheta=Twotheta -((Ac<0).*(As<0)*pi) +((Ac<0).*(~(As<0))*pi);theta=(twotheta/2);A=sqrt(Ac.^2 + As.^2);theta=theta*180/pi;%C=reshape(C,ROWS,COLS);%A=reshape(A,ROWS,COLs);%THETA=reshape(theta,ROWS,COLS);C=reshape(C,rows,cols);A=reshape(A,rows,cols);THETA=reshape(theta,rows,cols);clear Twotheta twotheta theta%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5ROWS=rows;COLS=cols;figure(1);hold offimagesc(C);axis(′image′,′off′,′xy′);%Xwide=COLS/0.7;Xwide=COLS/0.9;Ywide=round(Xwide*ROWS/COLS);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*ROWS/COLS;set(gcf,′position′,[10 300 Xwide Ywide]);colormap(gray(256)AA=(C-A);b=AA.*(AA>0);figure(2);hold offimagesc(b);axis(′image′,′off′,′xy′);%Xwide=COLS/0.7;Xwide=COLS/0.9;Ywide=round(Xwide*ROWS/COLS);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*ROWS/COLS;set(gcf,′position′,[10 300 Xwide Ywide]);colormap(gray(256))figure(3);hold offimagesc(A);axis(′image′,′off′,′xy′);%Xwide=COLS/0.7;Xwide=COLS/0.9;Ywide=round(Xwide*ROWS/COLS);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*ROWS/COLS;set(gcf,′position′,[10 300 Xwide Ywide]);colormap(gray(256))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear A Ac As C F I0 I45 I90 THETA AA bclearclose allpackload polarresl[ROWS,COLS]=size(C);c=reshape(C,1,ROWS*COLS);a=reshape(A,1,ROWS*COLS);ac=[c;a]boundeda=min(ac);boundedA=reshape(boundeda,ROWS,COLS);figure(1);   hold off   imagesc(C);   imagesc(A);   axis(′image′,′off′,′xy′);   %Xwide=COLS/0.7;   Xwide=COLS/0.9;   Ywide=round(Xwide*ROWS/COLS);   set(gcf,′position′,[55 109 Xwide Ywide]);   set(gca,′position′,
);   vec=get(gcf,′position′);   Ywide=vec(3)*ROWS/COLS;   set(gcf,′position′,[10 300 Xwide Ywide]);   colormap(gray(256))   figure(3+yiter);   hold offaxis(′image′,′off′,′xy′);   %Xwide=COLS/0.7;   Xwide=COLS/0.9;   Ywide=round(Xwide*ROWS/COLS);   set(gcf,′position′,[55 109 Xwide Ywide]);   set(gca,′position′,
);   vec=get(gcf,′position′);    Ywide=vec(3)*ROWS/COLS;   set(gcf,′position′,[10 300 Xwide Ywide]);   colormap(gray (256))   AA= (C-A);   b=AA.*(AA>0);   figure (4+yiter);   hold off   imagesc(b);   axis(′image′,′off′,′xy′);   %Xwide=COLS/0.7;   Xwide=COLS/0.9;   Ywide=round(Xwide*ROWS/COLS);   set(gcf,′position′,[55 109 Xwide Ywide]);   set(gca,′position′,
);   vec=get(gcf,′position′);   Ywide=vec(3)*ROWS/COLS;   set(gcf,′position′,[10 300 Xwide Ywide]);   colormap(gray(256))clearclose allpackload flatparams.matclear fused8 xlims ylims cparams%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%transit=40;uppersun=250;lowersun=205;left=97;right=600;bottom=200;roof=385;imagehight=roof-bottom+1;imagewidth=right-left+1;LeftspaceL=135;LeftspaceR=185;RightspaceL=315;RightspaceR=365;Leftspacewidth= LeftspaceR-LeftspaceL+1;Rightspacewidth=RightspaceR-RightspaceL+1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sig=0.5;invM=flatfun(bottomroof,leftright);SIG2=(sig.^2);cols=right-left+1;prof=profil2(cols,transit);inv_prof2=(1./prof).^2;FETHER2=ones(imagehight,1)*inv_prof2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%load coordpoll.mat %%%% LOADING THE COORDINATESbasefilename=′stepol_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[115];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MERGING IMAGES IN A BIG FRAMEstam=1;ystep=(2*stam)-1;frame_inds=1ystep(length(netframes)-1);k=1;for indframe=1length (frame_inds),%tic yframe_ind=frame_inds(indframe) if k=1   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 ′_1′ postfix];   shem2=[basefilename sframe1 ′_2′ postfix];   shem3=[basefilename sframe1 ′_3′ postfix];   rawim1=imread(shem1);   rawim2=imread(shem2);   rawim3=imread(shem3);   rawim=(double(flipud(rawim1))~1 ) + ( double(flipud(rawim2))-1 ) + (double(flipud(rawim3))-1);   basicim=(rawim/3);   Za=basicim.*flatfun;   Za=Za(bottomroof,leftright);   AL=NaN*ones(size(Za));AM=NaN*ones(size(Za));AR=NaN*ones(size(Za));   DELTAAL2=NaN*ones(size(Za)); DELTAAM2=NaN*ones(size(Za));DELTAAR2=NaN*ones(size(Za));    AL(,1LeftspaceL-1)=Za(,1LeftspaceL-1);   AM(,LeftspaceR+1RightspaceL-1)=Za(,LeftspaceR+1RightspaceL-1);   AR(,RightspaceR+1imagewidth)=Za(,RightspaceR+1imagewidth);   X=basicim(bottomroof,leftright);   Fethsatur=satursmooth7(X,uppersun,lowersun);   IFETHSATUR2=(1./Fethsatur).^2;   DELTAZa2=(SIG2*(invM.^2)).*IFETHSATUR2;% DELTAYa2=DELTAYa2.*FETHER2;% DELTAZa2=(SIG2*(invM.^2)).*(~dilsaturated) + sunerror.*dilsaturated;   DELTAAL2(,1LeftspaceL-1)=DELTAZa2(,1LeftspaceL-1);   DELTAAM2(,LeftspaceR+1RightspaceL~1)=DELTAZa2(,LeftspaceR+1RightspaceL-1);DELTAAR2(,RightspaceR+1imagewidth)=DELTAZa2(,RightspaceR+1imagewidth);   ALcell{1,1}=AL;   ALcell{2,1}=DELTAAL2;   ALcell{3,1}=centerx2(yframe_ind);   ALcell{4,1}=centery2(yframe_ind);   AMcell=ALcell;   AMcell{1,1}=AM;   AMcell{2,1}=DELTAAM2;   ARcell=ALcell;   ARcell{1,1}=AR;   ARcell{2,1}=DELTAAR2; else   ALcell=Lfusedcell;   AMcell=Mfusedcell;   ARcell=Rfusedcell; end yframe2=netframes(yframe_ind+ystep); sframe1=sprintf(′%d′,yframe2); shem1=[basefilename sframe1 ′_1′postfix]; shem2=[basefilename sframe1 ′_2′postfix]; shem3=[basefilename sframe1 ′_3′postfix]; rawim1=imread(shem1); rawim2=imread(shem2); rawim3=imread(shem3); rawim=(double(flipud(rawim1))-1 ) + ( double(flipud(rawim2))-1 ) + (double(flipud(rawim3))-1 ); basicim=(rawim/3); Zb=basicim.*flatfun; Zb=Zb(bottomroof,leftright); BL=NaN*ones(size(Zb));BM=NaN*ones(size(Zb));BR=NaN*ones(size(Zb)); DELTABL2=NaN*ones(size(Zb));DELTABM2=NaN*ones(size(Zb)); DELTABR2=NaN*ones(size(Zb)); BL(,1LeftspaceL-1)=Zb(,1LeftspaceL-1); BM(,LeftspaceR+1RightspaceL-1)=Zb(,LeftspaceR+1RightspaceL-1); BR(,RightspaceR+1imagewidth)=Zb(,RightspaceR+1imagewidth); X=basicim(bottomroof,leftright); Fethsatur=satursmooth7(X,uppersun,lowersun); IFETHSATUR2=(1./Fethsatur).^2; DELTAZb2=(SIG2*(invM.^2)).*IFETHSATUR2;% DELTAZb2=(SIG2*(invM.^2)).*(~dilsaturated) + sunerror.*dilsaturated;% DELTAZb2=DELTAZb2.*FETHER2; DELTABL2(,1LeftspaceL-1)=DELTAZb2(,1LeftspaceL-1); DELTABM2(,LeftspaceR+1RightspaceL~1)=DELTAZb2(,LeftspaceR+1RightspaceL-1); DELTABR2(,RightspaceR+1imagewidth)=DELTAZb2(,RightspaceR+1imagewidth); clear rawim basicim Fethsatur IFETHSATUR2BLcell{1,1}=BL; BLcell{2,1}=DELTABL2; BLcell{3,1}=centerx2{yframe_ind+ystep); BLcell{4,1}=centery2(yframe_ind+ystep); BMcell=BLcell; BMcell{1,1}=BM; BMcell{2,1}=DELTABM2; BRcell=BLcell; BRcell{1,1}=BR; BRcell{2,1}=DELTABR2; Lfusedcell=fuseml22(ALcell,BLcell); Mfusedcell=fuseml22 (AMcell,BMcell); Rfusedcell=fuseml22(ARcell,BRcell); k=k+1; clear ARcell BRcell AGcell BGcell ABcell BBcell%tocendLfused=double(Lfusedcell{1,1});Mfused=double(Mfusedcell{1,1});Rfused=double(Rfusedcell{1,1});greatcenterx=Lfusedcell{3,1};greatcentery=Lfusedcell{4,1};save fusecomponl Lfused Mfused Rfused greatcenterx greatcenteryclear Rfusedcell Gfusedcell Bfusedcellclear saturated logfused valint tikunx tikuny invM invm invm_squareclear invm_square dilsaturated flatcurve centery2 centerx2clear a b Ya Yb DELTAYb2 DELTAYa2 A B FETHER2 AB AG AR BB BG BRclear basefilename basephotoname greatcenterx greatcentery imagehightclear indframe inv_prof2 k lambda_pos left netframes postfix profclear rawname right roof sframel sframe2 shem1 shem2 stam transitclear yframel yframe2 yframe_ind ystep%clear bX CAMERESP invCAMERESP invCAMRESframe X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%yGAMMA=0.5;load fusecomponl[ROWS,COLS]=size(Lfused);fused(1ROWS,1COLS,1)=Lfused;fused(1ROWS,1COLS,2)=Mfused;fused(1ROWS,1COLS,3)=Rfused;mx=max(fused);fused=fused/mx;oldHSV=rgb2hsv(fused);luminance=oidHSV(,,3);mn=min(luminance);mx=max(luminance);newlumin=imadjust(luminance,[mn mx],
,yGAMMA);newHSV=oldHSV;newHSV(,,3)=newlumin;newfused=hsv2rgb(newHSV);mx=max(newfused);newfused=(newfused/mx).*(~(newfused<0));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5figure(1);hold offimagesc(newfused);axis(′image′,′off′,′xy′);%Xwide=COLS/0.7;Xwide=COLS/0.9;Ywide=round(Xwide*ROWS/COLS);set(gcf,′position′,[55 109 Xwide Ywide]);set(gca,′position′,
);vec=get(gcf,′position′);Ywide=vec(3)*ROWS/COLS;set(gcf,′position′,[10 300 Xwide Ywide]);newfused(,,1)=flipud(newfused(,,1));newfused(,,2)=flipud(newfused(,,2));newfused(,,3)=flipud(newfused(,,3));%imwrite(newfused,′rawmosNOccd2.jpg′,′jpeg′,′quality′,100);function fusedcell=fuseml22(oldcell,newcell)%% oldcell and newcell contain images to put in a single% coordinate system,and the vector-coordinates of the% image centers in (x,y).% fusedcell is the cell of the updated mosaic image% k is the number of images already fused%% merging=′average′ - The overlap area is an average of% all the images that exist in it,and the program% makes all weights be the same.% merging=′avdecay′ - The overlap area is the average of the% old image or mosaic,withthe new image.So,older% images in the old mosaic decay exponentially in time.% merging=′split′- the overlap is cut by half - each side% comes from a different image.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DECODING THE CELLS INTO THEIR CONSTITUENTSZold=oldcell{1,1}; % the 1st (old) imageDELTAZ2old=oldcell{2,1};oldcenterx=oldcell{3,1}; %coordinates of center of old imageoldcentery=oldcell{4,1};Znew=newcell{1,1}; % the 2nd (new) imageDELTAZ2new=newcell{2,1};newcenterx=newcell{3,1};newcentery=newcell{4,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; %the Left side of the new (2nd) imageRnew=newcenterx+(newcols-1)/2;Dnew=newcentery-(newrows-1)/2;Unew=newcentery+(newrows-1)/2;greatL=min([Lold,Lnew]); %the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find ((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find ((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find ((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find ((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;DELTAZoldgreat=NaN*zeros(size(greatX));DELTAZoldgreat(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=DELTAZ2old;%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find ((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find ((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find ((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;DELTAZnewgreat=NaN*zeros(size(greatX));DELTAZnewgreat(newLD_row-lnewLU_row,newLD_colnewRD_col)=DELTAZ2new;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TAKING CARE OF THE OVERLAP AREAnanold=uint8(isnan(Zoldgreat));nannew=uint8(isnan(Znewgreat));fused=NaN*zeros(size(greatX));%fused(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;%fused(oldLD_row-1oldLU_row,oldLD_cololdRD_col)=Zold;besdernew=find(~nannew);fused(besdernew)=Znewgreat(besdernew);besderold=find(-nanold);fused(besderold)=Zoldgreat(besderold);DELTA2fused=NaN*zeros(size(greatX));%DELTA2fused(newLD_row-1newLU_row,newLD_colnewRD_col)=DELTAZ2new;%DELTA2fused(oldLD_row-1oldLU_row,oldLD_coioldRD_col)=DELTAZ2old;DELTA2fused(besdernew)=DELTAZnewgreat(besdernew);DELTA2fused(besderold)=DELTAZoldgreat(besderold);%Sigmaoldgreat=sqrt(DELTAZoldgreat);%Sigmanewgreat=sqrt(DELTAZnewgreat);Sigma2oldgreat=DELTAZoldgreat;Sigma2newgreat=DELTAZnewgreat;nangreat=(nanold | nannew);E=~nangreat;F=find(E);if ~isempty(F) oldneto=NaN*ones(size(Zoldgreat)); newneto=NaN*ones(size(Znewgreat)); oldsigneto=oldneto;   newsigneto=newneto;   oldneto(F)=Zoldgreat(F);   newneto(F)=Znewgreat(F);   oldsig2neto(F)=Sigma2oldgreat(F);   newsig2neto(F)=Sigma2newgreat(F);   %size(oldneto(F))   %size(newneto(F))   %size(oldsig2neto(F)′)   %size(newsig2neto(F)′)   nominat=(oldneto(F)./oldsig2neto(F)′) + (newneto(F)./newsig2neto(F)′);   denomin=(1./oldsig2neto(F)′) + (1./newsig2neto(F)′);   sig2seam=1./denomin;   seamvalue=sig2seam.*nominat;   fused(F)=seamvalue;   DELTA2fused(F)=sigseam.^2;   DELTA2fused(F)=sig2seam;endfused=flipud(fused);DELTA2fused=flipud(DELTA2fused);%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT-COMBINED IMAGEgreatcenterx=(greatR+greatL)/2;greatcentery=(greatU+greatD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERfusedcell{1,1}=fused;fusedcell{2,1}=DELTA2fused;fusedcell{3,1}=greatcenterx;fusedcell{4,1}=greatcentery;function [small,delta2small]=mizoorml21(big,delta2big,levels,shape);newA=big;newdelta2A=delta2big;for iteration=1levels,  [reducedA,reducedelta2A]=reduceml21(newA,newdelta2A,shape);  newA=reducedA;  newdelta2A=reducedelta2A;endsmall=newA;delta2small=newdelta2A;function [small,delta2small]=mizoorml21_xy(big,delta2big,levels_x,levels_y,shape);newA=big;newdelta2A=delta2big;if ~(levels_x==levels_y)   for iteration=1min([levels_x,levels_y]),   [reducedA_x,reducedelta2A_x]=reduceml21_x{newA,newdelta2A,shape);   [reducedA,reducedelta2A]=reduceml21_y(reducedA_x,reducedelta2A_x,shape);  newA=reducedA;   newdelta2A=reducedelta2A;   end   [a,maindirection]=max([levels_x,levels_y]);   if maindirection==1 %%% x-direction is the main movement unceratinty  for iteration=levels_y+llevels_x,  [reducedA,reducedelta2A]=reduceml21_x(newA,newdelta2A,shape);  newA=reducedA;  newdelta2A=reducedelta2A;   end   else %%% y-direction is the main movementunceratinty   for iteration=levels_x+llevels_y,   [reducedA,reducedelta2A]=reduceml21_y(newA,newdelta2A,shape);   newA=reducedA;   newdelta2A=reducedelta2A;   endend else  levels=levels_x;  for iteration=1levels,   [reducedA_x,reducedelta2A_x]=reduceml21_x(newA,newdelta2A,shape);   [reducedA,reducedelta2A]=reduceml21_y(reducedA_x,reducedelta2A_x,shape);   newA=reducedA;   newdelta2A=reducedelta2A;   endendsmall=newA;delta2small=newdelta2A;clearpacklevels_x=3;levels_y=0;minlevel=0;yepsilon=eps;left=97;right=600;bottom=200;roof=385;imagehight=roof-bottom+1;imagewidth=right-left+1;LeftspaceL=135;LeftspaceR=185;RightspaceL=315;RightspaceR=365;Leftspacewidth=LeftspaceR-LeftspaceL+1;Rightspacewidth=RightspaceR-RightspaceL+1;ystep=1;%centerx2guess=100*ystep;%initializing the disparitycenterx2guess=115*ystep;%initializing the disparitycentery2guess=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%basefilename=′stepol_′;postfix=′.tif′;basephotoname=′photo′;rawname=′basicim′;netframes=[115];load flatparams.mat%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%spec=maspect2(basefilename,netframes,left,right,bottom,roof);sig=0.5;invM=flatfun(bottomroof,leftright);SIG2=(sig.^2); sun=250; sunerror=1/eps; SE=ones(3,5); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %h=fspecial(′prewitt′); centerx2=zeros(1,length (netframes)); centery2=zeros(1,length(netframes)); tikunx=zeros(1,length(netframes)); tikuny=zeros(1,length(netframes)); frame_inds=1ystep(length(netframes)-1); k=1; for indframe=1length(frame_inds-1),  %%%%%%%%%%%%%%%%%%%%%%%%%%%% READING AND PREPARING THE RAW IMAGES  yframe_ind=frame_inds(indframe);tic  if k==1   yframe1=netframes(yframe_ind);   sframe1=sprintf(′%d′,yframe1);   shem1=[basefilename sframe1 ′_1′postfix];   shem2=[basefilename sframe1 ′_2′postfix];   shem3=[basefilename sframe1 ′_3′postfix];   rawim1=imread(shem1);   rawim2=imread(shem2);   rawim3=imread(shem3);   rawim=(double(flipud(rawim1))-1) + (double(flipud(rawim2))-1) + (double(flipud(rawim3))-1 );   basicim=(rawim/3);   Za=basicim.*flatfun;   Za=Za(bottomroof,leftright);   X=basicim(bottomroof,leftright);   saturated=(X>sun);   dilsaturated=double(dilate(saturated,SE));   DELTAZa2=(SIG2*(invM.^2)).*(~dilsaturated) + sunerror.*dilsaturated;   DELTAZa2(,LeftspaceL LeftspaceR)=sunerror*ones(imagehight,Leftspacewidth);   DELTAZa2(,RightspaceLRightspaceR)=sunerror*ones(imagehight,Rightspacewidth);   centerx1=0;centery1=0;   holdcell{1,1}=Za;   holdcell{2,1}=DELTAZa2;   holdcell{3,1}=centerx1;   holdcell{4,1}=centeryl;   fusedcell=holdcell;  else   holdcell=fusedcell;  end  yframe2=netframes(yframe_ind+ystep)  sframe1=sprintf(′%d′,yframe2);shem1=[basefi1ename sframe1 ′_1′postfix]; shem2=[basefi1ename sframe1 ′_2′postfix]; shem3=[basefilename sframe1 ′_3′postfix]; rawim1=imread(shem1); rawim2=imread(shem2); rawim3=imread(shem3); rawim=(double(flipud(rawim1))-1) + (double(flipud(rawim2))-1) + (double(flipud(rawim3))-1); basicim=(rawim/3); Zb=basicim.*flatfun; Zb=Zb(bottomroof,leftright); X=basicim(bottomroof,leftright); saturated=(X>sun); dilsaturated=double(dilate(saturated,SE)); DELTAZb2=(SIG2*(invM.^2)).*(-dilsaturated) + sunerror.*dilsaturated; DELTAZb2(,LeftspaceL LeftspaceR)=sunerror*ones(imagehight,Leftspacewidth); DELTAZb2(,RightspaceLRightspaceR)=sunerror*ones(imagehight,Rightspacewidth); hnewcell{1,1}=Zb; hnewcell{2,1}=DELTAZb2; if k==1  hnewcell{3,1}=centerx2guess;  hnewcell{4,1}=centery2guess; else  hnewcell{3,1}=centerx2(old_yframe_ind)+centerx2guess;  hnewcell{4,1}=centery2(old_yframe_ind)+centery2guess; end newcell=hnewcell; clear rawim basicim %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DERIVING TRANSLATION BETWEEN IMAGES[tikunx(yframe_ind+ystep),tikuny(yframe_ind+ystep)]=gmultireml31(holdcell,hnewcell,levels_x,levels_y,minlevel);% centerx2(yframe_ind+ystep)=newcell{3,1} + tikunx(yframe_ind+ystep); centery2(yframe_ind+ystep)=newcell{4,1} + tikuny(yframe_ind+ystep); [tikunx(yframe_ind+ystep) tikuny(yframe_ind+ystep)] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUSING THE IMAGES newcell{3,1}=centerx2(yframe_ind+ystep); newcell{4,1}=centery2(yframe_ind+ystep); %%%%%%%%% NEED TO CORRECT THIS WITH THE NEW KIND OF FUSION fusedcell=fuseml22(fusedcell,newcell); old_yframe_ind=yframe_ind+ystep; k=k+1;tocsave yytemp.mat centerx2 centery2 tikunx tikuny %pause end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %save coordpol2.mat centerx2 centery2 tikunx tikuny fused=fusedcell{1,1}; deltaz2=fusedcell{2,1}; greatcenterx=fusedcell{3,1}; greatcentery=fusedcell{4,1}; %save fusedpoll.mat fused deltaz2 greatcenterx greatcentery %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear fusedcell image(fused) colormap(gray(256)) axis image axis xy set(gcf,′menubar′,′figure′) function prof=profil2(cols,transit) mn=eps^7; tmp=ones(1,cols); L=1transit; R=cols-L+1; %tmp(L)=mn+(1-mn) * (L-1)/(transit-1); %tmp(R)=mn+(1-mn) * (L-1)/(transit-1); tmp(L)=mn+(1-mn) * (1-cos(pi*(L-1)/(transit-1)))/2; tmp(R)=mn+(1-mn) * (1-cos(pi*(L-1)/(transit-1)))/2; prof=tmp; function [oldgreatcell,newgreatcell]=putbigm13(oldcell,newcell) % % oldcell and newcell contain images to put in a single % coordinate system,and the vector-coordinates of the % image centers in (x,y). % oldgreatcell and newgreatcell are the cells of the new large % frame with a common coorinate system,in each of which % the single images are embedded. These cells also contain the % coordinates of the center of the big frames. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DECODING THE CELLS INTO THEIR CONSTITUENTSZold =oldcell{1,1}; % the 1st (old) imageDELTAZ2old=oldcell{2,1};oldcenterx=oldcell{3,1}; %coordinates of center of old imageoldcentery=oldcell{4,1};Znew=newcell{1,1}; % the 2nd (new) imageDELTAZ2new=newcell{2,1};newcenterx=newcell{3,1};newcentery=newcell{4,1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[oldrows,oldcols]=size(Zold);Lold=oldcenterx-(oldcols-1)/2; %the Left side of the old imageRold=oldcenterx+(oldcols-1)/2; % the Right sideDold=oldcentery-(oldrows-1)/2; % the Down of the old imageUold=oldcentery+(oldrows-1)/2; % the Upper side of the old image[newrows,newcols]=size(Znew);Lnew=newcenterx-(newcols-1)/2; % the Left sideof the new (2nd) imageRnew=newcenterx+(newcols-1)/2;Dnew=newcentery-(newrows-1)/2;Unew=newcentery+(newrows-1)/2;greatL=min([Lold,Lnew]); % the Left side of the combined-greater imagegreatR=max([Rold,Rnew]);greatD=min([Dold,Dnew]);greatU=max([Uold,Unew]);spanx=greatR-greatL+1;spany=greatU-greatD+1;greatX=ones(spany,1)*(greatLgreatR); %the coordinate grid of the combined-greater imagegreatY=(greatDgreatU)′*ones(1,spanx);%%%%%%% PUTTING THE OLD IMAGE IN THE COMBINED GREAT-GRID[oldLDy,oldLDx]=find((greatX==Lold) (greatY==Dold));[oldLUy,oldLUx]=find((greatX==Lold) (greatY==Uold));[oldRDy,oldRDx]=find((greatX==Rold) (greatY==Dold));[oldRUy,oldRUx]=find((greatX==Rold) (greatY==Uold));[oldLD_row,oldLD_col]=xy2ij_brack([oldLDx,oldLDy],spany);[oldRD_row,oldRD_col]=xy2ij_brack([oldRDx,oldRDy],spany);[oldLU_row,oldLU_col]=xy2ij_brack([oldLUx,oldLUy],spany);[oldRU_row,oldRU_col]=xy2ij_brack([oldRUx,oldRUy],spany);Zoldgreat=NaN*zeros(size(greatX));Zoldgreat(oldLD_row-loldLU_row,oldLD_cololdRD_col)=Zold;Zoldgreat=flipud(Zoldgreat);DELTAZoldgreat=NaN*zeros(size(greatX));DELTAZoldgreat(oldLD_row-loldLU_row,oldLD_cololdRD_col)=DELTAZ2old;DELTAZoldgreat=flipud(DELTAZoldgreat);%%%%%%%%% PUTTING THE NEW IMAGE IN THE COMBINED GREAT-GRID[newLDy,newLDx]=find ((greatX==Lnew) (greatY==Dnew));[newLUy,newLUx]=find((greatX==Lnew) (greatY==Unew));[newRDy,newRDx]=find((greatX==Rnew) (greatY==Dnew));[newRUy,newRUx]=find((greatX==Rnew) (greatY==Unew));[newLD_row,newLD_col]=xy2ij_brack([newLDx,newLDy],spany);[newRD_row,newRD_col]=xy2ij_brack([newRDx,newRDy],spany);[newLU_row,newLU_col]=xy2ij_brack([newLUx,newLUy],spany);[newRU_row,newRU_col]=xy2ij_brack([newRUx,newRUy],spany);Znewgreat=NaN*zeros(size(greatX));Znewgreat(newLD_row-1newLU_row,newLD_colnewRD_col)=Znew;Znewgreat=flipud(Znewgreat);DELTAZnewgreat=NaN*zeros(size(greatX));DELTAZnewgreat(newLD_row-lnewLU_row,newLD_colnewRD_col)=DELTAZ2new;DELTAZnewgreat=flipud(DELTAZnewgreat);%%%%%%%%%%%%%%%%%% THE CENTRAL COORDINATE OF THE (WARPED) GREAT-COMBINED IMAGEgreatcenterx=(greatR+greatL)/2;greatcentery=(greatU+greatD)/2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PUTTING THE OUTPUT TO CELLS,FOR EASIER TRANSFERoldgreatcell{1,1}=Zoldgreat;oldgreatcell{2,1}=DELTAZoldgreat;oldgreatcell{3,1}=greatcenterx;oldgreatcell{4,1}=greatcentery;newgreatcell{1,1}=Znewgreat;newgreatcell{2,1}=DELTAZnewgreat;newgreatcell{3,1}=greatcenterx;newgreatcell{4,1}=greatcentery;function[small,smalldelta2]=reduceml21(big,bigdelta2,shape);gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25;gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=1./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(conv2(weightedbig,gs,shape),gs′,shape); %convolution in one lineblurdenomin=conv2(conv2(invsig2,gs,shape),gs′,shape);samples_row=1 2 size(blurnominat,1); %sample grid in one linesamples_col=1 2 size(blurnominat,2); %sample grid in one linesmallnominat=blurnominat(samples_row,samples_col); %sampling in one linesmalldenomin=blurdenomin(samples_row,samples_col); %sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig2;function[small,smalldelta2]=reduceml21_x(big,bigdelta2,shape)gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25;gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=1./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(weightedbig,gs,shape);%convolution in one lineblurdenomin=conv2(invsig2,gs,shape);samples_col=12size(blurnominat,2); %sample grid in one linesmallnominat=blurnominat(,samples_col);%sampling in one linesmalldenomin=blurdenomin(,samples_col);%sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig2=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig2;function[small,smalldelta2]=reduceml21_x(big,bigdelta2,shape)gs=zeros(5,1);a=0.4;gs(3)=a;gs(2)=0.25;gs(4)=gs(2);gs(1)=0.25-a/2;gs(5)=gs(1);%invsig=1./sqrt(bigdelta2);invsig2=1./bigdelta2;weightedbig=big.*invsig2;blurnominat=conv2(weightedbig,gs′,shape); %convolution in one lineblurdenomin=conv2(invsig2,gs′,shape);samples_row=1 2 size(blurnominat,1); %sample grid in one linesmallnominat=blurnominat(samples_row,); %sampling in one linesmalldenomin=blurdenomin(samples_row,); %sampling in one linesmall=smallnominat./smalldenomin;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%newsig2=1./smalldenomin;%smalldelta2=newsig.^2;smalldelta2=newsig;function Fethsatur=satursmooth7(X,uppersun,lowersun)SE=ones(3);mn=eps;saturated=uint8(X>uppersun);maybesatur=uint8(X>lowersun);s=max(saturated);if s==0  Fethsatur=ones(size(x));else  shulaim=uint8(1);  oldregion=uint8(saturated);  %k=0;  while~(shulaim==0)   %k=k+l   dilsaturated=(dilate(oldregion,SE));   newregion=(dilsaturated maybesatur);   added=newregion (~(oldregion));   shulaim=max(added);   oldregion=newregion;  end   tmp=ones(size(X));   fact=(1-mn)/(uppersun-(lowersun-1));   for level=lowersunuppersun,   badpixels=((X==level) newregion);   badvalue=mn+(uppersun-level)*fact;   lobeseder=find(badpixels);   if~isempty(lobeseder)  tmp(lobeseder)=badvalue;   end   end   lobeseder=find(saturated);   tmp(lobeseder)=mn;   Fethsatur=tmp;endfunction [enhanced,subenhanced,mnval,mxval]=stretchwin5(X,Y,photo,logphoto)XS=sort(X);YS=sort(Y);L=round(XS(2));R=round(XS(3));D=round(YS(2));U=round(YS(3));%LD LU RU RD LDXp=[L LRR L];Yp=[D UUD D];X0=(L+R)/2;Y0=(D+U)/2;spanx=size(photo,2);spany=size(photo,1);[LD_row,LD_col]=xy2ij_brack([L,D],spany);[RD_row,RD_col]=xy2ij_brack([R,D],spany);[LU_row,LU_col]=xy2ij_brack([L,U],spany);[RU_row,RU_col]=xy2ij_brack([R,U],spany);stam=flipud(photo);subarea=stam(LD_row-1LU_row,LD_colRD_col);mn=min(subarea);lemata=subarea-mn;mx=max(lemata);imx8=(2^8)/mx;subenhanced=uint8(round(lemata*imx8));mn=min(logphoto);lemata=double(logphoto)-double(mn);mx=max(lemata);imx8=(2^8)/mx;stam=flipud(lemata*imx8);stam(LD_row-1LU_row,LD_colRD_col)=subenhanced;enhanced=flipud(stam);mnval=min(subarea);mxval=max(subarea);% [I,J]=xy2ij_brack(XYcoord,N);% Converts x-y coordinates to ij coordinates% I is the row-number of the element,% J is the column-number of the element% N is the number of rows in the matrix.function [I,J]=xy2ij_brack(XYcoord,N);%XYcoordX=XYcoord(,1);Y=XYcoord(,2);%index in row-numberI=N-Y+1;%index in colsJ=X;%IJcoord=[I,J];]]>
權利要求
1.一種成像方法,其特徵在於包括用成像器作第一組測量以生成第一像值的第一步驟,第一組測量包括至少測量一次第一景區第一射線束的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束具有第一強度靈敏度特性,成像器對第一主射線射線束的強度具有第一動態範圍;用成像器作第二組測量以生成第二像值的第二步驟,第二組測量包括至少測量一次第一景區發射的第二射線束的強度,第二射線束在成像器參考幀內具有第二主射線,第二主射線不同於第一主射線,成像器對第二主射線射線束具有第二強度靈敏度特性,第二與第一強度靈敏度特性不同,而且成像器對第二主射線射線束的強度具有第二動態範圍;和對第一與第二像值作拼合操作而生成第三像值,第三像值與成像器對第一和第二射線束強度中至少一個的第三動態範圍相關,第三動態範圍大於成像器第一與第二動態範圍中的至少一個。
2.根據權利要求1的方法,其特徵在於第一射線束包括第一電磁射線束,第二射線束包括第二電磁射線束。
3.根據權利要求1的方法,其特徵在於還包括下列方法之一在第一與第二步驟之間相對第一景區轉動成像器;和在第一與第二步驟之間相對第一景區平移成像器。
4.根據權利要求1的方法,還包括校正成像器,其特徵在於,校正步驟包括用成像器測量第一組第一主射線射線束的強度以生成第一組校正測量值;通過測定下列值之一確定第一強度靈敏度特性的第一估值a)第一組校正測量值之和,和b)第一組校正測量值的平均值;用成像器測量第二組第二主射線射線束的強度,生成第二組校正測量值;和通過測定下列值之一確定第二強度靈敏度特性的第二估值a)第二組校正測量值之和,和b)第二組校正測量值的平均值。
5.根據權利要求1的方法,還包括校正成像器,其特徵在於,校正步驟包括用成像器作第三組測量而生成第四像值,第三組測量包括至少測量一次第三射線束的強度,第三射線束由第二景區發射並居於第一主射線;用成像器作第四組測量而生成第五像值,第四組測量包括至少測量一次第四射線束的強度,第四射線束由第二景區發射並且有第二主射線;和評估第一與第二強度靈敏度特性間的關係,評估步驟包括測定下列值之一a)第四與第五像值之差,和b)第四與第五像值之比。
6.根據權利要求1的方法,其特徵在於還包括用成像器作第三組測量而生成第四像值,第三組測量包括至少測量一次第三射線束至少一種選擇譜分量的強度,第三射線束由第一景區發射且在成像器參考幀內具有第三主射線,成像器對第三主射線射線束具有第一光譜靈敏度特性,第一光譜靈敏度特性包括具有第一波長靈敏度帶的帶通特性,第三射線束至少一個選擇的譜分量在第一波長靈敏度帶內有一波長;和用成像器作第四組測量而生成第五像值,第四組測量包括至少測量一次第四射線束至少一個選擇譜分量的強度,第四射線束由第一景區發射且在成像器參考幀內具有第四主射線,第四主射線不同於第三主射線,成像器對第四主射線射線束具有第二光譜靈敏度特性,第二光譜靈敏度特性包括具有第二波長靈敏度帶的帶通特性,第四射線束至少一個選擇的譜分量在第二波長靈敏度帶內有一波長,而第二與第一波長靈敏度帶不同。
7.根據權利要求1的方法,其特徵在於還包括用成像器作第三組測量而生成第四像值,第三組測量包括至少測量一次第一景區發射的第三射線束中至少一個選擇的偏振分量的強度,第三射線束在成像器參考幀內具有第三主射線,成像器對第三主射線射線束具有第一偏振靈敏度特性,第一偏振靈敏度特性包括對偏振角在第一角範圍外的射線分量降低的靈敏度,而第三射線束至少一個選擇的偏振分量在第一角範圍內有一偏振角;和用成像器作第四組測量而生成第五像值,第四組測量包括至少測量一次第四射線束至少一個所選偏振分量的強度,第四射線束由第一景區發射且在成像器參考幀內具有第四主射線,第四與第三主射線不同,成像器對第四主射線射線束具有第二偏振靈敏度特性,第二偏振靈敏度特性包括對偏振角在第二角範圍外的信號分量降低的靈敏度,第四射線束至少一個所選偏振分量在第二角範圍內有一偏振角,而第二與第一角範圍不同。
8.根據權利要求1的方法,其特徵在於還包括用成像器作第三組測量而生成第四像值,第三組測量包括至少測量一次第一景區發射的第三射線束的強度,第三射線束在成像器參考幀內具有第三主射線,成像器對第三主射線射線束具有第一聚焦特性,第一聚焦特性包括第一焦距;和用成像器作第四組測量而生成第五像值,第四組測量包括至少測量一次第一景區發射的第四射線束的強度,第四射線束在成像器參考幀內具有第四主射線,第四與第三主射線不同,成像器對第四主射線射線束具有第二聚焦特性,第二聚焦特性包括第二焦距,第二與第一焦距不同。
9.一種成像方法,其特徵在於包括用成像器作第一組測量而生成第一像值的第一步驟,第一組測量包括至少測量一次第一景區發射的第一射線束中至少一個所選偏振分量的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束具有第一偏振靈敏度特性,而第一偏振靈敏度特性包括對偏振角在第一角範圍外的信號分量降低的靈敏度,第一射線束至少一個所選偏振分量在第一角範圍內有一偏振角;用成像器作第二組測量而生成第二像值的第二步驟,第二組測量包括至少測量一次第一景區發射的第二射線束中至少一個所選偏振分量的強度,第二射線束在成像器參考幀內具有第二主射線,第二與第一主射線不同,成像器對第二主射線射線束具有第二偏振靈敏度特性,第二偏振靈敏度特性包括對偏振角在第二角範圍外的信號分量降低的靈敏度,第二射線束至少一個所選偏振分量在第二角範圍內有一偏振角,而第二與第一角範圍不同;移動成像器的第三步驟,包括下列之一在第一與第二步驟之間相對第一景區轉動成像器;和在第一與第二步驟之間相對第一景區平移成像器;和用第一與第二像值確定第一與第二射線束之一的偏振態。
10.根據權利要求9的方法,其特徵在於,第一射線束包括第一電磁射線束,第二射線束包括第二電磁射線束。
11.一種成像設備,其特徵在於包括用於作第一組測量而生成第一像值的第一成像器,第一組測量包括至少測量一次第一景區發射的第一射線束的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束具有第一強度靈敏度特性,而對第一主射線射線束的強度具有第一動態範圍;用於作第二組測量而生成第二像值的第二成像器,第二組測量包括至少測量一次第一景區發射的第二射線束的強度,第二射線束在成像器參考幀內具有第二主射線,第二與第一主射線不同,成像器對第二主射線射線束具有第二強度靈敏度特性,第二與第一強度靈敏度特性不同,成像器對第二主射線輻射信號組的強度具有第二動態範圍;和處理器,用於對第一和第二測量值作拼合操作而生成第三像值,第三像值與相對於第一和第二射線束強度中至少一個的第三動態範圍相關,第三動態範圍大於第一和第二動態範圍中的至少一個。
12.根據權利要求11的設備,其特徵在於第一射線束包括第一電磁射線束,第二射線束包括第二電磁射線束。
13.根據權利要求11的設備,其特徵在於第二成像器就是第一成像器,第一組測量不晚於第一時間,第二組測量不早於第二時間,第二時間晚於第一時間,而且該設備還包括下列結構之一在第一與第二時間之間相對於第一和第二景區轉動第一成像器的結構;和在第一與第二時間之間相對於第一和第二景區平移第一成像器的結構。
14.根據權利要求11的設備,其特徵在於還包括用於用成像器測量第一組第一主射線射線束的強度的處理器,生成第一組校正測量值;用於測定第一強度靈敏度特性的第一估值的處理器,包括下列之一a)測定第一組校正測量值之和的處理器,和b)測定第一組校正測量值平均值的處理器;用於用成像器測量第二組第二主射線射線束的強度的處理器,生成第二組校正測量值;和測定第二強度靈敏度特性的第二估值的處理器,包括下列之一a)測定第二組校正測量值之和的處理器,和b)測定第二組校正測量值平均值的處理器。
15.根據權利要求11的設備,其特徵在於還包括用於作第三組測量而生成第四像值的第三成像器,第三組測量包括至少測量一次第二景區發射的第三射線束的強度,第三射線束具有第一主射線;用於作第四組測量而生成第五像值的第四成像器,第四組測量包括至少測量一次第二景區發射的第四射線束的強度,第四射線束具有第二主射線;和用於評估第一與第二強度靈敏度特性的關係的處理器,包括下列之一a)測定第四與第五像值之差的處理器,和b)測定第四與第五像值之比的處理器。
16.根據權利要求11的設備,其特徵在於還包括用於作第三組測量而生成第四像值的第三成像器,第三組測量包括至少測量一次第一景區發射的第三射線束中至少一個所選譜分量的強度,第三射線束在成像器參考幀內具有第三主射線,成像器對第三主射線射線束具有第一光譜靈敏度特性,第一光譜靈敏度特性包括具有第一波長靈敏度帶的帶通特性第三射線束至少一個所選譜分量有一位於第一波長靈敏度帶內的波長;用於作第四組測量而生成第五像值的第四成像器,第四組測量包括至少測量一次第一景區發射的第四射線束中至少一個所選譜分量的強度,第四射線束在成像器參考幀內具有第四主射線,第四主射線與第三主射線不同,成像器對第四主射線射線束具有第二光譜靈敏度特性,第二光譜靈敏度特性包括具有第二波長靈敏度帶的帶通特性,第四射線束至少一個所選譜分琅具有位於第二波長靈敏度帶內的波長,第二與第一波長靈敏度帶不同。
17.根據權利要求11的設備,其特徵在於還包括用於作第三組測量而生成第四像值的第三成像器,第三組測量包括至少測量一次第一景區發射的第三射線束中至少一個所選偏振分量的強度,第三射線束在成像器參考幀內具有第三主射線,成像器對第三主射線射線束具有第一偏振靈敏度特性,第一偏振靈敏度特性包括對偏振角在第一角範圍外的信號分量降低的靈敏度,而第三射線束至少一個所選偏振分量的偏振角在第一角範圍內;用於作第四組測量而生成第五像值的第四成像器,第四組測量包括至少測量一次第一景區發射的第四射線束中至少一個所選偏振分量的強度,第四射線束在成像器參考幀內具有第四主射線,第四與第三主射線不同,成像器對第四主射線射線束具有第二偏振靈敏度特性,第二偏振靈敏度特性包括對偏振角在第二角範圍外的信號分量降低的靈敏度,第四射線束至少一個所選偏振分量的偏振角在第二角範圍內,第二與第一角範圍不同。
18.根據權利要求11的設備,其特徵在於還包括用於作第三組測量而生成第四像值的第三成像器,第三組測量包括至少測量一次第一景區發射的第三射線束的強度,第三射線束在成像器參考幀內具有第三主射線,而且成像器對第三主射線射線束具有第一聚焦特性,第一聚焦特性包括第一焦距;用於作第四組測量而生成第五像值的第四成像器,第四組測量包括至少測量一次第一景區發射的第四射線束的強度,第四射線束在成像器參考幀內具有第四主射線,第四與第三主射線不同,成像器對第四主射線射線束具有第二聚焦特性,第二聚焦特性包括第二焦距,第二與第一焦距不同。
19.一種成像設備,其特徵在於包括用於作第一組測量而生成第一像值的第一成像器,第一組測量包括至少測量一次第一景區發射的第一射線束中至少一個所選偏振分量的強度,第一射線束在成像器參考幀內具有第一主射線,成像器對第一主射線射線束參考幀內具有第一偏振靈敏度特性,第一偏振靈敏度特性包括對偏振角在第一角範圍外的信號分量減低的靈敏度,第一射線束至少一個所選偏振分量的偏振角在第一角範圍內;用於作第二組測量而生成第二像值的第二成像器,第二組測量包括至少測量一次第一景區發射的第二射線束中至少一個所選偏振分量的強度,第二射線束在成像器參考幀內具有第二主射線,第二與第一主射線不同,成像器對第二主射線射線束具有第二偏振靈敏度特性,第二偏振靈敏度特性包括對偏振角在第二角範圍外的信號分量減低的靈敏度,第二射線束至少一個所選偏振分量的偏振角在第二角範圍內,而第二與第一角範圍不同,其中第一組測量不晚於第一時間,第二組測量不早於第二時間,第二時間晚於第一時間;一種移動裝置,包括下列之一在第一與第二時間之間相對於第一和第二景區轉動第一成像器的結構,和在第一與第二時間之間相對於第一和第二景區平移第一成像器的結構;和應用第一和第二像值確定第一與第二射線束之一的偏振態的處理器。
20.根據權利要求19的設備,其特徵在於,第一射線束包括第一電磁射線束,第二射線束包括第二電磁射線束。
全文摘要
提供一種用攝像機或其他成像靈敏度特性在其視角內變化的成像器捕獲圖象的方法和設備。成像器的特性相對曝光、色靈敏度、偏振靈敏度、焦距和/或圖象檢測任一其他方面可以不均一。成像器可轉動或平移,以便捕獲被成像的不同景部。由於成像器在捕獲各景部時處於多個位置,各景部可用成像器靈敏度分布曲線的多個部分成像。
文檔編號G06T7/00GK1451230SQ01813115
公開日2003年10月22日 申請日期2001年7月23日 優先權日2000年7月21日
發明者Y·Y·謝希納, S·K·納亞爾 申請人:紐約市哥倫比亞大學託管會

同类文章

一種新型多功能組合攝影箱的製作方法

一種新型多功能組合攝影箱的製作方法【專利摘要】本實用新型公開了一種新型多功能組合攝影箱,包括敞開式箱體和前攝影蓋,在箱體頂部設有移動式光源盒,在箱體底部設有LED脫影板,LED脫影板放置在底板上;移動式光源盒包括上蓋,上蓋內設有光源,上蓋部設有磨沙透光片,磨沙透光片將光源封閉在上蓋內;所述LED脫影

壓縮模式圖樣重疊檢測方法與裝置與流程

本發明涉及通信領域,特別涉及一種壓縮模式圖樣重疊檢測方法與裝置。背景技術:在寬帶碼分多址(WCDMA,WidebandCodeDivisionMultipleAccess)系統頻分復用(FDD,FrequencyDivisionDuplex)模式下,為了進行異頻硬切換、FDD到時分復用(TDD,Ti

個性化檯曆的製作方法

專利名稱::個性化檯曆的製作方法技術領域::本實用新型涉及一種檯曆,尤其涉及一種既顯示月曆、又能插入照片的個性化檯曆,屬於生活文化藝術用品領域。背景技術::公知的立式檯曆每頁皆由月曆和畫面兩部分構成,這兩部分都是事先印刷好,固定而不能更換的。畫面或為風景,或為模特、明星。功能單一局限性較大。特別是畫

一種實現縮放的視頻解碼方法

專利名稱:一種實現縮放的視頻解碼方法技術領域:本發明涉及視頻信號處理領域,特別是一種實現縮放的視頻解碼方法。背景技術: Mpeg標準是由運動圖像專家組(Moving Picture Expert Group,MPEG)開發的用於視頻和音頻壓縮的一系列演進的標準。按照Mpeg標準,視頻圖像壓縮編碼後包

基於加熱模壓的纖維增強PBT複合材料成型工藝的製作方法

本發明涉及一種基於加熱模壓的纖維增強pbt複合材料成型工藝。背景技術:熱塑性複合材料與傳統熱固性複合材料相比其具有較好的韌性和抗衝擊性能,此外其還具有可回收利用等優點。熱塑性塑料在液態時流動能力差,使得其與纖維結合浸潤困難。環狀對苯二甲酸丁二醇酯(cbt)是一種環狀預聚物,該材料力學性能差不適合做纖

一種pe滾塑儲槽的製作方法

專利名稱:一種pe滾塑儲槽的製作方法技術領域:一種PE滾塑儲槽一、 技術領域 本實用新型涉及一種PE滾塑儲槽,主要用於化工、染料、醫藥、農藥、冶金、稀土、機械、電子、電力、環保、紡織、釀造、釀造、食品、給水、排水等行業儲存液體使用。二、 背景技術 目前,化工液體耐腐蝕貯運設備,普遍使用傳統的玻璃鋼容

釘的製作方法

專利名稱:釘的製作方法技術領域:本實用新型涉及一種釘,尤其涉及一種可提供方便拔除的鐵(鋼)釘。背景技術:考慮到廢木材回收後再加工利用作業的方便性與安全性,根據環保規定,廢木材的回收是必須將釘於廢木材上的鐵(鋼)釘拔除。如圖1、圖2所示,目前用以釘入木材的鐵(鋼)釘10主要是在一釘體11的一端形成一尖

直流氧噴裝置的製作方法

專利名稱:直流氧噴裝置的製作方法技術領域:本實用新型涉及ー種醫療器械,具體地說是ー種直流氧噴裝置。背景技術:臨床上的放療過程極易造成患者的局部皮膚損傷和炎症,被稱為「放射性皮炎」。目前對於放射性皮炎的主要治療措施是塗抹藥膏,而放射性皮炎患者多伴有局部疼痛,對於止痛,多是通過ロ服或靜脈注射進行止痛治療

新型熱網閥門操作手輪的製作方法

專利名稱:新型熱網閥門操作手輪的製作方法技術領域:新型熱網閥門操作手輪技術領域:本實用新型涉及一種新型熱網閥門操作手輪,屬於機械領域。背景技術::閥門作為流體控制裝置應用廣泛,手輪傳動的閥門使用比例佔90%以上。國家標準中提及手輪所起作用為傳動功能,不作為閥門的運輸、起吊裝置,不承受軸向力。現有閥門

用來自動讀取管狀容器所載識別碼的裝置的製作方法

專利名稱:用來自動讀取管狀容器所載識別碼的裝置的製作方法背景技術:1-本發明所屬領域本發明涉及一種用來自動讀取管狀容器所載識別碼的裝置,其中的管狀容器被放在循環於配送鏈上的文檔匣或託架裝置中。本發明特別適用於,然而並非僅僅專用於,對引入自動分析系統的血液樣本試管之類的自動識別。本發明還涉及專為實現讀