libift.c
2.28 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
//extern "C" {
#include "ift.h"
//}
// Prototypes
void NewDestroyScene(Scene *scn);
void NewDestroyScene(Scene *scn)
{
if(scn != NULL){
if (scn->data != NULL) free(scn->data);
if (scn->tby != NULL) free(scn->tby);
if (scn->tbz != NULL) free(scn->tbz);
free(scn);
}
}
Scene* EraseBackground(Scene *scn)
{
int otsu = (int)(Otsu3(scn)*0.8); // Get 80% of Otsu's threshold
int i,n = scn->xsize*scn->ysize*scn->zsize;
Scene *out = CopyScene(scn);
for (i=0; i < n; i++) {
if (out->data[i]<otsu) out->data[i]=0;
}
return out;
}
Scene* EraseSupport(Scene *scn)
{
// Check if RemoveBackground was already ran (workaround)
int c=0,i,n = scn->xsize*scn->ysize*scn->zsize;
Scene *scn2;
Scene *labels,*bin;
AdjRel3 *A;
int total[2000],p;
Voxel v;
int max;
Scene *out;
for (i=0; i < n; i++) {
if (scn->data[i]==0) c++; // count zero voxels
}
if (c<(n*0.15))
scn2 = EraseBackground(scn);
else
scn2 = CopyScene(scn);
// Get max connected component
A=Spheric(1.5);
bin = Threshold3(scn2,1,50000);
labels = LabelBinComp3(bin,A);
DestroyAdjRel3(&A);
// Count labels
for (i=0;i<2000;i++) total[i]=0;
for (v.z=0; v.z < scn->zsize; v.z++)
for (v.y=0; v.y < scn->ysize; v.y++)
for (v.x=0; v.x < scn->xsize; v.x++) {
p = labels->tbz[v.z] + labels->tby[v.y] + v.x;
if (labels->data[p]<2000) total[labels->data[p]]++;
if (bin->data[p]==0) total[labels->data[p]]=0; // exclude background
}
DestroyScene(&bin);
max=0;
for (i=0;i<2000;i++)
if (total[i]>total[max]) max=i;
// copy the maxcc to out
n = scn->xsize*scn->ysize*scn->zsize;
out = CreateScene(scn->xsize,scn->ysize,scn->zsize);
out->dx=scn->dx;
out->dy=scn->dy;
out->dz=scn->dz;
for (i=0; i < n; i++) {
if (labels->data[i]==max) out->data[i]=scn2->data[i];
}
DestroyScene(&labels);
DestroyScene(&scn2);
return out;
}
int ShiftScene(Scene *scn)
{
int n,i, min = MinimumValue3(scn);
if (min<0) {
n = scn->xsize*scn->ysize*scn->zsize;
for (i=0; i < n; i++)
scn->data[i] += 1024;
return 1;
}
return 0;
}
void UnShiftScene(Scene *scn, int flag)
{
int i,n;
if (flag!=0) {
n = scn->xsize*scn->ysize*scn->zsize;
for (i=0; i < n; i++)
scn->data[i] -= 1024;
}
}