354 lines
9.5 KiB
C
354 lines
9.5 KiB
C
/*---------------------------------------------------------------------------
|
|
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);
|
|
|
|
}
|