1 | /** file: imsmooth.c |
---|
2 | ** author: Andrea Vedaldi |
---|
3 | ** description: Smooth an image. |
---|
4 | **/ |
---|
5 | |
---|
6 | /* AUTORIGHTS |
---|
7 | Copyright (C) 2006 Andrea Vedaldi |
---|
8 | |
---|
9 | This file is part of VLUtil. |
---|
10 | |
---|
11 | VLUtil is free software; you can redistribute it and/or modify |
---|
12 | it under the terms of the GNU General Public License as published by |
---|
13 | the Free Software Foundation; either version 2, or (at your option) |
---|
14 | any later version. |
---|
15 | |
---|
16 | This program is distributed in the hope that it will be useful, |
---|
17 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
19 | GNU General Public License for more details. |
---|
20 | |
---|
21 | You should have received a copy of the GNU General Public License |
---|
22 | along with this program; if not, write to the Free Software Foundation, |
---|
23 | Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
---|
24 | */ |
---|
25 | |
---|
26 | #include"mex.h" |
---|
27 | #include<stdlib.h> |
---|
28 | #include<string.h> |
---|
29 | #include<math.h> |
---|
30 | #include<assert.h> |
---|
31 | |
---|
32 | #define greater(a,b) ((a) > (b)) |
---|
33 | #define min(a,b) (((a)<(b))?(a):(b)) |
---|
34 | #define max(a,b) (((a)>(b))?(a):(b)) |
---|
35 | |
---|
36 | #define PAD_BY_CONTINUITY |
---|
37 | |
---|
38 | const double win_factor = 1.5 ; |
---|
39 | const int nbins = 36 ; |
---|
40 | |
---|
41 | /* convolve and transopose */ |
---|
42 | void |
---|
43 | convolve(double* dst_pt, |
---|
44 | const double* src_pt, int M, int N, |
---|
45 | const double* filter_pt, int W) |
---|
46 | { |
---|
47 | int i,j ; |
---|
48 | for(j = 0 ; j < N ; ++j) { |
---|
49 | for(i = 0 ; i < M ; ++i) { |
---|
50 | double const* start = src_pt + max(i-W,0 ) ; |
---|
51 | double const* stop = src_pt + min(i+W,M-1) + 1 ; |
---|
52 | double const* g = filter_pt + max(0, W-i) ; |
---|
53 | double acc = 0.0 ; |
---|
54 | while(stop != start) acc += (*g++) * (*start++) ; |
---|
55 | *dst_pt = acc ; |
---|
56 | dst_pt += N ; |
---|
57 | } |
---|
58 | src_pt += M ; |
---|
59 | dst_pt -= M*N - 1 ; |
---|
60 | } |
---|
61 | } |
---|
62 | |
---|
63 | /* convolve and transopose and pad by continuity*/ |
---|
64 | void |
---|
65 | econvolve(double* dst_pt, |
---|
66 | const double* src_pt, int M, int N, |
---|
67 | const double* filter_pt, int W) |
---|
68 | { |
---|
69 | int i,j ; |
---|
70 | for(j = 0 ; j < N ; ++j) { |
---|
71 | for(i = 0 ; i < M ; ++i) { |
---|
72 | double acc = 0.0 ; |
---|
73 | double const* g = filter_pt ; |
---|
74 | double const* start = src_pt + (i-W) ; |
---|
75 | double const* stop ; |
---|
76 | double x ; |
---|
77 | |
---|
78 | /* beginning */ |
---|
79 | stop = src_pt + max(0, i-W) ; |
---|
80 | x = *stop ; |
---|
81 | while( start <= stop ) { acc += (*g++) * x ; start++ ; } |
---|
82 | |
---|
83 | /* middle */ |
---|
84 | stop = src_pt + min(M-1, i+W) ; |
---|
85 | while( start < stop ) acc += (*g++) * (*start++) ; |
---|
86 | |
---|
87 | /* end */ |
---|
88 | x = *start ; |
---|
89 | stop = src_pt + (i+W) ; |
---|
90 | while( start <= stop ) { acc += (*g++) * x ; start++ ; } |
---|
91 | |
---|
92 | /* save */ |
---|
93 | *dst_pt = acc ; |
---|
94 | dst_pt += N ; |
---|
95 | } |
---|
96 | /* next column */ |
---|
97 | src_pt += M ; |
---|
98 | dst_pt -= M*N - 1 ; |
---|
99 | } |
---|
100 | } |
---|
101 | |
---|
102 | void |
---|
103 | mexFunction(int nout, mxArray *out[], |
---|
104 | int nin, const mxArray *in[]) |
---|
105 | { |
---|
106 | int M,N,K,ndims ; |
---|
107 | int const *dims ; |
---|
108 | double* I_pt ; |
---|
109 | double* J_pt ; |
---|
110 | double s ; |
---|
111 | enum {I=0,S} ; |
---|
112 | enum {J=0} ; |
---|
113 | |
---|
114 | /* ------------------------------------------------------------------ |
---|
115 | ** Check the arguments |
---|
116 | ** --------------------------------------------------------------- */ |
---|
117 | if (nin != 2) { |
---|
118 | mexErrMsgTxt("Exactly two input arguments required."); |
---|
119 | } else if (nout > 1) { |
---|
120 | mexErrMsgTxt("Too many output arguments."); |
---|
121 | } |
---|
122 | |
---|
123 | if (!mxIsDouble(in[I]) || |
---|
124 | !mxIsDouble(in[S]) ) { |
---|
125 | mexErrMsgTxt("All arguments must be real.") ; |
---|
126 | } |
---|
127 | |
---|
128 | if(mxGetNumberOfDimensions(in[I]) > 3|| |
---|
129 | mxGetNumberOfDimensions(in[S]) > 2) { |
---|
130 | mexErrMsgTxt("I must be a two dimensional array and S a scalar.") ; |
---|
131 | } |
---|
132 | |
---|
133 | if(max(mxGetM(in[S]),mxGetN(in[S])) > 1) { |
---|
134 | mexErrMsgTxt("S must be a scalar.\n") ; |
---|
135 | } |
---|
136 | |
---|
137 | ndims = mxGetNumberOfDimensions(in[I]) ; |
---|
138 | dims = mxGetDimensions(in[I]) ; |
---|
139 | M = dims[0] ; |
---|
140 | N = dims[1] ; |
---|
141 | K = (ndims > 2) ? dims[2] : 1 ; |
---|
142 | |
---|
143 | out[J] = mxCreateNumericArray(ndims, dims, mxDOUBLE_CLASS, mxREAL) ; |
---|
144 | |
---|
145 | I_pt = mxGetPr(in[I]) ; |
---|
146 | J_pt = mxGetPr(out[J]) ; |
---|
147 | s = *mxGetPr(in[S]) ; |
---|
148 | |
---|
149 | /* ------------------------------------------------------------------ |
---|
150 | ** Do the job |
---|
151 | ** --------------------------------------------------------------- */ |
---|
152 | if(s > 0.01) { |
---|
153 | |
---|
154 | int W = (int) ceil(4*s) ; |
---|
155 | int j,k ; |
---|
156 | double* g0 = (double*) mxMalloc( (2*W+1)*sizeof(double) ) ; |
---|
157 | double* buffer = (double*) mxMalloc( M*N*sizeof(double) ) ; |
---|
158 | double acc=0.0 ; |
---|
159 | |
---|
160 | for(j = 0 ; j < 2*W+1 ; ++j) { |
---|
161 | g0[j] = exp(-0.5 * (j - W)*(j - W)/(s*s)) ; |
---|
162 | acc += g0[j] ; |
---|
163 | } |
---|
164 | for(j = 0 ; j < 2*W+1 ; ++j) { |
---|
165 | g0[j] /= acc ; |
---|
166 | } |
---|
167 | |
---|
168 | for(k = 0 ; k < K ; ++k) { |
---|
169 | #if defined PAD_BY_CONTINUITY |
---|
170 | |
---|
171 | econvolve(buffer, I_pt, M, N, g0, W) ; |
---|
172 | econvolve(J_pt, buffer, N, M, g0, W) ; |
---|
173 | |
---|
174 | #else |
---|
175 | |
---|
176 | convolve(buffer, I_pt, M, N, g0, W) ; |
---|
177 | convolve(J_pt, buffer, N, M, g0, W) ; |
---|
178 | |
---|
179 | #endif |
---|
180 | I_pt += M*N ; |
---|
181 | J_pt += M*N ; |
---|
182 | } |
---|
183 | |
---|
184 | mxFree(buffer) ; |
---|
185 | mxFree(g0) ; |
---|
186 | } else { |
---|
187 | memcpy(J_pt, I_pt, sizeof(double)*M*N) ; |
---|
188 | } |
---|
189 | } |
---|