面向爆炸壓力場的體繪制自動傳輸函數(shù)設(shè)計

陸曼君 解利軍 祁佳晨 王攀 姚成寶

收稿日期：2023-05-30；修回日期：2023-07-17

作者簡介：陸曼君（1997—），女，浙江杭州人，碩士研究生，主要研究方向為科學可視化；解利軍（1980—），男（通信作者），山西朔州人，副教授，碩導(dǎo)，博士，主要研究方向為科學可視化（Lijun_Xie@zju.edu.cn）；祁佳晨（1998—），男，浙江杭州人，碩士研究生，主要研究方向為計算機視覺；王攀（1983—），男，內(nèi)蒙古包頭人，助理研究員，博士，主要研究方向為高性能計算、虛擬現(xiàn)實與可視化；姚成寶（1984—），男，安徽廬江人，副研究員，博士，主要研究方向為爆炸與沖擊動力學數(shù)值模擬．

摘? 要：直接體繪制廣泛應(yīng)用于可視化領(lǐng)域，其中傳輸函數(shù)決定了最終繪制的效果。大規(guī)模爆炸數(shù)據(jù)屬于復(fù)雜的時變數(shù)據(jù)，其壓力場的取值范圍和分布會隨時間劇烈變化。為了實現(xiàn)對沖擊波傳播特征的連續(xù)跟蹤，提出了一種基于標量-時間梯度（scalar-temporal gradient，STG）二維空間的自動傳輸函數(shù)設(shè)計方法。結(jié)合爆炸過程中的物理規(guī)律，將STG空間劃分為多種狀態(tài)并總結(jié)了狀態(tài)之間的轉(zhuǎn)移規(guī)律。將沖擊波的傳播特征編碼為一段一定長度的狀態(tài)轉(zhuǎn)移鏈（state transitions chain，STC），稱為時變特征。傳輸函數(shù)中的顏色以時變特征為單位，基于完全狀態(tài)轉(zhuǎn)移鏈自動分配。傳輸函數(shù)中的不透明度采用基于壓強變化速率的分段非線性曲線，凸顯壓強變化劇烈的特征。選取了兩類典型的大規(guī)模爆炸場景對方法進行了驗證，繪制結(jié)果中都包含了豐富的明顯可區(qū)分的沖擊波傳播特征。結(jié)果表明，生成的傳輸函數(shù)在整個時間維度上都能穩(wěn)定、準確地捕捉?jīng)_擊波的傳播特征。

關(guān)鍵詞：直接體繪制；傳輸函數(shù)；時變數(shù)據(jù)；沖擊波；狀態(tài)轉(zhuǎn)移鏈

中圖分類號：TP391.41??? 文獻標志碼：A

文章編號：1001-3695（2024）02-043-0602-07

doi：10.19734/j.issn.1001-3695.2023.05.0284

Volume rendering automatic transfer function design for explosion pressure field

Lu Manjun1，Xie Lijun1，Qi Jiachen1，Wang Pan2，Yao Chengbao2

（1.School of Aeronautics & Astronautics，Zhejiang University，Hangzhou 310027，China；2.Northwest Institute of Nuclear Technology，Xian 710024，China）

Abstract：Direct volume rendering is widely used in scientific visualization.The transfer function largely determines the effect of rendering results.Explosion data is large-scale complex time-varying volume data.The value range and distribution of its pressure field will change drastically with time.In order to realize continuous tracking of shock wave propagation characteristics，this paper proposed an automatic transfer function design method based on scalar-temporal gradient （STG） 2D space.Combined with the physical laws in the explosion process，the method divided the STG space into multiple states and summarized the transition laws between states.It encoded the propagation characteristics of shock waves with state transition chain（STC） of a certain length，called time-varying features.The color assignment in the transfer function took time-varying features as the basic unit，which was done automatically through a complete state transition chain.The opacity in the transfer function used a piecewise non-linear curve based on the pressure change rate，highlighting the characteristics of severe pressure changes.The paper selected two typical large-scale explosion scenarios to verify the method.The rendering results all contain a wealth of clearly distinguishable shock wave propagation characteristics.The results show that the generated transfer function can stably and accurately capture the propagation characteristics of shock waves in the whole-time dimension.

Key words：direct volume rendering；transfer function；time-varying data；shock wave；state transition chain

0? 引言

從實際生活中采集或者計算模擬得到的數(shù)據(jù)包含大量時變數(shù)據(jù)，如氣象采集數(shù)據(jù)、流體數(shù)據(jù)等。面向時變數(shù)據(jù)的可視化研究是科學可視化的重要分支。爆炸壓力場數(shù)據(jù)屬于復(fù)雜的時變數(shù)據(jù)。一般的時變數(shù)據(jù)具有兩個特點：屬性取值范圍和分布會隨時間改變，以及數(shù)據(jù)規(guī)模大（包含時間維度）。相較于一般的時變數(shù)據(jù)，爆炸壓力場中壓強的取值和分布隨時間的變化更加劇烈，其壓強的取值在極短的時間內(nèi)能夠跨越幾個甚至十幾個數(shù)量級。這類爆炸數(shù)據(jù)的可視化難點在于無法保證特征定義在時間維度上的一致性和實現(xiàn)對特征的連續(xù)跟蹤。例如應(yīng)用直接體繪制展示爆炸特征的演變規(guī)律時，無法將劇烈變化的數(shù)據(jù)穩(wěn)定地映射到有限的顏色空間中。為了保證前后時間幀映射關(guān)系的一致性，就必然導(dǎo)致大量關(guān)鍵信息的丟失。如果為了盡可能保留關(guān)鍵信息而不考慮前后時間幀映射關(guān)系的一致性，就會產(chǎn)生劇烈的跳變。因此，穩(wěn)定追蹤爆炸壓力場的特征并進行可視化展示是一個重要的研究問題。

