# Towards Power Efficient DNN Accelerator Design on Reconfigurable Platform

Rourab Paul1, Sreetama Sarkar2, Suman Sau3, Sanghamitra Roy4, Koushik
Chakraborty5, Amlan Chakrabarti6 Computer Science & Engineering, Siksha O
Anusandhan, India1,
Electrical and Computer Engineering Technical University Munich, Germany2,
Computer Science & Information Technology, Siksha ’O’ Anusandhan, India3,
Dept. Electrical and Computer Engineering, Utah State University, Logan,
USA4,5,
School of IT, University of Calcutta, India6
mail<EMAIL_ADDRESS>
###### Abstract

The exponential emergence of Field Programmable Gate Array (FPGA) has
accelerated the research of hardware implementation of Deep Neural Network
(DNN). Among all DNN processors, domain specific architectures, such as,
Google’s Tensor Processor Unit (TPU) have outperformed conventional GPUs.
However, implementation of TPUs in reconfigurable hardware should emphasize
energy savings to serve the green computing requirement. Voltage scaling, a
popular approach towards energy savings, can be a bit critical in FPGA as it
may cause timing failure if not done in an appropriate way. In this work, we
present an ultra low power FPGA implementation of a TPU for edge applications.
We divide the systolic-array of a TPU into different FPGA partitions, where
each partition uses different near threshold (NTC) biasing voltages to run its
FPGA cores. The biasing voltage for each partition is roughly calculated by
the proposed static schemes. However, further calibration of biasing voltage
is done by the proposed runtime scheme. Four clustering algorithms based on
the minimum slack value of different design paths of Multiply Accumulates
(MACs) study the partitioning of FPGA. To overcome the timing failure caused
by NTC, the MACs which have higher minimum slack are placed in lower voltage
partitions and the MACs have lower minimum slack path are placed in higher
voltage partitions. The proposed architecture is simulated in a commercial
platform : $Vivado$ with Xilinx $Artix-7$ FPGA and academic platform VTR with
22nm, 45nm, 130nm FPGAs. The simulation results substantiate the
implementation of voltage scaled TPU in FPGAs and also justifies its power
efficiency.

###### Index Terms:

FPGA partition, Low Power, TPU, Voltage Scaling

##  I Introduction

The configurable logic block (CLB) and switch matrix of FPGAs are power-
hungry, which makes FPGAs energy inefficient compared to ASICs. Recently many
researchers [1], [2] have reported CPU-FPGA based hybrid data center
architectures which provide hardware acceleration facility for Deep Neural
Networks (DNNs). Despite power inefficiency, FPGA becomes popular in the
Cloud-Scale acceleration architecture due to its specialized hardware and the
economic benefits of homogeneity. Therefore, reducing power in FPGA for DNN
applications becomes a very relevant topic of research. B Salami et al. [3]
has studied the timing failure vs biasing voltage of DNN implementation in
FPGA. They have underscaled biasing voltage $V_{ccint}$ of the entire FPGA to
increase the power efficiency of Convolutional Neural Network (CNN)
accelerator by a factor of 3. A single $V_{ccint}$ for the entire FPGA might
not be the most power efficient solution. Partitioning an FPGA according to
the slacks and feeding different biasing voltages for different partitions can
cause further reduction of power for CNN implementations. In [4], the authors
have implemented a systolic array using near threshold (NTC) biasing voltage
in ASIC, which can predict the timing failure of multiplier-accumulators
(MACs) placed inside the systolic array of TPU. The prediction of timing
failure is based on $Razor$ flipflop [5]. Higher fluctuation of input bits
increases the possibility of timing failure in NTC condition. In [4], once the
timing failure of a MAC is predicted by its internal $Razor$ flipflop, the
biasing voltage of the MAC is boosted up.

Targeting FPGA based DNN applications [1], our work investigates voltage
scaling techniques of TPU in the FPGA platform. Different $V_{ccint}$ for each
of the MACs in a systolic array will be an absurd implementation for FPGA,
therefore this work partitions FPGA floor according to the minimum slack value
of design paths of MACs. Each partition consists of a group of MACs having
similar minimum slacks. Each partition is connected with different
$V_{ccint}$. The proposed methodology abstracts the synthesis timing report
from the $Vivado$ and $VTR$. tool. In a synthesized design, the $Vivado$ and
$VTR$ timing engine estimate the net delays of paths based on connectivity and
fanout. The clustering algorithms create clusters or groups based on the
minimum slack of MACs. The clusters consist of MACs which have lower minimum
slacks are placed in FPGA partitions with higher $V_{ccint}$ and the clusters
of MACs which have higher minimum slacks are placed in FPGA partitions with
lower $V_{ccint}$. Here the $V_{ccint}$ provides power to a FPGA core. The
tuning of $V_{ccint}$ with slack is done by unique $static-runtime$ strategy.
The circuit level challenges on the implementation of voltage scaling in FPGA
platform are beyond the present scope of our article. However, the feasibility
of implementing the necessary hardware for voltage scaling support is evident
considering the successful implementations in other ASIC technologies. As is
unavailable in current FPGAs we have simulated the design for the validation
of the claim. The contribution of the paper is as follows:

  * •

This paper proposes a new CAD flow to create voltage scaled TPU in FPGA based
platforms considering the trade off of circuit delay against biasing voltage.

  * •

The proposed algorithm divides the systolic array of TPU into different
partitions. Each partition will have different $V_{ccint}$. The $V_{ccint}$ in
different partitions is scaled against the different minimum slacks of
different MACs.

  * •

The calibration of $V_{ccint}$ of different partitions is done by the proposed
$runtime$ and $static$ schemes.

The organization of the article is as follows: Sec. II outlines our background
of FPGA environment. The working principle of Razor flipflop to detect runtime
timing failure is discussed in II-E. The methodology of the proposed work is
described in Sec. III. Sec. IV discusses the clustering algorithms. Result,
implementation and conclusion are organized in Sec. V and Sec. VI
respectively.

##  II Background: FPGA Environment

The proposed scheme has been simulated in the both commercial and academic CAD
tools. In our first approach, we have used the Xilinx $Vivado$ tool with
Artix-7 FPGA. For the sake of more accurate power data, the proposed voltage
scaled architecture is also simulated in VTR tool flow with 22nm, 45nm and
130nm academic FPGAs.

###  II-A Vivado Environment

A typical Xilinx FPGA in $Vivado$ environment has 3 conventional steps such as
synthesis, implementation and bit file generation whereas the adopted tool
flow of the proposed partitioned FPGA is divided into two environments: (i)
$Vivado$ Environment for synthesis, implementation and bit file generation and
(ii) Python Environment for clustering similar slacks. The entire tool flow is
shown in Fig. 1. The $Vivado$ environment is involved with 3 sub-steps stated
below:

####  II-A1 Synthesis

$Vivado$ synthesis process transforms register transistor logic (RTL) to gate
level representation. The synthesis process generates delays of all possible
paths of the design. The timing report of the synthesis process contains 12
information such as name of the path, slack value, level, high fanout, path
from, path to, total delay of path, logic delay, net delay, time requirement
source clock and destination clock as shown in Table I. It is to be noted that
the estimation of the slacks of each logic block is at a high level. The
actual timing behavior of the design depends on the net delays after placement
and routing.

TABLE I: A Fragment of Timing Report from Synthesis for 100 MgHz Clock

Name | Slack | Levels |  High Fanout | From | To |  Total Delay |  Logic Delay |  Net Delay |  Requir ement |  Source Clock |  Destination Clock
---|---|---|---|---|---|---|---|---|---|---|---
Path 1 | 5.34 | 8 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[16]/D | 4.37 | 2.80 | 1.57 | 10.00 | clk | clk
Path 2 | 5.49 | 8 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[15]/D | 4.40 | 2.83 | 1.57 | 10.00 | clk | clk
Path 3 | 5.52 | 9 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[14]/D | 4.36 | 2.89 | 1.47 | 10.00 | clk | clk
Path 4 | 5.59 | 9 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[13]/D | 4.30 | 2.83 | 1.47 | 10.00 | clk | clk
Path 5 | 5.78 | 7 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[12]/D | 4.10 | 2.54 | 1.57 | 10.00 | clk | clk
Path 6 | 5.83 | 7 | 8 | GEN_REG_I[0].GEN_REG_J[1].uut/prev_activ_reg[1]/C | GEN_REG_I[1].GEN_REG_J[1].uut/sig_mac_out_reg[11]/D | 4.05 | 2.49 | 1.57 | 10.00 | clk | clk

####  II-A2 Implementation

The $Vivado$ Implementation process is a timing-driven flow that transforms a
logical netlist and constraints (Xilinx Design Constraints format) into a
placed and routed design to make it ready for the bitstream generation
process. In our proposed tool flow, the logical netlist is provided by the
$Vivado$ synthesis process but the Xilinx Design Constraints (XDC) is
generated by a python script. The clustered MACs are considered for placing in
a specific location on the FPGA floor.

####  II-A3 Bit File Generation

Once the placement and routing are completed by the implementation process,
the flow generates bitstream of the systolic array. The Xilinx bitstream
generation program produces a bitstream for the Xilinx device configuration.
If there is any requirement of a processor, the design may include software
application data.

Figure 1: Vivado Tool Flow Figure 2: TPU Architecture

###  II-B VTR Environment

In a commercial CAD environment biasing voltage is fixed. The Verilog to
Routing (VTR) [6] tool is an open source academic CAD tool flow for FPGA
architecture which allows voltage scaling technology. The VTR contains 3
separate tools such as Odin II [7], ABC [8] and VPR [9].

####  II-B1 Synthesis

The synthesis process of the proposed VTR tool flow is processed by Odin II
and ABC. Odin II elaborates and synthesizes HDL into FPGA architectural
primitives like FFs, multipliers, adders. Thereafter the circuit logic is
handled by ABC to perform technology independent logic optimizations, then
technology maps the soft logic to LUTs. The information in the timing report
generated by ABC is similar to the Vivado synthesis report. The different
slack value of different design paths in synthesis report is used in cluster
algorithms.

####  II-B2 Implementation

The VPR [9] tool is a part of VTR flow which is used for the physical
implementation of the circuit on the target FPGA architecture along with
Synopsys Design Constraints File (sdc). In VTR flow, the logical netlist is
provided by Odin II and ABC synthesis process but the SDC is generated by a
python script. The clustered slack values generated by the python script are
considered for placing the logic paths in a specific location on the FPGA
floor. At the end, VPR analyzes the circuit implementation to generate area,
speed and power data, and a post-implementation netlist. Many commercial CAD
tools like Titan Flow uses Intel’s Quartus and Yosys used VPR for logic
synthesis, optimization and technology mapping.

There is a possibility that after the partitioning of the systolic array,
delays of design paths from the implementation process may differ from delays
of design paths from the synthesis process. If the minimum slack of any MAC is
changed due to this partitioning, the entire design needs to re-cluster based
on the new minimum slacks of MACs. Fig. 4 and Fig. 5 report the differecnes of
delays of 100 worst design paths of synthesis process and implementation
(after partition) process. Fig. 4 and Fig. 5 shows that partitioning process
does not effect design paths significantly. The little proportional changes in
the delays of design paths could not effect the minimum slack of MACs in each
partiontion. Therefore, partitioning based on minimum slack of MACs in each
partiontion will effective

Figure 3: VTR Tool Flow Figure 4: 100 worst Set up paths in Vivado Figure 5:
100 worst Hold paths in Vivado

###  II-C Python Environment

The contribution of the paper lies in augmenting the standard FPGA design tool
flow by incorporating a python-based environment, which consists of a script
to run three subsequent processes such as choice of $Clustering~{}Algorithms$,
$Cluster~{}Generation$ and $Constraint~{}Generation$.

####  II-C1 Choice of Clustering Algorithms

A clustering algorithm suited to the requirements is chosen at this step. As
stated in Sec. III, this paper investigates 4 commonly-used clustering
algorithms such as Hierarchical, K-means, Mean-shift and DBSCAN.

####  II-C2 Cluster Generation

We have assumed that the FPGA is divided into a few partitions and each
partition has a different biasing voltage $V_{ccint}$. The clustering
algorithms create few groups of MACs. The MACs having similar minimum slacks
form a group and they are placed in the same FPGA partition.

