C code used for simulating SPIM imaging of a Drosophila embryo


#include <iostream>

#include <cmath>

//#include <pthread.h>

 

#define UseInt 3 // Set to 0 for normal scan, otherwise to n for number of phases (SI)

#define Photons 2 // Set to 1 for single-photon excitation, 2 for 2-photon excitation

 

#define UnifVol 0 // Set to 1 for uniform point distribution, 0 for tip-enriched distribution

#define EmbryoLength 30 // Embryo length (in um)

#define EmbryoWidth 15 // Embryo width (in um)

#define EmbryoThick 5 // Embryo thickness (in um)

#define EmbryoNpoints 50000 // Number of points (fluorescent emitters) inside the embryo

 

#define DetectPSF_Xsize 50 // X size (pixels) of the detection PSF contained in the file Detect_PSF.dat

#define DetectPSF_Ysize 50 // Y size (pixels) of the detection PSF contained in the file Detect_PSF.dat

#define DetectPSF_Zsize 50 // Number of planes of the detection PSF contained in the file Detect_PSF.dat

#define DetectPSF_scale 0.1 // Pixel size (in um) of the detection PSF contained in the file Detect_PSF.dat

 

#define IllumPSF_Xsize 50 // X size (pixels) of the illumination PSF contained in the file Illum_PSF.dat

#define IllumPSF_Ysize 50 // Y size (pixels) of the illumination PSF contained in the file Illum_PSF.dat

#define IllumPSF_Zsize 300 // Number of planes of the illumination PSF contained in the file Illum_PSF.dat

#define IllumPSF_scale 0.1 // Pixel size (in um) of the illumination PSF contained in the file Illum_PSF.dat

#define IntensityScale 1E8 // Scale for the Intensity (a. u.)

 

#define CCD_Xsize 35 // X size (pixels) of the simulated CCD image

#define CCD_Ysize 300 // Y size (pixels) of the simulated CCD image

#define ImageScale 0.1 // Pixel size (um) of the simulated CCD image

 

#define BeamYpos 220 // Number of Y positions scanned by the illumination beam

#define BeamZpos 1 // Number of Z positions scanned by the illumination beam

#define BeamYstep 0.1 // Distance (in um) between two Y positions scanned by the beam

#define BeamZstep 0.1 // Distance (in um) between two Z positions scanned by the beam

#define ThrDist 2.1 // Threshold distance (in um) to consider a point hit by the illumination PSF

#define SIperiod 0.5 // Period (in um) of the SI grating

 

 

void ExtractZPlane (double *inputArray, int Xsize, int Ysize, int Zsize, int whichPlane, double *outPlane) {

// Function that extracts a Z plane from a 3D array

int i, j;

int offset = (Xsize) * (Ysize) * whichPlane;

 

for (j=0; j<Ysize; j++) {

for (i=0; i<Xsize; i++) {

 

outPlane [i+j*Xsize] = inputArray[i+j*Xsize+offset];

//printf("Indexing_in= %d , Value = %f\n", i+j*Xsize+offset, inputArray[i+j*Xsize+offset]);

//printf("Indexing_out= %d\n", i+j*Xsize);

}

}

 

 

}

 

 

void ExtractYPlane (double *inputArray, int Xsize, int Ysize, int Zsize, int whichPlane, double *outPlane) {

// Function that extracts a Y plane from a 3D array

int i, j;

int offset = (Xsize) * whichPlane;

 

for (j=0; j<Zsize; j++) {

for (i=0; i<Xsize; i++) {

 

outPlane [i+j*Xsize] = inputArray[i+j*Xsize*Ysize+offset];

//printf("Indexing_in= %d , Value = %f\n", i+j*Xsize*Ysize+offset, inputArray[i+j*Xsize*Ysize+offset]);

//printf("Indexing_out= %d\n", i+j*Xsize);

}

}

 

 

}

 

