diff --git a/build.gradle.kts b/build.gradle.kts index de4e3eea5..c31deb364 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -22,6 +22,7 @@ java { } repositories { + mavenCentral() if(project.properties["useMavenLocal"] == "true") { logger.warn("Using local Maven repository as source") mavenLocal() @@ -44,7 +45,7 @@ dependencies { exclude("org.lwjgl") } - val sceneryVersion = "0.11.2" + val sceneryVersion = "0.12.0" api("graphics.scenery:scenery:$sceneryVersion") { version { strictly(sceneryVersion) } exclude("org.biojava.thirdparty", "forester") diff --git a/src/main/java/sc/iview/commands/file/OpenDirOfTIFF.java b/src/main/java/sc/iview/commands/file/OpenDirOfTIFF.java new file mode 100644 index 000000000..12ad0ff6a --- /dev/null +++ b/src/main/java/sc/iview/commands/file/OpenDirOfTIFF.java @@ -0,0 +1,85 @@ +/*- + * #%L + * Scenery-backed 3D visualization package for ImageJ. + * %% + * Copyright (C) 2016 - 2021 SciView developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +package sc.iview.commands.file; + +import org.scijava.command.Command; +import org.scijava.io.IOService; +import org.scijava.log.LogService; +import org.scijava.plugin.Menu; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +import sc.iview.SciView; + +import java.io.File; +import java.io.IOException; + +import static sc.iview.commands.MenuWeights.FILE; +import static sc.iview.commands.MenuWeights.FILE_OPEN; + +/** + * Command to open a file in SciView + * + * @author Kyle Harrington + * + */ +@Plugin(type = Command.class, menuRoot = "SciView", // + menu = { @Menu(label = "File", weight = FILE), // + @Menu(label = "Open Directory of tif...", weight = FILE_OPEN) }) +public class OpenDirOfTIFF implements Command { + + @Parameter + private IOService io; + + @Parameter + private LogService log; + + @Parameter + private SciView sciView; + + // TODO: Find a more extensible way than hard-coding the extensions. + @Parameter(style = "directory") + private File file; + + @Parameter + private int onlyFirst = 0; + + @Override + public void run() { + try { + if(onlyFirst > 0) { + sciView.openDirTiff(file.toPath(), onlyFirst); + } else { + sciView.openDirTiff(file.toPath(), null); + } + } + catch (final IOException | IllegalArgumentException exc) { + log.error( exc ); + } + } +} diff --git a/src/main/kotlin/sc/iview/SciView.kt b/src/main/kotlin/sc/iview/SciView.kt index 5e7d92578..2b44d86f8 100644 --- a/src/main/kotlin/sc/iview/SciView.kt +++ b/src/main/kotlin/sc/iview/SciView.kt @@ -55,9 +55,9 @@ import graphics.scenery.utils.ExtractsNatives.Companion.getPlatform import graphics.scenery.utils.LogbackUtils import graphics.scenery.utils.SceneryPanel import graphics.scenery.utils.Statistics -import graphics.scenery.utils.extensions.times import graphics.scenery.volumes.Colormap import graphics.scenery.volumes.RAIVolume +import graphics.scenery.volumes.TransferFunction import graphics.scenery.volumes.Volume import graphics.scenery.volumes.Volume.Companion.fromXML import graphics.scenery.volumes.Volume.Companion.setupId @@ -85,6 +85,7 @@ import net.imglib2.type.numeric.integer.UnsignedByteType import net.imglib2.view.Views import org.joml.Quaternionf import org.joml.Vector3f +import org.joml.Vector4f import org.scijava.Context import org.scijava.`object`.ObjectService import org.scijava.display.Display @@ -100,6 +101,7 @@ import org.scijava.thread.ThreadService import org.scijava.util.ColorRGB import org.scijava.util.Colors import org.scijava.util.VersionUtils +import sc.iview.commands.demo.animation.ParticleDemo import sc.iview.commands.edit.InspectorInteractiveCommand import sc.iview.event.NodeActivatedEvent import sc.iview.event.NodeAddedEvent @@ -110,7 +112,6 @@ import sc.iview.ui.CustomPropertyUI import sc.iview.ui.MainWindow import sc.iview.ui.SwingMainWindow import sc.iview.ui.TaskManager -import ucar.units.ConversionException import java.awt.event.WindowListener import java.io.File import java.io.IOException @@ -118,6 +119,7 @@ import java.net.JarURLConnection import java.net.URL import java.nio.ByteBuffer import java.nio.FloatBuffer +import java.nio.file.Path import java.time.LocalDate import java.time.format.DateTimeFormatter import java.util.* @@ -134,6 +136,7 @@ import kotlin.concurrent.thread import javax.swing.JOptionPane import kotlin.math.cos import kotlin.math.sin +import kotlin.system.measureTimeMillis /** * Main SciView class. @@ -784,6 +787,25 @@ class SciView : SceneryBase, CalibratedRealInterval { } } + @Throws(IOException::class) + fun openDirTiff(source: Path, onlyFirst: Int? = null) + { + val v = Volume.fromPath(source, hub, onlyFirst) + v.name = "volume" + v.spatial().position = Vector3f(-3.0f, 10.0f, 0.0f) + v.colormap = Colormap.get("jet") + v.spatial().scale = Vector3f(15.0f, 15.0f,45.0f) + v.transferFunction = TransferFunction.ramp(0.05f, 0.8f) + v.metadata["animating"] = true + v.converterSetups.firstOrNull()?.setDisplayRange(0.0, 1500.0) + v.visible = true + + v.spatial().wantsComposeModel = true + v.spatial().updateWorld(true) + addNode(v) + } + + /** * Open a file specified by the source path. The file can be anything that SciView knows about: mesh, volume, point cloud * @param source string of a data source @@ -1697,28 +1719,41 @@ class SciView : SceneryBase, CalibratedRealInterval { val cam = scene.activeObserver as? DetachedHeadCamera ?: return var ti: TrackerInput? = null var hmdAdded = false - if (!hub.has(SceneryElement.HMDInput)) { - try { - val hmd = OpenVRHMD(false, true) - if (hmd.initializedAndWorking()) { - hub.add(SceneryElement.HMDInput, hmd) - ti = hmd - } else { - logger.warn("Could not initialise VR headset, just activating stereo rendering.") + + if (vrActive) { + + // VR activation logic + if (!hub.has(SceneryElement.HMDInput)) { + try { + val hmd = OpenVRHMD(false, true) + if (hmd.initializedAndWorking()) { + hub.add(SceneryElement.HMDInput, hmd) + ti = hmd + } else { + logger.warn("Could not initialise VR headset, just activating stereo rendering.") + } + } catch (e: Exception) { + logger.error("Could not add OpenVRHMD: $e") } hmdAdded = true - } catch (e: Exception) { - logger.error("Could not add OpenVRHMD: $e") + } else { + ti = hub.getWorkingHMD() + } + + // Set tracker on the DetachedHeadCamera + if (ti != null) { + cam.tracker = ti + logger.info("tracker set") } } else { - ti = hub.getWorkingHMD() - } - if (vrActive && ti != null) { - cam.tracker = ti - } else { + // VR deactivation logic + // Convert back to normal Camera + logger.info("Shutting down VR") cam.tracker = null } - renderer.pushMode = false + + // Enable push mode if VR is inactive, and the other way round + renderer.pushMode = !vrActive // we need to force reloading the renderer as the HMD might require device or instance extensions if (renderer is VulkanRenderer && hmdAdded) { @@ -1732,8 +1767,17 @@ class SciView : SceneryBase, CalibratedRealInterval { e.printStackTrace() } } - - renderer.toggleVR() + } + logger.debug("Replaced renderer.") + renderer.toggleVR() + // Cleanup HMD after VR has been toggled off + if (!vrActive) { + if (hub.has(SceneryElement.HMDInput)) { + val hmd = hub.get(SceneryElement.HMDInput) as? OpenVRHMD + hmd?.close() + // TODO hub.remove(hmd) + logger.debug("Closed HMD.") + } } } diff --git a/src/main/kotlin/sc/iview/commands/analysis/CellTrackingBase.kt b/src/main/kotlin/sc/iview/commands/analysis/CellTrackingBase.kt new file mode 100644 index 000000000..7bb7ad359 --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/CellTrackingBase.kt @@ -0,0 +1,471 @@ +package sc.iview.commands.analysis + +import graphics.scenery.DetachedHeadCamera +import graphics.scenery.Icosphere +import graphics.scenery.InstancedNode +import graphics.scenery.Mesh +import graphics.scenery.ShaderMaterial +import graphics.scenery.controls.OpenVRHMD +import graphics.scenery.controls.TrackerRole +import graphics.scenery.controls.behaviours.ControllerDrag +import graphics.scenery.primitives.Cylinder +import graphics.scenery.utils.MaybeIntersects +import graphics.scenery.utils.SystemHelpers +import graphics.scenery.utils.extensions.minus +import graphics.scenery.utils.lazyLogger +import graphics.scenery.volumes.RAIVolume +import graphics.scenery.volumes.Volume +import org.joml.Math +import org.joml.Matrix4f +import org.joml.Vector3f +import org.joml.Vector4f +import org.scijava.ui.behaviour.ClickBehaviour +import sc.iview.SciView +import sc.iview.commands.analysis.ConfirmableClickBehaviour +import sc.iview.commands.analysis.SpineMetadata +import sc.iview.commands.analysis.TimepointObserver +import java.io.BufferedWriter +import java.io.FileWriter +import java.nio.file.Path +import java.util.concurrent.atomic.AtomicInteger +import kotlin.concurrent.thread + +/** + * Base class for different VR cell tracking purposes. It includes functionality to add spines and edgehogs, + * as used by [EyeTracking], and registers controller bindings via [inputSetup]. It is possible to register observers + * that listen to timepoint changes with [registerObserver]. + * @param [sciview] The [sc.iview.SciView] instance to use + */ +open class CellTrackingBase( + open var sciview: SciView +) { + val logger by lazyLogger() + + lateinit var sessionId: String + lateinit var sessionDirectory: Path + + lateinit var hmd: OpenVRHMD + + val hedgehogs = Mesh() + val hedgehogIds = AtomicInteger(0) + lateinit var volume: Volume + + val referenceTarget = Icosphere(0.004f, 2) + + @Volatile var tracking = false + var playing = true + var direction = PlaybackDirection.Backward + var volumesPerSecond = 1 + var skipToNext = false + var skipToPrevious = false + + var volumeScaleFactor = 1.0f + + // determines whether the volume and hedgehogs should keep listening for updates or not + var cellTrackingActive: Boolean = false + + open var linkCreationCallback: ((HedgehogAnalysis.SpineGraphVertex) -> Unit)? = null + open var finalTrackCallback: (() -> Unit)? = null + + /** How to render the currently active hedgehog: all spines, only spines from the current time point, or none. */ + enum class HedgehogVisibility { Hidden, PerTimePoint, Visible } + + var hedgehogVisibility = HedgehogVisibility.Hidden + + enum class PlaybackDirection { + Forward, + Backward + } + + private val observers = mutableListOf() + + /** Registers a new observer that will get updated whenever the VR user triggers a timepoint update. */ + fun registerObserver(observer: TimepointObserver) { + observers.add(observer) + } + + /** Unregisters the timepoint observer. */ + fun unregisterObserver(observer: TimepointObserver) { + observers.remove(observer) + } + + private fun notifyObservers(timepoint: Int) { + observers.forEach { it.onTimePointChanged(timepoint) } + } + + /** Initialize new hedgehog, used for collecting spines (gaze rays) during eye tracking. */ + fun addHedgehog() { + logger.info("added hedgehog") + val hedgehog = Cylinder(0.005f, 1.0f, 16) + hedgehog.visible = false + hedgehog.setMaterial(ShaderMaterial.Companion.fromFiles("DeferredInstancedColor.frag", "DeferredInstancedColor.vert")) + val hedgehogInstanced = InstancedNode(hedgehog) + hedgehogInstanced.instancedProperties["ModelMatrix"] = { hedgehog.spatial().world} + hedgehogInstanced.instancedProperties["Metadata"] = { Vector4f(0.0f, 0.0f, 0.0f, 0.0f) } + hedgehogs.addChild(hedgehogInstanced) + } + + open fun inputSetup() + { + val cam = sciview.camera ?: throw IllegalStateException("Could not find camera") + + sciview.sceneryInputHandler?.let { handler -> + hashMapOf( + "move_forward_fast" to (TrackerRole.LeftHand to OpenVRHMD.OpenVRButton.Up), + "move_back_fast" to (TrackerRole.LeftHand to OpenVRHMD.OpenVRButton.Down), + "move_left_fast" to (TrackerRole.LeftHand to OpenVRHMD.OpenVRButton.Left), + "move_right_fast" to (TrackerRole.LeftHand to OpenVRHMD.OpenVRButton.Right)).forEach { (name, key) -> + handler.getBehaviour(name)?.let { b -> + hmd.addBehaviour(name, b) + hmd.addKeyBinding(name, key.first, key.second) + + } + } + } + + val toggleHedgehog = ClickBehaviour { _, _ -> + val current = HedgehogVisibility.entries.indexOf(hedgehogVisibility) + hedgehogVisibility = HedgehogVisibility.entries.get((current + 1) % 3) + + when (hedgehogVisibility) { + HedgehogVisibility.Hidden -> { + hedgehogs.visible = false + hedgehogs.runRecursive { it.visible = false } + cam.showMessage("Hedgehogs hidden", distance = 2f, size = 0.2f, centered = true) + } + + HedgehogVisibility.PerTimePoint -> { + hedgehogs.visible = true + cam.showMessage("Hedgehogs shown per timepoint", distance = 2f, size = 0.2f, centered = true) + } + + HedgehogVisibility.Visible -> { + hedgehogs.visible = true + cam.showMessage("Hedgehogs visible", distance = 2f, size = 0.2f, centered = true) + } + } + } + + val nextTimepoint = ClickBehaviour { _, _ -> + skipToNext = true + } + + val prevTimepoint = ClickBehaviour { _, _ -> + skipToPrevious = true + } + + val fasterOrScale = ClickBehaviour { _, _ -> + if (playing) { + volumesPerSecond = maxOf(minOf(volumesPerSecond + 1, 20), 1) + cam.showMessage("Speed: $volumesPerSecond vol/s", distance = 1.2f, size = 0.2f, centered = true) + } else { + volumeScaleFactor = minOf(volumeScaleFactor * 1.2f, 3.0f) + volume.spatial().scale = Vector3f(1.0f).mul(volumeScaleFactor) + } + } + + val slowerOrScale = ClickBehaviour { _, _ -> + if (playing) { + volumesPerSecond = maxOf(minOf(volumesPerSecond - 1, 20), 1) + cam.showMessage("Speed: $volumesPerSecond vol/s", distance = 2f, size = 0.2f, centered = true) + } else { + volumeScaleFactor = maxOf(volumeScaleFactor / 1.2f, 0.1f) + volume.spatial().scale = Vector3f(1.0f).mul(volumeScaleFactor) + } + } + + val playPause = ClickBehaviour { _, _ -> + playing = !playing + if (playing) { + cam.showMessage("Playing", distance = 2f, size = 0.2f, centered = true) + } else { + cam.showMessage("Paused", distance = 2f, size = 0.2f, centered = true) + } + } + + val move = ControllerDrag(TrackerRole.LeftHand, hmd) { volume } + + val deleteLastHedgehog = ConfirmableClickBehaviour( + armedAction = { timeout -> + cam.showMessage( + "Deleting last track, press again to confirm.", distance = 2f, size = 0.2f, + messageColor = Vector4f(1.0f, 1.0f, 1.0f, 1.0f), + backgroundColor = Vector4f(1.0f, 0.2f, 0.2f, 1.0f), + duration = timeout.toInt(), + centered = true + ) + + }, + confirmAction = { + hedgehogs.children.removeLast() + volume.children.last { it.name.startsWith("Track-") }?.let { lastTrack -> + volume.removeChild(lastTrack) + } + val hedgehogId = hedgehogIds.get() + val hedgehogFile = + sessionDirectory.resolve("Hedgehog_${hedgehogId}_${SystemHelpers.Companion.formatDateTime()}.csv") + .toFile() + val hedgehogFileWriter = BufferedWriter(FileWriter(hedgehogFile, true)) + hedgehogFileWriter.newLine() + hedgehogFileWriter.newLine() + hedgehogFileWriter.write("# WARNING: TRACK $hedgehogId IS INVALID\n") + hedgehogFileWriter.close() + + cam.showMessage( + "Last track deleted.", distance = 2f, size = 0.2f, + messageColor = Vector4f(1.0f, 0.2f, 0.2f, 1.0f), + backgroundColor = Vector4f(1.0f, 1.0f, 1.0f, 1.0f), + duration = 1000, + centered = true + ) + }) + + hmd.addBehaviour("playback_direction", ClickBehaviour { _, _ -> + direction = if (direction == PlaybackDirection.Forward) { + PlaybackDirection.Backward + } else { + PlaybackDirection.Forward + } + cam.showMessage("Playing: ${direction}", distance = 2f, centered = true) + }) + + val cellDivision = ClickBehaviour { _, _ -> + cam.showMessage("Adding cell division", distance = 2f, duration = 1000) + dumpHedgehog() + addHedgehog() + } + + hmd.addBehaviour("skip_to_next", nextTimepoint) + hmd.addBehaviour("skip_to_prev", prevTimepoint) + hmd.addBehaviour("faster_or_scale", fasterOrScale) + hmd.addBehaviour("slower_or_scale", slowerOrScale) + hmd.addBehaviour("play_pause", playPause) + hmd.addBehaviour("toggle_hedgehog", toggleHedgehog) + hmd.addBehaviour("delete_hedgehog", deleteLastHedgehog) + hmd.addBehaviour("trigger_move", move) + hmd.addBehaviour("cell_division", cellDivision) + + hmd.addKeyBinding("skip_to_next", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Right) + hmd.addKeyBinding("skip_to_prev", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Left) + hmd.addKeyBinding("faster_or_scale", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Up) + hmd.addKeyBinding("slower_or_scale", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Down) + hmd.addKeyBinding("play_pause", TrackerRole.LeftHand, OpenVRHMD.OpenVRButton.Menu) + hmd.addKeyBinding("toggle_hedgehog", TrackerRole.LeftHand, OpenVRHMD.OpenVRButton.Side) + hmd.addKeyBinding("delete_hedgehog", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Side) + hmd.addKeyBinding("playback_direction", TrackerRole.RightHand, OpenVRHMD.OpenVRButton.Menu) + hmd.addKeyBinding("cell_division", TrackerRole.LeftHand, OpenVRHMD.OpenVRButton.Trigger) + + hmd.allowRepeats += OpenVRHMD.OpenVRButton.Trigger to TrackerRole.LeftHand + logger.info("Registered VR controller bindings.") + + } + + /** + * Launches a thread that updates the volume time points, the hedgehog visibility and reference target color. + */ + fun launchUpdaterThread() { + thread { + while (!sciview.isInitialized) { + Thread.sleep(200) + } + + while (sciview.running && cellTrackingActive) { + if (playing || skipToNext || skipToPrevious) { + val oldTimepoint = volume.viewerState.currentTimepoint + if (skipToNext || playing) { + skipToNext = false + if (direction == PlaybackDirection.Forward) { + notifyObservers(oldTimepoint + 1) + } else { + notifyObservers(oldTimepoint - 1) + } + } else { + skipToPrevious = false + if (direction == PlaybackDirection.Forward) { + notifyObservers(oldTimepoint - 1) + } else { + notifyObservers(oldTimepoint + 1) + } + } + val newTimepoint = volume.viewerState.currentTimepoint + + + if (hedgehogs.visible) { + if (hedgehogVisibility == HedgehogVisibility.PerTimePoint) { + hedgehogs.children.forEach { hh -> + val hedgehog = hh as InstancedNode + hedgehog.instances.forEach { + if (it.metadata.isNotEmpty()) { + it.visible = + (it.metadata["spine"] as SpineMetadata).timepoint == volume.viewerState.currentTimepoint + } + } + } + } else { + hedgehogs.children.forEach { hh -> + val hedgehog = hh as InstancedNode + hedgehog.instances.forEach { it.visible = true } + } + } + } + + if (tracking && oldTimepoint == (volume.timepointCount - 1) && newTimepoint == 0) { + tracking = false + + referenceTarget.ifMaterial { diffuse = Vector3f(0.5f, 0.5f, 0.5f) } + sciview.camera!!.showMessage("Tracking deactivated.", distance = 1.2f, size = 0.2f) + dumpHedgehog() + } + } + + Thread.sleep((1000.0f / volumesPerSecond).toLong()) + } + logger.info("CellTracking updater thread has stopped.") + } + } + + /** Adds a single spine to the currently active hedgehog. */ + open fun addSpine(center: Vector3f, direction: Vector3f, volume: Volume, confidence: Float, timepoint: Int) { + val cam = sciview.camera as? DetachedHeadCamera ?: return + val sphere = volume.boundingBox?.getBoundingSphere() ?: return + + val sphereDirection = sphere.origin.minus(center) + val sphereDist = Math.sqrt(sphereDirection.x * sphereDirection.x + sphereDirection.y * sphereDirection.y + sphereDirection.z * sphereDirection.z) - sphere.radius + + val p1 = center + val temp = direction.mul(sphereDist + 2.0f * sphere.radius) + val p2 = Vector3f(center).add(temp) + + val spine = (hedgehogs.children.last() as InstancedNode).addInstance() + spine.spatial().orientBetweenPoints(p1, p2, true, true) + spine.visible = false + + val intersection = volume.spatial().intersectAABB(p1, (p2 - p1).normalize(), true) + + if (volume.boundingBox?.isInside(cam.spatial().position)!!) { + logger.info("Can't track inside the volume! Please move out of the volume and try again") + return + } + if(intersection is MaybeIntersects.Intersection) { + // get local entry and exit coordinates, and convert to UV coords + val localEntry = (intersection.relativeEntry) + val localExit = (intersection.relativeExit) + // TODO We dont need the local direction for grid traversal, but its still in the spine metadata for now + val localDirection = Vector3f(0f) + val (samples, samplePos) = volume.sampleRayGridTraversal(localEntry, localExit) ?: (null to null) + val volumeScale = (volume as RAIVolume).getVoxelScale() + + if (samples != null && samplePos != null) { + val metadata = SpineMetadata( + timepoint, + center, + direction, + intersection.distance, + localEntry, + localExit, + localDirection, + cam.headPosition, + cam.headOrientation, + cam.spatial().position, + confidence, + samples.map { it ?: 0.0f }, + samplePos.map { it?.mul(volumeScale) ?: Vector3f(0f) } + ) + val count = samples.filterNotNull().count { it > 0.2f } + + spine.metadata["spine"] = metadata + spine.instancedProperties["ModelMatrix"] = { spine.spatial().world } + // TODO: Show confidence as color for the spine + spine.instancedProperties["Metadata"] = + { Vector4f(confidence, timepoint.toFloat() / volume.timepointCount, count.toFloat(), 0.0f) } + } + } + } + + /** + * Dumps a given hedgehog including created tracks to a file. + * If [hedgehog] is null, the last created hedgehog will be used, otherwise the given one. + * If [hedgehog] is not null, the cell track will not be added to the scene. + */ + fun dumpHedgehog(){ + logger.info("dumping hedgehog...") + val lastHedgehog = hedgehogs.children.last() as InstancedNode + val hedgehogId = hedgehogIds.incrementAndGet() + + val hedgehogFile = sessionDirectory.resolve("Hedgehog_${hedgehogId}_${SystemHelpers.Companion.formatDateTime()}.csv").toFile() + val hedgehogFileWriter = hedgehogFile.bufferedWriter() + hedgehogFileWriter.write("Timepoint,Origin,Direction,LocalEntry,LocalExit,LocalDirection,HeadPosition,HeadOrientation,Position,Confidence,Samples\n") + + val trackFile = sessionDirectory.resolve("Tracks.tsv").toFile() + val trackFileWriter = BufferedWriter(FileWriter(trackFile, true)) + if(!trackFile.exists()) { + trackFile.createNewFile() + trackFileWriter.write("# BionicTracking cell track listing for ${sessionDirectory.fileName}\n") + trackFileWriter.write("# TIME\tX\tYt\t\tZ\tTRACK_ID\tPARENT_TRACK_ID\tSPOT\tLABEL\n") + } + + + val spines = lastHedgehog.instances.mapNotNull { spine -> + spine.metadata["spine"] as? SpineMetadata + } + + spines.forEach { metadata -> + hedgehogFileWriter.write("${metadata.timepoint};${metadata.origin};${metadata.direction};${metadata.localEntry};${metadata.localExit};${metadata.localDirection};${metadata.headPosition};${metadata.headOrientation};${metadata.position};${metadata.confidence};${metadata.samples.joinToString(";")}\n") + } + hedgehogFileWriter.close() + + val existingAnalysis = lastHedgehog.metadata["HedgehogAnalysis"] as? HedgehogAnalysis.Track + val track = if(existingAnalysis is HedgehogAnalysis.Track) { + existingAnalysis + } else { + val h = HedgehogAnalysis(spines, Matrix4f(volume.spatial().world)) + h.run() + } + + if(track == null) { + logger.warn("No track returned") + sciview.camera?.showMessage("No track returned", distance = 1.2f, size = 0.2f,messageColor = Vector4f( + 1.0f, + 0.0f, + 0.0f, + 1.0f + ) + ) + return + } + + lastHedgehog.metadata["HedgehogAnalysis"] = track + lastHedgehog.metadata["Spines"] = spines + + val parentId = 0 + val volumeDimensions = volume.getDimensions() + + trackFileWriter.newLine() + trackFileWriter.newLine() + trackFileWriter.write("# START OF TRACK $hedgehogId, child of $parentId\n") + if (linkCreationCallback != null && finalTrackCallback != null) { + track.points.windowed(2, 1).forEach { pair -> + linkCreationCallback?.let { it(pair[0].second) } + val p = Vector3f(pair[0].first).mul(Vector3f(volumeDimensions)) // direct product + val tp = pair[0].second.timepoint + trackFileWriter.write("$tp\t${p.x()}\t${p.y()}\t${p.z()}\t${hedgehogId}\t$parentId\t0\t0\n") + } + finalTrackCallback?.invoke() + } else { + track.points.windowed(2, 1).forEach { pair -> + val p = Vector3f(pair[0].first).mul(Vector3f(volumeDimensions)) // direct product + val tp = pair[0].second.timepoint + trackFileWriter.write("$tp\t${p.x()}\t${p.y()}\t${p.z()}\t${hedgehogId}\t$parentId\t0\t0\n") + } + } + trackFileWriter.close() + } + + /** + * Stops the current tracking environment and restore the original state. + * This method needs to be overridden. + */ + open fun stop() { + } + +} diff --git a/src/main/kotlin/sc/iview/commands/analysis/ConfirmableClickBehaviour.kt b/src/main/kotlin/sc/iview/commands/analysis/ConfirmableClickBehaviour.kt new file mode 100644 index 000000000..df09ed4be --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/ConfirmableClickBehaviour.kt @@ -0,0 +1,34 @@ +package sc.iview.commands.analysis + +import org.scijava.ui.behaviour.ClickBehaviour +import kotlin.concurrent.thread + +/** + * [org.scijava.ui.behaviour.ClickBehaviour] that waits [timeout] for confirmation by re-executing the behaviour. + * Executes [armedAction] on first invocation, and [confirmAction] on second invocation, if + * it happens within [timeout]. + * + * @author Ulrik Guenther + */ +class ConfirmableClickBehaviour(val armedAction: (Long) -> Any, val confirmAction: (Long) -> Any, var timeout: Long = 3000): + ClickBehaviour { + /** Whether the action is armed at the moment. Action becomes disarmed after [timeout]. */ + private var armed: Boolean = false + + /** + * Action fired at position [x]/[y]. Parameters not used in VR actions. + */ + override fun click(x : Int, y : Int) { + if(!armed) { + armed = true + armedAction.invoke(timeout) + + thread { + Thread.sleep(timeout) + armed = false + } + } else { + confirmAction.invoke(timeout) + } + } +} \ No newline at end of file diff --git a/src/main/kotlin/sc/iview/commands/analysis/EyeTracking.kt b/src/main/kotlin/sc/iview/commands/analysis/EyeTracking.kt new file mode 100644 index 000000000..0ad8f3bdd --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/EyeTracking.kt @@ -0,0 +1,346 @@ +package sc.iview.commands.analysis + +import graphics.scenery.BoundingGrid +import graphics.scenery.Box +import graphics.scenery.BufferUtils +import graphics.scenery.DetachedHeadCamera +import graphics.scenery.Icosphere +import graphics.scenery.Light +import graphics.scenery.Mesh +import graphics.scenery.PointLight +import graphics.scenery.attribute.material.Material +import graphics.scenery.controls.OpenVRHMD +import graphics.scenery.controls.TrackedDeviceType +import graphics.scenery.controls.TrackerRole +import graphics.scenery.controls.eyetracking.PupilEyeTracker +import graphics.scenery.primitives.Cylinder +import graphics.scenery.primitives.TextBoard +import graphics.scenery.textures.Texture +import graphics.scenery.utils.SystemHelpers +import graphics.scenery.utils.extensions.minus +import graphics.scenery.utils.extensions.xyz +import graphics.scenery.utils.extensions.xyzw +import graphics.scenery.volumes.Volume +import net.imglib2.type.numeric.integer.UnsignedByteType +import org.joml.Matrix4f +import org.joml.Vector2f +import org.joml.Vector3f +import org.joml.Vector3i +import org.joml.Vector4f +import org.scijava.ui.behaviour.ClickBehaviour +import sc.iview.SciView +import java.awt.image.DataBufferByte +import java.io.ByteArrayInputStream +import java.nio.file.Files +import java.nio.file.Paths +import javax.imageio.ImageIO +import kotlin.concurrent.thread +import kotlin.math.PI + +/** + * Tracking class used for communicating with eye trackers, tracking cells with them in a sciview VR environment. + * It calls the Hedgehog analysis on the eye tracking results and communicates the results to Mastodon via + * [linkCreationCallback], which is called on every spine graph vertex that is extracted, and + * [finalTrackCallback] which is called after all vertices of a track are iterated, giving Mastodon a chance to rebuild its tracks. + */ +class EyeTracking( + override var linkCreationCallback: ((HedgehogAnalysis.SpineGraphVertex) -> Unit)? = null, + override var finalTrackCallback: (() -> Unit)? = null, + sciview: SciView +): CellTrackingBase(sciview) { + + val pupilTracker = PupilEyeTracker( + calibrationType = PupilEyeTracker.CalibrationType.WorldSpace, + port = System.getProperty("PupilPort", "50020").toInt() + ) + val calibrationTarget = Icosphere(0.02f, 2) + val laser = Cylinder(0.005f, 0.2f, 10) + + val confidenceThreshold = 0.60f + + private lateinit var lightTetrahedron: List + private lateinit var debugBoard: TextBoard + + fun run() { + + sciview.toggleVRRendering() + cellTrackingActive = true + logger.info("VR mode has been toggled") + hmd = sciview.hub.getWorkingHMD() as? OpenVRHMD ?: throw IllegalStateException("Could not find headset") + sessionId = "BionicTracking-generated-${SystemHelpers.Companion.formatDateTime()}" + sessionDirectory = Files.createDirectory(Paths.get(System.getProperty("user.home"), "Desktop", sessionId)) + + referenceTarget.visible = false + referenceTarget.ifMaterial{ + roughness = 1.0f + metallic = 0.0f + diffuse = Vector3f(0.8f, 0.8f, 0.8f) + } + referenceTarget.name = "Reference Target" + sciview.camera?.addChild(referenceTarget) + + calibrationTarget.visible = false + calibrationTarget.material { + roughness = 1.0f + metallic = 0.0f + diffuse = Vector3f(1.0f, 1.0f, 1.0f) + } + calibrationTarget.name = "Calibration Target" + sciview.camera?.addChild(calibrationTarget) + + laser.visible = false + laser.ifMaterial{diffuse = Vector3f(1.0f, 1.0f, 1.0f) } + laser.name = "Laser" + sciview.addNode(laser) + + val shell = Box(Vector3f(20.0f, 20.0f, 20.0f), insideNormals = true) + shell.ifMaterial{ + cullingMode = Material.CullingMode.Front + diffuse = Vector3f(0.4f, 0.4f, 0.4f) + } + + shell.spatial().position = Vector3f(0.0f, 0.0f, 0.0f) + shell.name = "Shell" + sciview.addNode(shell) + + val volnodes = sciview.findNodes { node -> Volume::class.java.isAssignableFrom(node.javaClass) } + + val v = (volnodes.firstOrNull() as? Volume) + if(v == null) { + logger.warn("No volume found, bailing") + return + } else { + logger.info("found ${volnodes.size} volume nodes. Using the first one: ${volnodes.first()}") + volume = v + } + volume.visible = false + + val bb = BoundingGrid() + bb.node = volume + bb.visible = false + + sciview.addNode(hedgehogs) + + val eyeFrames = Mesh("eyeFrames") + val left = Box(Vector3f(1.0f, 1.0f, 0.001f)) + val right = Box(Vector3f(1.0f, 1.0f, 0.001f)) + left.spatial().position = Vector3f(-1.0f, 1.5f, 0.0f) + left.spatial().rotation = left.spatial().rotation.rotationZ(PI.toFloat()) + right.spatial().position = Vector3f(1.0f, 1.5f, 0.0f) + eyeFrames.addChild(left) + eyeFrames.addChild(right) + + sciview.addNode(eyeFrames) + + val pupilFrameLimit = 20 + var lastFrame = System.nanoTime() + + pupilTracker.subscribeFrames { eye, texture -> + if(System.nanoTime() - lastFrame < pupilFrameLimit*10e5) { + return@subscribeFrames + } + + val node = if(eye == 1) { + left + } else { + right + } + + val stream = ByteArrayInputStream(texture) + val image = ImageIO.read(stream) + val data = (image.raster.dataBuffer as DataBufferByte).data + + node.ifMaterial { + textures["diffuse"] = Texture( + Vector3i(image.width, image.height, 1), + 3, + UnsignedByteType(), + BufferUtils.Companion.allocateByteAndPut(data) + ) + } + + lastFrame = System.nanoTime() + } + + // TODO: Replace with cam.showMessage() + debugBoard = TextBoard() + debugBoard.name = "debugBoard" + debugBoard.spatial().scale = Vector3f(0.05f, 0.05f, 0.05f) + debugBoard.spatial().position = Vector3f(0.0f, -0.3f, -0.9f) + debugBoard.text = "" + debugBoard.visible = false + sciview.camera?.addChild(debugBoard) + + lightTetrahedron = Light.Companion.createLightTetrahedron( + Vector3f(0.0f, 0.0f, 0.0f), + spread = 5.0f, + radius = 15.0f, + intensity = 5.0f + ) + lightTetrahedron.forEach { sciview.addNode(it) } + + thread { + logger.info("Adding onDeviceConnect handlers") + hmd.events.onDeviceConnect.add { hmd, device, timestamp -> + logger.info("onDeviceConnect called, cam=${sciview.camera}") + if (device.type == TrackedDeviceType.Controller) { + logger.info("Got device ${device.name} at $timestamp") + device.model?.let { hmd.attachToNode(device, it, sciview.camera) } + } + } + } + thread { + logger.info("started thread for inputSetup") + inputSetup() + setupCalibration() + } + + launchUpdaterThread() + } + + private fun setupCalibration( + keybindingCalibration: Pair = (TrackerRole.RightHand to OpenVRHMD.OpenVRButton.Menu), + keybindingTracking: Pair = (TrackerRole.RightHand to OpenVRHMD.OpenVRButton.Trigger) + ) { + val startCalibration = ClickBehaviour { _, _ -> + thread { + val cam = sciview.camera as? DetachedHeadCamera ?: return@thread + pupilTracker.gazeConfidenceThreshold = confidenceThreshold + if (!pupilTracker.isCalibrated) { + logger.info("pupil is currently uncalibrated") + pupilTracker.onCalibrationInProgress = { + cam.showMessage( + "Crunching equations ...", + distance = 2f, + size = 0.2f, + messageColor = Vector4f(1.0f, 0.8f, 0.0f, 1.0f), + duration = 15000, + centered = true + ) + } + + pupilTracker.onCalibrationFailed = { + cam.showMessage( + "Calibration failed.", + distance = 2f, + size = 0.2f, + messageColor = Vector4f(1.0f, 0.0f, 0.0f, 1.0f), + centered = true + ) + } + + pupilTracker.onCalibrationSuccess = { + cam.showMessage( + "Calibration succeeded!", + distance = 2f, + size = 0.2f, + messageColor = Vector4f(0.0f, 1.0f, 0.0f, 1.0f), + centered = true + ) + + for (i in 0 until 20) { + referenceTarget.ifMaterial { diffuse = Vector3f(0.0f, 1.0f, 0.0f) } + Thread.sleep(100) + referenceTarget.ifMaterial { diffuse = Vector3f(0.8f, 0.8f, 0.8f) } + Thread.sleep(30) + } + + hmd.removeBehaviour("start_calibration") + hmd.removeKeyBinding("start_calibration") + + val toggleTracking = ClickBehaviour { _, _ -> + if (tracking) { + logger.info("deactivating tracking...") + referenceTarget.ifMaterial { diffuse = Vector3f(0.5f, 0.5f, 0.5f) } + cam.showMessage("Tracking deactivated.", distance = 2f, size = 0.2f, centered = true) + dumpHedgehog() + } else { + logger.info("activating tracking...") + addHedgehog() + referenceTarget.ifMaterial { diffuse = Vector3f(1.0f, 0.0f, 0.0f) } + cam.showMessage("Tracking active.", distance = 2f, size = 0.2f, centered = true) + } + tracking = !tracking + } + hmd.addBehaviour("toggle_tracking", toggleTracking) + hmd.addKeyBinding("toggle_tracking", keybindingTracking.first, keybindingTracking.second) + + volume.visible = true + playing = true + } + + pupilTracker.unsubscribeFrames() + sciview.deleteNode(sciview.find("eyeFrames")) + + logger.info("Starting eye tracker calibration") + cam.showMessage( + "Follow the white rabbit.", + distance = 2f, + size = 0.2f, + duration = 1500, + centered = true + ) + pupilTracker.calibrate( + cam, hmd, + generateReferenceData = true, + calibrationTarget = calibrationTarget + ) + + pupilTracker.onGazeReceived = when (pupilTracker.calibrationType) { + + PupilEyeTracker.CalibrationType.WorldSpace -> { gaze -> + if (gaze.confidence > confidenceThreshold) { + val p = gaze.gazePoint() + referenceTarget.visible = true + // Pupil has mm units, so we divide by 1000 here to get to scenery units + referenceTarget.spatial().position = p + (cam.children.find { it.name == "debugBoard" } as? TextBoard)?.text = + "${String.format("%.2f", p.x())}, ${ + String.format( + "%.2f", + p.y() + ) + }, ${String.format("%.2f", p.z())}" + + val headCenter = cam.spatial().viewportToWorld(Vector2f(0.0f, 0.0f)) + val pointWorld = Matrix4f(cam.spatial().world).transform(p.xyzw()).xyz() + val direction = (pointWorld - headCenter).normalize() + + if (tracking) { + addSpine( + headCenter, + direction, + volume, + gaze.confidence, + volume.viewerState.currentTimepoint + ) + } + } + } + } + logger.info("Calibration routine done.") + } + } + } + hmd.addBehaviour("start_calibration", startCalibration) + hmd.addKeyBinding("start_calibration", keybindingCalibration.first, keybindingCalibration.second) + } + + /** Toggles the VR rendering off, cleans up eyetracking-related scene objects and removes the light tetrahedron + * that was created for the calibration routine. */ + override fun stop() { + pupilTracker.unsubscribeFrames() + cellTrackingActive = false + logger.info("Stopped volume and hedgehog updater thread.") + lightTetrahedron.forEach { sciview.deleteNode(it) } + sciview.deleteNode(sciview.find("Shell")) + sciview.deleteNode(sciview.find("eyeFrames")) + listOf(referenceTarget, calibrationTarget, laser, debugBoard, hedgehogs).forEach { + sciview.deleteNode(it) + } + logger.info("Successfully cleaned up eye tracking environemt.") + sciview.toggleVRRendering() + logger.info("Shut down eye tracking environment and disabled VR.") + } + +} diff --git a/src/main/kotlin/sc/iview/commands/analysis/HedgehogAnalysis.kt b/src/main/kotlin/sc/iview/commands/analysis/HedgehogAnalysis.kt new file mode 100644 index 000000000..5364adc42 --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/HedgehogAnalysis.kt @@ -0,0 +1,435 @@ +package sc.iview.commands.analysis + +import org.joml.Vector3f +import org.joml.Matrix4f +import org.joml.Quaternionf +import graphics.scenery.utils.extensions.* +import graphics.scenery.utils.lazyLogger +import org.slf4j.LoggerFactory +import sc.iview.commands.analysis.SpineMetadata +import java.io.File +import kotlin.math.sqrt + +/** + * Performs analysis over a collection of eye-tracking spines (aka hedgehog). Extracts a list of local maxima from + * the sampled volume, removes statistical outliers and performs a graph optimization over the remaining maxima to + * extract the likeliest path of the cell. The companion object contains methods to load CSV files. + * @author Ulrik Günther + */ +class HedgehogAnalysis(val spines: List, val localToWorld: Matrix4f) { + + private val logger by lazyLogger() + + val timepoints = LinkedHashMap>() + + var avgConfidence = 0.0f + private set + var totalSampleCount = 0 + private set + + /** Data class for collecting track points, consisting of positions and a [SpineGraphVertex], and the averaged + * confidences of all spines. Returned by [kotlin.run]. */ + data class Track( + val points: List>, + val confidence: Float + ) + + init { + logger.info("Starting analysis with ${spines.size} spines") + + spines.forEach { spine -> + val timepoint = spine.timepoint + val current = timepoints[timepoint] + + if(current == null) { + timepoints[timepoint] = arrayListOf(spine) + } else { + current.add(spine) + } + + avgConfidence += spine.confidence + totalSampleCount++ + } + + avgConfidence /= totalSampleCount + } + + /** + * From a [list] of Floats, return both the index of local maxima, and their value, + * packaged nicely as a Pair + */ + private fun localMaxima(list: List): List> { + return list.windowed(3, 1).mapIndexed { index, l -> + val left = l[0] + val center = l[1] + val right = l[2] + + // we have a match at center + if (left < center && center > right) { + index * 1 + 1 to center + } else { + null + } + }.filterNotNull() + } + + /** Cell positions extracted from gaze analysis are collected in this data class together with other information + * such as the volume [value] at this point, and the [previous] and [next] vertices. */ + data class SpineGraphVertex(val timepoint: Int, + val position: Vector3f, + val worldPosition: Vector3f, + val index: Int, + val value: Float, + val metadata : SpineMetadata, + var previous: SpineGraphVertex? = null, + var next: SpineGraphVertex? = null) { + + fun distance(): Float { + val n = next + return if(n != null) { + val t = (n.worldPosition - this.worldPosition) + sqrt(t.x*t.x + t.y*t.y + t.z*t.z) + } else { + 0.0f + } + } + + fun drop() { + previous?.next = next + next?.previous = previous + } + + override fun toString() : String { + return "SpineGraphVertex for t=$timepoint, pos=$position,index=$index, worldPos=$worldPosition, value=$value" + } + } + + fun Iterable.stddev() = sqrt((this.map { (it - this.average()) * (it - this.average()) }.sum() / this.count())) + + fun Vector3f.toQuaternionf(forward: Vector3f = Vector3f(0.0f, 0.0f, -1.0f)): Quaternionf { + val cross = forward.cross(this) + val q = Quaternionf(cross.x(), cross.y(), cross.z(), this.dot(forward)) + + val x = sqrt((q.w + sqrt(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w)) / 2.0f) + + return Quaternionf(q.x/(2.0f * x), q.y/(2.0f * x), q.z/(2.0f * x), x) + } + + data class VertexWithDistance(val vertex: SpineGraphVertex, val distance: Float) + + fun run(): Track? { + + val startingThreshold = 0.002f + val localMaxThreshold = 0.001f + val zscoreThreshold = 2.0f + val removeTooFarThreshold = 5.0f + + if(timepoints.isEmpty()) { + return null + } + + + //step1: find the startingPoint by using startingthreshold + val startingPoint = timepoints.entries.firstOrNull { entry -> + entry.value.any { metadata -> metadata.samples.filterNotNull().any { it > startingThreshold } } + } ?: return null + + logger.info("Starting point is ${startingPoint.key}/${timepoints.size} (threshold=$startingThreshold)") + + // filter timepoints, remove all before the starting point + timepoints.filter { it.key > startingPoint.key } + .forEach { timepoints.remove(it.key) } + + logger.info("${timepoints.size} timepoints left") + + fun gaussSmoothing(samples: List, iterations: Int): List { + var smoothed = samples.toList() + val kernel = listOf(0.25f, 0.5f, 0.25f) + for (i in 0 until iterations) { + val newSmoothed = ArrayList(smoothed.size) + // Handle the first element + newSmoothed.add(smoothed[0] * 0.75f + smoothed[1] * 0.25f) + // Apply smoothing to the middle elements + for (j in 1 until smoothed.size - 1) { + val value = kernel[0] * smoothed[j-1] + kernel[1] * smoothed[j] + kernel[2] * smoothed[j+1] + newSmoothed.add(value) + } + // Handle the last element + newSmoothed.add(smoothed[smoothed.size - 2] * 0.25f + smoothed[smoothed.size - 1] * 0.75f) + + smoothed = newSmoothed + } + return smoothed + } + + //step2: find the maxIndices along the spine + // this will be a list of lists, where each entry in the first-level list + // corresponds to a time point, which then contains a list of vertices within that timepoint. + val candidates: List> = timepoints.map { tp -> + val vs = tp.value.mapIndexedNotNull { i, spine -> + // determine local maxima (and their indices) along the spine, aka, actual things the user might have + // seen when looking into the direction of the spine + val maxIndices = localMaxima(spine.samples.filterNotNull()) + logger.debug("Local maxima at ${tp.key}/$i are: ${maxIndices.joinToString(",")}") + + // if there actually are local maxima, generate a graph vertex for them with all the necessary metadata + if(maxIndices.isNotEmpty()) { + //maxIndices. +// filter the maxIndices which are too far away, which can be removed + //filter { it.first <1200}. + maxIndices.map { index -> + logger.debug("Generating vertex at index $index") + // get the position of the current index along the spine + val position = spine.samplePosList[index.first] + val worldPosition = localToWorld.transform((Vector3f(position)).xyzw()).xyz() + SpineGraphVertex(tp.key, + position, + worldPosition, + index.first, + index.second, + spine) + + } + } else { + null + } + } + vs + }.flatten() + + logger.info("SpineGraphVertices extracted") + + // step3: connect localMaximal points between 2 candidate spines according to the shortest path principle + // get the initial vertex, this one is assumed to always be in front, and have a local maximum - aka, what + // the user looks at first is assumed to be the actual cell they want to track + val initial = candidates.first().first { it.value > startingThreshold } + var current = initial + var shortestPath = candidates.drop(1).mapIndexedNotNull { time, vs -> + // calculate world-space distances between current point, and all candidate + // vertices, sorting them by distance + val vertices = vs + .filter { it.value > localMaxThreshold } + .map { vertex -> + val t = current.worldPosition - vertex.worldPosition + val distance = t.length() + VertexWithDistance(vertex, distance) + } + .sortedBy { it.distance } + + val closest = vertices.firstOrNull() + if(closest != null && closest.distance > 0) { + // create a linked list between current and closest vertices + current.next = closest.vertex + closest.vertex.previous = current + current = closest.vertex + current + } else { + null + } + }.toMutableList() + + // calculate average path lengths over all + val beforeCount = shortestPath.size + var avgPathLength = shortestPath.map { it.distance() }.average().toFloat() + var stdDevPathLength = shortestPath.map { it.distance() }.stddev().toFloat() + logger.info("Average path length=$avgPathLength, stddev=$stdDevPathLength") + + fun zScore(value: Float, m: Float, sd: Float) = ((value - m)/sd) + + //step4: if some path is longer than multiple average length, it should be removed + while (shortestPath.any { it.distance() >= removeTooFarThreshold * avgPathLength }) { + shortestPath = shortestPath.filter { it.distance() < removeTooFarThreshold * avgPathLength }.toMutableList() + shortestPath.windowed(3, 1, partialWindows = true).forEach { + // this reconnects the neighbors after the offending vertex has been removed + it.getOrNull(0)?.next = it.getOrNull(1) + it.getOrNull(1)?.previous = it.getOrNull(0) + it.getOrNull(1)?.next = it.getOrNull(2) + it.getOrNull(2)?.previous = it.getOrNull(1) + } + + } + + // recalculate statistics after offending vertex removal + avgPathLength = shortestPath.map { it.distance() }.average().toFloat() + stdDevPathLength = shortestPath.map { it.distance() }.stddev().toFloat() + + //step5: remove some vertices according to zscoreThreshold + var remaining = shortestPath.count { zScore(it.distance(), avgPathLength, stdDevPathLength) > zscoreThreshold } + logger.info("Iterating: ${shortestPath.size} vertices remaining, with $remaining failing z-score criterion") + while(remaining > 0) { + val outliers = shortestPath + .filter { zScore(it.distance(), avgPathLength, stdDevPathLength) > zscoreThreshold } + .map { + val idx = shortestPath.indexOf(it) + listOf(idx-1,idx,idx+1) + }.flatten() + + shortestPath = shortestPath.filterIndexed { index, _ -> index !in outliers }.toMutableList() + remaining = shortestPath.count { zScore(it.distance(), avgPathLength, stdDevPathLength) > zscoreThreshold } + + shortestPath.windowed(3, 1, partialWindows = true).forEach { + it.getOrNull(0)?.next = it.getOrNull(1) + it.getOrNull(1)?.previous = it.getOrNull(0) + it.getOrNull(1)?.next = it.getOrNull(2) + it.getOrNull(2)?.previous = it.getOrNull(1) + } + logger.info("Iterating: ${shortestPath.size} vertices remaining, with $remaining failing z-score criterion") + } + + val afterCount = shortestPath.size + logger.info("Pruned ${beforeCount - afterCount} vertices due to path length") + val singlePoints = shortestPath + .groupBy { it.timepoint } + .mapNotNull { vs -> vs.value.maxByOrNull{ it.metadata.confidence } } + .filter { + it.metadata.direction.dot(it.previous!!.metadata.direction) > 0.5f + } + + + logger.info("Returning ${singlePoints.size} points") + + + return Track(singlePoints.map { it.position to it}, avgConfidence) + } + + companion object { + private val logger by lazyLogger(System.getProperty("scenery.LogLevel", "info")) + + fun fromIncompleteCSV(csv: File, separator: String = ","): HedgehogAnalysis { + logger.info("Loading spines from incomplete CSV at ${csv.absolutePath}") + + val lines = csv.readLines() + val spines = ArrayList(lines.size) + + lines.drop(1).forEach { line -> + val tokens = line.split(separator) + val timepoint = tokens[0].toInt() + val confidence = tokens[1].toFloat() + val samples = tokens.subList(2, tokens.size - 1).map { it.toFloat() } + + val currentSpine = SpineMetadata( + timepoint, + Vector3f(0.0f), + Vector3f(0.0f), + 0.0f, + Vector3f(0.0f), + Vector3f(0.0f), + Vector3f(0.0f), + Vector3f(0.0f), + Quaternionf(), + Vector3f(0.0f), + confidence, + samples + ) + + spines.add(currentSpine) + } + + return HedgehogAnalysis(spines, Matrix4f()) + } + + private fun String.toVector3f(): Vector3f { + val array = this.replace("(", "").replace(")", "").trim().split(" ").filterNot { it == ""} + + if (array[0] == "+Inf" || array[0] == "-Inf") + return Vector3f(0.0f,0.0f,0.0f) + + return Vector3f(array[0].toFloat(),array[1].toFloat(),array[2].toFloat()) + } + + private fun String.toQuaternionf(): Quaternionf { + val array = this.replace("(", "").replace(")", "").trim().split(" ").filterNot { it == ""} + return Quaternionf(array[0].toFloat(), array[1].toFloat(), array[2].toFloat(), array[3].toFloat()) + } + fun fromCSVWithMatrix(csv: File, matrix4f: Matrix4f,separator: String = ";"): HedgehogAnalysis { + logger.info("Loading spines from complete CSV with Matrix at ${csv.absolutePath}") + + val lines = csv.readLines() + val spines = ArrayList(lines.size) + logger.info("lines number: " + lines.size) + lines.drop(1).forEach { line -> + val tokens = line.split(separator) + val timepoint = tokens[0].toInt() + val origin = tokens[1].toVector3f() + val direction = tokens[2].toVector3f() + val localEntry = tokens[3].toVector3f() + val localExit = tokens[4].toVector3f() + val localDirection = tokens[5].toVector3f() + val headPosition = tokens[6].toVector3f() + val headOrientation = tokens[7].toQuaternionf() + val position = tokens[8].toVector3f() + val confidence = tokens[9].toFloat() + val samples = tokens.subList(10, tokens.size - 1).map { it.toFloat() } + + val currentSpine = SpineMetadata( + timepoint, + origin, + direction, + 0.0f, + localEntry, + localExit, + localDirection, + headPosition, + headOrientation, + position, + confidence, + samples + ) + + spines.add(currentSpine) + } + + return HedgehogAnalysis(spines, matrix4f) + } + + fun fromCSV(csv: File, separator: String = ";"): HedgehogAnalysis { + logger.info("Loading spines from complete CSV at ${csv.absolutePath}") + + val lines = csv.readLines() + val spines = ArrayList(lines.size) + + lines.drop(1).forEach { line -> + val tokens = line.split(separator) + val timepoint = tokens[0].toInt() + val origin = tokens[1].toVector3f() + val direction = tokens[2].toVector3f() + val localEntry = tokens[3].toVector3f() + val localExit = tokens[4].toVector3f() + val localDirection = tokens[5].toVector3f() + val headPosition = tokens[6].toVector3f() + val headOrientation = tokens[7].toQuaternionf() + val position = tokens[8].toVector3f() + val confidence = tokens[9].toFloat() + val samples = tokens.subList(10, tokens.size - 1).map { it.toFloat() } + + val currentSpine = SpineMetadata( + timepoint, + origin, + direction, + 0.0f, + localEntry, + localExit, + localDirection, + headPosition, + headOrientation, + position, + confidence, + samples + ) + + spines.add(currentSpine) + } + + return HedgehogAnalysis(spines, Matrix4f()) + } + } +} + +fun main(filePath: String, args: Array) { + val logger = LoggerFactory.getLogger("HedgehogAnalysisMain") + // main should only be called for testing purposes + val file = File(filePath) + val analysis = HedgehogAnalysis.fromCSV(file) + val results = analysis.run() + logger.info("Results: \n$results") +} diff --git a/src/main/kotlin/sc/iview/commands/analysis/SpineMetadata.kt b/src/main/kotlin/sc/iview/commands/analysis/SpineMetadata.kt new file mode 100644 index 000000000..8417778c4 --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/SpineMetadata.kt @@ -0,0 +1,23 @@ +package sc.iview.commands.analysis + +import org.joml.Quaternionf +import org.joml.Vector3f + +/** + * Data class to store metadata for spines of the hedgehog. + */ +data class SpineMetadata( + val timepoint: Int, + val origin: Vector3f, + val direction: Vector3f, + val distance: Float, + val localEntry: Vector3f, + val localExit: Vector3f, + val localDirection: Vector3f, + val headPosition: Vector3f, + val headOrientation: Quaternionf, + val position: Vector3f, + val confidence: Float, + val samples: List, + val samplePosList: List = ArrayList() +) diff --git a/src/main/kotlin/sc/iview/commands/analysis/TimepointObserver.kt b/src/main/kotlin/sc/iview/commands/analysis/TimepointObserver.kt new file mode 100644 index 000000000..584cb5d8a --- /dev/null +++ b/src/main/kotlin/sc/iview/commands/analysis/TimepointObserver.kt @@ -0,0 +1,13 @@ +package sc.iview.commands.analysis + +/** + * Interface to allow subscription to timepoint updates, especially for updating sciview contents + * after a user triggered a timepoint change via controller input. + */ +interface TimepointObserver { + + /** + * Called when the timepoint was updated. + */ + fun onTimePointChanged(timepoint: Int) +} diff --git a/src/main/resources/sc/iview/commands/demo/animation/ParticleDemo.vert b/src/main/resources/sc/iview/commands/demo/animation/ParticleDemo.vert index 4a3d2b12f..1f612d645 100644 --- a/src/main/resources/sc/iview/commands/demo/animation/ParticleDemo.vert +++ b/src/main/resources/sc/iview/commands/demo/animation/ParticleDemo.vert @@ -58,7 +58,7 @@ mat4 mv; mv = (vrParameters.stereoEnabled ^ 1) * ViewMatrices[0] * iModelMatrix + (vrParameters.stereoEnabled * ViewMatrices[currentEye.eye] * iModelMatrix); projectionMatrix = (vrParameters.stereoEnabled ^ 1) * ProjectionMatrix + vrParameters.stereoEnabled * vrParameters.projectionMatrices[currentEye.eye]; - if(ubo.isBillboard > 0) { + if(ubo.isBillboard == 1) { mv[0][0] = 1.0f; mv[0][1] = .0f; mv[0][2] = .0f;