SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddsimageprofiles.c
Go to the documentation of this file.
1/**
2 * @file sddsimageprofiles.c
3 * @brief Analyze images stored as horizontal lines (one per column) based on the ideas of B-X Yang.
4 *
5 * @details
6 * This program processes image data stored in SDDS (Self Describing Data Set) format,
7 * allowing for various profile analyses such as center line, integrated, averaged, or peak profiles.
8 * It can handle background subtraction and define specific areas of interest within the image.
9 *
10 * @section Usage
11 * ```
12 * sddsimageprofiles [<inputfile>] [<outputfile>]
13 * [-pipe=[input][,output]]
14 * -columnPrefix=<prefix>
15 * [-profileType={x|y}]
16 * [-method={centerLine|integrated|averaged|peak}]
17 * [-background=<filename>|auto[,halfwidth=<value>][,keepNegative]]
18 * [-areaOfInterest=<rowStart>,<rowEnd>,<columnStart>,<columnEnd>]
19 * ```
20 *
21 * @section Options
22 * | Required | Description |
23 * |----------|-----------------------------------------------------------------------------|
24 * | `-columnPrefix` | Set the column prefix. |
25 *
26 * | Option | Description |
27 * |-------------------|--------------------------------------------------------------------------------------|
28 * | `-pipe` | Specify input and/or output via pipe. |
29 * | `-profileType` | Choose profile type: `x` or `y`. |
30 * | `-method` | Select the method for profile analysis: `centerLine`, `integrated`, `averaged`, `peak`.|
31 * | `-background` | Specify a background image file or estimate background from the image. |
32 * | `-areaOfInterest` | Define the area of interest within the image. |
33 *
34 * @copyright
35 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
36 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
37 *
38 * @license
39 * This file is distributed under the terms of the Software License Agreement
40 * found in the file LICENSE included with this distribution.
41 *
42 * @authors
43 * - R. Soliday
44 * - Xuesong Jiao
45 */
46
47#include "mdb.h"
48#include "SDDS.h"
49#include "scan.h"
50#include <float.h>
51
52#define BACKGROUND_AUTO 0x0001UL
53#define BACKGROUND_KEEP_NEGATIVE 0x0002UL
54#define AUTO_BACKGROUND_BINS 256
55
56/* Enumeration for option types */
57enum option_type {
58 SET_PIPE,
59 SET_PROFILETYPE,
60 SET_COLPREFIX,
61 SET_METHOD,
62 SET_AREAOFINTEREST,
63 SET_BACKGROUND,
64 N_OPTIONS,
65};
66
67char *option[N_OPTIONS] = {
68 "pipe", "profileType", "columnPrefix",
69 "method", "areaOfInterest", "background"};
70
71char *USAGE =
72 "sddsimageprofiles [<inputfile>] [<outputfile>]\n"
73 " [-pipe=[input][,output]]\n"
74 " -columnPrefix=<prefix>\n"
75 " [-profileType={x|y}]\n"
76 " [-method={centerLine|integrated|averaged|peak}]\n"
77 " [-background=<filename>|auto[,halfwidth=<value>][,keepNegative]]\n"
78 " [-areaOfInterest=<rowStart>,<rowEnd>,<columnStart>,<columnEnd>]\n"
79 "Options:\n"
80 " -pipe=[input][,output] Specify input and/or output via pipe.\n"
81 " -columnPrefix=<prefix> Set the column prefix.\n"
82 " -profileType={x|y} Choose profile type: 'x' or 'y'.\n"
83 " -method={centerLine|integrated|averaged|peak} Select the method for profile analysis.\n"
84 " -background=<filename> Specify a background image file.\n"
85 " -background=auto[,halfwidth=<value>][,keepNegative] Estimate background from the image, subtract it\n"
86 " before profiling, and clip negatives unless\n"
87 " keepNegative is given.\n"
88 " -areaOfInterest=<rowStart>,<rowEnd>,<columnStart>,<columnEnd> Define the area of interest.\n\n"
89 "Program by Robert Soliday. (\"" __DATE__ " " __TIME__ "\", SVN revision: " SVN_VERSION ")\n\n"
90 "-method:\n"
91 " If this option is not specified, it is a real profile.\n"
92 " If centerLine is specified, it will find the row with the\n"
93 " greatest integrated profile and display that line only.\n"
94 " If integrated is specified, it will sum all the profiles\n"
95 " together. If averaged is specified, it will divide the sum\n"
96 " of all the profiles by the number of profiles. If peak is\n"
97 " specified, it will find the peak point and display the profile\n"
98 " for that row.\n";
99
100typedef struct {
101 short *shortData;
102 unsigned short *ushortData;
103 int32_t *longData;
104 uint32_t *ulongData;
105 int64_t *long64Data;
106 uint64_t *ulong64Data;
107 float *floatData;
108 double *doubleData;
109} IMAGE_DATA;
110
111int xImageProfile(IMAGE_DATA *data, int32_t *type, int64_t rows, SDDS_DATASET *SDDS_dataset,
112 long method, int64_t x1, int64_t x2, long y1, long y2, long *colIndex, double *colIndex2);
113
114int yImageProfile(IMAGE_DATA *data, int32_t *type, int64_t rows, SDDS_DATASET *SDDS_dataset,
115 long method, int64_t x1, int64_t x2, long y1, long y2, long *colIndex, double *colIndex2);
116
117int64_t xPeakLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2);
118
119long yPeakLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2);
120
121int64_t xCenterLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2);
122
123long yCenterLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2);
124
125int64_t GetData(SDDS_DATASET *SDDS_orig, char *input, IMAGE_DATA **data, int32_t **type, long **colIndex,
126 double **colIndex2, char *colPrefix, long *validColumns);
127
128double getImageValue(const IMAGE_DATA *data, int32_t type, int64_t row);
129
130void convertImageDataToDouble(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns);
131
132double estimateBackgroundLevel(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns,
133 long backgroundHalfWidth);
134
135void subtractBackgroundLevel(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns,
136 double backgroundLevel, long clipNegative);
137
138int main(int argc, char **argv) {
139 SDDS_DATASET SDDS_dataset, SDDS_orig, SDDS_bg;
140 long i_arg;
141 SCANNED_ARG *s_arg;
142 IMAGE_DATA *data = NULL, *bg_data = NULL;
143 char *input = NULL, *output = NULL, *colPrefix = NULL, *background = NULL;
144 long profileType = 1, noWarnings = 0, tmpfile_used = 0, method = 0;
145 unsigned long pipeFlags = 0;
146 unsigned long backgroundFlags = 0;
147 int64_t rows, bg_rows, j;
148 long i, validColumns = 0, bg_validColumns = 0;
149 int32_t *type = NULL, *bg_type = NULL;
150 long backgroundHalfWidth = 1;
151 double backgroundLevel = 0;
152
153 int64_t rowStart = 1, rowEnd = 0;
154 long columnStart = 1, columnEnd = 0;
155
156 long *colIndex, *bg_colIndex;
157 double *colIndex2, *bg_colIndex2;
158
160 argc = scanargs(&s_arg, argc, argv);
161
162 if (argc < 3)
163 bomb(NULL, USAGE);
164
165 for (i_arg = 1; i_arg < argc; i_arg++) {
166 if (s_arg[i_arg].arg_type == OPTION) {
167 switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
168 case SET_PROFILETYPE:
169 if (s_arg[i_arg].n_items != 2)
170 SDDS_Bomb("invalid -profileType syntax");
171 if (strcmp("x", s_arg[i_arg].list[1]) == 0)
172 profileType = 1;
173 if (strcmp("y", s_arg[i_arg].list[1]) == 0)
174 profileType = 2;
175 break;
176 case SET_COLPREFIX:
177 if (s_arg[i_arg].n_items != 2)
178 SDDS_Bomb("invalid -columnPrefix syntax");
179 colPrefix = s_arg[i_arg].list[1];
180 break;
181 case SET_METHOD:
182 if (s_arg[i_arg].n_items != 2)
183 SDDS_Bomb("invalid -method syntax");
184 if ((strncasecmp("centralLine", s_arg[i_arg].list[1], strlen(s_arg[i_arg].list[1])) == 0) ||
185 (strncasecmp("centerLine", s_arg[i_arg].list[1], strlen(s_arg[i_arg].list[1])) == 0))
186 method = 1;
187 if (strncasecmp("integrated", s_arg[i_arg].list[1], strlen(s_arg[i_arg].list[1])) == 0)
188 method = 2;
189 if (strncasecmp("averaged", s_arg[i_arg].list[1], strlen(s_arg[i_arg].list[1])) == 0)
190 method = 3;
191 if (strncasecmp("peak", s_arg[i_arg].list[1], strlen(s_arg[i_arg].list[1])) == 0)
192 method = 4;
193 break;
194 case SET_AREAOFINTEREST:
195 if (s_arg[i_arg].n_items != 5)
196 SDDS_Bomb("invalid -areaOfInterest syntax");
197 if (sscanf(s_arg[i_arg].list[1], "%" SCNd64, &rowStart) != 1 || rowStart <= 0)
198 SDDS_Bomb("invalid -areaOfInterest syntax or value");
199 if (sscanf(s_arg[i_arg].list[2], "%" SCNd64, &rowEnd) != 1 || rowEnd <= 0)
200 SDDS_Bomb("invalid -areaOfInterest syntax or value");
201 if (sscanf(s_arg[i_arg].list[3], "%ld", &columnStart) != 1 || columnStart <= 0)
202 SDDS_Bomb("invalid -areaOfInterest syntax or value");
203 if (sscanf(s_arg[i_arg].list[4], "%ld", &columnEnd) != 1 || columnEnd <= 0)
204 SDDS_Bomb("invalid -areaOfInterest syntax or value");
205 break;
206 case SET_BACKGROUND:
207 if (s_arg[i_arg].n_items < 2)
208 SDDS_Bomb("invalid -background syntax");
209 if (strcasecmp(s_arg[i_arg].list[1], "auto") == 0) {
210 long items;
211 items = s_arg[i_arg].n_items - 1;
212 if (!scanItemList(&backgroundFlags, s_arg[i_arg].list + 1, &items, 0,
213 "auto", -1, NULL, 0, BACKGROUND_AUTO,
214 "halfwidth", SDDS_LONG, &backgroundHalfWidth, 1, 0,
215 "keepNegative", -1, NULL, 0, BACKGROUND_KEEP_NEGATIVE, NULL) ||
216 !(backgroundFlags & BACKGROUND_AUTO) ||
217 backgroundHalfWidth < 0)
218 SDDS_Bomb("invalid -background syntax/values");
219 } else {
220 if (s_arg[i_arg].n_items != 2)
221 SDDS_Bomb("invalid -background syntax");
222 background = s_arg[i_arg].list[1];
223 }
224 break;
225 case SET_PIPE:
226 if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
227 SDDS_Bomb("invalid -pipe syntax");
228 break;
229 default:
230 fprintf(stderr, "error: unknown switch: %s\n", s_arg[i_arg].list[0]);
231 exit(EXIT_FAILURE);
232 break;
233 }
234 } else {
235 if (input == NULL)
236 input = s_arg[i_arg].list[0];
237 else if (output == NULL)
238 output = s_arg[i_arg].list[0];
239 else
240 SDDS_Bomb("too many filenames");
241 }
242 }
243
244 if (colPrefix == NULL) {
245 fprintf(stderr, "error: missing columnPrefix\n");
246 exit(EXIT_FAILURE);
247 }
248
249 processFilenames("sddsimageprofiles", &input, &output, pipeFlags, noWarnings, &tmpfile_used);
250
251 /* Read in the image file */
252 rows = GetData(&SDDS_orig, input, &data, &type, &colIndex, &colIndex2, colPrefix, &validColumns);
253
254 if (rows < 0) {
255 fprintf(stderr, "error: no rows in image file\n");
256 exit(EXIT_FAILURE);
257 }
258
259 if (background != NULL) {
260 /* Read in the background image file */
261 bg_rows = GetData(&SDDS_bg, background, &bg_data, &bg_type, &bg_colIndex, &bg_colIndex2, colPrefix, &bg_validColumns);
262 if (rows != bg_rows) {
263 fprintf(stderr, "error: background has a different number of rows\n");
264 exit(EXIT_FAILURE);
265 }
266 if (validColumns != bg_validColumns) {
267 fprintf(stderr, "error: background has a different number of columns\n");
268 exit(EXIT_FAILURE);
269 }
270 /* Subtract the background from the image file */
271 for (i = 0; i < validColumns; i++) {
272 if (type[colIndex[i]] != bg_type[bg_colIndex[i]]) {
273 fprintf(stderr, "error: column types don't match with background image\n");
274 exit(EXIT_FAILURE);
275 }
276 if (colIndex2[i] != bg_colIndex2[i]) {
277 fprintf(stderr, "error: image rows don't match with background image\n");
278 exit(EXIT_FAILURE);
279 }
280 switch (type[colIndex[i]]) {
281 case SDDS_SHORT:
282 for (j = 0; j < rows; j++)
283 data[colIndex[i]].shortData[j] -= bg_data[bg_colIndex[i]].shortData[j];
284 break;
285 case SDDS_USHORT:
286 for (j = 0; j < rows; j++)
287 data[colIndex[i]].ushortData[j] -= bg_data[bg_colIndex[i]].ushortData[j];
288 break;
289 case SDDS_LONG:
290 for (j = 0; j < rows; j++)
291 data[colIndex[i]].longData[j] -= bg_data[bg_colIndex[i]].longData[j];
292 break;
293 case SDDS_ULONG:
294 for (j = 0; j < rows; j++)
295 data[colIndex[i]].ulongData[j] -= bg_data[bg_colIndex[i]].ulongData[j];
296 break;
297 case SDDS_LONG64:
298 for (j = 0; j < rows; j++)
299 data[colIndex[i]].long64Data[j] -= bg_data[bg_colIndex[i]].long64Data[j];
300 break;
301 case SDDS_ULONG64:
302 for (j = 0; j < rows; j++)
303 data[colIndex[i]].ulong64Data[j] -= bg_data[bg_colIndex[i]].ulong64Data[j];
304 break;
305 case SDDS_FLOAT:
306 for (j = 0; j < rows; j++)
307 data[colIndex[i]].floatData[j] -= bg_data[bg_colIndex[i]].floatData[j];
308 break;
309 case SDDS_DOUBLE:
310 for (j = 0; j < rows; j++)
311 data[colIndex[i]].doubleData[j] -= bg_data[bg_colIndex[i]].doubleData[j];
312 break;
313 default:
314 continue;
315 }
316 }
317 }
318
319 if (backgroundFlags & BACKGROUND_AUTO) {
320 convertImageDataToDouble(data, type, rows, colIndex, validColumns);
321 backgroundLevel = estimateBackgroundLevel(data, type, rows, colIndex, validColumns, backgroundHalfWidth);
322 if (!(pipeFlags & USE_STDOUT))
323 printf("Background level subtracted: %.15g\n", backgroundLevel);
324 subtractBackgroundLevel(data, type, rows, colIndex, validColumns, backgroundLevel,
325 !(backgroundFlags & BACKGROUND_KEEP_NEGATIVE));
326 }
327
328 /* Initialize the output file and define the columns */
329 if (!SDDS_InitializeOutput(&SDDS_dataset, SDDS_ASCII, 1, NULL, NULL, output)) {
330 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
331 exit(EXIT_FAILURE);
332 }
333
334 if (profileType == 1) {
335 if (SDDS_DefineParameter(&SDDS_dataset, "Zone", NULL, NULL, NULL, NULL, SDDS_STRING, 0) == -1) {
336 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
337 exit(EXIT_FAILURE);
338 }
339 if (SDDS_DefineColumn(&SDDS_dataset, "x", NULL, NULL, NULL, NULL, SDDS_LONG64, 0) == -1) {
340 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
341 exit(EXIT_FAILURE);
342 }
343 if (SDDS_DefineColumn(&SDDS_dataset, "y", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) == -1) {
344 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
345 exit(EXIT_FAILURE);
346 }
347 } else {
348 if (SDDS_DefineParameter(&SDDS_dataset, "Zone", NULL, NULL, NULL, NULL, SDDS_STRING, 0) == -1) {
349 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
350 exit(EXIT_FAILURE);
351 }
352 if (SDDS_DefineColumn(&SDDS_dataset, "x", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) == -1) {
353 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
354 exit(EXIT_FAILURE);
355 }
356 if (SDDS_DefineColumn(&SDDS_dataset, "y", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) == -1) {
357 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
358 exit(EXIT_FAILURE);
359 }
360 }
361 if (SDDS_WriteLayout(&SDDS_dataset) == 0) {
362 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
363 exit(EXIT_FAILURE);
364 }
365
366 if ((rowEnd > rows) || (rowEnd < rowStart))
367 rowEnd = rows;
368
369 if ((columnEnd > validColumns) || (columnEnd < columnStart))
370 columnEnd = validColumns;
371
372 if (profileType == 1) {
373 xImageProfile(data, type, rows, &SDDS_dataset, method, rowStart - 1, rowEnd, columnStart - 1, columnEnd, colIndex, colIndex2);
374 } else {
375 yImageProfile(data, type, rows, &SDDS_dataset, method, rowStart - 1, rowEnd, columnStart - 1, columnEnd, colIndex, colIndex2);
376 }
377
378 /* Close the output file */
379 if (SDDS_Terminate(&SDDS_dataset) != 1) {
380 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
381 exit(EXIT_FAILURE);
382 }
383
384 for (i = 0; i < validColumns; i++) {
385 switch (type[colIndex[i]]) {
386 case SDDS_SHORT:
387 free(data[colIndex[i]].shortData);
388 break;
389 case SDDS_USHORT:
390 free(data[colIndex[i]].ushortData);
391 break;
392 case SDDS_LONG:
393 free(data[colIndex[i]].longData);
394 break;
395 case SDDS_ULONG:
396 free(data[colIndex[i]].ulongData);
397 break;
398 case SDDS_LONG64:
399 free(data[colIndex[i]].long64Data);
400 break;
401 case SDDS_ULONG64:
402 free(data[colIndex[i]].ulong64Data);
403 break;
404 case SDDS_FLOAT:
405 free(data[colIndex[i]].floatData);
406 break;
407 case SDDS_DOUBLE:
408 free(data[colIndex[i]].doubleData);
409 break;
410 default:
411 continue;
412 }
413 }
414 free(data);
415 free(type);
416 if (background != NULL) {
417 for (i = 0; i < validColumns; i++) {
418 switch (bg_type[bg_colIndex[i]]) {
419 case SDDS_SHORT:
420 free(bg_data[bg_colIndex[i]].shortData);
421 break;
422 case SDDS_USHORT:
423 free(bg_data[bg_colIndex[i]].ushortData);
424 break;
425 case SDDS_LONG:
426 free(bg_data[bg_colIndex[i]].longData);
427 break;
428 case SDDS_ULONG:
429 free(bg_data[bg_colIndex[i]].ulongData);
430 break;
431 case SDDS_LONG64:
432 free(bg_data[bg_colIndex[i]].long64Data);
433 break;
434 case SDDS_ULONG64:
435 free(bg_data[bg_colIndex[i]].ulong64Data);
436 break;
437 case SDDS_FLOAT:
438 free(bg_data[bg_colIndex[i]].floatData);
439 break;
440 case SDDS_DOUBLE:
441 free(bg_data[bg_colIndex[i]].doubleData);
442 break;
443 default:
444 continue;
445 }
446 }
447 free(bg_data);
448 free(bg_type);
449 }
450
451 return EXIT_SUCCESS;
452}
453
454double getImageValue(const IMAGE_DATA *data, int32_t type, int64_t row) {
455 switch (type) {
456 case SDDS_SHORT:
457 return data->shortData[row];
458 case SDDS_USHORT:
459 return data->ushortData[row];
460 case SDDS_LONG:
461 return data->longData[row];
462 case SDDS_ULONG:
463 return data->ulongData[row];
464 case SDDS_LONG64:
465 return (double)data->long64Data[row];
466 case SDDS_ULONG64:
467 return (double)data->ulong64Data[row];
468 case SDDS_FLOAT:
469 return data->floatData[row];
470 case SDDS_DOUBLE:
471 return data->doubleData[row];
472 default:
473 return 0;
474 }
475}
476
477void convertImageDataToDouble(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns) {
478 double *doubleData;
479 int64_t row;
480 long i, column;
481
482 for (i = 0; i < validColumns; i++) {
483 column = colIndex[i];
484 if (type[column] == SDDS_DOUBLE)
485 continue;
486 if (!(doubleData = malloc(sizeof(*doubleData) * rows)))
487 SDDS_Bomb("memory allocation failure");
488 for (row = 0; row < rows; row++)
489 doubleData[row] = getImageValue(data + column, type[column], row);
490 switch (type[column]) {
491 case SDDS_SHORT:
492 free(data[column].shortData);
493 break;
494 case SDDS_USHORT:
495 free(data[column].ushortData);
496 break;
497 case SDDS_LONG:
498 free(data[column].longData);
499 break;
500 case SDDS_ULONG:
501 free(data[column].ulongData);
502 break;
503 case SDDS_LONG64:
504 free(data[column].long64Data);
505 break;
506 case SDDS_ULONG64:
507 free(data[column].ulong64Data);
508 break;
509 case SDDS_FLOAT:
510 free(data[column].floatData);
511 break;
512 default:
513 break;
514 }
515 data[column].doubleData = doubleData;
516 type[column] = SDDS_DOUBLE;
517 }
518}
519
520double estimateBackgroundLevel(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns,
521 long backgroundHalfWidth) {
522 double minValue, maxValue, binWidth, backgroundLevel, value;
523 double sum = 0;
524 int64_t row, pixelCount = 0;
525 long i, bin, modeBin = 0;
526 int64_t *histogram;
527 int64_t backgroundPixels = 0, maxCount = 0;
528
529 minValue = DBL_MAX;
530 maxValue = -DBL_MAX;
531 for (i = 0; i < validColumns; i++) {
532 for (row = 0; row < rows; row++) {
533 value = getImageValue(data + colIndex[i], type[colIndex[i]], row);
534 if (value < minValue)
535 minValue = value;
536 if (value > maxValue)
537 maxValue = value;
538 pixelCount++;
539 }
540 }
541 if (!pixelCount)
542 return 0;
543 if (minValue == maxValue)
544 return minValue;
545
546 if (!(histogram = calloc(AUTO_BACKGROUND_BINS, sizeof(*histogram))))
547 SDDS_Bomb("memory allocation failure");
548 binWidth = (maxValue - minValue) / AUTO_BACKGROUND_BINS;
549 if (binWidth <= 0) {
550 free(histogram);
551 return minValue;
552 }
553
554 for (i = 0; i < validColumns; i++) {
555 for (row = 0; row < rows; row++) {
556 value = getImageValue(data + colIndex[i], type[colIndex[i]], row);
557 bin = (long)((value - minValue) / binWidth);
558 if (bin < 0)
559 bin = 0;
560 else if (bin >= AUTO_BACKGROUND_BINS)
561 bin = AUTO_BACKGROUND_BINS - 1;
562 histogram[bin]++;
563 }
564 }
565
566 for (bin = 0; bin < AUTO_BACKGROUND_BINS; bin++) {
567 if (histogram[bin] > maxCount) {
568 maxCount = histogram[bin];
569 modeBin = bin;
570 }
571 }
572
573 for (i = 0; i < validColumns; i++) {
574 for (row = 0; row < rows; row++) {
575 value = getImageValue(data + colIndex[i], type[colIndex[i]], row);
576 bin = (long)((value - minValue) / binWidth);
577 if (bin < 0)
578 bin = 0;
579 else if (bin >= AUTO_BACKGROUND_BINS)
580 bin = AUTO_BACKGROUND_BINS - 1;
581 if (bin >= modeBin - backgroundHalfWidth && bin <= modeBin + backgroundHalfWidth) {
582 sum += value;
583 backgroundPixels++;
584 }
585 }
586 }
587 free(histogram);
588
589 if (backgroundPixels)
590 backgroundLevel = sum / backgroundPixels;
591 else
592 backgroundLevel = minValue + (modeBin + 0.5) * binWidth;
593
594 return backgroundLevel;
595}
596
597void subtractBackgroundLevel(IMAGE_DATA *data, int32_t *type, int64_t rows, long *colIndex, long validColumns,
598 double backgroundLevel, long clipNegative) {
599 double value;
600 int64_t row;
601 long i, column;
602
603 for (i = 0; i < validColumns; i++) {
604 column = colIndex[i];
605 if (type[column] != SDDS_DOUBLE)
606 continue;
607 for (row = 0; row < rows; row++) {
608 value = data[column].doubleData[row] - backgroundLevel;
609 if (clipNegative && value < 0)
610 value = 0;
611 data[column].doubleData[row] = value;
612 }
613 }
614}
615
616int xImageProfile(IMAGE_DATA *data, int32_t *type, int64_t rows, SDDS_DATASET *SDDS_dataset, long method,
617 int64_t x1, int64_t x2, long y1, long y2, long *colIndex, double *colIndex2) {
618 int64_t i, k = 0;
619 int j;
620 double val = 0;
621 long center;
622 int64_t *index;
623 double *values;
624 char value[100];
625 index = malloc(sizeof(int64_t) * rows);
626 values = malloc(sizeof(double) * rows);
627
628 if (method == 0) { /* Highest point */
629 for (i = x1; i < x2; i++) {
630 switch (type[colIndex[y1]]) {
631 case SDDS_SHORT:
632 val = data[colIndex[y1]].shortData[i];
633 break;
634 case SDDS_USHORT:
635 val = data[colIndex[y1]].ushortData[i];
636 break;
637 case SDDS_LONG:
638 val = data[colIndex[y1]].longData[i];
639 break;
640 case SDDS_ULONG:
641 val = data[colIndex[y1]].ulongData[i];
642 break;
643 case SDDS_LONG64:
644 val = data[colIndex[y1]].long64Data[i];
645 break;
646 case SDDS_ULONG64:
647 val = data[colIndex[y1]].ulong64Data[i];
648 break;
649 case SDDS_FLOAT:
650 val = data[colIndex[y1]].floatData[i];
651 break;
652 case SDDS_DOUBLE:
653 val = data[colIndex[y1]].doubleData[i];
654 break;
655 }
656 for (j = y1 + 1; j < y2; j++) {
657 switch (type[colIndex[j]]) {
658 case SDDS_SHORT:
659 if (val < data[colIndex[j]].shortData[i])
660 val = data[colIndex[j]].shortData[i];
661 break;
662 case SDDS_USHORT:
663 if (val < data[colIndex[j]].ushortData[i])
664 val = data[colIndex[j]].ushortData[i];
665 break;
666 case SDDS_LONG:
667 if (val < data[colIndex[j]].longData[i])
668 val = data[colIndex[j]].longData[i];
669 break;
670 case SDDS_ULONG:
671 if (val < data[colIndex[j]].ulongData[i])
672 val = data[colIndex[j]].ulongData[i];
673 break;
674 case SDDS_LONG64:
675 if (val < data[colIndex[j]].long64Data[i])
676 val = data[colIndex[j]].long64Data[i];
677 break;
678 case SDDS_ULONG64:
679 if (val < data[colIndex[j]].ulong64Data[i])
680 val = data[colIndex[j]].ulong64Data[i];
681 break;
682 case SDDS_FLOAT:
683 if (val < data[colIndex[j]].floatData[i])
684 val = data[colIndex[j]].floatData[i];
685 break;
686 case SDDS_DOUBLE:
687 if (val < data[colIndex[j]].doubleData[i])
688 val = data[colIndex[j]].doubleData[i];
689 break;
690 }
691 }
692 index[k] = i + 1;
693 values[k] = val;
694 k++;
695 }
696 } else if ((method == 1) || (method == 4)) { /* Center line or peak */
697 if (method == 4) {
698 center = yPeakLine(data, type, colIndex, x1, x2, y1, y2);
699 } else {
700 center = yCenterLine(data, type, colIndex, x1, x2, y1, y2);
701 }
702 for (i = x1; i < x2; i++) {
703 switch (type[colIndex[center]]) {
704 case SDDS_SHORT:
705 val = data[colIndex[center]].shortData[i];
706 break;
707 case SDDS_USHORT:
708 val = data[colIndex[center]].ushortData[i];
709 break;
710 case SDDS_LONG:
711 val = data[colIndex[center]].longData[i];
712 break;
713 case SDDS_ULONG:
714 val = data[colIndex[center]].ulongData[i];
715 break;
716 case SDDS_LONG64:
717 val = data[colIndex[center]].long64Data[i];
718 break;
719 case SDDS_ULONG64:
720 val = data[colIndex[center]].ulong64Data[i];
721 break;
722 case SDDS_FLOAT:
723 val = data[colIndex[center]].floatData[i];
724 break;
725 case SDDS_DOUBLE:
726 val = data[colIndex[center]].doubleData[i];
727 break;
728 }
729 index[k] = i + 1;
730 values[k] = val;
731 k++;
732 }
733 } else if ((method == 2) || (method == 3)) { /* Integrated or Averaged profile */
734 for (i = x1; i < x2; i++) {
735 val = 0;
736 for (j = y1; j < y2; j++) {
737 switch (type[colIndex[j]]) {
738 case SDDS_SHORT:
739 val += data[colIndex[j]].shortData[i];
740 break;
741 case SDDS_USHORT:
742 val += data[colIndex[j]].ushortData[i];
743 break;
744 case SDDS_LONG:
745 val += data[colIndex[j]].longData[i];
746 break;
747 case SDDS_ULONG:
748 val += data[colIndex[j]].ulongData[i];
749 break;
750 case SDDS_LONG64:
751 val += data[colIndex[j]].long64Data[i];
752 break;
753 case SDDS_ULONG64:
754 val += data[colIndex[j]].ulong64Data[i];
755 break;
756 case SDDS_FLOAT:
757 val += data[colIndex[j]].floatData[i];
758 break;
759 case SDDS_DOUBLE:
760 val += data[colIndex[j]].doubleData[i];
761 break;
762 }
763 }
764 index[k] = i + 1;
765 if (method == 2)
766 values[k] = val;
767 else
768 values[k] = val / (y2 - y1);
769 k++;
770 }
771 }
772
773 if (SDDS_StartPage(SDDS_dataset, x2 - x1) == 0) {
774 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
775 exit(EXIT_FAILURE);
776 }
777 sprintf(value, "(%" PRId64 ",%ld) x (%" PRId64 ",%ld)", x1 + 1, y1 + 1, x2, y2);
778 if (SDDS_SetParameters(SDDS_dataset, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "Zone", value, NULL) != 1) {
779 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
780 exit(EXIT_FAILURE);
781 }
782 if (SDDS_SetColumn(SDDS_dataset, SDDS_SET_BY_NAME, index, k, "x", NULL) != 1) {
783 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
784 exit(EXIT_FAILURE);
785 }
786 if (SDDS_SetColumn(SDDS_dataset, SDDS_SET_BY_NAME, values, k, "y", NULL) != 1) {
787 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
788 exit(EXIT_FAILURE);
789 }
790 if (SDDS_WritePage(SDDS_dataset) != 1) {
791 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
792 exit(EXIT_FAILURE);
793 }
794 free(index);
795 free(values);
796 return 0;
797}
798
799int yImageProfile(IMAGE_DATA *data, int32_t *type, int64_t rows, SDDS_DATASET *SDDS_dataset, long method,
800 int64_t x1, int64_t x2, long y1, long y2, long *colIndex, double *colIndex2) {
801 int i, k = 0;
802 int64_t j;
803 double val;
804 int64_t center;
805
806 double *index;
807 double *values;
808 char value[100];
809
810 index = malloc(sizeof(double) * (y2 - y1));
811 values = malloc(sizeof(double) * (y2 - y1));
812
813 if (method == 0) { /* Highest point */
814 for (i = y1; i < y2; i++) {
815 switch (type[colIndex[i]]) {
816 case SDDS_SHORT:
817 val = data[colIndex[i]].shortData[x1];
818 break;
819 case SDDS_USHORT:
820 val = data[colIndex[i]].ushortData[x1];
821 break;
822 case SDDS_LONG:
823 val = data[colIndex[i]].longData[x1];
824 break;
825 case SDDS_ULONG:
826 val = data[colIndex[i]].ulongData[x1];
827 break;
828 case SDDS_LONG64:
829 val = data[colIndex[i]].long64Data[x1];
830 break;
831 case SDDS_ULONG64:
832 val = data[colIndex[i]].ulong64Data[x1];
833 break;
834 case SDDS_FLOAT:
835 val = data[colIndex[i]].floatData[x1];
836 break;
837 case SDDS_DOUBLE:
838 val = data[colIndex[i]].doubleData[x1];
839 break;
840 default:
841 continue;
842 }
843 for (j = x1 + 1; j < x2; j++) {
844 switch (type[colIndex[i]]) {
845 case SDDS_SHORT:
846 if (val < data[colIndex[i]].shortData[j]) {
847 val = data[colIndex[i]].shortData[j];
848 }
849 break;
850 case SDDS_USHORT:
851 if (val < data[colIndex[i]].ushortData[j]) {
852 val = data[colIndex[i]].ushortData[j];
853 }
854 break;
855 case SDDS_LONG:
856 if (val < data[colIndex[i]].longData[j]) {
857 val = data[colIndex[i]].longData[j];
858 }
859 break;
860 case SDDS_ULONG:
861 if (val < data[colIndex[i]].ulongData[j]) {
862 val = data[colIndex[i]].ulongData[j];
863 }
864 break;
865 case SDDS_LONG64:
866 if (val < data[colIndex[i]].long64Data[j]) {
867 val = data[colIndex[i]].long64Data[j];
868 }
869 break;
870 case SDDS_ULONG64:
871 if (val < data[colIndex[i]].ulong64Data[j]) {
872 val = data[colIndex[i]].ulong64Data[j];
873 }
874 break;
875 case SDDS_FLOAT:
876 if (val < data[colIndex[i]].floatData[j]) {
877 val = data[colIndex[i]].floatData[j];
878 }
879 break;
880 case SDDS_DOUBLE:
881 if (val < data[colIndex[i]].doubleData[j]) {
882 val = data[colIndex[i]].doubleData[j];
883 }
884 break;
885 default:
886 continue;
887 }
888 }
889 index[k] = colIndex2[i];
890 values[k] = val;
891 k++;
892 }
893 } else if ((method == 1) || (method == 4)) { /* Center line or peak */
894 if (method == 4) {
895 center = xPeakLine(data, type, colIndex, x1, x2, y1, y2);
896 } else {
897 center = xCenterLine(data, type, colIndex, x1, x2, y1, y2);
898 }
899 for (i = y1; i < y2; i++) {
900 switch (type[colIndex[i]]) {
901 case SDDS_SHORT:
902 val = data[colIndex[i]].shortData[center];
903 break;
904 case SDDS_USHORT:
905 val = data[colIndex[i]].ushortData[center];
906 break;
907 case SDDS_LONG:
908 val = data[colIndex[i]].longData[center];
909 break;
910 case SDDS_ULONG:
911 val = data[colIndex[i]].ulongData[center];
912 break;
913 case SDDS_LONG64:
914 val = data[colIndex[i]].long64Data[center];
915 break;
916 case SDDS_ULONG64:
917 val = data[colIndex[i]].ulong64Data[center];
918 break;
919 case SDDS_FLOAT:
920 val = data[colIndex[i]].floatData[center];
921 break;
922 case SDDS_DOUBLE:
923 val = data[colIndex[i]].doubleData[center];
924 break;
925 default:
926 continue;
927 }
928 index[k] = colIndex2[i];
929 values[k] = val;
930 k++;
931 }
932 } else if ((method == 2) || (method == 3)) { /* Integrated or Averaged profile */
933 for (i = y1; i < y2; i++) {
934 val = 0;
935 switch (type[colIndex[i]]) {
936 case SDDS_SHORT:
937 for (j = x1; j < x2; j++)
938 val += data[colIndex[i]].shortData[j];
939 break;
940 case SDDS_USHORT:
941 for (j = x1; j < x2; j++)
942 val += data[colIndex[i]].ushortData[j];
943 break;
944 case SDDS_LONG:
945 for (j = x1; j < x2; j++)
946 val += data[colIndex[i]].longData[j];
947 break;
948 case SDDS_ULONG:
949 for (j = x1; j < x2; j++)
950 val += data[colIndex[i]].ulongData[j];
951 break;
952 case SDDS_LONG64:
953 for (j = x1; j < x2; j++)
954 val += data[colIndex[i]].long64Data[j];
955 break;
956 case SDDS_ULONG64:
957 for (j = x1; j < x2; j++)
958 val += data[colIndex[i]].ulong64Data[j];
959 break;
960 case SDDS_FLOAT:
961 for (j = x1; j < x2; j++)
962 val += data[colIndex[i]].floatData[j];
963 break;
964 case SDDS_DOUBLE:
965 for (j = x1; j < x2; j++)
966 val += data[colIndex[i]].doubleData[j];
967 break;
968 default:
969 continue;
970 }
971 index[k] = colIndex2[i];
972 if (method == 2)
973 values[k] = val;
974 else
975 values[k] = val / (x2 - x1);
976 k++;
977 }
978 }
979
980 if (SDDS_StartPage(SDDS_dataset, k) == 0) {
981 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
982 exit(EXIT_FAILURE);
983 }
984 sprintf(value, "(%" PRId64 ",%ld) x (%" PRId64 ",%ld)", x1 + 1, y1 + 1, x2, y2);
985 if (SDDS_SetParameters(SDDS_dataset, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "Zone", value, NULL) != 1) {
986 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
987 exit(EXIT_FAILURE);
988 }
989 if (SDDS_SetColumn(SDDS_dataset, SDDS_SET_BY_NAME, index, k, "y", NULL) != 1) {
990 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
991 exit(EXIT_FAILURE);
992 }
993 if (SDDS_SetColumn(SDDS_dataset, SDDS_SET_BY_NAME, values, k, "x", NULL) != 1) {
994 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
995 exit(EXIT_FAILURE);
996 }
997 if (SDDS_WritePage(SDDS_dataset) != 1) {
998 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
999 exit(EXIT_FAILURE);
1000 }
1001 free(index);
1002 free(values);
1003 return 0;
1004}
1005
1006int64_t xPeakLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2) {
1007 int64_t i, j;
1008 double maxValue = 0;
1009 int64_t index = 0;
1010
1011 index = x1;
1012 switch (type[colIndex[y1]]) {
1013 case SDDS_SHORT:
1014 maxValue = data[colIndex[y1]].shortData[x1];
1015 break;
1016 case SDDS_USHORT:
1017 maxValue = data[colIndex[y1]].ushortData[x1];
1018 break;
1019 case SDDS_LONG:
1020 maxValue = data[colIndex[y1]].longData[x1];
1021 break;
1022 case SDDS_ULONG:
1023 maxValue = data[colIndex[y1]].ulongData[x1];
1024 break;
1025 case SDDS_LONG64:
1026 maxValue = data[colIndex[y1]].long64Data[x1];
1027 break;
1028 case SDDS_ULONG64:
1029 maxValue = data[colIndex[y1]].ulong64Data[x1];
1030 break;
1031 case SDDS_FLOAT:
1032 maxValue = data[colIndex[y1]].floatData[x1];
1033 break;
1034 case SDDS_DOUBLE:
1035 maxValue = data[colIndex[y1]].doubleData[x1];
1036 break;
1037 }
1038 for (i = y1; i < y2; i++) {
1039 for (j = x1; j < x2; j++) {
1040 switch (type[colIndex[i]]) {
1041 case SDDS_SHORT:
1042 if (maxValue < data[colIndex[i]].shortData[j]) {
1043 maxValue = data[colIndex[i]].shortData[j];
1044 index = j;
1045 }
1046 break;
1047 case SDDS_USHORT:
1048 if (maxValue < data[colIndex[i]].ushortData[j]) {
1049 maxValue = data[colIndex[i]].ushortData[j];
1050 index = j;
1051 }
1052 break;
1053 case SDDS_LONG:
1054 if (maxValue < data[colIndex[i]].longData[j]) {
1055 maxValue = data[colIndex[i]].longData[j];
1056 index = j;
1057 }
1058 break;
1059 case SDDS_ULONG:
1060 if (maxValue < data[colIndex[i]].ulongData[j]) {
1061 maxValue = data[colIndex[i]].ulongData[j];
1062 index = j;
1063 }
1064 break;
1065 case SDDS_LONG64:
1066 if (maxValue < data[colIndex[i]].long64Data[j]) {
1067 maxValue = data[colIndex[i]].long64Data[j];
1068 index = j;
1069 }
1070 break;
1071 case SDDS_ULONG64:
1072 if (maxValue < data[colIndex[i]].ulong64Data[j]) {
1073 maxValue = data[colIndex[i]].ulong64Data[j];
1074 index = j;
1075 }
1076 break;
1077 case SDDS_FLOAT:
1078 if (maxValue < data[colIndex[i]].floatData[j]) {
1079 maxValue = data[colIndex[i]].floatData[j];
1080 index = j;
1081 }
1082 break;
1083 case SDDS_DOUBLE:
1084 if (maxValue < data[colIndex[i]].doubleData[j]) {
1085 maxValue = data[colIndex[i]].doubleData[j];
1086 index = j;
1087 }
1088 break;
1089 }
1090 }
1091 }
1092 return index;
1093}
1094
1095long yPeakLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2) {
1096 int i;
1097 int64_t j;
1098 double maxValue = 0;
1099 long index = 0;
1100
1101 index = y1;
1102 switch (type[colIndex[y1]]) {
1103 case SDDS_SHORT:
1104 maxValue = data[colIndex[y1]].shortData[x1];
1105 break;
1106 case SDDS_USHORT:
1107 maxValue = data[colIndex[y1]].ushortData[x1];
1108 break;
1109 case SDDS_LONG:
1110 maxValue = data[colIndex[y1]].longData[x1];
1111 break;
1112 case SDDS_ULONG:
1113 maxValue = data[colIndex[y1]].ulongData[x1];
1114 break;
1115 case SDDS_LONG64:
1116 maxValue = data[colIndex[y1]].long64Data[x1];
1117 break;
1118 case SDDS_ULONG64:
1119 maxValue = data[colIndex[y1]].ulong64Data[x1];
1120 break;
1121 case SDDS_FLOAT:
1122 maxValue = data[colIndex[y1]].floatData[x1];
1123 break;
1124 case SDDS_DOUBLE:
1125 maxValue = data[colIndex[y1]].doubleData[x1];
1126 break;
1127 }
1128 for (i = y1; i < y2; i++) {
1129 for (j = x1; j < x2; j++) {
1130 switch (type[colIndex[i]]) {
1131 case SDDS_SHORT:
1132 if (maxValue < data[colIndex[i]].shortData[j]) {
1133 maxValue = data[colIndex[i]].shortData[j];
1134 index = i;
1135 }
1136 break;
1137 case SDDS_USHORT:
1138 if (maxValue < data[colIndex[i]].ushortData[j]) {
1139 maxValue = data[colIndex[i]].ushortData[j];
1140 index = i;
1141 }
1142 break;
1143 case SDDS_LONG:
1144 if (maxValue < data[colIndex[i]].longData[j]) {
1145 maxValue = data[colIndex[i]].longData[j];
1146 index = i;
1147 }
1148 break;
1149 case SDDS_ULONG:
1150 if (maxValue < data[colIndex[i]].ulongData[j]) {
1151 maxValue = data[colIndex[i]].ulongData[j];
1152 index = i;
1153 }
1154 break;
1155 case SDDS_LONG64:
1156 if (maxValue < data[colIndex[i]].long64Data[j]) {
1157 maxValue = data[colIndex[i]].long64Data[j];
1158 index = i;
1159 }
1160 break;
1161 case SDDS_ULONG64:
1162 if (maxValue < data[colIndex[i]].ulong64Data[j]) {
1163 maxValue = data[colIndex[i]].ulong64Data[j];
1164 index = i;
1165 }
1166 break;
1167 case SDDS_FLOAT:
1168 if (maxValue < data[colIndex[i]].floatData[j]) {
1169 maxValue = data[colIndex[i]].floatData[j];
1170 index = i;
1171 }
1172 break;
1173 case SDDS_DOUBLE:
1174 if (maxValue < data[colIndex[i]].doubleData[j]) {
1175 maxValue = data[colIndex[i]].doubleData[j];
1176 index = i;
1177 }
1178 break;
1179 }
1180 }
1181 }
1182 return index;
1183}
1184
1185int64_t xCenterLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2) {
1186 int64_t i, j, start = 1;
1187 double val = 0, maxValue = 0;
1188 int64_t index = 0;
1189
1190 for (i = x1; i < x2; i++) {
1191 val = 0;
1192 for (j = y1; j < y2; j++) {
1193 switch (type[colIndex[j]]) {
1194 case SDDS_SHORT:
1195 val += data[colIndex[j]].shortData[i];
1196 break;
1197 case SDDS_USHORT:
1198 val += data[colIndex[j]].ushortData[i];
1199 break;
1200 case SDDS_LONG:
1201 val += data[colIndex[j]].longData[i];
1202 break;
1203 case SDDS_ULONG:
1204 val += data[colIndex[j]].ulongData[i];
1205 break;
1206 case SDDS_LONG64:
1207 val += data[colIndex[j]].long64Data[i];
1208 break;
1209 case SDDS_ULONG64:
1210 val += data[colIndex[j]].ulong64Data[i];
1211 break;
1212 case SDDS_FLOAT:
1213 val += data[colIndex[j]].floatData[i];
1214 break;
1215 case SDDS_DOUBLE:
1216 val += data[colIndex[j]].doubleData[i];
1217 break;
1218 }
1219 }
1220 if (start == 1) {
1221 index = i;
1222 maxValue = val;
1223 start = 0;
1224 } else {
1225 if (val > maxValue) {
1226 index = i;
1227 maxValue = val;
1228 }
1229 }
1230 }
1231 return index;
1232}
1233
1234long yCenterLine(IMAGE_DATA *data, int32_t *type, long *colIndex, int64_t x1, int64_t x2, long y1, long y2) {
1235 int i, start = 1;
1236 int64_t j;
1237 double val, maxValue = 0;
1238 long index = 0;
1239
1240 for (i = y1; i < y2; i++) {
1241 val = 0;
1242 for (j = x1; j < x2; j++)
1243 switch (type[colIndex[i]]) {
1244 case SDDS_SHORT:
1245 val += data[colIndex[i]].shortData[j];
1246 break;
1247 case SDDS_USHORT:
1248 val += data[colIndex[i]].ushortData[j];
1249 break;
1250 case SDDS_LONG:
1251 val += data[colIndex[i]].longData[j];
1252 break;
1253 case SDDS_ULONG:
1254 val += data[colIndex[i]].ulongData[j];
1255 break;
1256 case SDDS_LONG64:
1257 val += data[colIndex[i]].long64Data[j];
1258 break;
1259 case SDDS_ULONG64:
1260 val += data[colIndex[i]].ulong64Data[j];
1261 break;
1262 case SDDS_FLOAT:
1263 val += data[colIndex[i]].floatData[j];
1264 break;
1265 case SDDS_DOUBLE:
1266 val += data[colIndex[i]].doubleData[j];
1267 break;
1268 default:
1269 continue;
1270 }
1271 if (start == 1) {
1272 index = i;
1273 maxValue = val;
1274 start = 0;
1275 } else {
1276 if (val > maxValue) {
1277 index = i;
1278 maxValue = val;
1279 }
1280 }
1281 }
1282 return index;
1283}
1284
1285int64_t GetData(SDDS_DATASET *SDDS_orig, char *input, IMAGE_DATA **data, int32_t **type, long **colIndex,
1286 double **colIndex2, char *colPrefix, long *validColumns) {
1287 int64_t rows;
1288 long i, j, temp;
1289 double temp2;
1290 int32_t orig_column_names;
1291 char **orig_column_name;
1292
1293 /* Open image file */
1294 if (!SDDS_InitializeInput(SDDS_orig, input)) {
1295 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1296 exit(EXIT_FAILURE);
1297 }
1298
1299 /* Read column names */
1300 if (!(orig_column_name = SDDS_GetColumnNames(SDDS_orig, &orig_column_names))) {
1301 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1302 exit(EXIT_FAILURE);
1303 }
1304
1305 /* Allocate memory for image file data */
1306 *data = malloc(sizeof(IMAGE_DATA) * orig_column_names);
1307 *type = malloc(sizeof(int32_t) * orig_column_names);
1308
1309 /* Read page */
1310 if (SDDS_ReadPage(SDDS_orig) != 1) {
1311 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1312 exit(EXIT_FAILURE);
1313 }
1314
1315 /* Get number of rows */
1316 rows = SDDS_RowCount(SDDS_orig);
1317 *colIndex = malloc(sizeof(long) * orig_column_names);
1318 *colIndex2 = malloc(sizeof(double) * orig_column_names);
1319
1320 /* Read all numerical data from file */
1321 for (i = 0; i < orig_column_names; i++) {
1322 (*type)[i] = SDDS_GetNamedColumnType(SDDS_orig, orig_column_name[i]);
1323 if ((*type)[i] == 0) {
1324 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1325 exit(EXIT_FAILURE);
1326 }
1327
1328 switch ((*type)[i]) {
1329 case SDDS_SHORT:
1330 (*data)[i].shortData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1331 break;
1332 case SDDS_USHORT:
1333 (*data)[i].ushortData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1334 break;
1335 case SDDS_LONG:
1336 (*data)[i].longData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1337 break;
1338 case SDDS_ULONG:
1339 (*data)[i].ulongData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1340 break;
1341 case SDDS_LONG64:
1342 (*data)[i].long64Data = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1343 break;
1344 case SDDS_ULONG64:
1345 (*data)[i].ulong64Data = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1346 break;
1347 case SDDS_FLOAT:
1348 (*data)[i].floatData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1349 break;
1350 case SDDS_DOUBLE:
1351 (*data)[i].doubleData = SDDS_GetColumn(SDDS_orig, orig_column_name[i]);
1352 break;
1353 default:
1354 continue;
1355 }
1356 if (strncmp(colPrefix, orig_column_name[i], strlen(colPrefix)) != 0) {
1357 continue;
1358 }
1359 (*colIndex)[*validColumns] = i;
1360 *validColumns += 1;
1361 }
1362
1363 if (SDDS_Terminate(SDDS_orig) != 1) {
1364 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1365 exit(EXIT_FAILURE);
1366 }
1367
1368 if (*validColumns == 0) {
1369 fprintf(stderr, "error: no valid columns in image file\n");
1370 exit(EXIT_FAILURE);
1371 }
1372 /* For each valid column read the row from the column name */
1373 for (i = 0; i < *validColumns; i++) {
1374 (*colIndex2)[i] = atof(orig_column_name[(*colIndex)[i]] + strlen(colPrefix));
1375 }
1376 /* Sort columns by their row value ('row' is not related to the SDDS rows) */
1377 for (i = 0; i < *validColumns; i++) {
1378 for (j = i + 1; j < *validColumns; j++) {
1379 if ((*colIndex2)[j] < (*colIndex2)[i]) {
1380 temp = (*colIndex)[i];
1381 temp2 = (*colIndex2)[i];
1382 (*colIndex)[i] = (*colIndex)[j];
1383 (*colIndex2)[i] = (*colIndex2)[j];
1384 (*colIndex)[j] = temp;
1385 (*colIndex2)[j] = temp2;
1386 }
1387 }
1388 }
1389 return rows;
1390}
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
void * SDDS_GetColumn(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves a copy of the data for a specified column, including only rows marked as "of interest".
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:50
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
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.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
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.
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
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.
int32_t SDDS_GetNamedColumnType(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the data type of a column in the SDDS dataset by its name.
char ** SDDS_GetColumnNames(SDDS_DATASET *SDDS_dataset, int32_t *number)
Retrieves the names of all columns in the SDDS dataset.
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
#define SDDS_ULONG
Identifier for the unsigned 32-bit integer data type.
Definition SDDStypes.h:67
#define SDDS_FLOAT
Identifier for the float data type.
Definition SDDStypes.h:43
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
#define SDDS_ULONG64
Identifier for the unsigned 64-bit integer data type.
Definition SDDStypes.h:55
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
#define SDDS_SHORT
Identifier for the signed short integer data type.
Definition SDDStypes.h:73
#define SDDS_USHORT
Identifier for the unsigned short integer data type.
Definition SDDStypes.h:79
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
#define SDDS_LONG64
Identifier for the signed 64-bit integer data type.
Definition SDDStypes.h:49
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
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.
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:357
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:391
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.