特征追蹤的工作可以分為空間-時間并重方法、時間優(yōu)先方法和空間優(yōu)先方法三類［1］。空間-時間并重方法直接將時間作為新的維度加入，將時間序列數(shù)據(jù)看做一個四維數(shù)據(jù)。Kosara等人［2］將時間作為新的維度添加到傳統(tǒng)的直方圖中，提出了時間直方圖。Akiba等人［3］將時間直方圖引入了傳輸函數(shù)的設(shè)計過程中，實現(xiàn)了對整個時間序列的同時分類。這類可視化方法可以將面向靜態(tài)數(shù)據(jù)的可視化方法進一步推廣到時變數(shù)據(jù)的可視化工作中，統(tǒng)一提取特征。但是爆炸時變數(shù)據(jù)中包含成百上千個時間幀，同時分析和統(tǒng)一提取特征將面臨巨大的計算量，且結(jié)果容易受噪聲干擾產(chǎn)生波動［1］。

時間優(yōu)先方法將時變數(shù)據(jù)轉(zhuǎn)換為一個三維數(shù)組，每個元素都包含了一條時間相關(guān)的曲線。Fang等人［4］將時間活動曲線（time-activity curve，TAC）作為一類時間相關(guān)的屬性，計算不同TAC之間的相似度作為分類依據(jù)，定義并評估了三種不同的相似性度量方法。Lee等人［5］將TAC定義為時變特征，應(yīng)用動態(tài)時間規(guī)整（dynamic time warping，DTW）計算不同TAC之間的相似度。雖然此類方法適用于所有類型的時變數(shù)據(jù)，但是聚類的結(jié)果和數(shù)據(jù)本身包含的物理規(guī)律之間無法建立明確的對應(yīng)關(guān)系，且其聚類結(jié)果不穩(wěn)定，同一類物理特征在不同數(shù)據(jù)集中的聚類結(jié)果無法保證一致。

空間優(yōu)先方法首先從每一幀數(shù)據(jù)中提取空間特征，然后基于預(yù)設(shè)的匹配規(guī)則匹配相鄰時間幀的特征［6～9］。Silver等人［7］從每個時間幀中提取特征，基于空間重疊匹配特征，首先提出了通過有向無環(huán)圖（DAG）描述特征的演變情況。相對于上述兩種方法，空間優(yōu)先方法對特征的定義是基于空間的，具有物理直觀性，但其對特征的跟蹤建立在空間上的重疊，缺少對時變特征的描述能力。

爆炸壓力場數(shù)據(jù)中最重要的特征是沖擊波的傳播，其相關(guān)研究非常豐富［10～13］。在復(fù)雜環(huán)境中，沖擊波的傳播情況非常復(fù)雜，會產(chǎn)生反射、繞射、干涉等現(xiàn)象。雖然沖擊波的具體傳播過程會因為場景等因素存在差異，但是其中包含了一致的物理規(guī)律和共同特征。這些共同特征是爆炸可視化分析中的重點關(guān)注對象。雖然前述三類特征的定義和追蹤方法都可以應(yīng)用到爆炸壓力場可視化中，但是因為其對特征的定義是通用的，難以捕捉爆炸壓力場中特定的時變特征。所以，本文基于對爆炸壓力場的分析，提出了特定的捕捉和描述爆炸壓力場時變特征的方法。

特征演變規(guī)律可視化的方式主要分為抽象表示［6，14，15］和基于繪制的表示［16～18］兩種。前者主要通過圖、樹等數(shù)據(jù)結(jié)構(gòu)表示特征的演變規(guī)律；后者主要通過光學屬性的分配，以關(guān)鍵幀枚舉［16］、動畫［19］的形式展示特征的演變規(guī)律。基于繪制展示特征演變規(guī)律的方法更直觀，更好地利用了人類視覺感知顏色和捕捉變化的能力，在實際應(yīng)用中更為常見。其中，直接體繪制（DVR）方法直接對三維數(shù)據(jù)場進行數(shù)據(jù)分類和顏色合成，能實現(xiàn)對整個體的穿透式顯示，在爆炸壓力場可視化中廣受歡迎［20］。

直接體繪制的核心是傳輸函數(shù)的設(shè)計。傳輸函數(shù)決定了標量屬性值與光學屬性之間的映射關(guān)系。選取合適的傳輸函數(shù)能夠幫助用戶快速捕捉海量數(shù)據(jù)中的關(guān)鍵信息。傳輸函數(shù)的設(shè)計工作已經(jīng)非常豐富［21～23］，但是面向時變特征的傳輸函數(shù)設(shè)計工作依然存在很多挑戰(zhàn)，例如無法將劇烈變化的爆炸數(shù)據(jù)穩(wěn)定地映射到有限的顏色空間中。就筆者調(diào)研所知，目前暫時沒有針對短時間內(nèi)急劇變化的爆炸數(shù)據(jù)自動體繪制相關(guān)的研究工作。結(jié)合上述特征提取方法，本文提出了一種體繪制傳輸函數(shù)的自動設(shè)計方法。

本文主要的貢獻如下：

a）基于標量-時間梯度（scalar-temporal gradient，STG）二維空間分析和總結(jié)了沖擊波的傳播規(guī)律，將其定義為一段狀態(tài)轉(zhuǎn)移鏈（state transitions chain，STC），稱為時變特征。

b）提出一種面向爆炸壓力場數(shù)據(jù)的自動傳輸函數(shù)方法。顏色映射以時變特征為基本單位自動進行。不透明度映射曲線為基于壓強變化速率的分段非線性曲線。

