LaserProfile.cpp 10.7 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11
// ------------------------------------------------------------------------
// $RCSfile: LaserProfile.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.3.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: LaserProfile
//
// ------------------------------------------------------------------------
12

gsell's avatar
gsell committed
13 14
#include "Distribution/LaserProfile.h"
#include "Utility/Inform.h"
15
#include "Utilities/OpalException.h"
gsell's avatar
gsell committed
16

17
#include <boost/filesystem.hpp>
gsell's avatar
gsell committed
18

19 20 21
#include <fstream>
#include <iostream>
#include <cmath>
gsell's avatar
gsell committed
22

23
#define TESTLASEREMISSION
gsell's avatar
gsell committed
24

25
extern Inform *gmsg;
gsell's avatar
gsell committed
26

27 28 29 30 31 32 33 34 35 36 37
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){
gsell's avatar
gsell committed
38

39
    unsigned short *image = readFile(fileName, imageName, intensityCut);
gsell's avatar
gsell committed
40

41 42 43 44 45
    if (flags & FLIPX) flipX(image);
    if (flags & FLIPY) flipY(image);
    if (flags & ROTATE90) {
        swapXY(image);
        flipX(image);
gsell's avatar
gsell committed
46
    }
47 48 49
    if (flags & ROTATE180) {
        flipX(image);
        flipY(image);
gsell's avatar
gsell committed
50
    }
51 52 53
    if (flags & ROTATE270) {
        swapXY(image);
        flipY(image);
gsell's avatar
gsell committed
54 55
    }

56 57 58
#ifdef TESTLASEREMISSION
    saveOrigData(image);
#endif
gsell's avatar
gsell committed
59

60 61 62 63
    filterSpikes(image);
    normalizeProfileData(intensityCut, image);
    computeProfileStatistics(image);
    fillHistrogram(image);
gsell's avatar
gsell committed
64

65
    delete[] image;
gsell's avatar
gsell committed
66

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
    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");
gsell's avatar
gsell committed
91 92
    }

93 94
    hid_t h5 = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
    hid_t group = H5Gopen2 (h5, imageName.c_str(), H5P_DEFAULT);
gsell's avatar
gsell committed
95

96 97 98
    if (group < 0) {
        throw OpalException("LaserProfile::readFile",
                            "given image name '" + imageName + "' does not exist");
gsell's avatar
gsell committed
99 100
    }

101 102
    hsize_t dim[2];
    hid_t   dataSet = H5Dopen2 (group, "CameraData", H5P_DEFAULT);
gsell's avatar
gsell committed
103

104 105 106 107
    if (dataSet < 0) {
        throw OpalException("LaserProfile::readFile",
                            "data set 'CameraData' does not exist");
    }
gsell's avatar
gsell committed
108

109 110 111
    hid_t dataSetSpace = H5Dget_space(dataSet);
    H5Sget_simple_extent_dims(dataSetSpace, dim, NULL);
    hid_t filespace = H5Screate_simple(2, dim, NULL);
gsell's avatar
gsell committed
112

113 114
    sizeX_m        = dim[0];
    sizeY_m        = dim[1];
gsell's avatar
gsell committed
115

116 117
    hsize_t startHyperslab[]  = {0, 0};
    hsize_t blockCount[]  = {sizeX_m, sizeY_m};
gsell's avatar
gsell committed
118

119 120 121 122
    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);
gsell's avatar
gsell committed
123

124
    H5Dread(dataSet, H5T_NATIVE_USHORT, mem, filespace, H5P_DEFAULT, image);
gsell's avatar
gsell committed
125

126 127 128 129 130 131
    H5Sclose(mem);
    H5Sclose(filespace);
    H5Dclose(dataSetSpace);
    H5Dclose(dataSet);
    H5Gclose(group);
    H5Fclose(h5);
gsell's avatar
gsell committed
132

133
    return image;
gsell's avatar
gsell committed
134 135
}

