Roussos-Maragos Tensor-Driven Diffusion Interpolation
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
conv.c
Go to the documentation of this file.
1 
16 #include <string.h>
17 #include "conv.h"
18 
19 
21 #define CLAMP(X,A,B) (((X) < (A)) ? (A) : (((X) > (B)) ? (B) : (X)))
22 
23 
25 const filter NullFilter = {NULL, 0, 0};
26 
27 
40 void SampledConv1D(float *Dest, int DestStride, const float *Src,
41  int SrcStride, filter Filter, boundaryext Boundary, int N,
42  int nStart, int nStep, int nEnd)
43 {
44  const int SrcStrideStep = SrcStride*nStep;
45  const int LeftmostTap = 1 - Filter.Delay - Filter.Length;
46  const int StartInterior = CLAMP(-LeftmostTap, 0, N - 1);
47  const int EndInterior = (Filter.Delay < 0) ?
48  (N + Filter.Delay - 1) : (N - 1);
49  const float *SrcN, *SrcK;
50  float Accum;
51  int n, k;
52 
53 
54  if(nEnd < nStart || nStep <= 0 || N <= 0)
55  return;
56 
57  /* Handle the left boundary */
58  for(n = nStart; n < StartInterior; n += nStep, Dest += DestStride)
59  {
60  for(k = 0, Accum = 0; k < Filter.Length; k++)
61  Accum += Filter.Coeff[k]
62  * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
63 
64  *Dest = Accum;
65  }
66 
67  /* Compute the convolution on the interior of the signal:
68 
69  In the inner accumulation loop
70  SrcK = &inputdata[n - FilterDelay - k], k = FilterLength-1, ..., 0.
71 
72  The SrcN pointer is adjusted such that
73  SrcN = &inputdata[n + LeftmostTap].
74 
75  If n == StartInterior, then the loop starts with
76  n = -LeftmostTap, SrcN = &inputdata[0] if LeftmostTap <= 0
77  n = 0, SrcN = &inputdata[LeftmostTap] if LeftmostTap >= 0. */
78  SrcN = (LeftmostTap <= 0) ? Src : (Src + SrcStride*LeftmostTap);
79 
80  /* Adjust if n > StartInterior */
81  SrcN += SrcStride*(n - StartInterior);
82 
83  for(; n <= EndInterior; n += nStep, SrcN += SrcStrideStep, Dest += DestStride)
84  {
85  Accum = 0;
86  SrcK = SrcN;
87  k = Filter.Length;
88 
89  while(k)
90  {
91  Accum += Filter.Coeff[--k] * (*SrcK);
92  SrcK += SrcStride;
93  }
94 
95  *Dest = Accum;
96  }
97 
98  /* Handle the right boundary */
99  for(; n <= nEnd; n += nStep, Dest += DestStride)
100  {
101  for(k = 0, Accum = 0; k < Filter.Length; k++)
102  Accum += Filter.Coeff[k]
103  * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
104 
105  *Dest = Accum;
106  }
107 }
108 
109 
123 void SeparableConv2D(float *Dest, float *Buffer, const float *Src,
124  filter FilterX, filter FilterY, boundaryext Boundary,
125  int Width, int Height, int NumChannels)
126 {
127  const int NumPixels = Width*Height;
128  int i, Channel;
129 
130  for(Channel = 0; Channel < NumChannels; Channel++)
131  {
132  /* Filter Src horizontally and store the result in Buffer */
133  for(i = 0; i < Height; i++)
134  Conv1D(Buffer + Width*i, 1, Src + Width*i, 1,
135  FilterX, Boundary, Width);
136 
137  /* Filter Buffer vertically and store the result in Dest */
138  for(i = 0; i < Width; i++)
139  Conv1D(Dest + i, Width, Buffer + i, Width, FilterY,
140  Boundary, Height);
141 
142  Src += NumPixels;
143  Dest += NumPixels;
144  }
145 }
146 
147 
149 filter MakeFilter(float *Coeff, int Delay, int Length)
150 {
151  filter Filter;
152 
153  Filter.Coeff = Coeff;
154  Filter.Delay = Delay;
155  Filter.Length = Length;
156  return Filter;
157 }
158 
159 
161 filter AllocFilter(int Delay, int Length)
162 {
163  float *Coeff;
164 
165  if(Length > 0 && (Coeff = (float *)Malloc(sizeof(float)*Length)))
166  return MakeFilter(Coeff, Delay, Length);
167  else
168  return NullFilter;
169 }
170 
171 
173 int IsNullFilter(filter Filter)
174 {
175  return (Filter.Coeff == NULL) ? 1:0;
176 }
177 
178 
193 filter GaussianFilter(double Sigma, int R)
194 {
195  filter Filter = AllocFilter(-R, 2*R+1);
196 
197 
198  if(!IsNullFilter(Filter))
199  {
200  if(Sigma == 0)
201  Filter.Coeff[0] = 1;
202  else
203  {
204  float Sum;
205  int r;
206 
207  for(r = -R, Sum = 0; r <= R; r++)
208  {
209  Filter.Coeff[R + r] = (float)exp(-r*r/(2*Sigma*Sigma));
210  Sum += Filter.Coeff[R + r];
211  }
212 
213  for(r = -R; r <= R; r++)
214  Filter.Coeff[R + r] /= Sum;
215  }
216  }
217 
218  return Filter;
219 }
220 
221 
222 static float ZeroPaddedExtension(const float *Src, int Stride, int N, int n)
223 {
224  return (0 <= n && n < N) ? Src[Stride*n] : 0;
225 }
226 
227 
228 static float ConstantExtension(const float *Src, int Stride, int N, int n)
229 {
230  return Src[(n < 0) ? 0 : ((n >= N) ? Stride*(N - 1) : n)];
231 }
232 
233 
234 static float LinearExtension(const float *Src, int Stride, int N, int n)
235 {
236  if(0 <= n)
237  {
238  if(n < N)
239  return Src[Stride*n];
240  else if(N == 1)
241  return Src[0];
242  else
243  {
244  Src += Stride*(N - 1);
245  return Src[0] + (N - 1 - n)*(Src[-Stride] - Src[0]);
246  }
247  }
248  else if(N == 1)
249  return Src[0];
250  else
251  return Src[0] + n*(Src[Stride] - Src[0]);
252 }
253 
254 
255 static float PeriodicExtension(const float *Src, int Stride, int N, int n)
256 {
257  if(n < 0)
258  {
259  do
260  {
261  n += N;
262  }while(n < 0);
263  }
264  else if(n >= N)
265  {
266  do
267  {
268  n -= N;
269  }while(n >= N);
270  }
271 
272  return Src[Stride*n];
273 }
274 
275 
276 static float SymhExtension(const float *Src, int Stride, int N, int n)
277 {
278  while(1)
279  {
280  if(n < 0)
281  n = -1 - n;
282  else if(n >= N)
283  n = 2*N - 1 - n;
284  else
285  break;
286  }
287 
288  return Src[Stride*n];
289 }
290 
291 
292 static float SymwExtension(const float *Src, int Stride, int N, int n)
293 {
294  while(1)
295  {
296  if(n < 0)
297  n = -n;
298  else if(n >= N)
299  n = 2*(N - 1) - n;
300  else
301  break;
302  }
303 
304  return Src[Stride*n];
305 }
306 
307 
308 static float AsymhExtension(const float *Src, int Stride, int N, int n)
309 {
310  float Jump, Offset;
311 
312 
313  /* Use simple formulas for -N <= n <= 2*N - 1 */
314  if(0 <= n)
315  {
316  if(n < N)
317  return Src[Stride*n];
318  else if(n <= 2*N - 1)
319  return 3*Src[Stride*(N - 1)] - Src[Stride*(N - 2)]
320  - Src[Stride*(2*N - 1 - n)];
321  }
322  else if(-N <= n)
323  return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)];
324 
325  /* N == 1 is a special case */
326  if(N == 1)
327  return Src[0];
328 
329  /* General formula for extension at an arbitrary n */
330  Jump = 3*(Src[Stride*(N - 1)] - Src[0])
331  - (Src[Stride*(N - 2)] - Src[Stride]);
332  Offset = 0;
333 
334  if(n >= N)
335  {
336  do
337  {
338  Offset += Jump;
339  n -= 2*N;
340  }while(n >= N);
341  }
342  else
343  {
344  while(n < -N)
345  {
346  Offset -= Jump;
347  n += 2*N;
348  }
349  }
350 
351  if(n >= 0)
352  return Src[Stride*n] + Offset;
353  else
354  return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)] + Offset;
355 }
356 
357 
358 static float AsymwExtension(const float *Src, int Stride, int N, int n)
359 {
360  float Jump, Offset;
361 
362 
363  /* Use simple formulas for -N < n < 2*N - 1 */
364  if(0 <= n)
365  {
366  if(n < N)
367  return Src[Stride*n];
368  else if(n < 2*N - 1)
369  return 2*Src[Stride*(N - 1)] - Src[Stride*(2*(N - 1) - n)];
370  }
371  else if(-N < n)
372  return 2*Src[0] - Src[Stride*(-n)];
373 
374  /* N == 1 is a special case */
375  if(N == 1)
376  return Src[0];
377 
378  /* General formula for extension at an arbitrary n */
379  Jump = 2*(Src[Stride*(N - 1)] - Src[0]);
380  Offset = 0;
381 
382  if(n >= N)
383  {
384  do
385  {
386  Offset += Jump;
387  n -= 2*(N - 1);
388  }while(n >= N);
389  }
390  else
391  {
392  while(n <= -N)
393  {
394  Offset -= Jump;
395  n += 2*(N - 1);
396  }
397  }
398 
399  if(n >= 0)
400  return Src[Stride*n] + Offset;
401  else
402  return 2*Src[0] - Src[Stride*(-n)] + Offset;
403 }
404 
405 
421 boundaryext GetBoundaryExt(const char *Boundary)
422 {
423  if(!strcmp(Boundary, "zpd") || !strcmp(Boundary, "zero"))
424  return ZeroPaddedExtension;
425  else if(!strcmp(Boundary, "sp0") || !strcmp(Boundary, "const"))
426  return ConstantExtension;
427  else if(!strcmp(Boundary, "sp1") || !strcmp(Boundary, "linear"))
428  return LinearExtension;
429  else if(!strcmp(Boundary, "per") || !strcmp(Boundary, "periodic"))
430  return PeriodicExtension;
431  else if(!strcmp(Boundary, "sym")
432  || !strcmp(Boundary, "symh") || !strcmp(Boundary, "hsym"))
433  return SymhExtension;
434  else if(!strcmp(Boundary, "symw") || !strcmp(Boundary, "wsym"))
435  return SymwExtension;
436  else if(!strcmp(Boundary, "asym")
437  || !strcmp(Boundary, "asymh") || !strcmp(Boundary, "hasym"))
438  return AsymhExtension;
439  else if(!strcmp(Boundary, "asymw") || !strcmp(Boundary, "wasym"))
440  return AsymwExtension;
441  else
442  {
443  ErrorMessage("Unknown boundary extension \"%s\".\n", Boundary);
444  return NULL;
445  }
446 }