diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..7015d2de83655bf366d5115e417de29fbca8abbd
Binary files /dev/null and b/.DS_Store differ
diff --git a/manuscript/.DS_Store b/manuscript/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..6b7cae23c24b2371aa37887bf242ff81483ab255
Binary files /dev/null and b/manuscript/.DS_Store differ
diff --git a/manuscript/RSF.pdf b/manuscript/RSF.pdf
index 2893faf7b86a8006747abd276771f45114bbb2e2..16fe34c7fd817feb4ebd41cbb19085ba6c572887 100644
Binary files a/manuscript/RSF.pdf and b/manuscript/RSF.pdf differ
diff --git a/manuscript/RSF.tex b/manuscript/RSF.tex
index 0b5e4405a7a0981903993fdb2f3d84caf66ac15a..2f826004fe41618fc0755a9df3031f4f8a3a8b85 100644
--- a/manuscript/RSF.tex
+++ b/manuscript/RSF.tex
@@ -91,14 +91,16 @@ tl;dr RKHS make subsequent optimization learning easier to implement and more li
\subsection{Optical Coherence Tomography of the Retina}
-[Ronald?] -- ~ 4 paragraphs summarize impact of retinal OCT (Glaucoma, occulomics, etc.), background on VFMD, background on metadata, imaging physics.
-Key recent deep learning papers [Schuman et al.]. Other learning approaches [NNMF et al.].
--------
+The last few decades have ushered in revolutionary advancements in medical imaging which have helped elucidate relationships between biological structures and physiological function in health and disease. Combining computer-aided diagnosis with imaging data further augments physicians’ ability to make higher precision clinical decisions. Perhaps no better example of this could be made than with the impact that optical coherence tomography (OCT) has had on eye care. OCT is a non-invasive diagnostic imaging tool which employs principles of optical interferometry using a low-coherence light source to obtain cross-sectional digitally reconstructed images of biological tissue \cite{huang1991}. Often described as an optical analog to ultrasound, OCT images are generated from detection of electrical field produced from the echoed tissue and reference signals \cite{fujimoto2016development}. OCT datasets are comprised of multiple A-scans acquired in rapid succession to form cross-sectional images called B-Scans. These B-scans can be rapidly acuired in a raster pattern to obtain a volumetric image of the subject tissue. Current generations of spectral domain (SD) OCT have an axial resolution of 5-8 microns with a lateral resolution of 6-20 microns and a scanning rate of 100,000 A-scans per second, making them ideal for Capture of 3-D volumetric data in vivo. \cite{alexopoulos2022development}.
-Optical coherence tomography (OCT) is a non-invasive diagnostic imaging tool which employs principles of optical interferometry to obtain cross-sectional images from biological tissue \cite[]{huang1991}.
+Although predominantly used in clinical eye practice, OCT spans beyond ophthalmology. OCT has clinical applications in cardiology, otology, dermatology, and dentistry \cite{ali2021optical}. OCT’s ability to penetrate up to 3mm of tissue and produce high resolution images have provided use cases in diagnosing myocardial infarction with non-obtrusive coronary arteries \cite{reynolds2021coronary}, facilitating cochlear implants surgery \cite{starovoyt2019high}, diagnosing basal cell carcinoma \cite{chen2021evaluation}, and even dental cavity detection \cite{hariri2013estimation}. More recently, the sub-field of ophthalmic imaging known as Oculomics has shown the capability of using OCT to make prognostications of non-ocular systemic disorders like cardiovascular \cite{chan2023eyes} and neurological disease \cite{suh2023retina, lin2024individual}, as well as improved estimation of phenotypic age for predicting mortality \cite{nusinovici2022retinal}.
--------
+Mainly, this ability to safely and quickly acquire ocular images in a longitudinal and reproducible manner has revolutionized the way clinicians diagnose and manage blindness causing diseases such as glaucoma. Glaucoma is the global leader of irreversible blindness, with some projections indicating that the number of affected individuals could rise to 112 million by 2040 due to population aging \cite{tham2014global}. Glaucoma, characterized by persistent microstructural deterioriation of the optic nerve and retina, is often asymptomatic until the moderate to severe stages of disease. Therefore, OCT plays a critical role in mitigating disease through early detection of these microscopic structural changes. However, OCT alone is insufficient in determining disease severity. Glaucoma standard of care involves ophthalmic imaging with OCT in conjunction with functional assessment of the patient’s visual field (VF) through standardized automated perimetry \cite{chen2024deep}. One of the prominent ways to categorize disease severity and visual field defect is through visual field mean deviation (VFMD), which is a numerical estimation of light in decibels (dB) that an individual eye can perceive \cite{thirunavukarasu2024validated}. It is important to note, however, that structural changes in the retina are not linearly correlated with VFMD throughout the disease spectrum. Cross-sectional studies have demonstrated that there is little to no correlation between mean RNFL and VFMD until significant RNFL damage has occurred, at which point statistically significant associations with VFMD can be observed \cite{wollstein2012retinal}. Therefore, improving early detection of structural changes through ocular imaging are mission critical in the pursuit of vision preservation.
+
+Artificial intelligence (AI) applications using deep learners (DL) to detect, diagnose and predict disease progression have a potentially enormous impact on public health \cite{SCHUMAN2022e3}. The most common approaches look to classify glaucoma based on clinically accepted structural features such as the retinal nerve fiber layer (RNFL), ganglion cell inner plexiform layer (GCIPL) and optic nerve head (ONH) \cite{prahs2018oct, grewal2008artificial, shi2024rnflt2vec, berenguer2021automatic, christopher2018retinal}. Furthermore, incorporating features such as ganglion cell complex (GCC) thickness, ONH macrostructure, and RNFL reflectance maps into DL model training have been shown to significantly improve diagnostic accuracy \cite{tan2024hybrid}. Although there is no consensus on the exact clinical classification of glaucoma, the structure function relationship between OCT and VF are critical in understanding disease severity and progression. This reliability on VF for diagnosis is challenging because the VF test itself is challenging for patients. The test, which can last between 4-8 minutes per eye, requires patients to sit still and fixate on singular point without moving and click a button anytime a point of light is perceived in their field of vision. Test difficulty, combined with an older patient demographic, often makes VF test results more variable than OCT, which takes mere seconds to perform. Therefore, considerable efforts have been made towards predicting VF outcomes based on OCT. Previous studies have shown the ability to derive the spatial relationship between structure and function using OCT and VF point-by-point estimation \cite{chen2024deep} based on stochastic optimization DL models \cite{kingma2014adam}.
+
+Unsupervised feature agnostic approaches using OCT have shown promise in predicting disease progression measured by VF. Chen et al. \cite{chen2023segmentation}, describe a 3D based ResNet18 Convolutional neural network (CNN) capable of inferring pointwise VF sensitivities directly from segmentation-free OCT 3D volumes. Similarly, attention-guided network approaches have demonstrated improved glaucoma detection by taking OCT volumes and computing dual 3-D gradient class activation heatmaps to predict the Visual Field Index (VFI) estimation \cite{george2020attention}. DL approaches using OCT volumes from glaucoma subjects have identified 14 non-clinically defined surface shape patterns near the ONH which are capable of predicting specific VFMD loss rate with a best coefficient of determination of $r^2=0.37$ \cite{saini2022assessing}. These clinically agnostic DL approaches using OCT volumes have demonstrated the ability to identify previously undiscovered biomarkers for improving the prediction of glaucomatous functional defects, but, perhaps due to the insufficient confidence indices or inability to explicitly and metrically define the relationships, remain limited to research applications.
Key observation -- inherent anisotropy of imaging implements Frangi-like plate filter along the high resolution axis.
diff --git a/manuscript/rsf.bib b/manuscript/rsf.bib
index 5043d120da2bef08ee4a4bedc7e3272233a0ebe0..8dff073455a60e066c74d2f14c5bb841e65db877 100644
--- a/manuscript/rsf.bib
+++ b/manuscript/rsf.bib
@@ -261,3 +261,296 @@
year={1991},
doi={10.1126/science.1957169}
}
+
+@article{fujimoto2016development,
+ title={The development, commercialization, and impact of optical coherence tomography},
+ author={Fujimoto, James and Swanson, Eric},
+ journal={Investigative ophthalmology \& visual science},
+ volume={57},
+ number={9},
+ pages={OCT1--OCT13},
+ year={2016},
+ publisher={The Association for Research in Vision and Ophthalmology}
+}
+
+@article{fujimoto2000optical,
+ title={Optical coherence tomography: an emerging technology for biomedical imaging and optical biopsy},
+ author={Fujimoto, James G and Pitris, Costas and Boppart, Stephen A and Brezinski, Mark E},
+ journal={Neoplasia},
+ volume={2},
+ number={1-2},
+ pages={9--25},
+ year={2000},
+ publisher={Elsevier}
+}
+
+@article{alexopoulos2022development,
+ title={The development and clinical application of innovative optical ophthalmic imaging techniques},
+ author={Alexopoulos, Palaiologos and Madu, Chisom and Wollstein, Gadi and Schuman, Joel S},
+ journal={Frontiers in medicine},
+ volume={9},
+ pages={891369},
+ year={2022},
+ publisher={Frontiers Media SA}
+}
+
+@article{ali2021optical,
+ title={Optical coherence tomography’s current clinical medical and dental applications: a review},
+ author={Ali, Saqib and Gilani, Saqlain Bin Syed and Shabbir, Juzer and Almulhim, Khalid S and Bugshan, Amr and Farooq, Imran},
+ journal={F1000Research},
+ volume={10},
+ year={2021},
+ publisher={Faculty of 1000 Ltd}
+}
+
+@article{reynolds2021coronary,
+ title={Coronary optical coherence tomography and cardiac magnetic resonance imaging to determine underlying causes of myocardial infarction with nonobstructive coronary arteries in women},
+ author={Reynolds, Harmony R and Maehara, Akiko and Kwong, Raymond Y and Sedlak, Tara and Saw, Jacqueline and Smilowitz, Nathaniel R and Mahmud, Ehtisham and Wei, Janet and Marzo, Kevin and Matsumura, Mitsuaki and others},
+ journal={Circulation},
+ volume={143},
+ number={7},
+ pages={624--640},
+ year={2021},
+ publisher={Am Heart Assoc}
+}
+
+@article{starovoyt2019high,
+ title={High-resolution imaging of the human cochlea through the round window by means of optical coherence tomography},
+ author={Starovoyt, Anastasiya and Putzeys, Tristan and Wouters, Jan and Verhaert, Nicolas},
+ journal={Scientific reports},
+ volume={9},
+ number={1},
+ pages={14271},
+ year={2019},
+ publisher={Nature Publishing Group UK London}
+}
+
+@article{chen2021evaluation,
+ title={Evaluation of ultrahigh-resolution optical coherence tomography for basal cell carcinoma, seborrheic keratosis, and nevus},
+ author={Chen, Shufen and Xie, Fang and Hao, Tian and Xie, Jun and Wang, Xianghong and Chen, Si and Liu, Linbo and Li, Chengxin},
+ journal={Skin Research and Technology},
+ volume={27},
+ number={4},
+ pages={479--485},
+ year={2021},
+ publisher={Wiley Online Library}
+}
+
+@article{hariri2013estimation,
+ title={Estimation of the enamel and dentin mineral content from the refractive index},
+ author={Hariri, Ilnaz and Sadr, Alireza and Nakashima, Syozi and Shimada, Yasushi and Tagami, Junji and Sumi, Yasunori},
+ journal={Caries research},
+ volume={47},
+ number={1},
+ pages={18--26},
+ year={2013},
+ publisher={S. Karger AG}
+}
+
+@article{chan2023eyes,
+ title={Eyes as the windows into cardiovascular disease in the era of big data},
+ author={Chan, Yarn Kit and Cheng, Ching-Yu and Sabanayagam, Charumathi},
+ journal={Taiwan Journal of Ophthalmology},
+ volume={13},
+ number={2},
+ pages={151--167},
+ year={2023},
+ publisher={Medknow}
+}
+
+@article{suh2023retina,
+ title={Retina oculomics in neurodegenerative disease},
+ author={Suh, Alex and Ong, Joshua and Kamran, Sharif Amit and Waisberg, Ethan and Paladugu, Phani and Zaman, Nasif and Sarker, Prithul and Tavakkoli, Alireza and Lee, Andrew G},
+ journal={Annals of Biomedical Engineering},
+ volume={51},
+ number={12},
+ pages={2708--2721},
+ year={2023},
+ publisher={Springer}
+}
+
+@article{lin2024individual,
+ title={Individual prognostication of disease activity and disability worsening in multiple sclerosis with retinal layer thickness z scores},
+ author={Lin, Ting-Yi and Motamedi, Seyedamirhosein and Asseyer, Susanna and Chien, Claudia and Saidha, Shiv and Calabresi, Peter A and Fitzgerald, Kathryn C and Samadzadeh, Sara and Villoslada, Pablo and Llufriu, Sara and others},
+ journal={Neurology: Neuroimmunology \& Neuroinflammation},
+ volume={11},
+ number={5},
+ pages={e200269},
+ year={2024},
+ publisher={Lippincott Williams \& Wilkins Hagerstown, MD}
+}
+
+@article{nusinovici2022retinal,
+ title={Retinal photograph-based deep learning predicts biological age, and stratifies morbidity and mortality risk},
+ author={Nusinovici, Simon and Rim, Tyler Hyungtaek and Yu, Marco and Lee, Geunyoung and Tham, Yih-Chung and Cheung, Ning and Chong, Crystal Chun Yuen and Da Soh, Zhi and Thakur, Sahil and Lee, Chan Joo and others},
+ journal={Age and ageing},
+ volume={51},
+ number={4},
+ pages={afac065},
+ year={2022},
+ publisher={Oxford University Press}
+}
+
+@article{tham2014global,
+ title={Global prevalence of glaucoma and projections of glaucoma burden through 2040: a systematic review and meta-analysis},
+ author={Tham, Yih-Chung and Li, Xiang and Wong, Tien Y and Quigley, Harry A and Aung, Tin and Cheng, Ching-Yu},
+ journal={Ophthalmology},
+ volume={121},
+ number={11},
+ pages={2081--2090},
+ year={2014},
+ publisher={Elsevier}
+}
+
+@article{chen2024deep,
+ title={Deep-Learning-Based Group Pointwise Spatial Mapping of Structure to Function in Glaucoma},
+ author={Chen, Zhiqi and Ishikawa, Hiroshi and Wang, Yao and Wollstein, Gadi and Schuman, Joel S},
+ journal={Ophthalmology Science},
+ volume={4},
+ number={5},
+ pages={100523},
+ year={2024},
+ publisher={Elsevier}
+}
+
+@article{thirunavukarasu2024validated,
+ title={A validated web-application (GFDC) for automatic classification of glaucomatous visual field defects using Hodapp-Parrish-Anderson criteria},
+ author={Thirunavukarasu, Arun James and Jain, Nikhil and Sanghera, Rohan and Lattuada, Federico and Mahmood, Shathar and Economou, Anna and Yu, Helmut CY and Bourne, Rupert},
+ journal={npj Digital Medicine},
+ volume={7},
+ number={1},
+ pages={131},
+ year={2024},
+ publisher={Nature Publishing Group UK London}
+}
+
+@article{wollstein2012retinal,
+ title={Retinal nerve fibre layer and visual function loss in glaucoma: the tipping point},
+ author={Wollstein, Gadi and Kagemann, Larry and Bilonick, Richard A and Ishikawa, Hiroshi and Folio, Lindsey S and Gabriele, Michelle L and Ungar, Allison K and Duker, Jay S and Fujimoto, James G and Schuman, Joel S},
+ journal={British Journal of Ophthalmology},
+ volume={96},
+ number={1},
+ pages={47--52},
+ year={2012},
+ publisher={BMJ Publishing Group Ltd}
+}
+
+@article{SCHUMAN2022e3,
+title = {A Case for the Use of Artificial Intelligence in Glaucoma Assessment},
+journal = {Ophthalmology Glaucoma},
+volume = {5},
+number = {3},
+pages = {e3-e13},
+year = {2022},
+issn = {2589-4196},
+doi = {https://doi.org/10.1016/j.ogla.2021.12.003},
+url = {https://www.sciencedirect.com/science/article/pii/S2589419621002805},
+author = {Joel S. Schuman and Maria {De Los Angeles Ramos Cadena} and Rebecca McGee and Lama A. Al-Aswad and Felipe A. Medeiros and Michael Abramoff and Mark Blumenkranz and Emily Chew and Michael Chiang and Malvina Eydelman and David Myung and Carol Shields and Bhavna J. Antony and Tin Aung and Michael Boland and Tom Brunner and Robert T. Chang and Balwantray Chauhan and D. Hunter Cherwek and David Garway-Heath and Adrienne Graves and Jeffrey L. Goldberg and Minguang He and Naama Hammel and Donald Hood and Hiroshi Ishikawa and Chris Leung and Louis Pasquale and Harry A. Quigley and Calvin W. Roberts and Alan L. Robin and Elena Sturman and Remo Susanna and Jayme Vianna and Linda Zangwill},
+keywords = {Artificial intelligence, Deep learning, Glaucoma, Neural networks, OCT},
+abstract = {We hypothesize that artificial intelligence (AI) applied to relevant clinical testing in glaucoma has the potential to enhance the ability to detect glaucoma. This premise was discussed at the recent Collaborative Community on Ophthalmic Imaging meeting, “The Future of Artificial Intelligence–Enabled Ophthalmic Image Interpretation: Accelerating Innovation and Implementation Pathways,” held virtually September 3–4, 2020. The Collaborative Community on Ophthalmic Imaging (CCOI) is an independent self-governing consortium of stakeholders with broad international representation from academic institutions, government agencies, and the private sector whose mission is to act as a forum for the purpose of helping speed innovation in healthcare technology. It was 1 of the first 2 such organizations officially designated by the Food and Drug Administration in September 2019 in response to their announcement of the collaborative community program as a strategic priority for 2018–2020. Further information on the CCOI can be found online at their website (https://www.cc-oi.org/about). Artificial intelligence for glaucoma diagnosis would have high utility globally, because access to care is limited in many parts of the world and half of all people with glaucoma are unaware of their illness. The application of AI technology to glaucoma diagnosis has the potential to broadly increase access to care worldwide, in essence flattening the Earth by providing expert-level evaluation to individuals even in the most remote regions of the planet.}
+}
+
+@article{prahs2018oct,
+ title={OCT-based deep learning algorithm for the evaluation of treatment indication with anti-vascular endothelial growth factor medications},
+ author={Prahs, Philipp and Radeck, Viola and Mayer, Christian and Cvetkov, Yordan and Cvetkova, Nadezhda and Helbig, Horst and M{\"a}rker, David},
+ journal={Graefe's Archive for Clinical and Experimental Ophthalmology},
+ volume={256},
+ pages={91--98},
+ year={2018},
+ publisher={Springer}
+}
+
+@article{grewal2008artificial,
+ title={Artificial neural network-based glaucoma diagnosis using retinal nerve fiber layer analysis},
+ author={Grewal, DS and Jain, R and Grewal, SPS and Rihani, V},
+ journal={European journal of ophthalmology},
+ volume={18},
+ number={6},
+ pages={915--921},
+ year={2008},
+ publisher={SAGE Publications Sage UK: London, England}
+}
+
+@article{shi2024rnflt2vec,
+ title={RNFLT2Vec: Artifact-corrected representation learning for retinal nerve fiber layer thickness maps},
+ author={Shi, Min and Tian, Yu and Luo, Yan and Elze, Tobias and Wang, Mengyu},
+ journal={Medical Image Analysis},
+ volume={94},
+ pages={103110},
+ year={2024},
+ publisher={Elsevier}
+}
+
+@article{berenguer2021automatic,
+ title={Automatic segmentation of the retinal nerve fiber layer by means of mathematical morphology and deformable models in 2d optical coherence tomography imaging},
+ author={Berenguer-Vidal, Rafael and Verd{\'u}-Monedero, Rafael and Morales-S{\'a}nchez, Juan and Sell{\'e}s-Navarro, Inmaculada and Del Amor, Roc{\'\i}o and Garc{\'\i}a, Gabriel and Naranjo, Valery},
+ journal={Sensors},
+ volume={21},
+ number={23},
+ pages={8027},
+ year={2021},
+ publisher={MDPI}
+}
+
+@article{christopher2018retinal,
+ title={Retinal nerve fiber layer features identified by unsupervised machine learning on optical coherence tomography scans predict glaucoma progression},
+ author={Christopher, Mark and Belghith, Akram and Weinreb, Robert N and Bowd, Christopher and Goldbaum, Michael H and Saunders, Luke J and Medeiros, Felipe A and Zangwill, Linda M},
+ journal={Investigative ophthalmology \& visual science},
+ volume={59},
+ number={7},
+ pages={2748--2756},
+ year={2018},
+ publisher={The Association for Research in Vision and Ophthalmology}
+}
+
+@article{tan2024hybrid,
+ title={A Hybrid Deep Learning Classification of Perimetric Glaucoma Using Peripapillary Nerve Fiber Layer Reflectance and Other OCT Parameters from Three Anatomy Regions},
+ author={Tan, Ou and Greenfield, David S and Francis, Brian A and Varma, Rohit and Schuman, Joel S and Huang, David and Choi, Dongseok},
+ journal={ArXiv},
+ year={2024},
+ publisher={arXiv}
+}
+
+@article{kingma2014adam,
+ title={Adam: A method for stochastic optimization},
+ author={Kingma, Diederik P and Ba, Jimmy},
+ journal={arXiv preprint arXiv:1412.6980},
+ year={2014}
+}
+
+@article{chen2023segmentation,
+ title={Segmentation-free OCT-volume-based deep learning model improves pointwise visual field sensitivity estimation},
+ author={Chen, Zhiqi and Shemuelian, Eitan and Wollstein, Gadi and Wang, Yao and Ishikawa, Hiroshi and Schuman, Joel S},
+ journal={Translational Vision Science \& Technology},
+ volume={12},
+ number={6},
+ pages={28--28},
+ year={2023},
+ publisher={The Association for Research in Vision and Ophthalmology}
+}
+
+@article{george2020attention,
+ title={Attention-guided 3D-CNN framework for glaucoma detection and structural-functional association using volumetric images},
+ author={George, Yasmeen and Antony, Bhavna J and Ishikawa, Hiroshi and Wollstein, Gadi and Schuman, Joel S and Garnavi, Rahil},
+ journal={IEEE journal of biomedical and health informatics},
+ volume={24},
+ number={12},
+ pages={3421--3430},
+ year={2020},
+ publisher={IEEE}
+}
+
+@article{saini2022assessing,
+ title={Assessing surface shapes of the optic nerve head and peripapillary retinal nerve fiber layer in glaucoma with artificial intelligence},
+ author={Saini, Chhavi and Shen, Lucy Q and Pasquale, Louis R and Boland, Michael V and Friedman, David S and Zebardast, Nazlee and Fazli, Mojtaba and Li, Yangjiani and Eslami, Mohammad and Elze, Tobias and others},
+ journal={Ophthalmology Science},
+ volume={2},
+ number={3},
+ pages={100161},
+ year={2022},
+ publisher={Elsevier}
+}
+
+
+
+
diff --git a/src/+RSF/getImf.m b/src/+RSF/getImf.m
deleted file mode 100644
index d2d7887c4a7ae2d5a28571ba873f3b3a3e528e02..0000000000000000000000000000000000000000
--- a/src/+RSF/getImf.m
+++ /dev/null
@@ -1,40 +0,0 @@
-function imf = getImf(pData,radii,scanType)
-
-NGPU = HIP.Cuda.DeviceCount;
-p = ljsStartParallel(8*NGPU);
-
-imf = {};
-px = {};
-parfor pp = 1:height(pData)
- strDB = fullfile(pData.folder{pp},pData.filename{pp});
- [im] = SSF.loadImage(strDB,1);
- imf{pp} = RSF.imPreProcess(im,radii);
- p1 = imf{pp}(:);
- p1 = p1(find(p1));
- px{pp} = p1;
-end
-
-px = px{:};
-px_pos = px(px>0);
-px_neg = abs(px(px<0));
-clip_pos = [mean(px_pos)-std(px_pos), mean(px_pos)+std(px_pos)];
-clip_neg = [mean(px_neg)-std(px_neg), mean(px_neg)+std(px_neg)];
-
-for i = 1:length(imf)
- imp = imf{i};
- imp(imp<0)=0;
- imn = imf{i};
- imn(imn>0)=0;
- imn=abs(imn);
- imq_pos = SSF.quantize8(imp,clip_pos);
- imq_neg = SSF.quantize8(imn,clip_neg);
- if all(0==imq_pos)
- imf{i} = imq_neg;
- else
- % [0,127] for imq_neg, [128,255] imq_pos
- imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
- imq_neg = 0.5 .* imq_neg ;
- imf{i} = imq_neg + imq_pos;
- end
-end
-4;
diff --git a/src/.DS_Store b/src/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..e7b88337ea51a860d173a7fce6d85e277ebebcf7
Binary files /dev/null and b/src/.DS_Store differ
diff --git a/src/+RSF/checkMetaTiming.m b/src/MATLAB/+RSF/checkMetaTiming.m
similarity index 100%
rename from src/+RSF/checkMetaTiming.m
rename to src/MATLAB/+RSF/checkMetaTiming.m
diff --git a/src/+RSF/chooseImageSamples.m b/src/MATLAB/+RSF/chooseImageSamples.m
similarity index 100%
rename from src/+RSF/chooseImageSamples.m
rename to src/MATLAB/+RSF/chooseImageSamples.m
diff --git a/src/+RSF/clipRSF.m b/src/MATLAB/+RSF/clipRSF.m
similarity index 100%
rename from src/+RSF/clipRSF.m
rename to src/MATLAB/+RSF/clipRSF.m
diff --git a/src/+RSF/cohortOCT.m b/src/MATLAB/+RSF/cohortOCT.m
similarity index 100%
rename from src/+RSF/cohortOCT.m
rename to src/MATLAB/+RSF/cohortOCT.m
diff --git a/src/+RSF/corrOCT.m b/src/MATLAB/+RSF/corrOCT.m
similarity index 100%
rename from src/+RSF/corrOCT.m
rename to src/MATLAB/+RSF/corrOCT.m
diff --git a/src/+RSF/draw1RSF.m b/src/MATLAB/+RSF/draw1RSF.m
similarity index 100%
rename from src/+RSF/draw1RSF.m
rename to src/MATLAB/+RSF/draw1RSF.m
diff --git a/src/+RSF/drawOCT.m b/src/MATLAB/+RSF/drawOCT.m
similarity index 100%
rename from src/+RSF/drawOCT.m
rename to src/MATLAB/+RSF/drawOCT.m
diff --git a/src/+RSF/drawPatients.m b/src/MATLAB/+RSF/drawPatients.m
similarity index 100%
rename from src/+RSF/drawPatients.m
rename to src/MATLAB/+RSF/drawPatients.m
diff --git a/src/+RSF/drawProgressionStats.m b/src/MATLAB/+RSF/drawProgressionStats.m
similarity index 100%
rename from src/+RSF/drawProgressionStats.m
rename to src/MATLAB/+RSF/drawProgressionStats.m
diff --git a/src/+RSF/drawRSF.m b/src/MATLAB/+RSF/drawRSF.m
similarity index 100%
rename from src/+RSF/drawRSF.m
rename to src/MATLAB/+RSF/drawRSF.m
diff --git a/src/+RSF/drawRSF_lever.m b/src/MATLAB/+RSF/drawRSF_lever.m
similarity index 100%
rename from src/+RSF/drawRSF_lever.m
rename to src/MATLAB/+RSF/drawRSF_lever.m
diff --git a/src/+RSF/getDeltaTstats.m b/src/MATLAB/+RSF/getDeltaTstats.m
similarity index 100%
rename from src/+RSF/getDeltaTstats.m
rename to src/MATLAB/+RSF/getDeltaTstats.m
diff --git a/src/+RSF/getImRSF.m b/src/MATLAB/+RSF/getImRSF.m
similarity index 100%
rename from src/+RSF/getImRSF.m
rename to src/MATLAB/+RSF/getImRSF.m
diff --git a/src/MATLAB/+RSF/getImf.m b/src/MATLAB/+RSF/getImf.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d0465867b0b0b807dcb3d3675a59ea2a5a3ed85
--- /dev/null
+++ b/src/MATLAB/+RSF/getImf.m
@@ -0,0 +1,88 @@
+function imf = getImf(pData,radii,scanType)
+
+% radii is 1-element or 3-element cell array. [0,0,0] for unused channel
+% FLIF takes 1 or 3 UBYTE
+radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
+% radii = {[0.5,8,8]}; % plate
+
+% NGPU = HIP.Cuda.DeviceCount;
+% p = ljsStartParallel(8*NGPU);
+
+N_CUDA_PROCS = 48;
+startTarget = 1;
+imf = {};
+im1 = Composite();
+pp = Composite();
+while startTarget <= size(pData,1)
+ spmd
+ pp = spmdIndex - 1 + startTarget;
+ if spmdIndex <= N_CUDA_PROCS && pp <= size(pData,1)
+ strDB = fullfile(pData.folder{pp},pData.filename{pp});
+ [im] = SSF.loadImage(strDB,1);
+ % 3-channel plate,blob,0
+ for c = 1:length(radii)
+ if all(0==radii{c})
+ im1(:,:,:,c) = zeros(size(im));
+ else
+ im1(:,:,:,c) = RSF.imPreProcess(im,radii{c});
+ end
+ end
+ end
+ end
+ startTarget = startTarget + N_CUDA_PROCS;
+ for i = 1:N_CUDA_PROCS
+ if pp{i} <= size(pData,1)
+ imf{pp{i}} = im1{i};
+ end
+ end
+end
+
+% px = {}; % critical voxels
+% mx = []; % mean by stack
+% sx = []; % sigma by stack
+% for i = 1 : length(imf)
+% for c = 1:size(imf{i},4)
+% i1 = imf{i}(:,:,:,c);
+% px{i,c} = i1(find(i1));
+% mx(i,c) = mean(px{i,c});
+% sx(i,c) = std(px{i,c});
+% end
+% end
+
+% clipRange = [2.5,97.5];
+% % clipRange = [0.5,99.5];
+% for c = 1:size(imf{1},4)
+% p1 = vertcat(px{:,c});
+% px_pos = p1(p1>0);
+% px_neg = abs(p1(p1<0));
+%
+% clip_pos(c,:) = prctile(px_pos,clipRange);
+% clip_neg(c,:) = prctile(px_neg,clipRange);
+% end
+% 4;
+% for i = 1:length(imf)
+% for c = 1 : size(imf{i},4)
+% if all(0==imf{i}(:,:,c))
+% continue
+% end
+% imp = imf{i}(:,:,:,c);
+% imp(imp<0)=0;
+% imn = imf{i}(:,:,:,c);
+% imn(imn>0)=0;
+% imn=abs(imn);
+% imq_pos = SSF.quantize8(imp,clip_pos(c,:));
+% imq_neg = SSF.quantize8(imn,clip_neg(c,:));
+% if all(0==imq_pos)
+% imf{i}(:,:,:,c) = imq_neg;
+% else
+% % [0,127] for imq_neg, [128,255] imq_pos
+% imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
+% imq_neg = 0.5 .* imq_neg ;
+% imf{i}(:,:,:,c) = (imq_neg + imq_pos);
+% end
+%
+% end
+%
+% imf{i} = uint8(imf{i});
+% end
+% 4;
diff --git a/src/+RSF/getMetadata.m b/src/MATLAB/+RSF/getMetadata.m
similarity index 100%
rename from src/+RSF/getMetadata.m
rename to src/MATLAB/+RSF/getMetadata.m
diff --git a/src/+RSF/getPatientMetadata.m b/src/MATLAB/+RSF/getPatientMetadata.m
similarity index 100%
rename from src/+RSF/getPatientMetadata.m
rename to src/MATLAB/+RSF/getPatientMetadata.m
diff --git a/src/+RSF/imPreProcess.m b/src/MATLAB/+RSF/imPreProcess.m
similarity index 51%
rename from src/+RSF/imPreProcess.m
rename to src/MATLAB/+RSF/imPreProcess.m
index 0fbc9cdf46ea29f8164f7f8d1bba92f7d25a1b0d..62b3528f80a74edfb58521adead772e09fd4fcd5 100644
--- a/src/+RSF/imPreProcess.m
+++ b/src/MATLAB/+RSF/imPreProcess.m
@@ -1,28 +1,11 @@
function [imf,bwMask] = imPreProcess(im,radii)
-% AX = 1/sqrt(3) .* [15.38,1,1] ;
-% AX = 1/sqrt(3) .* [5,1,1] ;
-AX = 1/sqrt(3);
-% dx = 15.38;
-dx = 20;
-AX = 1/sqrt(3);
-% fun fact :: mat2gray here causes r2 decline
im = uint8(im); % values are already uint8, just wrapped in a double
-
-imp = double(0*im);
-imn = imp;
-for rr = 1:length(radii)
- sx = radii(rr);
- if 1 == length(sx)
- sx = [0.45,sx,sx];
- end
- sx = AX .* sx;
- cudaTarget = getCudaTarget;
- imf = HIP.LoG(im,sx,cudaTarget);
- [~,imp1,imn1] = LoG.separateLoG(imf);
- imp = max(imp,imp1);
- imn = max(imn,imn1);
-end
+AX = 1/sqrt(3);
+sx = AX .* radii;
+cudaTarget = getCudaTarget;
+imf = HIP.LoG(im,sx,cudaTarget);
+[~,imp,imn] = LoG.separateLoG(imf);
% threshold mask
bw = imbinarize(im);
diff --git a/src/+RSF/optimizeDeltaTstats.m b/src/MATLAB/+RSF/optimizeDeltaTstats.m
similarity index 100%
rename from src/+RSF/optimizeDeltaTstats.m
rename to src/MATLAB/+RSF/optimizeDeltaTstats.m
diff --git a/src/+RSF/pickMinVelocityStack.m b/src/MATLAB/+RSF/pickMinVelocityStack.m
similarity index 100%
rename from src/+RSF/pickMinVelocityStack.m
rename to src/MATLAB/+RSF/pickMinVelocityStack.m
diff --git a/src/+RSF/pickRSFradius.m b/src/MATLAB/+RSF/pickRSFradius.m
similarity index 100%
rename from src/+RSF/pickRSFradius.m
rename to src/MATLAB/+RSF/pickRSFradius.m
diff --git a/src/+RSF/platemetric.m b/src/MATLAB/+RSF/platemetric.m
similarity index 100%
rename from src/+RSF/platemetric.m
rename to src/MATLAB/+RSF/platemetric.m
diff --git a/src/+RSF/progressionStats.m b/src/MATLAB/+RSF/progressionStats.m
similarity index 100%
rename from src/+RSF/progressionStats.m
rename to src/MATLAB/+RSF/progressionStats.m
diff --git a/src/MATLAB/+RSF/quantizePair.m b/src/MATLAB/+RSF/quantizePair.m
new file mode 100644
index 0000000000000000000000000000000000000000..7bda6410bb87bbc71dd82d2e3a096cb961337db6
--- /dev/null
+++ b/src/MATLAB/+RSF/quantizePair.m
@@ -0,0 +1,26 @@
+% clipLimits = [minVal,maxVal];
+function [i1q,i2q] = quantizePair(i1,i2)
+
+if isa(i1,'uint8') && isa(i2,'uint8')
+ return
+end
+
+pix = [i1(find(i1));i2(find(i2))];
+
+% clip at 99.9%
+clipLimits = prctile(pix,[0.05,99.95]);
+
+i1q = applyClip(i1,clipLimits);
+i2q = applyClip(i2,clipLimits);
+
+4;
+
+function im = applyClip(im,clipLimits)
+im_in = im; % remember the zeros!
+im = max(im,clipLimits(1));
+im = min(im,clipLimits(end));
+im = (im - clipLimits(1)) ./ (clipLimits(end) - clipLimits(1));
+im = uint8(254 * im + 1); % 8 bit quantization [0,1]->[1,255]
+im(0 == im_in) = 0;
+4;
+
diff --git a/src/+RSF/thresholdRetinalOCT.m b/src/MATLAB/+RSF/thresholdRetinalOCT.m
similarity index 100%
rename from src/+RSF/thresholdRetinalOCT.m
rename to src/MATLAB/+RSF/thresholdRetinalOCT.m
diff --git a/src/MATLAB/AddSQLiteToPath.m b/src/MATLAB/AddSQLiteToPath.m
new file mode 100644
index 0000000000000000000000000000000000000000..9f696fff585765875fd59ff97153da4260bb9764
--- /dev/null
+++ b/src/MATLAB/AddSQLiteToPath.m
@@ -0,0 +1,26 @@
+function AddSQLiteToPath()
+ %CHECKJARPATH Summary of this function goes here
+ % Detailed explanation goes here
+
+ if isdeployed()
+ return
+ end
+ %% ensure that the sqlite jar file is on the path
+ dynamicPaths = javaclasspath('-dynamic');
+ bfIsLoaded = false;
+ if (~isempty(dynamicPaths))
+ for i=1:length(dynamicPaths)
+ [~,name,~] = fileparts(dynamicPaths{i});
+ if (strcmpi('sqlite-jdbc-3.21.0',name))
+ bfIsLoaded = true;
+ break
+ end
+ end
+ end
+
+ if (~bfIsLoaded)
+ curPath = mfilename('fullpath');
+ [pathstr,~,~] = fileparts(curPath);
+ javaaddpath(fullfile(pathstr,'sqlite-jdbc-3.21.0.jar'),'-end');
+ end
+end
\ No newline at end of file
diff --git a/src/MATLAB/metadata/getMetaOCT.m b/src/MATLAB/metadata/getMetaOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..acf9fc5352cd522b236a7feea224c9c1a2d8d672
--- /dev/null
+++ b/src/MATLAB/metadata/getMetaOCT.m
@@ -0,0 +1,31 @@
+function tblMeta = getMetaOCT(ROOT,eye,scanType)
+
+tblMeta = table();
+flist = dir(fullfile(ROOT,'**/*.LEVER'));
+mx = regexp({flist.name},'P(\d+)_(.+?)_(\d+-\d+-\d+)_(\d+-\d+-\d+)_(\w\w)_(.+?)_.*.LEVER','tokens');
+tblSkip = table();
+for ff = 1:length(flist)
+ nt = table();
+
+ nt.filename = {flist(ff).name};
+ nt.folder = {flist(ff).folder};
+ if isempty(mx{ff})
+ tblSkip = [tblSkip;nt];
+ continue;
+ end
+ nt.subjectID = str2double(mx{ff}{1}{1});
+ nt.scanType = mx{ff}{1}{2};
+ nt.dx = datetime(mx{ff}{1}{3},'InputFormat','MM-dd-yyyy');
+ nt.time = {mx{ff}{1}{4}};
+ nt.eye = mx{ff}{1}{5};
+ nt.scanID = {mx{ff}{1}{6}};
+
+
+ if ~strcmp(nt.eye,eye) || ~contains(lower(scanType),lower(nt.scanType))
+ continue;
+ end
+
+ tblMeta = [tblMeta;nt];
+end
+
+tblMeta = sortrows(tblMeta,{'subjectID','filename','dx'});
\ No newline at end of file
diff --git a/src/MATLAB/metadata/goLR.m b/src/MATLAB/metadata/goLR.m
new file mode 100644
index 0000000000000000000000000000000000000000..317a5b01b2b8a73e8d633db79d49573f886b10c5
--- /dev/null
+++ b/src/MATLAB/metadata/goLR.m
@@ -0,0 +1,53 @@
+tic
+ROOT = '/g/leverjs/Schuman_OCT/OCT/03-22-2024';
+
+PID = 'P10010';
+flist = dir(fullfile(ROOT,['**/' PID '*.LEVER']));
+mx = regexp({flist.name},'P(?<PID>\d+)_(?<scanType>.+)_(?<date>\d+-\d+-\d+)_(?<time>\d+-\d+-\d+)_(?<eye>\w\w)_(?<scanID>.+?)_.*.LEVER','names');
+tblMeta = struct2table(vertcat(mx{:}));
+tblMeta.PID = str2double(tblMeta.PID);
+tblMeta.folder = {flist.folder}';
+tblMeta.filename = {flist.name}';
+
+% tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+tblMeta = tblMeta(contains(tblMeta.scanType,'Macula'),:);
+% tblMeta = tblMeta(contains(tblMeta.eye,'OD'),:);
+p = ljsStartParallel(96);
+
+imf = RSF.getImf(tblMeta,[[0.5,8,8] ; [5,5,5]],'');
+d = [];
+parfor i = 1:length(imf)
+ dj = [];
+ for j = 1:length(imf)
+ dj(j) = SSF.ncd_ssf_volume(imf{i},imf{j});
+ end
+ d(i,:) = dj;
+end
+
+A = Cluster.Regularize(d);
+[~,Y] = Cluster.SpectralCluster(A,3);
+
+tblMeta.Y = Y;
+
+tS = tblMeta(strcmp(tblMeta.eye,'OD'),:);
+tD = tblMeta(strcmp(tblMeta.eye,'OS'),:);
+t1 = tS;
+t2 = tD;
+
+tO = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+tM = tblMeta(contains(tblMeta.scanType,'Macular'),:);
+t3 = tO;
+t4 = tM;
+
+clf;hold on
+plot3(t1.Y(:,1),t1.Y(:,2),t1.Y(:,3),'r*')
+plot3(t2.Y(:,1),t2.Y(:,2),t2.Y(:,3),'og')
+legend({'OS','OD'})
+xlabel('NCD1')
+ylabel('NCD2')
+zlabel('NCD3')
+
+plot3(t3.Y(:,1),t3.Y(:,2),t3.Y(:,3),'mx')
+plot3(t4.Y(:,1),t4.Y(:,2),t4.Y(:,3),'cs')
+legend({'ONH','macula'})
+toc
\ No newline at end of file
diff --git a/src/MATLAB/metadata/leftVright_3_3_3_d.mat b/src/MATLAB/metadata/leftVright_3_3_3_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e8c2b5dcc20f3bd5c2bbad88466a6fff5b2f11f0
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_3_3_3_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_5_5_5_d.mat b/src/MATLAB/metadata/leftVright_5_5_5_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..015fa8b8a1368c302731d40fb34fd36be5149c36
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_5_5_5_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_8_8_8_d.mat b/src/MATLAB/metadata/leftVright_8_8_8_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..57742ab1a6ee57568aa8f92a4a2c6b56c97dfde2
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_8_8_8_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat b/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat
new file mode 100644
index 0000000000000000000000000000000000000000..244895e9cd0412e37b85d92b4cd48b4334341e6f
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat differ
diff --git a/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat
new file mode 100644
index 0000000000000000000000000000000000000000..945f0f87208afcfcf7c35011e92f8721df52035e
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat differ
diff --git a/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat
new file mode 100644
index 0000000000000000000000000000000000000000..11aef006c4b53ee4302b9c5506b8795f61c7dec7
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat differ
diff --git a/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat b/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4b2f570003e28261a4524e026bcbf788a3b84e33
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat differ
diff --git a/src/MATLAB/ncdVSvfmd/getMetaOCT.m b/src/MATLAB/ncdVSvfmd/getMetaOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..acf9fc5352cd522b236a7feea224c9c1a2d8d672
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getMetaOCT.m
@@ -0,0 +1,31 @@
+function tblMeta = getMetaOCT(ROOT,eye,scanType)
+
+tblMeta = table();
+flist = dir(fullfile(ROOT,'**/*.LEVER'));
+mx = regexp({flist.name},'P(\d+)_(.+?)_(\d+-\d+-\d+)_(\d+-\d+-\d+)_(\w\w)_(.+?)_.*.LEVER','tokens');
+tblSkip = table();
+for ff = 1:length(flist)
+ nt = table();
+
+ nt.filename = {flist(ff).name};
+ nt.folder = {flist(ff).folder};
+ if isempty(mx{ff})
+ tblSkip = [tblSkip;nt];
+ continue;
+ end
+ nt.subjectID = str2double(mx{ff}{1}{1});
+ nt.scanType = mx{ff}{1}{2};
+ nt.dx = datetime(mx{ff}{1}{3},'InputFormat','MM-dd-yyyy');
+ nt.time = {mx{ff}{1}{4}};
+ nt.eye = mx{ff}{1}{5};
+ nt.scanID = {mx{ff}{1}{6}};
+
+
+ if ~strcmp(nt.eye,eye) || ~contains(lower(scanType),lower(nt.scanType))
+ continue;
+ end
+
+ tblMeta = [tblMeta;nt];
+end
+
+tblMeta = sortrows(tblMeta,{'subjectID','filename','dx'});
\ No newline at end of file
diff --git a/src/MATLAB/ncdVSvfmd/getStats.m b/src/MATLAB/ncdVSvfmd/getStats.m
new file mode 100644
index 0000000000000000000000000000000000000000..56b12eb542b7ecab65d31adc089757bd54c9397e
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getStats.m
@@ -0,0 +1,54 @@
+
+function tblVelo = getStats(dxx, pData, idxDel, tblMD)
+if ~exist('idxDel','var')
+ idxDel = [];
+end
+tblVelo = table();
+idxUnused = 1:height(pData);
+for i = 1:length(dxx)
+ for j = 1:size(dxx{i},1)
+ m = dxx{i}(j,2);
+ n = dxx{i}(j,3);
+ if ~isempty(intersect(m,idxDel)) || ~isempty(intersect(n,idxDel))
+ continue
+ end
+
+
+ if m~=n && pData.dx(m) == pData.dx(n) % same date pairs
+ continue
+ end
+%
+
+ dmd = pData.md(n) - pData.md(m);
+ nt = table();
+ nt.subjectID = pData.subjectID(m);
+ idxm = find(dxx{i}(:,2)==m & dxx{i}(:,3) == m);
+ idxn = find(dxx{i}(:,2)==n & dxx{i}(:,3) == n);
+ nt.X = [dxx{i}(idxm,1),dxx{i}(idxn,1),dxx{i}(j,1)];
+ nt.velocity = (dxx{i}(j,1));
+ nt.md = min(pData.md(n),pData.md(m));
+ nt.dmd = abs(dmd);
+ nt.mn = [m,n];
+ nt.ddx = pData.dx(n) - pData.dx(m);
+ % find sex
+ idx = find(tblMD.subjectID == pData.subjectID(m),1,'first');
+ nt.sex = tblMD.gender_full(idx);
+ if nt.dmd > 4
+ continue
+ end
+ if years(nt.ddx)>1
+ continue
+ end
+% if nt.velocity > 0.9 && nt.velocity < 0.975
+% continue
+% end
+ if m ~= n
+ idxUnused = setdiff(idxUnused,m);
+ idxUnused = setdiff(idxUnused,n);
+ end
+ tblVelo = [tblVelo;nt];
+ end
+end
+
+% tblVelo = trimSelfSimilar(tblVelo);
+4;
diff --git a/src/MATLAB/ncdVSvfmd/getVFMD.m b/src/MATLAB/ncdVSvfmd/getVFMD.m
new file mode 100644
index 0000000000000000000000000000000000000000..92e08e5386fdad126d4713b353c322c7cf318dec
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getVFMD.m
@@ -0,0 +1,67 @@
+function tblVFMD = getVFMD(eye)
+
+if ~exist('eye','var')
+ eye = 'OD';
+end
+fx = {'nyu VF qualified data great than -6 data export.xlsx',...
+ 'nyu VF qualified data from -12 to -6 data export.xlsx',...
+ 'nyu VF qualified data less than -12 data export.xlsx'};
+% read from 3 files
+tblVFMD = table();
+for ff = 1 : length(fx)
+
+ tx = readtable(fullfile("importMetadata/dataExcel/",fx{ff}));
+ % filter 1
+ idxVFMD = contains(lower(tx.qualification_status),'yes');
+ tx = tx(idxVFMD,:);
+ % filter 2
+ idxVFMD = contains(lower(tx.threshold),lower('Visual Field 24-2 Test Pattern'));
+ tx = tx(idxVFMD,:);
+ % filter 3
+ for i = 1:height(tx)
+ tx.operatorError(i) = sum([tx.false_positive_p(i), tx.false_negative_p(i), tx.fixation_loss_c_p(i)]);
+ end
+ idxVFMD = tx.operatorError < 33;
+ tx = tx(idxVFMD,:);
+
+ % filter by eye
+ idxVFMD = contains(lower(tx.Eye),lower(eye));
+ tx = tx(idxVFMD,:);
+
+ tblVFMD = [tblVFMD;tx];
+end
+
+tblVFMD.subjectID = cellfun(@(x) str2double(x),tblVFMD.project_id);
+tblVFMD.dx = tblVFMD.visit_date_date;
+tblVFMD = sortrows(tblVFMD,{'subjectID','dx'});
+tblVFMD = unique(tblVFMD,'rows');
+
+% filter by date ,patient
+[pC, p_ia, p_ic] = unique(tblVFMD.subjectID);
+for pp = 1:length(pC)
+ idxVFMD = tblVFMD.subjectID == pC(pp);
+ tbl1 = tblVFMD(idxVFMD,:);
+ dx = unique(tbl1.visit_date_date);
+ idxDelAll = [];
+ for dd = 1:length(dx)
+ d1 = tbl1(tbl1.visit_date_date == dx(dd),:);
+ if 1 == height(d1)
+ continue
+ end
+ % need to choose best visit
+ px = [];
+ for vv = 1:height(d1)
+ px(vv,:) = [d1.false_positive_p(vv),d1.false_negative_p(vv),d1.fixation_loss_c_p(vv)];
+ end
+ px = sum(abs(px),2);
+ [~,idxMin] = min(px);
+ idxDel = find(idxVFMD & tblVFMD.dx == dx(dd));
+ idxDel(idxMin) = []; % preserve the most likely
+ idxDelAll = [idxDelAll,idxDel];
+
+ 4;
+ end
+ tblVFMD(idxDelAll,:) = []; % remove the least likely
+end
+tblVFMD.mdString = tblVFMD.md;
+tblVFMD.md = cellfun(@str2double,tblVFMD.md);
diff --git a/src/MATLAB/ncdVSvfmd/goKernelCorr.m b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..992b050f2472db01dd10eaef12a4cf0dc0e3c587
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
@@ -0,0 +1,42 @@
+datetime
+kernelCorrStart = tic();
+eye = 'OD';
+ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+% scanType = 'Macular Cube 200x200';
+scanType = 'Optic Disc Cube 200x200';
+
+load('vfmd_meta.mat');
+[rho2,tblVelo, pData, dxx] = kernelCorr([0.5,8,8],tblVFMD,tblMeta,scanType,eye);
+[rmse1,rmse2] = regnetCorr(tblVelo)
+%
+% tblResults = table();
+% for rx = 1:length(radii)
+% rr = radii{rx};
+% [rho2,tblVelo, pData, dxx] = kernelCorr(rr,tblVFMD,tblMeta,scanType,eye);
+% nt = table();
+% nt.rr = {rr};
+% nt.tblVelo = {tblVelo};
+% nt.pData = {pData};
+% nt.dxx = {dxx};
+% nt.r2 = rho2;
+% if iscell(rr)
+% rstr = cellfun(@mat2str,rr,'UniformOutput',false);
+% rstr = horzcat(rstr{:});
+% else
+% rstr = mat2str(rr);
+% end
+% nt.rstr = {rstr};
+% tblResults = [tblResults;nt]
+% end
+%
+% for i = 1:height(tblResults)
+% [rmse1,rmse2] = regnetCorr(tblResults.tblVelo{i});
+% tblResults.rmse(i,:) = [rmse1,rmse2];
+% tblVelo = tblResults.tblVelo{i};
+% x = tblVelo.velocity;
+% y = tblVelo.dmd;
+% yt = log10(1+y);
+% [f,gofX]=fit(x,yt,'poly1');
+% [f,gof]=fit(x,y,'poly1');
+% tblResults.r2a(i,1:2) = [gof.adjrsquare,gofX.adjrsquare]
+% end
diff --git a/src/MATLAB/ncdVSvfmd/kernelCorr.m b/src/MATLAB/ncdVSvfmd/kernelCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..8320c86c516c154329814d3781f560a91e0f8978
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/kernelCorr.m
@@ -0,0 +1,93 @@
+function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,scanType,eye)
+
+kernelCorrStart = tic();
+
+mdInfo = [];
+pid = unique(tblVFMD.subjectID);
+for i = 1:length(pid)
+ idx = find(tblVFMD.subjectID == pid(i));
+ visits = unique(tblVFMD.dx(idx));
+ if length(visits)<2,continue,end
+ md1 = [];
+ md1.id = pid(i);
+ md1.dx = tblVFMD.dx(idx);
+ md1.md = tblVFMD.md(idx);
+ mdInfo = [mdInfo,md1];
+end
+
+for i = 1:length(mdInfo)
+ idx = find(tblMeta.subjectID == mdInfo(i).id);
+ [C,ia,ib] = intersect(mdInfo(i).dx,tblMeta.dx(idx));
+ mdInfo(i).dxOCT = C;
+ mdInfo(i).dx = mdInfo(i).dx(ia);
+ mdInfo(i).md = mdInfo(i).md(ia);
+end
+
+octLen = arrayfun(@(x) length(x.dxOCT),mdInfo);
+mdInfo(octLen<2) = [];
+
+pData = table();
+for i = 1:length(mdInfo)
+ for pj = 1:length(mdInfo(i).dx)
+ idx = find(tblMeta.subjectID == mdInfo(i).id & tblMeta.dx == mdInfo(i).dx(pj));
+ p1 = tblMeta(idx,:);
+ p1.md = repmat(mdInfo(i).md(pj),height(p1),1);
+ pData = [pData;p1];
+ end
+end
+
+% pData = pData(1:20,:); % shorten for test
+
+p=ljsStartParallel(96);
+
+imf = RSF.getImf(pData,radii,scanType);
+
+[C,ia,ic] = unique(pData.subjectID);
+
+
+workList = [];
+for px = 1:length(C)
+ work1 = [];
+ idx = find(ic==px);
+
+ work1.imf = imf(idx);
+
+ wl = [idx,idx];
+ wl = [wl;nchoosek(idx,2)];
+ work1.wl = wl;
+
+ idx0 = [1:length(idx)]';
+ wl0 = [idx0,idx0];
+ wl0 = [wl0;nchoosek(idx0,2)];
+ work1.wl0 = wl0;
+
+ workList = [workList,work1];
+end
+
+
+dxx = {};
+parfor wx = 1:length(workList)
+ for ii = 1:size(workList(wx).wl0,1)
+ m = workList(wx).wl0(ii,1) ;
+ n = workList(wx).wl0(ii,2);
+ im_n = workList(wx).imf{n};
+ im_m = workList(wx).imf{m};
+ d1 = NCD.ncd_ssf_volume(im_m,im_n);
+ if m==n
+ d2 = d1;
+ else
+ d2 = NCD.ncd_ssf_volume(im_n,im_m);
+ end
+ mm = workList(wx).wl(ii,1) ;
+ nn = workList(wx).wl(ii,2) ;
+ dxx{wx}(ii,:) = [min(d1,d2),mm,nn];
+ end
+end
+
+idxDel = [];
+
+tblVelo = getStats(dxx,pData,idxDel,tblVFMD);
+[rho,pCorr] = corr((tblVelo.velocity),(tblVelo.dmd));
+rho2 = rho^2
+toc(kernelCorrStart)
+4;
diff --git a/src/MATLAB/ncdVSvfmd/regnetCorr.m b/src/MATLAB/ncdVSvfmd/regnetCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..cbcc894a9766bacd4c4c9eb7c89df7f9639736c9
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/regnetCorr.m
@@ -0,0 +1,16 @@
+function [rmse1,rmse2] = regnetCorr(tblVelo)
+% tblVelo = getStats(dxx,pData,[],tblVFMD);
+% X = [tblVelo.velocity,tblVelo.md];
+t2 = tblVelo(tblVelo.dmd~=0,:);
+
+X1 = [t2.velocity];
+X2 = [t2.velocity,t2.md];
+T = t2.dmd;
+
+mdl1 = fitrnet(X1,T,'Activations','sigmoid','LayerSizes',[1]);
+cv1 = mdl1.crossval; rmse1 = cv1.kfoldLoss^0.5;
+mdl2 = fitrnet(X2,T,'Activations','sigmoid','LayerSizes',[10]);
+cv2 = mdl2.crossval; rmse2 = cv2.kfoldLoss^0.5;
+4;
+% mdl = fitrnet(X,T,"OptimizeHyperparameters","auto", "HyperparameterOptimizationOptions",struct("AcquisitionFunctionName","expected-improvement-plus"))
+
diff --git a/src/MATLAB/ncdVSvfmd/vfmd_meta.mat b/src/MATLAB/ncdVSvfmd/vfmd_meta.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9cffe29aed7df086bee08ddf99448c753115f65d
Binary files /dev/null and b/src/MATLAB/ncdVSvfmd/vfmd_meta.mat differ