136 137 138
void LaserProfile::flipX(unsigned short *image) {
    unsigned int col2 = sizeX_m - 1;
    for (unsigned int col1 = 0; col1 < col2; ++ col1, -- col2) {
gsell's avatar
gsell committed
139

140 141 142 143 144 145 146
        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 ++]);
        }
    }
}
gsell's avatar
gsell committed
147

148 149 150 151 152 153 154 155 156 157 158
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 --]);
        }
    }
gsell's avatar
gsell committed
159 160
}

161 162 163 164 165 166
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;
gsell's avatar
gsell committed
167

168 169 170
            copy[pixel2] = image[pixel1];
        }
    }
gsell's avatar
gsell committed
171

172 173
    for (unsigned int pixel = 0; pixel < sizeX_m * sizeY_m; ++ pixel) {
        image[pixel] = copy[pixel];
gsell's avatar
gsell committed
174 175
    }

176
    delete[] copy;
gsell's avatar
gsell committed
177

178
    std::swap(sizeX_m, sizeY_m);
gsell's avatar
gsell committed
179 180
}

181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
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;
                }
            }
gsell's avatar
gsell committed
205

206 207 208 209 210
            ++ idxNW; ++ idxN; ++ idxNE;
            ++ idxW;  ++ idx;  ++ idxE;
            ++ idxSW; ++ idxS; ++ idxSE;
        }
    }
gsell's avatar
gsell committed
211 212
}

213 214 215
void LaserProfile::normalizeProfileData(double intensityCut, unsigned short *image) {
    unsigned short int profileMax;
    getProfileMax(profileMax, image);
gsell's avatar
gsell committed
216

217 218 219 220
    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);
gsell's avatar
gsell committed
221

222 223 224 225 226
            val = std::max(0.0, val);
            image[pixel] = std::floor(val * profileMax + 0.5);
        }
    }
}
gsell's avatar
gsell committed
227

228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
void LaserProfile::computeProfileStatistics(unsigned short *image) {
    double totalMass = 0.0;
    centerMass_m = Vector_t(0.0);
    standardDeviation_m = Vector_t(0.0);

    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 = image[pixel];

            centerMass_m(0) += val * col;
            centerMass_m(1) += val * row;
            standardDeviation_m(0) += val * col*col;
            standardDeviation_m(1) += val * row*row;
            totalMass += val;
        }
gsell's avatar
gsell committed
244
    }
245 246 247 248 249 250

    centerMass_m /= totalMass;
    standardDeviation_m(0) = sqrt(standardDeviation_m(0) / totalMass -
                                  centerMass_m(0)*centerMass_m(0));
    standardDeviation_m(1) = sqrt(standardDeviation_m(1) / totalMass -
                                  centerMass_m(1)*centerMass_m(1));
gsell's avatar
gsell committed
251 252
}

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
void LaserProfile::fillHistrogram(unsigned short *image) {
    unsigned int histSizeX = std::floor(10 * standardDeviation_m(0) + 0.5);
    unsigned int histSizeY = std::floor(10 * standardDeviation_m(1) + 0.5);
    histSizeX += (1 - histSizeX % 2);
    histSizeY += (1 - histSizeY % 2);
    double histRangeX = histSizeX / standardDeviation_m(0) / 2;
    double histRangeY = histSizeY / standardDeviation_m(1) / 2;
    double binSizeX = histRangeX * 2 / histSizeX;
    double binSizeY = histRangeY * 2 / histSizeY;

    hist2d_m = gsl_histogram2d_alloc(histSizeX, histSizeY);
    gsl_histogram2d_set_ranges_uniform(hist2d_m,
                                       -histRangeX, histRangeX,
                                       -histRangeY, histRangeY);

    unsigned int pixel = 0;
    for (unsigned int col = 0; col < sizeX_m; ++ col) {
        double x = (col - centerMass_m(0)) / standardDeviation_m(0);
        x = std::floor(x / binSizeX + 0.5) * binSizeX;
        if (std::abs(x) > 5.0) {
            pixel += sizeY_m;
            continue;
        }
gsell's avatar
gsell committed
276 277


278 279 280 281
        for (unsigned int row = 0; row < sizeY_m; ++ row, ++ pixel) {
            double val = image[pixel];
            double y = (row - centerMass_m(1)) / standardDeviation_m(1);
            y = std::floor(y / binSizeY + 0.5) * binSizeY;
gsell's avatar
gsell committed
282

283
            if (std::abs(y) > 5.0) continue;
gsell's avatar
gsell committed
284

285 286 287
            gsl_histogram2d_accumulate(hist2d_m, x, y, val);
        }
    }
gsell's avatar
gsell committed
288 289
}