####  II-C3 Constraint Generation

Xilinx uses a constraint file format (XDC) to specify the coordinates of
different paths of the proposed systolic array. The XDC file is generated by
the python script.

###  II-D Clustering MACs based on their Minimum Slacks

The idea of voltage scaling for partitioned systolic array was initially based
on the slacks generated from the synthesis report. The slack based clustering
can group different similar design paths belonging to different MACs which may
be placed in the different physical locations of FPGA floor by the placement
and routing algorithm. For the slack based design path partitioning approach,
intervention of the proposed tool’s script is far more as compared to the
existing placement and routing process of existing EDA tools. As a result, the
timing parameters reported by the synthesis process are varied significantly
after the placement and routing process of existing EDA tools at
implementation level. For 4 partitions, 16 $\times$ 16 systolic array, the
Vivado tool generates 6.23 ns critical path. The same design gives 11.93 ns
critical path after placement routing which is almost two times the critical
path generated from the synthesis report. We have noticed that the placement
and routing process of slack based partitioning of 64 $\times$ 64 systolic
array takes 10 to 14 hours in i5, 8GB Linux platform. Later, instead of
clustering design paths based on slack, clustering is performed on MACs using
their minimum slack values. We find clustering MACs based on their minimum
slack is reasonable and better compare to the previous one for the below
reasons:

  * •

For the clustering of MACs based on their minimum slack, the intervention of
vendor’s technology dependent placement and routing algorithm is far more
compared to the previous idea. As a result, the critical path variation in
synthesis and implementation process is very minimum.

  * •

Placing all design paths in constraints file is much complicated compare to
placing entire MACs in constraint file.

  * •

The routing of wires on the FPGA floor is comparatively simpler for the MAC
clustering based on their minimum slacks.

###  II-E Razor Flipflop

Razor flipflop can be implemented in FPGA [5] by inserting a shadow flipflop
which is running by a delayed clock. We have assumed that a circuit register
$R$ is lying at the end of one or more timing paths originating from any of
the source registers. The shadow register $S$ samples the same data as R but
on a delayed clock $DCLK$ which is lagged by $T_{del}$ from the main clock
$CLK$. Any data that arrives after R samples but before S samples will cause a
discrepancy between the two registers that is detected by the error flag $F$.
In each MAC unit, this razor flipflop is placed. The multiplication and
addition process in each MAC is computed using the rising edge of $CLK$ and
$DLK$. The $CLK$ driven output of the multiplication and addition process of
each MAC is stored in $R$ register. The $DCLK$ driven output of multiplication
and addition process of each MAC is stored in shadow register $S$. The
inclusion of razor doubles the number of multiplier and adder required for the
systolic array but it can detect runtime failure occured in MACs due to the
near threshold biasing voltage. The timing diagram of the razor is shown in
Fig. 6.

Figure 6: Timing Diagram of Fault Detection

##  III Hybrid Configuration: Static & Runtime Schemes

To mitigate timing failure issue in the critical voltage region, we adopted
two sequential schemes such as (i)Static scheme which is involved with FPGA
partitioning and rough $V_{ccint_{i}}$ estimation depending on the FPGA
technology. (ii) Runtime scheme to calibrate suitable $V_{ccint_{i}}$ for each
partition of FPGA using $Razor$ flipflop. Each partition of FPGA consists a
group of MACs. All the groups of MACs form a systolic array of the TPU. Apart
from systolic array, TPU has memory to store active and weight inputs, PCI
interface, controlling circuitry etc. The architecture of TPU is shown in Fig.
2.

###  III-A Static Scheme

The proposed staic scheme works on $Vivado$ and $Python$ environments. As
shown in Fig. 1, synthesis is the first step of the proposed tool flow, which
takes a netlist of complex logic blocks (CLBs) of systolic array generated
from $Vivado$ tool. This netlist from the synthesis report is generated after
technology mapping and packing stages which contain time slacks of all the
possible paths of the systolic array. The proposed approach considers only
nodes along paths because (i) the nodes along the path have data dependencies,
which should be placed in the same FPGA partition even without considering the
voltage scaling [10]. (ii) The slack values of the nodes along paths are
usually close to each other. The second step of the proposed methodology is
involved with the choice of the clustering algorithm and cluster generation.
As stated in Sec. IV, the four clustering algorithms such as Hierarchy,
K-Mean, Mean-Shift and DBSCAN create multiple cluster of MACs with the paths
available in the synthesis report. Even for the same number of clusters,
different algorithms classify the data-points slightly differently.

The primary concern is to identify clusters of MACs, which can share the
minimum slacks available across the other MACs. Even for the same number of
clusters, different algorithms classify the data-points slightly differently.
Unlike K-means algorithm, the Hierarchical, Mean-Shift and DBSCAN do not need
the number of clusters to be specified beforehand. DBSCAN is found to perform
the best in this case as it groups together data-points close by, has a
reasonable time complexity and can also identify outliers. Hence, clustered
paths returned by DBSCAN are chosen for subsequent simulations.

Once the number of clusters is fixed we need to decide the voltage values of
different FPGA partitions. In Fig. 7., we illustrate 3 voltage regions in an
FPGA, which is also supported by the research work in [3]. The voltage below
FPGA crashing voltage $V_{crash}$ causes timing failure, which reduces the DNN
accuracy near to zero. The region between minimum voltage $V_{min}$ and
nominal voltage $V_{nom}$ is called guard band region where the DNN accuracy
will be 100% but power efficiency will be the least. In the critical region,
the closer the voltage is to $V_{crash}$, higher is the power efficiency and
lower the DNN accuracy. Similarly, if $V_{ccint}$ is closer to $V_{min}$ in
the critical region, the power efficiency decreases and DNN accuracy
increases. In our proposed architecture, we assume the operating voltage range
for the systolic array is $V_{crash}$ to $V_{min}$. If we have $n$ clusters
computed by the chosen clustering algorithm we need $n$ partitions in FPGA.
The primary $V_{ccint}$ estimation for each FPGA partition is computed by
Algorithm 1. In Xilinx FPGA, the coordinates of circuit components are
specified by two slice parameters $(X_{i},Y_{j})$. Each FPGA partition has
range of these coordinates. The clustered MACs are placed in same FPGA
partition by mentioning the slice parameters $(X_{i},Y_{j})$.

Figure 7: Voltage behaviour for $V_{ccint}$

In the third step of the proposed methodology, each clustered path computed by
the clustering algorithms is placed in a particular FPGA partition, which is
restricted by specific $X_{i},Y_{j}$ ranges. This restriction is done in the
xdc file during $Generate~{}Constraint~{}File$ process.

The rough $V_{ccint}$ calculation is done by Static Voltage Scaling algorithm
shown in Algorithm 1 which calculates a stepping voltage $V_{s}$ from
$V_{min}$ and $V_{crash}$. Thereafter, the $V_{ccint_{i}}$ of the $i^{th}$
partition is calculated based on the stepping voltage $V_{s}$. The Static
Voltage Scaling algorithm distributes $V_{ccint}$s for $i$ number of
partitions for the range started from $V_{crash}$ to $V_{min}$.
Algorithm 1 Static Voltage Scaling 0: $V_{ccint}$, $V_{min}$, $V_{crash}$ &
$n$ 1: $V_{s}$=$\frac{V_{min}-V_{crash}}{n}$ 2: $V_{l}=V_{crash}$ 3: for i=0
to n-1 do 4: $V_{ccint_{i}}$ = $\frac{V_{l}+V_{l}+V_{s}}{2}$ 5:
$V_{l}=V_{l}+V_{s}$ 6: end for Algorithm 2 Runtime Voltage Scaling 0:
$V_{ccint}$, $V_{s}$ 1: for i=0 to n-1 do 2: if $timing\\_fail-part-i==1$ then
3: $V_{ccint_{i}}$ = $V_{ccint_{i}}+V_{s}$ 4: else 5: $V_{ccint_{i}}$ =
$V_{ccint_{i}}-V_{s}$ 6: end if 7: end for

###  III-B Runtime Scheme

The $V_{ccint_{i}}$ of the $i^{th}$ FPGA partition calculated by Algorithm 1
is calibrated to $V_{ccint_{i}}$ pin of the $i^{th}$ FPGA partition. The
calculation of $V_{ccint_{i}}$ by Algorithm 1 is based on the number of
partitions $n$ and the critical voltage region $V_{min}-V_{crash}$ which
solely depends on the type of FPGA technology. However, the appropriate
$V_{ccint_{i}}$ of the $i^{th}$ FPGA partition should also depend on the
minimum slack values of MACs of that partition. At the static strategy we just
calculate a rough estimation of $V_{ccint_{i}}$ where as the runtime strategy
calibrates $V_{ccint_{i}}$ according to the runtime timing failure of the
systolic array. In the runtime scheme we used one of the most popular runtime
timing error detection scheme, $Razor$, which uses double sampling flipflop to
detect timing violation of pipeline stages. The $Razor$ flipflop is connected
with every MACs of the systolic array to indicate its the timing failure. Each
MAC has a timing failure flag which is controlled by the $Razor$ flipflop. If
any timing failure flag of any MAC placed in the $i^{th}$ FPGA partition is
high, the $V_{ccint_{i}}$ of that $i^{th}$ FPGA partition will be increased by
one step. If all the timing failure flags of all MACs placed in the $i^{th}$
FPGA partition is low, the $V_{ccint_{i}}$ of that $i^{th}$ FPGA partition
will be decreased by one step. Before starting the actual run of the proposed
systolic array, if we have trial run, all the $V_{ccint_{i}}$ of all
partitions will be tuned accurately by this $runtime$ process. The voltage
boosting circuit can be implemented externally following the technique
proposed in [11].

Figure 8: Example : Partitioned FPGA, n=4

In Fig. 8, we have shown that the cluster algorithm partitions the FPGA into 4
islands. The static scheme as stated in Sec. III-A calculates 4
$V_{ccint_{i}}$ such as $V_{ccint_{1}}$,$V_{ccint_{2}}$, $V_{ccint_{3}}$ and
$V_{ccint_{4}}$ for FPGA partition-1, partition-2, partition-3 and
partition-4, respectively. The power distribution unit distributes
$V_{ccint_{i}}$ such as $V_{ccint_{1}}$, $V_{ccint_{2}}$, $V_{ccint_{3}}$ and
$V_{ccint_{4}}$ to FPGA partition-1, partition-2, partition-3 and partition-4
respectively. Thereafter, the TPU circuit can be on and the runtime scheme
becomes functional. In Fig. 8, 4 FPGA partitions, partition-1, partition-2,
partition-3 and partition-4 have 4 flags form $Razor$ flipflops,
$timing\\_fail-part-1$, $timing\\_fail-part-2$, $timing\\_fail-part-3$ and
$timing\\_fail-part-4$ respectively to detect the timing failure of the
available partition of the FPGA. Each $timing\\_fail-part-i$ flag is ANDed
value of all error detection flag all MACs placed in the $i^{th}$ partition.
As shown in algorithm 2, if $i^{th}$ timing failure flag from $i^{th}$ FPGA
partitions becomes high, the power distribution network will step up
$V_{ccint_{i}}$ of that partition by $V_{s}$ else $V_{ccint_{i}}$ will be step
down by $V_{s}$.

Figure 9: Flow Diagram of Proposed Frame Work

##  IV Clustering Algorithms

We have investigated 4 clustering algorithms to group the MACs having similar
minimum slacks. Algorithms can be chosen based on the design requirements: if
we want to set a pre-defined number of clusters, or set hyperparameters to
automatically determine the number of clusters. Different algorithms work well
for different data distributions. Depending on our design requirements, we
choose among the following four algorithms:

###  IV-A Hierarchical

The hierarchical clustering [12] algorithm considers each data point as a
single cluster and measures distance between two clusters based on a chosen
distance measure (in this case, Euclidean distance). The two clusters that are
closest to each other are merged. The process is continued until all clusters
have been merged into a single cluster (root of the dendrogram). As shown in
fig. 10, the dendrogram is a tree-like structure used for visualizing the
hierarchy of clusters. The number of clusters can be decided from the
dendrogram. The hierarchical algorithm is computationally expensive for large
datasets, having a time complexity of $O(n^{3})$ where $n$ is the number of
data-points. As is evident from the dendrogram, the length of the branch
joining the last two clusters is the highest, indicating they are the most
dissimilar, followed by the third and fourth clusters. The result of
classifying the slack values into 2, 3 and 4 clusters is illustrated in the
fig. 11. Different clusters in Figs. 11, 12, 13 and 14 are indicated using
different colours.

