SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
sddsspotanalysis.c
Go to the documentation of this file.
1/**
2 * @file sddsspotanalysis.c
3 * @brief Analysis of spot images from SDDS files.
4 *
5 * @details
6 * This program analyzes spot images provided in the SDDS (Self Describing Data Sets) file format.
7 * It allows users to define regions of interest, apply filters, compute statistical parameters,
8 * and save analysis results to new SDDS files. The program also supports options for outputting
9 * visualizations of the processed spots.
10 *
11 * @section Usage
12 * ```
13 * sddsspotanalysis <inputfile> <outputfile>
14 * [-pipe[=in][,out]]
15 * [-ROI=[{xy}{01}value=<value>][,{xy}{01}parameter=<name>]]
16 * [-spotROIsize=[{xy}value=<value>][,{xy}parameter=<name>]]
17 * [-centerOn={{xy}centroid | {xy}peak} | {xy}Parameter=<name>}]
18 * [-imageColumns=<listOfNames>]
19 * [-xyz=<ix>,<iy>,<Intensity>]
20 * [-singleSpot]
21 * [-levels=[intensity=<integer>][,saturation=<integer>]]
22 * [-blankOut=[{xy}{01}value=<value>][,{xy}{01}parameter=<name>]]
23 * [-sizeLines=[{xy}value=<value>][,{xy}parameter=<name>]]
24 * [-background=[halfwidth=<value>][,symmetric][,antihalo][,antiloner[,lonerThreshold=<value>]]
25 * [-despike=[neighbors=<integer>][,passes=<integer>][,averageOf=<integer>][,threshold=<value>][,keep]]
26 * [-hotpixelFilter=level=<fraction>,distance=<integer>[,passes=<integer>]]
27 * [-clipNegative]
28 * [-spotImage=<filename>]
29 * [-majorOrder=row|column]
30 * ```
31 *
32 * @section Options
33 * | Optional | Description |
34 * |--------------------------------------|----------------------------------------------------------------------------------------|
35 * | `-pipe` | Use standard SDDS Toolkit pipe option. |
36 * | `-ROI` | Define a region of interest in pixel units. |
37 * | `-spotROIsize` | Specify the ROI size for analyzing spot properties. |
38 * | `-centerOn` | Center analysis on peak value, centroid, or specified parameter for both axes. |
39 * | `-imageColumns` | List column names containing image data. |
40 * | `-xyz` | Define columns for x, y indices, and intensity. |
41 * | `-singleSpot` | Retain only the most intense connected pixels, eliminating multiple spots. |
42 * | `-levels` | Set intensity levels and saturation thresholds. |
43 * | `-blankOut` | Define regions to blank out. |
44 * | `-sizeLines` | Specify number of lines for beam size analysis. |
45 * | `-background` | Configure background subtraction settings. |
46 * | `-despike` | Apply noise reduction with despiking options. |
47 * | `-hotpixelFilter` | Apply a hot-pixel filter to smooth intensity values. |
48 * | `-clipNegative` | Set negative intensity values to zero. |
49 * | `-spotImage` | Save images of spots for visualization. |
50 * | `-majorOrder` | Specify the output file order (row or column major). |
51 *
52 * @subsection Incompatibilities
53 * - `-imageColumns` and `-xyz` are mutually exclusive options.
54 * - For `-despike`:
55 * - `averageOf` must not exceed the number of neighbors.
56 *
57 * @copyright
58 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
59 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
60 *
61 * @license
62 * This file is distributed under the terms of the Software License Agreement
63 * found in the file LICENSE included with this distribution.
64 *
65 * @author
66 * L. Emery, M. Borland, R. Soliday
67 */
68
69#include "mdb.h"
70#include "scan.h"
71#include "SDDS.h"
72#include "SDDSutils.h"
73#include "fftpackC.h"
74
75/* Enumeration for option types */
76enum option_type {
77 SET_ROI,
78 SET_SPOTROISIZE,
79 SET_IMAGECOLUMNS,
80 SET_DESPIKE,
81 SET_PIPE,
82 SET_SIZELINES,
83 SET_LEVELS,
84 SET_BLANKOUT,
85 SET_BACKGROUND,
86 SET_SPOTIMAGE,
87 SET_SINGLESPOT,
88 SET_CENTERON,
89 SET_MAJOR_ORDER,
90 SET_CLIP_NEGATIVE,
91 SET_XYZ,
92 SET_HOTPIXELS,
93 N_OPTIONS
94};
95
96static char *option[N_OPTIONS] = {
97 "roi",
98 "spotroi",
99 "imagecolumns",
100 "despike",
101 "pipe",
102 "sizeLines",
103 "levels",
104 "blankout",
105 "background",
106 "spotimage",
107 "singlespot",
108 "centeron",
109 "majorOrder",
110 "clipnegative",
111 "xyz",
112 "hotpixelfilter"
113};
114
115char *USAGE1 = "Usage: sddsspotanalysis [<input>] [<output>] [-pipe[=in][,out]] \n\
116[-ROI=[{xy}{01}value=<value>][,{xy}{01}parameter=<name>]]\n\
117[-spotROIsize=[{xy}value=<value>][,{xy}parameter=<name>]]\n\
118[-centerOn={{xy}centroid | {xy}peak} | {xy}Parameter=<name>}]\n\
119{-imageColumns=<listOfNames> | -xyz=<ix>,<iy>,<Intensity>} [-singleSpot]\n\
120[-levels=[intensityLevels=<integer>][,saturationLevel=<integer>]]\n\
121[-blankOut=[{xy}{01}value=<value>][,{xy}{01}parameter=<name>]]\n\
122[-sizeLines=[{xy}value=<value>][,{xy}parameter=<name>]]\n\
123[-background=[halfwidth=<value>][,symmetric][,antihalo][,antiloner[,lonerThreshold=<value>]]\n\
124[-despike=[neighbors=<integer>][,passes=<integer>][,averageOf=<integer>][,threshold=<value>][,keep]]\n\
125[-hotpixelFilter=level=<fraction>,distance=<integer>[,passes=<integer>]]\n\
126[-clipNegative] [-spotImage=<filename>] [-majorOrder=row|column] \n\n";
127char *USAGE2 = "Options:\n"
128 " -pipe[=in][,out] Use the standard SDDS Toolkit pipe option.\n"
129 " -ROI Define the region of interest in pixel units.\n"
130 " All data outside this region is ignored.\n"
131 " -spotROIsize Specify the size in pixels of the ROI around the spot.\n"
132 " This ROI is used for computing spot properties.\n"
133 " -imagecolumns <list> List names of columns containing image data.\n"
134 " Wildcards are supported.\n"
135 " -xyz <ix>,<iy>,<Intensity> Specify columns for x and y pixel indices and intensity.\n"
136 " -singleSpot Eliminate multiple spots by retaining only the most intense connected pixels.\n"
137 " -centerOn <method> Center the analysis on the peak value, centroid, or a specified parameter for both x and y axes.\n"
138 " -levels intensityLevels=<int>, saturationLevel=<int>\n"
139 " Set intensity levels and saturation level.\n"
140 " -blankOut <parameters> Define regions to blank out based on pixel values or parameters.\n"
141 " -sizeLines <parameters> Number of lines to use for analyzing the beam size. Default is 3.\n"
142 " -background <parameters> Configure background subtraction with options like halfwidth, symmetric,\n"
143 " antihalo, antiloner, and lonerThreshold.\n"
144 " -despike <parameters> Apply despiking to smooth the data with options for neighbors, passes,\n"
145 " averageOf, threshold, and keep.\n"
146 " -hotpixelFilter <parameters> Apply a hot pixel filter with level, distance, and passes parameters.\n"
147 " -clipNegative Set negative pixel values to zero.\n"
148 " -spotImage <filename> Specify a file to save images of the spots for plotting with sddscontour.\n"
149 " -majorOrder <row|column> Define the output file order as either row or column.\n\n";
150
151char *USAGE3 = "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
152
153#define DESPIKE_AVERAGEOF 0x0001U
154#define DESPIKE_KEEP (DESPIKE_AVERAGEOF << 1)
155
156#define X0VALUE 0x0001U
157#define X1VALUE (X0VALUE << 1)
158#define Y0VALUE (X1VALUE << 1)
159#define Y1VALUE (Y0VALUE << 1)
160#define X0PARAM (Y1VALUE << 1)
161#define X1PARAM (X0PARAM << 1)
162#define Y0PARAM (X1PARAM << 1)
163#define Y1PARAM (Y0PARAM << 1)
164#define OPTGIVEN (Y1PARAM << 1)
165
166typedef struct
167{
168 double backgroundLevel, integratedSpotIntensity, peakSpotIntensity;
169 double saturationCount, backgroundKilledPositive, backgroundKilledNegative;
170 int32_t ROI[4], spotROI[4], spotCenter[2];
171 double spotSigma[2], spotRange50[2], spotRange65[2], spotRange80[2], spotCentroid[2];
172 double S11, S33, rmsCrossProduct, phi, majorAxis, minorAxis;
173} IMAGE_ANALYSIS;
174
175typedef struct
176{
177 unsigned long flags;
178 double fraction;
179 long passes, distance;
180} HOTPIXEL_SETTINGS;
181
182void replaceWithNearNeighbors(double **image, long iy0, long iy1, long ix0, long ix1, long iyc, long ixc, long distance);
183
184typedef struct
185{
186 unsigned long flags;
187 int32_t neighbors, passes, averageOf;
188 double threshold;
189} DESPIKE_SETTINGS;
190
191long SetUpOutputFile(SDDS_DATASET *SDDSout, char *outputFile, SDDS_DATASET *SDDSin, char ***copyParamName,
192 long *copyParamNames, short columnMajorOrder);
193long DetermineQuadLongValues(int32_t *ROI, unsigned long flags, char **parameter, SDDS_DATASET *SDDSin, long nx, long ny, long isROI);
194long DetermineDualLongValues(int32_t *spotROISize, unsigned long flags, char **parameter, SDDS_DATASET *SDDSin, long defaultValue);
195void BlankOutImageData(double **image, long nx, long ny, int32_t *region);
196long GetImageData(long *nx, double ***image, char **imageColumn, long ny, SDDS_DATASET *SDDSin);
197long GetXYZImageData(long *nx, long *ny, double ***image, char *ixName, char *iyName, char *intensityName, SDDS_DATASET *SDDSin);
198
199#define SYMMETRIC_BGREMOVAL 0x0001U
200#define ANTIHALO_BGREMOVAL 0x0002U
201#define REMOVE_LONERS 0x0004U
202#define SINGLE_SPOT 0x0008U
203#define XCENTER_PEAK 0x0010U
204#define YCENTER_PEAK 0x0020U
205#define XCENTER_CENTROID 0x0040U
206#define YCENTER_CENTROID 0x0080U
207#define CLIP_NEGATIVE 0x0100U
208#define XCENTER_PARAM 0x0200U
209#define YCENTER_PARAM 0x0400U
210
211long AnalyzeImageData(double **image, long nx, long ny, int32_t *ROI, int32_t *spotROISize,
212 int32_t *sizeLines, DESPIKE_SETTINGS *despikeSettings, HOTPIXEL_SETTINGS *hotpixelSettings,
213 long intensityLevels, long saturationLevel, long backgroundHalfWidth, long lonerThreshold,
214 long lonerPasses, unsigned long flags, IMAGE_ANALYSIS *analysisResults, SDDS_DATASET *SDDSspim,
215 double *centerValue);
216
217void DetermineBeamSizes(double *sigma, double *spotRange50, double *spotRange65, double *spotRange80,
218 double *lineBuffer, long i0, long i1);
219void DetermineBeamParameters(double **image, long *spotROI, long nx, long ny, double *S11, double *S33,
220 double *rmsCrossProduct, double *phi, double *majorAxis, double *minorAxis);
221
222int main(int argc, char **argv) {
223 SCANNED_ARG *scArg;
224 long iArg, imageColumns = 0, i, ip;
225 char **imageColumn = NULL;
226 char *ixName, *iyName, *IntensityName;
227 short xyzMode = 0;
228 unsigned long pipeFlags = 0, ROIFlags = 0, spotROIFlags = 0, blankOutFlags = 0;
229 int32_t intensityLevels = 256, saturationLevel = 254, backgroundHalfWidth = 3, lonerThreshold = 1, lonerPasses = 1;
230 long despike = 0;
231 DESPIKE_SETTINGS despikeSettings;
232 HOTPIXEL_SETTINGS hotpixelSettings;
233 short hotpixelFilter = 0;
234 char *inputFile = NULL, *outputFile = NULL;
235 int32_t ROI[4] = {-1, -1, -1, -1};
236 int32_t spotROISize[2] = {-1, -1};
237 char *ROIParameter[4] = {NULL, NULL, NULL, NULL};
238 int32_t blankOut[4] = {-1, -1, -1, -1};
239 char *blankOutParameter[4] = {NULL, NULL, NULL, NULL};
240 char *spotROISizeParameter[2] = {NULL, NULL};
241 SDDS_DATASET SDDSin, SDDSout;
242 long readStatus, nx, ny, outputRow, maxPages = 0;
243 double **image;
244 IMAGE_ANALYSIS analysisResults;
245 unsigned long sizeLinesFlags = 0, dummyFlags = 0, backgroundFlags = 0, centerFlags = 0, majorOrderFlag;
246 int32_t sizeLines[2] = {-1, -1};
247 char *sizeLinesParameter[2] = {NULL, NULL};
248 char *centerParameter[2] = {NULL, NULL};
249 double centerValue[2];
250 char **copyParamName;
251 long copyParamNames, copyBuffer[32];
252 char *spotImageFile = NULL;
253 short columnMajorOrder = -1;
254 SDDS_DATASET SDDSspim;
255
256 if (argc < 2) {
257 fprintf(stderr, "%s%s%s", USAGE1, USAGE2, USAGE3);
258 exit(EXIT_FAILURE);
259 }
260
261 ixName = iyName = IntensityName = NULL;
262
263 argc = scanargsg(&scArg, argc, argv);
264 for (iArg = 1; iArg < argc; iArg++) {
265 if (scArg[iArg].arg_type == OPTION) {
266 delete_chars(scArg[iArg].list[0], "_");
267 /* process options here */
268 switch (match_string(scArg[iArg].list[0], option, N_OPTIONS, 0)) {
269 case SET_MAJOR_ORDER:
270 majorOrderFlag = 0;
271 scArg[iArg].n_items--;
272 if (scArg[iArg].n_items > 0 &&
273 (!scanItemList(&majorOrderFlag, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
274 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
275 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
276 SDDS_Bomb("invalid -majorOrder syntax/values");
277 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
278 columnMajorOrder = 1;
279 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
280 columnMajorOrder = 0;
281 break;
282 case SET_IMAGECOLUMNS:
283 if (scArg[iArg].n_items < 2)
284 SDDS_Bomb("invalid -imagecolumns syntax");
285 if (imageColumns)
286 SDDS_Bomb("give one and only one -imagecolumns option");
287 imageColumns = scArg[iArg].n_items - 1;
288 imageColumn = tmalloc(sizeof(*imageColumn) * imageColumns);
289 for (i = 0; i < imageColumns; i++)
290 imageColumn[i] = scArg[iArg].list[i + 1];
291 xyzMode = 0;
292 break;
293 case SET_ROI:
294 if (ROIFlags & OPTGIVEN)
295 SDDS_Bomb("give -ROI only once");
296 ROIFlags = OPTGIVEN;
297 scArg[iArg].n_items--;
298 if (scArg[iArg].n_items < 1 ||
299 (!scanItemList(&ROIFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
300 "x0value", SDDS_LONG, ROI + 0, 1, X0VALUE,
301 "x1value", SDDS_LONG, ROI + 1, 1, X1VALUE,
302 "y0value", SDDS_LONG, ROI + 2, 1, Y0VALUE,
303 "y1value", SDDS_LONG, ROI + 3, 1, Y1VALUE,
304 "x0parameter", SDDS_STRING, ROIParameter + 0, 1, X0PARAM,
305 "x1parameter", SDDS_STRING, ROIParameter + 1, 1, X1PARAM,
306 "y0parameter", SDDS_STRING, ROIParameter + 2, 1, Y0PARAM,
307 "y1parameter", SDDS_STRING, ROIParameter + 3, 1, Y1PARAM,
308 NULL) ||
309 ((ROIFlags & X0VALUE) && ROI[0] < 0) || ((ROIFlags & X1VALUE) && ROI[1] < 0) ||
310 ((ROIFlags & Y0VALUE) && ROI[2] < 0) || ((ROIFlags & Y1VALUE) && ROI[3] < 0) ||
311 ((ROIFlags & X0PARAM) && (!ROIParameter[0] || !strlen(ROIParameter[0]))) ||
312 ((ROIFlags & X1PARAM) && (!ROIParameter[1] || !strlen(ROIParameter[1]))) ||
313 ((ROIFlags & Y0PARAM) && (!ROIParameter[2] || !strlen(ROIParameter[2]))) ||
314 ((ROIFlags & Y1PARAM) && (!ROIParameter[3] || !strlen(ROIParameter[3])))))
315 SDDS_Bomb("invalid -ROI syntax/values");
316 break;
317 case SET_BLANKOUT:
318 if (blankOutFlags & OPTGIVEN)
319 SDDS_Bomb("give -blankout only once");
320 blankOutFlags = OPTGIVEN;
321 scArg[iArg].n_items--;
322 if (scArg[iArg].n_items < 1 ||
323 (!scanItemList(&blankOutFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
324 "x0value", SDDS_LONG, blankOut + 0, 1, X0VALUE,
325 "x1value", SDDS_LONG, blankOut + 1, 1, X1VALUE,
326 "y0value", SDDS_LONG, blankOut + 2, 1, Y0VALUE,
327 "y1value", SDDS_LONG, blankOut + 3, 1, Y1VALUE,
328 "x0parameter", SDDS_STRING, blankOutParameter + 0, 1, X0PARAM,
329 "x1parameter", SDDS_STRING, blankOutParameter + 1, 1, X1PARAM,
330 "y0parameter", SDDS_STRING, blankOutParameter + 2, 1, Y0PARAM,
331 "y1parameter", SDDS_STRING, blankOutParameter + 3, 1, Y1PARAM,
332 NULL) ||
333 ((blankOutFlags & X0VALUE) && blankOut[0] < 0) || ((blankOutFlags & X1VALUE) && blankOut[1] < 0) ||
334 ((blankOutFlags & Y0VALUE) && blankOut[2] < 0) || ((blankOutFlags & Y1VALUE) && blankOut[3] < 0) ||
335 ((blankOutFlags & X0PARAM) && (!blankOutParameter[0] || !strlen(blankOutParameter[0]))) ||
336 ((blankOutFlags & X1PARAM) && (!blankOutParameter[1] || !strlen(blankOutParameter[1]))) ||
337 ((blankOutFlags & Y0PARAM) && (!blankOutParameter[2] || !strlen(blankOutParameter[2]))) ||
338 ((blankOutFlags & Y1PARAM) && (!blankOutParameter[3] || !strlen(blankOutParameter[3])))))
339 SDDS_Bomb("invalid -blankOut syntax/values");
340 break;
341 case SET_SPOTROISIZE:
342 if (spotROIFlags & OPTGIVEN)
343 SDDS_Bomb("give -spotROIsize only once");
344 spotROIFlags = OPTGIVEN;
345 scArg[iArg].n_items--;
346 if (scArg[iArg].n_items < 1 ||
347 (!scanItemList(&spotROIFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
348 "xvalue", SDDS_LONG, spotROISize + 0, 1, X0VALUE,
349 "yvalue", SDDS_LONG, spotROISize + 1, 1, Y0VALUE,
350 "xparameter", SDDS_STRING, spotROISizeParameter + 0, 1, X0PARAM,
351 "yparameter", SDDS_STRING, spotROISizeParameter + 1, 1, Y0PARAM,
352 NULL) ||
353 ((spotROIFlags & X0VALUE) && spotROISize[0] < 0) ||
354 ((spotROIFlags & Y0VALUE) && spotROISize[1] < 0) ||
355 ((spotROIFlags & X0PARAM) && (!spotROISizeParameter[0] || !strlen(spotROISizeParameter[0]))) ||
356 ((spotROIFlags & Y0PARAM) && (!spotROISizeParameter[1] || !strlen(spotROISizeParameter[1])))))
357 SDDS_Bomb("invalid -spotROIsize syntax/values");
358 break;
359 case SET_SIZELINES:
360 if (sizeLinesFlags & OPTGIVEN)
361 SDDS_Bomb("give -sizeLines only once");
362 sizeLinesFlags = OPTGIVEN;
363 scArg[iArg].n_items--;
364 if (scArg[iArg].n_items < 1 ||
365 (!scanItemList(&sizeLinesFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
366 "xvalue", SDDS_LONG, sizeLines + 0, 1, X0VALUE,
367 "yvalue", SDDS_LONG, sizeLines + 1, 1, Y0VALUE,
368 "xparameter", SDDS_STRING, sizeLinesParameter + 0, 1, X0PARAM,
369 "yparameter", SDDS_STRING, sizeLinesParameter + 1, 1, Y0PARAM,
370 NULL) ||
371 ((sizeLinesFlags & X0VALUE) && sizeLines[0] < 0) ||
372 ((sizeLinesFlags & Y0VALUE) && sizeLines[1] < 0) ||
373 ((sizeLinesFlags & X0PARAM) && (!sizeLinesParameter[0] || !strlen(sizeLinesParameter[0]))) ||
374 ((sizeLinesFlags & Y0PARAM) && (!sizeLinesParameter[1] || !strlen(sizeLinesParameter[1])))))
375 SDDS_Bomb("invalid -sizeLines syntax/values");
376 break;
377 case SET_DESPIKE:
378 scArg[iArg].n_items--;
379 despikeSettings.neighbors = 4;
380 despikeSettings.passes = 1;
381 despikeSettings.threshold = 0;
382 despikeSettings.averageOf = 2;
383 despikeSettings.flags = 0;
384 if (scArg[iArg].n_items > 0 &&
385 (!scanItemList(&despikeSettings.flags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
386 "neighbors", SDDS_LONG, &despikeSettings.neighbors, 1, 0,
387 "passes", SDDS_LONG, &despikeSettings.passes, 1, 0,
388 "averageOf", SDDS_LONG, &despikeSettings.averageOf, 1, DESPIKE_AVERAGEOF,
389 "threshold", SDDS_DOUBLE, &despikeSettings.threshold, 1, 0,
390 "keep", -1, NULL, 0, DESPIKE_KEEP, NULL) ||
391 despikeSettings.neighbors < 2 || despikeSettings.passes < 1 ||
392 despikeSettings.averageOf < 2 || despikeSettings.threshold < 0))
393 SDDS_Bomb("invalid -despike syntax/values");
394 if (!(despikeSettings.flags & DESPIKE_AVERAGEOF))
395 despikeSettings.averageOf = despikeSettings.neighbors;
396 if (despikeSettings.averageOf > despikeSettings.neighbors)
397 SDDS_Bomb("invalid -despike syntax/values: averageOf>neighbors");
398 despike = 1;
399 break;
400 case SET_HOTPIXELS:
401 scArg[iArg].n_items--;
402 hotpixelSettings.passes = 1;
403 hotpixelSettings.distance = 1;
404 hotpixelSettings.fraction = -1;
405 if (scArg[iArg].n_items > 0 &&
406 (!scanItemList(&hotpixelSettings.flags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
407 "fraction", SDDS_DOUBLE, &hotpixelSettings.fraction, 1, 0,
408 "passes", SDDS_LONG, &hotpixelSettings.passes, 1, 0,
409 "distance", SDDS_LONG, &hotpixelSettings.distance, 1, 0, NULL) ||
410 hotpixelSettings.fraction <= 0 || hotpixelSettings.passes < 1 || hotpixelSettings.distance < 1))
411 SDDS_Bomb("invalid -hotpixelFilter syntax/values");
412 hotpixelFilter = 1;
413 break;
414 case SET_PIPE:
415 if (!processPipeOption(scArg[iArg].list + 1, scArg[iArg].n_items - 1, &pipeFlags))
416 SDDS_Bomb("invalid -pipe syntax");
417 break;
418 case SET_LEVELS:
419 scArg[iArg].n_items--;
420 if (scArg[iArg].n_items < 1 ||
421 !scanItemList(&dummyFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
422 "intensityLevels", SDDS_LONG, &intensityLevels, 1, 0,
423 "saturationLevel", SDDS_LONG, &saturationLevel, 1, 0, NULL) ||
424 intensityLevels <= 10 || saturationLevel <= 0)
425 SDDS_Bomb("invalid -levels syntax/values");
426 break;
427 case SET_BACKGROUND:
428 scArg[iArg].n_items--;
429 if (scArg[iArg].n_items < 1 ||
430 !scanItemList(&backgroundFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
431 "halfwidth", SDDS_LONG, &backgroundHalfWidth, 1, 0,
432 "symmetric", -1, NULL, 0, SYMMETRIC_BGREMOVAL,
433 "antihalo", -1, NULL, 0, ANTIHALO_BGREMOVAL,
434 "antiloner", -1, NULL, 0, REMOVE_LONERS,
435 "lonerthreshold", SDDS_LONG, &lonerThreshold, 1, 0,
436 "lonerpasses", SDDS_LONG, &lonerPasses, 1, 0, NULL) ||
437 backgroundHalfWidth < 0)
438 SDDS_Bomb("invalid -background syntax/values");
439 break;
440 case SET_SINGLESPOT:
441 if (scArg[iArg].n_items != 1)
442 SDDS_Bomb("invalid -singleSpot syntax/values");
443 backgroundFlags |= SINGLE_SPOT;
444 break;
445 case SET_SPOTIMAGE:
446 if (scArg[iArg].n_items != 2 || !(spotImageFile = scArg[iArg].list[1]))
447 SDDS_Bomb("invalid -spotImage syntax/values");
448 break;
449 case SET_CLIP_NEGATIVE:
450 backgroundFlags |= CLIP_NEGATIVE;
451 break;
452 case SET_CENTERON:
453 scArg[iArg].n_items--;
454
455 if (scArg[iArg].n_items < 1 ||
456 !scanItemList(&centerFlags, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
457 "xcentroid", -1, NULL, 0, XCENTER_CENTROID,
458 "ycentroid", -1, NULL, 0, YCENTER_CENTROID,
459 "xpeak", -1, NULL, 0, XCENTER_PEAK,
460 "ypeak", -1, NULL, 0, YCENTER_PEAK,
461 "xparameter", SDDS_STRING, centerParameter + 0, 1, XCENTER_PARAM,
462 "yparameter", SDDS_STRING, centerParameter + 1, 1, YCENTER_PARAM, NULL) ||
463 bitsSet(centerFlags & (XCENTER_CENTROID + XCENTER_PEAK + XCENTER_PARAM)) != 1 ||
464 bitsSet(centerFlags & (YCENTER_CENTROID + YCENTER_PEAK + YCENTER_PARAM)) != 1)
465 SDDS_Bomb("invalid -centerOn syntax");
466 break;
467 case SET_XYZ:
468 if (scArg[iArg].n_items != 4 ||
469 !strlen(ixName = scArg[iArg].list[1]) ||
470 !strlen(iyName = scArg[iArg].list[2]) ||
471 !strlen(IntensityName = scArg[iArg].list[3]))
472 SDDS_Bomb("invalid -xyz syntax");
473 xyzMode = 1;
474 break;
475 default:
476 fprintf(stderr, "unknown option %s given\n", scArg[iArg].list[0]);
477 exit(EXIT_FAILURE);
478 break;
479 }
480 } else {
481 /* argument is input filename */
482 if (!inputFile) {
483 inputFile = scArg[iArg].list[0];
484 } else if (!outputFile) {
485 outputFile = scArg[iArg].list[0];
486 } else
487 SDDS_Bomb("too many filenames");
488 }
489 }
490
491 processFilenames("sddsspotanalysis", &inputFile, &outputFile, pipeFlags, 0, NULL);
492
493 if (!imageColumns && !IntensityName)
494 SDDS_Bomb("you must give either the -imageColumns or -xyz option");
495 if (imageColumns && IntensityName)
496 SDDS_Bomb("you must give either the -imageColumns or -xyz option, not both");
497
498 if (!SDDS_InitializeInput(&SDDSin, inputFile)) {
499 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
500 exit(EXIT_FAILURE);
501 }
502
503 if (spotImageFile &&
504 (!SDDS_InitializeOutput(&SDDSspim, SDDS_BINARY, 0, NULL, "sddsspotanalysis spot image", spotImageFile) ||
505 SDDS_DefineColumn(&SDDSspim, "ix", NULL, NULL, NULL, NULL, SDDS_SHORT, 0) < 0 ||
506 SDDS_DefineColumn(&SDDSspim, "iy", NULL, NULL, NULL, NULL, SDDS_SHORT, 0) < 0 ||
507 SDDS_DefineColumn(&SDDSspim, "Image", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
508 SDDS_DefineParameter(&SDDSspim, "nx", NULL, NULL, NULL, NULL, SDDS_SHORT, 0) < 0 ||
509 SDDS_DefineParameter(&SDDSspim, "ny", NULL, NULL, NULL, NULL, SDDS_SHORT, 0) < 0 ||
510 !SDDS_WriteLayout(&SDDSspim))) {
511 SDDS_SetError("Problem setting up spot image output file.");
512 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
513 exit(EXIT_FAILURE);
514 }
515
516 if (!xyzMode) {
517 if ((imageColumns = expandColumnPairNames(&SDDSin, &imageColumn, NULL, imageColumns, NULL, 0, FIND_NUMERIC_TYPE, 0)) <= 0) {
518 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
519 SDDS_Bomb("No quantities selected to analyze");
520 }
521 ny = imageColumns;
522 outputRow = 0;
523 if (!(image = calloc(ny, sizeof(*image))))
524 SDDS_Bomb("Memory allocation failure");
525 nx = 0;
526 } else
527 ny = nx = 0;
528
529 if (!SetUpOutputFile(&SDDSout, outputFile, &SDDSin, &copyParamName, &copyParamNames, columnMajorOrder) ||
530 !SDDS_StartPage(&SDDSout, maxPages = 10))
531 SDDS_Bomb("Problem setting up output file.");
532
533 while ((readStatus = SDDS_ReadPage(&SDDSin)) > 0) {
534 if (readStatus > maxPages) {
535 if (!SDDS_LengthenTable(&SDDSout, 10))
536 SDDS_Bomb("Problem lengthening output file.");
537 maxPages += 10;
538 }
539 /* Get image into array.
540 * N.B.: pixel (ix, iy) is accessed as image[iy][ix].
541 */
542 if (!xyzMode) {
543 if (!GetImageData(&nx, &image, imageColumn, ny, &SDDSin))
544 SDDS_Bomb("Problem getting image data.");
545 if (!nx)
546 continue;
547 } else {
548 if (!GetXYZImageData(&nx, &ny, &image, ixName, iyName, IntensityName, &SDDSin))
549 SDDS_Bomb("Problem getting image data.");
550 if (!nx || !ny)
551 continue;
552 }
553 if (!DetermineQuadLongValues(ROI, ROIFlags, ROIParameter, &SDDSin, nx, ny, 1)) {
554 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
555 SDDS_Bomb("Problem determining region of interest---verify that you gave sufficient information.");
556 }
557 if (blankOutFlags && !DetermineQuadLongValues(blankOut, blankOutFlags, blankOutParameter, &SDDSin, nx, ny, 0))
558 SDDS_Bomb("Problem determining blankout region---verify that you gave sufficient information.");
559 if (!DetermineDualLongValues(spotROISize, spotROIFlags, spotROISizeParameter, &SDDSin, 150))
560 SDDS_Bomb("Problem determine size of spot region of interest---verify that you gave sufficient information.");
561 if (!DetermineDualLongValues(sizeLines, sizeLinesFlags, sizeLinesParameter, &SDDSin, 3))
562 SDDS_Bomb("Problem determine size of number of lines to use for spot size computation---verify that you gave sufficient information.");
563 centerValue[0] = centerValue[1] = -1;
564 if (centerFlags & XCENTER_PARAM && !SDDS_GetParameterAsDouble(&SDDSin, centerParameter[0], centerValue + 0))
565 SDDS_Bomb("Problem getting center parameter value for x.");
566 if (centerFlags & YCENTER_PARAM && !SDDS_GetParameterAsDouble(&SDDSin, centerParameter[1], centerValue + 1))
567 SDDS_Bomb("Problem getting center parameter value for y.");
568 if (blankOutFlags)
569 BlankOutImageData(image, nx, ny, blankOut);
570 if (!AnalyzeImageData(image, nx, ny, ROI, spotROISize, sizeLines,
571 despike ? &despikeSettings : NULL, hotpixelFilter ? &hotpixelSettings : NULL, intensityLevels,
572 saturationLevel, backgroundHalfWidth, lonerThreshold, lonerPasses, backgroundFlags | centerFlags,
573 &analysisResults, spotImageFile ? &SDDSspim : NULL, centerValue))
574 continue;
575 for (ip = 0; ip < copyParamNames; ip++) {
576 if (!SDDS_GetParameter(&SDDSin, copyParamName[ip], copyBuffer))
577 SDDS_Bomb("Problem reading parameter data from input file.");
578 if (!SDDS_SetRowValues(&SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_REFERENCE, outputRow, copyParamName[ip], copyBuffer, NULL)) {
579 SDDS_Bomb("Problem copying parameter data from input file.");
580 }
581 }
582 if (!SDDS_SetRowValues(&SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, outputRow++,
583 "ImageIndex", readStatus - 1,
584 "xROILower", analysisResults.ROI[0],
585 "xROIUpper", analysisResults.ROI[1],
586 "yROILower", analysisResults.ROI[2],
587 "yROIUpper", analysisResults.ROI[3],
588 "xSpotROILower", analysisResults.spotROI[0],
589 "xSpotROIUpper", analysisResults.spotROI[1],
590 "ySpotROILower", analysisResults.spotROI[2],
591 "ySpotROIUpper", analysisResults.spotROI[3],
592 "xSpotCenter", analysisResults.spotCenter[0],
593 "ySpotCenter", analysisResults.spotCenter[1],
594 "xSpotCentroid", analysisResults.spotCentroid[0],
595 "ySpotCentroid", analysisResults.spotCentroid[1],
596 "xSpotSigma", analysisResults.spotSigma[0],
597 "xSpotRange50", analysisResults.spotRange50[0],
598 "xSpotRange65", analysisResults.spotRange65[0],
599 "xSpotRange80", analysisResults.spotRange80[0],
600 "ySpotSigma", analysisResults.spotSigma[1],
601 "ySpotRange50", analysisResults.spotRange50[1],
602 "ySpotRange65", analysisResults.spotRange65[1],
603 "ySpotRange80", analysisResults.spotRange80[1],
604 "BackgroundLevel", analysisResults.backgroundLevel,
605 "BackgroundKilledPositive", analysisResults.backgroundKilledPositive,
606 "BackgroundKilledNegative", analysisResults.backgroundKilledNegative,
607 "IntegratedSpotIntensity", analysisResults.integratedSpotIntensity,
608 "PeakSpotIntensity", analysisResults.peakSpotIntensity,
609 "SaturationCount", analysisResults.saturationCount,
610 "phi", analysisResults.phi,
611 "rmsCrossProduct", analysisResults.rmsCrossProduct,
612 "majorAxis", analysisResults.majorAxis,
613 "minorAxis", analysisResults.minorAxis,
614 "S11", analysisResults.S11,
615 "S33", analysisResults.S33,
616 "S13", analysisResults.rmsCrossProduct, NULL)) {
617 SDDS_SetError("Problem setting values into output file.");
618 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
619 }
620 }
621 if (!SDDS_WritePage(&SDDSout) || !SDDS_Terminate(&SDDSout) || !SDDS_Terminate(&SDDSin))
622 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
623 if (spotImageFile && !SDDS_Terminate(&SDDSspim))
624 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
625 for (i = 0; i < ny; i++)
626 free(image[i]);
627 if (!xyzMode) {
628 for (i = 0; i < ny; i++)
629 free(imageColumn[i]);
630 free(imageColumn);
631 }
632 for (i = 0; i < copyParamNames; i++)
633 free(copyParamName[i]);
634 free(copyParamName);
635 free(image);
636 return EXIT_SUCCESS;
637}
638
639long SetUpOutputFile(SDDS_DATASET *SDDSout, char *outputFile, SDDS_DATASET *SDDSin, char ***copyParamName,
640 long *copyParamNames, short columnMajorOrder) {
641 char **paramName;
642 int32_t paramNames;
643
644 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsspotanalysis output", outputFile) ||
645 SDDS_DefineColumn(SDDSout, "xROILower", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
646 SDDS_DefineColumn(SDDSout, "xROIUpper", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
647 SDDS_DefineColumn(SDDSout, "xSpotROILower", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
648 SDDS_DefineColumn(SDDSout, "xSpotROIUpper", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
649 SDDS_DefineColumn(SDDSout, "xSpotCenter", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
650 SDDS_DefineColumn(SDDSout, "xSpotCentroid", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
651 SDDS_DefineColumn(SDDSout, "xSpotSigma", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
652 SDDS_DefineColumn(SDDSout, "xSpotRange50", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
653 SDDS_DefineColumn(SDDSout, "xSpotRange65", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
654 SDDS_DefineColumn(SDDSout, "xSpotRange80", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
655 SDDS_DefineColumn(SDDSout, "yROILower", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
656 SDDS_DefineColumn(SDDSout, "yROIUpper", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
657 SDDS_DefineColumn(SDDSout, "ySpotROILower", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
658 SDDS_DefineColumn(SDDSout, "ySpotROIUpper", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
659 SDDS_DefineColumn(SDDSout, "ySpotCenter", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
660 SDDS_DefineColumn(SDDSout, "ySpotCentroid", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
661 SDDS_DefineColumn(SDDSout, "ySpotSigma", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
662 SDDS_DefineColumn(SDDSout, "ySpotRange50", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
663 SDDS_DefineColumn(SDDSout, "ySpotRange65", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
664 SDDS_DefineColumn(SDDSout, "ySpotRange80", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
665 SDDS_DefineColumn(SDDSout, "ImageIndex", NULL, NULL, NULL, NULL, SDDS_LONG, 0) < 0 ||
666 SDDS_DefineColumn(SDDSout, "BackgroundLevel", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
667 SDDS_DefineColumn(SDDSout, "BackgroundKilledNegative", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
668 SDDS_DefineColumn(SDDSout, "BackgroundKilledPositive", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
669 SDDS_DefineColumn(SDDSout, "IntegratedSpotIntensity", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
670 SDDS_DefineColumn(SDDSout, "PeakSpotIntensity", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
671 SDDS_DefineColumn(SDDSout, "SaturationCount", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
672 SDDS_DefineColumn(SDDSout, "phi", NULL, "degree", NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
673 SDDS_DefineColumn(SDDSout, "rmsCrossProduct", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
674 SDDS_DefineColumn(SDDSout, "majorAxis", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
675 SDDS_DefineColumn(SDDSout, "minorAxis", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
676 SDDS_DefineColumn(SDDSout, "S11", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
677 SDDS_DefineColumn(SDDSout, "S33", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0 ||
678 SDDS_DefineColumn(SDDSout, "S13", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0) {
679 SDDS_SetError("Problem setting up output file.");
680 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
681 exit(EXIT_FAILURE);
682 }
683 if (columnMajorOrder != -1)
684 SDDSout->layout.data_mode.column_major = columnMajorOrder;
685 else
686 SDDSout->layout.data_mode.column_major = SDDSin->layout.data_mode.column_major;
687
688 *copyParamNames = 0;
689 *copyParamName = NULL;
690 if ((paramName = SDDS_GetParameterNames(SDDSin, &paramNames)) && paramNames > 0) {
691 long ip;
692 if (!(*copyParamName = calloc(sizeof(**copyParamName), paramNames)))
693 SDDS_Bomb("memory allocation failure");
694 for (ip = 0; ip < paramNames; ip++) {
695 if (SDDS_GetColumnIndex(SDDSout, paramName[ip]) >= 0) {
696 free(paramName[ip]);
697 continue;
698 }
699 if (!SDDS_DefineColumnLikeParameter(SDDSout, SDDSin, paramName[ip], NULL)) {
700 SDDS_SetError("Problem setting up output file.");
701 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
702 exit(EXIT_FAILURE);
703 }
704 (*copyParamName)[*copyParamNames] = paramName[ip];
705 (*copyParamNames)++;
706 }
707 }
708 free(paramName);
709 if (!SDDS_WriteLayout(SDDSout)) {
710 SDDS_SetError("Problem setting up output file.");
711 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
712 exit(EXIT_FAILURE);
713 }
714 return 1;
715}
716
717long DetermineQuadLongValues(int32_t *ROI, unsigned long flags, char **parameter, SDDS_DATASET *SDDSin, long nx, long ny, long isROI) {
718 long i;
719 double value;
720 double defaultROI[4];
721 if (isROI) {
722 defaultROI[0] = defaultROI[2] = 0;
723 defaultROI[1] = nx - 1;
724 defaultROI[3] = ny - 1;
725 } else
726 defaultROI[0] = defaultROI[2] = defaultROI[1] = defaultROI[3] = -1;
727 for (i = 0; i < 4; i++) {
728 if (flags & (X0VALUE << i))
729 continue;
730 if (flags & (X0PARAM << i)) {
731 if (!SDDS_GetParameterAsDouble(SDDSin, parameter[i], &value)) {
732 SDDS_SetError("parameter is nonexistent or nonnumeric");
733 return 0;
734 }
735 ROI[i] = value + 0.5;
736 } else
737 ROI[i] = defaultROI[i];
738 }
739 if (ROI[0] < 0 || ROI[2] < 0) {
740 SDDS_SetError("lower bound of region less than zero");
741 return 0;
742 }
743 if (ROI[1] > nx - 1 || ROI[3] > ny - 1) {
744 SDDS_SetError("upper bound of region too large");
745 return 0;
746 }
747 if (ROI[0] >= ROI[1]) {
748 SDDS_SetError("x region has zero or negative width");
749 return 0;
750 }
751 if (ROI[2] >= ROI[3]) {
752 SDDS_SetError("y region has zero or negative width");
753 return 0;
754 }
755 return 1;
756}
757
758long DetermineDualLongValues(int32_t *returnValue, unsigned long flags, char **parameter, SDDS_DATASET *SDDSin, long defaultValue) {
759 long i;
760 double value;
761
762 for (i = 0; i < 2; i++) {
763 if (flags & (X0VALUE << (2 * i)))
764 continue;
765 if (flags & (X0PARAM << (2 * i))) {
766 if (!SDDS_GetParameterAsDouble(SDDSin, parameter[i], &value)) {
767 SDDS_SetError("parameter is nonexistent or nonnumeric");
768 return 0;
769 }
770 returnValue[i] = (int32_t)value;
771 } else
772 returnValue[i] = defaultValue;
773 }
774 if (returnValue[0] <= 0 || returnValue[1] <= 0) {
775 SDDS_SetError("determined value is nonpositive");
776 return 0;
777 }
778 return 1;
779}
780
781long GetImageData(long *nx, double ***image, char **imageColumn, long ny, SDDS_DATASET *SDDSin) {
782 /* Note that because each column is a horizontal line, the pixels are
783 * accessed as image[iy][ix]. I could swap it, but that would mean
784 * doubling the amount of memory used.
785 */
786 long i;
787 *nx = 0;
788 for (i = 0; i < ny; i++) {
789 if ((*image)[i]) {
790 free((*image)[i]);
791 (*image)[i] = NULL;
792 }
793 if (!((*image)[i] = SDDS_GetColumnInDoubles(SDDSin, imageColumn[i]))) {
794 SDDS_Bomb("Unable to get data from columns");
795 }
796 }
797 *nx = SDDS_RowCount(SDDSin);
798 return 1;
799}
800
801long GetXYZImageData(long *nx, long *ny, double ***image, char *ixName, char *iyName, char *intensityName, SDDS_DATASET *SDDSin) {
802 int32_t *ixData, *iyData;
803 double *intensityData;
804 long ixMin, iyMin, ixMax, iyMax, i, ix, iy, nRows;
805 char *ixIndexCheck, *iyIndexCheck;
806
807 if ((nRows = SDDS_RowCount(SDDSin)) <= 0 ||
808 !(ixData = SDDS_GetColumnInLong(SDDSin, ixName)) ||
809 !(iyData = SDDS_GetColumnInLong(SDDSin, iyName)) ||
810 !(intensityData = SDDS_GetColumnInDoubles(SDDSin, intensityName)))
811 return 0;
812
813 ixMax = iyMax = -(ixMin = iyMin = LONG_MAX);
814
815 for (i = 0; i < nRows; i++) {
816 if (ixMax < ixData[i])
817 ixMax = ixData[i];
818 if (ixMin > ixData[i])
819 ixMin = ixData[i];
820 if (iyMax < iyData[i])
821 iyMax = iyData[i];
822 if (iyMin > iyData[i])
823 iyMin = iyData[i];
824 }
825 if (ixMax == ixMin)
826 return 0;
827 if (iyMax == iyMin)
828 return 0;
829
830 ixIndexCheck = calloc((*nx = ixMax - ixMin + 1), sizeof(*ixIndexCheck));
831 iyIndexCheck = calloc((*ny = iyMax - iyMin + 1), sizeof(*iyIndexCheck));
832
833 *image = calloc(*ny, sizeof(**image));
834 for (iy = 0; iy < *ny; iy++)
835 (*image)[iy] = calloc(*nx, sizeof(***image));
836
837 for (i = 0; i < nRows; i++) {
838 ix = ixData[i] - ixMin;
839 iy = iyData[i] - iyMin;
840 ixIndexCheck[ix] = 1;
841 iyIndexCheck[iy] = 1;
842 (*image)[iy][ix] = intensityData[i];
843 }
844
845 for (i = 0; i <= (ixMax - ixMin); i++) { // Changed from < to <=
846 if (!ixIndexCheck[i]) {
847 fprintf(stderr, "Error: image file is missing some x index values\n");
848 /* Free allocated memory before returning */
849 for (long y = 0; y < *ny; y++) {
850 free((*image)[y]);
851 }
852 free(*image);
853 free(ixIndexCheck);
854 free(iyIndexCheck);
855 free(ixData);
856 free(iyData);
857 free(intensityData);
858 return 0;
859 }
860 }
861 for (i = 0; i <= (iyMax - iyMin); i++) { // Changed from < to <=
862 if (!iyIndexCheck[i]) {
863 fprintf(stderr, "Error: image file is missing some y index values\n");
864 /* Free allocated memory before returning */
865 for (long y = 0; y < *ny; y++) {
866 free((*image)[y]);
867 }
868 free(*image);
869 free(ixIndexCheck);
870 free(iyIndexCheck);
871 free(ixData);
872 free(iyData);
873 free(intensityData);
874 return 0;
875 }
876 }
877
878 free(ixIndexCheck);
879 free(iyIndexCheck);
880 free(ixData);
881 free(iyData);
882 free(intensityData);
883
884 return 1;
885}
886
887long AnalyzeImageData(double **image, long nx, long ny, int32_t *ROI, int32_t *spotROISize, int32_t *sizeLines,
888 DESPIKE_SETTINGS *despikeSettings, HOTPIXEL_SETTINGS *hotpixelSettings,
889 long intensityLevels, long saturationLevel, long backgroundHalfWidth, long lonerThreshold,
890 long lonerPasses, unsigned long flags, IMAGE_ANALYSIS *analysisResults,
891 SDDS_DATASET *SDDSspim, double *centerValue) {
892 long ix, iy, iy0, iy1, ix0, ix1, nxROI, nyROI;
893 double maxValue, minValue, background, value, sum1, sum2;
894 long ixCenter, iyCenter;
895 int64_t ixMax, ixMin;
896 int32_t i;
897 double *intensityHistogram, *lineBuffer = NULL, *x = NULL, *y = NULL;
898 long spotROI[4];
899 double mean, sum;
900
901 /* Steps in image analysis:
902 * . "Apply" the ROI
903 * . Remove hot pixels, if requested
904 * . Find the spot:
905 * optionally despike the image to remove noise
906 * scan each line for maximum until the maximum maximum is found
907 * . Determine and subtract the background:
908 * accumulate histogram of pixel intensities.
909 * find the mode
910 * . If requested, run "single-spot" filter.
911 * . Check that ixCenter and iyCenter are consistent with having the
912 * spot ROI inside the pixel map. If not, adjust.
913 * . Sum over the spot ROI and subtract background*nx*ny.
914 */
915
916 /* set the ROI */
917 if (ROI[0] >= 0 && ROI[1] >= 0 && ROI[0] < ROI[1]) {
918 ix0 = ROI[0];
919 ix1 = ROI[1];
920 if (ix1 >= nx)
921 ix1 = nx - 1;
922 } else {
923 ix0 = 0;
924 ix1 = nx - 1;
925 }
926 nxROI = ix1 - ix0 + 1;
927
928 if (ROI[2] >= 0 && ROI[3] >= 0 && ROI[2] < ROI[3]) {
929 iy0 = ROI[2];
930 iy1 = ROI[3];
931 if (iy1 >= ny)
932 iy1 = ny - 1;
933 } else {
934 iy0 = 0;
935 iy1 = ny - 1;
936 }
937 nyROI = iy1 - iy0 + 1;
938
939 if ((nyROI) < spotROISize[1] || (nxROI) < spotROISize[0])
940 SDDS_Bomb("spot ROI is larger than ROI");
941
942 /* Find the spot */
943 if (despikeSettings)
944 if (!(lineBuffer = SDDS_Malloc(sizeof(*lineBuffer) * (nx > ny ? nx : ny))))
945 SDDS_Bomb("memory allocation failure");
946 minValue = DBL_MAX;
947 maxValue = -DBL_MAX;
948 ixCenter = iyCenter = -1;
949
950 for (ix = ix0; ix <= ix1; ix++) {
951 for (iy = iy0; iy <= iy1; iy++) {
952 if ((flags & CLIP_NEGATIVE) && image[iy][ix] < 0)
953 image[iy][ix] = 0;
954 if (image[iy][ix] < 0 || image[iy][ix] >= intensityLevels)
955 SDDS_Bomb("image intensity outside intensity level range");
956 }
957 }
958
959 if (despikeSettings && (despikeSettings->flags & DESPIKE_KEEP)) {
960 /* despike vertical lines */
961 for (ix = ix0; ix <= ix1; ix++) {
962 for (iy = iy0; iy <= iy1; iy++)
963 lineBuffer[iy - iy0] = image[iy][ix];
964 despikeData(lineBuffer, nyROI, despikeSettings->neighbors, despikeSettings->passes, despikeSettings->averageOf,
965 despikeSettings->threshold, 0);
966 for (iy = iy0; iy <= iy1; iy++)
967 image[iy][ix] = lineBuffer[iy - iy0];
968 }
969 }
970
971 /* despike horizontal lines, also find point of max intensity */
972 for (iy = iy0; iy <= iy1; iy++) {
973 if (despikeSettings) {
974 memcpy(lineBuffer, image[iy] + ix0, nxROI * sizeof(*lineBuffer));
975 despikeData(lineBuffer, nxROI, despikeSettings->neighbors, despikeSettings->passes, despikeSettings->averageOf,
976 despikeSettings->threshold, 0);
977 if (despikeSettings->flags & DESPIKE_KEEP)
978 memcpy(image[iy] + ix0, lineBuffer, nxROI * sizeof(*lineBuffer));
979 } else
980 lineBuffer = image[iy] + ix0;
981 index_min_max(&ixMin, &ixMax, lineBuffer, nxROI);
982 if (lineBuffer[ixMax] > maxValue) {
983 maxValue = lineBuffer[ixMax];
984 ixCenter = ixMax + ix0; /* since lineBuffer copy starts at ix=ix0 */
985 iyCenter = iy;
986 }
987 if (lineBuffer[ixMin] < minValue)
988 minValue = lineBuffer[ixMin];
989 }
990 if (despikeSettings) {
991 free(lineBuffer);
992 lineBuffer = NULL;
993 }
994 if (flags & XCENTER_PARAM)
995 ixCenter = centerValue[0];
996 else
997 centerValue[0] = ixCenter;
998 if (flags & YCENTER_PARAM)
999 iyCenter = centerValue[1];
1000 else
1001 centerValue[1] = iyCenter;
1002 if (ixCenter == -1 || iyCenter == -1)
1003 return 0;
1004
1005 /* Determine the background */
1006 if (!(intensityHistogram = calloc(intensityLevels, sizeof(*intensityHistogram))))
1007 SDDS_Bomb("memory allocation failure");
1008 for (iy = iy0; iy <= iy1; iy++)
1009 make_histogram(intensityHistogram, intensityLevels, -0.5, intensityLevels + 0.5, image[iy] + ix0, nxROI, iy == iy0);
1010 index_min_max(&ixMin, &ixMax, intensityHistogram, intensityLevels);
1011 sum1 = sum2 = 0;
1012 for (i = ixMax - backgroundHalfWidth; i <= ixMax + backgroundHalfWidth; i++) {
1013 if (i < 0 || i >= intensityLevels)
1014 continue;
1015 sum1 += intensityHistogram[i];
1016 sum2 += intensityHistogram[i] * i;
1017 }
1018 if (sum1)
1019 background = sum2 / sum1;
1020 else
1021 background = ixMax;
1022 if (background < 0)
1023 background = 0;
1024 free(intensityHistogram);
1025
1026 if (hotpixelSettings) {
1027 long min, max, pass;
1028 /* remove hot pixels:
1029 * 1. Histogram the intensity levels
1030 * 2. Identify pixels with levels that are >defined fraction below max
1031 * 3. Replace identified intensities with average of near neighbors
1032 */
1033 pass = hotpixelSettings->passes;
1034 while (pass--) {
1035 max = -(min = LONG_MAX);
1036 for (iy = iy0; iy <= iy1; iy++) {
1037 for (ix = ix0; ix <= ix1; ix++) {
1038 if (image[iy][ix] > max)
1039 max = image[iy][ix];
1040 if (image[iy][ix] < min)
1041 min = image[iy][ix];
1042 }
1043 }
1044 if (min >= max)
1045 SDDS_Bomb("Can't apply hotpixel filter (min=max). Are data non-integer?");
1046 for (iy = iy0; iy <= iy1; iy++)
1047 for (ix = ix0; ix <= ix1; ix++)
1048 if ((image[iy][ix] - min) > (max - min) * hotpixelSettings->fraction)
1049 replaceWithNearNeighbors(image, iy0, iy1, ix0, ix1, iy, ix, hotpixelSettings->distance);
1050 }
1051 }
1052
1053 /* Compute new ROI for spot only */
1054 if (flags & XCENTER_CENTROID) {
1055 mean = sum = 0;
1056 for (ix = ix0; ix <= ix1; ix++)
1057 for (iy = iy0; iy <= iy1; iy++) {
1058 sum += image[iy][ix];
1059 mean += image[iy][ix] * ix;
1060 }
1061 mean /= sum;
1062 spotROI[0] = mean - spotROISize[0] / 2;
1063 } else {
1064 spotROI[0] = ixCenter - spotROISize[0] / 2;
1065 }
1066 spotROI[1] = spotROI[0] + spotROISize[0] - 1;
1067 if (spotROI[0] < ix0) {
1068 spotROI[0] = ix0;
1069 spotROI[1] = ix0 + spotROISize[0] - 1;
1070 } else if (spotROI[1] > ix1) {
1071 spotROI[1] = ix1;
1072 spotROI[0] = ix1 - spotROISize[0] + 1;
1073 }
1074 if (spotROI[0] < ix0 || spotROI[1] > ix1)
1075 SDDS_Bomb("spot ROI is larger than ROI for x"); /* shouldn't ever see this */
1076
1077 if (flags & XCENTER_CENTROID) {
1078 mean = sum = 0;
1079 for (ix = ix0; ix <= ix1; ix++)
1080 for (iy = iy0; iy <= iy1; iy++) {
1081 sum += image[iy][ix];
1082 mean += image[iy][ix] * iy;
1083 }
1084 mean /= sum;
1085 spotROI[2] = mean - spotROISize[1] / 2;
1086 } else {
1087 spotROI[2] = iyCenter - spotROISize[1] / 2;
1088 }
1089
1090 spotROI[3] = spotROI[2] + spotROISize[1] - 1;
1091 if (spotROI[2] < iy0) {
1092 spotROI[2] = iy0;
1093 spotROI[3] = iy0 + spotROISize[1] - 1;
1094 } else if (spotROI[3] > iy1) {
1095 spotROI[3] = iy1;
1096 spotROI[2] = iy1 - spotROISize[1] + 1;
1097 }
1098 if (spotROI[2] < iy0 || spotROI[3] > iy1)
1099 SDDS_Bomb("spot ROI is larger than ROI for y"); /* shouldn't ever see this */
1100
1101 /* perform background removal */
1102 analysisResults->saturationCount = 0;
1103 analysisResults->integratedSpotIntensity = 0;
1104 analysisResults->backgroundKilledNegative = 0;
1105 analysisResults->backgroundKilledPositive = 0;
1106 saturationLevel -= (long)background; // Cast to long to avoid implicit conversion warning
1107 for (ix = spotROI[0]; ix <= spotROI[1]; ix++) {
1108 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1109 if ((value = image[iy][ix] - background) <= 0) {
1110 analysisResults->backgroundKilledNegative += 1;
1111 value = 0;
1112 }
1113 image[iy][ix] = value;
1114 }
1115 }
1116 if (flags & SYMMETRIC_BGREMOVAL) {
1117 for (ix = spotROI[0]; ix <= spotROI[1]; ix++) {
1118 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1119 long ox, oy, goodPixels;
1120 long ox0, ox1, oy0, oy1;
1121 if (image[iy][ix] > 0 && image[iy][ix] <= backgroundHalfWidth) {
1122 /* if no more than one of the pixels around this pixel are >backgroundHalfWidth, then
1123 * set this pixel to zero
1124 */
1125 ox0 = (ix == spotROI[0] ? 0 : -1);
1126 ox1 = (ix == spotROI[1] ? 0 : 1);
1127 oy0 = (iy == spotROI[2] ? 0 : -1);
1128 oy1 = (iy == spotROI[3] ? 0 : 1);
1129 /* the pixel itself is counted but guaranteed to be subtracted off, too */
1130 /* this avoids an extra test in the loop */
1131 goodPixels = (ox1 - ox0 + 1) * (oy1 - oy0 + 1);
1132 for (ox = ox0; ox <= ox1; ox++) {
1133 for (oy = oy0; oy <= oy1; oy++) {
1134 if (image[iy + oy][ix + ox] <= backgroundHalfWidth)
1135 goodPixels--;
1136 }
1137 }
1138 if (goodPixels <= 1) {
1139 analysisResults->backgroundKilledPositive += 1;
1140 image[iy][ix] = 0;
1141 }
1142 }
1143 }
1144 }
1145 }
1146 if (flags & REMOVE_LONERS) {
1147 while (lonerPasses-- > 0) {
1148 for (ix = spotROI[0]; ix <= spotROI[1]; ix++) {
1149 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1150 long ox, oy, company;
1151 long ox0, ox1, oy0, oy1;
1152 if (image[iy][ix] > 0) {
1153 /* if none of the pixels around this pixel are nonzero, then set it to zero.
1154 */
1155 ox0 = (ix == spotROI[0] ? 0 : -1);
1156 ox1 = (ix == spotROI[1] ? 0 : 1);
1157 oy0 = (iy == spotROI[2] ? 0 : -1);
1158 oy1 = (iy == spotROI[3] ? 0 : 1);
1159 /* the pixel itself is counted but guaranteed to be subtracted off, too */
1160 /* this avoids an extra test in the loop */
1161 company = (ox1 - ox0 + 1) * (oy1 - oy0 + 1);
1162 for (ox = ox0; ox <= ox1; ox++) {
1163 for (oy = oy0; oy <= oy1; oy++) {
1164 if (image[iy + oy][ix + ox] == 0)
1165 company--;
1166 }
1167 }
1168 if (company <= lonerThreshold) {
1169 analysisResults->backgroundKilledPositive += 1;
1170 image[iy][ix] = 0;
1171 }
1172 }
1173 }
1174 }
1175 }
1176 }
1177 if (flags & ANTIHALO_BGREMOVAL) {
1178 long tryCount; // Renamed from 'try' to avoid confusion with C++ keyword
1179 for (tryCount = 0; tryCount < 2; tryCount++) {
1180 for (ix = spotROI[0]; ix <= spotROI[1]; ix++) {
1181 for (iy = spotROI[2]; iy < spotROI[3]; iy++) {
1182 if (image[iy][ix] > backgroundHalfWidth || image[iy + 1][ix] > backgroundHalfWidth)
1183 break;
1184 if (image[iy][ix]) {
1185 image[iy][ix] = 0;
1186 analysisResults->backgroundKilledPositive += 1;
1187 }
1188 }
1189 if (iy != spotROI[3])
1190 for (iy = spotROI[3]; iy > spotROI[2]; iy--) {
1191 if (image[iy][ix] > backgroundHalfWidth || image[iy - 1][ix] > backgroundHalfWidth)
1192 break;
1193 if (image[iy][ix]) {
1194 image[iy][ix] = 0;
1195 analysisResults->backgroundKilledPositive += 1;
1196 }
1197 }
1198 }
1199 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1200 for (ix = spotROI[0]; ix < spotROI[1]; ix++) {
1201 if (image[iy][ix] > backgroundHalfWidth || image[iy][ix + 1] > backgroundHalfWidth)
1202 break;
1203 if (image[iy][ix]) {
1204 image[iy][ix] = 0;
1205 analysisResults->backgroundKilledPositive += 1;
1206 }
1207 }
1208 if (ix != spotROI[1])
1209 for (ix = spotROI[1]; ix > spotROI[0]; ix--) {
1210 if (image[iy][ix] > backgroundHalfWidth || image[iy][ix - 1] > backgroundHalfWidth)
1211 break;
1212 if (image[iy][ix]) {
1213 image[iy][ix] = 0;
1214 analysisResults->backgroundKilledPositive += 1;
1215 }
1216 }
1217 }
1218 }
1219 }
1220 if (flags & SINGLE_SPOT) {
1221 short **connected, changed;
1222 double maxVal;
1223 long ix_max = 0, iy_max = 0;
1224
1225 connected = tmalloc(sizeof(*connected) * nx);
1226 for (ix = 0; ix < nx; ix++)
1227 connected[ix] = tmalloc(sizeof(**connected) * ny);
1228 maxVal = -HUGE_VAL;
1229 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1230 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1231 connected[ix][iy] = 0;
1232 if (image[iy][ix] > maxVal) { // Changed from image[ix][iy] to image[iy][ix]
1233 ix_max = ix;
1234 iy_max = iy;
1235 maxVal = image[iy][ix];
1236 }
1237 }
1238 connected[ix_max][iy_max] = 1;
1239
1240 do {
1241 changed = 0;
1242 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1243 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1244 if (!image[iy][ix] || connected[ix][iy])
1245 continue;
1246 if (ix > spotROI[0] && connected[ix - 1][iy]) {
1247 changed += (connected[ix][iy] = 1);
1248 continue;
1249 }
1250 if (ix < spotROI[1] && connected[ix + 1][iy]) {
1251 changed += (connected[ix][iy] = 1);
1252 continue;
1253 }
1254 if (iy > spotROI[2] && connected[ix][iy - 1]) {
1255 changed += (connected[ix][iy] = 1);
1256 continue;
1257 }
1258 if (iy < spotROI[3] && connected[ix][iy + 1]) {
1259 changed += (connected[ix][iy] = 1);
1260 continue;
1261 }
1262 }
1263 for (ix = spotROI[1]; ix >= spotROI[0]; ix--)
1264 for (iy = spotROI[3]; iy >= spotROI[2]; iy--) {
1265 if (!image[iy][ix] || connected[ix][iy])
1266 continue;
1267 if (ix > spotROI[0] && connected[ix - 1][iy]) {
1268 changed += (connected[ix][iy] = 1);
1269 continue;
1270 }
1271 if (ix < spotROI[1] && connected[ix + 1][iy]) {
1272 changed += (connected[ix][iy] = 1);
1273 continue;
1274 }
1275 if (iy > spotROI[2] && connected[ix][iy - 1]) {
1276 changed += (connected[ix][iy] = 1);
1277 continue;
1278 }
1279 if (iy < spotROI[3] && connected[ix][iy + 1]) {
1280 changed += (connected[ix][iy] = 1);
1281 continue;
1282 }
1283 }
1284 } while (changed);
1285
1286 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1287 for (iy = spotROI[2]; iy <= spotROI[3]; iy++)
1288 if (!connected[ix][iy])
1289 image[iy][ix] = 0;
1290
1291 /* Free connected memory */
1292 for (ix = 0; ix < nx; ix++) {
1293 free(connected[ix]);
1294 }
1295 free(connected);
1296 }
1297
1298 /* check for saturation */
1299 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1300 for (iy = spotROI[2]; iy <= spotROI[3]; iy++)
1301 if (image[iy][ix] > saturationLevel)
1302 analysisResults->saturationCount += 1;
1303
1304 /* find the spot intensity and centroids */
1305 analysisResults->spotCentroid[0] = analysisResults->spotCentroid[1] = 0;
1306 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1307 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1308 analysisResults->integratedSpotIntensity += image[iy][ix];
1309 analysisResults->spotCentroid[0] += image[iy][ix] * ix;
1310 analysisResults->spotCentroid[1] += image[iy][ix] * iy;
1311 }
1312 if (analysisResults->integratedSpotIntensity)
1313 for (i = 0; i < 2; i++)
1314 analysisResults->spotCentroid[i] /= analysisResults->integratedSpotIntensity;
1315
1316 /* find the spot size in y using central lines around the peak */
1317 if (!(lineBuffer = calloc(ny, sizeof(*lineBuffer))))
1318 SDDS_Bomb("memory allocation failure");
1319 if (!(x = calloc(ny, sizeof(*x))))
1320 SDDS_Bomb("memory allocation failure");
1321 if (!(y = calloc(ny, sizeof(*y))))
1322 SDDS_Bomb("memory allocation failure");
1323 for (ix = ixCenter - sizeLines[0] / 2; ix <= ixCenter + sizeLines[0] / 2; ix++) {
1324 if (ix < ix0 || ix > ix1)
1325 continue;
1326 for (iy = spotROI[2]; iy <= spotROI[3]; iy++)
1327 lineBuffer[iy] += image[iy][ix];
1328 }
1329 for (iy = spotROI[2]; iy <= spotROI[3]; iy++)
1330 y[iy] = lineBuffer[iy];
1331
1332 DetermineBeamSizes(&analysisResults->spotSigma[1], &analysisResults->spotRange50[1], &analysisResults->spotRange65[1],
1333 &analysisResults->spotRange80[1], lineBuffer, spotROI[2], spotROI[3]);
1334 free(lineBuffer);
1335 /* find the spot size in x using central lines around the peak */
1336 if (!(lineBuffer = calloc(nx, sizeof(*lineBuffer))))
1337 SDDS_Bomb("memory allocation failure");
1338 for (iy = iyCenter - sizeLines[1] / 2; iy <= iyCenter + sizeLines[1] / 2; iy++) {
1339 if (iy < iy0 || iy > iy1)
1340 continue;
1341 for (ix = spotROI[0]; ix <= spotROI[1]; ix++)
1342 lineBuffer[ix] += image[iy][ix];
1343 }
1344 DetermineBeamSizes(&analysisResults->spotSigma[0], &analysisResults->spotRange50[0], &analysisResults->spotRange65[0],
1345 &analysisResults->spotRange80[0], lineBuffer, spotROI[0], spotROI[1]);
1346 free(lineBuffer);
1347 DetermineBeamParameters(image, spotROI, nx, ny, &analysisResults->S11, &analysisResults->S33, &analysisResults->rmsCrossProduct,
1348 &analysisResults->phi, &analysisResults->majorAxis, &analysisResults->minorAxis);
1349 /* put results in the structure for return */
1350 analysisResults->peakSpotIntensity = maxValue - background;
1351 analysisResults->spotCenter[0] = ixCenter;
1352 analysisResults->spotCenter[1] = iyCenter;
1353 analysisResults->backgroundLevel = background;
1354 analysisResults->ROI[0] = ix0;
1355 analysisResults->ROI[1] = ix1;
1356 analysisResults->ROI[2] = iy0;
1357 analysisResults->ROI[3] = iy1;
1358 for (i = 0; i < 4; i++)
1359 analysisResults->spotROI[i] = spotROI[i];
1360 if (x)
1361 free(x);
1362 if (y)
1363 free(y);
1364
1365 if (SDDSspim) {
1366 long i_row;
1367 if (!SDDS_StartPage(SDDSspim, (spotROI[1] - spotROI[0] + 1) * (spotROI[3] - spotROI[2] + 1)))
1368 SDDS_Bomb("Problem starting page for spot image output file.");
1369 if (!SDDS_SetParameters(SDDSspim, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, "nx", (short)(spotROI[1] - spotROI[0] + 1),
1370 "ny", (short)(spotROI[3] - spotROI[2] + 1), NULL))
1371 SDDS_Bomb("Problem setting parameter values for spot image output file.");
1372 i_row = 0;
1373 for (iy = spotROI[2]; iy <= spotROI[3]; iy++) {
1374 for (ix = spotROI[0]; ix <= spotROI[1]; ix++) {
1375 if (!SDDS_SetRowValues(SDDSspim, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, i_row++, "ix", (short)ix,
1376 "iy", (short)iy, "Image", image[iy][ix], NULL)) {
1377 SDDS_Bomb("Problem setting row values for spot image output file.");
1378 }
1379 }
1380 }
1381 if (!SDDS_WritePage(SDDSspim)) {
1382 SDDS_Bomb("Problem writing page for spot image output file.");
1383 }
1384 }
1385 return 1;
1386}
1387
1388void DetermineBeamSizes(double *sigma, double *Range50, double *Range65, double *Range80, double *lineBuffer, long i0, long i1) {
1389 double centroid, sum;
1390 long i, j;
1391 double pLevel[6] = {.10, .175, .25, .75, .825, .90};
1392 double pValue[6];
1393
1394 centroid = sum = 0;
1395 *sigma = *Range80 = *Range65 = *Range50 = 0;
1396 for (i = i0; i <= i1; i++) {
1397 sum += lineBuffer[i];
1398 centroid += lineBuffer[i] * i;
1399 }
1400 if (sum) {
1401 centroid = centroid / sum;
1402 for (i = i0; i <= i1; i++)
1403 *sigma += lineBuffer[i] * sqr(i - centroid);
1404 *sigma = sqrt(*sigma / sum);
1405
1406 /* integrate the intensity */
1407 for (i = i0 + 1; i <= i1; i++)
1408 lineBuffer[i] += lineBuffer[i - 1];
1409 if (lineBuffer[i1]) {
1410 for (i = i0; i <= i1; i++)
1411 lineBuffer[i] /= lineBuffer[i1];
1412 i = i0 + 1;
1413 for (j = 0; j < 6; j++) {
1414 pValue[j] = 0;
1415 while (i <= i1 && lineBuffer[i] < pLevel[j])
1416 i++;
1417 if (i > i1) { // Added check to prevent out-of-bounds access
1418 pValue[j] = i1;
1419 } else if (lineBuffer[i] == lineBuffer[i - 1])
1420 pValue[j] = i - 0.5;
1421 else
1422 pValue[j] = i - (lineBuffer[i] - pLevel[j]) / (lineBuffer[i] - lineBuffer[i - 1]);
1423 }
1424 *Range80 = pValue[5] - pValue[0];
1425 *Range65 = pValue[4] - pValue[1];
1426 *Range50 = pValue[3] - pValue[2];
1427 }
1428 }
1429}
1430
1431void BlankOutImageData(double **image, long nx, long ny, int32_t *region) {
1432 long ix, iy, count = 0;
1433 for (ix = region[0]; ix <= region[1]; ix++)
1434 for (iy = region[2]; iy <= region[3]; iy++, count++)
1435 image[iy][ix] = 0;
1436}
1437
1438void DetermineBeamParameters(double **image, long *spotROI, long nx, long ny, double *S11, double *S33,
1439 double *rmsCrossProduct, double *phi, double *majorAxis, double *minorAxis) {
1440 long i, j;
1441 long x1, x2, y1, y2;
1442 double imageArea = 0, x2Ave = 0, y2Ave = 0, xyAve = 0, dominator = 0, xcentroid = 0, ycentroid = 0;
1443
1444 x1 = spotROI[0];
1445 x2 = spotROI[1];
1446 y1 = spotROI[2];
1447 y2 = spotROI[3];
1448
1449 /*calcuate the area of image(x,y) in a square defined by (x1,y1), (x1,y2), (x2,y2) and (x2,y1) */
1450 for (i = y1; i <= y2; i++)
1451 for (j = x1; j <= x2; j++) {
1452 imageArea += image[i][j];
1453 xcentroid += image[i][j] * j;
1454 ycentroid += image[i][j] * i;
1455 }
1456 if (imageArea == 0) {
1457 *rmsCrossProduct = *majorAxis = *minorAxis = DBL_MAX;
1458 } else {
1459 xcentroid = xcentroid / imageArea;
1460 ycentroid = ycentroid / imageArea;
1461 for (i = y1; i <= y2; i++)
1462 for (j = x1; j <= x2; j++) {
1463 x2Ave += sqr(j - xcentroid) * image[i][j];
1464 y2Ave += sqr(i - ycentroid) * image[i][j];
1465 xyAve += (i - ycentroid) * (j - xcentroid) * image[i][j];
1466 }
1467 x2Ave = x2Ave / imageArea;
1468 y2Ave = y2Ave / imageArea;
1469 xyAve = xyAve / imageArea;
1470 dominator = x2Ave * y2Ave - xyAve * xyAve;
1471 *S11 = x2Ave;
1472 *S33 = y2Ave;
1473 *rmsCrossProduct = xyAve;
1474 *phi = 0.5 * atan2(2 * xyAve, x2Ave - y2Ave) / PI * 180;
1475 if ((x2Ave + y2Ave - sqrt(sqr(x2Ave - y2Ave) + 4 * sqr(xyAve))) != 0) {
1476 *majorAxis = sqrt(2 * dominator / (x2Ave + y2Ave - sqrt(sqr(x2Ave - y2Ave) + 4 * sqr(xyAve))));
1477 } else {
1478 *majorAxis = DBL_MAX;
1479 }
1480 if ((x2Ave + y2Ave + sqrt(sqr(x2Ave - y2Ave) + 4 * sqr(xyAve))) != 0) {
1481 *minorAxis = sqrt(2 * dominator / (x2Ave + y2Ave + sqrt(sqr(x2Ave - y2Ave) + 4 * sqr(xyAve))));
1482 } else {
1483 *minorAxis = DBL_MAX;
1484 }
1485 }
1486}
1487
1488void replaceWithNearNeighbors(double **image, long iy0, long iy1, long ix0, long ix1, long iyc, long ixc, long distance) {
1489 long ix, iy, nnn;
1490 double sum;
1491
1492 if ((iyc - distance) > iy0)
1493 iy0 = iyc - distance;
1494 if ((iyc + distance) < iy1)
1495 iy1 = iyc + distance;
1496
1497 if ((ixc - distance) > ix0)
1498 ix0 = ixc - distance;
1499 if ((ixc + distance) < ix1)
1500 ix1 = ixc + distance;
1501
1502 sum = 0;
1503 nnn = 0;
1504 for (iy = iy0; iy <= iy1; iy++)
1505 for (ix = ix0; ix <= ix1; ix++) {
1506 if (ix == ixc && iy == iyc)
1507 continue;
1508 sum += image[iy][ix];
1509 nnn++;
1510 }
1511
1512 if (nnn > 0)
1513 image[iyc][ixc] = sum / nnn;
1514}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_LengthenTable
int32_t SDDS_LengthenTable(SDDS_DATASET *SDDS_dataset, int64_t n_additional_rows)
Definition SDDS_dataprep.c:297
SDDS_SetRowValues
int32_t SDDS_SetRowValues(SDDS_DATASET *SDDS_dataset, int32_t mode, int64_t row,...)
Definition SDDS_dataprep.c:988
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:64
SDDS_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:369
SDDS_GetParameterAsDouble
double * SDDS_GetParameterAsDouble(SDDS_DATASET *SDDS_dataset, char *parameter_name, double *memory)
Retrieves the value of a specified parameter as a double from the current data table of an SDDS datas...
Definition SDDS_extract.c:2782
SDDS_GetColumnInLong
int32_t * SDDS_GetColumnInLong(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of 32-bit integers,...
Definition SDDS_extract.c:1157
SDDS_GetParameter
void * SDDS_GetParameter(SDDS_DATASET *SDDS_dataset, char *parameter_name, void *memory)
Retrieves the value of a specified parameter from the current data table of a data set.
Definition SDDS_extract.c:2481
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_InitializeOutput
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
Definition SDDS_output.c:601
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
SDDS_DefineColumn
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
Definition SDDS_output.c:1715
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
SDDS_DefineParameter
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
Definition SDDS_output.c:1472
SDDS_DefineColumnLikeParameter
int32_t SDDS_DefineColumnLikeParameter(SDDS_DATASET *target, SDDS_DATASET *source, char *name, char *newName)
Defines a column in the target dataset based on a parameter definition from the source dataset.
Definition SDDS_transfer.c:264
SDDS_SetError
void SDDS_SetError(char *error_text)
Records an error message in the SDDS error stack.
Definition SDDS_utils.c:379
SDDS_GetParameterNames
char ** SDDS_GetParameterNames(SDDS_DATASET *SDDS_dataset, int32_t *number)
Retrieves the names of all parameters in the SDDS dataset.
Definition SDDS_utils.c:2501
SDDS_GetColumnIndex
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
Definition SDDS_utils.c:1309
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_Malloc
void * SDDS_Malloc(size_t size)
Allocates memory of a specified size.
Definition SDDS_utils.c:639
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
SDDS_STRING
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
SDDS_SHORT
#define SDDS_SHORT
Identifier for the signed short integer data type.
Definition SDDStypes.h:73
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDSutils.h
Utility functions for SDDS dataset manipulation and string array operations.
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
bitsSet
long bitsSet(unsigned long data)
Counts the number of set bits (1s) in the given data.
Definition binary.c:52
delete_chars
char * delete_chars(char *s, char *t)
Removes all occurrences of characters found in string t from string s.
Definition delete_chars.c:28
index_min_max
int index_min_max(int64_t *imin, int64_t *imax, double *list, int64_t n)
Finds the indices of the minimum and maximum values in a list of doubles.
Definition findMinMax.c:116
make_histogram
long make_histogram(double *hist, long n_bins, double lo, double hi, double *data, int64_t n_pts, long new_start)
Compiles a histogram from data points.
Definition makeHistogram.c:33
match_string
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
Definition match_string.c:34
scanargsg
int scanargsg(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:163
processPipeOption
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:356
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390
scanItemList
long scanItemList(unsigned long *flags, char **item, long *items, unsigned long mode,...)
Scans a list of items and assigns values based on provided keywords and types.
Definition scanitemlist.c:41
despikeData
long despikeData(double *data, long rows, long neighbors, long passes, long averageOf, double threshold, long countLimit)
Remove spikes from a data array by comparing each point to its neighbors.
Definition smooth.c:86
DESPIKE_SETTINGS
Definition sddsspotanalysis.c:185
HOTPIXEL_SETTINGS
Definition sddsspotanalysis.c:176
IMAGE_ANALYSIS
Definition sddsspotanalysis.c:167
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0