1? 基于標量-時間梯度二維空間的傳輸函數(shù)設(shè)計

應(yīng)用本文提出的傳輸函數(shù)獲取繪制結(jié)果的流程如圖1所示。圖中展示了從數(shù)據(jù)輸入到獲取繪制結(jié)果的整個流程。其中存在兩條支線，分別對應(yīng)顏色和不透明度的映射過程。傳輸函數(shù)設(shè)計的流程如圖2所示。實線箭頭指示了設(shè)計過程中主要的步驟，虛線箭頭指示了兩個主要步驟之間具體的設(shè)計步驟。

1.1? 狀態(tài)定義

原始爆炸壓力場數(shù)據(jù)只包含壓強屬性，本節(jié)首先在壓強屬性的基礎(chǔ)上計算了時間梯度。時間梯度的符號和絕對值分別表示壓強的變化趨勢和變化速率，可以描述壓強在相鄰時間幀的變化情況。應(yīng)用前向差分計算得到除第一幀外所有時間幀的時間梯度屬性：

GiT=Pi-Pi-1Δt（1）

其中：Δt表示相鄰時間幀之間的時間間隔；GiT、Pi分別表示時間幀i對應(yīng)的時間梯度GT和壓強P。在模擬數(shù)據(jù)中，相鄰時間幀的時間間隔通常設(shè)置為固定值。為了加速計算，本文在實際計算中使用幀序差1替代。

構(gòu)建STG二維空間，以標量屬性為橫坐標、時間梯度為縱坐標。對于除第一幀外的每一幀數(shù)據(jù)，將所有數(shù)據(jù)點按照其對應(yīng)的壓強和時間梯度映射到STG二維空間中，得到一張STG散點圖。

STG散點圖可以用于分析屬性的分布規(guī)律和變化速率，所以本節(jié)基于STG散點圖分析和編碼爆炸壓力場的基本變化特征。以城市爆炸數(shù)據(jù)為例（屬于近地爆炸，可視化結(jié)果見2.1.3節(jié)），圖3展示了整個爆炸過程中沖擊波傳播處于不同階段的三張散點圖。圖（a）對應(yīng)最外層沖擊波與建筑群相遇前的階段，散點分布規(guī)則且密集；圖（b）對應(yīng)最外層沖擊波在建筑群中傳播的階段，散點分布散亂；圖（c）對應(yīng)最外層沖擊波離開建筑群后的階段，散點的分布情況介于圖（a）（b）。三張散點圖中均能清晰地觀察到由散點聚集形成的左端點為（u，0）、斜率為1的射線（壓強初始值為u，城市爆炸數(shù)據(jù)中u=101301）。所有的散點均分布在坐標空間的右下區(qū)域（圖4中藍色區(qū)域，下文稱其為目標區(qū)域），區(qū)域的上界為GT=P，左邊界為P=0。

通過分析各類爆炸數(shù)據(jù)后，可以將目標區(qū)域劃分為如圖4所示的八個子區(qū)域，分別對應(yīng)八種不同的狀態(tài)，具體的劃分依據(jù)如表1所示。除S、A外的所有狀態(tài)對應(yīng)的數(shù)值范圍是否包含邊界值對結(jié)果的影響可以忽略。狀態(tài)S對應(yīng)的區(qū)域只包含一個坐標點（u，0），是所有數(shù)據(jù)點的初始狀態(tài)。狀態(tài)A是以（u，0）為左端點的射線。射線對應(yīng)的數(shù)據(jù)點在當前時間幀初次受到?jīng)_擊波的影響，位于最外層沖擊波所在區(qū)域。因為數(shù)據(jù)點受到?jīng)_擊波影響后壓強的增長幅度不相同，就形成了一條斜線。根據(jù)式（1）可知，該射線的斜率為1/Δt。只有位于斜率為1/Δt的直線上的數(shù)據(jù)點，其前一幀的壓強取值相等且取值與這條特殊直線的截距b在數(shù)值上成比例關(guān)系。

b=GT-PΔt（2）

本文在實際計算中Δt取1，所以射線斜率為1。因為壓力場非負的特性，所以任何數(shù)據(jù)點只能映射到目標區(qū)域內(nèi)。

1.2? 時變特征定義

數(shù)據(jù)點所處狀態(tài)會隨時間發(fā)生轉(zhuǎn)換，但是這種轉(zhuǎn)換并不是任意的，必須符合爆炸的物理規(guī)律。當數(shù)據(jù)點受到?jīng)_擊波的影響，壓強增長達到頂峰后，通常會開始下降，即狀態(tài)A、E通常轉(zhuǎn)換到狀態(tài)F。數(shù)據(jù)點在受到最外層沖擊波的影響后，還會受到繞射波或反射波的影響，轉(zhuǎn)換到狀態(tài)C、E或D。狀態(tài)之間所有可能的轉(zhuǎn)換關(guān)系使用圖的形式連接，形成了一張有向圖，如圖5所示。數(shù)據(jù)模擬中壓力場初始值相同，所有數(shù)據(jù)點的初始狀態(tài)必定為S。數(shù)據(jù)點最初受到?jīng)_擊波的影響，壓強呈增長趨勢，狀態(tài)S只能轉(zhuǎn)移到狀態(tài)A。圖5實線框中描述了一組（3個）狀態(tài)之間的轉(zhuǎn)移規(guī)律，綠色箭頭指示了當前狀態(tài)可能轉(zhuǎn)移到的一組狀態(tài)。其中組內(nèi)狀態(tài)要轉(zhuǎn)換到組外狀態(tài)一定會經(jīng)過灰色的狀態(tài)（參見電子版）。一個數(shù)據(jù)點在整個生命周期中的所有狀態(tài)可以按時間順序連接起來形成一條單鏈，除去未受到?jīng)_擊波影響的轉(zhuǎn)換關(guān)系S→S后，得到一條狀態(tài)轉(zhuǎn)移鏈。

