Initial revision
This commit is contained in:
cvs
353
integrate.c
Normal file
353
integrate.c
Normal file
@ -0,0 +1,353 @@
|
||||
/*---------------------------------------------------------------------------
|
||||
P E A K I N T E G R A T E
|
||||
|
||||
Scan Peak integration along the scheme described by Grant & Gabe in
|
||||
J. Appl. Cryst (1978), 11, 114-120
|
||||
|
||||
|
||||
Mark Koennecke, March 1998
|
||||
|
||||
--------------------------------------------------------------------------*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <assert.h>
|
||||
#include <math.h>
|
||||
#include "fortify.h"
|
||||
#include "integrate.h"
|
||||
|
||||
/* define PRINT if you want more info about the peak limits and what is
|
||||
in between
|
||||
|
||||
#define PRINT 1
|
||||
*/
|
||||
|
||||
#define PROBABILITY 0.01
|
||||
#define LEFT -1
|
||||
#define RIGHT 1
|
||||
|
||||
#define ABS(x) (x < 0 ? -(x) : (x))
|
||||
|
||||
/*----------------- a private data structure to keep track of things ----*/
|
||||
|
||||
typedef struct {
|
||||
long *lCounts;
|
||||
int iCounts;
|
||||
int m;
|
||||
int iLeft;
|
||||
int iMaxCount;
|
||||
int iRight;
|
||||
} IntegData, *pIntegData;
|
||||
|
||||
/*-------------------------------------------------------------------------*/
|
||||
static int FindMaxCount(long lCounts[], int iCounts)
|
||||
{
|
||||
int i, iPtr;
|
||||
long l = -100000;
|
||||
|
||||
for(i = 0, iPtr = 0; i < iCounts;i++)
|
||||
{
|
||||
if(lCounts[i] > l)
|
||||
{
|
||||
l = lCounts[i];
|
||||
iPtr = i;
|
||||
}
|
||||
}
|
||||
return iPtr;
|
||||
}
|
||||
/*------------------------------------------------------------------------
|
||||
Finds the index where the counts are less then half peak height. Assumes
|
||||
that iMaxCount has been set to a proper value beforehand. This helps to
|
||||
solve a problem: Given very broad peaks, the window condition would be
|
||||
fulfilled when starting for background counts from the top. The solution
|
||||
is to search for the half width point and start searchnig for background
|
||||
from there. Saves calculation time as well. FindHalf returns an index on
|
||||
success, else - 1 if it goes beyond limits!
|
||||
-----------------------------------------------------------------------*/
|
||||
static int FindHalf(pIntegData self, int iSign)
|
||||
{
|
||||
int i;
|
||||
long lHalf;
|
||||
|
||||
assert(self);
|
||||
|
||||
lHalf = self->lCounts[self->iMaxCount]/2;
|
||||
i = self->iMaxCount;
|
||||
while(self->lCounts[i] > lHalf)
|
||||
{
|
||||
i += iSign;
|
||||
if( (i < 0) || (i > self->iCounts) )
|
||||
{
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
return i;
|
||||
}
|
||||
/*--------------------------------------------------------------------------
|
||||
FindWindowHits find the number of hits in the window defined by
|
||||
0.67 sigma of the value at iPos and m. Returns -1 if insufficient data
|
||||
is available to do the trick. iSign can be -1 for search to left or
|
||||
1 for search to right.
|
||||
---------------------------------------------------------------------------*/
|
||||
static int FindWindowHits(pIntegData self, int iPos, int iSign)
|
||||
{
|
||||
float fWinMin, fWinMax, fVar;
|
||||
int i, iTest, iCount;
|
||||
|
||||
assert(self);
|
||||
assert( (iSign == -1) || (iSign == 1) );
|
||||
|
||||
/* check limits of the show */
|
||||
iTest = iPos + iSign * self->m;
|
||||
if( (iTest < 0) || (iTest >= self->iCounts) )
|
||||
{
|
||||
return -1;
|
||||
}
|
||||
|
||||
/* calculate limits */
|
||||
if(self->lCounts[iPos] > 0)
|
||||
{
|
||||
fVar = sqrt((double)self->lCounts[iPos]);
|
||||
fWinMin = (float)self->lCounts[iPos] - 0.67 * fVar;
|
||||
fWinMax = (float)self->lCounts[iPos] + 0.67 * fVar;
|
||||
}
|
||||
else
|
||||
{
|
||||
fWinMin = 0.;
|
||||
fWinMax = 1.;
|
||||
}
|
||||
|
||||
/* count the window hits */
|
||||
for(i = 0, iCount = 0; i < self->m; i++)
|
||||
{
|
||||
iTest = iPos + iSign * i;
|
||||
if( (self->lCounts[iTest] >= fWinMin) &&
|
||||
(self->lCounts[iTest] <= fWinMax) )
|
||||
{
|
||||
iCount++;
|
||||
}
|
||||
}
|
||||
return iCount;
|
||||
}
|
||||
/*----------------------------------------------------------------------------
|
||||
AdvanceBackground advances in the background until the predefined
|
||||
Probability has been reached. The probability is defined at the top
|
||||
of the file. iSign and errors are the same as above.
|
||||
----------------------------------------------------------------------------*/
|
||||
static int AdvanceBackground(pIntegData self, int iPos, int iSign)
|
||||
{
|
||||
float fProb;
|
||||
int iRet, iTest;
|
||||
|
||||
assert(self);
|
||||
assert( (iSign == 1) || (iSign == -1) );
|
||||
|
||||
/* initialize, calculate initial Probability */
|
||||
iRet = FindWindowHits(self,iPos,iSign);
|
||||
if(iRet < 0)
|
||||
{
|
||||
return iRet;
|
||||
}
|
||||
fProb = (self->m - iRet)/self->m;
|
||||
iTest = iPos;
|
||||
|
||||
/* loop further on */
|
||||
while(fProb > PROBABILITY)
|
||||
{
|
||||
iTest += iSign;
|
||||
iRet = FindWindowHits(self,iTest,iSign);
|
||||
if(iRet < 0)
|
||||
{
|
||||
return iRet;
|
||||
}
|
||||
fProb *= (self->m - iRet)/self->m;
|
||||
}
|
||||
|
||||
if( (iTest < 0) || (iTest >= self->iCounts) )
|
||||
{
|
||||
return -1;
|
||||
}
|
||||
else
|
||||
{
|
||||
return iTest;
|
||||
}
|
||||
}
|
||||
/*--------------------------------------------------------------------------
|
||||
FindLimit assumes, that iMaxCount has been set to something sensible
|
||||
in advance!
|
||||
--------------------------------------------------------------------------*/
|
||||
static int FindLimit(pIntegData self, int iStart, int iSign)
|
||||
{
|
||||
int i, iHalf, iHit, iLeft, iTest;
|
||||
|
||||
assert(self);
|
||||
assert( (iSign == 1) || (iSign == -1) );
|
||||
|
||||
iHit = 0;
|
||||
iHalf = self->m/2;
|
||||
i = iStart;
|
||||
|
||||
/* loop until iHalf is reached */
|
||||
while((iHit = FindWindowHits(self,i,iSign)) < iHalf)
|
||||
{
|
||||
i += iSign;
|
||||
if(iHit < 0)
|
||||
{
|
||||
if(iSign == LEFT)
|
||||
{
|
||||
return INTEGLEFT;
|
||||
}
|
||||
else
|
||||
{
|
||||
return INTEGRIGHT;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* step further into the backgound */
|
||||
iTest = AdvanceBackground(self,i,iSign);
|
||||
if(iTest < 0)
|
||||
{
|
||||
if(iSign == LEFT)
|
||||
{
|
||||
return INTEGLEFT;
|
||||
}
|
||||
else
|
||||
{
|
||||
return INTEGRIGHT;
|
||||
}
|
||||
}
|
||||
|
||||
/* store it away and return success */
|
||||
if(iSign == LEFT)
|
||||
{
|
||||
self->iLeft = iTest;
|
||||
}
|
||||
else
|
||||
{
|
||||
self->iRight = iTest;
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
/*--------------------------------------------------------------------------
|
||||
This routine does the actual integration once the limits have been
|
||||
determined.
|
||||
--------------------------------------------------------------------------*/
|
||||
static int DoIntegrate(pIntegData self, float *fSum, float *fVariance)
|
||||
{
|
||||
long lLeftBack, lRightBack, lBackMean;
|
||||
int i,iRes = 1;
|
||||
float fScan, fBackVar,fLeft, fRight, fN, fNP, fNB, fSumm;
|
||||
float fMLeft, fMRight;
|
||||
|
||||
assert(self);
|
||||
assert( (self->iLeft != 0) && (self->iRight != 0) );
|
||||
|
||||
|
||||
/* sum the left background */
|
||||
for(i = 0, lLeftBack = 0; i < self->iLeft; i++)
|
||||
{
|
||||
lLeftBack += self->lCounts[i];
|
||||
}
|
||||
lLeftBack = lLeftBack;
|
||||
|
||||
/* sum the right background */
|
||||
for(i = self->iRight, lRightBack = 0; i < self->iCounts; i++)
|
||||
{
|
||||
lRightBack += self->lCounts[i];
|
||||
}
|
||||
lRightBack = lRightBack;
|
||||
|
||||
|
||||
/* sum all of it */
|
||||
fScan = 0.;
|
||||
for(i = 0; i < self->iCounts; i++)
|
||||
{
|
||||
fScan += self->lCounts[i];
|
||||
}
|
||||
|
||||
/* sum the peak */
|
||||
fSumm = 0.;
|
||||
for(i = self->iLeft; i < self->iRight +1; i++)
|
||||
{
|
||||
fSumm += (float)self->lCounts[i];
|
||||
}
|
||||
|
||||
/* calculate the values */
|
||||
fLeft = (float)self->iLeft;
|
||||
fRight = (float)self->iRight;
|
||||
fN = (float)self->iCounts;
|
||||
|
||||
fNP = fRight - fLeft + 1.;
|
||||
fNB = fLeft + (fN - fRight);
|
||||
|
||||
*fSum = fSumm - (fNP/fNB)*((float)lLeftBack + (float)lRightBack);
|
||||
*fVariance = sqrt(fSumm + (fNP/fNB)*(fNP/fNB)*((float)lLeftBack
|
||||
+ (float)lRightBack));
|
||||
|
||||
/* check the background */
|
||||
fMLeft = lLeftBack/fLeft;
|
||||
fMRight = lRightBack/fRight;
|
||||
if(ABS(fMLeft - fMRight) > ((fMLeft + fMRight)/2) )
|
||||
{
|
||||
iRes = INTEGFUNNYBACK;
|
||||
}
|
||||
|
||||
return iRes;
|
||||
}
|
||||
/*---------------------------------------------------------------------------*/
|
||||
int GabePeakIntegrate(int m, int iCounts, long lCounts[],
|
||||
float *fIntens, float *fVariance)
|
||||
{
|
||||
IntegData self;
|
||||
int iTest,i, iStart;
|
||||
|
||||
/* check input */
|
||||
assert(iCounts > 0);
|
||||
assert(lCounts);
|
||||
|
||||
|
||||
/* initialize Data Structure */
|
||||
self.iCounts = iCounts;
|
||||
self.lCounts = lCounts; /* no need to duplicate data */
|
||||
self.iLeft = 0;
|
||||
self.iRight = 0;
|
||||
self.iMaxCount = 0;
|
||||
self.m = m;
|
||||
|
||||
/* find maximum */
|
||||
self.iMaxCount = FindMaxCount(self.lCounts, self.iCounts);
|
||||
|
||||
/* find Left Limit */
|
||||
iStart = FindHalf(&self,LEFT);
|
||||
if(iStart < 0)
|
||||
{
|
||||
return INTEGLEFT;
|
||||
}
|
||||
iTest = FindLimit(&self,iStart,LEFT);
|
||||
if(iTest < 0)
|
||||
{
|
||||
return iTest;
|
||||
}
|
||||
|
||||
/* find right limit */
|
||||
iStart = FindHalf(&self,RIGHT);
|
||||
if(iStart < 0)
|
||||
{
|
||||
return INTEGRIGHT;
|
||||
}
|
||||
iTest = FindLimit(&self,iStart,RIGHT);
|
||||
if(iTest < 0)
|
||||
{
|
||||
return iTest;
|
||||
}
|
||||
#ifdef PRINT
|
||||
printf("Limits: %d %d\n",self.iLeft, self.iRight);
|
||||
for(i = self.iLeft; i < self.iRight; i++)
|
||||
{
|
||||
printf(" %d\n",self.lCounts[i]);
|
||||
}
|
||||
#endif
|
||||
/* integrate */
|
||||
return DoIntegrate(&self,fIntens, fVariance);
|
||||
|
||||
}
|
||||
Reference in New Issue
Block a user