Prefilter

/*********************** CubicBSpline Prefilter ****************************/
327		/* This function creates a prefilter in the frequency domain. It can be
328		   viewed as a compliment to a traditional bicubic filter. */
329		
330		void create_prefilter(ImBuf *in, ImBuf *out)
331		{
332		    int x, y;
333		    
334		    // A row of float coefficients.
335		    float *coeff = out->rect_float;
336		    
337		    if (in == NULL || in->rect_float == NULL || out == NULL ||
338		        out->rect_float == NULL || in->x <= 2 || in->y <=2)
339		        return;
340		    
341		    // Initialize
342		    memcpy(coeff, in->rect_float, in->x * in->y * 4 * sizeof(float));
343		
344		    // X Recursion
345		    for(y = 0; y < in->y; y++)
346		        prefilter_doRecursion(&coeff[y * in->x * 4], in->x, 4);
347		    // Y Recursion
348		    for(x = 0; x < in->x; x++)
349		        prefilter_doRecursion(&coeff[x * 4], in->y, in->x * 4);
350		        
351		    // Copy to the out ImBuf.
352		    memcpy(out->rect_float, coeff, (in->x * in->y * 4 * sizeof(float)));
353		}
354		
355		// Get the initial causual coefficient for a given run of floats at stride count
356		// of floats.
357		void prefilter_initInitialCausal(float *coeffIn, float *sum, unsigned int length, unsigned int stride)
358		{
359		    const float Zp = sqrt(3.0f)-2.0f;
360		    const float lambda = (1.0f - Zp) * (1.0f - 1.0f / Zp);
361		    unsigned int count = 0, horizon = MIN2(12, length) * 4;
362		    float Zn = Zp;
363		    
364		    sum[0] = coeffIn[0]; // R
365		    sum[1] = coeffIn[1]; // G
366		    sum[2] = coeffIn[2]; // B
367		    sum[3] = coeffIn[3]; // A
368		    
369		    for (count = 0; count < horizon; count += stride) {
370		        sum[0] += Zn * coeffIn[count];    // R
371		        sum[1] += Zn * coeffIn[count+1];  // G
372		        sum[2] += Zn * coeffIn[count+2];  // B
373		        sum[3] += Zn * coeffIn[count+3];  // A
374		        
375		        Zn *= Zp;
376		    }
377		    
378		    coeffIn[0] = (lambda * sum[0]); // R
379		    coeffIn[1] = (lambda * sum[1]); // G
380		    coeffIn[2] = (lambda * sum[2]); // B
381		    coeffIn[3] = (lambda * sum[3]); // A
382		}
383		
384		// Set the initial anti-causual coefficient at the end of the given run of
385		// floats at stride count of floats.
386		void prefilter_initInitialAntiCausal(float *coeffIn)
387		{
388		    const float Zp  = sqrt(3.0f)-2.0f;
389		    const float iZp = (Zp / (Zp - 1.0f));
390		    
391		    coeffIn[0] = iZp * coeffIn[0]; //R
392		    coeffIn[1] = iZp * coeffIn[1]; //G
393		    coeffIn[2] = iZp * coeffIn[2]; //B
394		    coeffIn[3] = iZp * coeffIn[3]; //A
395		}
396		
397		
398		void prefilter_doRecursion(float *coeffIn, unsigned int length, unsigned int stride)
399		{
400		    const float Zp  = sqrt(3.0f)-2.0f;
401		    const float lambda = (1.0f - Zp) * (1.0f - 1.0f / Zp);
402		    float prevcoeff[4] = { 0, 0, 0, 0 };
403		    float sum[4] = { 0, 0, 0, 0 };
404		    int count = 0;
405		    unsigned int total_floats = length * 4;
406		    
407		    prefilter_initInitialCausal(coeffIn, sum, length, 4);
408		    
409		    prevcoeff[0] = coeffIn[0]; // R
410		    prevcoeff[1] = coeffIn[1]; // G
411		    prevcoeff[2] = coeffIn[2]; // B
412		    prevcoeff[3] = coeffIn[3]; // A
413		    
414		    for (count = (1 * stride); count < total_floats; count += stride) {
415		        coeffIn[count]     = prevcoeff[0] = (lambda * coeffIn[count])     + (Zp * prevcoeff[0]); // R
416		        coeffIn[count + 1] = prevcoeff[1] = (lambda * coeffIn[count + 1]) + (Zp * prevcoeff[1]); // G
417		        coeffIn[count + 2] = prevcoeff[2] = (lambda * coeffIn[count + 2]) + (Zp * prevcoeff[2]); // B
418		        coeffIn[count + 3] = prevcoeff[3] = (lambda * coeffIn[count + 3]) + (Zp * prevcoeff[3]); // A
419		    }
420		    
421		    count -= stride;
422		    prefilter_initInitialAntiCausal(&coeffIn[count]);
423		   
424		    prevcoeff[0] = coeffIn[count];     // R
425		    prevcoeff[1] = coeffIn[count + 1]; // G
426		    prevcoeff[2] = coeffIn[count + 2]; // B
427		    prevcoeff[3] = coeffIn[count + 3]; // A
428		
429		    for (count -= stride; count >= 0; count -= stride) {
430		        coeffIn[count]     = prevcoeff[0] = Zp * (prevcoeff[0] - coeffIn[count]);     // R
431		        coeffIn[count + 1] = prevcoeff[1] = Zp * (prevcoeff[1] - coeffIn[count + 1]); // G
432		        coeffIn[count + 2] = prevcoeff[2] = Zp * (prevcoeff[2] - coeffIn[count + 2]); // B
433		        coeffIn[count + 3] = prevcoeff[3] = Zp * (prevcoeff[3] - coeffIn[count + 3]); // A
434		    }
435		}