Code sample for an image rectification algorithm in C++

gistfile1.cpp

// MatrixImagingLibrary.cpp [fragment]                Alejandro Carrillo (July 2013)

// At the Neuro-Electronic Research Foundation in Belgium, I wrote image
// rectification software for their custom multiphoton microscope.
// The microscope functions by emitting photons and measuring the fluorescence
// of neurons in mice that have genetically encoded calcium indicators.
// In order to sweep the laser across the microscope subject at a maximally high
// frequency, the laser mirror is driven via harmonic resonance.
// The consequence, however, is that the laser sweep is sinusoidal. The following
// functions concern rectifying the sine-warped image that is generated from the
// raw datapoints.
//
// Note that, while the side-sweeping action of the mirror is sinusoidal in nature,
// the slower downwards sweep from one row to another is linear with time.


// ------------ VARIABLES ----------------


int *SineBinningTable; // Contains mapping from element G[i] to H[j]
int mhz;               // Side-sweep scan frequency.


// ------------ FUNCTIONS ----------------


/**
 * Generates a lookup table to map sine-warped datapoints to a normalized image.
 *
 * Suppose we have pixel line G containing n raw datapoints. G[i] is remapped to
 * H[j] by letting j = SineBinningTable[i]. This function generates SineBinningTable.
 *
 * dsratio: the down sampling ratio, if we want to lower the resolution of our output.
 *
 * TODO: Make this work even if user changes scan rate in the middle of experimenting.
 */
void MatroxImagingLibrary::GenerateBinningTable(int dsratio) {
    int lineWidthOut; // Width of sine-corrected image
    double deltaPhase;
    int nSamplesPerCycle;// Width of raw image. Equivalent to one scan period.
    std::ofstream myfile;
    myfile.open ("binning_tbl.txt"); // For debugging


    if (mhz == 20) {
        nSamplesPerCycle = 2308;
    } else if (mhz == 80) {
        nSamplesPerCycle = 9000;
    } else {
        printf("Failed to make binning table! Scan frequency was neither 80 MHz nor 20 MHz");
        return;
    }

    assert(dsratio == 1 || dsratio == 2 || dsratio == 4); // dsratio should be 1, 2, or 4!

    switch (nSamplesPerCycle) {
        case 9000: // 80 MHz
            lineWidthOut = 2880/dsratio; // Chosen so no sample discarded at FOV center
            break;
        case 2308: // 20 MHz
            lineWidthOut = 720/dsratio;
            break;
        default:
            printf("Generating Binning Table failed!");
            break;
    }


    // phase increment per sample, where nSamplesPerCycle == one period.
    deltaPhase = 2.0*M_PI / nSamplesPerCycle;
    SineBinningTable = new int[nSamplesPerCycle / 2]; // Integer division (OK)

    // Populate binning table
    double phi;//phases, incremented from 0 to pi
    for (int i = 0; i < nSamplesPerCycle / 2; i++)
    {
        phi = i * deltaPhase; // current phase angle
        SineBinningTable[i] = (int) ((1 - cos(phi))*(lineWidthOut-1)/2.0) + 1;
        //printf("Binning Table value:\t%d\n",SineBinningTable[i]);
        myfile << SineBinningTable[i];
        myfile << "\n";
    }

    nBinnedWidth = lineWidthOut;

    printf("Generated Binning Table.\n");
    myfile.close();
}

/**
 * Transforms pixels from space *src to *dest via the line mapping in SineBinningTable.
 * Each raw line in *src actually contains two periods of a sine wave. This is
 * because the frame grabber sweeps back and forth in a full cycle before making a
 * new line in the raw output image. Thus each line must be broken into two before
 * the datapoints in each line are then re-mapped according to SineBinningTable.
 *
 * Once raw pixels have been remapped to *dest, pixels that have been binned into the
 * same index are averaged.
 *
 * srcHeight: Number of raw frame grabber lines in the image
 * srcWidth:  Total number of datapoints in the raw line, which is actually two lines in one
 * *src:      Source 2D array. Rows are sine-warped; columns are not.
 * *dest:     Output array.
 */
void MatroxImagingLibrary::BinningTransform(int srcHeight, int srcWidth, byte *src, byte *dest) {

    printf("In BinningTransform!\n");
    int numInputCycles, numOutputLines;
    int delay; // Inter-line delay
    int nSamplesPerInputLine;
	long int dest_ind1, dest_ind2;
    int *pos; // A buffer cursor, with indices to the start of each row.
    unsigned short count;


    numInputCycles=srcHeight;          // number of input cycles
    numOutputLines = srcHeight * 2;    // number of output lines
    nSamplesPerInputLine = srcWidth/2; // Since each line of length srcWidth is actually two real image lines

    pos = new int[numInputCycles];

    for (int i = 0; i < numInputCycles;i++) {
        pos[i] = i*srcWidth;
    }


    count = 0;
    delay = 1; //TODO: Make this an adjustable paramater!
    int srcval1, srcval2, destval1, destval2;

    std::fill_n(dest,nBinnedWidth*numOutputLines,0); // Initialize frame buffer to all zeros

    // loop over raw samples
    for(int j=0;j<nSamplesPerInputLine-1;j++) {

        //loop over lines
        for(int i=0;i<numInputCycles;i++) {
            // bin samples odd lines (left-to-right part of the sidesweep)
			dest_ind1 = SineBinningTable[j]+2*nBinnedWidth*i;
            srcval1 = src[pos[i]+j+delay];
            dest[dest_ind1] += srcval1;

            // bin samples even lines (right-to-left part of the sidesweep)
			dest_ind2 = (nBinnedWidth-(SineBinningTable[j]+1))+2*i*nBinnedWidth+1+(SineBinningTable[j]*2);
			srcval2 = src[pos[i] + srcWidth -1 - j];
            dest[dest_ind2] += srcval2;

		}

        count = count + 1;

        // compute average for each bin
        if (j + 1 < nSamplesPerInputLine && SineBinningTable[j+1]>SineBinningTable[j] && count > 0)
        {
            if (count>1)
            {
                 for(int i=0;i<numInputCycles;i++)
                 {
					 dest_ind1 = SineBinningTable[j]+2*nBinnedWidth*i;
					 dest_ind2 = (nBinnedWidth-(SineBinningTable[j]+1))+2*i*nBinnedWidth+1+(SineBinningTable[j]*2);

                     dest[dest_ind1] = dest[dest_ind1] / count;
                     dest[dest_ind2] = dest[dest_ind2] / count;
                 }
            }
             count = 0;
        }
    }

    int j = nSamplesPerInputLine-1;

    // Just makes sure to average the binned contents of the last column
    // Not exactly sure why I needed this special case outside of the nested for loops above -- I didn't write a comment
    // over the summer about this block...but the code is tested and works.
    if (count>1)
    {
        for(int i=0;i<numInputCycles;i++)
        {
					 dest_ind1 = SineBinningTable[j]+2*nBinnedWidth*i;
					 dest_ind2 = (nBinnedWidth-(SineBinningTable[j]+1))+2*i*nBinnedWidth+1+(SineBinningTable[j]*2);
                     dest[dest_ind1] = dest[dest_ind1] / count;
                     dest[dest_ind2] = dest[dest_ind2] / count;
         }
     }


    delete [] pos;
    return;
}