An Integrated Deep Learning and Statistical Framework for Whole-Network Gene--Environment Association with Leaf Vascular Architecture

# An Integrated Deep Learning and Statistical Framework for Whole-Network Gene–Environment Association with Leaf Vascular Architecture
Source: [https://arxiv.org/html/2607.22763](https://arxiv.org/html/2607.22763)
Geran ZhaoYangsheng WangDepartment of Mathematics and Statistics, Binghamton University, Binghamton, New York 13902, USAXiaotian DaiDepartment of Mathematics, Illinois State University, Normal, Illinois 61790, USAGuifang FuCorresponding author: gfu@binghamton\.eduDepartment of Mathematics and Statistics, Binghamton University, Binghamton, New York 13902, USA

###### Abstract

Leaf veins exhibit remarkable diversity in architecture and patterning, yet existing gene–environment association studies have primarily quantified leaf venation using a small collection of low\-dimensional summary traits, thereby discarding most of the structural information contained in the original images\. We propose an integrated deep learning and statistical framework that, to the best of our knowledge, is the first to investigate gene–environment associations by quantifying leaf vascular architecture as a whole\-network\. The proposed framework achieves four methodological advances\. First, it represents the complete leaf vascular architecture as a whole\-network image phenotype, preserving both fine\-scale local details and global topological structures\. Second, it fine\-tunes the deep learning\-based Edge Detection with Transformers \(EDTER\) model to accurately extract whole\-network leaf vascular architecture from RGB images by jointly learning local and global contextual features\. Compared with traditional image\-processing techniques, EDTER reduces the need for specialized image preparation while maintaining robust performance on low\-quality images with weak vein contrast\. Third, it constructs a new annotated leaf image database by integrating edge maps generated by DiffusionEdge with the Berkeley Segmentation Database \(BSDS500\), providing a bridge between existing benchmark datasets and user\-generated image collections for domain adaptation and specialized edge\-detection tasks\. Fourth, it applies Semiparametric Sparse Canonical Correlation Analysis \(SSCCA\) to perform variable selection and model associations between repeatedly measured high\-dimensional Bivariate image responses and high\-dimensional predictors while simultaneously accommodating sparse, zero\-inflated data represented by edge maps through a truncated latent Gaussian copula model\. Two simulation studies demonstrate the performance of the proposed framework under increasing levels of complexity\. Application to a real*Populus*dataset identifies three significant gene–geography interactions associated with leaf vascular architecture, providing new biological insights and establishing a broadly applicable methodological framework for high\-dimensional complex image phenotypes\.

Keywords:Edge Detection, Leaf Vascular Architecture, Deep Learning, Semiparametric Sparse Canonical Correlation Analysis, Genotype–Phenotype Associations

## 1Introduction

Leaf veins exhibit exceptional diversity in architecture and patterning, representing one of the most striking examples of vascular network variation in nature\[[24](https://arxiv.org/html/2607.22763#bib.bib43)\]\. Veins are typically organized hierarchically into first\-order \(midvein\), second\-order \(lateral veins\), third\-order, and additional subtle minor veins\[[24](https://arxiv.org/html/2607.22763#bib.bib43)\]\. This diversity is manifested in leaf function, hydraulic transport, mechanical support, and evolutionary adaptation\. Specifically, leaf veins serve as the primary pathways for water transport, thereby influencing hydraulic function, stomatal behavior, maximum photosynthetic rate, drought tolerance, and hydraulic safety margins\[[3](https://arxiv.org/html/2607.22763#bib.bib27)\]\. Moreover, leaf veins contribute to the mechanical toughness of the lamina by serving as the structural barriers that resist fracture\[[17](https://arxiv.org/html/2607.22763#bib.bib3)\]\. Leaf venation is also linked to broader ecological adaptations, including herbivory resistance and variation in leaf lifespan\[[19](https://arxiv.org/html/2607.22763#bib.bib59)\]\. Existing studies have demonstrated that genes play a critical role in regulating leaf vascular architecture, which is highly heritable and provides a unique opportunity to investigate natural variation and uncover genotype–phenotype relationships\[[23](https://arxiv.org/html/2607.22763#bib.bib53),[18](https://arxiv.org/html/2607.22763#bib.bib61)\]\.

Understanding how genes inlfuence leaf vascular architecture requires not only precise quantification of leaf venation but also accurate downstream statistical modeling, both of which involve high\-dimensional phenotypic and genetic data\. Leaf vein extraction is particularly challenging because vein structures often appear as fine, subtle, and low\-contrast internal edges that are difficult to detect accurately\. Existing studies have quantified leaf venation using four broad categories of traits\[[23](https://arxiv.org/html/2607.22763#bib.bib53),[18](https://arxiv.org/html/2607.22763#bib.bib61)\]: structural, topological, geometric, and mechanical traits\. Structural traits, such as vein density \(VD\), midvein area \(MVA\), vein number \(VN\), vein diameter, and hierarchical vein orders, primarily determine water transport efficiency\[[24](https://arxiv.org/html/2607.22763#bib.bib43)\]\. Topological traits, including areole density, loopiness, and connectivity, influence hydraulic redundancy and damage resilience\[[2](https://arxiv.org/html/2607.22763#bib.bib31)\]\. Geometric traits, such as vein length, orientation, and branching angles, describe the spatial organization of leaf venation, and influence leaf folding patterns and hydraulic pathways\[[21](https://arxiv.org/html/2607.22763#bib.bib38)\]\. Mechanical and anatomical traits, including bundle sheath dimensions, sclerenchyma investment, and leaf toughness, contribute to leaf durability, herbivory resistance, and leaf lifespan\[[20](https://arxiv.org/html/2607.22763#bib.bib39)\]\. However, as low\-dimensional summaries extracted from the original images, they provided only a partial representation of the vascular architecture, thereby discarding most of the rich structural information contained in the full venascular architecture\.

Several well\-established tools have been developed to extract leaf vein traits from images, including phenoVein, LIMANI, LEAF GUI, and NEFI\[[21](https://arxiv.org/html/2607.22763#bib.bib38),[5](https://arxiv.org/html/2607.22763#bib.bib37),[4](https://arxiv.org/html/2607.22763#bib.bib49),[6](https://arxiv.org/html/2607.22763#bib.bib46)\]\. These tools extracted low\-dimensional venation traits such as vein length, vein density, branching points, and areole area, etc\. They primarily relied on image\-processing techniques involving segmentation, skeletonization, and graph reconstruction rather than learning\-based approaches\. Consequently, their performance was often sensitive to image quality and performed best on images with high vein contrast\. Segmentation errors can lead to disconnected veins, false branches, or inaccurate trait extraction\. Moreover, some methods required specialized image preparation\. For example, LIMANI was developed for chemically cleared leaves, in which pigments and soft tissues were chemically removed to enhance the recognition of the leaf vein traits\.

Few learning\-based approaches have been developed for leaf venation analysis\.Lagergrenet al\.\[[12](https://arxiv.org/html/2607.22763#bib.bib62)\]employed a few\-shot convolutional neural network \(CNN\) to segment vein structure from high\-resolution RGB leaf images\. However, they still subsequently extracted low\-dimensional summary traits, such as vein length and vein density, from the segmented vein images using RhizoVision Explorer \(RVE\), rather than directly modeling the full venation architecture\.

After low\-dimensional leaf vein traits have been extracted, downstream gene–vein association analyses in the existing literature primarily relied on mixed linear models \(MLM\), fixed and random model circulating probability unification \(FarmCPU\), and Bayesian\-information and linkage\-disequilibrium iteratively nested keyway \(BLINK\)\[[32](https://arxiv.org/html/2607.22763#bib.bib20),[15](https://arxiv.org/html/2607.22763#bib.bib52),[11](https://arxiv.org/html/2607.22763#bib.bib54)\]\.Rishmawiet al\.\[[23](https://arxiv.org/html/2607.22763#bib.bib53)\]employed MLM and Markov models to identify candidate genes associated with VD, areole number, and the number of vein endpoints in*Arabidopsis thaliana*\.Narawatthanaet al\.\[[18](https://arxiv.org/html/2607.22763#bib.bib61)\]applied MLM, FarmCPU, BLINK, and haplotype analysis to identify genetic associations with the distance between minor veins \(IVD\) in rice\. While these statistical approaches successfully identified gene–vein associations, they primarily focused on low\-dimensional traits, and were not directly applicable to high\-dimensional image traits containing the whole\-network\. Furthermore, these approaches were largely restricted by parametric structures or normality assumptions\.

In this article, we propose an integrated deep learning and statistical framework that achieves four methodological advances driven by gene–vein association analysis\. First, we quantify the complete leaf vascular architecture as a whole\-network phenotype rather than a small collection of low\-dimensional summary traits\. Unlike conventional approaches, the proposed framework preserves the complete structural information of the vascular network, including both fine\-scale local details and global network topology, and creates a new framework for investigating gene–environment associations with the leaf vascular whole\-network phenotype\.

Second, we fine\-tune the Edge Detection with Transformers \(EDTER\), a deep learning\-based approach that overcomes several limitations of traditional handcrafted image\-processing techniques by reducing the need for specialized image preparation and maintaining robust performance on low\-quality images with weak vein contrast\. EDTER takes RGB images of the leaf lamina as input and outputs a pixel\-level edge map of the same spatial dimension by jointly learning global and local contextual features\[[22](https://arxiv.org/html/2607.22763#bib.bib75)\]\. CNN\-based methods, such as Holistically\-Nested Edge Detection \(HED\) and Richer Convolutional Features \(RCF\), improved upon traditional image\-processing approaches by learning multi\-scale hierarchical feature representations\[[29](https://arxiv.org/html/2607.22763#bib.bib72),[16](https://arxiv.org/html/2607.22763#bib.bib73)\]\. However, CNNs primarily relied on local receptive fields and therefore had limited ability to capture long\-range spatial dependencies\. This limitation is particularly relevant for leaf vein extraction, where vein networks may extend across large spatial regions and exhibit complex global connectivity patterns\. Transformer\-based edge detection methods address this challenge through self\-attention mechanisms that directly model long\-range dependencies and global contextual information\.

Third, the widely used edge\-detection benchmark, the Berkeley Segmentation Database \(BSDS500\)\[[1](https://arxiv.org/html/2607.22763#bib.bib68)\], contains relatively few leaf images, which may limit the ability of the pretrained EDTER models to accurately identify fine leaf veins\. To address this limitation, we construct a new leaf image database and generate high\-quality annotated edge maps using DiffusionEdge, a state\-of\-the\-art edge\-detection framework based on Diffusion Probabilistic Models\[[30](https://arxiv.org/html/2607.22763#bib.bib78)\]\. Although the edge maps generated by DiffusionEdge are not manually annotated ground\-truth labels, they provide an efficient and scalable source of pseudo\-ground\-truth supervision\. This strategy offers a practical alternative to labor\-intensive manual annotation while enabling rapid construction of domain\-specific training datasets\. More importantly, it establishes a bridge between existing benchmark datasets and user\-generated image collections, facilitating domain adaptation and fine\-tuning for specialized edge\-detection tasks\.

Fourth, instead of relying on parametric statistical models for low\-dimensional trait associations, we apply Semiparametric Sparse Canonical Correlation Analysis \(SSCCA\) to model the relationships between repeatedly measured high\-dimensional Bivariate images and high\-dimensional predictors\. The edge maps are high\-dimensional images in which the majority of pixels correspond to background and only a small fraction represent vein structures\. After normalizing pixel intensities to the interval\[0,1\]\[0,1\], most background pixels take values of zero or near zero, resulting in highly sparse and zero\-inflated image responses\. SSCCA incorporates a truncated latent Gaussian copula model to estimate the correlation structure of high\-dimensional zero\-inflated data\[[31](https://arxiv.org/html/2607.22763#bib.bib79)\]\.

We conduct two simulation studies with increasing levels of complexity to evaluate the performance of the proposed framework\. We further apply the proposed framework to a real*Populus*dataset consisting of 100 leaves, each photographed from both the front and back sides, resulting in 200 RGB leaf images, together with 118 genetic\- and geography\-related predictor variables\. Because the front and back sides of the same leaf share the same genetic and geographic information, the corresponding edge maps are treated as repeatedly measured Bivariate images for each subject\. To the best of our knowledge, this is the first study to employ deep learning techniques to represent leaf vascular architecture as a whole\-network image phenotype for gene–environment association analysis, and the identified associations provide new biological insights into genotype–phenotype relationships\. The proposed framework provides a general methodology for the analysis of broader network\-structured and image\-based data\.

## 2Methods

### 2\.1The Structure of EDTER

Edge detection methods are commonly designed to identify object boundaries and internal edges in images\. EDTER is a Transformer\-based edge detector that combines global contextual information with local fine\-grained cues through a two\-stage architecture\. Specifically, Stage I is designed to learn global contextual information by employing a global transformer encoder and a global bidirectional multi\-level aggregation \(BiMLA\) decoder to generate high\-resolution feature representations\. Stage II focuses on learning short\-range local cues and applies a local transformer encoder and a local BiMLA module to produce pixel\-level feature maps\. The outputs from both stages are integrated through a Feature Fusion Module \(FFM\) to predict the final edge maps\[[22](https://arxiv.org/html/2607.22763#bib.bib75)\]\. See Figure[1](https://arxiv.org/html/2607.22763#S2.F1)for a schematic overview of EDTER\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/EDTER_Architecture.png)Figure 1:A schematic overview of EDTER\.El​o​c​a​lE\_\{local\}denotes the final edge map produced by EDTER through the fusion of global and local feature representations\.#### 2\.1\.1Transformer Encoder

In the EDTER framework, the transformer encoder adopts the Vision Transformer \(ViT\) architecture, which represents an image as a sequence of embedded image patches and processes them through multiple transformer blocks\[[7](https://arxiv.org/html/2607.22763#bib.bib76)\]\. For Stage I, the input RGB imageI∈ℝH×W×3I\\in\\mathbb\{R\}^\{H\\times W\\times 3\}is partitioned into a sequence of non\-overlapping16×1616\\times 16image patches, whereHHandWWdenote the height and width of the image\. For Stage II, the same input image is first divided into four non\-overlapping windows, each of sizeH2×W2\\frac\{H\}\{2\}\\times\\frac\{W\}\{2\}\. Each window is then partitioned into a sequence of non\-overlapping8×88\\times 8image patches\. The resulting patches are flattened and projected into a latent embedding space through a learnable linear projection\. Positional embeddings are subsequently added to preserve spatial information\.

The transformer encoder consists of a sequence of transformer blocks, each composed of Layer Normalizations \(LayerNorm\), a Multi\-head Self\-Attention \(MSA\) module, residual connections, and a Multi\-Layer Perceptron \(MLP\)\. LetMMdenote the number of attention heads andLLdenote the number of transformer blocks\. After patch embedding and positional embedding, the initial token embedding is denoted byZ0Z\_\{0\}\. LetZℓZ\_\{\\ell\}denote the output representation of theℓ\\ellth transformer block\. The transformer blocks are connected recursively throughZℓ=Blockℓ⁡\(Zℓ−1\),ℓ=1,…,L,Z\_\{\\ell\}=\\operatorname\{Block\}\_\{\\ell\}\(Z\_\{\\ell\-1\}\),~\\ell=1,\\ldots,L,whereZℓ∈ℝN×JZ\_\{\\ell\}\\in\\mathbb\{R\}^\{N\\times J\},NNis the number of image patches \(tokens\), andJJis the embedding dimension of each token\. Here,Blockℓ⁡\(⋅\)\\operatorname\{Block\}\_\{\\ell\}\(\\cdot\)denotes the completeℓ\\ellth transformer block, consisting of LayerNorm, an MSA, residual connections, and an MLP\.

Within theℓ\\ellth transformer block, Layer Normalization is first applied to the input,Z~ℓ−1=LayerNorm⁡\(Zℓ−1\)\.\\widetilde\{Z\}\_\{\\ell\-1\}=\\operatorname\{LayerNorm\}\(Z\_\{\\ell\-1\}\)\.The normalized token embeddings serve as the common input to allMMattention heads\. For themmth attention head \(m=1,…,Mm=1,\\ldots,M\), the query, key, and value matrices are computed as

Qℓ\(m\)\\displaystyle Q\_\{\\ell\}^\{\(m\)\}=Z~ℓ−1​WQ,ℓ\(m\),\\displaystyle=\\widetilde\{Z\}\_\{\\ell\-1\}W\_\{Q,\\ell\}^\{\(m\)\},Kℓ\(m\)\\displaystyle K\_\{\\ell\}^\{\(m\)\}=Z~ℓ−1​WK,ℓ\(m\),\\displaystyle=\\widetilde\{Z\}\_\{\\ell\-1\}W\_\{K,\\ell\}^\{\(m\)\},Vℓ\(m\)\\displaystyle V\_\{\\ell\}^\{\(m\)\}=Z~ℓ−1​WV,ℓ\(m\),\\displaystyle=\\widetilde\{Z\}\_\{\\ell\-1\}W\_\{V,\\ell\}^\{\(m\)\},whereWQ,ℓ\(m\),WK,ℓ\(m\),WV,ℓ\(m\)∈ℝJ×dW\_\{Q,\\ell\}^\{\(m\)\},W\_\{K,\\ell\}^\{\(m\)\},W\_\{V,\\ell\}^\{\(m\)\}\\in\\mathbb\{R\}^\{J\\times d\}are learnable projection matrices for each of the attention head in theℓ\\ellth transformer block, anddddenotes the dimension of each attention head\.

The output of themmth attention head in theℓ\\ellth transformer block is computed as

Oℓ\(m\)=softmax⁡\(Qℓ\(m\)​\(Kℓ\(m\)\)⊤d\)​Vℓ\(m\),O\_\{\\ell\}^\{\(m\)\}=\\operatorname\{softmax\}\\\!\\left\(\\frac\{Q\_\{\\ell\}^\{\(m\)\}\(K\_\{\\ell\}^\{\(m\)\}\)^\{\\top\}\}\{\\sqrt\{d\}\}\\right\)V\_\{\\ell\}^\{\(m\)\},which are concatenated across allMMattention heads and projected back to the embedding dimension,

Oℓ=\[Oℓ\(1\),Oℓ\(2\),…,Oℓ\(M\)\]​WO,ℓ,O\_\{\\ell\}=\\left\[O\_\{\\ell\}^\{\(1\)\},O\_\{\\ell\}^\{\(2\)\},\\ldots,O\_\{\\ell\}^\{\(M\)\}\\right\]W\_\{O,\\ell\},whereWO,ℓ∈ℝM​d×JW\_\{O,\\ell\}\\in\\mathbb\{R\}^\{Md\\times J\}is the projection matrix of theℓ\\ellth MSA\.

Note thatOℓO\_\{\\ell\}is the output of a MSA module and it is subsequently processed through residual connections, LayerNorm, and an MLP to produce the output representationZℓZ\_\{\\ell\}of theℓ\\ellth transformer block\. In this study, we setM=16M=16\. The global transformer encoder consists ofL=24L=24transformer blocks, producing the sequence of output representations\{Z1,…,Z24\}\\\{Z\_\{1\},\\ldots,Z\_\{24\}\\\}, whereas the local transformer encoder consists ofL=12L=12transformer blocks, producing\{T1,…,T12\}\.\\\{T\_\{1\},\\ldots,T\_\{12\}\\\}\.

#### 2\.1\.2Bi\-directional Multi\-Level Aggregation Decoder

The global transformer encoder produces 24 transformer blocks, which are evenly divided into four groups\. The last block from each group, namelyZ6Z\_\{6\},Z12Z\_\{12\},Z18Z\_\{18\}, andZ24Z\_\{24\}, is selected as the input to the BiMLA decoder\. Each output is reshaped from anN×JN\\times Jmatrix into a three\-dimensional feature map of sizeH/16×W/16×JH/16\\times W/16\\times J, restoring the spatial arrangement of the image patches\. A subsequent1×11\\times 1convolution aligns the feature channels before feature aggregation\. As illustrated in Figure[2](https://arxiv.org/html/2607.22763#S2.F2), the top\-down path propagates high\-level semantic information through successive feature fusion operations, whereas the bottom\-up path propagates low\-level spatial details while preserving fine edge information\. Together, these complementary paths effectively aggregate multi\-scale contextual information\.

To recover spatial resolution, deconvolution layers with4×44\\times 4and16×1616\\times 16kernels are employed to upsample feature maps from different scales\. Each deconvolution layer is followed by Batch Normalization \(BN\) and a ReLU activation function\. The upsampled feature maps are subsequently concatenated and refined by three consecutive3×33\\times 3convolution layers followed by a final1×11\\times 1convolution layer, producing the global edge feature mapfg​l​o​b​a​lf\_\{global\}\.

Figure[2](https://arxiv.org/html/2607.22763#S2.F2)illustrates a schematic overview of the global BiMLA decoder\. The local BiMLA decoder follows the same overall design, except that the3×33\\times 3convolution layers are replaced by1×11\\times 1convolution layers to reduce the introduction of artificial edges while preserving local edge details\. The local BiMLA takesT3T\_\{3\},T6T\_\{6\},T9T\_\{9\}, andT12T\_\{12\}as its inputs and outputs edge feature mapfl​o​c​a​lf\_\{local\}\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/BiMLA_figure.png)Figure 2:A schematic overview of the Bi\-directional Multi\-Level Aggregation Decoder\.
#### 2\.1\.3Feature Fusion Module

To integrate the complementary information extracted from the global and local transformer encoders, EDTER employs a Feature Fusion Module, illustrated in Figure[3](https://arxiv.org/html/2607.22763#S2.F3)\. The FFM is built upon a Spatial Feature Transform \(SFT\) layer\[[27](https://arxiv.org/html/2607.22763#bib.bib74)\], followed by convolutional, Batch Normalization, and ReLU operations\. The global feature mapfg​l​o​b​a​lf\_\{global\}is first processed by two parallel convolution branches to generate two sets of spatially adaptive modulation parameters\. These parameters provide contextual guidance for refining the local feature representation and are applied to the local feature mapfl​o​c​a​lf\_\{local\}through element\-wise multiplication followed by element\-wise addition\. The resulting fused feature map is then processed by an additional3×33\\times 3convolutional layer, followed by BN and ReLU, before being passed to the local decision head to predict the refined edge map\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/FFM.png)Figure 3:A schematic overview of the Feature Fusion Module\.

### 2\.2Fine\-tuning with Enlarged New Leaf Database

We construct a new leaf image database consisting of 324 selected healthy leaves from the Mendeley Leaf Image Database\[[25](https://arxiv.org/html/2607.22763#bib.bib81)\]\. However, these images do not contain manually annotated edge maps, which are required as ground truth labels for fine\-tuning EDTER\. Unlike deterministic edge\-detection methods, DiffusionEdge formulates edge detection as a conditional image generation task, progressively refining edge maps from Gaussian noise under the guidance of the original RGB image\. Figure[4](https://arxiv.org/html/2607.22763#S2.F4)presents four examples from the new leaf database together with the corresponding edge maps annotated by DiffusionEdge, which exhibit sharp and thin vein structures that visually resemble the manually annotated edge maps provided in BSDS500\.

To improve the model’s adaptation to leaf venation images, we combine this dataset with BSDS500 and utilize the annotated edge maps as supervisory labels for fine\-tuning EDTER\[[1](https://arxiv.org/html/2607.22763#bib.bib68)\]\. We initialize the model with pre\-trained weights provided by EDTER and further train it for 20,000 more iterations for stage I and 40,000 more iterations for stage II\. The fine\-tuned EDTER model is subsequently applied to the leaf images in the real dataset to generate venation edge maps for downstream association analysis\.

Figure[5](https://arxiv.org/html/2607.22763#S2.F5)compares the edge maps generated by the pretrained and fine\-tuned EDTER models\. The fine\-tuned EDTER model \(right panel\) preserves substantially more minor veins and fine\-scale branching structures, resulting in a more complete representation of the whole\-network\. In addition, vein segments of the right panel appear more continuous and better connected than the left panel, reducing fragmentation in the extracted venation patterns\. These improvements are particularly evident in regions with dense tertiary and quaternary venation, where many faint vein segments become clearly noticeable after fine\-tuning\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/leaf_data_example.jpg)Figure 4:Four randomly selected examples from the newly constructed leaf image dataset, together with the corresponding RGB images and edge maps annotated by DiffusionEdge\.![Refer to caption](https://arxiv.org/html/2607.22763v1/Fine_tuned_comparision_v2.jpg)Figure 5:Comparison of the EDTER model before and after fine\-tuning\. The left panel illustrates the output of the pre\-trained EDTER model, whereas the right panel shows the output of the fine\-tuned EDTER model\.
### 2\.3Semiparametric Sparse Canonical Correlation Analysis

Let𝐀\(1\),𝐀\(2\)∈ℝH×W\\mathbf\{A\}^\{\(1\)\},\\mathbf\{A\}^\{\(2\)\}\\in\\mathbb\{R\}^\{H\\times W\}denote the edge maps extracted by EDTER from the front and back sides of the same leaf, respectively\. The original edge maps produced by EDTER assign pixel values close to 0 to vein structures and values close to 1 to background regions\. To improve interpretability of the downstream statistical analysis, we invert the pixel intensities so that vein pixels take values close to 1, whereas background pixels take values close to 0\. Throughout this article,𝐀\(1\)\\mathbf\{A\}^\{\(1\)\}and𝐀\(2\)\\mathbf\{A\}^\{\(2\)\}refer to the inverted edge maps\.

Let𝐘\(1\),𝐘\(2\)∈ℝp2/2\\mathbf\{Y\}^\{\(1\)\},\\mathbf\{Y\}^\{\(2\)\}\\in\\mathbb\{R\}^\{p\_\{2\}/2\}denote the vectorized representations of𝐀\(1\)\\mathbf\{A\}^\{\(1\)\}and𝐀\(2\)\\mathbf\{A\}^\{\(2\)\}, respectively, obtained by concatenating all image pixels into vectors, wherep2/2=H×Wp\_\{2\}/2=H\\times W\. We define the Bivariate image response as

𝐘=\[𝐘\(1\)𝐘\(2\)\]∈ℝp2\.\\mathbf\{Y\}=\\begin\{bmatrix\}\\mathbf\{Y\}^\{\(1\)\}\\\\ \\mathbf\{Y\}^\{\(2\)\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{p\_\{2\}\}\.Let𝐗=\(X1,X2,…,Xp1\)∈ℝp1\\mathbf\{X\}=\(X\_\{1\},X\_\{2\},\\ldots,X\_\{p\_\{1\}\}\)\\in\\mathbb\{R\}^\{p\_\{1\}\}denote the high\-dimensional predictor vector\.

Classical canonical correlation analysis \(CCA\) is a widely used statistical method for modeling the association between two sets of vectors\[[10](https://arxiv.org/html/2607.22763#bib.bib65)\]\. CCA seeks the canonical weight vectors𝝎x\\boldsymbol\{\\omega\}\_\{x\}and𝝎y\\boldsymbol\{\\omega\}\_\{y\}that maximize the correlation between the linear combinations𝝎x⊤​𝐗\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\mathbf\{X\}and𝝎y⊤​𝐘\\boldsymbol\{\\omega\}\_\{y\}^\{\\top\}\\mathbf\{Y\}by solving

max𝝎x,𝝎y⁡𝝎x⊤​𝚺x​y​𝝎ys\.t\.𝝎x⊤​𝚺x​𝝎x=1,𝝎y⊤​𝚺y​𝝎y=1,\\max\_\{\\boldsymbol\{\\omega\}\_\{x\},\\boldsymbol\{\\omega\}\_\{y\}\}\\ \\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\boldsymbol\{\\Sigma\}\_\{xy\}\\boldsymbol\{\\omega\}\_\{y\}\\quad\\text\{s\.t\.\}\\quad\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\boldsymbol\{\\Sigma\}\_\{x\}\\boldsymbol\{\\omega\}\_\{x\}=1,\\qquad\\boldsymbol\{\\omega\}\_\{y\}^\{\\top\}\\boldsymbol\{\\Sigma\}\_\{y\}\\boldsymbol\{\\omega\}\_\{y\}=1,\(1\)where𝚺x=cov⁡\(𝐗\)\\boldsymbol\{\\Sigma\}\_\{x\}=\\operatorname\{cov\}\(\\mathbf\{X\}\),𝚺y=cov⁡\(𝐘\)\\boldsymbol\{\\Sigma\}\_\{y\}=\\operatorname\{cov\}\(\\mathbf\{Y\}\), and𝚺x​y=cov⁡\(𝐗,𝐘\)\\boldsymbol\{\\Sigma\}\_\{xy\}=\\operatorname\{cov\}\(\\mathbf\{X\},\\mathbf\{Y\}\)\. Classical CCA is not directly applicable to high\-dimensional settings because the sample covariance matrices become singular or ill\-conditioned, resulting in unstable or non\-unique solutions\. To overcome this limitation, numerous sparse CCA methods have been proposed by incorporating anℓ1\\ell\_\{1\}penalty into the canonical weight vectors\[[28](https://arxiv.org/html/2607.22763#bib.bib66)\]\. However, these methods are primarily developed for continuous data and do not explicitly account for the highly sparse and zero\-inflated characteristics of edge\-map images\. Therefore, we employ the Semiparametric Sparse CCA, which is specifically designed to accommodate high\-dimensional zero\-inflated sparse data\[[31](https://arxiv.org/html/2607.22763#bib.bib79)\]\. Rather than relying on the conventional sample covariance matrix, SSCCA estimates the correlation structure under a truncated latent Gaussian copula model\.

Definition\[[31](https://arxiv.org/html/2607.22763#bib.bib79)\]\(Truncated Latent Gaussian Copula model\)\.The image \(or Bivariate image\) response𝐘\\mathbf\{Y\}is assumed to follow a truncated latent Gaussian copula model if there exists a latent random vector

𝐔=\(U1,…,Up2\)⊤∼NPN​\(𝟎,𝚺,f\),\\mathbf\{U\}=\(U\_\{1\},\\ldots,U\_\{p\_\{2\}\}\)^\{\\top\}\\sim\\mathrm\{NPN\}\(\\mathbf\{0\},\\boldsymbol\{\\Sigma\},f\),such thatYj=I​\(Uj\>Cj\)​Uj,j=1,…,p2,Y\_\{j\}=I\(U\_\{j\}\>C\_\{j\}\)U\_\{j\},~j=1,\\ldots,p\_\{2\},whereI​\(⋅\)I\(\\cdot\)denotes the indicator function,𝐂=\(C1,…,Cp2\)⊤\\mathbf\{C\}=\(C\_\{1\},\\ldots,C\_\{p\_\{2\}\}\)^\{\\top\}is a vector of truncation constants, and𝚺\\boldsymbol\{\\Sigma\}is the latent correlation matrix\. Here,NPN\\mathrm\{NPN\}denotes the nonparanormal distribution, andf=\{fj,j=1,…,p2\}f=\\\{f\_\{j\},~j=1,\\ldots,p\_\{2\}\\\}represents a collection of monotonically increasing transformation functions\. Consequently,𝐘∼TLNPN​\(𝟎,𝚺,f,𝐂\),\\mathbf\{Y\}\\sim\\mathrm\{TLNPN\}\(\\mathbf\{0\},\\boldsymbol\{\\Sigma\},f,\\mathbf\{C\}\),whereTLNPN\\mathrm\{TLNPN\}denotes the truncated latent nonparanormal distribution\.

To estimate the latent correlation matrix between the zero\-inflated response𝐘\\mathbf\{Y\}and the continuous predictor𝐗\\mathbf\{X\}, SSCCA employs a bridge function that relates Kendall’sτ\\tauto the latent correlation coefficient\. Specifically, the latent correlation betweenYjY\_\{j\}andXkX\_\{k\}is given byΣj​k=F−1​\(τj​k\),\\Sigma\_\{jk\}=F^\{\-1\}\(\\tau\_\{jk\}\),whereF​\(⋅\)F\(\\cdot\)denotes the bridge function, whose explicit form for this mixed data type is given below\[[31](https://arxiv.org/html/2607.22763#bib.bib79)\],

F​\(Σj​k;Δj\)=−2​Φ2​\(−Δj,0;12\)\+4​Φ3​\(−Δj,0,0;Σ3′\),F\(\\Sigma\_\{jk\};\\Delta\_\{j\}\)=\-2\\Phi\_\{2\}\\left\(\-\\Delta\_\{j\},0;\\frac\{1\}\{\\sqrt\{2\}\}\\right\)\+4\\Phi\_\{3\}\(\-\\Delta\_\{j\},0,0;\\Sigma^\{\\prime\}\_\{3\}\),and

Σ3′=\(11/2Σj​k/21/21Σj​kΣj​k/2Σj​k1\)\.\\Sigma^\{\\prime\}\_\{3\}=\\begin\{pmatrix\}1&1/\\sqrt\{2\}&\\Sigma\_\{jk\}/\\sqrt\{2\}\\\\ 1/\\sqrt\{2\}&1&\\Sigma\_\{jk\}\\\\ \\Sigma\_\{jk\}/\\sqrt\{2\}&\\Sigma\_\{jk\}&1\\end\{pmatrix\}\.Here,Φd​\(⋯,Σd′\)\\Phi\_\{d\}\(\\cdots,\\Sigma^\{\\prime\}\_\{d\}\)denotes the cdf of the standarddd\-variate normal distribution with correlation matrixΣd′\\Sigma^\{\\prime\}\_\{d\}\. The latent thresholdΔj=fj​\(Cj\)\\Delta\_\{j\}=f\_\{j\}\(C\_\{j\}\)is estimated byΔ^j=Φ−1​\(1−π^j\)\\widehat\{\\Delta\}\_\{j\}=\\Phi^\{\-1\}\(1\-\\widehat\{\\pi\}\_\{j\}\), whereπ^j\\widehat\{\\pi\}\_\{j\}is the proportion of non\-zero observations inYjY\_\{j\}\.

Letnndenote the sample size and let\{\(Xi​k,Yi​j\)\}i=1n\\\{\(X\_\{ik\},Y\_\{ij\}\)\\\}\_\{i=1\}^\{n\}be independent and identically distributed observations forYjY\_\{j\}andXkX\_\{k\}\. The sample Kendall’sτ\\tauis defined as

τ^j​k=2n​\(n−1\)​∑1≤i<i′≤nsign⁡\(Yi​j−Yi′​j\)​sign⁡\(Xi​k−Xi′​k\)\.\\widehat\{\\tau\}\_\{jk\}=\\frac\{2\}\{n\(n\-1\)\}\\sum\_\{1\\leq i<i^\{\\prime\}\\leq n\}\\operatorname\{sign\}\(Y\_\{ij\}\-Y\_\{i^\{\\prime\}j\}\)\\operatorname\{sign\}\(X\_\{ik\}\-X\_\{i^\{\\prime\}k\}\)\.Its population counterpart satisfiesτj​k=𝔼​\(τ^j​k\)=F​\(Σj​k,Δj\)\.\\tau\_\{jk\}=\\mathbb\{E\}\(\\widehat\{\\tau\}\_\{jk\}\)=F\(\\Sigma\_\{jk\},\\Delta\_\{j\}\)\.Accordingly, SSCCA estimates the population latent correlationΣj​k\\Sigma\_\{jk\}byr^j​k\\widehat\{r\}\_\{jk\}, wherer^j​k=F−1​\(τ^j​k\)\\widehat\{r\}\_\{jk\}=F^\{\-1\}\(\\widehat\{\\tau\}\_\{jk\}\)is obtained by replacing the population Kendall’sτj​k\\tau\_\{jk\}with its sample counterpartτ^j​k\\widehat\{\\tau\}\_\{jk\}through the inverse bridge function\. The resulting rank\-based correlation estimator is𝐑^=\(r^j​k\)\.\\widehat\{\\mathbf\{R\}\}=\(\\widehat\{r\}\_\{jk\}\)\.

To ensure positive definiteness,𝐑^\\widehat\{\\mathbf\{R\}\}is regularized as

𝐑~=\(1−μ\)​𝐑^PSD\+μ​𝐈,\\widetilde\{\\mathbf\{R\}\}=\(1\-\\mu\)\\widehat\{\\mathbf\{R\}\}^\{\\,\\mathrm\{PSD\}\}\+\\mu\\mathbf\{I\},where𝐈\\mathbf\{I\}is the identity matrix,μ\\muis a small positive threshold \(set to 0\.01 in this article\), and𝐑^PSD\\widehat\{\\mathbf\{R\}\}^\{\\,\\mathrm\{PSD\}\}denotes the projection of𝐑^\\widehat\{\\mathbf\{R\}\}onto the cone of positive semidefinite matrices\[[8](https://arxiv.org/html/2607.22763#bib.bib80),[31](https://arxiv.org/html/2607.22763#bib.bib79)\]\.

Finally, the SSCCA optimization problem is formulated as

min𝝎x,𝝎y\{−𝝎x⊤​𝐑~x​y​𝝎y\+λ1‖𝝎x∥1\+λ2​‖𝝎y‖1\}s\.t\.𝝎x⊤​𝐑~x​𝝎x≤1,𝝎y⊤​𝐑~y​𝝎y≤1\.\\begin\{split\}\\min\_\{\\boldsymbol\{\\omega\}\_\{x\},\\boldsymbol\{\\omega\}\_\{y\}\}\\;&\\left\\\{\-\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{xy\}\\boldsymbol\{\\omega\}\_\{y\}\+\\lambda\_\{1\}\\\|\\boldsymbol\{\\omega\}\_\{x\}\\\|\_\{1\}\+\\lambda\_\{2\}\\\|\\boldsymbol\{\\omega\}\_\{y\}\\\|\_\{1\}\\right\\\}\\\\ \\text\{s\.t\.\}\\quad&\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{x\}\\boldsymbol\{\\omega\}\_\{x\}\\leq 1,\\qquad\\boldsymbol\{\\omega\}\_\{y\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{y\}\\boldsymbol\{\\omega\}\_\{y\}\\leq 1\.\\end\{split\}\(2\)
This constrained optimization problem is equivalent to solving

𝝎~x=arg⁡min𝝎x⁡\{12​𝝎x⊤​𝐑~x​𝝎x−𝝎x⊤​𝐑~x​y​𝝎y\+λ1​‖𝝎x‖1\}\.\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}=\\arg\\min\_\{\\boldsymbol\{\\omega\}\_\{x\}\}\\left\\\{\\frac\{1\}\{2\}\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{x\}\\boldsymbol\{\\omega\}\_\{x\}\-\\boldsymbol\{\\omega\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{xy\}\\boldsymbol\{\\omega\}\_\{y\}\+\\lambda\_\{1\}\\\|\\boldsymbol\{\\omega\}\_\{x\}\\\|\_\{1\}\\right\\\}\.
The solution is then normalized by setting

𝝎^x=\{𝟎,𝝎~x=𝟎,𝝎~x\(𝝎~x⊤​𝐑~x​𝝎~x\)1/2,𝝎~x≠𝟎\.\\widehat\{\\boldsymbol\{\\omega\}\}\_\{x\}=\\begin\{cases\}\\mathbf\{0\},&\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}=\\mathbf\{0\},\\\\\[5\.69054pt\] \\displaystyle\\frac\{\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\}\{\\left\(\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{x\}\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\\right\)^\{1/2\}\},&\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\\neq\\mathbf\{0\}\.\\end\{cases\}
Another attractive property of SSCCA is its ability to perform variable selection through anℓ1\\ell\_\{1\}\(LASSO\) penalty, which shrinks the coefficients of non\-informative𝐗\\mathbf\{X\}toward zero while retaining only the predictors most strongly associated with𝐘\\mathbf\{Y\}\. The tuning parametersλ1\\lambda\_\{1\}andλ2\\lambda\_\{2\}are selected using the Bayesian Information Criterion \(BIC\) proposed byYoonet al\.\[[31](https://arxiv.org/html/2607.22763#bib.bib79)\], defined as

BIC=g​\(𝝎~x\)\+d​f𝝎~x​log⁡nn,\\mathrm\{BIC\}=g\(\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\)\+df\_\{\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\}\\frac\{\\log n\}\{n\},\(3\)whereg​\(𝝎~x\)=𝝎~x⊤​𝐑~x​𝝎~x−2​𝝎~x⊤​𝐑~x​y​𝝎y\+𝝎y⊤​𝐑~y​𝝎y,g\(\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\)=\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{x\}\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\-2\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{xy\}\\boldsymbol\{\\omega\}\_\{y\}\+\\boldsymbol\{\\omega\}\_\{y\}^\{\\top\}\\widetilde\{\\mathbf\{R\}\}\_\{y\}\\boldsymbol\{\\omega\}\_\{y\},andd​f𝝎~xdf\_\{\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\}denotes the number of nonzero entries in𝝎~x\\widetilde\{\\boldsymbol\{\\omega\}\}\_\{x\}\. By the symmetry of the objective function,𝝎y\\boldsymbol\{\\omega\}\_\{y\}can be estimated analogously\. The SSCCA model is implemented using themixedCCARpackage\.

## 3Simulation

We design two simulation scenarios of increasing complexity to evaluate the accuracy of SSCCA to identify the true predictors associated with the simulated leaf vein image data\. The following four evaluation metrics are adopted\[[13](https://arxiv.org/html/2607.22763#bib.bib69)\]\.

- •𝒮\\mathcal\{S\}is defined as the average number of predictors selected by the model across 100 replications\.
- •𝒫i\\mathcal\{P\}\_\{i\}is defined as the proportion of each individual true predictor being selected by the model across 100 replications\.
- •𝒫s\\mathcal\{P\}\_\{s\}is defined as the proportion of all true predictors being simultaneously selected by the model across 100 replications\.
- •ℱ​𝒮​ℛ\\mathcal\{FSR\}is defined as the false selection rate that represents the proportion of noise predictors that are incorrectly selected by the model across 100 replications\.

### 3\.1Simulation 1

We setn=200n=200andp1=50p\_\{1\}=50\. The predictor vectors𝐗\\mathbf\{X\}are generated from𝐗∼𝒩​\(𝟎,𝚺x\)\\mathbf\{X\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\boldsymbol\{\\Sigma\}\_\{x\}\), where𝚺x=\(σk1​k2\)\\boldsymbol\{\\Sigma\}\_\{x\}=\(\\sigma\_\{k\_\{1\}k\_\{2\}\}\)is defined byσk1​k2=ρ\|k1−k2\|,ρ=0\.8\.\\sigma\_\{k\_\{1\}k\_\{2\}\}=\\rho^\{\|k\_\{1\}\-k\_\{2\}\|\},~\\rho=0\.8\.To mimic the categorical nature of genetic markers, each component of𝐗\\mathbf\{X\}is discretized into three categorical levels\. LetXkX\_\{k\}denote thekkth predictor variable\. We define the set of true predictors asS∗=\{X6,X13,X21,X34,X41\},S^\{\*\}=\\\{X\_\{6\},X\_\{13\},X\_\{21\},X\_\{34\},X\_\{41\}\\\},with corresponding coefficientsβ6=1,β13=−2,β21=−1,β34=1,β41=0\.5\.\\beta\_\{6\}=1,~\\beta\_\{13\}=\-2,~\\beta\_\{21\}=\-1,~\\beta\_\{34\}=1,~\\beta\_\{41\}=0\.5\.All remaining coefficients are set to zero, that is,βk=0,k∉S∗\.\\beta\_\{k\}=0,~k\\notin S^\{\*\}\.

To generate the leaf vein image response data, we select one edge map extracted by EDTER from the real dataset as a template image, whose pixel values are denoted byTs,tT\_\{s,t\},s=1,…,Hs=1,\\ldots,Handt=1,…,Wt=1,\\ldots,W\. For subjectii, the simulated image𝐀i=\(Ai,s,t\)\\mathbf\{A\}\_\{i\}=\(A\_\{i,s,t\}\)is generated by perturbing the template image with a predictor\-dependent signal and random noise according to

Ai,s,t=Ts,t\+γ​\(∑k=150βk​Xi​k\)\+ϵi,s,t,A\_\{i,s,t\}=T\_\{s,t\}\+\\gamma\\Big\(\\sum\_\{k=1\}^\{50\}\\beta\_\{k\}X\_\{ik\}\\Big\)\+\\epsilon\_\{i,s,t\},whereϵi,s,t\\epsilon\_\{i,s,t\}are i\.i\.d\.𝒩​\(0,1\)\\mathcal\{N\}\(0,1\),γ=10\\gamma=10\. Finally, the pixel intensities of𝐀i\\mathbf\{A\}\_\{i\}are truncated to the rangeAi,s,t∈\[0,255\]A\_\{i,s,t\}\\in\[0,255\]to satisfy the valid intensity range of grayscale images\.

Figure[6](https://arxiv.org/html/2607.22763#S3.F6)presents an example of an image generated under Simulation 1\. In this setting, the predictor\-dependent signal is added uniformly across all pixels, affecting both background regions and vein structures\. As shown in Figure[6](https://arxiv.org/html/2607.22763#S3.F6), the resulting images are heavily contaminated by noise, making the underlying minor veins difficult to distinguish from the background\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/sim_1_v2.jpg)Figure 6:Illustration of Simulation 1\. The left panel shows the edge map of the reference leaf, whereas the right panel presents a randomly selected simulated leaf generated from the reference edge map\.Table 1:Results of Simulation 1\.𝒮\\mathcal\{S\}denotes the average model size,𝒫i\\mathcal\{P\}\_\{i\}the individual success rate,𝒫s\\mathcal\{P\}\_\{s\}the simultaneous success rate,ℱ​𝒮​ℛ\\mathcal\{FSR\}the false selection rate, and SD the sample standard deviation computed from 100 simulation replicates\.The results of Simulation 1 are summarized in Table[1](https://arxiv.org/html/2607.22763#S3.T1)\. The average model size,𝒮\\mathcal\{S\}, selected by SSCCA is 5\.72, which is very close to the true number of active predictors, indicating high variable selection accuracy\. The individual success rate,𝒫i\\mathcal\{P\}\_\{i\}, for each of the five true predictors is remarkably high\. In particular,X6X\_\{6\},X13X\_\{13\},X34X\_\{34\}, andX41X\_\{41\}are selected 100% with SD = 0\. The simultaneous success rate,𝒫s\\mathcal\{P\}\_\{s\}, reaches 99%, indicating near\-perfect discovery of all true signals simultaneously\. Meanwhile, the false selection rate,ℱ​𝒮​ℛ\\mathcal\{FSR\}, remains low at 0\.016, demonstrating effective control of false discoveries\.

### 3\.2Simulation 2

To increase the complexity level, Simulation 2 restricts the predictor effects to vein regions only\. Since vein pixels typically have higher intensity values than background pixels, we introduce a thresholdα\\alphaand allow the five true predictors to influence only those pixels with intensities exceedingα\\alpha\(in this article, we setα=80\\alpha=80\)\. Specifically, the response images are generated according to

Ai,s,t=Ts,t\+γ​\(∑k=150βk​Xi​k\)​I​\(Ts,t\>α\)\+ϵi,s,t,A\_\{i,s,t\}=T\_\{s,t\}\+\\gamma\\Big\(\\sum\_\{k=1\}^\{50\}\\beta\_\{k\}X\_\{ik\}\\Big\)I\(T\_\{s,t\}\>\\alpha\)\+\\epsilon\_\{i,s,t\},\(4\)whereI​\(⋅\)I\(\\cdot\)denotes the indicator function,ϵi,s,t\\epsilon\_\{i,s,t\}are i\.i\.d\.𝒩​\(0,1\)\\mathcal\{N\}\(0,1\), and the resulting pixel intensities are truncated to the interval\[0,255\]\[0,255\]\. All remaining parameters and coefficients are set to be identical to those in Simulation 1\.

The left panel of Figure[7](https://arxiv.org/html/2607.22763#S3.F7)displays the original template leaf edge mapTs,tT\_\{s,t\}, whereas the right panel shows the added signal \(i\.e\., the second and third term in equation \([4](https://arxiv.org/html/2607.22763#S3.E4)\)\), which enables us to visualize the regions where the predictor\-related signals and random noise are introduced, as well as their relative intensities\.

Table[2](https://arxiv.org/html/2607.22763#S3.T2)summarizes the results of Simulation 2\. The average model size,𝒮=5\.55\\mathcal\{S\}=5\.55, remains very close to the true number of active predictors, indicating accurate variable selection performance\. Compared with Simulation 1, the association signal is restricted to a small subset of pixels located on or near the vein structures, making signal discovery substantially more challenging\. But SSCCA still achieves a91%91\\%simultaneous success rate while maintaining a low false selection rate ofℱ​𝒮​ℛ=0\.021\\mathcal\{FSR\}=0\.021\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/sim_2_v2.jpg)Figure 7:Illustration of Simulation 2\. The left panel shows the edge map of the reference leaf, whereas the right panel illustrates the spatial locations and intensities of the simulated predictor signal and random noise added to the reference edge map\.Table 2:Results of Simulation 2\.𝒮\\mathcal\{S\}denotes the average model size,𝒫i\\mathcal\{P\}\_\{i\}the individual success rate,𝒫s\\mathcal\{P\}\_\{s\}the simultaneous success rate,ℱ​𝒮​ℛ\\mathcal\{FSR\}the false selection rate, and SD the sample standard deviation computed from 100 simulation replicates\.

## 4Real Data Analysis

The real data analyzed in this study were collected from a population of*Populus szechuanica*var\.*tibetica*, a member of the*Tacamahaca*section of the genus*Populus*\. Native to the Tibetan Plateau, this species is distributed across Gansu, Shaanxi, Sichuan, Xizang, and Yunnan provinces of China\. The sampled trees span a broad geographical range, with latitudes from28∘​58\.337′​N28^\{\\circ\}58\.337^\{\\prime\}\\mathrm\{N\}to29∘​58\.915′​N29^\{\\circ\}58\.915^\{\\prime\}\\mathrm\{N\}, longitudes from88∘​52\.283′​E88^\{\\circ\}52\.283^\{\\prime\}\\mathrm\{E\}to94∘​47\.903′​E94^\{\\circ\}47\.903^\{\\prime\}\\mathrm\{E\}, and elevations ranging from 2540 to 4081 meters\[[9](https://arxiv.org/html/2607.22763#bib.bib42)\]\. The leaves of*P\. szechuanica*var\.*tibetica*exhibit substantial variation in both morphology and vascular architecture, making them an ideal system for investigating the genetic and environmental factors regulating leaf vascular network\.

For each of thePopulustrees under study, one leaf was randomly selected and photographed from both the front and back sides\. To characterize the vascular whole\-network of each leaf, we apply the fine\-tuned EDTER model to the RGB images and generate corresponding edge maps\. Figure[8](https://arxiv.org/html/2607.22763#S4.F8)presents two representative leaf samples together with their edge maps\. The extracted edge maps preserve the whole\-network\. In some cases, the edge maps generated from the back side capture substantially more vein details than those obtained from the front side, justifying the inclusion of both sides in the analysis\.

To reduce computational cost, each Bivariate image is resized from its original dimension of900×600×2900\\times 600\\times 2to60×40×260\\times 40\\times 2\. The predictor set contains 118 standardized predictors, including 28 categorical genetic markers, 3 continuous geographic variables \(latitude, longitude, and elevation\), 84 gene–geography interaction terms, and 3 geography–geography interaction terms\. We then apply SSCCA to investigate the associations between the Bivariate edge\-map response images and the 118 predictor variables\.

![Refer to caption](https://arxiv.org/html/2607.22763v1/real3.jpg)Figure 8:Two randomly selected leaves from the real*Populus*dataset\. Each column corresponds to one leaf, with the top and bottom panels showing the front and back sides, respectively\. For each side, the original RGB image and the corresponding edge map generated by the fine\-tuned EDTER model are displayed side by side\.Among the 31 main effects and 87 interaction terms, SSCCA selects three variables as nonzero canonical coefficients, all of which correspond to gene–geography interaction effects:*GCPM\_1053\-1*×\\timeslongitude,*GCPM\_1036\-1*×\\timeselevation, and*GCPM\_1131*×\\timeslatitude\. These results suggest that gene–environment interactions play a very important role in influencing the vascular whole\-network of*Populus*leaves\.

In addition to identifying new scientific findings, our results also reconfirm some previous studies that were based on low\-dimensional vein traits\. For example, several leaf venation traits were reported to vary with elevation in woody plants, with trees at higher elevations generally exhibiting lower vein density but greater vein thickness and volume\[[26](https://arxiv.org/html/2607.22763#bib.bib90)\]\. Similarly,Zhuet al\.\[[33](https://arxiv.org/html/2607.22763#bib.bib88)\]found that minor vein density in natural populations of oriental oak decreased significantly with increasing latitude\. A study of three herbaceous species further reported that leaf major vein thickness increased with elevation\[[14](https://arxiv.org/html/2607.22763#bib.bib91)\]\.

## 5Conclusion

The proposed framework provides four methodological advances over existing approaches for gene–environment association studies of leaf vascular architecture\. By leveraging state\-of\-the\-art deep learning techniques, the proposed framework preserves the complete whole\-network of the leaf vascular architecture that cannot be adequately captured by a limited set of summary traits\. In addition, the transformer\-based EDTER model reduces the dependence on specialized image preparation and demonstrates robust performance in extracting vein structures from images with weak contrast\. From a statistical perspective, the SSCCA model enables the joint analysis of repeatedly measured high\-dimensional Bivariate image responses and high\-dimensional predictors, while simultaneously performing variable selection and accommodating the sparse and zero\-inflated nature of edge\-map images\.

The limitation of the real dataset analyzed in this study is its relatively small sample size and low marker resolution, which may limit the ability to identify specific gene names\. Nevertheless, as a methodology\-oriented study, our primary objective is to develop and evaluate a new framework that is readily applicable to larger datasets with high\-dimensional image responses and predictors\. Beyond leaf venation, the proposed framework can be extended to a broad range of network\-structured and image\-based phenotypes, including root system architecture, human vascular networks, neuroimaging, and transportation networks\. It is particularly well suited for applications where manually annotated ground\-truth labels are scarce or unavailable\.

Although SSCCA is designed for high\-dimensional predictors and responses, its computational efficiency may deteriorate when applied to ultra\-high\-dimensional image data due to memory limitations, computational burden, and increased processing time\. To address this challenge, we reduce the image dimension through resizing prior to the downstream association analysis\. We experimentally evaluate several image resolutions, including180×120180\\times 120,120×80120\\times 80,90×6090\\times 60, and60×4060\\times 40, through simulation studies\. The results indicate that the60×4060\\times 40resolution achieves the most favorable balance between computational efficiency and statistical accuracy\. Consequently, all analyses presented in this article, including both simulation studies and the real data application, are conducted using images resized to60×4060\\times 40\. Compared with alternative dimension\-reduction approaches such as principal component analysis \(PCA\), image resizing preserves the complete image information, facilitating biological interpretation\. While pixel level details may be smoothed during resizing, the overall venation architecture remains well preserved, making resizing a practical and interpretable strategy for handling ultrahigh\-dimensional image responses\.

## 6Competing interests

The authors declare that they have no competing interests\.

## References

- \[1\]P\. Arbelaez, M\. Maire, C\. Fowlkes, and J\. Malik\(2010\)Contour detection and hierarchical image segmentation\.IEEE transactions on pattern analysis and machine intelligence33\(5\),pp\. 898–916\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p8.1),[§2\.2](https://arxiv.org/html/2607.22763#S2.SS2.p2.1)\.
- \[2\]B\. Blonder, C\. Violle, L\. P\. Bentley, and B\. J\. Enquist\(2011\)Venation networks and the origin of the leaf economics spectrum\.Ecology letters14\(2\),pp\. 91–100\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p2.1)\.
- \[3\]T\. J\. Brodribb, T\. S\. Feild, and L\. Sack\(2010\)Viewing leaf structure and evolution from a hydraulic perspective\.Functional Plant Biology37\(6\),pp\. 488–498\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1)\.
- \[4\]J\. Bühler, L\. Rishmawi, D\. Pflugfelder, G\. Huber, H\. Scharr, M\. Hülskamp, M\. Koornneef, U\. Schurr, and S\. Jahnke\(2015\)PhenoVein—a tool for leaf vein segmentation and analysis\.Plant physiology169\(4\),pp\. 2359–2370\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p3.1)\.
- \[5\]S\. Dhondt, D\. Van Haerenborgh, C\. Van Cauwenbergh, R\. M\. Merks, W\. Philips, G\. T\. Beemster, and D\. Inzé\(2012\)Quantitative analysis of venation patterns of arabidopsis leaves by supervised image analysis\.The Plant Journal69\(3\),pp\. 553–563\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p3.1)\.
- \[6\]M\. Dirnberger, T\. Kehl, and A\. Neumann\(2015\)NEFI: network extraction from images\.Scientific reports5\(1\),pp\. 15669\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p3.1)\.
- \[7\]A\. Dosovitskiy\(2020\)An image is worth 16x16 words: transformers for image recognition at scale\.arXiv preprint arXiv:2010\.11929\.Cited by:[§2\.1\.1](https://arxiv.org/html/2607.22763#S2.SS1.SSS1.p1.6)\.
- \[8\]J\. Fan, H\. Liu, Y\. Ning, and H\. Zou\(2017\)High dimensional semiparametric latent graphical model for mixed data\.Journal of the Royal Statistical Society Series B: Statistical Methodology79\(2\),pp\. 405–421\.Cited by:[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p7.5)\.
- \[9\]G\. Fu, W\. Bo, X\. Pang, Z\. Wang, L\. Chen, Y\. Song, Z\. Zhang, J\. Li, and R\. Wu\(2013\)Mapping shape quantitative trait loci using a radius\-centroid\-contour model\.Heredity110\(6\),pp\. 511–519\.Cited by:[§4](https://arxiv.org/html/2607.22763#S4.p1.4)\.
- \[10\]H\. Hotelling\(1992\)Relations between two sets of variates\.InBreakthroughs in statistics: methodology and distribution,pp\. 162–190\.Cited by:[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p3.4)\.
- \[11\]M\. Huang, X\. Liu, Y\. Zhou, R\. M\. Summers, and Z\. Zhang\(2019\)BLINK: a package for the next level of genome\-wide association studies with both individuals and markers in the millions\.Gigascience8\(2\),pp\. giy154\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p5.1)\.
- \[12\]J\. Lagergren, M\. Pavicic, H\. B\. Chhetri, L\. M\. York, D\. Hyatt, D\. Kainer, E\. M\. Rutter, K\. Flores, J\. Bailey\-Bale, M\. Klein,et al\.\(2023\)Few\-shot learning enables population\-scale analysis of leaf traits in populus trichocarpa\.Plant Phenomics5,pp\. 0072\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p4.1)\.
- \[13\]R\. Li, W\. Zhong, and L\. Zhu\(2012\)Feature screening via distance correlation learning\.Journal of the American Statistical Association107\(499\),pp\. 1129–1139\.Cited by:[§3](https://arxiv.org/html/2607.22763#S3.p1.1)\.
- \[14\]W\. Liu, L\. Zheng, and D\. Qi\(2020\)Variation in leaf traits at different altitudes reflects the adaptive strategy of plants to environmental changes\. ecol evol 10: 8166–8175\.Cited by:[§4](https://arxiv.org/html/2607.22763#S4.p5.1)\.
- \[15\]X\. Liu, M\. Huang, B\. Fan, E\. S\. Buckler, and Z\. Zhang\(2016\)Iterative usage of fixed and random effect models for powerful and efficient genome\-wide association studies\.PLoS genetics12\(2\),pp\. e1005767\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p5.1)\.
- \[16\]Y\. Liu, M\. Cheng, X\. Hu, K\. Wang, and X\. Bai\(2017\)Richer convolutional features for edge detection\.InProceedings of the IEEE conference on computer vision and pattern recognition,pp\. 3000–3009\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p7.1)\.
- \[17\]P\. Lucas, M\. Choong, H\. Tan, I\. Turner, and A\. Berrick\(1991\)The fracture toughness of the leaf of the dicotyledon calophyllum inophyllum l\.\(guttiferae\)\.Philosophical Transactions of the Royal Society of London\. Series B: Biological Sciences334\(1269\),pp\. 95–106\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1)\.
- \[18\]S\. Narawatthana, Y\. Phansenee, B\. Thammasamisorn, and P\. Vejchasarn\(2023\)Multi\-model genome\-wide association studies of leaf anatomical traits and vein architecture in rice\.Frontiers in Plant Science14,pp\. 1107718\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1),[§1](https://arxiv.org/html/2607.22763#S1.p2.1),[§1](https://arxiv.org/html/2607.22763#S1.p5.1)\.
- \[19\]A\. Nardini\(2022\)Hard and tough: the coordination between leaf mechanical resistance and drought tolerance\.Flora288,pp\. 152023\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1)\.
- \[20\]Y\. Onoda, M\. Westoby, P\. B\. Adler, A\. M\. Choong, F\. J\. Clissold, J\. H\. Cornelissen, S\. Díaz, N\. J\. Dominy, A\. Elgart, L\. Enrico,et al\.\(2011\)Global patterns of leaf mechanical properties\.Ecology letters14\(3\),pp\. 301–312\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p2.1)\.
- \[21\]C\. A\. Price, O\. Symonova, Y\. Mileyko, T\. Hilley, and J\. S\. Weitz\(2011\)Leaf extraction and analysis framework graphical user interface: segmenting and analyzing the structure of leaf veins and areoles\.Plant Physiology155\(1\),pp\. 236–245\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p2.1),[§1](https://arxiv.org/html/2607.22763#S1.p3.1)\.
- \[22\]M\. Pu, Y\. Huang, Y\. Liu, Q\. Guan, and H\. Ling\(2022\)Edter: edge detection with transformer\.InProceedings of the IEEE/CVF conference on computer vision and pattern recognition,pp\. 1402–1412\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p7.1),[§2\.1](https://arxiv.org/html/2607.22763#S2.SS1.p1.1)\.
- \[23\]L\. Rishmawi, J\. Bühler, B\. Jaegle, M\. Hülskamp, and M\. Koornneef\(2017\)Quantitative trait loci controlling leaf venation in arabidopsis\.Plant, cell & environment40\(8\),pp\. 1429–1441\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1),[§1](https://arxiv.org/html/2607.22763#S1.p2.1),[§1](https://arxiv.org/html/2607.22763#S1.p5.1)\.
- \[24\]L\. Sack and C\. Scoffoni\(2013\)Leaf venation: structure, function, development, evolution, ecology and applications in the past, present and future\.New phytologist198\(4\),pp\. 983–1000\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p1.1),[§1](https://arxiv.org/html/2607.22763#S1.p2.1)\.
- \[25\]S\. C\. Siddharth, U\. Singh, A\. Kaul, and S\. Jain\(2019\)A database of leaf images: practice towards plant conservation with plant pathology\.Mendeley Data\.Cited by:[§2\.2](https://arxiv.org/html/2607.22763#S2.SS2.p1.1)\.
- \[26\]R\. Wang, H\. Chen, X\. Liu, Z\. Wang, J\. Wen, and S\. Zhang\(2020\)Plant phylogeny and growth form as drivers of the altitudinal variation in woody leaf vein traits\.Frontiers in Plant Science10,pp\. 1735\.Cited by:[§4](https://arxiv.org/html/2607.22763#S4.p5.1)\.
- \[27\]X\. Wang, K\. Yu, C\. Dong, and C\. C\. Loy\(2018\)Recovering realistic texture in image super\-resolution by deep spatial feature transform\.InProceedings of the IEEE conference on computer vision and pattern recognition,pp\. 606–615\.Cited by:[§2\.1\.3](https://arxiv.org/html/2607.22763#S2.SS1.SSS3.p1.3)\.
- \[28\]D\. M\. Witten, R\. Tibshirani, and T\. Hastie\(2009\)A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis\.Biostatistics10\(3\),pp\. 515–534\.Cited by:[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p3.8)\.
- \[29\]S\. Xie and Z\. Tu\(2015\)Holistically\-nested edge detection\.InProceedings of the IEEE international conference on computer vision,pp\. 1395–1403\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p7.1)\.
- \[30\]Y\. Ye, K\. Xu, Y\. Huang, R\. Yi, and Z\. Cai\(2024\)Diffusionedge: diffusion probabilistic model for crisp edge detection\.InProceedings of the AAAI conference on artificial intelligence,Vol\.38,pp\. 6675–6683\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p8.1)\.
- \[31\]G\. Yoon, R\. J\. Carroll, and I\. Gaynanova\(2020\)Sparse semiparametric canonical correlation analysis for data of mixed types\.Biometrika107\(3\),pp\. 609–625\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p9.1),[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p11.5),[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p3.8),[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p4.1.1),[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p5.7),[§2\.3](https://arxiv.org/html/2607.22763#S2.SS3.p7.5)\.
- \[32\]J\. Yu, G\. Pressoir, W\. H\. Briggs, I\. Vroh Bi, M\. Yamasaki, J\. F\. Doebley, M\. D\. McMullen, B\. S\. Gaut, D\. M\. Nielsen, J\. B\. Holland,et al\.\(2006\)A unified mixed\-model method for association mapping that accounts for multiple levels of relatedness\.Nature genetics38\(2\),pp\. 203–208\.Cited by:[§1](https://arxiv.org/html/2607.22763#S1.p5.1)\.
- \[33\]Y\. Zhu, H\. Kang, Q\. Xie, Z\. Wang, S\. Yin, and C\. Liu\(2012\)Pattern of leaf vein density and climate relationship of quercus variabilis populations remains unchanged with environmental changes\.Trees26\(2\),pp\. 597–607\.Cited by:[§4](https://arxiv.org/html/2607.22763#S4.p5.1)\.