Figure 10: Dendrogram

Figure 11: Hierarchical Cluster of Slack of Systolic Array $16\times 16$: a)
#clusters=2, b) #clusters=3, c) #clusters=4

###  IV-B K-Means Clustering

K-Means Clustering can cluster data into a predefined number of groups ($k$).
At the beginning, $k$ cluster centers are randomly initialized [13]. The
algorithm computes the distance between each data-point and the cluster-
centers and assigns data-points to the cluster whose center is closest to it.
The cluster centers are then recomputed as the mean of the data-points
belonging to that cluster. The process is repeated for a predefined number of
steps or until cluster centers do not change significantly. The K-Means
Clustering is simple, fast, and its time complexity is $O(n)$. Fig. 12
illustrates the results of applying K-Means clustering algorithm on the
minimum slack values of a $16\times 16$ Systolic Array (256 MACs) for 3, 4 and
5 clusters.

Figure 12: K Means Cluster of Slack of Systolic Array$16\times 16$: a)
#clusters=3, b) #clusters=4, c) #clusters=5

###  IV-C Mean-Shift Clustering

Mean Shift Clustering [14] is based on the idea of Kernel Density Estimation
(KDE). KDE assumes that the data points are generated from an underlying
distribution and tries to estimate the distribution by assigning a kernel to
each data point. The most commonly used kernel is the Gaussian or RBF kernel.
The mean-shift algorithm is designed in a way that the points iteratively
climb the KDE surface and are shifted to the nearest KDE peaks. It starts with
a randomly selected point as the center of the RBF kernel. Thereafter, it
proceeds by moving the kernel towards regions of higher density by shifting
the center of the kernel to the mean of the points within the window (hence
the algorithm is termed mean-shift). This is continued until shifting the
kernel no longer includes more points. This algorithm does not need the number
of clusters to be specified beforehand, but it is computationally expensive
compared to K-Means (time complexity of $O(n*log(n))$ in lower dimension for
sklearn implementation). The selection of the window size/radius $r$ can be
non-trivial and plays a key-role in the success of the algorithm. Setting the
radius as 0.4 for the slack values of a $16\times 16$ Systolic array, yields 4
clusters as observed in the Fig. 13.
Figure 13: Mean-Shift Clustering of Slack of Systolic Array$16\times 16$
Figure 14: DB Scan Clustering of Slack of Systolic Array $16\times 16$

###  IV-D DBSCAN

The DBSCAN algorithm has two important hyperparameters, based on which it
determines the number of clusters [15], epsilon: The maximum distance between
two samples for one to be considered as in the neighborhood of the other and
minpoints: The number of samples in a neighborhood for a point to be
considered as a core point. At each step, a data-point that has not been
visited before is taken. If there are more data-points than $minpoints$ within
its $epsilon$ radius, all the data-points are marked as belonging to a
cluster, otherwise the first point is marked as noise. For all points in the
newly-formed cluster, points within their ‘epsilon’ neighborhood are checked
and labeled as either belonging to a cluster or noise. The process is
continued until all data-points have been labeled. The greatest advantage of
DBSCAN is that it can identify outliers as noise, unlike other algorithms
which throw all points into a cluster even if one data point is significantly
different from the rest. The time complexity of this algorithm is $O(n)$ for
reasonable $epsilon$. This algorithm is not effective for clusters with
varying density since $epsilon$ and $minpoints$ are different for different
clusters.

##  V Implementation and Result

As mentioned in Sec. II, the two proposed tool flow has 2 environments. The
clustering algorithms for both $Vivado$ and VTR are implemented in Python
using the Scikit-learn library. The $synthesis$, $implementation$ and
$bit~{}file~{}generation$ of $Vivado$ flow is done by the board support
package of $Artix-7$ FPGA. The $synthesis$ and $implementation$ of $VTR$ flow
is done by the board support package of 22nm, 45nm and 130nm academic FPGAs.
As shown in and Fig. 9, the cluster algorithm generates $P$ no. of partition,
and the dimension of each partition $(n\times m)$. The static scheme generates
biasing voltages:

$P\times(n\times
m)\\{V_{ccint_{1}},V_{ccint_{2}},...,V_{ccint_{i}},...,V_{ccint_{p}}\\}$ (1)

for $P$ no.of partitions. The runtime scheme calibrate the biasing voltage
according to the timing failure detected by Razor placed in every MAC. Runtime
scheme gives the final set of biasing voltages:

$\begin{split}V_{ccint_{1}}+C1.V_{s},V_{ccint_{2}}+C2.V_{s},...,V_{ccint_{i}}+Ci.V_{s},...,\\\
V_{ccint_{p}}+Cp.V_{s}\end{split}$ (2)

Here $C_{1}$, $C_{2}$, .. $C_{p}$ are integers starts from $0$ to any positive
value

TABLE II: Comparison of Dynamic Power (mw) for $Vivado$ and $VTR$ flow

Schemes under 25 ∘ Ambient Temperature &100MHz Clock |  Dimension of Systolic Array |  Partition No. |  $V_{ccint_{i}}$ volt |  Vivado 28nm Artix-7 |  VTR 22nm |  VTR 45nm |  VTR 130nm
---|---|---|---|---|---|---|---
Without Voltage Scaling |  $16\times 16$ |  NA |  1.00 |  408 | 269 | 387 | 1543
| $8\times 8$ |  partition-1 | 0.96 |  |  |  |
Voltage | $8\times 8$ |  partition-2 | 0.97 | 382 | 263 | 380 | 1531
Scaled | $8\times 8$ |  partition-3 | 0.98 |  |  |  |
| $8\times 8$ |  partition-4 | 0.99 |  |  |  |
% of Reduction | 6.37 | 1.86 | 1.8 | 0.7
Without Voltage Scaling |  $32\times 32$ |  NA |  1.00 |  1538 | 1072 | 1549 | 6172
| $16\times 16$ |  partition-1 | 0.96 |  |  |  |
Voltage | $16\times 16$ |  partition-2 | 0.97 | 1404 | 1051 | 1520 | 6125
Scaled | $16\times 16$ |  partition-3 | 0.98 |  |  |  |
| $16\times 16$ |  partition-4 | 0.99 |  |  |  |
% of Reduction | 6.76 | 1.95 | 1.87 | 0.76
Without Voltage Scaling |  $64\times 64$ |  NA |  1.00 |  5920 | 4284 | 6200 | 24693
| $32\times 32$ |  partition-1 | 0.96 |  |  |  |
Voltage | $32\times 32$ |  partition-2 | 0.97 | 5534 | 4205 | 6090 | 24503
Scaled | $32\times 32$ |  partition-3 | 0.98 |  |  |  |
| $32\times 32$ |  partition-4 | 0.99 |  |  |  |
% of Reduction | 6.52 | 1.84 | 1.77 | 0.77
Without Voltage Scaling |  $64\times 64$ |  NA |  0.9 |  not supported | 3965 | 5798 | 23961
| $32\times 32$ |  partition-1 | 0.7 |  |  |  |
Voltage | $32\times 32$ |  partition-2 | 0.8 | not | 3818 | 5656 | 23631
Scaled | $32\times 32$ |  partition-3 | 0.9 | supported |  |  |
| $32\times 32$ |  partition-4 | 1.00 |  |  |  |
% of Reduction | - | 3.7 | 2.4 | 1.37

###  V-A Implementational Challenges

The proposed design could not be implemented as none of the present-day FPGA
devices support variable voltage scaling in the different logic partitions.
The implementation issues of power distribution unit with multiple $V_{ccint}$
in different partitions are beyond the scope of our paper. However, we
consider, the implementation of voltage scaling technology in ASIC [4]
establishes the feasibility of implementation of voltage scaling technology in
FPGA.

Figure 15: Comparison of Dynamic Power (mw) of Various Variance of 64 $\times$
64 Systolic Array on 22nm and 45nm  Figure 16: Comparison of Dynamic Power
(mw) of Various Variance of 64 $\times$ 64 Systolic Array on 130nm

###  V-B Our Validation Strategy

To validate the claim of the proposal, we have simulated the proposed scheme
using VTR and Vivado flow. We have designed a 3 Systolic arrays with the
dimension of $16\times 16$, $32\times 32$, and $64\times 64$. Let us take the
example of $16\times 16$ systolic array where $16\times 16=256$ MACs are
placed in the FPGA. As shown in Fig. 12(b), the K means clustering algorithm
mentioned in Sec. IV divides $16\times 16$ systolic array into 4 partitions:
$partition-1$, $partition-2$, $partition-3$ and $partition-4$. Though the size
of the partitions in Fig. 12 are not same, for sake of simplicity of
implementation we have assumed the same partition size $(8\times 8)$. As the
current $Vivado$ tool does not allow simulating the design in critical voltage
region, our $16\times 16$ systolic array is tested in the guardband region.
Due to the unavailability of multiple $V_{ccint}$ support in single a FPGA
dice, our design has been implemented in one partition at a time. Therefore,
the power measurement of 4 partitions is also done separately where each
partition is considered as an individual circuit. Though $VTR$ allows design
in critical regions, for sake of better comparative study we have also used
the same voltage ranges used in $Vivado$.

###  V-C Results

The guardband region for $Artix-7$ FPGA is 0.95 volt to 1.00 volt. For this
example $n=4$, $V_{min}=V_{nom}=1.00~{}volt$, $V_{crash}=V_{min}=0.95~{}volt$,
therefor $V_{b}=0.05~{}volt$. Algorithm 1 calculates the $V_{ccint_{i}}$ of
the 4 FPGA partitions of this design which are : $V_{ccint_{1}}=0.956$ for
$partition-1$, $V_{ccint_{2}}=0.968$ for $partition-2$, $V_{ccint_{3}}=0.985$
for $partition-3$ and $V_{ccint_{4}}=0.993$ for $partition-4$. We have
observed that when the partial sums are moved to the bottom rows of systolic
array, the timing error increases significantly [4]. In this example the MACs
of bottom rows have less minimum slacks, which should be placed in
$partition-3$ and $partition-4$ where $V_{ccint_{i}}$ is more compared to the
existing $V_{ccint_{i}}$. The MACs of upper rows should have more minimum
slacks which should be placed in $partition-1$ and $partition-2$ where
$V_{ccint_{i}}$ is less compared to the existing $V_{ccint_{i}}$. As shown in
Fig 8, the clustering algorithm divides the $16\times 16$ systolic array into
four $8\times 8$ systolic array partitions and each partition has $8\times
8=64$ MACs. The top-left partition-1 consists of a $8\times 8$ systolic array
which has $V_{ccint_{1}}=0.956\approx 0.96$. Similarly, top-right partition-2
has $V_{ccint_{2}}=0.968\approx 0.97$, bottom-left partition-3 has
$V_{ccint_{3}}=0.985\approx 0.98$ and bottom right partition-4 has
$V_{ccint_{4}}=0.993\approx 0.99$. Table II shows the dynamic power
consumption of $16\times 16$, $32\times 32$ and $64\times 64$ systolic arrays
with 4 partitions. Table II shows the adoption of voltage scaling technology
reduces 6.97% to 6.76% dynamic power consumption for $Vivado$ commercial FPGA
and 0.7% to 1.95% VTR academic FPGA. Due to the limited range of voltage in
Vivado, these are the lower bounds on improvements. As we can drop biasing
voltage closer to NTC for next generation FPGAs, we expect the improvement to
improve substantially. In the $4^{th}$ instant, VTR allows more lower biasing
voltages which reduces the dynamic power consumption 3.7%, 2.4% and 1.37% for
22nm, 45nm and 130nm respectively.

