322 lines
8.2 KiB
C
322 lines
8.2 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 = 0, 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];
|
|
}
|
|
|
|
/* sum the right background */
|
|
for (i = self->iRight, lRightBack = 0; i < self->iCounts; i++) {
|
|
lRightBack += self->lCounts[i];
|
|
}
|
|
|
|
|
|
/* 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);
|
|
|
|
}
|