The
Design and Implementation
of a
Lens Flare Filter

 

flare: noun — a sudden blaze of bright light

What follows is the description of the algorithm used to create the image above. So, how does it work? As you can see in the image, the light is concentrated around spikes (or spokes, if you will). The first thing we have to do is to describe the flare. In the demonstration program, flare.c (see end of page), the array spike_angles contains the angles expressed in radians of the spikes. The array spike_lengths contains the lengths of the spikes. So both arrays have the same size. We also need to specify the center of the flare. In flare.c, this is done with CENTERX and CENTERY. Now we are ready to render the flare. We calculate the color for each pixel (x,y) by following these steps:

  1. Find between which spikes the pixel (x,y) lies. We can do that with the C function atan2(y,x) (yes, that's (y,x) and NOT (x,y)). Given a pixel (x,y), atan2(y,x) returns the arctangent of y/x in the range -π to π radians. So, if y is negative, atan2 will return a negative value. Once we have the arctangent of pixel (x,y), we can search in the array spike_angles for the first spike with an angle smaller than the arctangent of y/x. Let's call the first spike A and the second spike B (the second spike is simply the next spike in the array).

  2. Calculate the distance from (x,y) to the spikes. To calculate the distance from a point to a line, we need five points: the point (x,y) and two lines, (x1,y1)-(x2,y2) and (x3,y3)-(x4,y4). The points (x1,y1) and (x3,y3) are the center (CENTERX,CENTERY). Point (x2,y2) must lie on spike A: (x1+cos(green angle),y1+sin(green angle)). Point (x4,y4) must lie on spike B: (x3+cos(red angle),y3+sin(red angle)). The picture below will make things clear. The formula to calculate the distance can be found in the function distance_point_spike(). Once we have the distance from pixel (x,y) to a spike, we convert it to range [0,1]. How do we do that? A spike “emits” light: the further away from the spike, the darker the color of the pixel. Suppose a spike has a “width” of 16: if the distance from pixel (x,y) to that spike is 20 then the color of that spike will not contribute to the color of the pixel. So, converting the distance to range [0,1] is nothing more than dividing the distance by the width of the spike. The color of pixel (x,y) can then be calculated by summing the light contributed by spike A and B: color(x,y) = (distance_A × color_flare) + (distance_B × color_flare).

    If we would have a flare with four spikes with a width of 16, we would end up with this image:

    As you can see, the joints could be better. We can enhance them by simply squaring the distances. This way, the light emitted by a spike will decrease in intensity quadraticly instead of linearly:

    Here is the new image:

  3. Up to now, we didn't take the spike lengths into account. As you can see in the image at the top of the page, the light intensity of a pixel decreases when it's further away from the center. We could use the same trick as we did with the distance from a pixel to a spike: calculating the distance and dividing it by a maximum length to convert it to range [0,1]. But the big question is: “What is that maximum length?” A flare doesn't have a radius: the spikes can have different lengths. So, the answer is very simple: we will have to invent our own maximum length. Now, suppose that spike A has a length of 4, spike B a length of 10 and that the angle between the two spikes is 30°. Also suppose that the angle from spike A to line (x,y)-(CENTERX,CENTERY) is 3° (the blue line between spike A and B). We know that the maximum length is 4 if pixel (x,y) is on spike A (a pixel is on spike A if the angle between the blue line and spike A is 0; remember that we can find the angle of the blue line using atan2) and 10 if it's on spike B (a pixel is on spike B if the angle of the blue line minus the angle of spike A equals the angle between spike A and spike B). And now we can use interpolation to find the maximum length for pixel (x,y) on the blue line:

     30° ->   6 (10 - 4)
      1° -> 1/5 (6 / 30)
      3° -> 3/5 (3 * 1/5)
    

    Thus we must divide the length of vector (x,y) by 4.6 (= 4 + 3/5; we have to add 4 because that is the minimum value). After dividing the length of vector (x,y) by that maximum length, we square the result again for esthetical reasons. Then we subtract that value from 1 (because the light is brighter towards the center) and multiply the result with the color value obtained in step 2. In case of a flare with four spikes, we would get this image:

    Believe it or not, but we're done. When we render a flare with 56 spikes and random angles and lengths using this algorithm, we obtain this image:

    We have seen one method to render a flare. Now it's up to you to add color, halos and other effects!

Below is the source code of flare.c. Experiment with different values for MAX_SPIKES and SPIKEWIDTH.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#define WIDTH         400
#define HEIGHT        400
#define MAX_SPIKES    56
#define RADIUS        190
#define CENTERX       200
#define CENTERY       200
#define SPIKEWIDTH    16

int find_spike_index(double spike_angles[], int x, int y) {
    int i;
    double a=atan2(y-CENTERY,x-CENTERX);
    if (a<0) a+=M_PI*2.0;              /* range [0-2*PI] */
    for(i=0;i<MAX_SPIKES;i++) {
        if (a<=spike_angles[i])
            break;
    }
    if (i==0)
        return MAX_SPIKES-1;
    return i-1;
}

double distance_point_spike(double spike_angle, int x, int y) {
    double x0=CENTERX;
    double y0=CENTERY;
    double x1=x0+cos(spike_angle);
    double y1=y0+sin(spike_angle);
    double d=fabs((x1-x0)*(y0-y)-(x0-x)*(y1-y0));
    d/=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1- y0));
    d/=SPIKEWIDTH;                            /* range [0-1] */
    if (d>1.0) d=1.0;                         /* maximum 1.0 */
    d=1.0-d;                       /* 1.0 is on spike, not 0 */
    return d;
}