In Fig. 15 and Fig. 16, we have shown dynamic power consumption of different
variance of 64 $\times$ 64 systolic array on 22nm, 45nm and 130nm academic
FPGAs using VTR flow. Fig. 15 and Fig. 16 show that variation of 3 parameters
such as number of partition $P$, biasing voltage $V_{ccint_{i}}$ of each
partition and dimension of each FPGA partition $n\times m$ change the dynamic
power consumption of 64 $\times$ 64 systolic array by 18%, 21% and 39% for
22nm, 45nm and 130nm academic FPGAs respectively. Here the number of cluster
or partition $P$ and dimension of each partition $n\times m$ will be
calculated by cluster algorithms. The biasing voltage $V_{ccint_{i}}$ of each
FPGA partition are roughly calculated by static scheme and further calibration
of accurate $V_{ccint_{i}}$ is done by the runtime scheme. Such significant
effect of $P$, $n\times m$ and $V_{ccint_{i}}$ on dynamic power consumption
shows cluster algorithm, static and runtime schemes are very crucial steps of
proposed frame. As an example, in Fig. 15 and 16, a name of one variance of
systolic array is $4\times(32\times 32)\\{0.8,1.0,1.2,1.3\\}$ (the rightmost
bar at Fig. 16), here $P$=4, $n\times m$ =$32\times 32$ and biasing voltages
are of 4 partitions are 0.8 volt, 1.0 volt, 1.2 volt and 1.3 volt. In Fig. 15
and 16, the $V_{ccint_{i}}$ for 130 nm varies from threshold voltage 0.7 volt
to 1.3 volt whereas for 22 nm and 45 nm $V_{ccint_{i}}$ varies from 0.5 volt
to 1.2 volt. Though the threshold voltage of 45nm is 0.5 volt and for 22nm it
is 0.45 volt, for comparative purposes, in both cases we have measured it form
0.5 volt. It is known that the dynamic power reduces by the square of the
supply voltage $V_{ccint_{i}}$. Also, the $2\times(32\times 64)\\{0.5,0.6\\}$
variant of systolic array implemented in 22nm and 45nm technology have maximum
number of MACs which are running with minimum $V_{ccint_{i}}$ as compared to
other reported variants as shown in Fig. 15. Thus, the aforementioned variant
consumes minimum dynamic power as compared to other variants reported in Fig.
15 . Going by the same reasoning, the $2\times(32\times 64)\\{0.7,0.8\\}$
variant in 135 technology consumes minimum dynamic power when compared to the
other variants as reported in Fig 16. The minimum voltage step of the power
supply [11] is considered as 0.1 volt. We observe that the timing reports of
16 $\times$ 16, 32 $\times$ 32 and 64 $\times$ 64 systolic arrays before
partitioning and after partitioning, which shows very insignificant effects on
delay in wires and placement and routing difficulties. Hence, the re-
clustering process is not required for the aforementioned systolic arrays.

##  VI Conclusion

FPGA is becoming popular in DNN based configurable cloud because of its
efficiency and manageability but immoderate power consumption is a growing
concern for present FPGA technology. A lot of effort has been made to reduce
the power consumption by using multiple biasing voltages in FPGA. This paper
proposes a systolic array where the MACs are placed in different partitions of
FPGA based on the minimum slacks of different MACs. Each partition of the FPGA
uses different biasing voltage $V_{ccint}$. The proposed $runtime$ and
$static$ schemes can tune appropriate $V_{ccint}$ with the group MACs having
similar minimum slacks placed in the partitions. The experimental results show
that the voltage scaled systolic array can reduce power consumption. In future
we will address two points such as (i) Improvement of $V_{ccint}$ calibration
by grouping input sequences with similar delay characteristics to predict
future timing failures. (ii) Study the tradeoff between the DNN accuracy
estimated in terms of timing failures with the no. of partitions and that
between no. of partitions and dynamic power. The same partition based voltage
scaling can be used for other high performance hardware accelerator to reduce
the power consumption.

## References

  * [1] A. M. Caulfield and et. al.  A cloud-scale acceleration architecture.  In 2016 49th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO), pages 1–13, 2016.
  * [2] A. Putnam and et al.  A reconfigurable fabric for accelerating large-scale datacenter services.  IEEE Micro, 35(3):10–22, 2015.
  * [3] B. Salami, E. B. Onural, I. E. Yuksel, F. Koc, O. Ergin, A. Cristal Kestelman, O. Unsal, H. Sarbazi-Azad, and O. Mutlu.  An experimental study of reduced-voltage operation in modern fpgas for neural network acceleration.  In 2020 50th Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), pages 138–149, 2020.
  * [4] P. Pandey, P. Basu, K. Chakraborty, and S. Roy.  Greentpu: Improving timing error resilience of a near-threshold tensor processing unit.  In 2019 56th ACM/IEEE Design Automation Conference (DAC), pages 1–6, 2019.
  * [5] D. Ernst, Nam Sung Kim, S. Das, S. Pant, R. Rao, Toan Pham, C. Ziesler, D. Blaauw, T. Austin, K. Flautner, and T. Mudge.  Razor: a low-power pipeline based on circuit-level timing speculation.  In Proceedings. 36th Annual IEEE/ACM International Symposium on Microarchitecture, 2003. MICRO-36., pages 7–18, 2003.
  * [6] Kevin E. Murray, Oleg Petelin, Sheng Zhong, Jia Min Wang, Mohamed Eldafrawy, Jean-Philippe Legault, Eugene Sha, Aaron G. Graham, Jean Wu, Matthew J. P. Walker, Hanqing Zeng, Panagiotis Patros, Jason Luu, Kenneth B. Kent, and Vaughn Betz.  Vtr 8: High-performance cad and customizable fpga architecture modelling.  ACM Trans. Reconfigurable Technol. Syst., 13(2), May 2020.
  * [7] P. Jamieson, K. B. Kent, F. Gharibian, and L. Shannon.  Odin ii - an open-source verilog hdl synthesis tool for cad research.  In 2010 18th IEEE Annual International Symposium on Field-Programmable Custom Computing Machines, pages 149–156, 2010.
  * [8] Berkley Logic Synthesis and Verification Group. 2018.  Abc: A system for sequential synthesis and verification.  In Revison, 2018.
  * [9] Jason Luu, Ian Kuon, Peter Jamieson, Ted Campbell, Andy Ye, Wei Mark Fang, Kenneth Kent, and Jonathan Rose.  Vpr 5.0: Fpga cad and architecture exploration tools with single-driver routing, heterogeneity and process scaling.  ACM Trans. Reconfigurable Technol. Syst., 4(4), December 2011.
  * [10] R. Mukherjee and Seda Ogrenci Memik.  Realizing low power fpgas : A design partitioning algorithm for voltage scaling and a comparative evaluation of voltage scaling techniques for fpgas.  2005\.
  * [11] T. N. Miller, X. Pan, R. Thomas, N. Sedaghati, and R. Teodorescu.  Booster: Reactive core acceleration for mitigating the effects of process variation and application imbalance in low-voltage chips.  In IEEE International Symposium on High-Performance Comp Architecture, pages 1–12, 2012.
  * [12] Stanford University.  Hierarchical agglomerative clustering.  2008 Cambridge University Press, 2008.
  * [13] David Arthur and Sergei Vassilvitskii.  K-means++: the advantages of careful seeding.  In In Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms, 2007.
  * [14] D. Comaniciu and P. Meer.  Mean shift: a robust approach toward feature space analysis.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5):603–619, 2002.
  * [15] Martin Ester, Hans-Peter Kriegel, Jörg Sander, and Xiaowei Xu.  A density-based algorithm for discovering clusters in large spatial databases with noise.  In Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, KDD’96, page 226–231. AAAI Press, 1996\.


# Robust Tumor Segmentation with Hyperspectral Imaging
and Graph Neural Networks

Mayar Lotfy1,2 Anna Alperovich3 Tommaso Giannantonio3 Björn Barz3 Xiaohan
Zhang1 Felix Holm2 Nassir Navab2 Felix Boehm4 Carolin Schwamborn4, Thomas K.
Hoffmann4, Patrick J. Schuler4
1Carl Zeiss Meditec AG, Oberkochen, Germany
2CAMP, Technical University of Munich, Garching, Germany
3Carl Zeiss AG, Corporate Research & Technology, Oberkochen, Germany
4Dept. of Otorhinolaryngology, University Hospital Ulm, Ulm, Germany

###### Abstract

Segmenting the boundary between tumor and healthy tissue during surgical
cancer resection poses a significant challenge. In recent years, Hyperspectral
Imaging (HSI) combined with Machine Learning (ML) has emerged as a promising
solution. However, due to the extensive information contained within the
spectral domain, most ML approaches primarily classify individual HSI
(super-)pixels, or tiles, without taking into account their spatial context.
In this paper, we propose an improved methodology that leverages the spatial
context of tiles for more robust and smoother segmentation. To address the
irregular shapes of tiles, we utilize Graph Neural Networks (GNNs) to
propagate context information across neighboring regions. The features for
each tile within the graph are extracted using a Convolutional Neural Network
(CNN), which is trained simultaneously with the subsequent GNN. Moreover, we
incorporate local image quality metrics into the loss function to enhance the
training procedure’s robustness against low-quality regions in the training
images. We demonstrate the superiority of our proposed method using a clinical
ex vivo dataset consisting of 51 HSI images from 30 patients. Despite the
limited dataset, the GNN-based model significantly outperforms context-
agnostic approaches, accurately distinguishing between healthy and tumor
tissues, even in images from previously unseen patients. Furthermore, we show
that our carefully designed loss function, accounting for local image quality,
results in additional improvements. Our findings demonstrate that context-
aware GNN algorithms can robustly find tumor demarcations on HSI images,
ultimately contributing to better surgery success and patient outcome.

Figure 1: Overview of the segmentation methods we consider. The top section
(a, blue arrows) describes the baseline CNN method, while the bottom section
(b, orange arrows) represents the CNN+GNN approach. Beginning from the left,
we represent the full HSI data cube for one image (the grayscale shows
luminosity at a single representative wavelength) with exemplary SLIC tiling
overlaid in yellow. The numbers on the $x,y,\lambda$ axes represent the
dimensions. Along the (a) path we show a representative cropped and zero-
padded tile, the architecture of the CNN classifier model (each block
represents a convolutional layer with feature dimensions at the bottom,
spatial dimensions at the left, and kernel dimensions near the blue
rectangle), its possible predictions, and their aggregation to image
segmentation. Along the (b) path, we show the graph representation of the full
HSI image, incorporating as node features the CNN embeddings, the GNN model
(each block represents a GAT layer with numbers of input and output channels
and attention heads), and the resulting node predictions, which can be
directly interpreted as image segmentation result.

##  1 Introduction

Head and Neck (H&N) cancer is the sixth most common cancer worldwide [46]. It
is located in the upper aerodigestive tract, which comprises the oral and
nasal cavities, pharynx, and larynx. The majority of these cancers are
squamous cell carcinoma (SCC) and they are often lethal, with a total number
of 529,000 new cases yearly resulting in 350,000 cancer-related deaths [15,
16]. Currently, the primary treatment remains surgery, whose outcome is
compromised by inadequate cancer margin detection in 10% to 20% of the cases.
If residues of cancerous tissue remain, recurrence rates after surgery can be
as high as 55% [39].

In the typical surgical workflow, doctors use pre-operative imaging, such as
CT or MRI, to plan tumor resection. In addition to that, during the surgical
procedure, surgeons rely strongly on prior experience and visual cues, as well
as intra-operative pathologist consultations (IPCs) that provide guidance
towards the total resection of the cancer tissue. However, IPC is time
consuming and may not fully reflect the extent and the delineation of the
lesion due to limitations in tissue sampling and preparation. Thus, tumor
segmentation solutions based on optical imaging combined with machine learning
recently emerged as an attractive addition to the classical pathological
testing [65]. This approach holds the promise of non-invasive and real-time
tumor segmentation during the surgery.

Hyperspectral imaging (HSI) [7] is an optical modality that combines the
benefits of imaging and spectroscopy. It covers a continuous part of the
visible (and often IR) spectrum and images the whole field of view at hundreds
of narrow wavelength bands, thus recording significantly more information than
regular broadband imaging, such as RGB. HSI is an emerging technique in the
biomedical laboratory and in the operating room [37, 7], as it is non-
invasive, non-ionizing, and label-free, i.e., does not require a contrast
agent. The information provided by the spectral properties of biological
tissues can be used to distinguish them [30].

