Generate gradients using sobel operator (#280)

* Sobel gradients * Apply bilateral image filtering to data values before calculating gradients

Author: mlavik1

Assets/Editor/ImportSettingsEditorWindow.cs

@@ -43,6 +43,14 @@ private void OnGUI()
  PlayerPrefs.SetInt("ClampHounsfield", clampHounsfield ? 1 : 0);
  EditorGUILayout.EndHorizontal();
 
+ EditorGUILayout.BeginHorizontal();
+ EditorGUILayout.LabelField("Default gradient computation method");
+ GradientType gradientType = GradientType.CentralDifference;
+ Enum.TryParse(PlayerPrefs.GetString("DefaultGradientType"), out gradientType);
+ gradientType = (GradientType)EditorGUILayout.EnumPopup(gradientType, "");
+ PlayerPrefs.SetString("DefaultGradientType", gradientType.ToString());
+ EditorGUILayout.EndHorizontal();
+
  EditorGUILayout.Space();
  EditorGUILayout.Space();
 

Assets/Editor/VolumeRenderedObjectCustomInspector.cs

@@ -112,6 +112,14 @@ public override void OnInspectorGUI()
  else
  volrendObj.SetLightingEnabled(false);
 
+ GradientType oldGradientType = volrendObj.GetGradientType();
+ GradientType newGradientType = (GradientType)EditorGUILayout.EnumPopup("Gradient", oldGradientType);
+
+ if (newGradientType != oldGradientType)
+ {
+ volrendObj.SetGradientTypeAsync(newGradientType, new ProgressHandler(this));
+ }
+
  if (volrendObj.GetLightingEnabled() || volrendObj.GetRenderMode() == RenderMode.IsosurfaceRendering)
  {
  LightSource oldLightSource = volrendObj.GetLightSource();

Assets/Scripts/Gradients/CentralDifferenceGradientComputator.cs

@@ -0,0 +1,27 @@
+using itk.simple;
+using openDicom.Image;
+using System;
+using System.Runtime.InteropServices;
+using UnityEngine;
+
+namespace UnityVolumeRendering
+{
+ public class CentralDifferenceGradientComputator : GradientComputator
+ {
+ public CentralDifferenceGradientComputator(VolumeDataset dataset, bool smootheDataValues) : base(dataset, smootheDataValues)
+ {
+ }
+
+ public override Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange)
+ {
+ float x1 = data[Math.Min(x + 1, dimX - 1) + y * dimX + z * (dimX * dimY)] - minValue;
+ float x2 = data[Math.Max(x - 1, 0) + y * dimX + z * (dimX * dimY)] - minValue;
+ float y1 = data[x + Math.Min(y + 1, dimY - 1) * dimX + z * (dimX * dimY)] - minValue;
+ float y2 = data[x + Math.Max(y - 1, 0) * dimX + z * (dimX * dimY)] - minValue;
+ float z1 = data[x + y * dimX + Math.Min(z + 1, dimZ - 1) * (dimX * dimY)] - minValue;
+ float z2 = data[x + y * dimX + Math.Max(z - 1, 0) * (dimX * dimY)] - minValue;
+
+ return new Vector3((x2 - x1) / maxRange, (y2 - y1) / maxRange, (z2 - z1) / maxRange);
+ }
+ }
+}

Assets/Scripts/Gradients/GradientComputator.cs

@@ -0,0 +1,71 @@
+using itk.simple;
+using System;
+using System.Runtime.InteropServices;
+using UnityEngine;
+
+namespace UnityVolumeRendering
+{
+ public abstract class GradientComputator
+ {
+ protected float[] data;
+ protected int dimX, dimY, dimZ;
+
+ public GradientComputator(VolumeDataset dataset, bool smootheDataValues)
+ {
+ this.data = dataset.data;
+ this.dimX = dataset.dimX;
+ this.dimY = dataset.dimY;
+ this.dimZ = dataset.dimZ;
+ 
+ if (smootheDataValues)
+ {
+#if UVR_USE_SIMPLEITK
+ Image image = new Image((uint)dimX, (uint)dimY, (uint)dimZ, PixelIDValueEnum.sitkFloat32);
+
+ for (uint z = 0; z < dimZ; z++)
+ {
+ for (uint y = 0; y < dimY; y++)
+ {
+ for (uint x = 0; x < dimX; x++)
+ {
+ float value = data[x + y * dimX + z * (dimX * dimY)];
+ image.SetPixelAsFloat(new VectorUInt32() { x, y, z }, value);
+ }
+ }
+ }
+
+ BilateralImageFilter filter = new BilateralImageFilter();
+ image = filter.Execute(image);
+
+ this.data = new float[data.Length];
+ IntPtr imgBuffer = image.GetBufferAsFloat();
+ Marshal.Copy(imgBuffer, data, 0, data.Length);
+#else
+ Debug.LogWarning("SimpleITK is required to generate smooth gradients.");
+#endif
+ }
+ }
+
+ public abstract Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange);
+ }
+
+ public class GradientComputatorFactory
+ {
+ public static GradientComputator CreateGradientComputator(VolumeDataset dataset, GradientType gradientType)
+ {
+ switch (gradientType)
+ {
+ case GradientType.CentralDifference:
+ return new CentralDifferenceGradientComputator(dataset, false);
+ case GradientType.SmoothedCentralDifference:
+ return new CentralDifferenceGradientComputator(dataset, true);
+ case GradientType.Sobel:
+ return new SobelGradientComputator(dataset, false);
+ case GradientType.SmoothedSobel:
+ return new SobelGradientComputator(dataset, true);
+ default:
+ throw new NotImplementedException();
+ }
+ }
+ }
+}