double vectorlen(int x, int y) {
    return sqrt((CENTERX-x)*(CENTERX-x)+(CENTERY-y)*(CENTERY-y));
}

double maxvectorlen(int x, int y, double a1, double a2,
                    double l1, double l2) {
    double a, vlen;
    a=atan2(y-CENTERY,x-CENTERX);
    if (a<0) a+=M_PI*2.0;                  /* range [0-2*PI] */
    if (a2<a1) a2+=2*M_PI;    /* in case spike 2 has index 0 */
    if (a<a1) a+=2*M_PI;
    vlen=l1+((a-a1)/(a2-a1))*(l2-l1);
    return vlen;
}

int main(void) {
    unsigned char *image;
    FILE *f;
    unsigned short ws, hs;
    int i, x, y, index1, index2;
    double spike_angles[MAX_SPIKES];
    double spike_lengths[MAX_SPIKES];
    double d1, d2, dr, dg, db, l, ml;
    for(i=0;i<MAX_SPIKES;i++) {
        d1=(2.0*M_PI)/MAX_SPIKES;
        d2=(((double)rand())/RAND_MAX);
        spike_angles[i]=(d1*i)+(d1/2.0)*d2;
        spike_lengths[i]=((((double)rand())/RAND_MAX)+0.2)*RADIUS;
        if (spike_lengths[i]>RADIUS) spike_lengths[i]=RADIUS;
    }
    image=(unsigned char*)malloc(sizeof(unsigned char)*WIDTH*HEIGHT*3);
    memset(image,0,sizeof(unsigned char)*WIDTH*HEIGHT*3);
    x=0;y=0;
    for(i=0;i<WIDTH*HEIGHT;i++) {
        if (!(x==CENTERX && y==CENTERY)) {   /* prevent domain error */
            index1=find_spike_index(spike_angles,x,y);
            d1=distance_point_spike(spike_angles[index1],x,y);
            d1*=d1;                              /* linear is boring */
            index2=index1+1;
            if (index2>=MAX_SPIKES) index2=0;
            d2=distance_point_spike(spike_angles[index2],x,y);
            d2*=d2;                              /* linear is boring */
            l=vectorlen(x,y);
            ml=maxvectorlen(x,y,spike_angles[index1],
                            spike_angles[index2],spike_lengths[index1],
                            spike_lengths[index2]);
            l/=ml;
            if (l>1.0) l=1.0;
            l=1.0-l;
            l*=l;
        } else {
            d1=1.0;                      /* full intensity at center */
            d2=1.0;
            l=1.0;
        }
        dr=(d1*255.0)+(d2*255.0);
        dg=(d1*255.0)+(d2*255.0);
        db=(d1*255.0)+(d2*255.0);
        dr*=l;
        dg*=l;
        db*=l;
        if (dr>255.0) dr=255.0;
        if (dg>255.0) dg=255.0;
        if (db>255.0) db=255.0;
        image[i*3+0]=(unsigned char)dr;
        image[i*3+1]=(unsigned char)dg;
        image[i*3+2]=(unsigned char)db;
        x++;
        if (x%WIDTH==0) {x=0;y++;}
    }
    f=fopen("flare.tga","wb");
    ws = WIDTH;
    hs = HEIGHT;
    fwrite("\0\0\2\0\0\0\0\0\0\0\0\0", sizeof(char), 12, f);
    fwrite(&ws, sizeof(unsigned short), 1, f);
    fwrite(&hs, sizeof(unsigned short), 1, f);
    fwrite("\x18\x20", sizeof(char), 2, f);
    for (i=0;i<WIDTH*HEIGHT;i++) {
        fputc(image[i*3+2],f);
        fputc(image[i*3+1],f);
        fputc(image[i*3+0],f);
    }
    fclose(f);
    return 0;
}


Download flare.zip

 

Back