總結(jié)了不同爆炸場景和不同空間位置的數(shù)據(jù)點對應(yīng)的STC的特點后，本文將沖擊波的傳播特征編碼為一段長度為2的STC，稱為時變特征。雖然從單個狀態(tài)中可以獲知當前數(shù)據(jù)點的壓強取值在相鄰時間幀內(nèi)的基本變化趨勢和變化速率，但這只是瞬時特點。使用一段STC來定義特征，能夠追蹤壓強在一段時間內(nèi)的變化特點和描述爆炸過程中包含的物理規(guī)律。例如，S→A→E表示數(shù)據(jù)點的壓強從初始值開始連續(xù)兩個時間幀內(nèi)保持增長趨勢，對應(yīng)于數(shù)據(jù)點初次受到?jīng)_擊波影響并在兩幀內(nèi)均受到最外層沖擊波影響的情況。雖然STC越長所能獲取的信息就會越多，但是其中包含的變化特征會越復(fù)雜，導(dǎo)致視覺混亂和分析困難。綜合考慮以上因素，本文將長度為2的STC定義為時變特征（除去S→S）。排除S→S，符合圖5中轉(zhuǎn)換規(guī)律的時變特征總共有22種，這個數(shù)量的特征類型適合人類觀察和分析。

1.3? 傳輸函數(shù)

在傳輸函數(shù)設(shè)計過程中，包含過多手動設(shè)置的內(nèi)容，需要用戶對爆炸沖擊波的產(chǎn)生和傳播、傳輸函數(shù)的原理等有一定的知識儲備。本節(jié)提出了基于STG域的自動傳輸函數(shù)設(shè)計方法，主要分為顏色傳輸函數(shù)和不透明度傳輸函數(shù)兩部分。

1.3.1? 顏色傳輸函數(shù)

顏色傳輸函數(shù)將時變特征作為顏色分配的基本單位。顏色分配按時變特征中包含的狀態(tài)在完全狀態(tài)轉(zhuǎn)移鏈（complete state transition chain，CSTC）中的位置信息自動進行。CSTC將所有的狀態(tài)按照一定的轉(zhuǎn)移規(guī)則串聯(lián)在了一起。狀態(tài)在其中的位置信息將作為分配顏色的關(guān)鍵依據(jù)。數(shù)據(jù)點在壓強變化緩慢時，會長時間處于同一狀態(tài)或在鄰近狀態(tài)（目標區(qū)域中相鄰）之間轉(zhuǎn)移，在受到?jīng)_擊波的影響時會轉(zhuǎn)移到目標區(qū)域中相隔較遠的狀態(tài)。CSTC中按照盡可能減少時間梯度符號的反轉(zhuǎn)次數(shù)和轉(zhuǎn)移到鄰近狀態(tài)的原則選擇后一狀態(tài)，保證時間梯度符號一致的鄰近狀態(tài)對應(yīng)的顏色也相近。圖6展示了目標區(qū)域?qū)?yīng)的CSTC，并按照狀態(tài)在其中的先后順序標注了序號。

顏色傳輸函數(shù)采用HSV顏色空間，按照CSTC中狀態(tài)之間的有序關(guān)系，分別基于當前狀態(tài)和前一狀態(tài)確定色調(diào)（H）、飽和度（S）的取值。具體的顏色傳輸函數(shù)如圖7所示。為了讓繪制結(jié)果更加明亮，設(shè)定明度（V）為固定值1.0（V的取值為0.0～1.0）。相較于RGB模型，HSV模型與人類解釋顏色的方式一致，三個參數(shù)作用明確、相關(guān)性低。時變特征中包含了數(shù)據(jù)點在前一時間幀所處的狀態(tài)和當前時間幀所處的狀態(tài)。時變特征的當前狀態(tài)有7種（不包含狀態(tài)S）。首先按照當前狀態(tài)確定色調(diào)值，基于CSTC中的順序?qū)?種狀態(tài)線性映射到色調(diào)的0°～240°。然后按照前一狀態(tài)確定飽和度的取值，進一步區(qū)分當前狀態(tài)一致的一組時變特征。圖5中狀態(tài)的入度有三種情況： 1（當前狀態(tài)為A），4（當前狀態(tài)為E、F、G），3（當前狀態(tài)為B、C、D）。飽和度的分配策略按照入度的取值也分為三種，依據(jù)序號的大小關(guān)系線性映射到0.4～1.0。飽和度過低會導(dǎo)致顏色偏白，所以均設(shè)置了0.4的偏移值。

1.3.2? 不透明度傳輸函數(shù)

在爆炸數(shù)據(jù)中，最感興趣的部分為壓強變化劇烈的區(qū)域，壓強變化平緩的區(qū)域通常為不感興趣的部分。本節(jié)引入新變量φ來衡量壓強在相鄰時間幀之間的變化速率。

φ=arctan（GTP）（3）

通過反正切函數(shù)將φ的取值映射到了一個有限范圍內(nèi)，值域為（-π/2，π/2）。因為P的非負性，φ的實際取值為（-π/2，π/4）。不透明度傳輸函數(shù)設(shè)計為

α（φ）=σ2-σ1（θ1-π/2）2（φ+θ1）2+σ1-π2<φ<-θ10-θ1≤φ≤θ2σ2-σ1（π/4-θ2）2（φ-θ2）2+σ1θ2<φ<π4

θ1∈（0，π2），θ2∈（0，π4），σ1、σ2∈[0，1]（4）

