mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-04-23 15:00:02 +02:00
f6e76145c1
* removed Makefile for moench and integrated the build in CMake * broke out tiff reading and writing to its own library * moved tiff includes to include/sls * moved tiffio source to src * removed incorrectly used bps * cleanup and tests for tiffio * removed using namespace std from header * some fixing for moench04 * Program for offline processing renamed Co-authored-by: Anna Bergamaschi <anna.bergamaschi@psi.ch>
400 lines
14 KiB
C++
400 lines
14 KiB
C++
// SPDX-License-Identifier: LGPL-3.0-or-other
|
|
// Copyright (C) 2021 Contributors to the SLS Detector Package
|
|
#ifndef ETA_INTERPOLATION_RANDOMBINS_H
|
|
#define ETA_INTERPOLATION_RANDOMBINS_H
|
|
|
|
#include "sls/tiffIO.h"
|
|
//#include "etaInterpolationBase.h"
|
|
#include "etaInterpolationPosXY.h"
|
|
#include <algorithm>
|
|
#include <cstdlib>
|
|
//#include <math>
|
|
#include <cmath> // std::abs
|
|
|
|
#define PI 3.14159265
|
|
#define TWOPI 2. * PI
|
|
|
|
using namespace std;
|
|
|
|
class etaInterpolationRandomBins : public etaInterpolationPosXY {
|
|
|
|
private:
|
|
double calcDiff(double avg, float *hx, float *hy) {
|
|
double p_tot = 0;
|
|
double diff = 0;
|
|
double bsize = 1. / nSubPixels;
|
|
|
|
for (int ipx = 0; ipx < nSubPixels; ipx++) {
|
|
for (int ipy = 0; ipy < nSubPixels; ipy++) {
|
|
p_tot = 0;
|
|
for (int ibx = 0; ibx < nbeta; ibx++) {
|
|
for (int iby = 0; iby < nbeta; iby++) {
|
|
if (hx[ibx + iby * nbeta] >= ((ipx)*bsize) &&
|
|
hx[ibx + iby * nbeta] < ((ipx + 1) * bsize) &&
|
|
hy[ibx + iby * nbeta] >= ((ipy)*bsize) &&
|
|
hy[ibx + iby * nbeta] < ((ipy + 1) * bsize)) {
|
|
p_tot += heta[ibx + iby * nbeta];
|
|
}
|
|
}
|
|
}
|
|
|
|
// cout << p_tot << " \t ";
|
|
|
|
diff += (p_tot - avg) * (p_tot - avg);
|
|
}
|
|
// cout << "\n";
|
|
}
|
|
return diff;
|
|
}
|
|
|
|
double iterate(float *newhhx, float *newhhy, double avg) {
|
|
|
|
double bsize = 1. / nSubPixels;
|
|
|
|
double hy[nbeta]; // profile y
|
|
double hx[nbeta]; // profile x
|
|
double hix[nbeta]; // integral of projection x
|
|
double hiy[nbeta]; // integral of projection y
|
|
|
|
double tot_eta_x = 0;
|
|
double tot_eta_y = 0;
|
|
|
|
int p0;
|
|
int vx[(nSubPixels + 1) * (nSubPixels + 1)],
|
|
vy[(nSubPixels + 1) * (nSubPixels + 1)];
|
|
|
|
int arrx[nSubPixels + 1], arry[nSubPixels + 1];
|
|
|
|
int bad = 1;
|
|
|
|
int isby, isbx;
|
|
int ii = 0;
|
|
|
|
// using default comparison (operator <):
|
|
// std::sort (myvector.begin(), myvector.begin()+4); //(12 32
|
|
// 45 71)26 80 53 33
|
|
|
|
for (isby = 0; isby < (nSubPixels + 1) / 2 + 1; isby++) {
|
|
|
|
for (isbx = 0; isbx < (nSubPixels + 1) / 2 + 1; isbx++) {
|
|
p0 = isby * (nSubPixels + 1) + isbx;
|
|
// for (int iv=0; iv<(nSubPixels+1)*(nSubPixels+1); iv++) {
|
|
if (isbx == 0) {
|
|
vy[p0] = isby * nbeta / nSubPixels;
|
|
vx[p0] = 0;
|
|
} else if (isby == 0) {
|
|
vy[p0] = 0;
|
|
vx[p0] = isbx * nbeta / nSubPixels;
|
|
} else {
|
|
vy[p0] = rand() % (nbeta / 2);
|
|
vx[p0] = rand() % (nbeta / 2);
|
|
if (nSubPixels % 2 == 0 && isbx == nSubPixels / 2)
|
|
vx[p0] = nbeta / 2;
|
|
if (nSubPixels % 2 == 0 && isby == nSubPixels / 2)
|
|
vy[p0] = nbeta / 2;
|
|
}
|
|
// cout << "(" << vx[p0] << " , " << vy[p0] << " ) \t" ;
|
|
// }
|
|
}
|
|
// cout << endl;
|
|
}
|
|
// cout << "rand" << endl;
|
|
|
|
while (bad) {
|
|
|
|
for (isby = 0; isby < (nSubPixels + 1) / 2 + 1; isby++) {
|
|
|
|
for (isbx = 0; isbx < (nSubPixels + 1) / 2 + 1; isbx++) {
|
|
arrx[isbx] = vx[isby * (nSubPixels + 1) + isbx];
|
|
arry[isbx] = vy[isbx * (nSubPixels + 1) + isby];
|
|
// cout << isbx << " " << arrx[isbx] << " " << isby << " "
|
|
// << arry[isbx] << endl;
|
|
}
|
|
|
|
sort(arrx, arrx + (nSubPixels + 1) / 2 + 1);
|
|
sort(arry, arry + (nSubPixels + 1) / 2 + 1);
|
|
|
|
// cout << "*****"<< endl;
|
|
// cout << endl;
|
|
|
|
for (int isbx = 0; isbx < (nSubPixels + 1) / 2 + 1; isbx++) {
|
|
vx[isby * (nSubPixels + 1) + isbx] = arrx[isbx];
|
|
vy[isbx * (nSubPixels + 1) + isby] = arry[isbx];
|
|
|
|
vx[(nSubPixels - isby) * (nSubPixels + 1) +
|
|
(nSubPixels - isbx)] = nbeta - arrx[isbx];
|
|
vy[(nSubPixels - isbx) * (nSubPixels + 1) +
|
|
(nSubPixels - isby)] = nbeta - arry[isbx];
|
|
|
|
vx[isby * (nSubPixels + 1) + (nSubPixels - isbx)] =
|
|
nbeta - arrx[isbx];
|
|
vy[isbx * (nSubPixels + 1) + (nSubPixels - isby)] =
|
|
arry[isbx];
|
|
|
|
vx[(nSubPixels - isby) * (nSubPixels + 1) + (isbx)] =
|
|
arrx[isbx];
|
|
vy[(nSubPixels - isbx) * (nSubPixels + 1) + (isby)] =
|
|
nbeta - arry[isbx];
|
|
}
|
|
}
|
|
|
|
/* for (isby=0; isby<nSubPixels+1; isby++) { */
|
|
|
|
/* for (isbx=0; isbx<nSubPixels+1; isbx++) { */
|
|
|
|
/* cout << "("<< vx[isby*(nSubPixels+1)+isbx] << " " <<
|
|
* vy[isby*(nSubPixels+1)+isbx] << ")\t";//<< endl; */
|
|
/* } */
|
|
/* cout << endl; */
|
|
/* } */
|
|
|
|
bad = 0;
|
|
for (isby = 1; isby < (nSubPixels + 1) / 2 + 1; isby++) {
|
|
|
|
for (isbx = 1; isbx < (nSubPixels + 1) / 2 + 1; isbx++) {
|
|
|
|
if (heta[vx[isby * (nSubPixels + 1) + isbx] +
|
|
vy[isby * (nSubPixels + 1) + isbx] * nbeta] <
|
|
avg * (nSubPixels * nSubPixels) / (nbeta * nbeta)) {
|
|
// cout << ii << " " << isbx << " " << isby << " "
|
|
//<< vx[isby*(nSubPixels+1)+isbx] << " " <<
|
|
//vy[isby*(nSubPixels+1)+isbx] << " " <<
|
|
//heta[vx[isby*(nSubPixels+1)+isbx]+vy[isby*(nSubPixels+1)+isbx]*nbeta]
|
|
//<< endl;
|
|
if (nSubPixels % 2 == 0 && isbx == nSubPixels / 2)
|
|
;
|
|
else
|
|
vx[isby * (nSubPixels + 1) + isbx] =
|
|
rand() % (nbeta / 2);
|
|
|
|
if (nSubPixels % 2 == 0 && isbx == nSubPixels / 2)
|
|
;
|
|
else
|
|
vy[isby * (nSubPixels + 1) + isbx] =
|
|
rand() % (nbeta / 2);
|
|
|
|
if (bad == 0)
|
|
ii++;
|
|
|
|
bad = 1;
|
|
// break;
|
|
}
|
|
}
|
|
// if (bad) break;
|
|
}
|
|
// cout << "sort" << endl;
|
|
}
|
|
|
|
cout << ii << " sub iteractions "
|
|
<< avg * (nSubPixels * nSubPixels) / (nbeta * nbeta) << endl;
|
|
|
|
double m, q;
|
|
int in_quad;
|
|
int p[4];
|
|
int p1x, p2x, p1y, p2y;
|
|
// cout << nbeta << endl;
|
|
double angle;
|
|
double dtheta, theta1, theta2;
|
|
|
|
for (int ibx = 0; ibx < nbeta; ibx++) {
|
|
|
|
for (int iby = 0; iby < nbeta; iby++) {
|
|
|
|
in_quad = 0;
|
|
|
|
/* if (ibx==0) */
|
|
/* isbx=0; */
|
|
/* else */
|
|
/* isbx= (newhhx[ibx-1+iby*nbeta])/bsize-1; */
|
|
/* if (isbx<0) isbx=0; */
|
|
/* if (isbx>nSubPixels-1) isbx=nSubPixels-1; */
|
|
|
|
/* if (iby==0) */
|
|
/* isby=0; */
|
|
/* else */
|
|
/* isby= (newhhx[ibx+(iby-1)*nbeta])/bsize-1; */
|
|
|
|
/* if (isby<0) isbx=0; */
|
|
/* if (isby>nSubPixels-1) isby=nSubPixels-1; */
|
|
/* // cout << isbx << " " << isby << endl; */
|
|
|
|
for (isby = 0; isby < nSubPixels; isby++) {
|
|
|
|
for (isbx = 0; isbx < nSubPixels; isbx++) {
|
|
|
|
// cout << ibx << " " << iby << " " << isbx << " " <<
|
|
// isby << endl;
|
|
p[0] = isby * (nSubPixels + 1) + isbx;
|
|
p[1] = isby * (nSubPixels + 1) + isbx + 1;
|
|
p[2] = (isby + 1) * (nSubPixels + 1) + isbx + 1;
|
|
p[3] = (isby + 1) * (nSubPixels + 1) + isbx;
|
|
|
|
angle = 0;
|
|
for (int i = 0; i < 4; i++) {
|
|
p1x = vx[p[i]] - ibx;
|
|
p1y = vy[p[i]] - iby;
|
|
p2x = vx[p[(i + 1) % 4]] - ibx;
|
|
p2y = vy[p[(i + 1) % 4]] - iby;
|
|
theta1 = atan2(p1y, p1x);
|
|
theta2 = atan2(p2y, p2x);
|
|
dtheta = theta2 - theta1;
|
|
|
|
while (dtheta > PI)
|
|
dtheta -= TWOPI;
|
|
while (dtheta < -PI)
|
|
dtheta += TWOPI;
|
|
|
|
angle += dtheta;
|
|
}
|
|
|
|
if (abs((double)angle) < PI)
|
|
in_quad = 0;
|
|
else
|
|
in_quad = 1;
|
|
|
|
if (in_quad) {
|
|
newhhx[ibx + iby * nbeta] = bsize * ((double)isbx);
|
|
newhhy[ibx + iby * nbeta] = bsize * ((double)isby);
|
|
break;
|
|
}
|
|
}
|
|
if (in_quad)
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
// cout << "hist" << endl;
|
|
return calcDiff(avg, newhhx, newhhy);
|
|
}
|
|
|
|
public:
|
|
etaInterpolationRandomBins(int nx = 400, int ny = 400, int ns = 25,
|
|
int nb = -1, double emin = 1, double emax = 0)
|
|
: etaInterpolationPosXY(nx, ny, ns, nb, emin, emax){};
|
|
|
|
etaInterpolationRandomBins(etaInterpolationRandomBins *orig)
|
|
: etaInterpolationPosXY(orig){};
|
|
|
|
virtual etaInterpolationRandomBins *Clone() {
|
|
|
|
return new etaInterpolationRandomBins(this);
|
|
};
|
|
|
|
virtual void prepareInterpolation(int &ok) {
|
|
ok = 1;
|
|
cout << "Adaptive bins" << endl;
|
|
|
|
///*Eta Distribution Rebinning*///
|
|
double bsize = 1. / nSubPixels; // precision
|
|
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels
|
|
// = "<<nSubPixels<<endl;
|
|
double tot_eta = 0;
|
|
for (int ip = 0; ip < nbeta * nbeta; ip++)
|
|
tot_eta += heta[ip];
|
|
if (tot_eta <= 0) {
|
|
ok = 0;
|
|
return;
|
|
};
|
|
|
|
int ii = 0;
|
|
|
|
int nint = 1000;
|
|
|
|
double thr = 1. / ((double)nSubPixels);
|
|
double avg = tot_eta / ((double)(nSubPixels * nSubPixels));
|
|
cout << "total eta entries is :" << tot_eta << " avg: " << avg << endl;
|
|
cout << "Start " << endl;
|
|
double old_diff = -1, new_diff = -1;
|
|
// cout << " diff= " << new_diff << endl;
|
|
|
|
etaInterpolationPosXY::prepareInterpolation(ok);
|
|
|
|
old_diff = calcDiff(avg, hhx, hhy);
|
|
cout << " diff= " << old_diff << endl;
|
|
|
|
int iint = 0;
|
|
float *newhhx = new float[nbeta * nbeta]; // profile x
|
|
float *newhhy = new float[nbeta * nbeta]; // profile y
|
|
int igood = 0, ibad = 0;
|
|
#ifdef SAVE_ALL
|
|
int etabins = nbeta;
|
|
float *etah = new float[nbeta * nbeta];
|
|
char tit[1000];
|
|
#endif
|
|
|
|
while (iint < nint) {
|
|
|
|
cout << "Iteration " << iint << endl;
|
|
new_diff = iterate(newhhx, newhhy, avg);
|
|
// new_diff=calcDiff(avg, newhhx, newhhy);
|
|
cout << " diff= " << new_diff << " ( " << old_diff << " ) " << endl;
|
|
|
|
/* #ifdef SAVE_ALL */
|
|
/* for (int ii=0; ii<etabins*etabins; ii++) { */
|
|
/* etah[ii]=newhhx[ii]; */
|
|
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii <<
|
|
* etah[ii] << endl; */
|
|
/* } */
|
|
/* sprintf(tit,"/scratch/randeta_hhx_%d.tiff",iint); */
|
|
/* WriteToTiff(etah, tit, etabins, etabins); */
|
|
|
|
/* for (int ii=0; ii<etabins*etabins; ii++) { */
|
|
/* etah[ii]=newhhy[ii]; */
|
|
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii <<
|
|
* etah[ii] << endl; */
|
|
/* } */
|
|
/* sprintf(tit,"/scratch/randeta_hhy_%d.tiff",iint); */
|
|
/* WriteToTiff(etah, tit, etabins, etabins); */
|
|
/* #endif */
|
|
|
|
if (new_diff < old_diff) {
|
|
|
|
cout << "******************** GOOD! ***********************"
|
|
<< endl;
|
|
delete[] hhx;
|
|
delete[] hhy;
|
|
igood++;
|
|
hhx = newhhx;
|
|
hhy = newhhy;
|
|
newhhx = new float[nbeta * nbeta]; // profile x */
|
|
newhhy = new float[nbeta * nbeta]; // profile y */
|
|
old_diff = new_diff;
|
|
} else
|
|
ibad++;
|
|
|
|
iint++;
|
|
}
|
|
delete[] newhhx;
|
|
delete[] newhhy;
|
|
|
|
cout << "performed " << iint << " iterations of which " << igood
|
|
<< " positive " << endl;
|
|
|
|
/* #ifdef SAVE_ALL */
|
|
|
|
/* for (int ii=0; ii<etabins*etabins; ii++) { */
|
|
/* etah[ii]=hhx[ii]; */
|
|
/* } */
|
|
/* sprintf(tit,"/scratch/eta_hhx_%d.tiff",id); */
|
|
/* WriteToTiff(etah, tit, etabins, etabins); */
|
|
|
|
/* for (int ii=0; ii<etabins*etabins; ii++) { */
|
|
/* etah[ii]=hhy[ii]; */
|
|
/* } */
|
|
/* sprintf(tit,"/scratch/eta_hhy_%d.tiff",id); */
|
|
/* WriteToTiff(etah, tit, etabins, etabins); */
|
|
|
|
/* for (int ii=0; ii<etabins*etabins; ii++) { */
|
|
/* etah[ii]=heta[ii]; */
|
|
/* } */
|
|
/* sprintf(tit,"/scratch/eta_%d.tiff",id); */
|
|
/* WriteToTiff(etah, tit, etabins, etabins); */
|
|
/* delete [] etah; */
|
|
/* #endif */
|
|
return;
|
|
}
|
|
};
|
|
|
|
#endif
|