Skip to content

Commit fbeee24

Browse files
mlavik1mlavik1
authored
Generate gradients using sobel operator (#280)
* Sobel gradients * Apply bilateral image filtering to data values before calculating gradients
1 parent 62f9fea commit fbeee24

9 files changed

+276
-27
lines changed

‎Assets/Editor/ImportSettingsEditorWindow.cs

+8
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,14 @@ private void OnGUI()
4343
PlayerPrefs.SetInt("ClampHounsfield", clampHounsfield ? 1 : 0);
4444
EditorGUILayout.EndHorizontal();
4545

46+
EditorGUILayout.BeginHorizontal();
47+
EditorGUILayout.LabelField("Default gradient computation method");
48+
GradientType gradientType = GradientType.CentralDifference;
49+
Enum.TryParse(PlayerPrefs.GetString("DefaultGradientType"), out gradientType);
50+
gradientType = (GradientType)EditorGUILayout.EnumPopup(gradientType, "");
51+
PlayerPrefs.SetString("DefaultGradientType", gradientType.ToString());
52+
EditorGUILayout.EndHorizontal();
53+
4654
EditorGUILayout.Space();
4755
EditorGUILayout.Space();
4856

‎Assets/Editor/VolumeRenderedObjectCustomInspector.cs

+8
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,14 @@ public override void OnInspectorGUI()
112112
else
113113
volrendObj.SetLightingEnabled(false);
114114

115+
GradientType oldGradientType = volrendObj.GetGradientType();
116+
GradientType newGradientType = (GradientType)EditorGUILayout.EnumPopup("Gradient", oldGradientType);
117+
118+
if (newGradientType != oldGradientType)
119+
{
120+
volrendObj.SetGradientTypeAsync(newGradientType, new ProgressHandler(this));
121+
}
122+
115123
if (volrendObj.GetLightingEnabled() || volrendObj.GetRenderMode() == RenderMode.IsosurfaceRendering)
116124
{
117125
LightSource oldLightSource = volrendObj.GetLightSource();

‎Assets/Scripts/Gradients/CentralDifferenceGradientComputator.cs

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
using itk.simple;
2+
using openDicom.Image;
3+
using System;
4+
using System.Runtime.InteropServices;
5+
using UnityEngine;
6+
7+
namespace UnityVolumeRendering
8+
{
9+
public class CentralDifferenceGradientComputator : GradientComputator
10+
{
11+
public CentralDifferenceGradientComputator(VolumeDataset dataset, bool smootheDataValues) : base(dataset, smootheDataValues)
12+
{
13+
}
14+
15+
public override Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange)
16+
{
17+
float x1 = data[Math.Min(x + 1, dimX - 1) + y * dimX + z * (dimX * dimY)] - minValue;
18+
float x2 = data[Math.Max(x - 1, 0) + y * dimX + z * (dimX * dimY)] - minValue;
19+
float y1 = data[x + Math.Min(y + 1, dimY - 1) * dimX + z * (dimX * dimY)] - minValue;
20+
float y2 = data[x + Math.Max(y - 1, 0) * dimX + z * (dimX * dimY)] - minValue;
21+
float z1 = data[x + y * dimX + Math.Min(z + 1, dimZ - 1) * (dimX * dimY)] - minValue;
22+
float z2 = data[x + y * dimX + Math.Max(z - 1, 0) * (dimX * dimY)] - minValue;
23+
24+
return new Vector3((x2 - x1) / maxRange, (y2 - y1) / maxRange, (z2 - z1) / maxRange);
25+
}
26+
}
27+
}

‎Assets/Scripts/Gradients/GradientComputator.cs

+71
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
using itk.simple;
2+
using System;
3+
using System.Runtime.InteropServices;
4+
using UnityEngine;
5+
6+
namespace UnityVolumeRendering
7+
{
8+
public abstract class GradientComputator
9+
{
10+
protected float[] data;
11+
protected int dimX, dimY, dimZ;
12+
13+
public GradientComputator(VolumeDataset dataset, bool smootheDataValues)
14+
{
15+
this.data = dataset.data;
16+
this.dimX = dataset.dimX;
17+
this.dimY = dataset.dimY;
18+
this.dimZ = dataset.dimZ;
19+
20+
if (smootheDataValues)
21+
{
22+
#if UVR_USE_SIMPLEITK
23+
Image image = new Image((uint)dimX, (uint)dimY, (uint)dimZ, PixelIDValueEnum.sitkFloat32);
24+
25+
for (uint z = 0; z < dimZ; z++)
26+
{
27+
for (uint y = 0; y < dimY; y++)
28+
{
29+
for (uint x = 0; x < dimX; x++)
30+
{
31+
float value = data[x + y * dimX + z * (dimX * dimY)];
32+
image.SetPixelAsFloat(new VectorUInt32() { x, y, z }, value);
33+
}
34+
}
35+
}
36+
37+
BilateralImageFilter filter = new BilateralImageFilter();
38+
image = filter.Execute(image);
39+
40+
this.data = new float[data.Length];
41+
IntPtr imgBuffer = image.GetBufferAsFloat();
42+
Marshal.Copy(imgBuffer, data, 0, data.Length);
43+
#else
44+
Debug.LogWarning("SimpleITK is required to generate smooth gradients.");
45+
#endif
46+
}
47+
}
48+
49+
public abstract Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange);
50+
}
51+
52+
public class GradientComputatorFactory
53+
{
54+
public static GradientComputator CreateGradientComputator(VolumeDataset dataset, GradientType gradientType)
55+
{
56+
switch (gradientType)
57+
{
58+
case GradientType.CentralDifference:
59+
return new CentralDifferenceGradientComputator(dataset, false);
60+
case GradientType.SmoothedCentralDifference:
61+
return new CentralDifferenceGradientComputator(dataset, true);
62+
case GradientType.Sobel:
63+
return new SobelGradientComputator(dataset, false);
64+
case GradientType.SmoothedSobel:
65+
return new SobelGradientComputator(dataset, true);
66+
default:
67+
throw new NotImplementedException();
68+
}
69+
}
70+
}
71+
}