Assets/Scripts/Gradients/GradientType.cs

@@ -0,0 +1,23 @@
+using System;
+using UnityEngine;
+
+namespace UnityVolumeRendering
+{
+ public enum GradientType
+ {
+ CentralDifference,
+ SmoothedCentralDifference,
+ Sobel,
+ SmoothedSobel
+ }
+
+ public class GradientTypeUtils
+ {
+ public static GradientType GetDefaultGradientType()
+ {
+ GradientType gradientType = GradientType.CentralDifference;
+ Enum.TryParse(PlayerPrefs.GetString("DefaultGradientType"), out gradientType);
+ return gradientType;
+ }
+ }
+}

Assets/Scripts/Gradients/SobelGradientComputator.cs

@@ -0,0 +1,73 @@
+using System;
+using System.Drawing;
+using UnityEngine;
+
+namespace UnityVolumeRendering
+{
+ public class SobelGradientComputator : GradientComputator
+ {
+ public SobelGradientComputator(VolumeDataset dataset, bool smootheDataValues) : base(dataset, smootheDataValues)
+ {
+ }
+
+ private static readonly float[,,] kernelx = new float[,,]
+ {
+ { { -1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } },
+ { { -2, 0, 2 }, { -4, 0, 4 }, { -2, 0, 2 } },
+ { { -1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } }
+ };
+
+ private static readonly float[,,] kernely = new float[,,]
+ {
+ { { -1, -2, -1 }, { 0, 0, 0 }, { 1, 2, 1 } },
+ { { -2, -4, -2 }, { 0, 0, 0 }, { 2, 4, 2 } },
+ { { -1, -2, -1 }, { 0, 0, 0 }, { 1, 2, 1 } }
+ };
+
+ private static readonly float[,,] kernelz = new float[,,]
+ {
+ { { -1, -2, -1 }, { -2, -4, -2 }, { -1, -2, -1 } },
+ { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } },
+ { { 1, 2, 1 }, { 2, 4, 2 }, { 1, 2, 1 } }
+ };
+
+ private float GetData(int x, int y, int z)
+ {
+ return data[x + y * dimX + z * (dimX * dimY)];
+ }
+
+ private Vector3 ConvolveWithKernels(int x, int y, int z)
+ {
+ Vector3 result = Vector3.zero;
+ for (int iz = 0; iz <= 2; iz++)
+ {
+ for (int iy = 0; iy <= 2; iy++)
+ {
+ for (int ix = 0; ix <= 2; ix++)
+ {
+ float dataValue = GetData(x + ix - 1, y + iy - 1, z + iz - 1);
+ result.x += kernelx[iz, iy, ix] * dataValue;
+ result.y += kernely[iz, iy, ix] * dataValue;
+ result.z += kernelz[iz, iy, ix] * dataValue;
+ }
+ }
+ }
+ return result;
+ }
+
+ public override Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange)
+ {
+ // TODO
+ if (x < 2 || y < 2 || z < 2 || x > dimX - 3 || y > dimY - 3 || z > dimZ - 3)
+ {
+ return Vector3.zero;
+ }
+
+ Vector3 gradient = ConvolveWithKernels(x, y, z);
+
+ float divident = maxRange * 12;
+
+ return new Vector3(gradient.x / divident, gradient.y / divident, gradient.z / divident);
+ }
+ }
+}

