apps/ccam/lasercalc.c

Go to the documentation of this file.

/*
 * lasercalc.c
 *
 */
#include <fcntl.h>  /*open*/
#include <unistd.h> /* close */
#include <string.h> /* memset */

#include <stdio.h>
#include <stdlib.h>
#include <sys/mman.h>           /* mmap */
#include <math.h>

#include "lasercalc.h"
#define D(x)


//   MD(fprintf(stderr,"width- %d, height %d, quality - %d, contrast - %d, color - %d, orient - %d, depth - %d\r\n", ImageWidth,ImageHeight, Quality, Contrast, Color, bayerOrient, Depth));

/*
struct beam_params {
// primary sums
  MYPRECISION                   SQ;             // Sum of all pixels values
  MYPRECISION                   SxQ;    // Sum of (pixels values multiplied by x)
  MYPRECISION                   SyQ;    // Sum of (pixels values multiplied by y)
  MYPRECISION                   SxxQ;   // Sum of (pixels values multiplied by x*x)
  MYPRECISION                   SyyQ;   // Sum of (pixels values multiplied by y*x)
  MYPRECISION                   SxyQ;   // Sum of (pixels values multiplied by x*y)
//momentums
  MYPRECISION                   xc;             // center x
  MYPRECISION                   yc;             // center y
  MYPRECISION                   s2x;    // (sigma_x)^2
  MYPRECISION                   s2y;    // (sigma_y)^2
  MYPRECISION                   sxy;    //  (sigma_xy)^2
//results
  MYPRECISION                   dx;             // diamerter_x
  MYPRECISION                   dy;             // diamerter_y
  MYPRECISION                   d;              // diamerter
  MYPRECISION                   az;             // ellipse azimuth
  MYPRECISION                   dxe;    // ellipse width x
  MYPRECISION                   dye;    // ellipse width y
  MYPRECISION                   ar;             // aspect ratio
};

*/
int calculateBeamParams(const char * dmaFileName,
                                        int     width,
                                        int height,
                                        struct beam_params * beam
                                        ) {
   int  i,j,k,l,y;
   int  d,q;
   int  DMALongsPerRow;
   int  data_fd;
   long long iSQ=0;
   long long iSxQ=0;
   long long iSyQ=0;
   long long iSxxQ=0;
   long long iSyyQ=0;
   long long iSxyQ=0;
   long long d2;
   int   eps;

   unsigned long        * DMAData;  // opened through mmap
   unsigned long        scols [MAX_IMAGE_WIDTH];        // sum of each column
   unsigned long        scolsy[MAX_IMAGE_WIDTH];        // sum of each column: (pixel value*y)
   unsigned long        srows [MAX_IMAGE_HEIGHT];       // sum of each row

   MD(fprintf(stderr,"dmaFileName- %s,width - %d, height %d\r\n",dmaFileName,width,height));
   if ( (width>MAX_IMAGE_WIDTH) ||
        (width<=0) ||
            (height>MAX_IMAGE_HEIGHT) ||
        (height<=0) ) return -1;
   i=width/3;
   DMALongsPerRow = (3*i == width)? i: i+1;      
//try open dmaFileName
   if ((data_fd = open(dmaFileName, O_RDONLY))==-1) return -2; // file open error
 MD(fprintf(stderr,"\r\n"));
        DMAData = (unsigned long *) mmap(0, (DMALongsPerRow*height)<<2, PROT_READ, MAP_OPTIONS, data_fd, 0);
 MD(fprintf(stderr,"\r\n"));

        if((int)(DMAData) == -1) { // Error in mmap
                close(data_fd);
                return -3;
        }
 MD(fprintf(stderr,"\r\n"));
 // prepare arrays for integration
    memset (scols, 0,sizeof(scols ));
    memset (scolsy,0,sizeof(scolsy));
    memset (srows, 0,sizeof(srows));
        for (i=0;i<height;i++) {
          l=DMALongsPerRow*i;
          d=DMAData[l++];
          k=0;
          y=height-i;
          for (j=0;j<width;j++) {
             q=(d>>k)&0x3ff;
                 scols [j]+=q;
                 scolsy[j]+=q*y;
                 srows [y]+=q;
             if ((k+=10)>20) {
                   k=0;
                   d=DMAData[l++];
                 }

          }
        }
// calculate long long (64-bit) integer sums
        for (j=0;j<width;j++) {
                iSQ+= scols[j];
                d2=((long long) j) * scols[j];
                iSxQ+=d2;
                iSxxQ+=d2 * j;
                iSxyQ+=((long long) j) * scolsy[j];
        }
        for (i=0;i<height;i++) {
                d2=((long long) i) * srows[i];
                iSyQ+=d2;
                iSyyQ+=d2 * i;
        }
// convert them to floating point format (see if single or double precision)
        beam->SQ=       iSQ;
        beam->SxQ=      iSxQ;
        beam->SyQ=      iSyQ;
        beam->SxxQ=     iSxxQ;
        beam->SyyQ=     iSyyQ;
        beam->SxyQ=     iSxyQ;
// now - center and second order moments
        beam->xc=       beam->SxQ/beam->SQ;
        beam->yc=       beam->SyQ/beam->SQ;
        beam->s2x=      beam->SxxQ/beam->SQ - (beam->xc*beam->xc);
        beam->s2y=      beam->SyyQ/beam->SQ - (beam->yc*beam->yc);
        beam->sxy=      beam->SxyQ/beam->SQ - (beam->xc*beam->yc);
// results
        beam->dx=       4*sqrt(beam->s2x);
        beam->dy=       4*sqrt(beam->s2y);
        beam->d=                4*sqrt((beam->s2x+beam->s2y)/2);
        beam->az=       0.5*atan2(2*beam->sxy,beam->s2x-beam->s2y);
        
        eps= (beam->s2x > beam->s2y) ? 1: -1;
        beam->dxe=      2.0*sqrt(2.0)*sqrt(beam->s2x+beam->s2y+\
                eps*sqrt((beam->s2x-beam->s2y)*(beam->s2x-beam->s2y) +4*(beam->sxy*beam->sxy)));
        beam->dye=      2.0*sqrt(2.0)*sqrt(beam->s2x+beam->s2y-\
                eps*sqrt((beam->s2x-beam->s2y)*(beam->s2x-beam->s2y) +4*(beam->sxy*beam->sxy)));
        beam->ar=       beam->dye/beam->dxe;
    return 0; // all OK
}

---