φ的取值分為三個區(qū)間，通過θ1、θ2兩個參數(shù)調(diào)節(jié)邊界。θ1、θ2為分段曲線中兩個可交互設(shè)置的分界點。不透明度的非零取值為［σ1，σ2］，σ1、σ2分別為交互設(shè)置的最小和最大不透明度取值。

不透明度傳輸函數(shù)的示意圖如圖8所示，展示了不透明度α與變量φ的關(guān)系曲線。圖中：①和③段曲線均為開口朝上的拋物線，同一區(qū)間內(nèi)|φ|越小，曲線就越平緩，反之越陡峭；②段曲線為α取值恒為0的線段，|φ|越小，表示壓強的變化越緩慢。不透明度傳輸函數(shù)通過①③段映射曲線凸顯了變化劇烈的特征并弱化了變化緩慢的特征，通過②段映射曲線過濾了整個區(qū)間內(nèi)壓強變化最緩慢的特征。

2? 實驗結(jié)果

在配備有AMD Ryzen 5 3500X 6-Core Processor（3.59 GHz）和 NVIDIA GeForce GTX 1660 SUPER的Windows 10電腦上，使用Qt編寫了面向爆炸數(shù)據(jù)的交互界面，并驗證了基于STG二維空間的傳輸函數(shù)設(shè)計方法的有效性。實驗部分選擇了兩類典型場景的爆炸數(shù)據(jù)：a）維度為501×51×501，相鄰時間幀壓強差值最高達7個數(shù)量級的水下爆炸數(shù)據(jù)（在水中設(shè)置了9個爆炸點）；b）維度為902×202×802，相鄰時間幀的壓強差值最高達9個數(shù)量級的城市爆炸數(shù)據(jù)。

2.1? 結(jié)果分析

在分析爆炸數(shù)據(jù)的過程中，為了探究各種不同的特征、方便調(diào)節(jié)參數(shù)和顯示關(guān)鍵信息，基于QT編寫了分析數(shù)據(jù)的軟件，如圖9所示。界面主要分為三個部分：第一個部分為數(shù)據(jù)的讀取、信息的顯示和交互信息的設(shè)置； 第二部分為散點圖和圖元選中區(qū)域的顯示； 第三部分為繪制窗口，繪制的結(jié)果對應(yīng)散點圖中圖元選中區(qū)域。用戶可以選擇不同的圖元形狀，并通過鼠標移動圖元或者拖拽控制點修改圖元的大小和形狀，在繪制窗口實時觀察選中區(qū)域?qū)?yīng)的繪制結(jié)果。

2.1.1? 算法實現(xiàn)

本文應(yīng)用光線投射（ray casting）算法［24］實現(xiàn)體繪制，主要對傳輸函數(shù)部分進行了修改。傳輸函數(shù)的實現(xiàn)步驟如下：

輸入：采樣點在前一幀和當前幀的壓強和時間梯度。

輸出：采樣點對應(yīng)的顏色（h，s，v）和不透明度α。

step1? 初始化顏色為（0.0，0.0，0.0），不透明度為0.0。

step2? 確定時變特征。

step2.1? 根據(jù)表1，確定前一時間幀的狀態(tài)lastState和當前時間幀的狀態(tài)curState。

step2.2? 串聯(lián)前一狀態(tài)和當前狀態(tài)，得到時變特征lastState→curState。

step2.3? 若時變特征為S→S，直接返回初始顏色和不透明度。

step3? 計算顏色。

step3.1? 亮度v設(shè)置為1.0。基于時變特征查詢圖7得到飽和度s和色相h的取值。

step4? 計算不透明度α。

step4.1? 根據(jù)式（3）計算得到當前時間幀的壓強變化速率φ。

step4.2? 將φ代入式（4），計算得到對應(yīng)的不透明度α。

在傳輸函數(shù)的執(zhí)行過程中，執(zhí)行時間主要分布在step2.1、step3.1中的圖表查詢和step4.1、step4.2中的公式計算兩部分。其中圖表查詢部分為3次索引表的直接查詢，其查詢時間復(fù)雜度均為O（1）；計算部分為式（3）（4）的直接計算，其復(fù)雜度也為O（1）。因此本文提出的傳輸函數(shù)的整體計算復(fù)雜度為O（1），與其他同類型傳輸函數(shù)［25，26］相比，計算復(fù)雜度一致。在現(xiàn)代GPU中，整個算法耗時極少。

2.1.2? 時變特征解析

在繪制結(jié)果中，顏色和時變特征之間直觀的對應(yīng)關(guān)系可以幫助分析不同爆炸場景中包含的物理規(guī)律。圖10展示了兩類爆炸數(shù)據(jù)繪制結(jié)果中部分顏色和時變特征的對應(yīng)關(guān)系，以11種時變特征為例（總共有22種，如圖7所示）。圖中序號1、2、9標注的時變特征在空間上總是相鄰的。在爆炸中心向外的方向上，這三種時變特征的次序總是9、2、1。在最外層沖擊波的影響下，數(shù)據(jù)點的壓強首先會迅速上升。因為沖擊波存在厚度，該區(qū)域壓強持續(xù)上升，然后迅速下降。圖10中上面兩張子圖對應(yīng)城市爆炸場景，爆炸后沖擊波近似球狀向四周傳播。其中向地面?zhèn)鞑サ臎_擊波發(fā)生反射，傳播方向發(fā)生了變化，形成了左上圖中序號4、5對應(yīng)的時變特征。沖擊波在建筑群中傳播時，沖擊波近地的部分與建筑群發(fā)生碰撞，形成了右上圖中序號3、5、11對應(yīng)的時變特征。下面兩張子圖對應(yīng)水下爆炸場景，不同爆炸點形成的沖擊波發(fā)生交匯，形成了右下圖中序號5、10對應(yīng)的時變特征。因為繪制結(jié)果中存在顏色疊加和光照效果，所以紅框內(nèi)的顏色（參見電子版）和標注的顏色會存在小范圍偏差，這種誤差可忽略。

