sicspsi/autowin.c

/**
 * This is a routine which calculates automatic windows for RITA's
 * monochromatic window mode. In this mode data from 9 analyzer blades 
 * is displayed on the detector. RITA needs the windows in which to 
 * sum the contributions of each detector blade. The algorithm to find 
 * them goes like this:
 * 1) the HM is summed in y, resulting in a 1D histogram
 * 2) In the 1D histogram peaks are located.
 * 3) For each peak from 2, the COG is calculated to get a good peak
 *    center
 * 4) The window is the peak center +- 9 bixels.
 * 
 * copyright: see file COPYRIGHT
 * 
 * Mark Koennecke, July 2009
 */
#include <math.h>
#include <sics.h>
#include <sicsobj.h>
#include <HistMem.h>
#include <HistMem.i>
#include <HistDriv.i>
#include <lld.h>
#include <dynstring.h>
#include <sicshipadaba.h>

extern float nintf(float x);
/*------------------------------------------------------------------------*/
static HistInt *sumY(pHistMem hm, SConnection *pCon)
{
	int dim[10], i;
	char command[30];
	HistInt *data = NULL;
	
	GetHistDim(hm,dim,&i);
	snprintf(command,30,"sum:0:0:%d", dim[1]);
	data = hm->pDriv->SubSample(hm->pDriv, pCon, 0, command);
	/*
	for(i = 0; i < dim[0]; i++){
		printf(" %d %d\n", i, data[i] );
	}
	*/
	return data;
}
/*-----------------------------------------------------------------------*/
static int findPeaks(HistInt *data, int length){
	int result, i, j;
	float sum, avg;
	
	result = LLDcreate(sizeof(int));
	
	for(i = 3; i < length - 3; i++){
		sum = 0.;
		for(j = i -3; j < i +3; j++){
			if(j != 4){
				sum += data[j];
			}
		}
		avg = sum/6.;
		if(data[i] > 8.*sqrt(avg)){
			if(data[i+1] < data[i] && data[i] > data[i-1]){
				LLDnodeAppendFrom(result,&i);
			}
		}
	}
	
	return result;
}
/*------------------------------------------------------------------------*/
static float calcCOG(HistInt *data, int length, int startPos)
{
	float threshold, sum = .0, xsum = .0, cog;
	int i;
	
	threshold = data[startPos] *  .33;
	
	/* sum to left */
	i = startPos;
	while(i > 0 && data[i] > threshold){
		sum += data[i];
		xsum += data[i]*i;
		i--;
	}
	
	/* sum to right */
	i = startPos+1;
	while(i < length && data[i] > threshold){
		sum += data[i];
		xsum += data[i]*i;
		i++;
	}
	
	cog = xsum/sum;
	return cog;
}
/*-----------------------------------------------------------------------*/
static void printPeaks(SConnection *pCon, int peakList, 
		HistInt *data, int length)
{
  int status, startPos, oldStartPos = 0;
  float cog;
  pDynString buf;
  char buffer[512];

  buf = CreateDynString(256,256);
  if(buf == NULL){
    SCWrite(pCon,"ERROR: out of memory printing peaks", eError);
    return;
  }
  status = LLDnodePtr2First(peakList);
  while(status > 0){
    LLDnodeDataTo(peakList,&startPos);
    cog = calcCOG(data, length, startPos);
    startPos = (int)nintf(cog);
    if((startPos - oldStartPos) > 2){
      snprintf(buffer,512,"%d,%d,%d",startPos - 5, startPos + 5, startPos);
      DynStringConcatLine(buf,buffer);
      oldStartPos = startPos;
    }
    status = LLDnodePtr2Next(peakList);
  }
  SCWrite(pCon,GetCharArray(buf), eValue);
  DeleteDynString(buf);
}
/*------------------------------------------------------------------------*/
static int CalcCmd(pSICSOBJ self, SConnection *pCon, pHdb commandNode,
                  pHdb par[], int nPar)
{
	
	HistInt *data = NULL;
	pHistMem hm = NULL;
	int peakList, dim[10];
	char hmName[30];
	
	if(nPar < 1){
		strcpy(hmName,"hm");
	} else {
		strlcpy(hmName,par[0]->value.v.text,30);
	}
	
	hm  = FindCommandData(pServ->pSics, 
			hmName, "HistMem");
	if(hm == NULL || hm->pDriv->SubSample == NULL){
		SCPrintf(pCon,eError,"ERROR: hm %s not found", hmName);
		return 0;
	}

	data = sumY(hm,pCon);
	if(data == NULL){
		SCWrite(pCon,"ERROR: subsampling failed", eError);
		return 0;
	}
	
	GetHistDim(hm,dim,&peakList);
	
	peakList = findPeaks(data, dim[0]);
	
	printPeaks(pCon,peakList, data, dim[0]);
	free(data);
	
	return 1;
}
/*----------------------------------------------------------------*/
int MakeRitaWin(SConnection *pCon, SicsInterp *pSics, void *pData, 
		int argc, char *argv[])
{
	pSICSOBJ pNew = NULL;
	pHdb cmd, node;
	
	pNew = SetupSICSOBJ(pCon,pSics,pData,argc,argv);
	if(pNew == NULL){
		return 0;
	}
	
	cmd = AddSICSHdbPar(pNew->objectNode, "calc", usSpy, 
			MakeSICSFunc(CalcCmd));
	SetHdbProperty(cmd,"type","command");
	SetHdbProperty(cmd,"priv","spy");
	node = MakeSICSHdbPar("hm",usSpy,MakeHdbText("banana"));
	AddHipadabaChild(cmd,node,NULL);
	
	return 1;
}