HSI has been used on biopsies, histopathological and fluorometric analysis,
and disease biology [37, 14, 43]. In medical computer vision research, it was
used for tumor/healthy tissue differentiation of multiple cancer types [24,
11, 21].

These recent studies relied on the assumption that the spectral domain should
contain most of the necessary information for tumor differentiation. Thus,
hyperspectral images are decomposed into single pixels or small spatial tiles
that are analyzed independently. This has two significant drawbacks: firstly,
it destroys any spatial context information, and secondly, the tiles typically
have varying data quality, which may bias the results.

#### Contributions

In this work, we present an ML method for overcoming the two aforementioned
limitations of the current state of the art. By introducing additional context
and application information to the ML model, we improve the prediction quality
for the task of tumor/healthy tissue segmentation on unseen samples.

Firstly, we introduce spatial context information to the model by utilizing
Graph Neural Networks (GNNs). We connect hyperspectral tiles based on their
spatial proximity into a graph and feed it to a GNN, whereby image
segmentation is seen as a graph node classification task. Our novel approach
extracts tile-specific features with a CNN encoder and uses them as the node
features for the GNN (see Fig. 1). An end-to-end training strategy ensures
that task-specific information is extracted simultaneously on various scales,
optimally combining spatial and spectral information.

Secondly, we investigate whether the local quality of HSI data influences the
segmentation results by changing the training set from only spatially and
spectrally homogeneous tiles to all available tiles. We confirm that tiles of
a lower quality harm network performance. We then derive a custom loss
function that weights the contribution from each spatial tile based on custom
quality metrics. Thus, information of the lower-quality data is still used,
but it does not contribute as much as good-quality data to the final model.

Our quantitative evaluations show that all GNN models outperform CNN-based
counterparts, improving accuracy for all classes. Qualitative results
illustrate that predictions from GNN-based models are smoother and do not
suffer from isolated classification errors. Furthermore, our novel quality-
aware loss function exhibits superior performance compared to the baseline,
especially on the tumor class.

The rest of this paper is structured as follows: after discussing related work
in Sec. 2, we present our methods in Sec. 3, describe the results in Sec. 4,
and conclude in Sec. 5.

##  2 Related work

Here we first summarize the state of the art of medical HSI data segmentation,
then we review GNN applications to medical non-HSI data and to non-medical HSI
data. To the best of our knowledge, no previous work used GNNs on medical HSI
data.

###  2.1 Semantic segmentation of medical HSI data

Different methods exist to maximise information extraction for HSI-based
tissue segmentation [45]. In the simplest cases, a physical inverse model can
be used [34]. More often however, a forward model is lacking, and data-driven
ML models are required [20], either classical [40, 2] or based on deep
learning [35, 32, 8].

Different strategies are possible for the combinations of spatial and spectral
information. The simplest single-pixel classification relies entirely on the
spectral domain [2]; alternatively, spectral and spatial domains can be
combined, e.g., by engineering spectral-spatial features in classical ML [28].
Most deep learning approaches have been based on multi-channel 2D CNNs,
applied independently to macropixels (patches) to improve signal-to-noise and
utilize at least part of the spatial information. This was the strategy
adopted by the HELICoiD project [12], which developed an intraoperative
demonstrator for in vivo brain tumor segmentation for surgical support. Their
final dataset consisted of 36 images from 22 gliobastoma patients [13]. A
combination of classical and deep learning methods (2D CNNs) based on single
pixels and patches achieved accuracies $>98\%$ [11, 25, 38, 22]; patch-wise
fully 3D CNNs were also used [25, 38]. The same patch-based models were used
for in vitro segmentation [44] and other medical applications, including H&N
cancer (SCC) segmentation [26, 23, 25]. A similar pixel-wise model was applied
by Urbanos et al. [59] on an independent HSI brain cancer dataset, and patch-
wise using SLIC [1] aggregation by Giannantonio et al.[21] on low-grade glioma
data.

While the typical HSI data sizes and data scarcity can make full-image models
unwieldy, full-image segmentation approaches were also proposed. Trajanovski
et. al [57, 58] segmented an HSI dataset of tongue cancer from 14 patients
using U-Net [49] based on patches of size $256\times 256$. Cervantes et al.
[5] segmented HSI data from hepatic (18 images) and tyrhoid (21 images)
surgeries with U-Net on patches of size $8\times 8$. Garifullin et al. [19]
segmented 55 retinal HSI images using SegNet [3] and Dense-FCN [31]. Wang et
al. [63] defined the 3D Hyper-Net encoder-decoder model to segment HSI
melanoma images. Yun et al. [70] presented the SpecTr Spectral Transformer
model, by introducing a transformer into the encoding part of a U-Net, on a
dataset of 514 cholangiocarcinoma HSI images. Seidlitz et al. [52] compared
pixel- vs. patch-based vs. full-image organ segmentation on 506 HSI animal
images.

###  2.2 GNN segmentation of medical (non-HSI) data

Graph Neural Networks (GNNs) [66, 6] are deep learning models that operate on
data in graph form [4]; they have been successfully applied to diverse tasks
in recent years, since the introduction of the Graph Convolutional Network
model [33] and its later extensions such as GraphSAGE [27] and the Graph
Attention Network (GAT) [60].

Any image can be seen as a graph, with the pixels as the nodes and the edges
drawn, e.g., to the neighbours. This idea is attractive in case of large
medical data, such as HSI but also volumetric data such as MRI or CT scans,
which can be first aggregated into super-pixels (-voxels), and then processed
as a graph with GNNs. The benefits are twofold: circumventing the cost of
processing the raw data cube with CNNs, and incorporating domain knowledge in
the definition of nodes and edges.

This approach has been applied to several medical image segmentation task. Yan
et al. [67] used GNNs (ChebNets) to segment brain tissue from MRI data. Juarez
et al. [18] combined U-Nets with GNNs to segment airways from CT scans.
Saueressig et al. [50, 51] performed brain tumor segmentation from MRI images,
obtaining best results with GraphSAGE. Duan et al. [10] segmented eye fundus
images (RGB) using a combination of CNN and GNN features extracted with a GAT
model. Patel et al. [47] also segmented volumetric MRI data using GAT. Gaggion
et al. [17] defined the CNN+GNN encoder-decoder model HybridGNet for landmark-
based segmentation of chest X-ray data. Shao et al. [56] tracked liver tumor
in X-ray data combining GNNs with a biomechanical model.

###  2.3 GNN segmentation of non-medical HSI data

Outside the medical domain, multiple GNN models have been developed for HSI
image segmentation, mostly taylored to remote sensing and land coverage
classification.

Qin et al. [48] defined the spectral-spatial GCN that combines spectral and
spatial similarity in the edge definition. Shahraki et al. [55] combined a 1D
CNN feature extractor with GCNs. Sha et al. [54] employed edge convolutions in
semi-supervised mode, while Hong et al. [29] presented the miniGCN model that
subdivides nodes in mini-batches. Wan et al. [61] introduced the multi-scale
dynamic GCN (MDGCN) model that considers the connection to multi-hop
neighbors; later from this model the context-aware dynamic GCN was derived
[62], where the adjacency matrix is dynamically refined. Several further
variants were introduced [36], including defining a graph over the whole HSI
image [41] and using GATs [53, 64] and GraphSAGE [68]. Zhang et al. [71] later
introduced the Global Random Graph Convolution Network (GR-GCN), which reduces
graph size by random sampling. Ding et al. [9] proposed the multi-scale
receptive fields graph attention network (MRGAT), which combines CNNs with
GATs on multiple scales. Yao et al. [69] presented a deep hybrid multi-graph
neural network (DHMG) that employs a spectral filter to extract features and
an auto-regressive average for noise suppression.

##  3 Methodology

###  3.1 Data acquisition

Figure 2: Laboratory setup at the ENT department of University Hospital Ulm:
an IMEC HSI snapscan camera connected to a ZEISS OPMI® PENTERO® microscope for
ex vivo HSI and RGB data acquisition.

We follow the configuration decribed in Ref. [21] for HSI data acquisition. In
brief, ex vivo HSI image acquisition was performed at the ENT department of
University Hospital Ulm, employing an IMEC Snapscan VNIR hyperspectral camera
in conjunction with a ZEISS OPMI® PENTERO® 900 surgical microscope through its
optical port.

The HSI setup shown in Fig. 2 is prepared before the start of tumor surgery.
Following tumor resection, the ex vivo tumor specimen undergoes a saline
solution rinse to eliminate residual blood from its surface. Next, the sample
is promptly taken to the laboratory-based HSI setup for data acquisition. For
each specimen, a minimum of two HSI images are typically captured: an initial
panoramic image covering the entire sample area to facilitate subsequent
annotation of tissue categories, and an image with a reduced working distance
and higher magnification to scrutinize the tumor center and margin in greater
detail. A consistent lighting scheme and region of interest are maintained by
capturing an RGB photograph prior to each HSI image acquisition using the
PENTERO microscope. The selection of microscope working distance and
magnification values is determined based on the dimensions of the respective
specimens. A constant light intensity of 50% is used on the PENTERO microscope
for all HSI measurements.

Post-operation, surgeons annotate diverse tissue classes within RGB images
based on their prior knowledge from pre- and intra-surgical diagnostics. The
RGB images and annotations are subsequently registered with their HSI
counterparts.

###  3.2 Data characterization and preparation

(a) RGB image

(b) Annotated image

Figure 3: RGB image of one exemplary resected tissue sample (a) and its
corresponding labels (b) annotated by medical experts, depicting in red,
yellow, blue, green, and magenta the classes tumor, tumor margin, healthy,
muscle and background respectively.

Our dataset consists of HSI images of otolaryngological tumor of SCC type. It
comprises a total of 51 annotated images from 30 patients with resolution
$1600\times 1600$ pixels. We remove two images, as based on medical assessment
they were found to display a different tumor type, thus leaving us with 49
images for analysis. Every HSI image has $N_{\lambda}=109$ wavelength channels
in the visible spectrum between approximately 468 nm – 790 nm. Four classes
were annotated by medical experts namely, tumor, tumor margin, healthy and
muscle. We additionally annotated a background class using the same tool for
extending model robustness. For all experiments, we only consider the three
classes tumor, healthy, and background, as the muscle class is poorly sampled
and tumor margin has large annotation uncertainty. We show in Fig. 3 one
example RGB image and its annotations.

In order to boost signal-to-noise and support the segmentation task, we
aggregate pixels into tiles (macro-pixels) by using a modified version of the
SLIC algorithm [1], which aggregates pixels to minimize a combination of
spatial and color distances between them. We set the color distance to be the
Spectral Angle Mapper (SAM) distance, which for two pixels $i,j$ of spectra
$\mathbf{s}_{i},\mathbf{s}_{j}$ is given by

$\mathrm{SAM}(\mathbf{s}_{i},\mathbf{s}_{j})=\cos^{-1}\left(\frac{\mathbf{s}_{i}\cdot\mathbf{s}_{j}}{\lVert\mathbf{s}_{i}\rVert\,\lVert\mathbf{s}_{j}\rVert}\right)\,.$
(1)

Alternatively, the $L_{2}$ distance can be used, as given by

$L_{2}(\mathbf{s}_{i},\mathbf{s}_{j})=\sqrt{\sum_{l=1}^{N_{\lambda}}(s_{il}-s_{jl})^{2}}\,.$
(2)

We set the number SLIC tiles to contain $\sim 200$ pixels each.

|
---|---
|

(a) RGB image

(b) Average intensity

(c) Average $L_{2}$

(d) Average SAM

Figure 4: Tile quality metrics of the representative tissue sample shown in
panel (a). The remaining panels display the three tile quality metrics we
consider: average tile intensity $I$ (panel b), and the average intra-tile
spectral distances from the mean $L_{2},\mathrm{SAM}$ in panels (c, d)
respectively.