2.1.3? 時變特征過濾

為了過濾不感興趣的時變特征，幫助用戶更清晰地觀察反射波和繞射波在建筑群中的傳播規(guī)律，不透明度傳輸函數(shù)中提供了θ1、θ2兩個可交互調(diào)節(jié)的參數(shù)。當最外層沖擊波在建筑群中傳播時，與建筑發(fā)生碰撞形成的反射波和繞射波會互相疊加和追趕，向與傳播方向相反的空中擴散。這導(dǎo)致從空中向下的角度無法觀察到反射波和繞射波在建筑群中的傳播情況。反射發(fā)生在建筑物附近，形成的反射波擴散到空中后，能量已經(jīng)發(fā)生了衰減。不透明度傳輸函數(shù)基于壓強變化速率過濾空中的沖擊波，通過調(diào)節(jié)θ1、θ2設(shè)置過濾范圍。圖11展示了6組不同θ1、θ2取值對應(yīng)的城市爆炸繪制結(jié)果，從左到右θ1增大，從上到下θ2增大。圖（a）完全無法觀察到建筑群附近的沖擊波傳播情況；隨著θ1值增大，繪制結(jié)果中冷色對應(yīng)的區(qū)域明顯減少，如圖（a）（b）之間；隨著θ2值增大，繪制結(jié)果中暖色對應(yīng)的區(qū)域明顯減少，如圖（b）（c）之間。

2.2? 值域上區(qū)分特征的能力

本文方法可以在局部取值范圍內(nèi)清晰地區(qū)分多種特征。傳統(tǒng)的可視化方法通常將具有相近標量值的區(qū)域定義為相同或相似特征，導(dǎo)致局部取值范圍對應(yīng)的繪制結(jié)果中顏色總是趨于一致，難以進行進一步的區(qū)分。數(shù)據(jù)可視化中普遍使用的分析工具有Voreen［25］和ParaView［26］，兩個平臺都提供了一些可選的自動傳輸函數(shù)方法。本節(jié)選取了基于標量的一維全局線性傳輸函數(shù)BS-GLTF（基于標量的一維全局線性傳輸函數(shù)）和全局非線性傳輸函數(shù)BS-GNLTF（基于標量的一維全局非線性傳輸函數(shù)）兩種自動傳輸函數(shù)方法作為對比，驗證本文方法在值域上區(qū)分特征的能力。

一維全局線性傳輸函數(shù)作為最簡單的模式，廣泛應(yīng)用在數(shù)據(jù)可視化中。首先獲取壓強在整個爆炸過程中的取值范圍，將壓強線性映射到［0，1］。

w=P-PminPmax-Pmin（5）

其中：Pmax和Pmin分別為壓強的全局最大值和最小值。然后，基于w值設(shè)置不透明度和顏色的色調(diào)。

α=wρc.h=240（1-w）（6）

其中：α和c分別表示不透明度和顏色；ρ為預(yù)設(shè)的不透明度最大值；h表示色調(diào)。飽和度和明度設(shè)置為固定值。

圖像增強領(lǐng)域為了增強低照明條件下捕獲的數(shù)字圖像，會首先采用非線性映射關(guān)系改變強度的全局分布［27］。傳輸函數(shù)設(shè)計與圖像增強的目標存在一致性。BS-GNLTF方法對應(yīng)的映射曲線如圖12所示。

水下爆炸數(shù)據(jù)的壓強變化最高跨越了7個數(shù)量級。相較于整個壓強取值范圍，±105的取值范圍可映射的顏色范圍非常小。基于水下爆炸數(shù)據(jù)，圖13展示了兩組不同的壓強范圍對應(yīng)的繪制結(jié)果。BS-GLTF對應(yīng)的繪制結(jié)果整體不透明度偏低，顏色趨于一致。相較于前者，BS-GNLTF對應(yīng)的繪制結(jié)果在不透明度上有了更明顯的區(qū)分。但是，這兩種方法對應(yīng)的繪制結(jié)果均顏色單一，無法通過顏色對特征進行區(qū)分。圖（a）中本文方法對應(yīng)的繪制結(jié)果，中心區(qū)域包含了4種以上明顯可區(qū)分的顏色； 圖（b）中對應(yīng)的繪制結(jié)果，從內(nèi)到外包含了2、3種明顯可區(qū)分的顏色。不同顏色的區(qū)域之間界限清晰，不存在顏色上的混淆。本文方法得到的繪制結(jié)果，無論是在數(shù)據(jù)點分布密集、特征豐富的壓強范圍（1×105±105內(nèi)），還是在數(shù)據(jù)點分布稀疏、包含特征較少的壓強范圍（2×106±105內(nèi)），均能通過顏色清晰地辨別多種特征。

2.3? 時域上區(qū)分特征的能力

爆炸數(shù)據(jù)屬于時變數(shù)據(jù)，屬性取值范圍和分布會隨時間發(fā)生變化。本文設(shè)計了一種與具體屬性取值范圍和分布無關(guān)的全局傳輸函數(shù)，保證了時變特征與顏色的對應(yīng)關(guān)系在時間維度上的一致性。本節(jié)選取了同樣結(jié)合時間相關(guān)信息的BTG-GNLTF方法（基于時間梯度的一維全局非線性傳輸函數(shù)）作為對比，驗證本文傳輸函數(shù)方法在時域上區(qū)分特征和展現(xiàn)特征演變規(guī)律的有效性。