290 291
void LaserProfile::setupRNG() {
    unsigned int histSizeX = hist2d_m->nx, histSizeY = hist2d_m->ny;
gsell's avatar
gsell committed
292

293
    gsl_rng_env_setup();
gsell's avatar
gsell committed
294

295 296
    const gsl_rng_type *T = gsl_rng_default;
    rng_m = gsl_rng_alloc(T);
gsell's avatar
gsell committed
297

298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
    pdf_m = gsl_histogram2d_pdf_alloc(histSizeX, histSizeY);
    gsl_histogram2d_pdf_init(pdf_m, hist2d_m);
}

void LaserProfile::printInfo() {
    INFOMSG(level3 << "* **********************************************************************************\n");
    INFOMSG("* LASER PROFILE \n");
    INFOMSG("* size = " << sizeX_m << " x " << sizeY_m << " pixels \n");
    INFOMSG("* center of mass: x = " << centerMass_m(0) << ", y = " << centerMass_m(1) << "\n");
    INFOMSG("* standard deviation: x = " << standardDeviation_m(0) << ", y = " << standardDeviation_m(1) << "\n");
    INFOMSG("* **********************************************************************************" << endl);
}

void LaserProfile::saveOrigData(unsigned short *image) {
    std::ofstream out("data/originalLaserProfile.dat");
gsell's avatar
gsell committed
313

314 315 316 317 318 319 320 321 322 323 324 325 326 327
    unsigned int index = 0;
    for (unsigned int i = 0; i < sizeX_m; ++ i) {
        for (unsigned int j = 0; j < sizeY_m; ++ j) {
            out << image[index++] << "\t";
        }
        out << std::endl;
    }
}

void LaserProfile::saveHistogram() {
    FILE  *fh = fopen("data/LaserHistogram.dat", "w");
    gsl_histogram2d_fprintf(fh, hist2d_m, "%g", "%g");
    fclose(fh);
}
gsell's avatar
gsell committed
328

329 330 331
void LaserProfile::sampleDist() {
    std::ofstream fh("data/LaserEmissionSampled.dat");
    double x, y;
gsell's avatar
gsell committed
332

333 334 335 336
    for(int i = 0; i < 1000000; i++) {
        getXY(x, y);
        fh << x << "\t" << y << "\n";
    }
gsell's avatar
gsell committed
337

338 339
    fh.close();
}
gsell's avatar
gsell committed
340

341 342 343 344
void LaserProfile::getXY(double &x, double &y) {
    double u = gsl_rng_uniform(rng_m);
    double v = gsl_rng_uniform(rng_m);
    gsl_histogram2d_pdf_sample(pdf_m, u, v, &x, &y);
gsell's avatar
gsell committed
345 346
}

347 348 349 350 351 352 353 354 355 356
void LaserProfile::getProfileMax(unsigned short &maxIntensity, unsigned short *image) {
    unsigned int numberPixels = sizeX_m * sizeY_m;

    maxIntensity = 0;

    for(unsigned int i = 0; i < numberPixels; i++) {
        if(image[i] > maxIntensity)
            maxIntensity = image[i];
    }
}