We consider as possible proxies of each tile’s $T_{i}$ quality its average
pixel intensity $I_{i}$ (to isolate under-exposed and saturated regions) and
the quantities $L_{2,i}$, $\mathrm{SAM}_{i}$ that are obtained by taking the
average distances between the spectrum at each pixel $\mathbf{s}_{j}$ $\forall
j\in T_{i}$ and the tile’s average spectrum $\overline{\mathbf{s}}_{i}$, i.e.,
$L_{2,i}\equiv\frac{1}{|T_{i}|}\sum_{j\in
T_{i}}L_{2}(\mathbf{s}_{j},\overline{\mathbf{s}}_{i})$ (large spectral
distances indicate low uniformity). (We drop the $i$ indices henceforth for
simplicity.) As depicted in Fig. 4, tile quality across an image displays
significant variations. Therefore, to build a dataset with only high-quality
tiles, we follow the approach by Ref. [21] and discard the tiles that have
poor quality by applying threshold cuts on the three mentioned metrics: low
spectral uniformity (we discard the upper 75% percentile on both average SAM
and $L_{2}$), extremely high or low mean intensity (we discard the lowest and
highest 10% percentiles in $I$), as well as tiles with non-uniform class
labels. This results in a dataset consisting of high-quality tiles only.

The total amount of tiles in our dataset is 36 693, of which 26 190 are high-
quality tiles. We split the data into training+validation and testing samples
patient-wise, retaining all images from six patients for testing, i.e., we use
39 images for training and validation and 10 images for testing. We divide
accordingly both full and high-quality dataset into three sets containing 65%,
16.5%, and 18.5% of the tiles for the train, val, test respectively.

###  3.3 Tissue segmentation methods

We consider the two main segmentation methods shown in Fig. 1, CNN-based and
CNN+GNN-based, both of which can be extended with our quality-weighted loss
function.

#### CNN method

The CNN-based method closely follows the approach of Ref. [21]. Here we
isolate the hyperspectral SLIC tiles and spatially zero-pad them to a common
shape of $48\times 48\times 109$ ($x\times y\times\lambda$), so that they can
be passed to a CNN model for classification. We then aggregate the CNN
predictions into a composite image segmentation result.

The CNN model uses a similar encoder architecture to Ref. [21]. We first apply
a channel compression layer, implemented with a $1\times 1$ convolutional
kernel with stride 1 and no padding, to reduce the number of spectral
dimensions from $N_{\lambda}=109$ to $N_{\mathrm{ch}}=12$. We then apply four
2D convolutional blocks bringing the number of features from $N_{\mathrm{ch}}$
to $(N_{f}/2,N_{f},2N_{f},4N_{f})$ respectively, where we take $N_{f}=12$. The
kernel sizes are respectively $(3,4,4,4)$, the strides $(1,2,2,2)$ and padding
$(1,1,1,1)$. The resulting features of dimensions $48\times 6\times 6$ are
then average-pooled, thus yielding for each tile $i$ an encoded representation
$\mathbf{z}_{i}\in\mathbb{R}^{48}$, which we pass to a fully connected
classification head. We apply a LeakyReLU activation and BatchNorm after each
convolutional layer.

During training we apply the following data augmentation: a 50% dropout before
the classification layer, a random shift along spatial axes, a random pixel
brightness shift, and random rotation, rescaling, and Gaussian blurring.

Our baseline strategy consists of training only on the high-quality tiles
(training and validation splits) obtained with the procedure described in Sec.
3.2, using a regular cross-entropy loss function summed over all tiles. At
test time, we always evaluate the model on all annotated image areas (test
split, whereby tiles are generated on the test data on the fly, so the exact
number will vary).

#### CNN with quality-weighted loss

While the exclusion of low-quality tiles at train time generally makes the
models more robust on the good-quality areas, in any realistic surgical
setting it would be desirable to have meaningful predictions over the full
image, including lower-quality areas.

To address this requirement, we incorporate all tiles to the training and
validation sets while introducing a quality-weighted loss function to
downweight low-quality tiles. We define the loss weighting factor $w_{i}$ for
each tile $T_{i}$ as

$w_{i}=w_{I}(I_{i})\cdot w_{L_{2}}(L_{2,i})\cdot
w_{\mathrm{SAM}}(\mathrm{SAM}_{i})\,,$ (3)

where the quality metrics $I$, $L_{2}$, $\mathrm{SAM}$ are defined in Sec.
3.2. We choose for the three components the forms:

$\displaystyle w_{I}(I)$ $\displaystyle=$
$\displaystyle\begin{cases}6I&\text{if}\>\>\>\>I\leq 0.167\\\
1&\text{if}\>\>\>\>0.167<I<0.5\\\ 2-2I&\text{if}\>\>\>\>I\geq 0.5\end{cases}$
$\displaystyle w_{L_{2}}(L_{2})$ $\displaystyle=$ $\displaystyle 0.7^{L_{2}}$
$\displaystyle w_{\mathrm{SAM}}(\mathrm{SAM})$ $\displaystyle=$ $\displaystyle
0.4^{\mathrm{SAM}}\,.$ (4)

With this strategy, the model during training also uses the low-quality tiles,
albeit with less weight, so that information from those regions is not lost.
At test time, we expect better generalization, especially in the regions with
high density of blood vessels, low light, and close to saturated areas.

We compare the results obtained by training on all tiles with the weighted
loss with the model trained on good tiles only, as well as with the same model
trained on all tiles with a regular cross-entropy loss.

Table 1: Quantitative comparison of model performance. We report accuracy,
$F_{1}$-score, and IoU on each class (H, T, B refer to healthy, tumor and
background respectively) and macro-averaged across all classes (Avg);
additionally, we also report the ROC AUC for the healthy/tumor binary
classification task. The top half refers to the CNN-only method, the bottom
half describes the CNN+GNN method. Models trained on good tiles only, on all
tiles, and on all tiles with a quality-weighted loss are denoted with _g, _a,
_aW respectively. The best results for each class are presented in bold. The
CNN_g model is equivalent to Ref. [21].

| Accuracy |  $F_{1}$-score | IoU | AUC
---|---|---|---|---
Model | H | T | B | Avg. | H | T | B | Avg. | H | T | B | Avg. | H vs. T
CNN_g [21] | 0.83 | 0.68 | 0.86 | 0.79 | 0.79 | 0.63 | 0.83 | 0.75 | 0.68 | 0.48 | 0.75 | 0.64 | 0.86
CNN_a | 0.85 | 0.76 | 0.94 | 0.85 | 0.82 | 0.69 | 0.92 | 0.81 | 0.71 | 0.56 | 0.71 | 0.66 | 0.88
CNN_aW | 0.84 | 0.80 | 0.97 | 0.87 | 0.85 | 0.71 | 0.88 | 0.81 | 0.74 | 0.58 | 0.73 | 0.68 | 0.91
GNN_g | 0.78 | 0.92 | 1.00 | 0.90 | 0.84 | 0.77 | 0.91 | 0.84 | 0.73 | 0.65 | 0.78 | 0.72 | 0.92
GNN_a | 0.94 | 0.64 | 0.99 | 0.86 | 0.85 | 0.67 | 0.99 | 0.84 | 0.75 | 0.50 | 0.84 | 0.69 | 0.92
GNN_aW | 0.87 | 0.83 | 0.98 | 0.89 | 0.87 | 0.74 | 0.95 | 0.85 | 0.78 | 0.60 | 0.82 | 0.73 | 0.95

#### CNN + GNN method

We improve over the CNN baseline by incorporating the spatial context of the
tiles using GNNs. For every HSI image, we construct an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertices $\mathcal{V}$ and edges
$\mathcal{E}$. Here $N=|\mathcal{V}|$ is the number of nodes, and
$A\in\mathbb{R}^{N\times N}$ is the adjacency matrix. In this context, the
vertices represent the tiles, and the edges their relational structure.

In order to characterize the nodes, we need to aggregate the tiles’ spectral
information into a matrix $X\in\mathbb{R}^{N\times F}$, where $F$ is the
dimensionality of the aggregated features per tile. To encode intra-tile
structure in both spatial and spectral domains, we choose to use our baseline
CNN model as a feature extractor, thus setting $X=Z\in\mathbb{R}^{N\times
48}$, where we aggregate the latent space representation of all tiles into one
matrix $Z$. Differently from the baseline CNN model, in this case we do not
apply batch normalization to the CNN, and we use for the final dropout layer
of the CNN model a 30% probability, as the GNN is less prone to overfitting.

We define the adjacency matrix $A$ by simply selecting for each tile its $k$
nearest neighbors with a kNN algorithm, using $k=2$; we find this small number
of neighbors to be effective in combating overfitting, which we observed in
increasing amounts for $k\in\\{3,5\\}$.

The graphs thus defined are then passed to a Graph Attention Network (GAT)
[60] with 64 hidden channels, 3 hidden GAT layer and 3 attention heads. After
each GAT layer, a ReLU activation is applied, and before the last layer we
apply a 30% dropout during training. While we experimented with different GNN
choices, we found the GAT model to be suitable for our use-case as opposed to
simpler graph-convolutional layers such as the GCN [33]. The attention
mechanism provides the flexibility of applying different weights for different
nodes in each tile’s neighborhood, which is likely to be beneficial, e.g.,
near tissue class boundary regions. We note however that a full optimization
of the GNN architecture is beyond the scope of this work.

During training, in addition to the CNN-level augmentation described above, we
apply graph-level data augmentation by randomly jittering the spatial
coordinates associated with each and by randomly dropping up to 30% of the
nodes. We use per-node cross-entropy as the loss function. While it is
possible to train the CNN+GNN model end-to-end in principle, we found that in
order to reduce overfitting, it is beneficial to not let gradients flow from
the GNN head into the CNN backbone during backpropagation, and instead
minimize the combined loss
$L_{\mathrm{comb}}=L_{\mathrm{CNN}}+L_{\mathrm{GNN}}$, where both losses are
cross-entropy.

As for the CNN case, our baseline strategy is to first train the model using
graphs constructed from the good-quality tiles only; at test time we always
evaluate the model on all annotated image areas of the 10 left-out images.
Then we include all tiles in the training by introducing the quality-weighted
loss. Also here, we will compare the results obtained by training on graphs
constructed from all tiles using weighted and unweighted losses with the
baseline model trained on graphs constructed from good tiles only.

##  4 Results

We trained all models described in Sec. 3.3 on a Linux machine with 12 CPU
cores, 64 GB RAM, and an NVIDIA GeForce RTX™ 3090 GPU (24 GB RAM).

###  4.1 Quantitative evaluation

We summarize in Tab. 1 the quantitative results of all experiments described
in Sec. 3.3 and evaluated on the test set that includes all annotated image
areas. Here we report the accuracy, $F_{1}$-score and IoU [42] metrics per
class (H, T, B for healthy, tumor, and background respectively) and macro-
averaged across all classes. Additionally, we report the integrated area under
the ROC curve (AUC) for the tumor/healthy binary classification problem, where
we excluded the background class, since it is irrelevant for the clinical
application. The ROC curves for all models are shown in Fig. 5.

In the analysis of the results, the closest previous work we can compare to is
the method by Ref. [21], which is equivalent to our CNN_g model. Other CNN-
based models build on top of the baseline CNN_g model.

Figure 5: ROC curves for each model for the binary tumor/healthy
classification task (background class was removed). The thick lines
corresponding to GNN models are always above the thin lines describing their
CNN counterparts. While the training on all tiles (_a models) with unweighted
loss is not always beneficial, the inclusion of a weighted loss (_aW models)
yields a significant improvement in all cases. The CNN_g model corresponds to
the approach by Ref. [21].

While the CNN models already reach nontrivial accuracy, surpassing 80% in most
cases, the GNN method leads to an overall superior performance. All best
results of Tab. 1 lie in the GNN section, and a like-to-like comparison shows
that the GNN approach improves average accuracy by 11 percent points (pp) for
training on good-quality tiles compared to the CNN_g baseline, by 1 pp for all
tiles, and by 2 pp in the weighted loss case.

Models trained on all tiles without weighted loss function yield contradicting
results: they show a better performance on the healthy class, but in some
cases yield degraded results on the tumor class. This behaviour can be
explained by the decreased quality of the dataset, where individual tiles
contain heterogeneous spectral information and suffer from high noise and
signal saturation. The introduction of a weighted loss improves the results
compared to the use of all tiles without weights, by rebalancing the
performance between classes. The healthy class quality decreases, while the
tumor class quality increases, leading to an average accuracy improvement of 2
pp (3 pp) for the CNN (GNN) methods. Quantitative results according to the
different metrics are largely consistent: GNN methods are superior, and the
usage of a weighted loss improves most metrics, although from the point of
view of a correct detection of the tumor class, the GNN model trained on good
tiles only is ahead.

