1 /*
2 * ColourBrightness
3 *
4 * Based on the program PTStitcher by Helmut Dersch.
5 *
6 * It is intended to duplicate the functionality of original program
7 *
8 * Dec 2005
9 *
10 * This program is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU General Public
12 * License as published by the Free Software Foundation; either
13 * version 2 of the License, or (at your option) any later version.
14 *
15 * This software is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with this software; see the file COPYING. If not, a copy
22 * can be downloaded from http://www.gnu.org/licenses/gpl.html, or
23 * obtained by writing to the Free Software Foundation, Inc.,
24 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 *
26 *
27 * Author: Daniel M German dmgerman at uvic doooot ca
28 *
29 */
30
31
32 #include "filter.h"
33 #include "PTcommon.h"
34 #include "ColourBrightness.h"
35 #include "pt_stdint.h"
36
37 #include "tiffio.h"
38 #include <assert.h>
39 #include <math.h>
40
41 #include "pttiff.h"
42
43 #ifdef WIN32
44 #ifdef _MSC_VER
45
46 // MSVC doesn't support round()
47 //#define round(x) ( (int) (x+0.5) )
48 #define round(x) (int)(x)
49 #endif // _MSC_VER
50
51 // MSVC wants htons() to be a library function
52 // here we define it as a macro instead
53 #undef htons
54 #undef htonl
55 // byte reordering macros -- avoids loading sockets lib on Windows
56 #define LITTLE_ENDIAN // change if your Windows box is from Mars
57 #if defined(BIG_ENDIAN) && !defined(LITTLE_ENDIAN)
58 #define htons(A) (A)
59 #define htonl(A) (A)
60 #define ntohs(A) (A)
61 #define ntohl(A) (A)
62 #elif defined(LITTLE_ENDIAN) && !defined(BIG_ENDIAN)
63 #define htons(A) (((uint16)(A) & 0xff00) >> 8 | ((uint16)(A) & 0x00ff) << 8 )
64 #define htonl(A) ((((uint32)(A) & 0xff000000) >> 24) | \
65 ((uint32)(A) & 0x00ff0000) >> 8 | \
66 ((uint32)(A) & 0x0000ff00) << 8 | \
67 ((uint32)(A) & 0x000000ff) << 24)
68 #define ntohs htons
69 #define ntohl htohl
70 #else
71 #error "Either BIG_ENDIAN or LITTLE_ENDIAN must be #defined, but not both."
72 #endif
73
74 #endif //def WIN32
75
76 #if defined(__DragonFly__)
77 #include <sys/endian.h>
78 #include <arpa/inet.h>
79 #endif
80
81 FILE *debugFile = 0;
82
83
84 #ifdef __TESTING__
85 #include <dmalloc.h>
86 #endif
87
88
89 #define IDX_RED 0
90 #define IDX_GREEN 1
91 #define IDX_BLUE 2
92 #define IDX_INTENSITY 3
93 #define IDX_SATURATION 4
94 #define IDX_HUE 5
95
96
InitializeMagnolia(int numberImages,int size,calla_function parm2)97 magnolia_struct *InitializeMagnolia(int numberImages, int size, calla_function parm2)
98 {
99 int j;
100 int i;
101 magnolia_struct *magnolia = NULL;
102 double *ptrDouble;
103 double var16;
104 int ecx;
105
106 if ((magnolia = malloc(numberImages * sizeof(magnolia_struct) )) == 0) {
107
108 return 0;
109
110 }
111
112 var16 = (size -1 ) / 255.0; /// shouldn't this be a 256?
113
114 for (i=0; i < numberImages; i++) {
115
116 magnolia[i].components = size;
117
118 magnolia[i].function = parm2;
119
120 for (j=0; j<6 ; j++) {
121
122 if ((ptrDouble = calloc(size, sizeof(double))) == NULL) {
123 /// @TODO clean up memory leak
124 return NULL;
125 }
126
127 assert( magnolia[i].components == size);
128
129 for (ecx = 0; ecx < size; ecx ++) {
130
131 ptrDouble[ecx] = ecx * var16;
132
133 } // if
134
135 magnolia[i].fieldx04[j] = ptrDouble;
136
137 } // for (j=0; j<6; j++) {
138
139 } // end for (i=0; i < numberImages; i++) {
140
141 return magnolia;
142
143 }
144
RemapPoint(int value,double mapTable[])145 int RemapPoint(int value, double mapTable[])
146 {
147
148 double delta;
149
150 double deltaPrev;
151 double deltaNext;
152
153 double var48;
154
155 double tablePrevValue;
156 double tableNextValue;
157 int nextValueInt;
158 int prevValueInt;
159
160 int edx;
161
162 // Find previous and next. Extrapolate if necessary
163 #if 0
164 int i;
165 printf("Colour map\n");
166 for (i=0;i < 255; i++) {
167 printf("%d: %7.3f ", i, mapTable[i]);
168 }
169 #endif
170 if ( value == 0 ) {
171
172 tablePrevValue = 2* mapTable[0] - mapTable[1]; // XXXXXXXXXXXXXXXXX This just looks backwards!
173
174 } else {
175
176 tablePrevValue = mapTable[value-1];
177
178 }
179
180 if ( value == 0xff ) {
181
182 tableNextValue = 2 * mapTable[255] - mapTable[254];
183
184
185 } else {
186
187 tableNextValue = mapTable[value+1];
188
189 }
190
191 if (fabs(tableNextValue - tablePrevValue) <= 2.0) {
192 // if the difference |f(value - 1) -f(value+1) is too small
193
194 if ((int)mapTable[value] == 0xff ) {
195 return 0xff;
196 }
197
198 delta = mapTable[value] - (int)mapTable[value];
199
200 // chance of being one colour or the next.
201 // of n points, n * delta = base and n(1-*delta) will be base + 1
202 // this guarantees that the distribution is faithfull. clever.
203
204 if ( delta * RAND_MAX < rand() ) {
205 // if (C0 == 1)
206 return (int)mapTable[value];
207
208 } else {
209
210 return ((int)mapTable[value]) + 1;
211
212 }
213 assert(0); // nothing should reach here
214
215 } // if (value ... 2.0) {
216
217 // THIS CASE IS WHEN THE TANGENT is > 1
218
219
220 nextValueInt = (int)tableNextValue;
221
222 if ( (int)tableNextValue > 0xff ) {
223
224 nextValueInt = 0xff;
225
226 }
227
228 prevValueInt = (int)tablePrevValue;
229
230 if ( (int)tablePrevValue < tablePrevValue ) {
231 prevValueInt ++;
232 }
233
234 // prevValueInt == ceiling(tablePrevValue)
235
236 if ( prevValueInt < 0 ) {
237 prevValueInt = 0;
238 }
239
240
241 deltaPrev = mapTable[value] - tablePrevValue;
242 deltaNext = tableNextValue - mapTable[value];
243
244 edx = prevValueInt;
245
246 var48 = 0;
247
248 // if ( %edx > %nextValueInt ) /// [(int)tablePrevValue] > [(int)tableNextValue]
249 while ( edx <= nextValueInt ) { /// [(int)tablePrevValue] > [(int)tableNextValue]
250
251 if (edx < mapTable[value]) {
252 var48 += (edx - tablePrevValue)/deltaPrev;
253 } else {
254 var48 += (tableNextValue - edx)/deltaNext;
255 }
256 edx ++;
257 } // while...
258
259 // where did edx = ebx go?
260
261 var48 = (rand() * var48) / RAND_MAX;
262
263 edx = prevValueInt;
264
265 while ( edx <= nextValueInt ) {
266
267 if (edx < mapTable[value]) {
268 var48 -= (edx - tablePrevValue) / deltaPrev;
269 } else { // if ( %ah == 0x1 ) {
270 var48 -= (tableNextValue - edx) / deltaNext;
271 } // if ( %ah == 0x1 )
272
273 if ( var48 < 0 ) {
274 return edx;
275 }
276 edx ++;
277
278 } // while ( %edx <= %nextValueInt )
279
280 return nextValueInt;
281
282 }
283
284
RemapDouble(double value,double mapTable[])285 double RemapDouble(double value, double mapTable[])
286 {
287
288 double delta, deltaY, tableNextValue;
289 int valueInt;
290
291 if (!(value >=0.0 && value <= 255.0)) {
292 printf("Wrong value %f\n", value);
293 }
294
295 assert(value >=0.0 && value <= 255.0);
296
297
298 valueInt = (int)value;
299
300 assert(valueInt >=0 && valueInt <= 255);
301
302 if ( value == 0xff ) {
303 tableNextValue = 2 * mapTable[255] - mapTable[254];
304 } else {
305 tableNextValue = mapTable[valueInt+1];
306 }
307 deltaY = tableNextValue - mapTable[valueInt];
308 assert(deltaY>=0);
309
310 delta = (value - valueInt) * deltaY;
311
312 return mapTable[valueInt] + delta;
313
314 }
315
316
317
Unknown47(unsigned char parm0,unsigned char parm1,unsigned char parm2)318 unsigned char Unknown47(unsigned char parm0, unsigned char parm1, unsigned char parm2)
319 {
320
321 int eax = parm1;
322 int edx = parm2;
323 int ecx = parm0;
324
325
326 ecx = ecx * 3;
327 ecx = ecx + 2 * eax - 256;
328 ecx = (ecx + 2 * edx - 256) * 2 /3;
329
330 if (ecx >= 0) {
331 if (ecx > 0xff)
332 return 0xff;
333 else
334 return ecx;
335 } else {
336 return 0;
337 }
338 }
339
340
341
DisplayHistogramsError(int numberHistograms,histograms_struct * ptrHistograms)342 void DisplayHistogramsError(int numberHistograms, histograms_struct * ptrHistograms)
343 {
344 int index;
345 histogram_type * ptrOtherHistograms;
346 histogram_type * ptrBaseHistograms;
347 int currentColour;
348 int i;
349
350 for (index = 0; index < numberHistograms; index++ ) {
351
352 // if the number of overlapping pixels is less than 1k then it ignores the images
353
354 if ( ptrHistograms[index].overlappingPixels > 999 ) {
355
356 fprintf(debugFile, "Histogram %d Images %d %d, %d Pixels: ", index ,
357 ptrHistograms[index].baseImageNumber,
358 ptrHistograms[index].otherImageNumber,
359 ptrHistograms[index].overlappingPixels);
360
361 ptrBaseHistograms = &(ptrHistograms[index].ptrBaseHistograms);
362 ptrOtherHistograms = &(ptrHistograms[index].ptrOtherHistograms);
363
364 for (currentColour = 0; currentColour < 6; currentColour++) {
365
366 double sum = 0;
367 int *var192;
368 int *var200;
369
370 var192 = (*ptrBaseHistograms)[currentColour];
371
372 var200 = (*ptrOtherHistograms)[currentColour];
373
374 for (i =0; i < 0x100; i++) {
375
376 int diff = var192[i] - var200[i];
377
378 sum += diff * diff;
379
380 }
381
382 fprintf(debugFile, " %g", sum * 1.0 /ptrHistograms[index].overlappingPixels);
383
384 } //for (currentColour = 0; currentColour < 3; currentColour++ {
385
386 fprintf(debugFile, "\n");
387
388 } // if ( > 999)
389 } // for (index = 0; index < numberHistograms; index++ ) {
390
391 }
392
OutputPhotoshopCurve(FILE * output,int size,double * curve)393 int OutputPhotoshopCurve(FILE *output, int size, double *curve)
394 {
395 uint16_t shortValue;
396 uint16_t x;
397 int i;
398
399 // so far we only support size == 256
400
401 // fprintf(stderr, "size %d\n", size);
402
403 assert(size == 256);
404
405 // I am not sure why I can't write more than 14 points to the curves file.
406 // it might have to do with the way a spline is computed
407 // The algorithm currently writes the first point, 12, and the last point
408
409
410 shortValue = (uint16_t) htons(14);
411
412 if (fwrite(&shortValue, 2, 1, output) != 1) {
413 goto error;
414 }
415
416
417 // for (i = 0; i< size; i+=16) {
418 for (i = 0; i< size; i+=20) {
419 int temp;
420 uint16_t y;
421 uint16_t x;
422
423 // fprintf(stderr, "value of curve %i: %10.5f\n", i, curve[i]);
424
425 temp = round(curve[i]);
426
427 // be paranoic
428 assert(temp >= 0 && temp <= 255);
429
430 y = (uint16_t) htons(temp);
431 x = (uint16_t) htons(i);
432
433 // For some reason y is first in the output
434 if (fwrite(&y, 2, 1, output) != 1 ||
435 fwrite(&x, 2, 1, output) != 1 ) {
436 goto error;
437 }
438 }
439
440 // Write the very last point
441 x=htons(255);
442 if (fwrite(&x, 2, 1, output) != 1 ||
443 fwrite(&x, 2, 1, output) != 1 ) {
444 goto error;
445 }
446
447 return 1;
448
449 error:
450 PrintError("Error writing to curves file");
451 return 0;
452
453 }
454
OutputPhotoshopFlatArbitraryMap(FILE * output)455 int OutputPhotoshopFlatArbitraryMap(FILE *output)
456 {
457 unsigned int i;
458 for (i = 0; i< 256; i++) {
459 unsigned char c;
460 c = i;
461 if (fputc(c, output) != c) {
462 PrintError("Error writing to curves file");
463 return 0;
464 }
465 }
466 return 1;
467 }
468
OutputPhotoshopArbitraryMap(FILE * output,int size,double * curve)469 int OutputPhotoshopArbitraryMap(FILE *output, int size, double *curve)
470 {
471 int i;
472
473 // so far we only support size == 0xff
474
475 assert(size == 256);
476
477 // for (i = 0; i< size; i+=16) {
478 for (i = 0; i< size; i++) {
479 unsigned int temp;
480 unsigned int value;
481
482 // fprintf(stderr, "value of curve %i: %10.5f\n", i, curve[i]);
483
484 temp = round(curve[i]);
485
486 // be paranoic
487 assert(temp >= 0 && temp <= 255);
488 value = temp;
489
490 // For some reason y is first in the output
491 if (fputc(value, output) != value) {
492 PrintError("Error writing to curves file");
493 return 0;
494 }
495 }
496
497 return 1;
498 }
499
500
OutputEmptyPhotoshopCurve(FILE * output)501 int OutputEmptyPhotoshopCurve(FILE *output)
502 {
503 // The empty curve is 2 entries: 0,0, and 255,255
504 // It is easier to write it in one go
505
506 #define EMPTY_CURVE "\x00\x02\x00\x00\x00\x00\x00\xff\x00\xff"
507 #define LENGTH_EMPTY_CURVE_DATA 10
508
509 // Make sure it is the right size
510 // Remember that strings are padded with \0 at the end
511
512 if (fwrite(EMPTY_CURVE, LENGTH_EMPTY_CURVE_DATA, 1, output) != 1) {
513 PrintError("Error writing to curves file");
514 return 0;
515 }
516
517 return 1;
518 }
519
OutputCurves(int index,magnolia_struct * curves,int typeOfCorrection,fullPath * panoFileName,int curveType)520 int OutputCurves(int index, magnolia_struct *curves, int typeOfCorrection,
521 fullPath *panoFileName, int curveType)
522 {
523 char outputFileName[512];
524 char temp[12];
525 int i;
526 FILE *output;
527 char *curveExtension[2] = {
528 ".amp",
529 ".acv"};
530
531 #define PHOTOSHOP_CURVES_MAGIC_NUMBER "\x00\x04\x00\x05"
532
533 // TODO: panotool usually creates a temp file, and then renames it to the
534 // actual name. This is not done yet
535
536 strncpy(outputFileName, panoFileName->name, 500);
537 sprintf(temp, "%04d", index);
538 strcat(outputFileName, temp);
539
540 panoReplaceExt(outputFileName, curveExtension[curveType-1]);
541
542 // fprintf(stderr, "Creating output file %s\n", outputFileName);
543
544 if ((output = fopen(outputFileName, "w+")) == NULL) {
545 PrintError("Unable to create output curves file %s", outputFileName);
546 return 0;
547 }
548
549 switch (curveType) {
550 case CB_OUTPUT_CURVE_SMOOTH:
551 if (fwrite(PHOTOSHOP_CURVES_MAGIC_NUMBER, 4, 1, output) != 1)
552 goto error;
553
554 // Output main curve
555 if (OutputEmptyPhotoshopCurve(output) == 0)
556 goto error;
557
558
559 // Output Red, Green and Blue Channels
560 for (i=0;i<3;i++) {
561 if (OutputPhotoshopCurve(output, curves->components, curves->fieldx04[i]) == 0)
562 goto error;
563 }
564 // it needs an empty curve at the end, and I donot know why
565 if (OutputEmptyPhotoshopCurve(output) == 0) {
566 PrintError("Unable to create output curves file %s", outputFileName);
567 return 0;
568 }
569 break;
570 case CB_OUTPUT_CURVE_ARBITRARY:
571 if (OutputPhotoshopFlatArbitraryMap(output) == 0)
572 goto error;
573
574 // Output Red, Green and Blue Channels
575 for (i=0;i<3;i++) {
576 if (OutputPhotoshopArbitraryMap(output, curves->components, curves->fieldx04[i]) == 0)
577 goto error;
578 }
579 break;
580 }
581
582 fclose(output);
583 return 1;
584
585
586 error:
587 PrintError("Unable to output curves file %s", outputFileName);
588 return 0;
589 }
590
591
592
ColourBrightness(fullPath * fullPathImages,fullPath * outputFullPathImages,int counterImages,int indexReferenceImage,int parm3,int createCurvesType)593 void ColourBrightness( fullPath *fullPathImages, fullPath *outputFullPathImages,
594 int counterImages, int indexReferenceImage, int parm3, int createCurvesType)
595 {
596
597 histograms_struct * ptrHistograms2;
598 int numberHistograms;
599 int index;
600 calla_struct calla;
601 char string[128];
602 int i;
603 extern FILE *debugFile; // 0x8054600
604
605 numberHistograms = ((counterImages-1) * counterImages)/2;
606
607 if(debugFile)
608 {
609 fclose(debugFile);
610 debugFile = 0;
611 }
612
613 debugFile = fopen("Debug.txt", "w");
614 // debugFile = stderr;
615
616 fprintf(debugFile, "Entering function \"colourbrightness\" with %d images, nfix =%d\n", counterImages, indexReferenceImage);
617
618 calla.ptrHistograms = ReadHistograms(fullPathImages, counterImages);
619
620 if ( calla.ptrHistograms == 0 )
621 return ;
622
623 fprintf(debugFile, "\nQuality before optimization:\n");
624
625 // printf("\nQuality before optimization:\n");
626
627 index = 0;
628
629 DisplayHistogramsError(numberHistograms, calla.ptrHistograms);
630
631 ///////////////////////////////////////////
632
633 calla.fullPathImages = fullPathImages;
634
635 calla.numberImages = counterImages;
636
637 calla.indexReferenceImage = indexReferenceImage; //
638
639 calla.magnolia = InitializeMagnolia(counterImages, 0x100, MapFunction);
640
641 if (calla.magnolia == 0 )
642 return ;
643
644 if (ComputeColourBrightnessCorrection(&calla) == 0)
645 return ;
646
647 fprintf(debugFile, "\nResults of Optimization:");
648
649 if (createCurvesType !=0) {
650 fprintf(debugFile, "\nResults of Optimization:");
651
652 for (index=0; index < counterImages; index++ ) {
653 if (OutputCurves(index, &(calla.magnolia[index]), parm3, &(outputFullPathImages[index]), createCurvesType) == 0) {
654 PrintError("Error creating curves files");
655 return ;
656 }
657 }
658 return ;
659 }
660
661 // printf("\nQuality before optimization: 12 3\n");
662
663 for (index=0; index < counterImages; index++ ) {
664 int i;
665 int histo;
666 magnolia_struct *magnolias;
667
668 magnolias = calla.magnolia; // ecx
669
670 for (histo = 0; histo< 6; histo++) {
671 switch (histo) {
672 case 0:
673 fprintf(debugFile, "\nImage %d:\nRed Channel: ", index);
674 break;
675 case 1:
676 fprintf(debugFile, "\nImage %d:\nGreen Channel: ", index);
677 break;
678 case 2:
679 fprintf(debugFile, "\nImage %d:\nBlue Channel: ", index);
680 break;
681 case IDX_INTENSITY:
682 fprintf(debugFile, "\nImage %d:\nIntensity: ", index);
683 break;
684 case IDX_SATURATION:
685 fprintf(debugFile, "\nImage %d:\nSauration: ", index);
686 break;
687 case IDX_HUE:
688 fprintf(debugFile, "\nImage %d:\nHue: ", index);
689 break;
690 default:
691 assert(0);
692 }
693
694 for (i = 0; i < magnolias[index].components; i++) {
695
696 double *array;
697
698 array = magnolias[index].fieldx04[histo];
699
700 fprintf(debugFile, "%g ", array[i]);
701
702 } // while (( %edi < (%ecx+ebx) ) {
703 }
704
705 } // for (index=0; index < counterImages; index++ ) {
706
707 if ( ptQuietFlag == 0 ) {
708 switch (parm3) {
709 case 0:
710 Progress(_initProgress, "Adjusting Colour and Brightness");
711 break;
712 case 1:
713 Progress(_initProgress, "Adjusting Brightness");
714 break;
715 case 2:
716 Progress(_initProgress, "Adjusting Saturation");
717 break;
718 default:
719 assert(0);
720 }
721 }
722 index = 0;
723
724 for (index = 0; index <counterImages; index++ ) {
725
726 sprintf(string, "%d", index * 100 / counterImages);
727
728 if ( ptQuietFlag == 0 ) {
729
730 if (Progress(_setProgress, string) == 0)
731 return ;
732
733 } //if ( ptQuietFlag == $0x0 )
734
735
736 if (strcmp(fullPathImages[index].name,
737 outputFullPathImages[index].name) != 0 ||
738 index != indexReferenceImage ) {
739
740
741 // Do the correction if the input and output filename are different
742 // or if it is not the reference image
743
744 /// Due to laiziness we are processing the reference image in order to copy it from input to output
745 // We should avoid it, but it is not a big deal
746
747 //printf("To correct index %d %d indexReferenceImage\n", index, indexReferenceImage);
748 if (CorrectFileColourBrightness(&fullPathImages[index], &outputFullPathImages[index],
749 &calla.magnolia[index], parm3) != 0) {
750
751 PrintError("Error creating colour corrected file\n");
752 return;
753 }
754 }
755 } // for (index = 0; index <counterImages; index++ ) {
756
757
758 ptrHistograms2 = ReadHistograms(outputFullPathImages, counterImages);
759
760 fprintf(debugFile, "\nQuality after optimization:\n");
761
762 DisplayHistogramsError(numberHistograms, ptrHistograms2);
763
764 FreeHistograms(calla.ptrHistograms, numberHistograms);
765 FreeHistograms(ptrHistograms2, numberHistograms);
766
767 for (i = 0; i <counterImages; i++ ) {
768
769 for (index = 0; index < 6; index ++)
770 free(calla.magnolia[i].fieldx04[index]);
771
772 }
773 free(calla.magnolia);
774
775 }
776
FreeHistograms(histograms_struct * ptrHistograms,int count)777 void FreeHistograms(histograms_struct *ptrHistograms, int count)
778 {
779 int i;
780 int index;
781 histograms_struct *savedPtrHistograms = ptrHistograms;
782
783 for (i = 0; i < count; i++ ) {
784 for (index = 0; index < 6; index ++) {
785 free(ptrHistograms->ptrBaseHistograms[index]);
786 free(ptrHistograms->ptrOtherHistograms[index]);
787 }
788 ptrHistograms++;
789 } // for (i = 0; i < numberHistograms; i++ ) {
790 free(savedPtrHistograms);
791 }
792
793
CorrectFileColourBrightness(fullPath * inPath,fullPath * outPath,magnolia_struct * magnolia,int parm3)794 int CorrectFileColourBrightness(fullPath *inPath, fullPath *outPath, magnolia_struct *magnolia, int parm3)
795 {
796 Image image;
797 char tempString[512];
798 if (panoTiffRead (&image, inPath->name) == 0) {
799 sprintf(tempString, "Could not read TIFF file %s", inPath->name);
800 PrintError(tempString);
801 return -1;
802 }
803
804 CorrectImageColourBrigthness(&image, magnolia, parm3);
805
806 if (panoTiffWrite(&image, outPath->name) == 0) {
807 PrintError("Could not read TIFF file %s", inPath->name);
808 panoImageDispose(&image);
809 return -1;
810 }
811
812 panoImageDispose(&image);
813 return(0);
814
815 }
816
817
FindNextCandidate(int candidates[],calla_struct * calla)818 int FindNextCandidate(int candidates[], calla_struct *calla)
819 {
820
821 // Find image not yet corrected with the maximum overlap over the
822 // corrected area
823
824 // The algorithm is simple.
825
826 // For each daisy
827 // if one of the images is in, but no the other
828 // edi[not in image] += overlap
829
830 // return image with largest overlap
831 // if not found return -1
832
833 int *overlapping;
834 int i;
835 int max;
836 int overlappingPixels;
837 int baseImage;
838 int otherImage;
839 int returnValue;
840 int numberDaisies;
841
842 histograms_struct *ptrHistograms;
843
844 // Candidates is an array of boolean. If true the image has been processed
845
846 // A daisy exists for each pair of images
847 // It contains information about their ovelap
848
849 numberDaisies = calla->numberImages * (calla->numberImages -1)/2;
850
851 if ((overlapping = malloc(calla->numberImages * sizeof(int))) == 0) {
852 PrintError("Not enough memory\n");
853 return -1; // not more memory
854 }
855
856 for (i = 0; i < calla->numberImages; i++) {
857
858 // Number of pixels ovelapping for this image with corrected area
859 overlapping[i] = 0;
860
861 } // for (i = 0; i < calla->numberImages; i++) {
862
863
864 for (i = 0; i < numberDaisies; i++) {
865
866 ptrHistograms = calla->ptrHistograms;
867
868 // How many pixels overlap
869 overlappingPixels = ptrHistograms[i].overlappingPixels;
870 baseImage = ptrHistograms[i].baseImageNumber;
871 otherImage = ptrHistograms[i].otherImageNumber;
872
873 assert(baseImage < calla->numberImages);
874 assert(otherImage < calla->numberImages);
875 assert(baseImage >= 0);
876 assert(otherImage >= 0);
877 assert(baseImage != otherImage);
878
879 // fprintf(stderr, "Overlap %2d:%2d:%7d\n", baseImage, otherImage, overlappingPixels);
880
881 if ( overlappingPixels > 1000 ) {
882
883 // WE ONLY Process images with an overlap of at least 1000 pixels
884
885 if (candidates[baseImage] == 0 ||
886 candidates[otherImage] != 0 ) {
887
888 // baseImageNumber not in OR other in
889
890 // Here we have four alternatives:
891
892 // 1. baseImageNumber not AND other not <- we need to skip this case
893 // 2. baseImageNumber not AND other in
894
895 // 3. baseImageNumber in AND other in <- we need to skip this case
896 // 4. baseImageNumber not AND other in
897
898
899 // so the condition is: not baseImageNumber and otherIn
900
901
902 if (candidates[otherImage] != 0 &&
903 candidates[baseImage] == 0) {
904
905 overlapping[baseImage]+=overlappingPixels;
906 }
907
908 } else {
909
910 // This code is executed if:
911
912 // NOT (baseImageNumber == 0 OR otherImageNumber == 1) =>
913 // NOT (baseImageNumber == 0 ) AND NOT (otherImageNumber == 1) =>
914 // baseImageNumber in AND otherImageNumber is not in
915
916 // if the base is in, AND not the other, then add it to the other
917
918 overlapping[otherImage]+=overlappingPixels;
919
920 }
921
922 } // if ( overlappingPixels > 1000 ) {
923
924
925 } // for (i = 0; i < numberDaisies; i++) {
926
927
928 returnValue = -1;
929
930 // Find maximum
931 max = 0;
932 for (i =0; i < calla->numberImages ; i++ ) {
933
934 // fprintf(stderr, "Overlap Image %2d:%7d\n", i, overlapping[i]);
935
936 if ( overlapping[i] > max ) {
937
938 max = overlapping[i];
939
940 returnValue = i;
941
942 }
943
944 } // for (i =0; i < calla->numberImages ; i++ ) {
945
946 free(overlapping);
947
948 // We need to assert the value
949
950 if(returnValue >= 0) {
951 assert(returnValue < calla->numberImages);
952 assert(candidates[returnValue] == 0);
953 }
954
955
956 return returnValue;
957
958
959 }
960
961
962
ComputeColourBrightnessCorrection(calla_struct * calla)963 int ComputeColourBrightnessCorrection(calla_struct *calla)
964 {
965
966
967 double *remappedSourceHistogram = 0;
968 double *accumToCorrectHistogram = 0;
969 double *accumSourceHistogram = 0;
970 int *processedImages = 0;
971
972 int currentImageNumber;
973 int channel;
974 int numberIntersections;
975 histograms_struct *currentHistogram = 0;
976 int j;
977
978 int **ptrHistogram;
979 int *array;
980 int i;
981
982 // How many overlaps do we expect
983
984 numberIntersections = ((calla->numberImages - 1) * calla->numberImages)/2;
985
986 // Keep an array of booleans that keeps track of what we have done
987 processedImages = calloc(calla->numberImages, sizeof(int));
988
989 // Histograms for source and to correct
990 accumToCorrectHistogram = malloc(0x100 * sizeof(double));
991 accumSourceHistogram = malloc(0x100 * sizeof(double));
992
993 // Histogram of the to-correct when it has been remapped
994 remappedSourceHistogram = malloc(0x100 * sizeof(double));
995
996
997 if ( processedImages == 0
998 || accumToCorrectHistogram == 0
999 || accumSourceHistogram == 0
1000 || remappedSourceHistogram == 0 )
1001 {
1002 if ( processedImages != 0 )
1003 free(processedImages);
1004
1005 if ( remappedSourceHistogram != 0 )
1006 free(remappedSourceHistogram);
1007
1008 if ( accumToCorrectHistogram != 0 )
1009 free(accumToCorrectHistogram);
1010
1011 if ( accumSourceHistogram != 0 )
1012 free(accumSourceHistogram);
1013
1014 return 0;
1015 }
1016
1017 // Mark starting image as done
1018 processedImages[calla->indexReferenceImage] = 1;
1019
1020 // We keep repeating this loop while there are "candidates"
1021 // A candidate is the image with the largest overlap with the already corrected image
1022
1023 while ((currentImageNumber = FindNextCandidate(processedImages, calla)) != -1) {
1024
1025 // We have a candidate image: currentImageNumber
1026
1027 // fprintf(stderr, "We are processing image %d\n", currentImageNumber);
1028
1029 // make sure it is a valid image number and it has not been processed
1030 assert(currentImageNumber >= 0);
1031 assert(currentImageNumber < calla->numberImages);
1032 assert(processedImages[currentImageNumber] == 0);
1033
1034 // For every channel it does the correction independently.
1035
1036 for (channel = 0; channel < 6; channel++) {
1037
1038 int currentIntersection;
1039
1040 // clean the accum histograms
1041 for (i =0; i <= 0xff; i ++) {
1042
1043 accumSourceHistogram[i] = 0;
1044 accumToCorrectHistogram[i] = 0;
1045
1046 }
1047
1048 // for each intersection between 2 images
1049
1050 for (currentIntersection = 0; currentIntersection < numberIntersections ; currentIntersection++) {
1051
1052 currentHistogram = &(calla->ptrHistograms[currentIntersection]);
1053
1054
1055 // DO it only if the overlap is bigger than 1k pixels
1056
1057 if ( currentHistogram->overlappingPixels > 1000 ) { // it is not consistent XXX
1058
1059 if ( currentHistogram->baseImageNumber == currentImageNumber &&
1060 processedImages[currentHistogram->otherImageNumber] != 0 ) {
1061
1062 // this means the otherImage has been already processed
1063 // but not baseImageNumber
1064
1065 // REMAP histogram according to current mapping function
1066
1067 // fprintf(stderr, "Original histogram Channel %d [%d,Base %d]: \n", channel, currentImageNumber, currentHistogram->otherImageNumber);
1068 // for (j = 0; j<0x100; j++) {
1069 // fprintf(stderr, " %d:%d", j, currentHistogram->ptrOtherHistograms[channel][j]);
1070 // }
1071 // fprintf(stderr, "\n");
1072
1073
1074
1075 // Remap histogram of other image
1076
1077 RemapHistogram(currentHistogram->ptrOtherHistograms[channel],
1078 remappedSourceHistogram,
1079 &(calla->magnolia[currentHistogram->otherImageNumber]), channel);
1080
1081 // fprintf(stderr, "\n");
1082
1083 // fprintf(stderr, "Remapped histogram :\n");
1084
1085 //for (j = 0; j<0x100; j++) {
1086 //fprintf(stderr, " %d:%g", j, remappedSourceHistogram[j]);
1087 //}
1088 //fprintf(stderr, "\n");
1089
1090 // Add source Histogram to the Accumulated
1091
1092 for (j = 0; ( j <= 0xff ); j ++) {
1093
1094 accumSourceHistogram[j] += remappedSourceHistogram[j];
1095
1096 }
1097
1098 // Add target histogram to its accumulated
1099 // This guarantees that both accum histograms have the same total
1100 // number of points
1101
1102
1103 ptrHistogram = calla->ptrHistograms[currentIntersection].ptrBaseHistograms;
1104
1105 for (j = 0; j <= 0xff; j ++) {
1106
1107 array = ptrHistogram[channel];
1108 accumToCorrectHistogram[j] += array[j];
1109
1110 }
1111
1112 continue;
1113
1114
1115 } else { //
1116
1117 assert(currentHistogram == &calla->ptrHistograms[currentIntersection]);
1118
1119 if (currentHistogram->otherImageNumber == currentImageNumber &&
1120 processedImages[currentHistogram->baseImageNumber] != 0 ) {
1121
1122 // MIRROR version of the code above.
1123 // In this case baseImageNumber is done
1124 // and otherImageNumber is to be done
1125
1126 //fprintf(stderr, "Original histogram Channel %d [Base %d,%d]: \n", channel, currentImageNumber, currentHistogram->baseImageNumber);
1127
1128 //for (j = 0; j<0x100; j++) {
1129 //fprintf(stderr, " %d:%d", j, currentHistogram->ptrBaseHistograms[channel][j]);
1130 //}
1131 //fprintf(stderr, "\n");
1132
1133 // Remap source histogram
1134 RemapHistogram(currentHistogram->ptrBaseHistograms[channel],
1135 remappedSourceHistogram,
1136 &(calla->magnolia[currentHistogram->baseImageNumber]),
1137 channel);
1138
1139 //fprintf(stderr, "Remapped histogram: \n");
1140
1141 //for (j = 0; j<0x100; j++) {
1142 // fprintf(stderr, " %d:%g", j, remappedSourceHistogram[j]);
1143 //}
1144 //fprintf(stderr, "\n");
1145
1146 // Add it to the accumulated
1147 for (j = 0; j <= 0xff; j ++) {
1148
1149 accumSourceHistogram[j] += remappedSourceHistogram[j];
1150
1151 } // for (j = 0; j <= 0xff; j ++) {
1152
1153 ptrHistogram = currentHistogram->ptrOtherHistograms;
1154
1155
1156 // add target histogram to its accumulated
1157 // This guarantees that both accum histograms have the same total
1158 // number of points
1159 for (j = 0; j <= 0xff; j ++) {
1160
1161 array = ptrHistogram[channel];
1162 accumToCorrectHistogram[j] += array[j];
1163
1164 } // for (j = 0; j <= 0xff; j ++) {
1165
1166 } // end of if
1167
1168 } // end of else
1169 } // end of if
1170
1171 } // for (currentIntersection = 0; %currentIntersection < numberIntersections ; currentIntersection++) {
1172
1173
1174 ComputeAdjustmentCurve(accumToCorrectHistogram,
1175 accumSourceHistogram,
1176 (calla->magnolia[currentImageNumber].fieldx04)[channel]);
1177
1178 } // for (channel = 0; var < 6; var++) {
1179
1180
1181
1182 processedImages[currentImageNumber] = 1;
1183
1184 } // while (Unknown43 (...) != -1);
1185
1186 for (i =0 ; i< calla->numberImages; i++) {
1187 // are all the images opimized?
1188 assert(processedImages[i]);
1189 }
1190
1191
1192 free(processedImages);
1193
1194 free(remappedSourceHistogram);
1195
1196 free(accumToCorrectHistogram);
1197
1198 free(accumSourceHistogram);
1199
1200 return 1;
1201
1202 }
1203
1204
1205
1206
1207 // Returns an array of (n * n-1) /2 daisies, with the histograms
1208
ReadHistograms(fullPath * fullPathImages,int numberImages)1209 histograms_struct *ReadHistograms (fullPath *fullPathImages, int numberImages)
1210 {
1211
1212 int value;
1213
1214 unsigned char *ptrCurrentLineBuffer;
1215 unsigned char *ptrCurrentPixelLineBuffer;
1216 histograms_struct *ptrHistograms; // << arrays of n * (n-1)/2 daisies
1217 int bytesPerLine;
1218 int bitsPerPixel;
1219 unsigned char *imagesDataBuffer; // numberOfImages * bytesPerLine
1220 int bytesPerPixel;
1221 int currentPixel;
1222 int currentRow;
1223 int otherImage;
1224 int currentImage;
1225 TIFF **ptrTIFFs = NULL;
1226 uint16 samplesPerPixel;
1227 uint16 bitsPerSample;
1228 uint32 imageLength;
1229 uint32 imageWidth;
1230 char tempString[512];
1231 char tempString2[512];
1232 int *ptrInt;
1233 histograms_struct * currentHistogram;
1234 histograms_struct * saveReturnValue;
1235
1236 CropInfo *crop_info = NULL;
1237 CropInfo *crop_info_array = NULL;
1238
1239 int temp;
1240
1241 int i;
1242
1243 int totalPixels = 0;
1244
1245 // esi = fullPathImages;
1246
1247 // (n * n-1)/2
1248
1249 if ( numberImages == 0 )
1250 return 0;
1251
1252 if ( ptQuietFlag == 0 ) {
1253
1254 Progress(_initProgress, "Reading Histograms");
1255
1256 }
1257
1258 saveReturnValue = ptrHistograms = calloc(numberImages * (numberImages-1)/2, sizeof(histograms_struct)); // Allocates one per every intersection n * (n-1) /2
1259
1260 if ( ptrHistograms == NULL )
1261 return 0;
1262
1263 ptrTIFFs = calloc(numberImages , sizeof(TIFF*));
1264 crop_info_array = (CropInfo *)calloc(numberImages, sizeof(CropInfo));
1265
1266 if ( ptrTIFFs == NULL || crop_info_array == NULL )
1267 {
1268 saveReturnValue = 0;
1269 goto Exit;
1270 }
1271
1272 currentImage = 0;
1273
1274 // Open TIFFs
1275 for (currentImage = 0; currentImage < numberImages ; currentImage ++ ) {
1276
1277 if (GetFullPath(&fullPathImages[currentImage],tempString) != 0) {
1278 PrintError("Could not get filename");
1279 saveReturnValue = 0;
1280 goto Exit;
1281 }
1282
1283 if ((ptrTIFFs[currentImage] = TIFFOpen(tempString, "r")) == NULL) {
1284 sprintf(tempString2, "Could not open TIFF file [%s]", tempString);
1285 PrintError(tempString2);
1286 saveReturnValue = 0;
1287 goto Exit;
1288 }
1289
1290 getCropInformationFromTiff(ptrTIFFs[currentImage], &(crop_info_array[currentImage]));
1291
1292 } // for loop
1293
1294 // Set defaults for all images (assumes they have all the same
1295 //TIFFGetField(ptrTIFFs[0], TIFFTAG_IMAGEWIDTH, &imageWidth);
1296 //TIFFGetField(ptrTIFFs[0], TIFFTAG_IMAGELENGTH, &imageLength);
1297 imageWidth = crop_info_array[0].full_width;
1298 imageLength = crop_info_array[0].full_height;
1299
1300 TIFFGetField(ptrTIFFs[0], TIFFTAG_BITSPERSAMPLE, &bitsPerSample);
1301 TIFFGetField(ptrTIFFs[0], TIFFTAG_SAMPLESPERPIXEL, &samplesPerPixel); // 0x11c
1302
1303 bitsPerPixel = samplesPerPixel * bitsPerSample;
1304
1305 bytesPerPixel = (bitsPerPixel + 7 ) / 8;
1306
1307 //bytesPerLine = TIFFScanlineSize(ptrTIFFs[0]);
1308 bytesPerLine = imageWidth * bytesPerPixel;
1309
1310 imagesDataBuffer = calloc(numberImages, bytesPerLine);
1311
1312 if ( imagesDataBuffer == 0 ) {
1313 PrintError("Not enough memory");
1314 saveReturnValue = 0;
1315 goto Exit;
1316 }
1317
1318 currentHistogram = ptrHistograms;
1319 currentImage = 0;
1320
1321 // Initializes the data structure
1322
1323 for (currentImage = 0; currentImage < numberImages; currentImage++) {
1324
1325 for (otherImage = currentImage + 1; otherImage < numberImages ; otherImage++, currentHistogram++) {
1326
1327 currentHistogram->overlappingPixels = 0;
1328
1329 currentHistogram->bytesPerSample = (bitsPerSample+7)/8;
1330
1331 currentHistogram->numberDifferentValues = 0x100;
1332
1333 currentHistogram->baseImageNumber = currentImage;
1334
1335 currentHistogram->otherImageNumber = otherImage;
1336
1337
1338 for (i = 0 ; i < 6 ; i++) {
1339
1340 if ((currentHistogram->ptrBaseHistograms[i] = calloc(currentHistogram->numberDifferentValues, sizeof(int))) == NULL)
1341 {
1342 saveReturnValue = 0;
1343 goto Exit;
1344 }
1345 if ((currentHistogram->ptrOtherHistograms[i] = calloc(currentHistogram->numberDifferentValues,sizeof(int))) == NULL)
1346 {
1347 saveReturnValue = 0;
1348 goto Exit;
1349 }
1350 } // for
1351
1352 } //for (otherImage = currentImage + 1; otherImage < numberImages ; otherImage++, currentHistogram++)
1353
1354 } // for (currentImage = 0; currentImage < numberImages; currentImage++) {
1355
1356 // fprintf(stderr,"Width %d Length %d BytesPerPixel %d per line%d\n", imageWidth, imageLength, bytesPerPixel, bytesPerLine);
1357
1358 for (currentRow = 0; currentRow < (int) imageLength; currentRow ++) {
1359
1360 if (currentRow * 2 == (int)(currentRow / 5.0) * 10) {
1361
1362 // This probably means every 2 percent update progress
1363
1364 sprintf(tempString, "%d", (int)(currentRow * 100/imageLength));
1365
1366 if ( ptQuietFlag == 0 ) {
1367 if (Progress(_setProgress, tempString) == 0) {
1368
1369 for (currentImage = 0 ; currentImage < numberImages ; currentImage++) {
1370
1371 TIFFClose(ptrTIFFs[currentImage]);
1372
1373 } //for
1374 saveReturnValue = 0;
1375 goto Exit;
1376
1377 } // progresss
1378
1379 } // quiet flag
1380
1381 }// if ( %ecx == %eax ) {
1382
1383
1384 ptrCurrentLineBuffer = imagesDataBuffer;
1385
1386 // Read the current row from each image
1387 for (currentImage=0; currentImage < numberImages; currentImage++) {
1388
1389 //fill buffer with empty space (zeros)
1390 memset(ptrCurrentLineBuffer, 0, bytesPerLine);
1391
1392 //need to handle cropped tiffs carefully here...
1393 crop_info = crop_info_array + currentImage;
1394
1395 if (currentRow>=crop_info->y_offset && currentRow < crop_info->y_offset + crop_info->cropped_height) {
1396 TIFFReadScanline(ptrTIFFs[currentImage], ptrCurrentLineBuffer + crop_info->x_offset * bytesPerPixel, currentRow - crop_info->y_offset, 0);
1397 }
1398
1399 RGBAtoARGB(ptrCurrentLineBuffer, imageWidth, bitsPerPixel);
1400
1401 ptrCurrentLineBuffer+= bytesPerLine;
1402
1403 } // while (currentImage < numberImages)
1404
1405 ptrCurrentPixelLineBuffer = imagesDataBuffer;
1406
1407 //for each pixel in the current line...
1408
1409 for (currentPixel = 0; currentPixel < (int)imageWidth; currentPixel++, ptrCurrentPixelLineBuffer+= bytesPerPixel ) {
1410
1411 unsigned char *ptrPixel;
1412 unsigned char *ptrOtherPixel;
1413
1414 // We process each currentHistogram
1415 currentHistogram = ptrHistograms;
1416
1417 currentImage = 0;
1418
1419 ptrPixel = ptrCurrentPixelLineBuffer;
1420 assert(ptrPixel < imagesDataBuffer + numberImages * bytesPerLine);
1421
1422 // for each pixel of each line of each image...
1423 for (currentImage = 0; currentImage < numberImages ; currentImage++, ptrPixel+= bytesPerLine) {
1424
1425 assert(ptrPixel < imagesDataBuffer + numberImages * bytesPerLine);
1426
1427 ptrOtherPixel = ptrPixel + bytesPerLine;
1428
1429 // for each pixel of each line of each current image, and images 'above' the current
1430
1431 for (otherImage = currentImage + 1; otherImage < numberImages; otherImage++, ptrOtherPixel += bytesPerLine) {
1432
1433 assert(ptrOtherPixel < imagesDataBuffer + numberImages * bytesPerLine);
1434
1435 assert(ptrPixel < ptrOtherPixel);
1436 assert(((int)(ptrOtherPixel - ptrPixel)) % bytesPerLine == 0);
1437
1438 /* Only process if the alpha channel is == 0xff in both pixels*/
1439
1440 if (0xff == ptrPixel[0] && 0xff == ptrOtherPixel[0] ) {
1441
1442 totalPixels ++;
1443 currentHistogram->overlappingPixels++;
1444
1445 // esi == ptrPixel
1446 // ebx == ptrOtherPixel
1447
1448
1449 // This seems to record the frequency of every pixels of every channel one of 3 channels of every image
1450 for (i = 0; i < 3; i++) {
1451
1452 // First byte is alpha channel
1453
1454 value = (unsigned char)ptrPixel[i+1];
1455 assert(value >= 0 && value <= 0xff);
1456
1457 ptrInt = currentHistogram->ptrBaseHistograms[i];
1458
1459 ptrInt[value] ++;
1460
1461 value = (unsigned char)ptrOtherPixel[i + 1];
1462 assert(value >= 0 && value <= 0xff);
1463
1464 ptrInt = currentHistogram->ptrOtherHistograms[i] ; // eax = ptrInt
1465
1466 ptrInt[value] ++;
1467
1468 }
1469
1470 // compute the other 6 histograms
1471 temp = panoColourComputeIntensity(ptrPixel[1], ptrPixel[2], ptrPixel[3]);
1472 assert(temp >= 0 && temp <= 0xff);
1473 ptrInt = currentHistogram->ptrBaseHistograms[IDX_INTENSITY];
1474 ptrInt[temp] ++;
1475
1476 temp = panoColourComputeIntensity(ptrOtherPixel[1], ptrOtherPixel[2], ptrOtherPixel[3]);
1477 assert(temp >= 0 && temp <= 0xff);
1478 ptrInt = currentHistogram->ptrOtherHistograms[IDX_INTENSITY];
1479 ptrInt[temp] ++;
1480
1481 temp = panoColourComputeSaturation(ptrPixel[1], ptrPixel[2], ptrPixel[3]);
1482 assert(temp >= 0 && temp <= 0xff);
1483 ptrInt = currentHistogram->ptrBaseHistograms[IDX_SATURATION];
1484 ptrInt[temp] ++;
1485
1486 temp = panoColourComputeSaturation(ptrOtherPixel[1], ptrOtherPixel[2], ptrOtherPixel[3]);
1487 assert(temp >= 0 && temp <= 0xff);
1488 ptrInt = currentHistogram->ptrOtherHistograms[IDX_SATURATION];
1489 ptrInt[temp] ++;
1490
1491 temp = panoColourComputeHue(ptrPixel[1], ptrPixel[2], ptrPixel[3]);
1492 assert(temp >= 0 && temp <= 0xff);
1493 ptrInt = currentHistogram->ptrBaseHistograms[IDX_HUE];
1494 ptrInt[temp] ++;
1495
1496 temp = panoColourComputeHue(ptrOtherPixel[1], ptrOtherPixel[2], ptrOtherPixel[3]);
1497 assert(temp >= 0 && temp <= 0xff);
1498 ptrInt = currentHistogram->ptrOtherHistograms[IDX_HUE];
1499 ptrInt[temp] ++;
1500
1501 } // if ( 0 != *ptrPixel and 0 != *ptrOtherPixel ) {
1502
1503
1504 currentHistogram++; //edi = edi + 0x44; //68 ?? again, this should be edi ++;
1505
1506 } // for (otherImage = currentImage + 1; otherImage < numberImages; otherImage++ ) {
1507
1508 }// for (currentImage = 0; currentImage >=numberImages ;...
1509
1510 } // for (currentPixel = 0; currentPixel < imageWidth; currentPixel++... ) {
1511
1512 } //for (currentRow = 0; currentRow < imageLength; currentRow ++) {
1513
1514
1515 // for (i = 0; i< numberImages * (numberImages-1)/2; i++) {
1516 // fprintf(stderr, "Histogram %d Images %d %d, %d Pixels\n", i ,
1517 // ptrHistograms[i].baseImageNumber,
1518 // ptrHistograms[i].otherImageNumber,
1519 // ptrHistograms[i].overlappingPixels);
1520 // totalPixels2 += ptrHistograms[i].overlappingPixels;
1521 // }
1522
1523 // assert(totalPixels2 == totalPixels);
1524
1525 for (currentImage=0; currentImage < numberImages; currentImage++) {
1526
1527 TIFFClose(ptrTIFFs[currentImage]);
1528
1529 }
1530
1531 Exit:
1532 free(ptrTIFFs);
1533 free(imagesDataBuffer);
1534 free(crop_info_array);
1535
1536 return(saveReturnValue);
1537 }
1538
1539
MinDouble(double a,double b,double c)1540 double MinDouble(double a, double b, double c)
1541 {
1542 double min = 0;
1543 if (a < b) {
1544 min = a;
1545 } else {
1546 min = b;
1547 }
1548 if (c < min)
1549 return c;
1550 else
1551 return min;
1552 }
1553
MaxDouble(double a,double b,double c)1554 double MaxDouble(double a, double b, double c)
1555 {
1556 double max = 0;
1557 if (a > b) {
1558 max = a;
1559 } else {
1560 max = b;
1561 }
1562 if (c > max)
1563 return c;
1564 else
1565 return max;
1566 }
1567
panoColourRGBtoHSV(int R,int G,int B,double * pH,double * pS,double * pV)1568 void panoColourRGBtoHSV(int R, int G, int B, double *pH, double *pS, double *pV)
1569 {
1570 double r = ( R / 255.0 );
1571 double g = ( G / 255.0 );
1572 double b = ( B / 255.0 );
1573 double V, H,S;
1574
1575 double min, max, delta;
1576
1577 min = MinDouble( r, g, b );
1578 max = MaxDouble( r, g, b );
1579 V = max; // v
1580
1581 delta = max - min;
1582
1583
1584 if( max != 0 )
1585 S = delta / max; // s
1586 else {
1587 // r = g = b = 0 // s = 0, v is undefined
1588 S = 0;
1589 H = 0;
1590
1591 *pH = H;
1592 *pS = S;
1593 *pV = V;
1594 return;
1595 }
1596
1597 if (delta == 0) {
1598 H = 0;
1599 } else {
1600
1601 if( r == max )
1602 H = ( g - b ) / delta; // between yellow & magenta
1603 else if( g == max )
1604 H = 2 + ( b - r ) / delta; // between cyan & yellow
1605 else
1606 H = 4 + ( r - g ) / delta; // between magenta & cyan
1607
1608 H *= 60; // degrees
1609 if( H < 0 )
1610 H += 360;
1611 }
1612 *pH = H;
1613 *pS = S;
1614 *pV = V;
1615 }
1616
panoColourHSVtoRGB(double H,double S,double V,int * pR,int * pG,int * pB)1617 void panoColourHSVtoRGB(double H, double S, double V, int *pR, int *pG, int *pB)
1618 {
1619
1620 int i;
1621 double f, p, q, t;
1622 double R, G, B;
1623
1624 if( fabs(S) < 0.000001 ) {
1625 // achromatic (grey)
1626 *pR = (int)(*pG = *pB = V * 255);
1627 return;
1628 }
1629
1630 H /= 60; // sector 0 to 5
1631 i = (int)H;
1632
1633
1634 f = H - i; // factorial part of h
1635 p = V * ( 1 - S );
1636 q = V * ( 1 - S * f );
1637 t = V * ( 1 - S * ( 1 - f ) );
1638
1639 switch( i ) {
1640 case 0:
1641 R = V;
1642 G = t;
1643 B = p;
1644 break;
1645 case 1:
1646 R = q;
1647 G = V;
1648 B = p;
1649 break;
1650 case 2:
1651 R = p;
1652 G = V;
1653 B = t;
1654 break;
1655 case 3:
1656 R = p;
1657 G = q;
1658 B = V;
1659 break;
1660 case 4:
1661 R = t;
1662 G = p;
1663 B = V;
1664 break;
1665 default: // case 5:
1666 R = V;
1667 G = p;
1668 B = q;
1669 break;
1670 }
1671 *pR = (int)(R * 255);
1672 *pG = (int)(G * 255);
1673 *pB = (int)(B * 255);
1674 }
1675
1676
1677
1678
panoColourComputeHue(unsigned char red,unsigned char green,unsigned char blue)1679 unsigned char panoColourComputeHue (unsigned char red, unsigned char green, unsigned char blue)
1680 {
1681 // Compute Hue
1682 double H, S, V;
1683 int h;
1684
1685 panoColourRGBtoHSV(red, green, blue, &H, &S, &V);
1686
1687 // Normalize to 0..255 because the mapping tables are configured that way...
1688
1689 H = H * (255.0/ 360);
1690 h = (int)H;
1691 assert(h >=0 && h <= 255);
1692
1693 return h;
1694 }
1695
1696
panoColourComputeIntensity(unsigned char red,unsigned char green,unsigned char blue)1697 unsigned char panoColourComputeIntensity (unsigned char red, unsigned char green, unsigned char blue)
1698 {
1699
1700 unsigned temp = (red + green + blue)/3;
1701
1702 assert(temp >= 0 && temp <= 255);
1703
1704 return temp;
1705 }
1706
panoColourComputeSaturation(unsigned char red,unsigned char green,unsigned char blue)1707 unsigned char panoColourComputeSaturation (unsigned char red, unsigned char green, unsigned char blue)
1708 {
1709 double H, S, V;
1710 int h;
1711
1712 panoColourRGBtoHSV(red, green, blue, &H, &S, &V);
1713
1714 // Normalize to 0..255 because the mapping tables are configured that way...
1715
1716 S = S * 255.0;
1717 h = (int)S;
1718 assert(h >=0 && h <= 255);
1719
1720 return h;
1721 }
1722
1723
1724
1725 //int ComputeColourBrightnessCorrection(calla_struct *calla)
1726 //{
1727 //
1728 //
1729 // // var32 0xffffffe0(%ebp) var32
1730 // // var28 0xffffffe4(%ebp) double *var28
1731 // // var24 0xffffffe8(%ebp) var24
1732 // // var20 0xffffffec(%ebp) int *correctedImages
1733 // // var16 0xfffffff0(%ebp) currentImageNumber
1734 // // var12 0xfffffff4(%ebp) channel
1735 // // var08 0xfffffff8(%ebp) int numberIntersections
1736 //
1737 //
1738 // numberIntersections = (calla->numberImages - 1) * calla->numberImages;
1739 //
1740 //
1741 // correctedImages = calloc(calla->numberImages, sizeof(int));
1742 //
1743 //
1744 // edi = malloc(0x100 * sizeof(double));
1745 // var24 = malloc(0x100 * sizeof(double));
1746 // var28 = malloc(0x100 * sizeof(double));
1747 //
1748 //
1749 // if ( correctedImages == 0 )
1750 // return 0;
1751 //
1752 //
1753 // if ( var24 == 0 )
1754 // return 0;
1755 //
1756 //
1757 // if ( edi == 0 )
1758 // return 0;
1759 //
1760 //
1761 //
1762 // eax = calla->indexReferenceImage;
1763 //
1764 //
1765 // correctedImages[calla->indexReferenceImage] = 1;
1766 //
1767 // while ((currentImageNumber = Unknown43(correctedImages, calla)) != -1) {
1768 //
1769 // assert(currentImageNumber > 0);
1770 // assert(currentImageNumber < calla->numberImages);
1771 // assert(correctedImages[currentImageNumber] == 0);
1772 //
1773 // channel = 0;
1774 //
1775 // for (channel = 0; var < 6; var++) {
1776 //
1777 //
1778 // for (ecx =0; ecx < 0xff; ecx ++) {
1779 //
1780 // edi[ecx] = 0;
1781 // var24[ecx] = 0;
1782 //
1783 // }
1784 //
1785 //
1786 // // for each daisy records (histograms)
1787 // for (esi = 0; esi < numberIntersections ; esi++) {
1788 //
1789 // currentHistogram = &calla->ptrHistograms[esi];
1790 //
1791 //
1792 // if ( currentHistogram->overlappingPixels > 1000 ) { // it is not consistent XXX
1793 //
1794 // if ( currentHistogram->baseImageNumber == currentImageNumber &&
1795 //
1796 // correctedImages[currentHistogram->otherImageNumber] != 0 ) {
1797 //
1798 //
1799 // // Do something with the already corrected histogram
1800 // Unknown37(currentHistogram->ptrOtherHistograms[channel],
1801 // var28,
1802 // calla->magnolia[currentHistogram->otherImageNumber], channel);
1803 //
1804 //
1805 // for (ecx = 0; ( %ecx <= $0xff ); ecx ++) {
1806 //
1807 // edi[ecx] += var28[ecx];
1808 //
1809 // }
1810 //
1811 // // Now process the current histogram
1812 //
1813 // ptrHistogram = calla->ptrHistograms[esi].ptrBaseHistograms;
1814 //
1815 // for (ecx = 0; ecx <= 0xff; ecx ++) {
1816 //
1817 // array = ptrHistogram[channel];
1818 // var24[ecx] += array[ecx];
1819 //
1820 // }
1821 //
1822 // } else { // if ( 0xc(%ecx, edx, 1) == %ebx ) {
1823 //
1824 //
1825 // currentHistogram = &calla->ptrHistograms[esi];
1826 //
1827 //
1828 // if (currentHistogram->otherImageNumber == currentImageNumber ) {
1829 //
1830 // if ( correctedImages[currentHistogram->baseImageNumber] != 0 ) {
1831 //
1832 // Unknown37(currentHistogram->ptrBaseHistograms[channel],
1833 // var28,
1834 // calla->magnolia[currentHistogram->baseImageNumber], channel);
1835 //
1836 // for (ecx = 0; ecx <= 0xff; ecx ++) {
1837 //
1838 // edi[ecx] += var28[ecx];
1839 //
1840 //
1841 // } // for (ecx = 0; ecx <= 0xff; ecx ++) {
1842 //
1843 //
1844 // currentHistogram = calla->ptrHistograms[esi];
1845 //
1846 // ptrHistogram = currentHistogram->ptrOtherHistogram;
1847 //
1848 //
1849 // for (ecx = 0; ecx <= 0xff; ecx ++) {
1850 //
1851 // array = ptrHistogram[channel];
1852 //
1853 // var24[ecx] += array[ecx];
1854 //
1855 //
1856 // } // for (ecx = 0; ecx <= 0xff; ecx ++) {
1857 //
1858 // } // if ( (%ebx,eax,4) != 0 ) {
1859 //
1860 // } // if ( 0x10(%ecx,edx,1) == %ebx ) {
1861 // } // end of else
1862 // } //if ( (%ecx, edx, 1) > $0x3e8 ) {
1863 //
1864 // } // for (esi = 0; %esi < numberIntersections ; esi++) {
1865 //
1866 // // This should be responsible for updating the current correction table
1867 // // for the currentimage
1868 //
1869 // Unknown41(var24,
1870 // edi,
1871 // (calla->magnolia[currentImageNumber].fieldx04)[channel]);
1872 //
1873 //
1874 //
1875 // } // for (channel = 0; var < 6; var++) {
1876 //
1877 //
1878 // correctedImages[currentImageNumber] = 1;
1879 //
1880 // } // while (Unknown43 (...) != -1);
1881 //
1882 // for (i =0 ; i< calla->numberImages; i++) {
1883 // // are all the images opimized?
1884 // assert(correctedImages[i]);
1885 // }
1886 //
1887 // free(var20intar);
1888 //
1889 // free(var28);
1890 //
1891 // free(var24);
1892 //
1893 // free(edi);
1894 //
1895 // return 0;
1896 //
1897 //
1898 //}
1899 //
1900 //
1901
AssertEqualCurves(double * first,double * second,int size)1902 void AssertEqualCurves(double *first, double *second, int size)
1903 {
1904 int i;
1905 for (i=0;i<size;i++) {
1906 assert(first[i] == second[i]);
1907 }
1908 }
1909
1910
1911
1912
CorrectImageColourBrigthness(Image * image,magnolia_struct * magnolia,int parm3)1913 void CorrectImageColourBrigthness(Image *image, magnolia_struct *magnolia, int parm3)
1914 {
1915
1916 int currentRow;
1917 int currentPixel;
1918 unsigned char *pixel;
1919 int edi;
1920 double * (mappingCurves[6]);
1921 unsigned char *ptrPixel;
1922 int channel;
1923 int level;
1924
1925 for (channel = 0; channel < 6; channel++ ) {
1926
1927 if ((mappingCurves[channel] = calloc(0x100 , sizeof(double))) == NULL) {
1928 PrintError("Not enough memory\n");
1929 return ;
1930 }
1931
1932 }
1933
1934
1935 // fprintf(stderr, "Colour correcting image\n");
1936
1937 edi = 0;
1938
1939 for (level = 0; level < 0x100; level++) {
1940
1941 for (channel = 0; channel< 6; channel ++ ) {
1942 double *ptr;
1943 // takes an array of 0x100 doubles, a double between 0 and 0x100, and number of doubles
1944 // Remaps the correction functions
1945
1946 ptr = mappingCurves[channel];
1947
1948 ptr[level] = (*magnolia->function)(magnolia->fieldx04[channel], level, magnolia->components);
1949
1950 }
1951
1952 } // for (edi = 0; edi < 0x100; edi++) {
1953
1954 pixel = *image->data;
1955
1956 switch (parm3) {
1957
1958 case 1:
1959
1960 printf("**************************Bright\n");
1961
1962 currentRow = 0;
1963
1964 while ( currentRow < image->height) {
1965
1966 currentPixel = 0;
1967
1968 ptrPixel = pixel;
1969
1970 while ( currentPixel < image->width ) {
1971
1972
1973 if (*ptrPixel != 0 ) {
1974
1975 int R, G, B;
1976 double H, S, I;
1977
1978 R = ptrPixel[1];
1979 G = ptrPixel[2];
1980 B = ptrPixel[3];
1981
1982 panoColourRGBtoHSV(R, G, B, &H, &S, &I);
1983
1984 assert(H >= 0.0 && H<=360.0);
1985 assert(S >= 0.0 && S<=1.0);
1986 assert(I >= 0.0 && I<=1.0);
1987
1988 // Now remap the intensity
1989
1990
1991 I = RemapDouble(I * 255.0, mappingCurves[IDX_INTENSITY]) / 255.0;
1992
1993 panoColourHSVtoRGB(H, S, I, &R, &G, &B);
1994
1995 if (R < 0 || R > 255 || G < 0 || G > 255 || B < 0 || B > 255) {
1996 printf("Value of R G B %d %d %d\n", R, G, B);
1997 assert(0);
1998 }
1999 //assert(R >= 0 && R <= 255);
2000 //assert(G >= 0 && G <= 255);
2001 //assert(B >= 0 && B <= 255);
2002
2003 ptrPixel[1] = R;
2004 ptrPixel[2] = G;
2005 ptrPixel[3] = B;
2006
2007
2008 #if 0
2009
2010
2011 ecx = ptrPixel[1] + edx;
2012
2013 if (ecx < 0) {
2014 ptrPixel[1] = 0;
2015 } else if (ecx > 0xff) {
2016 ptrPixel[1] = 0xff;
2017 } else {
2018 ptrPixel[1] = ecx;
2019 }
2020
2021 ecx = ptrPixel[2] + edx;
2022
2023 if (ecx < 0) {
2024 ptrPixel[2] = 0;
2025 } else if (ecx > 0xff) {
2026 ptrPixel[2] = 0xff;
2027 } else {
2028 ptrPixel[2] = ecx;
2029 }
2030
2031 ecx = ptrPixel[3] + edx;
2032
2033 if (ecx < 0) {
2034 ptrPixel[3] = 0;
2035 } else if (ecx > 0xff) {
2036 ptrPixel[3] = 0xff;
2037 } else {
2038 ptrPixel[3] = ecx;
2039 }
2040 #endif
2041
2042 }
2043
2044 currentPixel++;
2045 ptrPixel+=4; // advance an entire pixel
2046
2047 } // while ( currentPixel < image->width )
2048
2049 currentRow++;
2050
2051 pixel += image->bytesPerLine;
2052
2053 } //while ( currentRowx < image->height) {
2054 break;
2055
2056
2057 case 0: //case of switch
2058
2059 // Normal colour correction: hue + intensity
2060
2061 ptrPixel = *(image->data);
2062
2063 for (currentRow = 0; currentRow < image->height; currentRow++) {
2064
2065 for (currentPixel = 0 ; currentPixel < image->width ; currentPixel++) {
2066
2067 if ( ptrPixel[0] != 0 ) {
2068
2069 // If we have data... do each of the channels
2070
2071
2072 for (channel = 0; channel < 3 ; channel ++) {
2073
2074 ptrPixel[channel+1] = RemapPoint(ptrPixel[channel+1], mappingCurves[channel]);
2075
2076 }
2077
2078 } // if ( (ptrPixel) != 0 ) {
2079
2080 ptrPixel +=4;
2081
2082 } // for (currentPixel = 0 ; currentPixel < imageWidth ; currentPixel++) {
2083 }
2084 break;
2085
2086 case 2: // case of switch parm3
2087
2088 // Colour-only correction: saturation
2089
2090 currentRow = 0;
2091 ptrPixel = *(image->data);
2092
2093 for (currentRow = 0; currentRow < image->height; currentRow++) {
2094
2095 for (currentPixel = 0; currentPixel < image->width; currentPixel++) {
2096
2097 if (ptrPixel[0] != 0) {
2098
2099 int R, G, B;
2100 double H, S, I;
2101
2102 R = ptrPixel[1];
2103 G = ptrPixel[2];
2104 B = ptrPixel[3];
2105
2106 panoColourRGBtoHSV(R, G, B, &H, &S, &I);
2107
2108 assert(H >= 0.0 && H<=360.0);
2109 assert(S >= 0.0 && S<=1.0);
2110 assert(I >= 0.0 && I<=1.0);
2111
2112 // Now remap the intensity
2113
2114 // First normalize it to 255
2115 H /= 360.0;
2116 H *= 255.0;
2117
2118 // I find that changing hue is not a very good idea. it tends to create a very strong colour cast
2119
2120 //H = RemapDouble(H, mappingCurves[IDX_HUE]);
2121
2122 H /= 255.0;
2123 H *= 360.0;
2124
2125 S = RemapPoint(S * 255.0, mappingCurves[IDX_SATURATION]) / 255.0;
2126
2127 assert(S >= 0.0 && S <= 1.0);
2128 assert(H >= 0.0 && S <= 360);
2129
2130 panoColourHSVtoRGB(H, S, I, &R, &G, &B);
2131
2132 assert(R >= 0 && R <= 255);
2133 assert(G >= 0 && G <= 255);
2134 assert(B >= 0 && B <= 255);
2135
2136 ptrPixel[1] = R;
2137 ptrPixel[2] = G;
2138 ptrPixel[3] = B;
2139
2140 #if 0
2141
2142 // THIS IS THE OLD CODE
2143
2144 int ebx,ecx;
2145
2146 int ebx, var49, var56, esi;
2147
2148 ebx = panoColourComputeIntensity(ptrPixel[1], ptrPixel[2], ptrPixel[3]);
2149
2150 var49 = RemapPoint( panoColourComputeSaturation(ptrPixel[1], ptrPixel[2], ptrPixel[3]),
2151 mappingCurves[IDX_SATURATION]); //var08
2152
2153 var56 = RemapPoint(panoColourComputeHue(ptrPixel[1], ptrPixel[2], ptrPixel[3]),
2154 mappingCurves[IDX_HUE]); // var04
2155 esi = var49;
2156
2157 ptrPixel[1] = Unknown47(ebx, var49, var56);
2158
2159 ptrPixel[2] = Unknown48(ebx, var49, var56);
2160
2161 ptrPixel[3] = Unknown49(ebx, var49, var56);
2162
2163 #endif
2164
2165 } // if (ptrPixel[0] != 0) {
2166
2167 ptrPixel += 4;
2168
2169 } // for (currentPixel = 0; currentPixel < image->width; currentPixel++) {
2170
2171 } // for (currentRow = 0; currentRow < image->height; currentRow++) {
2172
2173 break;
2174
2175 }
2176
2177
2178 for (edi = 0;edi < 6; edi++) {
2179
2180 free(mappingCurves[edi]);
2181
2182 }
2183
2184
2185 }//end of function
2186
2187
2188
MapFunction(double parm[],double x,int n)2189 double MapFunction(double parm[], double x, int n)
2190 {
2191 double x_1;
2192 int e;
2193 double result;
2194 /*
2195
2196 Assume we have a curve defined from [0 to n-1] but it is
2197 discretized. So we actually have an array p[] of n elements that
2198 defines this curve (called it p). The curve grows monotonically,
2199 from 0 to n-1.
2200
2201
2202 We also have the same curve, but "stretched" from [0 to 255]. Call it
2203 f. We have a point x in this range. Now we need to find f(x)
2204
2205 So we need to find the corresponding x_1 in the domain of p.
2206
2207 But x_1 might fall in between 2 different integers e and e+1. We
2208 assume the cuve is a straight line between these two points:
2209
2210 e = floor(x_1);
2211
2212 The result y is:
2213
2214 y = p(x_1) = ((x_1 - e ) * (p[e+1] - p[e])) + p[e]
2215
2216 */
2217
2218 x_1 = (x * 255.0) / ( n - 1 );
2219
2220 e = floor(x_1);
2221
2222 if ( e < 0 ) {
2223
2224 result = parm[0];
2225
2226 } else if ( e < n - 1 ) {
2227 assert(e < n);
2228 assert(e >= 0);
2229
2230 result = ((x_1 - e) * (parm[e+1] - parm[e])) + parm[e];
2231
2232 assert( result >= parm[e]);
2233
2234 } else {
2235
2236 result = parm[n-1];
2237 }
2238
2239
2240 if (result >= 0x100) {
2241 // THIS CODE IS JUST FOR DEBUGGING PURPOSES
2242
2243 int i;
2244 fprintf(stderr, "Result %g Value %d Array: ", result, n);
2245
2246 for (i=0; i<= 0xff; i++) {
2247 fprintf(stderr, "%d: %g ", i, parm[i]);
2248 }
2249 fprintf(stderr, "\n");
2250 assert(0);
2251 }
2252
2253
2254 return result;
2255
2256
2257 // should return a double
2258
2259 } // end of function Unknown33
2260
2261
2262
2263
2264
2265
Unknown48(unsigned char parm0,unsigned char parm1,unsigned char parm2)2266 unsigned char Unknown48(unsigned char parm0, unsigned char parm1, unsigned char parm2)
2267 {
2268 // (3a - 4b + 512 + 2c - 256) * (2/3)
2269 // => (3a - 4b + 2c + 256) * (2/3)
2270 // if a = b = c
2271 // => (3a -4a + 2a + 256) * (2/3 => (a +256) * 2/3
2272
2273
2274 int ecx;
2275
2276 ecx = parm0 * 3 - (parm1 * 4 - 512);
2277
2278 ecx = (ecx + 2 * parm2 - 256);
2279
2280 ecx = (ecx * 2)/3;
2281
2282 if (ecx < 0) {
2283 return 0;
2284 } else {
2285 if (ecx > 0xff)
2286 return 0xff;
2287 else
2288 return ecx;
2289 }
2290 }
2291
2292
2293
Unknown49(unsigned char parm0,unsigned char parm1,unsigned char parm2)2294 unsigned char Unknown49(unsigned char parm0, unsigned char parm1, unsigned char parm2)
2295 {
2296 return Unknown48(parm0, parm1, parm2) ;
2297 }
2298
2299
2300
RemapHistogram(int * histogram,double * remappedHistogram,magnolia_struct * magnolia,int channel)2301 void RemapHistogram(int *histogram, double *remappedHistogram, magnolia_struct *magnolia, int channel)
2302 {
2303
2304 int index;
2305 double mappingFunction[256];
2306
2307 double prevValue;
2308 double nextValue;
2309
2310 int value;
2311
2312 double delta;
2313
2314 double sumR, sumH;
2315
2316 double contribution;
2317
2318 // fprintf(stderr, "Doubles remappedHistogram: \n");
2319
2320 // Compute the mapping function.
2321
2322 // fprintf(stderr, "Doubles array: \n");
2323 for (index = 0; index < 0x100 ; index++) {
2324
2325 mappingFunction[index] = (*magnolia->function)(magnolia->fieldx04[channel], (double)index, magnolia->components);
2326
2327 // DEBUGGING code
2328 if ((int)mappingFunction[index] < 0 || (int)mappingFunction[index] > 0xff) {
2329 fprintf(stderr, "error %d %g\n", index, mappingFunction[index]);
2330 assert(0);
2331 }
2332
2333 } //
2334
2335 // Initialize the remapped histogram
2336 for (index = 0; index < 0x100; index ++) {
2337
2338 remappedHistogram[index] = 0;
2339
2340 } //if ( %index <= $0xff )
2341
2342 index = 0;
2343
2344 sumR = 0;
2345 sumH = 0;
2346 for (index = 0; index < 0x100; index++) {
2347 sumH += histogram[index];
2348 }
2349 if (sumR != sumH) {
2350 // printf("Starting*** Sum in histogram: %f\n", sumH);
2351 }
2352
2353
2354 prevValue = 0.0;
2355 nextValue = 0.0;
2356
2357 for (index = 0; index < 0x100; index++) {
2358
2359 {
2360 int i;
2361 double sumR = 0;
2362 double sumH = 0;
2363
2364 for (i = 0; i < index; i++) {
2365 sumH += histogram[i];
2366 }
2367 for (i = 0; i < 0x100; i++) {
2368 sumR += remappedHistogram[i];
2369 }
2370 if (fabs(sumR - sumH) > 0.00001) {
2371 printf("****B********** Sum in histograms: %d R %f H %f, difference %g\n", index, sumR, sumH, sumH-sumR);
2372 assert(0);
2373 }
2374 // printf("******* Sum in histograms: %d remapped %f original %f\n", index, sumR, sumH);
2375 }
2376
2377 if (index == 0 ) {
2378
2379 prevValue = mappingFunction[1] - 2 * mappingFunction[0];
2380 // assert((mappingFunction[1] - 2 * mappingFunction[0]) >= 0);
2381
2382 } else { // if ( %index == 0 ) {
2383
2384 prevValue = mappingFunction[index - 1]; // makes sense, it would be undefined for index == 0
2385
2386 }// if ( %index == 0 ) {
2387
2388 if ( index == 0xff ) {
2389 // Extrapolate another value, a[xff] + delta (a[xff] - a[xff-1])
2390 nextValue = 2 * mappingFunction[0xff] - mappingFunction[0xff-1];
2391
2392 } else { // if ( %index == $0xff ) {
2393
2394 nextValue = mappingFunction[index + 1];
2395
2396 } // if ( %index == $0xff ) {
2397
2398 // ***************************************************
2399
2400 if ( (int)mappingFunction[index] == 0xff ) {
2401 // REPEATED FROM AAAAAAA
2402 remappedHistogram[255] = histogram[index] + remappedHistogram[255];
2403 continue;
2404
2405 } // if ( (int) mappingFunction[index] == $0xff ) {
2406
2407 if (fabs(nextValue - prevValue) <= 2.0) { // remember to negate as part of the if jump less
2408
2409 // RESET stack
2410 // remove 2 values from the stack
2411
2412 // fprintf(stderr, "PTdouble[%d] = %g\n", index, mappingFunction[index]);
2413 assert(mappingFunction[index] >= 0 && mappingFunction[index] <= 0xff);
2414
2415
2416 } else { // if (top stack (and pop) ?? 2.0) { // remember to negate as part of the if jump less
2417
2418
2419 int ecx;
2420 int var2072;
2421 int edx;
2422
2423 //////////////////////////////////////////////////////////
2424 double st_0;
2425
2426 double checkSum = histogram[index];
2427
2428 double deltaPrev;
2429 double deltaNext;
2430
2431 ecx = (int)nextValue;
2432
2433 if ((int)nextValue > 0xff ) {
2434 ecx = 0xff;
2435 }
2436
2437 edx = (int)prevValue;//TOP
2438 if (edx < prevValue) {
2439 edx ++;
2440 } //if
2441 assert(edx == ceil(prevValue));
2442
2443 if ( edx < 0 ) {
2444 edx = 0;
2445 } // if ( %edx < 0 ) {
2446
2447
2448 st_0 = 0;
2449
2450 deltaPrev = (mappingFunction[index] - prevValue);
2451 deltaNext = (nextValue - mappingFunction[index]);
2452
2453 for (var2072 = edx ; var2072 <= ecx; var2072++ ) {
2454
2455 if (var2072 < mappingFunction[index]) { // I AM IN DOUBT of this one
2456
2457 if ( deltaPrev != 0.0 ) {
2458 assert(mappingFunction[index] - prevValue > 0);
2459 //
2460 st_0 += (var2072 - prevValue)/deltaPrev;
2461
2462 }
2463 } else { // if (var2072 < mappingFunction[index])
2464
2465 if ( 0.0 != deltaNext ) {
2466 assert(nextValue - mappingFunction[index] > 0);
2467 st_0 += (nextValue- var2072) / deltaNext;
2468 }
2469 } // if (var2072 < mappingFunction[index])
2470 } // for loop
2471
2472 assert(st_0 != 0.0);
2473 if (0.0 != st_0 ) {
2474
2475 for (var2072 = edx; var2072 <= ecx; var2072++) {
2476
2477 if (var2072 < mappingFunction[index]) {
2478
2479 if (deltaPrev != 0.0) {
2480
2481 double contribution = (histogram[index] * (var2072 - prevValue)) / (deltaPrev * st_0);
2482
2483 remappedHistogram[var2072] += contribution;
2484
2485 checkSum -= contribution;
2486
2487 }
2488
2489 } else {
2490
2491 if (deltaNext != 0.0) {
2492
2493 double contribution = (histogram[index] * (nextValue - var2072) ) / (deltaNext * st_0);
2494
2495 //((nextValue - var2072) / deltaNext) * histogram[index])/st_0;
2496
2497 // ????????????? it was remappedHistogram[edi]
2498 remappedHistogram[var2072] += contribution;
2499
2500 checkSum -= contribution;
2501
2502 } //if
2503 } //if
2504 } // for (var2072 = edx; var2072 <= ecx; var2072++) {
2505
2506
2507 if (checkSum > 0) {
2508 // printf("Missing------------------------->%f\n", checkSum);
2509 remappedHistogram[index] += checkSum; // perhaps this should be spreaded over several points XXXXXXXXXXXXX
2510 checkSum = 0;
2511 }
2512
2513 continue; //???
2514 }
2515
2516 assert(0);
2517 assert(checkSum == 0.0);
2518 assert(st_0 == 0);
2519
2520 //////////////////////////////////////////////////////////
2521
2522 } //// if (top stack (and pop) ?? 2.0) { // remember to negate as part of the if jump less
2523
2524
2525 value = (int) mappingFunction[index];
2526
2527 delta = mappingFunction[index] - (int) mappingFunction[index];
2528
2529 // delta determines the rounding error. 1-delta*histogram gets
2530 // added to the current value and
2531 // delta *histogram to the next value.
2532
2533 assert(value < 255);
2534
2535 contribution = (1 - delta) *histogram[index];
2536 // contribution2 = delta * histogram[index];
2537
2538 remappedHistogram[value] += contribution;
2539 remappedHistogram[value+1] += histogram[index] - contribution;
2540
2541
2542
2543 } // for (index = 0; index < 0x100; index++) {
2544
2545 {
2546 double sumR = 0, sumH = 0;
2547 for (index = 0; index < 0x100; index++) {
2548 sumR += remappedHistogram[index];
2549 sumH += histogram[index];
2550 }
2551 if (fabs(sumR - sumH) > 0.00001) {
2552 printf("F************* Sum in histograms: %f %f\n", sumR, sumH);
2553 assert(0);
2554 }
2555 }
2556
2557
2558 }
2559
2560
ComputeAdjustmentCurve(double * sourceHistogram,double * referenceHistogram,double * curve)2561 void ComputeAdjustmentCurve(double *sourceHistogram, double *referenceHistogram, double *curve)
2562 {
2563 // Unknown41
2564
2565 double copyReferenceHistogram[0x100];
2566 double copySourceHist[0x100];
2567 double otherArray[0x100];
2568 int i,j;
2569
2570 double sum2 = 0.0;
2571 double sum = 0.0;
2572 // double total = 0;
2573
2574 // FIRST, make sure we parameters are sound
2575
2576 // Should the two histograms have the same accumated total?
2577
2578 double total1 = 0;
2579 double total2 = 0;
2580 for (i=0;i<0x100;i++) {
2581 if (sourceHistogram[i] < 0) {
2582
2583 printf("I am going to crash %f\n", sourceHistogram[i]);
2584 }
2585 if (referenceHistogram[i] < 0) {
2586 int j;
2587 for (j=0;j<0x100;j++) {
2588 printf("I am going to crash %f ", referenceHistogram[j]);
2589 }
2590 printf("I am going to crash at i %d %f ", i, referenceHistogram[i]);
2591 printf("\n");
2592 }
2593 assert(sourceHistogram[i]>= 0);
2594 assert(referenceHistogram[i]>= 0);
2595 total1 += sourceHistogram[i];
2596 total2 += referenceHistogram[i];
2597
2598 }
2599 // if (total1 != total2) {
2600 // fprintf(stderr, "%g %g difference %g\n", total1, total2, total1 - total2);
2601 // }
2602
2603 // assert(total1 == total2);
2604
2605
2606 // FIRST COPY THE histograms to a temporal place.
2607
2608
2609 memcpy(copySourceHist, sourceHistogram, 0x200 * 4); // Check this one ???
2610 memcpy(copyReferenceHistogram, referenceHistogram, 0x200 * 4); // Check this one ???
2611
2612
2613 //// Debugging messages
2614
2615 // fprintf(stderr, "41->Histogram source: ");
2616 // for (i=0; i<0x100;i++) {
2617 // fprintf(stderr, " %d:%g ", i, copySourceHist[i]);
2618 // assert(copySourceHist[i] == sourceHistogram[i]);
2619 // total += copySourceHist[i];
2620 // }
2621 // fprintf(stderr, "\nTotal: %g\n", total);
2622 //
2623 // total = 0;
2624 // fprintf(stderr, "41->Histogram Reference : ");
2625 // for (i=0; i<0x100;i++) {
2626 // fprintf(stderr, " %d:%g ", i, copyReferenceHistogram[i]);
2627 // assert(copyReferenceHistogram[i] == referenceHistogram[i]);
2628 // total += copyReferenceHistogram[i];
2629 // }
2630 // fprintf(stderr, "\nTotal: %g\n", total);
2631 //
2632
2633 for (i = 0; i <= 0xff; i ++ ) {
2634
2635 double st_0 = copySourceHist[i];
2636
2637
2638 for (j = 0; j <= 0xff; j ++) {
2639
2640 if ( st_0 == 0.0 ) {
2641
2642 otherArray[j] = 0;
2643
2644 } else if (st_0 < copyReferenceHistogram[j] ) {
2645
2646 otherArray[j] = st_0;
2647 copyReferenceHistogram[j] -= st_0;
2648
2649 st_0 = 0.0;
2650
2651 } else {
2652 //
2653 otherArray[j] = copyReferenceHistogram[j];
2654
2655 st_0 -= copyReferenceHistogram[j];
2656
2657 copyReferenceHistogram[j] = 0.0;
2658
2659 }
2660
2661 }// for (j = 0; j <= 0xff; j ++) {
2662
2663 sum = 0.0;
2664 for (j = 0; j <= 0xff;j++) {
2665
2666 sum += otherArray[j];
2667
2668 }// if ( %j <= $0xff )
2669
2670 if ( sum == 0.0 ) {
2671 // Any value with sum 0 should be interpolated
2672 // between 2 other values
2673 // This is done at the end. In the meantime, mark it
2674 // with a -1
2675 if ( i == 0 ) {
2676
2677 curve[0] = 0;
2678
2679 continue;
2680
2681 } else {
2682
2683 if ( i == 0xff ) {
2684
2685 curve[0xff] = 255.0;
2686
2687 continue;
2688
2689 }
2690 curve[i] = -1.0;
2691
2692 }
2693 assert( curve[i] == -1.0);
2694 } else { // sum != 0.0
2695
2696 assert(sum != 0.0);
2697 // fprintf(stderr, "Value of sum %g %g \n", sum, otherArray[0]);
2698 sum2 = 0.0;
2699 for (j = 0; j <= 0xff; j ++ ) {
2700
2701 sum2 += otherArray[j] * j;
2702
2703 }
2704 // fprintf(stderr, "Value of sum2 %g\n", sum2);
2705
2706 curve[i] = sum2/sum;
2707
2708 } //else { //if ( ?? ) {
2709
2710
2711 } // for (i = 0; i <= 0xff; i ++ ) {
2712
2713
2714 // LET Us verify the results so far
2715 for (i=0;i<0x100;i++) {
2716 assert(curve[i] == -1 || curve[i] >= 0);
2717 assert(curve[i] < 0x100);
2718 }
2719
2720 // fprintf(stderr, "41->Magnolia: ");
2721 // for (i=0; i<0x100;i++) {
2722 // fprintf(stderr, " %d:%g ", i, curve[i]);
2723 // }
2724 // fprintf(stderr, "\n");
2725 //
2726 // // Check if the curve got any negative values.
2727
2728 for (i = 1; i <= 0xfe; i ++) {
2729
2730 if ( curve[i] == -1.0 ) {
2731
2732 // if the computed value was negative, then
2733 // interpolate between its 2 non-negative neighbors
2734
2735 int j;
2736
2737 j = i +1;
2738
2739 if ( j <= 0xff ) {
2740
2741 while ( curve[j] == -1.0 ) {
2742 ++j;
2743 if (j > 0xff) {
2744 break; //
2745 }
2746 }
2747 } // if ( %j <= $0xff )
2748
2749 assert(curve[j] >= 0);
2750 assert(curve[i-1] >= 0);
2751
2752 curve[i] = curve[i - 1] + (curve[j] - curve[i-1])/(j + 1 - i) ;
2753
2754 }
2755
2756 } // for (i = 1; i <= $0xfe; i ++) {
2757
2758 // LET Us verify the result
2759 for (i=0;i<0x100;i++) {
2760 assert(curve[i] >= 0);
2761 assert(curve[i] < 0x100);
2762 }
2763
2764
2765 // fprintf(stderr, "41->Magnolia: ");
2766 // for (i=0; i<0x100;i++) {
2767 // fprintf(stderr, " %d:%g ", i, curve[i]);
2768 // }
2769 // fprintf(stderr, "\n");
2770 //
2771 // exit(0);
2772
2773
2774 }
2775
2776 #ifdef OLDVERSION
2777
Unknown41(double * sourceHistogram,double * targetHistogram,double * magnoliaArray)2778 void Unknown41(double *sourceHistogram, double *targetHistogram, double *magnoliaArray)
2779 {
2780
2781
2782 double copyTargetHist[0x100];
2783 double copySourceHist[0x100];
2784 double otherArray[0x100];
2785 double contribution;
2786
2787 double sum;
2788
2789 int i,j;
2790
2791 int edx;
2792
2793
2794 memcpy(copySourceHist, sourceHistogram, 0x100 * sizeof(double)); // Check this one ???
2795
2796 memcpy(copyTargetHist, targetHistogram, 0x100 * sizeof(double)); // Check this one ???
2797
2798 for (i = 0; i <= 0xff; i ++ ) {
2799
2800
2801 contribution = copySourceHist[i];
2802
2803 for (j = 0; j <= 0xff; j ++) {
2804
2805 if ( copySourceHist[i] == 0.0 ) {
2806
2807 otherArray[j] = 0;
2808
2809 } else { // if ( ?? ) {
2810
2811 if (copySourceHist[i] < copyTargetHist[j] ) {
2812
2813 copyTargetHist[j] -= copySourceHist[i];
2814
2815 otherArray[j] = copySourceHist[i];
2816
2817 contribution = 0.0;
2818
2819 } else {
2820
2821 otherArray[j] = copyTargetHist[j];
2822
2823
2824 contribution = copySourceHist[i] - copyTargetHist[j];
2825
2826 copyTargetHist[j] = 0.0;
2827
2828 }
2829 } // end of else if ( ?? ) {
2830
2831
2832 }// for (j = 0; j <= 0xff; j ++) {
2833
2834
2835 sum = 0.0;
2836
2837
2838 for (edx = 0; edx <= 0xff;edx++) {
2839
2840 sum += otherArray[edx];
2841
2842 }// if ( %edx <= $0xff )
2843
2844 if ( sum == 0.0 ) {
2845
2846 contribution = 0.0;
2847
2848 if ( i == 0 ) {
2849
2850 magnoliaArray[0] = 0;
2851
2852 continue;
2853
2854 } else {
2855
2856 if (i == 0xff ) {
2857 magnoliaArray[0x7f8/8] = 255.0;
2858 continue;
2859
2860 }
2861
2862 contribution = -1;
2863
2864 }
2865
2866 assert(contribution == -1);
2867
2868 } else { //if ( ?? ) {
2869
2870
2871 sum = 0.0;
2872 for (edx = 0; edx <= 0xff; edx ++ ) {
2873 sum+= edx * otherArray[edx];
2874 }
2875 contribution /= sum;
2876
2877 }
2878
2879 magnoliaArray[i] = contribution;
2880
2881 } // for (i = 0; i <= 0xff; i ++ ) {
2882
2883
2884 for (edx = 1; edx <= 0xfe; edx ++) {
2885 int ecx;
2886
2887 if ( magnoliaArray[edx] == -1.0 ) {
2888
2889 ecx = edx +1;
2890
2891 if ( ecx <= 0xff ) {
2892
2893 while ( magnoliaArray[ecx] == -1.0 ) {
2894 if (++ecx > 0xff) {
2895 break;
2896 }
2897 }
2898 } //
2899
2900 magnoliaArray[edx] = magnoliaArray[edx-1] + (ecx - edx -1)/(magnoliaArray[ecx] - magnoliaArray[edx-1]);
2901
2902
2903 }//
2904
2905 } // for (edx = 1; edx <= $0xfe; edx ++) {
2906
2907 }
2908 #endif
2909
2910
2911
2912