Image Interpolation with Contour Stencils

conv.c

Go to the documentation of this file.
00001 
00016 #include <string.h>
00017 #include "conv.h"
00018 
00019 
00021 #define CLAMP(X,A,B)    (((X) < (A)) ? (A) : (((X) > (B)) ? (B) : (X)))
00022 
00023 
00025 const filter NullFilter = {NULL, 0, 0};
00026 
00027 
00040 void SampledConv1D(float *Dest, int DestStride, const float *Src,
00041     int SrcStride, filter Filter, boundaryext Boundary, int N, 
00042     int nStart, int nStep, int nEnd)
00043 {    
00044     const int SrcStrideStep = SrcStride*nStep;
00045     const int LeftmostTap = 1 - Filter.Delay - Filter.Length;    
00046     const int StartInterior = CLAMP(-LeftmostTap, 0, N - 1);
00047     const int EndInterior = (Filter.Delay <  0) ? 
00048                             (N + Filter.Delay - 1) : (N - 1);    
00049     const float *SrcN, *SrcK;
00050     float Accum;
00051     int n, k;
00052 
00053     
00054     if(nEnd < nStart || nStep <= 0 || N <= 0)
00055         return;
00056     
00057     /* Handle the left boundary */
00058     for(n = nStart; n < StartInterior; n += nStep, Dest += DestStride)
00059     {
00060         for(k = 0, Accum = 0; k < Filter.Length; k++)
00061             Accum += Filter.Coeff[k] 
00062                 * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
00063         
00064         *Dest = Accum;
00065     }
00066     
00067     /* Compute the convolution on the interior of the signal:
00068     
00069        In the inner accumulation loop
00070        SrcK = &inputdata[n - FilterDelay - k],  k = FilterLength-1, ..., 0.
00071        
00072        The SrcN pointer is adjusted such that
00073        SrcN = &inputdata[n + LeftmostTap]. 
00074        
00075        If n == StartInterior, then the loop starts with
00076           n = -LeftmostTap, SrcN = &inputdata[0]  if LeftmostTap <= 0
00077           n = 0, SrcN = &inputdata[LeftmostTap]   if LeftmostTap >= 0. */
00078     SrcN = (LeftmostTap <= 0) ? Src : (Src + SrcStride*LeftmostTap);
00079     
00080     /* Adjust if n > StartInterior */
00081     SrcN += SrcStride*(n - StartInterior);
00082     
00083     for(; n <= EndInterior; n += nStep, SrcN += SrcStrideStep, Dest += DestStride)
00084     {
00085         Accum = 0;
00086         SrcK = SrcN;
00087         k = Filter.Length;
00088         
00089         while(k)
00090         {
00091             Accum += Filter.Coeff[--k] * (*SrcK);
00092             SrcK += SrcStride;
00093         }
00094         
00095         *Dest = Accum;
00096     }
00097     
00098     /* Handle the right boundary */
00099     for(; n <= nEnd; n += nStep, Dest += DestStride)
00100     {
00101         for(k = 0, Accum = 0; k < Filter.Length; k++)
00102             Accum += Filter.Coeff[k]
00103                 * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
00104         
00105         *Dest = Accum;
00106     }
00107 }
00108 
00109 
00123 void SeparableConv2D(float *Dest, float *Buffer, const float *Src,
00124     filter FilterX, filter FilterY, boundaryext Boundary, 
00125     int Width, int Height, int NumChannels)
00126 {
00127     const int NumPixels = Width*Height;
00128     int i, Channel;
00129     
00130     for(Channel = 0; Channel < NumChannels; Channel++)
00131     {
00132         /* Filter Src horizontally and store the result in Buffer */
00133         for(i = 0; i < Height; i++)
00134             Conv1D(Buffer + Width*i, 1, Src + Width*i, 1, 
00135                 FilterX, Boundary, Width);
00136         
00137         /* Filter Buffer vertically and store the result in Dest */
00138         for(i = 0; i < Width; i++)
00139             Conv1D(Dest + i, Width, Buffer + i, Width, FilterY, 
00140                 Boundary, Height);
00141             
00142         Src += NumPixels;
00143         Dest += NumPixels;
00144     }
00145 }
00146 
00147 
00149 filter MakeFilter(float *Coeff, int Delay, int Length)
00150 {
00151     filter Filter;
00152     
00153     Filter.Coeff = Coeff;
00154     Filter.Delay = Delay;
00155     Filter.Length = Length;
00156     return Filter;
00157 }
00158 
00159 
00161 filter AllocFilter(int Delay, int Length)
00162 {
00163     float *Coeff;
00164     
00165     if(Length > 0 && (Coeff = (float *)Malloc(sizeof(float)*Length)))
00166         return MakeFilter(Coeff, Delay, Length);
00167     else
00168         return NullFilter;
00169 }
00170 
00171 
00173 int IsNullFilter(filter Filter)
00174 {
00175     return (Filter.Coeff == NULL) ? 1:0;
00176 }
00177 
00178 
00193 filter GaussianFilter(double Sigma, int R)
00194 {
00195     filter Filter = AllocFilter(-R, 2*R+1);
00196     
00197 
00198     if(!IsNullFilter(Filter))
00199     {
00200         if(Sigma == 0)
00201             Filter.Coeff[0] = 1;
00202         else
00203         {
00204             float Sum;
00205             int r;
00206             
00207             for(r = -R, Sum = 0; r <= R; r++)
00208             {
00209                 Filter.Coeff[R + r] = (float)exp(-r*r/(2*Sigma*Sigma));
00210                 Sum += Filter.Coeff[R + r];
00211             }
00212             
00213             for(r = -R; r <= R; r++)
00214                 Filter.Coeff[R + r] /= Sum;
00215         }
00216     }
00217     
00218     return Filter;
00219 }
00220 
00221 
00222 static float ZeroPaddedExtension(const float *Src, int Stride, int N, int n)
00223 {
00224     return (0 <= n && n < N) ? Src[Stride*n] : 0;
00225 }
00226 
00227 
00228 static float ConstantExtension(const float *Src, int Stride, int N, int n)
00229 {
00230     return Src[(n < 0) ? 0 : ((n >= N) ? Stride*(N - 1) : n)];
00231 }
00232 
00233 
00234 static float LinearExtension(const float *Src, int Stride, int N, int n)
00235 {
00236     if(0 <= n)
00237     {
00238         if(n < N)
00239             return Src[Stride*n];
00240         else if(N == 1)
00241             return Src[0];
00242         else
00243         {
00244             Src += Stride*(N - 1);
00245             return Src[0] + (N - 1 - n)*(Src[-Stride] - Src[0]);
00246         }
00247     }
00248     else if(N == 1)
00249         return Src[0];
00250     else
00251         return Src[0] + n*(Src[Stride] - Src[0]);
00252 }
00253 
00254 
00255 static float PeriodicExtension(const float *Src, int Stride, int N, int n)
00256 {
00257     if(n < 0)
00258     {
00259         do
00260         {
00261             n += N;
00262         }while(n < 0);
00263     }
00264     else if(n >= N)
00265     {
00266         do
00267         {
00268             n -= N;
00269         }while(n >= N);
00270     }
00271     
00272     return Src[Stride*n];
00273 }
00274 
00275 
00276 static float SymhExtension(const float *Src, int Stride, int N, int n)
00277 {
00278     while(1)
00279     {
00280         if(n < 0)
00281             n = -1 - n;
00282         else if(n >= N)
00283             n = 2*N - 1 - n;
00284         else
00285             break;
00286     }
00287     
00288     return Src[Stride*n];
00289 }
00290 
00291 
00292 static float SymwExtension(const float *Src, int Stride, int N, int n)
00293 {
00294     while(1)
00295     {
00296         if(n < 0)
00297             n = -n;
00298         else if(n >= N)
00299             n = 2*(N - 1) - n;
00300         else
00301             break;
00302     }
00303     
00304     return Src[Stride*n];
00305 }
00306 
00307 
00308 static float AsymhExtension(const float *Src, int Stride, int N, int n)
00309 {
00310     float Jump, Offset;
00311     
00312     
00313     /* Use simple formulas for -N <= n <= 2*N - 1 */
00314     if(0 <= n)
00315     {
00316         if(n < N)
00317             return Src[Stride*n];        
00318         else if(n <= 2*N - 1)
00319             return 3*Src[Stride*(N - 1)] - Src[Stride*(N - 2)]
00320                 - Src[Stride*(2*N - 1 - n)];
00321     }
00322     else if(-N <= n)
00323         return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)];
00324     
00325     /* N == 1 is a special case */
00326     if(N == 1)
00327         return Src[0];
00328     
00329     /* General formula for extension at an arbitrary n */
00330     Jump = 3*(Src[Stride*(N - 1)] - Src[0]) 
00331         - (Src[Stride*(N - 2)] - Src[Stride]);
00332     Offset = 0;        
00333     
00334     if(n >= N)
00335     {
00336         do
00337         {
00338             Offset += Jump;
00339             n -= 2*N;
00340         }while(n >= N);
00341     }
00342     else
00343     {
00344         while(n < -N)
00345         {
00346             Offset -= Jump;
00347             n += 2*N;
00348         }
00349     }
00350     
00351     if(n >= 0)
00352         return Src[Stride*n] + Offset;
00353     else
00354         return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)] + Offset;
00355 }
00356 
00357 
00358 static float AsymwExtension(const float *Src, int Stride, int N, int n)
00359 {
00360     float Jump, Offset;
00361     
00362     
00363     /* Use simple formulas for -N < n < 2*N - 1 */
00364     if(0 <= n)
00365     {
00366         if(n < N)
00367             return Src[Stride*n];        
00368         else if(n < 2*N - 1)
00369             return 2*Src[Stride*(N - 1)] - Src[Stride*(2*(N - 1) - n)];
00370     }
00371     else if(-N < n)
00372         return 2*Src[0] - Src[Stride*(-n)];
00373     
00374     /* N == 1 is a special case */
00375     if(N == 1)
00376         return Src[0];
00377     
00378     /* General formula for extension at an arbitrary n */
00379     Jump = 2*(Src[Stride*(N - 1)] - Src[0]);
00380     Offset = 0;        
00381     
00382     if(n >= N)
00383     {
00384         do
00385         {
00386             Offset += Jump;
00387             n -= 2*(N - 1);
00388         }while(n >= N);
00389     }
00390     else
00391     {
00392         while(n <= -N)
00393         {
00394             Offset -= Jump;
00395             n += 2*(N - 1);
00396         }
00397     }
00398     
00399     if(n >= 0)
00400         return Src[Stride*n] + Offset;
00401     else
00402         return 2*Src[0] - Src[Stride*(-n)] + Offset;
00403 }
00404 
00405 
00420 boundaryext GetBoundaryExt(const char *Boundary)
00421 {
00422     if(!strcmp(Boundary, "zpd") || !strcmp(Boundary, "zero"))
00423         return ZeroPaddedExtension;
00424     else if(!strcmp(Boundary, "sp0") || !strcmp(Boundary, "const"))
00425         return ConstantExtension;
00426     else if(!strcmp(Boundary, "sp1") || !strcmp(Boundary, "linear"))
00427         return LinearExtension;
00428     else if(!strcmp(Boundary, "per") || !strcmp(Boundary, "periodic"))
00429         return PeriodicExtension;
00430     else if(!strcmp(Boundary, "sym") 
00431         || !strcmp(Boundary, "symh") || !strcmp(Boundary, "hsym"))
00432         return SymhExtension;
00433     else if(!strcmp(Boundary, "symw") || !strcmp(Boundary, "wsym"))
00434         return SymwExtension;
00435     else if(!strcmp(Boundary, "asym") 
00436         || !strcmp(Boundary, "asymh") || !strcmp(Boundary, "hasym"))
00437         return AsymhExtension;
00438     else if(!strcmp(Boundary, "asymw") || !strcmp(Boundary, "wasym"))
00439         return AsymwExtension;
00440     else
00441     {
00442         ErrorMessage("Unknown boundary extension \"%s\".\n", Boundary);
00443         return NULL;
00444     }
00445 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines