Mercurial > forge

comparison main/image/cordflt2.cc @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 82264140e58b
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 // Copyright (C) 2000 Teemu Ikonen
2 //
3 // This program is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU General Public License
5 // as published by the Free Software Foundation; either version 2
6 // of the License, or (at your option) any later version.
7 //
8 // This program is distributed in the hope that it will be useful, but
9 // WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with this program; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
16
17 #include <octave/oct.h>
18
19 #ifdef HAVE_OCTAVE_20
20 typedef Matrix boolMatrix;
21 #define bool_matrix_value matrix_value
22 #endif
23
24 #define SWAP(a, b) { SWAP_temp = (a); (a)=(b); (b) = SWAP_temp; }
25
26 // Template function for comparison
27 // ET is the type of the matrix element
28 template <class ET>
29 inline bool compare(const ET a, const ET b)
30 {
31 if(a > b)
32 return 1;
33 else
34 return 0;
35 }
36
37 // Explicit template function for complex compare
38 template <> inline bool compare<Complex>(const Complex a, const Complex b)
39 {
40 double anorm2 = a.real() * a.real() + a.imag() * a.imag();
41 double bnorm2 = b.real() * b.real() + b.imag() * b.imag();
42
43 if( anorm2 > bnorm2 ) {
44 return 1;
45 } else {
46 return 0;
47 }
48 }
49
50 // select nth largest member from the array values
51 // Partitioning algorithm, see Numerical recipes chap. 8.5
52 template <class ET>
53 ET selnth(ET *vals, int len, int nth)
54 {
55 ET SWAP_temp;
56 ET hinge;
57 int l, r, mid, i, j;
58
59 l = 0;
60 r = len - 1;
61 for(;;) {
62 // if partition size is 1 or two, then sort and return
63 if(r <= l+1) {
64 if(r == l+1 && compare<ET>(vals[l], vals[r])) {
65 SWAP(vals[l], vals[r]);
66 }
67 return vals[nth];
68 } else {
69 mid = (l+r) >> 1;
70 SWAP(vals[mid], vals[l+1]);
71 // choose median of l, mid, r to be the hinge element
72 // and set up sentinels in the borders (order l, l+1 and r)
73 if(compare<ET>(vals[l], vals[r])) {
74 SWAP(vals[l], vals[r]);
75 }
76 if(compare<ET>(vals[l+1], vals[r])) {
77 SWAP(vals[l+1], vals[r]);
78 }
79 if(compare<ET>(vals[l], vals[l+1])) {
80 SWAP(vals[l], vals[l+1]);
81 }
82 i = l+1;
83 j = r;
84 hinge = vals[l+1];
85 for(;;) {
86 do i++; while(compare<ET>(hinge, vals[i]));
87 do j--; while(compare<ET>(vals[j], hinge));
88 if(i > j)
89 break;
90 SWAP(vals[i], vals[j]);
91 }
92 vals[l+1] = vals[j];
93 vals[j] = hinge;
94 if(j >= nth)
95 r = j - 1;
96 if(j <= nth)
97 l = i;
98 }
99 }
100 }
101
102 // Template function for doing the actual filtering
103 // MT is the type of the matrix to be filtered (Matrix or ComplexMatrix)
104 // ET is the type of the element of the matrix (double or Complex)
105 template <class MT, class ET>
106 octave_value_list do_filtering(MT A, int nth, boolMatrix dom, MT S)
107 {
108 int i, j, c, d;
109
110 int len = 0;
111 for(j = 0; j < dom.columns(); j++) {
112 for(i = 0; i < dom.rows(); i++) {
113 if(dom.elem(i,j))
114 len++;
115 }
116 }
117 if(nth > len - 1) {
118 warning("nth should be less than number of non-zero values in domain");
119 warning("setting nth to largest possible value\n");
120 nth = len - 1;
121 }
122 if(nth < 0) {
123 warning("nth should be non-negative, setting to 1\n");
124 nth = 0; // nth is a c-index
125 }
126
127 int rowoffset = (dom.columns() + 1)/2 - 1;
128 int coloffset = (dom.rows() + 1)/2 - 1;
129
130 //outputs
131 octave_value_list out;
132 const int origx = A.columns() - dom.columns()+1;
133 const int origy = A.rows() - dom.rows()+1;
134 MT retval = MT(origy, origx);
135
136 int *offsets = new int[len];
137 ET *values = new ET[len];
138 ET *adds = new ET[len];
139
140 c = 0;
141 d = A.rows();
142 for(j = 0; j < dom.columns(); j++) {
143 for(i = 0; i < dom.rows(); i++) {
144 if(dom.elem(i,j)) {
145 offsets[c] = (i - coloffset) + (j - rowoffset)*d;
146 adds[c] = S.elem(i,j);
147 c++;
148 }
149 }
150 }
151
152 ET *data = A.fortran_vec();
153 int base = coloffset + A.rows()*rowoffset;
154 for(j = 0; j < retval.columns(); j++) {
155 for(i = 0; i < retval.rows(); i++) {
156 for(c = 0; c < len; c++) {
157 values[c] = data[base + offsets[c]] + adds[c];
158 }
159 base++;
160 retval(i, j) = selnth(values, len, nth);
161 }
162 base += dom.rows() - 1;
163 }
164
165 out(0) = octave_value(retval);
166
167 return out;
168 }
169
170 // instantiate template functions
171 template inline bool compare<double>(const double, const double);
172 template double selnth(double *, int, int);
173 template Complex selnth(Complex *, int, int);
174 template octave_value_list do_filtering<Matrix, double>(Matrix, int, boolMatrix, Matrix);
175 // g++ is broken, explicit instantiation of specialized template function
176 // confuses the compiler.
177 //template int compare<Complex>(const Complex, const Complex);
178 template octave_value_list do_filtering<ComplexMatrix, Complex>(ComplexMatrix, int, boolMatrix, ComplexMatrix);
179
180 DEFUN_DLD(cordflt2, args, ,
181 "function retval = cordflt2(A, nth, domain, S)
182
183 Implementation of two-dimensional ordered filtering. User interface
184 in ordfilt2.m
185 ")
186 {
187 if(args.length() != 4) {
188 print_usage ("ordfilt2");
189 return octave_value_list();
190 }
191
192 // nth is an index to an array, thus - 1
193 int nth = (int) (args(1).vector_value())(0) - 1;
194 boolMatrix dom = args(2).bool_matrix_value();
195
196 octave_value_list retval;
197
198 if(args(0).is_real_matrix()) {
199 Matrix A = args(0).matrix_value();
200 Matrix S = args(3).matrix_value();
201 retval = do_filtering<Matrix, double>(A, nth, dom, S);
202 }
203 else if(args(0).is_complex_matrix()) {
204 ComplexMatrix A = args(0).complex_matrix_value();
205 ComplexMatrix S = args(3).complex_matrix_value();
206 retval = do_filtering<ComplexMatrix, Complex>(A, nth, dom, S);
207 }
208 else {
209 error("A should be real or complex matrix\n");
210 return octave_value_list();
211 }
212
213 return retval;
214
215 }