Akiba等人［3］提出了基于時間梯度的時間直方圖，用于識別和區(qū)分隨時間變化高度活躍的區(qū)域。本文方法在BS-GNLTF的基礎(chǔ)上用時間梯度替換標量屬性，設(shè)計了一種分段式的非線性傳輸函數(shù)。選取五個關(guān)鍵值，即0、最大值曲線中的最小值gs min和最大值gmax、最小值曲線中的最大值gs max和最小值gmin，作為傳輸函數(shù)的分界點，具體的映射曲線如圖14所示。

本文方法保證了時變特征和光學屬性的映射關(guān)系在時間維度上的一致性，實現(xiàn)了對時變特征的連續(xù)跟蹤。基于水下爆炸數(shù)據(jù)，圖15展示了BS-GNLTF、BTG-GNLTF和本文方法對應(yīng)的繪制結(jié)果，均只截取了單個爆炸位置的繪制結(jié)果。圖（a）中幀序1對應(yīng)的繪制結(jié)果可以明顯觀察到橙色的環(huán)狀特征，即最外層沖擊波，但是在幀序3和6對應(yīng)的繪制結(jié)果中，已經(jīng)無法捕捉到這個特征。圖（b）的三幀繪制結(jié)果均能觀察到顏色區(qū)別于中心區(qū)域的環(huán)狀特征，但是在幀序3和6的繪制結(jié)果中，該環(huán)狀特征的顏色和鄰近區(qū)域非常相近，且環(huán)狀特征和光學屬性的對應(yīng)關(guān)系發(fā)生了明顯的變化。最外層沖擊波通常攜帶著巨大的能量向四周傳播，但是這種能量會隨著傳播距離的增加發(fā)生衰減。BS-GNLTF和BTG-GNLTF都與屬性的具體取值范圍和分布緊密相關(guān)，無法在時域上連續(xù)跟蹤沖擊波的傳播特征。圖（c）中時變特征與顏色的對應(yīng)關(guān)系未隨時間發(fā)生變化，可以通過顏色清晰地觀察到時變特征的演變規(guī)律。

隨著傳播距離的增加，最外層沖擊波攜帶的能量會迅速衰減，極大地增加了中后期捕捉特征的難度。本文方法在沖擊波傳播的中后期仍然能準確地捕捉?jīng)_擊波的傳播特征。基于城市爆炸數(shù)據(jù)，圖16展示了BTG-GNLTF和本文方法對應(yīng)的兩組處于3個不同沖擊波傳播階段的繪制結(jié)果。BTG-GNLTF基于時間梯度的全局取值范圍進行顏色映射，整體呈冷色。圖（a）中幀序13對應(yīng)的繪制結(jié)果可以觀察到最外層沖擊波和與地面碰撞后形成的反射波。圖（a）中幀序52、91對應(yīng)的繪制結(jié)果只能觀察到前者，并且，在幀序91對應(yīng)的結(jié)果中，最外層沖擊波已經(jīng)非常稀薄。圖（b）中三幀繪制結(jié)果均能夠清晰地觀察到最外層沖擊波和殘留在建筑群中的繞射波和反射波。

分析了兩份爆炸數(shù)據(jù)對應(yīng)的繪制結(jié)果后得出結(jié)論：本文方法在時域上具備區(qū)分重要特征和穩(wěn)定跟蹤沖擊波傳播特征的能力。

3? 結(jié)束語

本文面向大規(guī)模爆炸數(shù)據(jù)提出了一種基于STG空間的自動傳輸函數(shù)設(shè)計方法。爆炸數(shù)據(jù)對應(yīng)的壓力場取值范圍和分布會隨時間發(fā)生劇烈變化，極大地增加了連續(xù)跟蹤沖擊波傳播特征的難度。為此，本文將沖擊波傳播特征編碼為一段狀態(tài)轉(zhuǎn)移鏈，提出了基于時變特征的顏色傳輸函數(shù)和基于壓強變化率的不透明度傳輸函數(shù)。在不透明度傳輸函數(shù)中設(shè)置了可交互的參數(shù)θ1、θ2，過濾不感興趣的時變特征。本文方法與壓力場的具體取值范圍和分布無關(guān)，因此它可以在局部壓強范圍內(nèi)區(qū)分多種特征和連續(xù)跟蹤沖擊波的傳播特征。將本文方法應(yīng)用在兩種經(jīng)典的爆炸場景中，均獲得了期望的效果。

本文方法暫時只應(yīng)用在了爆炸數(shù)據(jù)的可視化中，但是計算流體力學、燃燒等領(lǐng)域的計算模擬數(shù)據(jù)也屬于復(fù)雜的時變數(shù)據(jù)，具有與爆炸數(shù)據(jù)相同的數(shù)據(jù)特點。在未來的工作中，筆者希望將本文方法推廣應(yīng)用到更多種類的大規(guī)模復(fù)雜時變數(shù)據(jù)的可視化工作中。

參考文獻：

［1］Bai Zhihui，Tao Yubo，Lin Hai.Time-varying volume visualization：a survey［J］.Journal of Visualization，2020，23（5）：745-761.

［2］Kosara R，Bendix F，Hauser H.Time histograms for large，time-dependent data ［C］// Proc of the 6th Joint Eurographics-IEEE TCVG conference on Visualization.Goslar：Eurographics Association，2004：45-54.

［3］Akiba H，F(xiàn)out N，Ma K L.Simultaneous classification of time-varying volume data based on the time histogram［C］//Proc of the 8th Joint Eurographics/IEEE VGTC Conference on Visualization.Goslar：Eurographics Association，2006：171-178.