Assets/Scripts/VolumeData/VolumeDataset.cs

@@ -113,7 +113,7 @@ public Texture3D GetGradientTexture()
  {
  if (gradientTexture == null)
  {
- gradientTexture = AsyncHelper.RunSync<Texture3D>(() => CreateGradientTextureInternalAsync(NullProgressHandler.instance));
+ gradientTexture = AsyncHelper.RunSync<Texture3D>(() => CreateGradientTextureInternalAsync(GradientTypeUtils.GetDefaultGradientType(), NullProgressHandler.instance));
  return gradientTexture;
  }
  else
@@ -122,6 +122,29 @@ public Texture3D GetGradientTexture()
  }
  }
 
+ public async Task<Texture3D> RegenerateGradientTextureAsync(GradientType gradientType, IProgressHandler progressHandler = null)
+ {
+ await createGradientTextureLock.WaitAsync();
+ try
+ {
+ if (progressHandler == null)
+ progressHandler = new NullProgressHandler();
+ try
+ {
+ gradientTexture = await CreateGradientTextureInternalAsync(gradientType, progressHandler != null ? progressHandler : NullProgressHandler.instance);
+ }
+ catch (System.Exception exception)
+ {
+ Debug.LogException(exception);
+ }
+ }
+ finally
+ {
+ createGradientTextureLock.Release();
+ }
+ return gradientTexture;
+ }
+
  /// <summary>
  /// Gets the gradient texture, containing the gradient values (direction of change) of the dataset.
  /// Will create the gradient texture if it does not exist, without blocking the main thread.
