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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
/***************************************************************************
filter_morph.h - description
-------------------
begin : Fr Nov 29 2019
copyright : (C) 2019-2021 by Thies Jochimsen
email : thies@jochimsen.de
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#ifndef FILTER_MORPH_H
#define FILTER_MORPH_H
#include <odindata/filter_step.h>
#include <odindata/utils.h> // for on_grid(...)
enum morphOp { erode=0, dilate };
template <morphOp morphop>
class FilterMorph: public FilterStep {
LDRfloat radius;
STD_string label() const {if(morphop==erode) return "erode"; else return "dilate";}
STD_string description() const {return label()+" image using spherical kernel as structuring element";}
bool process(Data<float,4>& data, Protocol& prot) const {
Log<Filter> odinlog(c_label(),"process");
Range all=Range::all();
TinyVector<int,4> shape=data.shape();
TinyVector<int,3> spatshape(shape(sliceDim), shape(phaseDim), shape(readDim));
TinyVector<float,3> voxel_spacing;
voxel_spacing(0)=FileFormat::voxel_extent(prot.geometry, sliceDirection, data.extent(sliceDim));
voxel_spacing(1)=FileFormat::voxel_extent(prot.geometry, phaseDirection, data.extent(phaseDim));
voxel_spacing(2)=FileFormat::voxel_extent(prot.geometry, readDirection, data.extent(readDim));
ODINLOG(odinlog,normalDebug) << "voxel_spacing=" << voxel_spacing << STD_endl;
TinyVector<int,3> enclbox;
enclbox=radius/voxel_spacing+1;
ODINLOG(odinlog,normalDebug) << "enclbox=" << enclbox << STD_endl;
STD_vector<TinyVector<int,3> > neighboffset;
for(int k=-enclbox(0); k<=enclbox(0); k++) {
for(int j=-enclbox(1); j<=enclbox(1); j++) {
for(int i=-enclbox(2); i<=enclbox(2); i++) {
TinyVector<float,3> distvec(voxel_spacing(0)*k, voxel_spacing(1)*j, voxel_spacing(2)*i);
float dist=sqrt(sum(distvec*distvec));
if(dist<=radius) {
neighboffset.push_back(TinyVector<int,3>( k, j, i));
}
}
}
}
unsigned int numof_neigb=neighboffset.size();
ODINLOG(odinlog,normalDebug) << "numof_neigb=" << numof_neigb << STD_endl;
for(unsigned int i=0; i<numof_neigb; i++) {
ODINLOG(odinlog,normalDebug) << "neighboffset[" << i << "]=" << neighboffset[i] << STD_endl;
}
Data<float,3> involdata(spatshape);
Data<unsigned int,3> neighbcount(spatshape); neighbcount=0;
Data<float,3> outmask(spatshape); outmask=0.0;
TinyVector<int,3> index, neighbindex;
for(int irep=0; irep<shape(timeDim); irep++) {
involdata(all,all,all)=data(irep,all,all,all);
for(unsigned int i=0; i<involdata.size(); i++) {
index=involdata.create_index(i);
if(involdata(index)!=0.0) {
for(unsigned int j=0; j<numof_neigb; j++) {
neighbindex=index+neighboffset[j];
if(on_grid<3>(spatshape, neighbindex)) {
neighbcount(neighbindex)++;
}
}
}
}
for(unsigned int i=0; i<involdata.size(); i++) {
index=involdata.create_index(i);
unsigned int count=neighbcount(index);
if(morphop==erode && count==numof_neigb) outmask(index)=1.0;
if(morphop==dilate && count>0) outmask(index)=1.0;
}
data(irep,all,all,all)=outmask(all,all,all);
}
return true;
}
FilterStep* allocate() const {return new FilterMorph<morphop>();}
void init() {
radius.set_unit(ODIN_SPAT_UNIT).set_description("radius of kernel");
append_arg(radius,"radius");
}
};
///////////////////////////////////////////////////////////////////////////
#endif
|