‎Assets/Scripts/Gradients/GradientType.cs

+23
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
using System;
2+
using UnityEngine;
3+
4+
namespace UnityVolumeRendering
5+
{
6+
public enum GradientType
7+
{
8+
CentralDifference,
9+
SmoothedCentralDifference,
10+
Sobel,
11+
SmoothedSobel
12+
}
13+
14+
public class GradientTypeUtils
15+
{
16+
public static GradientType GetDefaultGradientType()
17+
{
18+
GradientType gradientType = GradientType.CentralDifference;
19+
Enum.TryParse(PlayerPrefs.GetString("DefaultGradientType"), out gradientType);
20+
return gradientType;
21+
}
22+
}
23+
}

‎Assets/Scripts/Gradients/SobelGradientComputator.cs

+73
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
using System;
2+
using System.Drawing;
3+
using UnityEngine;
4+
5+
namespace UnityVolumeRendering
6+
{
7+
public class SobelGradientComputator : GradientComputator
8+
{
9+
public SobelGradientComputator(VolumeDataset dataset, bool smootheDataValues) : base(dataset, smootheDataValues)
10+
{
11+
}
12+
13+
private static readonly float[,,] kernelx = new float[,,]
14+
{
15+
{ { -1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } },
16+
{ { -2, 0, 2 }, { -4, 0, 4 }, { -2, 0, 2 } },
17+
{ { -1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } }
18+
};
19+
20+
private static readonly float[,,] kernely = new float[,,]
21+
{
22+
{ { -1, -2, -1 }, { 0, 0, 0 }, { 1, 2, 1 } },
23+
{ { -2, -4, -2 }, { 0, 0, 0 }, { 2, 4, 2 } },
24+
{ { -1, -2, -1 }, { 0, 0, 0 }, { 1, 2, 1 } }
25+
};
26+
27+
private static readonly float[,,] kernelz = new float[,,]
28+
{
29+
{ { -1, -2, -1 }, { -2, -4, -2 }, { -1, -2, -1 } },
30+
{ { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } },
31+
{ { 1, 2, 1 }, { 2, 4, 2 }, { 1, 2, 1 } }
32+
};
33+
34+
private float GetData(int x, int y, int z)
35+
{
36+
return data[x + y * dimX + z * (dimX * dimY)];
37+
}
38+
39+
private Vector3 ConvolveWithKernels(int x, int y, int z)
40+
{
41+
Vector3 result = Vector3.zero;
42+
for (int iz = 0; iz <= 2; iz++)
43+
{
44+
for (int iy = 0; iy <= 2; iy++)
45+
{
46+
for (int ix = 0; ix <= 2; ix++)
47+
{
48+
float dataValue = GetData(x + ix - 1, y + iy - 1, z + iz - 1);
49+
result.x += kernelx[iz, iy, ix] * dataValue;
50+
result.y += kernely[iz, iy, ix] * dataValue;
51+
result.z += kernelz[iz, iy, ix] * dataValue;
52+
}
53+
}
54+
}
55+
return result;
56+
}
57+
58+
public override Vector3 ComputeGradient(int x, int y, int z, float minValue, float maxRange)
59+
{
60+
// TODO
61+
if (x < 2 || y < 2 || z < 2 || x > dimX - 3 || y > dimY - 3 || z > dimZ - 3)
62+
{
63+
return Vector3.zero;
64+
}
65+
66+
Vector3 gradient = ConvolveWithKernels(x, y, z);
67+
68+
float divident = maxRange * 12;
69+
70+
return new Vector3(gradient.x / divident, gradient.y / divident, gradient.z / divident);
71+
}
72+
}
73+
}