void ExtractXPlane (double *inputArray, int Xsize, int Ysize, int Zsize, int whichPlane, double *outPlane) {

// Function that extracts a X plane from a 3D array

int i, j;

int offset = whichPlane;

 

for (j=0; j<Zsize; j++) {

for (i=0; i<Ysize; i++) {

 

outPlane [i+j*Ysize] = inputArray[i*Xsize+j*Xsize*Ysize+offset];

//printf("Indexing_in= %d , Value = %f\n", i*Xsize+j*Xsize*Ysize+offset, inputArray[i*Xsize+j*Xsize*Ysize+offset]);

//printf("Indexing_out= %d\n", i+j*Xsize);

}

}

 

 

}

 

 

 

int main (int argc, char * const argv[]) {

 

int i=0;

int j=0;

int k=0;

int l=0;

int m=0;

 

 

FILE *outFile_Debug;

FILE *inFile_DetectPSF;

FILE *inFile_IllumPSF;

FILE *outFile_CCDsim;

 

 

// Adding EmbryoNpoints randomnly inside an ellipsoid of major axis EmbryoLength and minor axis EmbryoWidth, EmbryoThick (using hit-or-miss)

 

double Embryo[EmbryoNpoints][3];

double TheX;

double TheY;

double TheZ;

double TheR1;

double TheR2;

 

outFile_Debug = fopen("Embryo.dat", "wb");

 

 

if (UnifVol == 1) {

// Uniform distribution

while(1){

 

TheX = (double) rand()/RAND_MAX - 0.5;

TheY = (double) rand()/RAND_MAX - 0.5;

TheZ = (double) rand()/RAND_MAX - 0.5;

 

if (sqrt(TheX*TheX+TheY*TheY+TheZ*TheZ) < (0.5)) {

TheX = TheX * EmbryoLength;

TheY = TheY * EmbryoWidth;

TheZ = TheZ * EmbryoThick;

 

Embryo[i][0] = TheX;

Embryo[i][1] = TheY;

Embryo[i][2] = TheZ;

 

i ++;

 

fwrite(&TheX, sizeof(double), 1, outFile_Debug);

fwrite(&TheY, sizeof(double), 1, outFile_Debug);

fwrite(&TheZ, sizeof(double), 1, outFile_Debug);

 

 

 

}

 

if (i == EmbryoNpoints) {

break;

}

 

}

}else{

 

// Tip-enriched placement

 

for (i=0; i<EmbryoNpoints; i++) {

 

TheX = (double) rand()/RAND_MAX * EmbryoLength - EmbryoLength/2;

 

TheR1 = EmbryoWidth/2 * sqrt (1-TheX*TheX/(EmbryoLength/2*EmbryoLength/2));

 

while (1) {

 

TheY = (double) rand()/RAND_MAX * EmbryoWidth - EmbryoWidth/2;

if(sqrt(TheY*TheY)<TheR1){

break;

}

}

 

TheR2 = EmbryoThick/2 * sqrt (1 - TheX*TheX/(EmbryoLength/2*EmbryoLength/2)) *sqrt (1 - TheY*TheY/(EmbryoWidth/2*EmbryoWidth/2));

while (1) {

 

TheZ = (double) rand()/RAND_MAX * EmbryoThick - EmbryoThick/2;

if(sqrt(TheZ*TheZ)<TheR2){

break;

}

}

 

Embryo[i][0] = TheX;

Embryo[i][1] = TheY;

Embryo[i][2] = TheZ;

 

fwrite(&TheX, sizeof(double), 1, outFile_Debug);

fwrite(&TheY, sizeof(double), 1, outFile_Debug);

fwrite(&TheZ, sizeof(double), 1, outFile_Debug);

}

 

 

}

 

// Saving the cartesian coordinates of the EmbryoNpoints to a file

fclose (outFile_Debug);

 

double TheDist;

double TheInt;

double Attenuation;

int DeFocus;

int Xoffset;

int Yoffset;

 

double SIphase;

int ThisPhase = 0;

 

 

double SIphaseOffset = 2*M_PI/UseInt*ThisPhase;

 

 

int TheN=0;

 

double Detect_PSF[DetectPSF_Zsize][DetectPSF_Ysize][DetectPSF_Xsize];

double Illum_PSF[IllumPSF_Zsize][IllumPSF_Ysize][IllumPSF_Xsize];

 

// Loading the detection PSF from the file Detect_PSF.dat (must be in the same directory of the program)

inFile_DetectPSF = fopen("Detect_PSF.dat", "rb");

 

for (k=0; k<DetectPSF_Zsize; k++) {

for (j=0; j<DetectPSF_Ysize; j++){

for (i=0; i<DetectPSF_Xsize; i++) {

 

fread(&Detect_PSF[k][j][i], sizeof(double), 1, inFile_DetectPSF);

}

}

 

}

 

fclose(inFile_DetectPSF);

 

// Loading the illumination PSF from the file Illum_PSF.dat (must be in the same directory of the program)

inFile_IllumPSF = fopen("Illum_PSF.dat", "rb");

 

for (k=0; k<IllumPSF_Zsize; k++) {

for (j=0; j<IllumPSF_Ysize; j++){

for (i=0; i<IllumPSF_Xsize; i++) {

 

fread(&Illum_PSF[k][j][i], sizeof(double), 1, inFile_IllumPSF);

if (Photons==2){

// Squaring (needed for 2-ph)

Illum_PSF[k][j][i] = Illum_PSF[k][j][i] * Illum_PSF[k][j][i];

}

}

}

 

}

 

fclose(inFile_IllumPSF);

 

//outFile_Debug = fopen("Embryo_N.dat", "wb");

 

 

// Allocating ThisImage (contains the simulated CCD image) and initializing to 0

double ThisImage[CCD_Ysize][CCD_Xsize];

 

for (j=0; j<CCD_Ysize; j++){

for (i=0; i<CCD_Xsize; i++){

 

ThisImage[j][i] = 0;

}

}

 

 

 

 

char fileName[19];

 

 

 

// Allocating ThisDotImage (contains the Detect_PSF plane that will be copied to ThisImage)

double ThisDotImage[DetectPSF_Ysize][DetectPSF_Xsize];

 

 

for(ThisPhase=0; ThisPhase<UseInt; ThisPhase++){

 

sprintf(fileName, "CCD_sim_phase%d.dat\0", ThisPhase);

 

outFile_CCDsim = fopen(fileName, "wb");

 

SIphaseOffset = 2*M_PI/UseInt*ThisPhase;

 

// Scanning the beam in the YZ plane

for (j=0; j<BeamZpos; j++) {

 

for (i=0; i<BeamYpos; i++) {

 

// Position of the beam

TheY = (-BeamYpos/2 + i)*BeamYstep;

TheZ = (-BeamZpos/2 + j)*BeamZstep;

 

//printf("Moving the beam to Y = %f Z = %f\n", TheY, TheZ);

 

// Checking how many points fall within ThrDist from the beam position

for (k=0; k<EmbryoNpoints; k++) {

 

// Calculating the YZ distance between the k-th Embryo point and the beam

TheDist=sqrt((TheY-Embryo[k][1])*(TheY-Embryo[k][1])+(TheZ-Embryo[k][2])*(TheZ-Embryo[k][2]));

 

if (TheDist < ThrDist){

// The k-th point was within ThrDist

//printf("Point Embryo[%d] at X = %f Y = %f Z =%f is %f from beam at Y = %f Z = %f\n", k , Embryo[k][0], Embryo[k][1], Embryo[k][2], TheDist, TheY, TheZ);

 

// Calculating the offset in pixel of the k-th point (needed to place the point on the CCD image)

// Yoffset on the CCD image corresponds to the plane of the Illum_PSF

// Xoffset on the CCD image corresponds to the Y distance between the beam position and the k-th point

Yoffset = (int) ((Embryo[k][0])/ImageScale+CCD_Ysize/2-DetectPSF_Ysize/2);

Xoffset = (int) ((TheY-Embryo[k][1])/ImageScale+CCD_Xsize/2-DetectPSF_Xsize/2);

 

// Calculating which plane of the Detect_PSF to extract (Z distance between beam position and k-th point)

DeFocus = (int) ((TheZ-Embryo[k][2])/DetectPSF_scale+DetectPSF_Zsize/2);

 

// Calculating the intensity of the illumination PSF at the k-th point location

 

if(UseInt==0) {

// Normal imaging

// Attenuation factor (???)

//Attenuation =  IntensityScale*(((double) (Yoffset))/CCD_Ysize);

//printf("Attenuation: %E\n", ((double) Yoffset)/CCD_Ysize);

//printf("Indexes: k=%d j=%d i=%d\n",(int) ((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale), (int) ((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2), (int) ((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2) );

Attenuation = IntensityScale;

}else {

// SI

SIphase = 2*M_PI*TheY/SIperiod;

Attenuation = IntensityScale * (1 + cos(SIphase+SIphaseOffset));

}

 

/*

// Debug: checking that the offsets make sense

if (((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale) < 0) {

printf("Z at k=%d is %E!!!\n", k, ((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale));

return -2;

}

 

if ( (int) ((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale) > IllumPSF_Zsize) {

printf("Z at k=%d is %E!!!\n", k, ((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale));

return -3;

}

 

if ((int) ((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2) < 0) {

printf("Y at k=%d is %E!!!\n", k, ((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2));

return -4;

}

 

if ((int)((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2) > IllumPSF_Ysize) {

printf("Y at k=%d is %E!!!\n", k, ((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2));

return -5;

}

 

if ((int)((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2) < 0) {

printf("X at k= %d is %E!!!\n", k, ((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2));

printf("TheZ = %E, TheDist = %E\n", TheZ, TheDist);

printf("Embryo = %E, %E, %E\n", Embryo[k][0], Embryo[k][1], Embryo[k][2]);

printf("IllumSize = %d IllumScale = %f\n", IllumPSF_Xsize, IllumPSF_scale);

return -6;

}

 

if ((int)((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2) > IllumPSF_Xsize) {

printf("X at k=%d  is %E!!!\n", k, ((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2));

return -7;

}

*/

 

TheInt = Attenuation * Illum_PSF[(int) ((Embryo[k][0]+EmbryoLength/2)/IllumPSF_scale)][(int) ((TheY-Embryo[k][1])/IllumPSF_scale+IllumPSF_Ysize/2)][(int) ((TheZ-Embryo[k][2])/IllumPSF_scale+IllumPSF_Xsize/2)];

 

// Extracting the de-focused Detect_PSF plane into ThisDotImage

ExtractZPlane (&Detect_PSF[0][0][0], DetectPSF_Xsize, DetectPSF_Ysize, DetectPSF_Zsize, DeFocus, &ThisDotImage[0][0]);

 

// Putting ThisDotImage into ThisImage

for (l=0; l<DetectPSF_Ysize; l++){

 

for (m=0; m<DetectPSF_Xsize; m++){

 

if (((l+Yoffset)>0)&&((l+Yoffset)<CCD_Ysize)&&((m+Xoffset)>0)&&((m+Xoffset)<CCD_Xsize) ) {

ThisImage[l+Yoffset][m+Xoffset] += TheInt*ThisDotImage[l][m];

 

}

}

}

 

// Adding 1 to the number of points hit by the beam

TheN++;

}

}

 

 

if(TheN>0){

printf("Found %d points at beam location Y = %.3f, Z = %.3f\n", TheN, TheY, TheZ);

printf("Phase %d: SI intensity is %.3f at phi=%.3f\n", ThisPhase, 1 + cos(SIphase+SIphaseOffset), SIphaseOffset);

 

}

 

//fwrite(&TheN, sizeof(int), 1, outFile_Debug);

 

TheN=0;

 

// Saving ThisImage to the output file and re-initializing it to 0

for(k=0; k<CCD_Ysize; k++){

for(l=0; l<CCD_Xsize; l++){

fwrite(&ThisImage[k][l], sizeof(double), 1, outFile_CCDsim);

ThisImage[k][l] = 0;

}

}

 

}

 

 

}

 

// Closing the output files

//fclose (outFile_Debug);

fclose(outFile_CCDsim);

 

}

 

 

 

 

}