The ROC curve of Fig. 5, together with the reported AUC, provides a good
single-number summary of each model quality for clinical applications. We
observe once again that the GNN models are superior, that the addition of a
weighted loss is beneficial, and specifically that each GNN model is ahead of
its CNN counterpart. Our best model, GNN_aW, is according to this metric 9
percent points ahead of the baseline from Ref. [21], CNN_g, with an evident
gap across all threshold levels in the ROC curve.

(a) GT

(b) CNN_g good tiles

(c) CNN_a all tiles

(d) CNN_aW weigh

(e) GNN_g good tiles

(f) GNN_a all tiles

(g) GNN_aW weigh

Figure 6: Qualitative comparison of the inference results on the test set for
the different approaches described in Sec. 3.3 and quantitatively summarized
in Tab. 1. Each row shows the inference on one test example; column (a) shows
the ground truth, and the remaining columns display the results with each
model. The grayscale images in the background represent luminosity in an
intermediate wavelength channels for illustration, and the overlayed colors
indicate the tissue classes: red, green and blue denote tumor, healthy and
background respectively.

###  4.2 Qualitative evaluation

We next compare the results qualitatively on four representative images from
the test set in Fig. 6. We show for each example the ground truth followed by
the segmentation predictions using each of the models of Tab. 1.

Both CNN and GNN models succeed in producing a qualitatively correct result in
most cases, such as the first three examples. The last example is a more
challenging sample, where we observe a substantial fraction of false
positives, which are however reduced when using more advanced approaches
(columns d, g). For all example results, the CNN models (columns b-d) are as
expected noisier, as they cannot make use of contextual information. The GNN
models (columns e-g), on the contrary, leverage context to achieve a smoother
prediction.

Models trained on good tiles only (columns b, e) feature as expected a high
error rate on lower-quality areas that were not seen during training, e.g., in
over-/under-exposed patches and in noisy areas. As also seen quantitatively in
Tab. 1, GNN models trained on all tiles (column f) display a degraded
performance on the tumor class prediction. A more robust prediction of this
class is then restored once the weighted loss is introduced (columns d, g for
CNN and GNN respectively). The qualitative evaluation confirms that our best
model, GNN_aW, brings a significant improvement compared with the baseline
from Ref. [21], CNN_g.

##  5 Conclusions and future work

In this work, we presented a novel approach for HSI-based tumor segmentation
that leverages recent advances in the field of Graph Neural Networks (GNNs) to
fully exploit HSI information in both spectral and spatial domains. We used
for the first time a composite CNN+GNN model for tile-based medical HSI image
segmentation. After the creation of macropixels (tiles) with an adapted
version of the SLIC algorithm [1], we defined for each HSI image a graph
object whose nodes are the tiles and whose edges connect neighboring tiles.
The node features are defined by applying a CNN feature extractor to the tile
data cube, thus maximising the extraction of intra-tile spatial and spectral
information. Global spatial information is in turn optimally used by passing
the graph through a GAT [60] graph model for node classification.

We evaluated this approach on a dataset of 51 HSI ex vivo images from 30
patients. We compared the segmentation results of the GNN models with a
baseline CNN by Ref. [21]. We additionally introduced a loss weighting scheme
based on local image quality that aims at recovering good segmentation results
also on lower-quality regions.

We found that all models, including the CNN baseline from Ref. [21], achieved
good accuracy ($>80\%$) on a test dataset consisting of images from unseen
patients. However, as expected, GNN models do make better use of spatial
information, thus leading to more accurate and uniform results compared with
the baseline: for models trained on good-quality image regions only, GNNs
improved the average accuracy with respect to the baseline from 79% to 90%; a
qualitative analysis also confirmed the superior performance and better
spatial uniformity of the GNN results.

We also demonstrated that our new quality-aware loss weighting scheme improves
the robustness of the predictions for both CNN and GNN models: the average
accuracy of CNN models improves from 79% to 87% with the weighted loss, and
while the overall average accuracy of GNN models does not improve, its
$F_{1}$-score does, as well as its qualitative spatial uniformity.

Our work shows that GNNs together with a quality-aware loss weighting strategy
lead to robust tissue segmentation that generalizes to unseen patients. The
next step for future research will be to apply this method to in vivo and
intra-operative settings, thus paving the way for improved surgery workflow
and patient outcome.