‎Assets/Scripts/VolumeData/VolumeDataset.cs

+32-26
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,7 @@ public Texture3D GetGradientTexture()
113113
{
114114
if (gradientTexture == null)
115115
{
116-
gradientTexture = AsyncHelper.RunSync<Texture3D>(() => CreateGradientTextureInternalAsync(NullProgressHandler.instance));
116+
gradientTexture = AsyncHelper.RunSync<Texture3D>(() => CreateGradientTextureInternalAsync(GradientTypeUtils.GetDefaultGradientType(), NullProgressHandler.instance));
117117
return gradientTexture;
118118
}
119119
else
@@ -122,6 +122,29 @@ public Texture3D GetGradientTexture()
122122
}
123123
}
124124

125+
public async Task<Texture3D> RegenerateGradientTextureAsync(GradientType gradientType, IProgressHandler progressHandler = null)
126+
{
127+
await createGradientTextureLock.WaitAsync();
128+
try
129+
{
130+
if (progressHandler == null)
131+
progressHandler = new NullProgressHandler();
132+
try
133+
{
134+
gradientTexture = await CreateGradientTextureInternalAsync(gradientType, progressHandler != null ? progressHandler : NullProgressHandler.instance);
135+
}
136+
catch (System.Exception exception)
137+
{
138+
Debug.LogException(exception);
139+
}
140+
}
141+
finally
142+
{
143+
createGradientTextureLock.Release();
144+
}
145+
return gradientTexture;
146+
}
147+
125148
/// <summary>
126149
/// Gets the gradient texture, containing the gradient values (direction of change) of the dataset.
127150
/// 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
132155
{
133156
if (gradientTexture == null)
134157
{
135-
await createGradientTextureLock.WaitAsync();
136-
try
137-
{
138-
if (progressHandler == null)
139-
progressHandler = new NullProgressHandler();
140-
gradientTexture = await CreateGradientTextureInternalAsync(progressHandler != null ? progressHandler : NullProgressHandler.instance);
141-
}
142-
finally
143-
{
144-
createGradientTextureLock.Release();
145-
}
158+
gradientTexture = await RegenerateGradientTextureAsync(GradientTypeUtils.GetDefaultGradientType(), progressHandler);
146159
}
147160
return gradientTexture;
148161
}
@@ -332,7 +345,7 @@ await Task.Run(() => {
332345
return texture;
333346
}
334347

335-
private async Task<Texture3D> CreateGradientTextureInternalAsync(IProgressHandler progressHandler)
348+
private async Task<Texture3D> CreateGradientTextureInternalAsync(GradientType gradientType, IProgressHandler progressHandler)
336349
{
337350
Debug.Log("Async gradient generation. Hold on.");
338351

@@ -364,6 +377,8 @@ await Task.Run(() => {
364377
Texture3D textureTmp = new Texture3D(dimX, dimY, dimZ, texformat, false);
365378
textureTmp.wrapMode = TextureWrapMode.Clamp;
366379

380+
GradientComputator gradientComputator = GradientComputatorFactory.CreateGradientComputator(this, gradientType);
381+
367382
for (int x = 0; x < dimX; x++)
368383
{
369384
progressHandler.ReportProgress(x, dimX, "Calculating gradients for slice");
@@ -372,7 +387,7 @@ await Task.Run(() => {
372387
for (int z = 0; z < dimZ; z++)
373388
{
374389
int iData = x + y * dimX + z * (dimX * dimY);
375-
Vector3 grad = GetGrad(x, y, z, minValue, maxRange);
390+
Vector3 grad = gradientComputator.ComputeGradient(x, y, z, minValue, maxRange);
376391

377392
textureTmp.SetPixel(x, y, z, new Color(grad.x, grad.y, grad.z, (float)(data[iData] - minValue) / maxRange));
378393
}
@@ -390,6 +405,8 @@ await Task.Run(() => {
390405

391406
progressHandler.StartStage(0.6f, "Creating gradient texture");
392407
await Task.Run(() => {
408+
GradientComputator gradientComputator = GradientComputatorFactory.CreateGradientComputator(this, gradientType);
409+
393410
for (int z = 0; z < dimZ; z++)
394411
{
395412
progressHandler.ReportProgress(z, dimZ, "Calculating gradients for slice");
@@ -398,7 +415,7 @@ await Task.Run(() => {
398415
for (int x = 0; x < dimX; x++)
399416
{
400417
int iData = x + y * dimX + z * (dimX * dimY);
401-
Vector3 grad = GetGrad(x, y, z, minValue, maxRange);
418+
Vector3 grad = gradientComputator.ComputeGradient(x, y, z, minValue, maxRange);
402419

403420
cols[iData] = new Color(grad.x, grad.y, grad.z, (float)(data[iData] - minValue) / maxRange);
404421
}
@@ -418,17 +435,6 @@ await Task.Run(() => {
418435
return texture;
419436

420437
}
421-
public Vector3 GetGrad(int x, int y, int z, float minValue, float maxRange)
422-
{
423-
float x1 = data[Math.Min(x + 1, dimX - 1) + y * dimX + z * (dimX * dimY)] - minValue;
424-
float x2 = data[Math.Max(x - 1, 0) + y * dimX + z * (dimX * dimY)] - minValue;
425-
float y1 = data[x + Math.Min(y + 1, dimY - 1) * dimX + z * (dimX * dimY)] - minValue;
426-
float y2 = data[x + Math.Max(y - 1, 0) * dimX + z * (dimX * dimY)] - minValue;
427-
float z1 = data[x + y * dimX + Math.Min(z + 1, dimZ - 1) * (dimX * dimY)] - minValue;
428-
float z2 = data[x + y * dimX + Math.Max(z - 1, 0) * (dimX * dimY)] - minValue;
429-
430-
return new Vector3((x2 - x1) / maxRange, (y2 - y1) / maxRange, (z2 - z1) / maxRange);
431-
}
432438

433439
public float GetAvgerageVoxelValues(int x, int y, int z)
434440
{

‎Assets/Scripts/VolumeObject/VolumeObjectFactory.cs

+2
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ public static VolumeRenderedObject CreateObject(VolumeDataset dataset)
1010
{
1111
GameObject outerObject = new GameObject("VolumeRenderedObject_" + dataset.datasetName);
1212
VolumeRenderedObject volObj = outerObject.AddComponent<VolumeRenderedObject>();
13+
volObj.SetGradientType(GradientTypeUtils.GetDefaultGradientType());
1314

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

2931
GameObject meshContainer = GameObject.Instantiate((GameObject)Resources.Load("VolumeContainer"));
3032
volObj.volumeContainerObject = meshContainer;

0 commit comments

Comments
 (0)