［4］Fang Zhe，Mller T，Hamarneh G，et al.Visualization and exploration of time-varying medical image data sets ［C］// Proc of Graphics Interface.New York：ACM Press，2007：281-288.

［5］Lee T Y，Shen H W.Visualizing time-varying features with TAC-based distance fields［C］//Proc of IEEE Pacific Visualization Symposium.Piscataway，NJ：IEEE Press，2009：1-8.

［6］Bai Zhihui，Tao Yubo，Lin Hai.FeatureFlow：exploring feature evolution for time-varying volume data［J］.Journal of Visualization，2019，22（5）：927-940.

［7］Silver D，Wang X.Tracking and visualizing turbulent 3D features［J］.IEEE Trans on Visualization and Computer Graphics，1997，3（2）：129-141.

［8］Muelder C，Ma K L.Interactive feature extraction and tracking by utilizing region coherency［C］//Proc of IEEE Pacific Visualization Symposium.Piscataway，NJ：IEEE Press，2009：17-24.

［9］Soler M，Plainchault M，Conche B，et al.Lifted Wasserstein matcher for fast and robust topology tracking ［C］//Proc of the 8th IEEE Symposium on Large Data Analysis and Visualization.Piscataway，NJ：IEEE Press，2018：23-33.

［10］姚成寶，王宏亮，浦錫鋒，等.空中強爆炸沖擊波地面反射規(guī)律數(shù)值模擬研究［J］.爆炸與沖擊，2019，39（11）：24-31.（Yao Chengbao，Wang Hongliang，Pu Xifeng，et al.Numerical simulation of intense blast wave reflected on rigid ground［J］.Explosion and Shock Waves，2019，39（11）：24-31.）

［11］Soukup J，Klimenda F，Skocˇilas J，et al.Finite element modelling of shock wave propagation over obstacles［J］.Manufacturing Techno-logy，2019，19（3）：499-507.

［12］Wang Yan，Wang Hua，Cui Cunyan，et al.Investigating different grounds effects on shock wave propagation resulting from near-ground explosion［J］.Applied Sciences，2019，9（17）：3639.

［13］曹濤，孫浩，周游，等.近地爆炸沖擊波傳播特性數(shù)值模擬與應(yīng)用［J］.兵器裝備工程學報，2020，41（12）：187-191.（Cao Tao，Sun Hao，Zhou You，et al.Numerical simulation and application of propagation characteristics of shock wave near ground explosion［J］.Journal of Ordnance Equipment Engineering，2020，41（12）：187-191.）

［14］Bremer P T，Weber G，Pascucci V，et al.Analyzing and tracking bur-ning structures in lean premixed hydrogen flames［J］.IEEE Trans on Visualization and Computer Graphics，2010，16（2）：248-260.

［15］Saikia H，Weinkauf T.Global feature tracking and similarity estimation in time-dependent scalar fields ［J］.Computer Graphics Forum，2017，36（3）：1-11.

［16］Dutta S，Shen H W.Distribution driven extraction and tracking of features for time-varying data analysis ［J］.IEEE Trans on Visualization and Computer Graphics，2016，22（1）：837-846.

［17］Wang K C，Wei T H，Shareef N，et al.Ray-based exploration of large time-varying volume data using per-ray proxy distributions［J］.IEEE Trans on Visualization and Computer Graphics，2020，26（11）：3299-3313.

［18］劉力.面向時變體數(shù)據(jù)的特征可視化方法［J］.中國圖象圖形學報，2022，27（4）：1302-1313.（Liu Li.A feature visualization method for time-varying volume data［J］.Journal of Image and Graphics，2022，27（4）：1302-1313.）

［19］Akiba H，Wang C，Ma K L.AniViz：a template-based animation tool for volume visualization ［J］.IEEE Computer Graphics and Applications，2010，30（5）：61-71.

［20］Elvins T T.A survey of algorithms for volume visualization［J］.ACM SIGGRAPH Computer Graphics，1992，26（3）：194-201.

［21］ereda P，Vilanova A，Gerritsen F A.Automating transfer function design for volume rendering using hierarchical clustering of material boundaries［C］//Proc of the 8th Joint Eurographics/IEEE VGTC Conference on Visualization.Goslar：Eurographics Association，2006：243-250.

［22］朱奭，常晉義.一種改進的基于 CUDA 的紋理映射和光線投射結(jié)合的體繪制算法 ［J］.計算機應(yīng)用研究，2015，32（6）：1884-1887.（Zhu Shi，Chang Jinyi.Improved algorithm of volume rendering combined texture mapping with ray casting based on CUDA［J］.Application Research of Computers，2015，32（6）：1884-1887.）

［23］Cai L，Nguyen B P，Chui C K，et al.A two-level clustering approach for multidimensional transfer function specification in volume visua-lization［J］.The Visual Computer，2017，33（2）：163-177.

［24］Kruger J，Westermann R.Acceleration techniques for GPU-based vo-lume rendering［C］//Proc of the 14th IEEE Visualization.Washington DC：IEEE Computer Society，2003：38.

［25］Meyer-Spradow J，Ropinski T，Mensmann J，et al.Voreen：a rapid-prototyping environment for ray-casting-based volume visualizations ［J］.IEEE Computer Graphics and Applications，2009，29（6）：6-13.

［26］Ahrens J，Geveci B，Law C.ParaView：an end-user tool for large data visualization ［M］//Visualization Handbook.[S.l]：Butterworth-Heineman，2005：717-731.

［27］Tao Li，Asari V K.Adaptive and integrated neighborhood-dependent approach for nonlinear enhancement of color images ［J］.Journal of Electronic Imaging，2005，14（4）：043006.