## References

  * [1] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi, Pascal Fua, and Sabine Süsstrunk.  Slic superpixels compared to state-of-the-art superpixel methods.  IEEE transactions on pattern analysis and machine intelligence, 34(11):2274–2282, 2012.
  * [2] Hamed Akbari, Yukio Kosugi, Kazuyuki Kojima, and Naofumi Tanaka.  Wavelet-based compression and segmentation of hyperspectral images in surgery.  In Medical Imaging and Augmented Reality: 4th International Workshop Tokyo, Japan, August 1-2, 2008 Proceedings 4, pages 142–149. Springer, 2008.
  * [3] Vijay Badrinarayanan, Alex Kendall, and Roberto Cipolla.  Segnet: A deep convolutional encoder-decoder architecture for image segmentation.  IEEE transactions on pattern analysis and machine intelligence, 39(12):2481–2495, 2017.
  * [4] John Adrian Bondy.  Graph theory with applications.  1982\.
  * [5] Fernando Cervantes-Sanchez, Marianne Maktabi, Hannes Köhler, Robert Sucher, Nada Rayes, Juan Gabriel Avina-Cervantes, Ivan Cruz-Aceves, and Claire Chalopin.  Automatic tissue segmentation of hyperspectral images in liver and head neck surgeries using machine learning.  Artif Intell Surg, 1:22–37, 2021.
  * [6] Xiaojun Chang, Pengzhen Ren, Pengfei Xu, Zhihui Li, Xiaojiang Chen, and Alex Hauptmann.  A comprehensive survey of scene graphs: Generation and application.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 45(1):1–26, 2021.
  * [7] Neil T Clancy, Geoffrey Jones, Lena Maier-Hein, Daniel S Elson, and Danail Stoyanov.  Surgical spectral imaging.  Medical image analysis, 63:101699, 2020.
  * [8] Rong Cui, He Yu, Tingfa Xu, Xiaoxue Xing, Xiaorui Cao, Kang Yan, and Jiexi Chen.  Deep learning in medical hyperspectral images: A review.  Sensors, 22(24):9790, 2022.
  * [9] Yao Ding, Zhili Zhang, Xiaofeng Zhao, Danfeng Hong, Wei Cai, Nengjun Yang, and Bei Wang.  Multi-scale receptive fields: Graph attention neural network for hyperspectral image classification.  Expert Systems with Applications, 223:119858, 2023.
  * [10] Sixu Duan, Pu Huang, Min Chen, Ting Wang, Xiaolei Sun, Meirong Chen, Xueyuan Dong, Zekun Jiang, and Dengwang Li.  Semi-supervised classification of fundus images combined with cnn and gcn.  Journal of Applied Clinical Medical Physics, 23(12):e13746, 2022\.
  * [11] Himar Fabelo, Martin Halicek, Samuel Ortega, Maysam Shahedi, Adam Szolna, Juan F Piñeiro, Coralia Sosa, Aruma J O’Shanahan, Sara Bisshopp, Carlos Espino, et al.  Deep learning-based framework for in vivo identification of glioblastoma tumor using hyperspectral images of human brain.  Sensors, 19(4):920, 2019.
  * [12] Himar Fabelo, Samuel Ortega, Silvester Kabwama, Gustavo M Callico, Diederik Bulters, Adam Szolna, Juan F Pineiro, and Roberto Sarmiento.  Helicoid project: A new use of hyperspectral imaging for brain cancer detection in real-time during neurosurgical operations.  In Hyperspectral Imaging Sensors: Innovative Applications and Sensor Standards 2016, volume 9860, page 986002. SPIE, 2016.
  * [13] Himar Fabelo, Samuel Ortega, Adam Szolna, Diederik Bulters, Juan F Piñeiro, Silvester Kabwama, Aruma JO’Shanahan, Harry Bulstrode, Sara Bisshopp, B Ravi Kiran, et al.  In-vivo hyperspectral human brain image database for brain cancer detection.  IEEE Access, 7:39098–39116, 2019.
  * [14] Baowei Fei.  Hyperspectral imaging in medical applications.  In Data Handling in Science and Technology, volume 32, pages 523–565. Elsevier, 2020.
  * [15] Jacques Ferlay, Isabelle Soerjomataram, Rajesh Dikshit, Sultan Eser, Colin Mathers, Marise Rebelo, Donald Maxwell Parkin, David Forman, and Freddie Bray.  Cancer incidence and mortality worldwide: sources, methods and major patterns in globocan 2012.  International journal of cancer, 136(5):E359–E386, 2015.
  * [16] Miranda M Fidler, Freddie Bray, Salvatore Vaccarella, and Isabelle Soerjomataram.  Assessing global transitions in human development and colorectal cancer incidence.  International journal of cancer, 140(12):2709–2715, 2017.
  * [17] Nicolás Gaggion, Lucas Mansilla, Candelaria Mosquera, Diego H Milone, and Enzo Ferrante.  Improving anatomical plausibility in medical image segmentation via hybrid graph neural networks: applications to chest x-ray analysis.  IEEE Transactions on Medical Imaging, 42(2):546–556, 2022.
  * [18] Antonio Garcia-Uceda Juarez, Raghavendra Selvan, Zaigham Saghir, and Marleen de Bruijne.  A joint 3d unet-graph neural network-based method for airway segmentation from chest cts.  In Machine Learning in Medical Imaging: 10th International Workshop, MLMI 2019, Held in Conjunction with MICCAI 2019, Shenzhen, China, October 13, 2019, Proceedings 10, pages 583–591. Springer, 2019.
  * [19] Azat Garifullin, Peeter Kööbi, Pasi Ylitepsa, Kati Ådjers, Markku Hauta-Kasari, Hannu Uusitalo, and Lasse Lensu.  Hyperspectral image segmentation of retinal vasculature, optic disc and macula.  In 2018 Digital Image Computing: Techniques and Applications (DICTA), pages 1–5. IEEE, 2018.
  * [20] Pedram Ghamisi, Javier Plaza, Yushi Chen, Jun Li, and Antonio J Plaza.  Advanced spectral classifiers for hyperspectral images: A review.  IEEE Geoscience and Remote Sensing Magazine, 5(1):8–32, 2017.
  * [21] Tommaso Giannantonio, Anna Alperovich, Piercosimo Semeraro, Manfredo Atzori, Xiaohan Zhang, Christoph Hauger, Alexander Freytag, Siri Luthman, Roeland Vandebriel, Murali Jayapala, et al.  Intra-operative brain tumor detection with deep learning-optimized hyperspectral imaging.  In Optical Biopsy XXI: Toward Real-Time Spectroscopic Imaging and Diagnosis, volume 12373, pages 74–92. SPIE, 2023.
  * [22] Martin Halicek, James D Dormer, James V Little, Amy Y Chen, and Baowei Fei.  Tumor detection of the thyroid and salivary glands using hyperspectral imaging and deep learning.  Biomedical Optics Express, 11(3):1383–1400, 2020.
  * [23] Martin Halicek, James D Dormer, James V Little, Amy Y Chen, Larry Myers, Baran D Sumer, and Baowei Fei.  Hyperspectral imaging of head and neck squamous cell carcinoma for cancer margin detection in surgical specimens from 102 patients using deep learning.  Cancers, 11(9):1367, 2019.
  * [24] Martin Halicek, Himar Fabelo, Samuel Ortega, Gustavo M Callico, and Baowei Fei.  In-vivo and ex-vivo tissue analysis through hyperspectral imaging techniques: revealing the invisible features of cancer.  Cancers, 11(6):756, 2019.
  * [25] Martin Halicek, James V Little, Xu Wang, Mihir Patel, Christopher C Griffith, Amy Y Chen, and Baowei Fei.  Tumor margin classification of head and neck cancer using hyperspectral imaging and convolutional neural networks.  In Medical Imaging 2018: Image-Guided Procedures, Robotic Interventions, and Modeling, volume 10576, pages 17–27. SPIE, 2018.
  * [26] Martin Halicek, Guolan Lu, James V Little, Xu Wang, Mihir Patel, Christopher C Griffith, Mark W El-Deiry, Amy Y Chen, and Baowei Fei.  Deep convolutional neural networks for classifying head and neck cancer using hyperspectral imaging.  Journal of biomedical optics, 22(6):060503–060503, 2017.
  * [27] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  Advances in neural information processing systems, 30, 2017.
  * [28] Lin He, Jun Li, Chenying Liu, and Shutao Li.  Recent advances on spectral–spatial hyperspectral image classification: An overview and new guidelines.  IEEE Transactions on Geoscience and Remote Sensing, 56(3):1579–1597, 2017.
  * [29] Danfeng Hong, Lianru Gao, Jing Yao, Bing Zhang, Antonio Plaza, and Jocelyn Chanussot.  Graph convolutional networks for hyperspectral image classification.  IEEE Transactions on Geoscience and Remote Sensing, 59(7):5966–5978, 2020.
  * [30] Steven L Jacques.  Optical properties of biological tissues: a review.  Physics in Medicine & Biology, 58(11):R37, 2013.
  * [31] Simon Jégou, Michal Drozdzal, David Vazquez, Adriana Romero, and Yoshua Bengio.  The one hundred layers tiramisu: Fully convolutional densenets for semantic segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pages 11–19, 2017.
  * [32] Uzair Khan, Sidike Paheding, Colin P Elkin, and Vijaya Kumar Devabhaktuni.  Trends in deep learning for medical hyperspectral image analysis.  IEEE Access, 9:79534–79548, 2021.
  * [33] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907, 2016.
  * [34] Raquel Leon, Sofia H Gelado, Himar Fabelo, Samuel Ortega, Laura Quintana, Adam Szolna, Juan F Piñeiro, Francisco Balea-Fernandez, Jesus Morera, Bernardino Clavo, et al.  Hyperspectral imaging for in-vivo/ex-vivo tissue analysis of human brain cancer.  In Medical Imaging 2022: Image-Guided Procedures, Robotic Interventions, and Modeling, volume 12034, pages 525–534. SPIE, 2022.
  * [35] Shutao Li, Weiwei Song, Leyuan Fang, Yushi Chen, Pedram Ghamisi, and Jon Atli Benediktsson.  Deep learning for hyperspectral image classification: An overview.  IEEE Transactions on Geoscience and Remote Sensing, 57(9):6690–6709, 2019.
  * [36] Bing Liu, Kuiliang Gao, Anzhu Yu, Wenyue Guo, Ruirui Wang, and Xibing Zuo.  Semisupervised graph convolutional network for hyperspectral image classification.  Journal of Applied Remote Sensing, 14(2):026516–026516, 2020.
  * [37] Guolan Lu and Baowei Fei.  Medical hyperspectral imaging: a review.  Journal of biomedical optics, 19(1):010901–010901, 2014.
  * [38] Francesca Manni, Fons van der Sommen, Himar Fabelo, Svitlana Zinger, Caifeng Shan, Erik Edström, Adrian Elmi-Terander, Samuel Ortega, Gustavo Marrero Callicó, and Peter HN de With.  Hyperspectral imaging for glioblastoma surgery: Improving tumor identification using a deep spectral-spatial approach.  Sensors, 20(23):6955, 2020.
  * [39] Shanthi Marur and Arlene A Forastiere.  Head and neck squamous cell carcinoma: update on epidemiology, diagnosis, and treatment.  In Mayo Clinic Proceedings, volume 91, pages 386–396. Elsevier, 2016.
  * [40] Sara Moccia, Sebastian J Wirkert, Hannes Kenngott, Anant S Vemuri, Martin Apitz, Benjamin Mayer, Elena De Momi, Leonardo S Mattos, and Lena Maier-Hein.  Uncertainty-aware organ classification for surgical data science applications in laparoscopy.  IEEE Transactions on Biomedical Engineering, 65(11):2649–2659, 2018\.
  * [41] Lichao Mou, Xiaoqiang Lu, Xuelong Li, and Xiao Xiang Zhu.  Nonlocal graph convolutional networks for hyperspectral image classification.  IEEE Transactions on Geoscience and Remote Sensing, 58(12):8246–8257, 2020.
  * [42] Dominik Müller, Iñaki Soto-Rey, and Frank Kramer.  Towards a guideline for evaluation metrics in medical image segmentation.  BMC Research Notes, 15(1):1–8, 2022.
  * [43] Samuel Ortega, Martin Halicek, Himar Fabelo, Gustavo M Callico, and Baowei Fei.  Hyperspectral and multispectral imaging in digital and computational pathology: a systematic review.  Biomedical Optics Express, 11(6):3195–3233, 2020.
  * [44] Samuel Ortega, Martin Halicek, Himar Fabelo, Rafael Camacho, María de la Luz Plaza, Fred Godtliebsen, Gustavo M. Callicó, and Baowei Fei.  Hyperspectral imaging for the detection of glioblastoma tumor cells in h&e slides using convolutional neural networks.  Sensors, 20(7):1911, 2020.
  * [45] Samuel Ortega, Martin Halicek, Himar Fabelo, Eduardo Quevedo, Baowei Fei, and Gustavo Marrero Callico.  Information extraction techniques in hyperspectral imaging biomedical applications.  In Multimedia Information Retrieval. IntechOpen, 2020.
  * [46] D Max Parkin, Freddie Bray, J Ferlay, and Paola Pisani.  Global cancer statistics, 2002.  CA: a cancer journal for clinicians, 55(2):74–108, 2005.
  * [47] Dhrumil Patel, Dhruv Patel, Rudra Saxena, and Thangarajah Akilan.  Multi-class brain tumor segmentation using graph attention network.  arXiv preprint arXiv:2302.05598, 2023.
  * [48] Anyong Qin, Zhaowei Shang, Jinyu Tian, Yulong Wang, Taiping Zhang, and Yuan Yan Tang.  Spectral–spatial graph convolutional networks for semisupervised hyperspectral image classification.  IEEE Geoscience and Remote Sensing Letters, 16(2):241–245, 2018\.
  * [49] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * [50] Camillo Saueressig, Adam Berkley, Elliot Kang, Reshma Munbodh, and Ritambhara Singh.  Exploring graph-based neural networks for automatic brain tumor segmentation.  In From Data to Models and Back: 9th International Symposium, DataMod 2020, Virtual Event, October 20, 2020, Revised Selected Papers 9, pages 18–37. Springer, 2021.
  * [51] Camillo Saueressig, Adam Berkley, Reshma Munbodh, and Ritambhara Singh.  A joint graph and image convolution network for automatic brain tumor segmentation.  In International MICCAI Brainlesion Workshop, pages 356–365. Springer, 2021.
  * [52] Silvia Seidlitz, Jan Sellner, Jan Odenthal, Berkin Özdemir, Alexander Studier-Fischer, Samuel Knödler, Leonardo Ayala, Tim J Adler, Hannes G Kenngott, Minu Tizabi, et al.  Robust deep learning-based semantic organ segmentation in hyperspectral images.  Medical Image Analysis, 80:102488, 2022.
  * [53] Anshu Sha, Bin Wang, Xiaofeng Wu, and Liming Zhang.  Semisupervised classification for hyperspectral images using graph attention networks.  IEEE Geoscience and Remote Sensing Letters, 18(1):157–161, 2020\.
  * [54] Anshu Sha, Bin Wang, Xiaofeng Wu, Liming Zhang, Bo Hu, and Jian Qiu Zhang.  Semi-supervised classification for hyperspectral images using edge-conditioned graph convolutional networks.  In IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium, pages 2690–2693. IEEE, 2019.
  * [55] Farideh Foroozandeh Shahraki and Saurabh Prasad.  Graph convolutional neural networks for hyperspectral data classification.  In 2018 IEEE global conference on signal and information processing (GlobalSIP), pages 968–972. IEEE, 2018.
  * [56] Hua-Chieh Shao, Jing Wang, Ti Bai, Jaehee Chun, Justin C Park, Steve Jiang, and You Zhang.  Real-time liver tumor localization via a single x-ray projection using deep graph neural network-assisted biomechanical modeling.  Physics in Medicine & Biology, 67(11):115009, 2022.
  * [57] Stojan Trajanovski, Caifeng Shan, Pim JC Weijtmans, Susan G Brouwer de Koning, and Theo JM Ruers.  Tumor semantic segmentation in hyperspectral images using deep learning.  In International Conference on Medical Imaging with Deep Learning–Extended Abstract Track, 2019.
  * [58] Stojan Trajanovski, Caifeng Shan, Pim JC Weijtmans, Susan G Brouwer de Koning, and Theo JM Ruers.  Tongue tumor detection in hyperspectral images using deep learning semantic segmentation.  IEEE transactions on biomedical engineering, 68(4):1330–1340, 2020\.
  * [59] Gemma Urbanos, Alberto Martín, Guillermo Vázquez, Marta Villanueva, Manuel Villa, Luis Jimenez-Roldan, Miguel Chavarrías, Alfonso Lagares, Eduardo Juárez, and César Sanz.  Supervised machine learning methods and hyperspectral imaging techniques jointly applied for brain cancer classification.  Sensors, 21(11):3827, 2021.
  * [60] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * [61] Sheng Wan, Chen Gong, Ping Zhong, Bo Du, Lefei Zhang, and Jian Yang.  Multiscale dynamic graph convolutional network for hyperspectral image classification.  IEEE Transactions on Geoscience and Remote Sensing, 58(5):3162–3177, 2019.
  * [62] Sheng Wan, Chen Gong, Ping Zhong, Shirui Pan, Guangyu Li, and Jian Yang.  Hyperspectral image classification with context-aware dynamic graph convolutional network.  IEEE Transactions on Geoscience and Remote Sensing, 59(1):597–612, 2020.
  * [63] Qian Wang, Li Sun, Yan Wang, Mei Zhou, Menghan Hu, Jiangang Chen, Ying Wen, and Qingli Li.  Identification of melanoma from hyperspectral pathology image using 3d convolutional networks.  IEEE Transactions on Medical Imaging, 40(1):218–227, 2020.
  * [64] Tinghuai Wang, Guangming Wang, Kuan Eeik Tan, and Donghui Tan.  Spectral pyramid graph attention network for hyperspectral image classification.  arXiv preprint arXiv:2001.07108, 2020.
  * [65] Xi Wang and Bin-bin Li.  Deep learning in head and neck tumor multiomics diagnosis and analysis: review of the literature.  Frontiers in Genetics, 12:624820, 2021.
  * [66] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  IEEE transactions on neural networks and learning systems, 32(1):4–24, 2020.
  * [67] Zhang Yan, Kong Youyong, Wu Jiasong, Gouenou Coatrieux, and Shu Huazhong.  Brain tissue segmentation based on graph convolutional networks.  In 2019 IEEE International Conference on Image Processing (ICIP), pages 1470–1474. IEEE, 2019.
  * [68] Pan Yang, Lei Tong, Bin Qian, Zheng Gao, Jing Yu, and Chuangbai Xiao.  Hyperspectral image classification with spectral and spatial graph using inductive representation learning network.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 14:791–800, 2020.
  * [69] Ding Yao, Zhang Zhi-li, Zhao Xiao-feng, Cai Wei, He Fang, Cai Yao-ming, and Wei-Wei Cai.  Deep hybrid: multi-graph neural network collaboration for hyperspectral image classification.  Defence Technology, 23:164–176, 2023.
  * [70] Boxiang Yun, Yan Wang, Jieneng Chen, Huiyu Wang, Wei Shen, and Qingli Li.  Spectr: Spectral transformer for hyperspectral pathology image segmentation.  arXiv preprint arXiv:2103.03604, 2021.
  * [71] Chaozi Zhang, Jianli Wang, and Kainan Yao.  Global random graph convolution network for hyperspectral image classification.  Remote Sensing, 13(12):2285, 2021.