@@ -132,17 +155,7 @@ public async Task<Texture3D> GetGradientTextureAsync(IProgressHandler progressHa
  {
  if (gradientTexture == null)
  {
- await createGradientTextureLock.WaitAsync();
- try
- {
- if (progressHandler == null)
- progressHandler = new NullProgressHandler();
- gradientTexture = await CreateGradientTextureInternalAsync(progressHandler != null ? progressHandler : NullProgressHandler.instance);
- }
- finally
- {
- createGradientTextureLock.Release();
- }
+ gradientTexture = await RegenerateGradientTextureAsync(GradientTypeUtils.GetDefaultGradientType(), progressHandler);
  }
  return gradientTexture;
  }
@@ -332,7 +345,7 @@ await Task.Run(() => {
  return texture;
  }
 
- private async Task<Texture3D> CreateGradientTextureInternalAsync(IProgressHandler progressHandler)
+ private async Task<Texture3D> CreateGradientTextureInternalAsync(GradientType gradientType, IProgressHandler progressHandler)
  {
  Debug.Log("Async gradient generation. Hold on.");
 
@@ -364,6 +377,8 @@ await Task.Run(() => {
  Texture3D textureTmp = new Texture3D(dimX, dimY, dimZ, texformat, false);
  textureTmp.wrapMode = TextureWrapMode.Clamp;
 
+ GradientComputator gradientComputator = GradientComputatorFactory.CreateGradientComputator(this, gradientType);
+
  for (int x = 0; x < dimX; x++)
  {
  progressHandler.ReportProgress(x, dimX, "Calculating gradients for slice");
@@ -372,7 +387,7 @@ await Task.Run(() => {
  for (int z = 0; z < dimZ; z++)
  {
  int iData = x + y * dimX + z * (dimX * dimY);
- Vector3 grad = GetGrad(x, y, z, minValue, maxRange);
+ Vector3 grad = gradientComputator.ComputeGradient(x, y, z, minValue, maxRange);
 
  textureTmp.SetPixel(x, y, z, new Color(grad.x, grad.y, grad.z, (float)(data[iData] - minValue) / maxRange));
  }
@@ -390,6 +405,8 @@ await Task.Run(() => {
 
  progressHandler.StartStage(0.6f, "Creating gradient texture");
  await Task.Run(() => {
+ GradientComputator gradientComputator = GradientComputatorFactory.CreateGradientComputator(this, gradientType);
+
  for (int z = 0; z < dimZ; z++)
  {
  progressHandler.ReportProgress(z, dimZ, "Calculating gradients for slice");
@@ -398,7 +415,7 @@ await Task.Run(() => {
  for (int x = 0; x < dimX; x++)
  {
  int iData = x + y * dimX + z * (dimX * dimY);
- Vector3 grad = GetGrad(x, y, z, minValue, maxRange);
+ Vector3 grad = gradientComputator.ComputeGradient(x, y, z, minValue, maxRange);
 
  cols[iData] = new Color(grad.x, grad.y, grad.z, (float)(data[iData] - minValue) / maxRange);
  }
@@ -418,17 +435,6 @@ await Task.Run(() => {
  return texture;
 
  }
- public Vector3 GetGrad(int x, int y, int z, float minValue, float maxRange)
- {
- float x1 = data[Math.Min(x + 1, dimX - 1) + y * dimX + z * (dimX * dimY)] - minValue;
- float x2 = data[Math.Max(x - 1, 0) + y * dimX + z * (dimX * dimY)] - minValue;
- float y1 = data[x + Math.Min(y + 1, dimY - 1) * dimX + z * (dimX * dimY)] - minValue;
- float y2 = data[x + Math.Max(y - 1, 0) * dimX + z * (dimX * dimY)] - minValue;
- float z1 = data[x + y * dimX + Math.Min(z + 1, dimZ - 1) * (dimX * dimY)] - minValue;
- float z2 = data[x + y * dimX + Math.Max(z - 1, 0) * (dimX * dimY)] - minValue;
-
- return new Vector3((x2 - x1) / maxRange, (y2 - y1) / maxRange, (z2 - z1) / maxRange);
- }
 
  public float GetAvgerageVoxelValues(int x, int y, int z)
  {

Assets/Scripts/VolumeObject/VolumeObjectFactory.cs

@@ -10,6 +10,7 @@ public static VolumeRenderedObject CreateObject(VolumeDataset dataset)
  {
  GameObject outerObject = new GameObject("VolumeRenderedObject_" + dataset.datasetName);
  VolumeRenderedObject volObj = outerObject.AddComponent<VolumeRenderedObject>();
+ volObj.SetGradientType(GradientTypeUtils.GetDefaultGradientType());
 
  GameObject meshContainer = GameObject.Instantiate((GameObject)Resources.Load("VolumeContainer"));
  volObj.volumeContainerObject = meshContainer;
@@ -25,6 +26,7 @@ public static async Task<VolumeRenderedObject> CreateObjectAsync(VolumeDataset d
  {
  GameObject outerObject = new GameObject("VolumeRenderedObject_" + dataset.datasetName);
  VolumeRenderedObject volObj = outerObject.AddComponent<VolumeRenderedObject>();
+ volObj.SetGradientType(GradientTypeUtils.GetDefaultGradientType());
 
  GameObject meshContainer = GameObject.Instantiate((GameObject)Resources.Load("VolumeContainer"));
  volObj.volumeContainerObject = meshContainer;