Skip to content

GitLab

Commit 5a1927a3 authored by kraus's avatar kraus
Browse files

- detect number of different particle types automaticly 

- new laser profile class with automatic detection of center of mass of profile and possiblity to flip and rotate
parent af1a4084
Hide whitespace changes
Inline Side-by-side
Showing with 384 additions and 314 deletions
classic/5.0/src/Algorithms/PartBunch.cpp
View file @ 5a1927a3
......@@ -622,10 +622,10 @@ void PartBunch::calcGammas() {
reduce(pInBin, pInBin, OpAddAssign());
if(pInBin != 0) {
bingamma_m[i] /= pInBin;
INFOMSG(level2 << "Bin " << i << " gamma = " << setw(8) << scientific << setprecision(5) << bingamma_m[i] << "; NpInBin= " << setw(8) << setfill(' ') << pInBin << endl);
INFOMSG(level2 << "Bin " << std::setw(3) << i << " gamma = " << setw(8) << scientific << setprecision(5) << bingamma_m[i] << "; NpInBin= " << setw(8) << setfill(' ') << pInBin << endl);
} else {
bingamma_m[i] = 1.0;
INFOMSG(level2 << "Bin " << i << " has no particles " << endl);
INFOMSG(level2 << "Bin " << std::setw(3) << i << " has no particles " << endl);
}
s += pInBin;
}
......
src/Algorithms/ParallelTTracker.cpp
View file @ 5a1927a3
......@@ -2094,7 +2094,7 @@ void ParallelTTracker::computeExternalFields() {
}
if(ne > 0)
msg << level2 << "* Deleted " << ne << " particles, "
msg << level1 << "* Deleted " << ne << " particles, "
<< "remaining " << itsBunch->getTotalNum() << " particles" << endl;
}
......@@ -2847,4 +2847,4 @@ void ParallelTTracker::freeDeviceMemory() {
// mode:c
// c-basic-offset: 4
// indent-tabs-mode:nil
// End:
// End:
\ No newline at end of file
src/Distribution/Distribution.cpp
View file @ 5a1927a3
......@@ -134,6 +134,11 @@ namespace AttributesT
LASERPROFFN,
IMAGENAME,
INTENSITYCUT,
FLIPX,
FLIPY,
ROTATE90,
ROTATE180,
ROTATE270,
NPDARKCUR,
INWARDMARGIN,
EINITHR,
......@@ -847,8 +852,8 @@ Inform &Distribution::printInfo(Inform &os) const {
if (emitting_m) {
os << "* Distribution is emitted. " << endl;
os << "* Emission time = " << tEmission_m << " [sec]" << endl;
os << "* Time per bin = " << tEmission_m/numberOfEnergyBins_m << " [sec]" << endl;
os << "* Bin delta t = " << tBin_m << " [sec]" << endl;
os << "* Time per bin = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
os << "* Delta t during emission = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
os << "* " << endl;
PrintEmissionModel(os);
os << "* " << endl;
......@@ -1553,8 +1558,8 @@ void Distribution::CreateOpalE(Beam *beam,
SetDistParametersGauss(beam->getMass());
break;
case DistrTypeT::GUNGAUSSFLATTOPTH:
INFOMSG("GUNGAUSSFLATTOPTH"<<endl)
SetDistParametersFlattop(beam->getMass());
INFOMSG("GUNGAUSSFLATTOPTH"<<endl);
SetDistParametersFlattop(beam->getMass());
beamEnergy = Attributes::getReal(itsAttr[AttributesT::EKIN]);
break;
default:
......@@ -1592,9 +1597,7 @@ void Distribution::CreateOpalT(PartBunch &beam,
std::vector<Distribution *> addedDistributions,
size_t &numberOfParticles,
bool scan) {
// This is PC from BEAM
avrgpz_m = OpalData::getInstance()->getP0(); // beam.getP()/beam.getM();
addedDistributions_m = addedDistributions;
CreateOpalT(beam, numberOfParticles, scan);
}
......@@ -1871,7 +1874,7 @@ size_t Distribution::EmitParticles(PartBunch &beam, double eZ) {
currentSampleBin_m++;
if (currentSampleBin_m == numberOfSampleBins_m) {
INFOMSG("*Bin number: "
INFOMSG(level3 << "* Bin number: "
<< currentEnergyBin_m
<< " has emitted all particles (new emit)." << endl);
currentSampleBin_m = 0;
......@@ -2286,7 +2289,7 @@ void Distribution::GenerateBinomial(size_t numberOfParticles) {
yDist_m.push_back(x[1]);
pyDist_m.push_back(p[1]);
tOrZDist_m.push_back(x[2]);
pzDist_m.push_back(avrgpz_m *(1+p[2]));
pzDist_m.push_back(avrgpz_m * (1 + p[2]));
}
}
......@@ -2302,7 +2305,7 @@ void Distribution::GenerateFlattopLaserProfile(size_t numberOfParticles) {
double y = 0.0;
double py = 0.0;
laserProfile_m->GetXY(&x, &y);
laserProfile_m->getXY(x, y);
// Save to each processor in turn.
saveProcessor++;
......@@ -2310,9 +2313,9 @@ void Distribution::GenerateFlattopLaserProfile(size_t numberOfParticles) {
saveProcessor = 0;
if (Ippl::myNode() == saveProcessor) {
xDist_m.push_back(x);
xDist_m.push_back(x * sigmaR_m[0]);
pxDist_m.push_back(px);
yDist_m.push_back(y);
yDist_m.push_back(y * sigmaR_m[1]);
pyDist_m.push_back(py);
}
}
......@@ -2580,12 +2583,12 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
throw OpalException("Distribution::GenerateGaussZ",
"oops... something wrong with GSL matvec");
}
xDist_m.push_back( sigmaR_m[0]*gsl_vector_get(ry, 0));
pxDist_m.push_back( sigmaP_m[0]*gsl_vector_get(ry, 1));
yDist_m.push_back( sigmaR_m[1]*gsl_vector_get(ry, 2));
pyDist_m.push_back( sigmaP_m[1]*gsl_vector_get(ry, 3));
tOrZDist_m.push_back( sigmaR_m[2]*gsl_vector_get(ry, 4));
pzDist_m.push_back(avrgpz_m +(sigmaP_m[2]*gsl_vector_get(ry, 5)));
xDist_m.push_back( sigmaR_m[0] * gsl_vector_get(ry, 0));
pxDist_m.push_back( sigmaP_m[0] * gsl_vector_get(ry, 1));
yDist_m.push_back( sigmaR_m[1] * gsl_vector_get(ry, 2));
pyDist_m.push_back( sigmaP_m[1] * gsl_vector_get(ry, 3));
tOrZDist_m.push_back( sigmaR_m[2] * gsl_vector_get(ry, 4));
pzDist_m.push_back(avrgpz_m + sigmaP_m[2] * gsl_vector_get(ry, 5));
}
}
......@@ -3656,6 +3659,16 @@ void Distribution::SetAttributes() {
= Attributes::makeReal("INTENSITYCUT", "For background subtraction of laser "
"image profile, in % of max intensity.",
0.0);
itsAttr[AttributesT::FLIPX]
= Attributes::makeBool("FLIPX", "Flip laser profile horizontally", false);
itsAttr[AttributesT::FLIPY]
= Attributes::makeBool("FLIPY", "Flip laser profile vertically", false);
itsAttr[AttributesT::ROTATE90]
= Attributes::makeBool("ROTATE90", "Rotate laser profile 90 degrees counter clockwise", false);
itsAttr[AttributesT::ROTATE180]
= Attributes::makeBool("ROTATE180", "Rotate laser profile 180 degrees counter clockwise", false);
itsAttr[AttributesT::ROTATE270]
= Attributes::makeBool("ROTATE270", "Rotate laser profile 270 degrees counter clockwise", false);
// Dark current and field emission model parameters.
itsAttr[AttributesT::NPDARKCUR]
......@@ -3878,7 +3891,6 @@ void Distribution::SetEmissionTime(double &maxT, double &minT) {
case DistrTypeT::GUNGAUSSFLATTOPTH:
tEmission_m = tPulseLengthFWHM_m + (cutoffR_m[2] - sqrt(2.0 * log(2.0)))
* (sigmaTRise_m + sigmaTFall_m);
INFOMSG("tEm= " << tEmission_m << endl);
break;
default:
......@@ -4068,9 +4080,17 @@ void Distribution::SetDistParametersFlattop(double massIneV) {
if (!(laserProfileFileName_m == std::string(""))) {
laserImageName_m = Attributes::getString(itsAttr[AttributesT::IMAGENAME]);
laserIntensityCut_m = std::abs(Attributes::getReal(itsAttr[AttributesT::INTENSITYCUT]));
short flags = 0;
if (Attributes::getBool(itsAttr[AttributesT::FLIPX])) flags |= LaserProfile::FLIPX;
if (Attributes::getBool(itsAttr[AttributesT::FLIPY])) flags |= LaserProfile::FLIPY;
if (Attributes::getBool(itsAttr[AttributesT::ROTATE90])) flags |= LaserProfile::ROTATE90;
if (Attributes::getBool(itsAttr[AttributesT::ROTATE180])) flags |= LaserProfile::ROTATE180;
if (Attributes::getBool(itsAttr[AttributesT::ROTATE270])) flags |= LaserProfile::ROTATE270;
laserProfile_m = new LaserProfile(laserProfileFileName_m,
laserImageName_m,
laserIntensityCut_m);
laserIntensityCut_m,
flags);
}
// Legacy for ASTRAFLATTOPTH.
......
src/Distribution/LaserProfile.cpp
View file @ 5a1927a3
......@@ -9,319 +9,348 @@
// Class: LaserProfile
//
// ------------------------------------------------------------------------
#include "Distribution/LaserProfile.h"
#include "Utility/Inform.h"
#include "Utilities/OpalException.h"
using namespace std;
void LaserProfile::ReadFile(string fn, string imagestr, double cut) {
Inform m("LaserProfile::ReadFile ");
hid_t h5 = H5Fopen(fn.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
hid_t gr = H5Gopen2 (h5, imagestr.c_str(), H5P_DEFAULT);
#include <boost/filesystem.hpp>
hid_t dsetBit = H5Dopen2 (gr, "CameraBits", H5P_DEFAULT);
dsetBit = H5Dopen2 (gr, "CameraGain", H5P_DEFAULT);
#include <fstream>
#include <iostream>
#include <cmath>
double bitData[] = { -1 };
H5Dread(dsetBit, H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, bitData);
int bits = static_cast<int>(bitData[0]);
#define TESTLASEREMISSION
extern Inform *gmsg;
hsize_t dim[2];
hid_t dset = H5Dopen2 (gr, "CameraData", H5P_DEFAULT);
LaserProfile::LaserProfile(const std::string &fileName,
const std::string &imageName,
double intensityCut,
short flags):
sizeX_m(0),
sizeY_m(0),
hist2d_m(NULL),
rng_m(NULL),
pdf_m(NULL),
centerMass_m(0.0),
standardDeviation_m(0.0){
hid_t filespace = H5Dget_space(dset);
H5Sget_simple_extent_dims(filespace, dim, NULL);
filespace = H5Screate_simple(2, dim, NULL);
/*
filespace = H5Screate_simple(3, dim, NULL);
*/
unsigned short *image = readFile(fileName, imageName, intensityCut);
sizeX_m = dim[0];
sizeY_m = dim[1];
/*
int nImage = dim[0];
sizeX_m = dim[1];
sizeY_m = dim[2];
*/
// selection of hyperslab in hdf5 file:
hsize_t start[] = {0, 0}; // Start of hyperslab
hsize_t count[] = {sizeX_m, sizeY_m}; // Block count
/*
hsize_t start[] = {imageNumber,0,0}; // Start of hyperslab
hsize_t count[] = {imageNumber,sizeX_m,sizeY_m}; // Block count
*/
H5Sselect_hyperslab(filespace, H5S_SELECT_SET, start, NULL, count, NULL);
unsigned short int *image = new unsigned short int [sizeX_m*sizeY_m];
hid_t mem = H5Screate_simple(2, count, NULL);
/*
hid_t mem = H5Screate_simple(3, count, NULL);
*/
H5Dread(dset, H5T_NATIVE_USHORT, mem, filespace, H5P_DEFAULT, image);
m << "# Read image done .... sizeX= " << sizeX_m << " sizeY= " << sizeY_m << " Bits= " << bits << endl;
int pixel = 0;
hsize_t col = 1;
hsize_t row = 1;
hist2d_m = gsl_histogram2d_alloc(sizeX_m, sizeY_m);
gsl_histogram2d_set_ranges_uniform(hist2d_m, 0.0, 1.0, 0.0, 1.0); // all bins are set to 0 too
cout << "sizeX_m=" << sizeX_m << "; sizeY_m=" << sizeY_m << endl;
bool doBackGroundCut_m = 1;
unsigned short int profileMax_m;
if(doBackGroundCut_m)
GetProfileMax(&profileMax_m, image);
// exit ();
/*
for (col=1; col<=sizeX_m; col++) {
for (row=1; row<=sizeY_m; row++) {
// int val = (image[pixel++] >> shift); // lowest bits are insignificant
int val = (image[pixel++]);
// if ((200 < col < 300) && (200 < row < 300))
// gsl_histogram2d_accumulate (hist2d_m, col/sizeX_m, row/sizeY_m, val);
if ( (700<col)&&(col<1000)&&(1<row)&&(row<300) ) {
gsl_histogram2d_accumulate (hist2d_m, double(col)/sizeX_m, double(row)/sizeY_m, val);
if( ((col==5)||(col==250)||(col==500)||(col==800))&&((row==100)||(row==300)||(row==500)||(row==700)) ) {
cout << col << " " << row << " val=" << val << endl;
cout << "col/sizeX_m=" << double(col)/sizeX_m << "; row/sizeY_m=" << double(row)/sizeY_m << endl;
}
}
if (flags & FLIPX) flipX(image);
if (flags & FLIPY) flipY(image);
if (flags & ROTATE90) {
swapXY(image);
flipX(image);
}
if (flags & ROTATE180) {
flipX(image);
flipY(image);
}
cout << "xmax=" << gsl_histogram2d_xmax(hist2d_m) << "; xmin=" << gsl_histogram2d_xmin(hist2d_m) << endl;
cout << "ymax=" << gsl_histogram2d_ymax(hist2d_m) << "; ymin=" << gsl_histogram2d_ymin(hist2d_m) << endl;
cout << "max_val=" << gsl_histogram2d_max_val(hist2d_m) << "; min_val=" << gsl_histogram2d_min_val(hist2d_m) << endl;
cout << "xmean=" << gsl_histogram2d_xmean(hist2d_m) << "; xsigma=" << gsl_histogram2d_xsigma(hist2d_m) << endl;
cout << "ymean=" << gsl_histogram2d_ymean(hist2d_m) << "; ysigma=" << gsl_histogram2d_ysigma(hist2d_m) << endl;
cout << "sum=" << gsl_histogram2d_sum(hist2d_m) << endl;
*/
/*
if (doBackGroundCut_m) {
GetProfileMax();
double wkjdfkjdf=0.14;
for (col=1; col<=sizeX_m; col++) {
for (row=1; row<=sizeY_m; row++) {
val =
// Background cut starts here
if ( (val*1.0)/(MAX*1.0) < wkjdfkjdf )
val = 0;
if ( (val*1.0)/(profileMax_m*1.0) >= wkjdfkjdf )
val = ((val*1.0)/profileMax_m - wkjdfkjdf)/(1-wkjdfkjdf) * profileMax_m;
// end of cut
if ( (700<col)&&(col<1000)&&(1<row)&&(row<300) ) {
gsl_histogram2d_accumulate (hist2d_m, double(col)/sizeX_m, double(row)/sizeY_m, val);
if( ((col==5)||(col==250)||(col==500)||(col==800))&&((row==100)||(row==300)||(row==500)||(row==700)) ) {
cout << col << " " << row << " val=" << val << endl;
cout << "col/sizeX_m=" << double(col)/sizeX_m << "; row/sizeY_m=" << double(row)/sizeY_m << endl;
}
}
}
}
cout << "profileMax = " << profileMax << endl;
if (flags & ROTATE270) {
swapXY(image);
flipY(image);
}
*/
/*
else {
*/
double insCut = 0.45;
// double MAX = gsl_histogram2d_max_val(hist2d_m); // MAX =Maximum im Laserprofil
// cout << "profileMax_m = " << profileMax_m << endl;
for(col = sizeX_m; col >= 1; col--) {
for(row = 1; row <= sizeY_m; row++) {
//
// if (pixel == 0)
// cout << "pixel == 0; image[pixel] = " << image[pixel] << "; image[pixel+1] = " << image[pixel+1] << endl;
// int val = (image[pixel++] >> shift); // lowest bits are insignificant
double val;
/*
if (pixel == 0)
{ val = (image[pixel++]);
cout << "val = " << val << "; image[pixel] = " << image[pixel] << endl; }
else
val = (image[pixel++]);
*/
val = (image[pixel++]);
// Background cut starts here
if((val * 1.0) / (profileMax_m * 1.0) < insCut)
val = 0;
if((val * 1.0) / (profileMax_m * 1.0) >= insCut)
val = ((val * 1.0) / profileMax_m - insCut) / (1 - insCut) * profileMax_m;
// end of cut
if((700 < col) && (col < 1000) && (1 < row) && (row < 300)) {
gsl_histogram2d_accumulate(hist2d_m, double(col) / sizeX_m, double(row) / sizeY_m, val);
// if( ((col==5)||(col==250)||(col==500)||(col==800))&&((row==100)||(row==300)||(row==500)||(row==700)) ) {
// cout << col << " " << row << " val=" << val << endl;
// cout << "col/sizeX_m=" << double(col)/sizeX_m << "; row/sizeY_m=" << double(row)/sizeY_m << endl;
// }
}
#ifdef TESTLASEREMISSION
saveOrigData(image);
#endif
filterSpikes(image);
normalizeProfileData(intensityCut, image);
computeProfileStatistics(image);
fillHistrogram(image);
delete[] image;
// if ((200 < col < 300) && (200 < row < 300))
// gsl_histogram2d_accumulate (hist2d_m, col/sizeX_m, row/sizeY_m, val);
}
setupRNG();
#ifdef TESTLASEREMISSION
saveHistogram();
sampleDist();
#endif
printInfo();
}
LaserProfile::~LaserProfile() {
gsl_histogram2d_pdf_free(pdf_m);
gsl_histogram2d_free(hist2d_m);
gsl_rng_free(rng_m);
}
unsigned short * LaserProfile::readFile(const std::string &fileName,
const std::string &imageName,
double insCut) {
namespace fs = boost::filesystem;
if (!fs::exists(fileName)) {
throw OpalException("LaserProfile::readFile",
"given file '" + fileName + "' does not exist");
}
/*
hid_t h5 = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
hid_t group = H5Gopen2 (h5, imageName.c_str(), H5P_DEFAULT);
if (group < 0) {
throw OpalException("LaserProfile::readFile",
"given image name '" + imageName + "' does not exist");
}
*/
hsize_t dim[2];
hid_t dataSet = H5Dopen2 (group, "CameraData", H5P_DEFAULT);
if (dataSet < 0) {
throw OpalException("LaserProfile::readFile",
"data set 'CameraData' does not exist");
}
SaveDist();
hid_t dataSetSpace = H5Dget_space(dataSet);
H5Sget_simple_extent_dims(dataSetSpace, dim, NULL);
hid_t filespace = H5Screate_simple(2, dim, NULL);
gsl_rng_env_setup();
sizeX_m = dim[0];
sizeY_m = dim[1];
const gsl_rng_type *T = gsl_rng_default;
r_m = gsl_rng_alloc(T);
hsize_t startHyperslab[] = {0, 0};
hsize_t blockCount[] = {sizeX_m, sizeY_m};
pdf_m = gsl_histogram2d_pdf_alloc(sizeX_m, sizeY_m);
gsl_histogram2d_pdf_init(pdf_m, hist2d_m);
H5Sselect_hyperslab(filespace, H5S_SELECT_SET, startHyperslab, NULL, blockCount, NULL);
unsigned short *image = new unsigned short [sizeX_m * sizeY_m];
hid_t mem = H5Screate_simple(2, blockCount, NULL);
SampleDist();
H5Dread(dataSet, H5T_NATIVE_USHORT, mem, filespace, H5P_DEFAULT, image);
double aa, bb;
GetXY(&aa, &bb);
cout << "x aa=" << aa << "; y bb=" << bb << endl;
H5Sclose(mem);
H5Sclose(filespace);
H5Dclose(dataSetSpace);
H5Dclose(dataSet);
H5Gclose(group);
H5Fclose(h5);
return image;
}
void LaserProfile::flipX(unsigned short *image) {
unsigned int col2 = sizeX_m - 1;
for (unsigned int col1 = 0; col1 < col2; ++ col1, -- col2) {
unsigned int pixel1 = col1 * sizeY_m;
unsigned int pixel2 = col2 * sizeY_m;
for (unsigned int row = 0; row < sizeY_m; ++ row) {
std::swap(image[pixel1 ++], image[pixel2 ++]);
}
}
}
/*-----*/
void LaserProfile::SaveDist() {
FILE *fh = fopen("dist2dhist.dat", "w");
gsl_histogram2d_fprintf(fh, hist2d_m, "%g", "%g");
fclose(fh);
void LaserProfile::flipY(unsigned short *image) {
for (unsigned int col = 0; col < sizeX_m; ++ col) {
unsigned int pixel1 = col * sizeY_m;
unsigned int pixel2 = (col + 1) * sizeY_m - 1;
unsigned int row2 = sizeY_m - 1;
for (unsigned int row1 = 0; row1 < row2; ++ row1, -- row2) {
std::swap(image[pixel1 ++], image[pixel2 --]);
}
}
}
void LaserProfile::swapXY(unsigned short *image) {
unsigned short *copy = new unsigned short [sizeX_m * sizeY_m];
for (unsigned int col = 0; col < sizeX_m; ++ col) {
for (unsigned int row = 0; row < sizeY_m; ++ row) {
unsigned int pixel1 = col * sizeY_m + row;
unsigned int pixel2 = row * sizeX_m + col;
void LaserProfile::SampleDist() {
FILE *fh = fopen("distsampled.dat", "w");
double x, y;
copy[pixel2] = image[pixel1];
}
}
for(int i = 0; i < 10000; i++) {
GetXY(&x, &y);
fprintf(fh, "%g %g \n", x, y);
for (unsigned int pixel = 0; pixel < sizeX_m * sizeY_m; ++ pixel) {
image[pixel] = copy[pixel];
}
delete[] copy;
fclose(fh);
std::swap(sizeX_m, sizeY_m);
}
void LaserProfile::GetXY(double *s_x, double *s_y) {
double x, y;
double u = gsl_rng_uniform(r_m);
double v = gsl_rng_uniform(r_m);
gsl_histogram2d_pdf_sample(pdf_m, u, v, &x, &y);
*s_x = x;
*s_y = y;
void LaserProfile::filterSpikes(unsigned short *image) {
hsize_t idx = sizeY_m + 1;
hsize_t idxN = idx - sizeY_m;
hsize_t idxNW = idxN - 1, idxNE = idxN + 1;
hsize_t idxW = idx - 1, idxE = idx + 1;
hsize_t idxS = idx + sizeY_m;
hsize_t idxSW = idxS - 1, idxSE = idxS + 1;
for (hsize_t i = 1; i < sizeX_m - 1; ++ i) {
for (hsize_t j = 1; j < sizeY_m - 1; ++ j) {
if (image[idx] > 0) {
if (image[idxNW] == 0 &&
image[idxN] == 0 &&
image[idxNE] == 0 &&
image[idxW] == 0 &&
image[idxE] == 0 &&
image[idxSW] == 0 &&
image[idxS] == 0 &&
image[idxSE] == 0) {
image[idx] = 0;
}
}
// *s_x = x*sizeX_m;
// *s_y = y*sizeY_m;
++ idxNW; ++ idxN; ++ idxNE;
++ idxW; ++ idx; ++ idxE;
++ idxSW; ++ idxS; ++ idxSE;
}
}
}
void LaserProfile::GetProfileMax(unsigned short int *profileMax_m, unsigned short int image[]) {
void LaserProfile::normalizeProfileData(double intensityCut, unsigned short *image) {
unsigned short int profileMax;
getProfileMax(profileMax, image);
int image_len = sizeX_m * sizeY_m;
//*profileMax_m = image[i]; // start with max = first element
unsigned int pixel = 0;
for (unsigned int col = 0; col < sizeX_m; ++ col) {
for(unsigned int row = 0; row < sizeY_m; ++ row, ++ pixel) {
double val = (double(image[pixel]) / profileMax - intensityCut) / (1.0 - intensityCut);
*profileMax_m = 0;
val = std::max(0.0, val);
image[pixel] = std::floor(val * profileMax + 0.5);
}
}
}
cout << " profileMax_m = " << profileMax_m << " ; *profileMax_m = " << *profileMax_m << endl;
for(int i = 0; i < image_len; i++) {
if(image[i] > *profileMax_m)
*profileMax_m = image[i];
void LaserProfile::computeProfileStatistics(unsigned short *image) {
double totalMass = 0.0;
centerMass_m = Vector_t(0.0);
standardDeviation_m = Vector_t(0.0);