I fixed the bug myself. Therefor I changed the files image.cpp and 
multiblend.cpp. I have these files in the attachment. But the file 
multiblend.cpp contains additionally my prior changes that enable the 
syntax @file for reading the image filenames from a file.

Furthermore, I fixed a memory leak bug in the file multiblend.cpp (the 
objects of class Image that are created via new Image(...) were not freed 
with delete at the end of main()).
Florian Königstein schrieb am Freitag, 17. Dezember 2021 um 04:48:07 UTC+1:

> Thank you for the update. I took it, but kept my changes in the file 
> multiblend.cpp so that I can use the syntax @file for reading the image 
> file names from a file.
>
> I took a large panorama project with 1351 images and 273471 control points 
> and preprocessed the images with nona.exe . Nona procuced at least one 
> image that has no content, but only transparency (the size is 250 x 152 
> pixels, GIMP says 10333 x 2583).
>
> Multiblend crashes when reading this image in the function   void 
> Image::Read()   in the file image.cpp . When running it with the debugger, 
> I see the following exception:
> Run-Time Check Failure #3 - The variable 'right' is being used without 
> being initialized (line 331 in file image.cpp  :   if (right < left) 
> std::swap(left, right); ).
> Probably in the preceding for-loops the if-interrogation never evaluate to 
> true so that 'right' is uninitialized.
> I attach the file with that the error occurs.
>
>
> Monkey schrieb am Sonntag, 12. Dezember 2021 um 15:57:15 UTC+1:
>
>> I finally got around to packaging up the files with the last few 
>> bugfixes. It doesn't include list file processing though, which is still on 
>> the to-do list.
>>
>> https://horman.net/multiblend/
>>
>> On Wednesday, 20 October 2021 at 15:51:01 UTC+1 [email protected] wrote:
>>
>>> @Florian K¨nigstein, @Monkey, do you plane to release a RC5 with the 
>>> fixes you identified ? (limited number of images because of windows args in 
>>> console, and silent crashes of Multiblend?)
>>> I almost have the same problems with big panoramas while charging tiff 
>>> files or blending: Multiblend stops with no errors.
>>>
>>> I also found that multiblend reads the geotiff tag, but in order to make 
>>> it work, I have to modify the georeference informations with Y=-Y.
>>> I'm trying to give you the exact lines of code to change.
>>>
>>> This tool is a game changer in th image assembly solutions!
>>> Best regards,
>>>
>>> Le dimanche 4 juillet 2021 à 20:36:19 UTC+2, Florian Königstein a écrit :
>>>
>>>> Very fine, it works now. Thank you for fixing.
>>>> For this panorama Multiblend took about 40 minutes. The final tiff file 
>>>> has about 14 GBytes. Windows could not open it (probably too large), but 
>>>> with GIMP I could view it.
>>>>
>>>> Monkey schrieb am Sonntag, 4. Juli 2021 um 12:13:05 UTC+2:
>>>>
>>>>> The most probable culprit is line 167 in pyramid.h which should be 
>>>>> changed to:
>>>>>
>>>>>     dst_p += *(size_t)*(sy - (level ? 1 : 0)) * pitch;
>>>>>
>>>>> Line 169 should also be changed to:
>>>>>
>>>>>     p_pt += *(size_t)*m128_pitch * sy;
>>>>>
>>>>> although it is probably not the cause of the problem in this case.
>>>>> On Sunday, 4 July 2021 at 05:48:50 UTC+1 Florian Königstein wrote:
>>>>>
>>>>>>
>>>>>> With the Debugger I see that it now crashed in the function  
>>>>>>  template <typename F> void Pyramid::OutPlanar8 _OP_  in the line
>>>>>>           _mm_storeu_si128(dst_pp_m++, pixels);
>>>>>> dst_pp_m points to an unreadable position. However *dst_p is readable 
>>>>>> (above dst_pp_m = (__m128i*)dst_p; was executed).
>>>>>>
>>>>>> Monkey schrieb am Samstag, 3. Juli 2021 um 21:51:58 UTC+2:
>>>>>>
>>>>>>> Actually never mind, that probably isn't it... I'll keep looking.
>>>>>>>
>>>>>>> On Saturday, 3 July 2021 at 20:46:48 UTC+1 Monkey wrote:
>>>>>>>
>>>>>>>> Line 683 of image.cpp - please replace
>>>>>>>>
>>>>>>>>     size_t channel_bytes = (width * height) << (bpp >> 4);
>>>>>>>>
>>>>>>>> with
>>>>>>>>
>>>>>>>>     size_t channel_bytes = ((size_t)width * height) << (bpp >> 4);
>>>>>>>>
>>>>>>>> I would test this myself but my one attempt caused my computer to 
>>>>>>>> freeze and need a hard reset, which I hate doing!
>>>>>>>> On Saturday, 3 July 2021 at 20:06:38 UTC+1 Florian Königstein wrote:
>>>>>>>>
>>>>>>>>> I changed it and ran it in a Release build because running the 
>>>>>>>>> Debug version takes much more time. The release version stopped with 
>>>>>>>>> the 
>>>>>>>>> following console output:
>>>>>>>>> Warning: test1191347.tif is fully obscured by other images
>>>>>>>>> Warning: test1191349.tif is fully obscured by other images
>>>>>>>>> Warning: test1191350.tif is fully obscured by other images
>>>>>>>>> Shrinking masks...
>>>>>>>>> Blending...
>>>>>>>>>
>>>>>>>>> So, no meesage "writing .. output file". And test119.tif (the 
>>>>>>>>> output file) contains 0 bytes. Probably the Debugger would raise an 
>>>>>>>>> error 
>>>>>>>>> at the same line. At night I will run the new version with the 
>>>>>>>>> debugger.
>>>>>>>>> Monkey schrieb am Samstag, 3. Juli 2021 um 20:16:46 UTC+2:
>>>>>>>>>
>>>>>>>>>> Ah, how foolish of me, I think I failed to account for the 
>>>>>>>>>> required number of bytes overflowing a 32-bit value.
>>>>>>>>>>
>>>>>>>>>> Does it work if you change line 31 of pyramid.cpp from
>>>>>>>>>>
>>>>>>>>>>         size_t bytes = pitch * ((height + y_shift + 3) & ~3) * 
>>>>>>>>>> sizeof(float);
>>>>>>>>>>
>>>>>>>>>> to
>>>>>>>>>>
>>>>>>>>>>         size_t bytes = (size_t)pitch * ((height + y_shift + 3) & 
>>>>>>>>>> ~3) * sizeof(float);
>>>>>>>>>>
>>>>>>>>>> ?
>>>>>>>>>> On Saturday, 3 July 2021 at 15:26:41 UTC+1 Florian Königstein 
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> I tested it so that more memory than I have RAM is used.
>>>>>>>>>>> But Multiblend crashed when stitching 1351 photos to an output 
>>>>>>>>>>> size of 88291 x 61567 Pixels.
>>>>>>>>>>>
>>>>>>>>>>> I started it again with the Visual Studio 2019 Debugger.
>>>>>>>>>>>
>>>>>>>>>>> The input options were Multiblend --argfile=test.txt with 
>>>>>>>>>>> test.txt containing
>>>>>>>>>>> --bigtiff
>>>>>>>>>>> --all-threads
>>>>>>>>>>> -f88291x61567+18125+0
>>>>>>>>>>> --compression=LZW
>>>>>>>>>>> -o
>>>>>>>>>>> test119.tif
>>>>>>>>>>> test1190000.tif
>>>>>>>>>>> test1190001.tif
>>>>>>>>>>> test1190002.tif
>>>>>>>>>>> ...
>>>>>>>>>>> test1191350.tif
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> The console output was
>>>>>>>>>>> ignoring Enblend option -f
>>>>>>>>>>>
>>>>>>>>>>> Multiblend v2.0.0 (c) 2021 David Horman        
>>>>>>>>>>> http://horman.net/multiblend/
>>>>>>>>>>>
>>>>>>>>>>> ----------------------------------------------------------------------------
>>>>>>>>>>> Processing test1190000.tif...
>>>>>>>>>>> Processing test1190001.tif...
>>>>>>>>>>> Processing test1190002.tif...
>>>>>>>>>>> ...
>>>>>>>>>>> Processing test1191350.tif...
>>>>>>>>>>>
>>>>>>>>>>> 88291 x 61567, 11 levels, 8 bpp
>>>>>>>>>>>
>>>>>>>>>>> Seaming...
>>>>>>>>>>> Warning: test1190001.tif is fully obscured by other images
>>>>>>>>>>> Warning: test1190002.tif is fully obscured by other images
>>>>>>>>>>>    .... some other images obscured ....
>>>>>>>>>>> Warning: test1191350.tif is fully obscured by other images
>>>>>>>>>>> Shrinking masks...
>>>>>>>>>>> Blending...
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> The crash was in the function CompositeLine() in the following 
>>>>>>>>>>> line:
>>>>>>>>>>>        if (i == 0) memset(&output_p[x], 0, mask_count << 2);
>>>>>>>>>>> The pointer output_p is 0x0000026321b85980, "Unable to read 
>>>>>>>>>>> memory". mask_count is 88291.
>>>>>>>>>>>
>>>>>>>>>>> In the calling function main() and the lambda function the 
>>>>>>>>>>> variables were:
>>>>>>>>>>>    i = 0
>>>>>>>>>>>    l = 0
>>>>>>>>>>> - in_level {width=12324 height=5502 pitch=12328 ...} 
>>>>>>>>>>> Pyramid::Level
>>>>>>>>>>> width 12324 int
>>>>>>>>>>> height 5502 int
>>>>>>>>>>> pitch 12328 int
>>>>>>>>>>> m128_pitch 3082 int
>>>>>>>>>>> bytes 271413248 unsigned __int64
>>>>>>>>>>> + data 0x0000025ed3e23080 {0.00000000} float *
>>>>>>>>>>> x 75967 int
>>>>>>>>>>> y 0 int
>>>>>>>>>>> x_shift true bool
>>>>>>>>>>> y_shift false bool
>>>>>>>>>>> upper_x_shift false bool
>>>>>>>>>>> upper_m128_pitch 0 int
>>>>>>>>>>> - bands { size=3 } std::vector<int,std::allocator<int>>
>>>>>>>>>>> [capacity] 3 __int64
>>>>>>>>>>> + [allocator] allocator 
>>>>>>>>>>> std::_Compressed_pair<std::allocator<int>,std::_Vector_val<std::_Simple_types<int>>,1>
>>>>>>>>>>> [0] 0 int
>>>>>>>>>>> [1] 2748 int
>>>>>>>>>>> [2] 5502 int
>>>>>>>>>>> + [Raw View] {_Mypair=allocator } 
>>>>>>>>>>> std::vector<int,std::allocator<int>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> - out_level {width=88291 height=61567 pitch=88296 ...} 
>>>>>>>>>>> Pyramid::Level
>>>>>>>>>>> width 88291 int
>>>>>>>>>>> height 61567 int
>>>>>>>>>>> pitch 88296 int
>>>>>>>>>>> m128_pitch 22074 int
>>>>>>>>>>> bytes 4564963328 unsigned __int64
>>>>>>>>>>> + data 0x0000026099c20000 {0.00000000} float *
>>>>>>>>>>> x 0 int
>>>>>>>>>>> y 0 int
>>>>>>>>>>> x_shift false bool
>>>>>>>>>>> y_shift false bool
>>>>>>>>>>> upper_x_shift false bool
>>>>>>>>>>> upper_m128_pitch 0 int
>>>>>>>>>>> - bands { size=3 } std::vector<int,std::allocator<int>>
>>>>>>>>>>> [capacity] 3 __int64
>>>>>>>>>>> + [allocator] allocator 
>>>>>>>>>>> std::_Compressed_pair<std::allocator<int>,std::_Vector_val<std::_Simple_types<int>>,1>
>>>>>>>>>>> [0] 0 int
>>>>>>>>>>> [1] 30780 int
>>>>>>>>>>> [2] 61567 int
>>>>>>>>>>> + [Raw View] {_Mypair=allocator } 
>>>>>>>>>>> std::vector<int,std::allocator<int>>
>>>>>>>>>>>
>>>>>>>>>>> x_offset = 75967
>>>>>>>>>>> y_offset = 0
>>>>>>>>>>> sy = 30780
>>>>>>>>>>> ey = 61567
>>>>>>>>>>> y = 30780
>>>>>>>>>>> The input files have a total size of 90 GBytes, but if nothing 
>>>>>>>>>>> helps, I could upload them somewhere.
>>>>>>>>>>>
>>>>>>>>>>> On Thursday I stitched the same project, but it was scaled with 
>>>>>>>>>>> nona so that the output became 36534 x 25476 Pixels. Multiblend 
>>>>>>>>>>> used only 
>>>>>>>>>>> the RAM (about 24 GBytes). It ran without problems.
>>>>>>>>>>>
>>>>>>>>>>> Maybe you can provocate this error when simulating less RAM e.g. 
>>>>>>>>>>> with a virtual machine and using less images.
>>>>>>>>>>> Monkey schrieb am Donnerstag, 1. Juli 2021 um 20:57:38 UTC+2:
>>>>>>>>>>>
>>>>>>>>>>>> Don't forget that the 64-bit version of Multiblend will use 
>>>>>>>>>>>> disk space (system temp or specified directory) if there isn't 
>>>>>>>>>>>> enough RAM.
>>>>>>>>>>>>
>>>>>>>>>>>> On Thursday, 1 July 2021 at 16:17:11 UTC+1 Florian Königstein 
>>>>>>>>>>>> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> For me it does not work if I use  multiblend.exe @test.txt  
>>>>>>>>>>>>> instead of  multiblend.exe --argfile=test.txt on Windows.
>>>>>>>>>>>>> AFAIK on Linux you can do this with some other syntax (I'm not 
>>>>>>>>>>>>> so familiar with Linux). But building something like an --argfile 
>>>>>>>>>>>>> option 
>>>>>>>>>>>>> into Multiblend has the other advantage that is works OS 
>>>>>>>>>>>>> independent.
>>>>>>>>>>>>>
>>>>>>>>>>>>> I have stitched nearly 1 GPixels with Multiblend in about 6 
>>>>>>>>>>>>> minutes. It's super fast. Thanks @ Monkey !
>>>>>>>>>>>>> The maximum memory usage was about 24 GBytes. My photos would 
>>>>>>>>>>>>> allow to stitch it to about 10 GPixels, but due to my RAM (64 
>>>>>>>>>>>>> GBytes) I 
>>>>>>>>>>>>> will probably only be able to stitch about 2 - 2.5 GPixels.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Monkey schrieb am Donnerstag, 1. Juli 2021 um 12:04:02 UTC+2:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> Thanks Florian, that's a great suggestion and I'll 
>>>>>>>>>>>>>> incorporate into the source distribution at some point. Out of 
>>>>>>>>>>>>>> interest, 
>>>>>>>>>>>>>> how long did the blend take? Was the final pixel count?
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Thursday, 1 July 2021 at 10:16:06 UTC+1 
>>>>>>>>>>>>>> [email protected] wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> AFAIK if you pass the parameter @filename to a program on ms 
>>>>>>>>>>>>>>> windows the contents of the file "filename" is used as 
>>>>>>>>>>>>>>> command-line 
>>>>>>>>>>>>>>> parameters. Thw last time I tried if the parameters are read 
>>>>>>>>>>>>>>> from a file 
>>>>>>>>>>>>>>> the maximum length was higher than the 256 bytes the limit was 
>>>>>>>>>>>>>>> at back then.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Am 1. Juli 2021 07:04:12 MESZ schrieb "Florian Königstein" <
>>>>>>>>>>>>>>> [email protected]>:
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> I tried to stitch a panorama with 1350 images with 
>>>>>>>>>>>>>>>> multiblend. It didn't work in Windows because the command line 
>>>>>>>>>>>>>>>> where all 
>>>>>>>>>>>>>>>> the image filenames are listed was longer than 32768 
>>>>>>>>>>>>>>>> characters. At least 
>>>>>>>>>>>>>>>> in Windows the limit is 32768 (or maybe one less).
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> I suggest adding the possibility to read the command line 
>>>>>>>>>>>>>>>> arguments from a file.
>>>>>>>>>>>>>>>> I changed multiblend.cpp so that you can add a command line 
>>>>>>>>>>>>>>>> option  --argfile filename  or  --argfile=filename . After 
>>>>>>>>>>>>>>>> this no further 
>>>>>>>>>>>>>>>> arguments may follow in the command line, but each line in the 
>>>>>>>>>>>>>>>> file 
>>>>>>>>>>>>>>>> "filename" counts as another argument, e.g. call
>>>>>>>>>>>>>>>> multiblend.exe" --argfile=test.txt
>>>>>>>>>>>>>>>> with test.txt containing e.g.
>>>>>>>>>>>>>>>> --compression=LZW
>>>>>>>>>>>>>>>> -o
>>>>>>>>>>>>>>>> test110.tif
>>>>>>>>>>>>>>>> --
>>>>>>>>>>>>>>>> test1100000.tif
>>>>>>>>>>>>>>>> test1100001.tif
>>>>>>>>>>>>>>>> ...
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> In the attachment I have the modified version of 
>>>>>>>>>>>>>>>> multiblend.cpp.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Maybe Hugin and HuginExecutor could be changed so that the 
>>>>>>>>>>>>>>>> arguments are written in a file if they are many.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Florian
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> [email protected] schrieb am Sonntag, 13. Juni 2021 um 
>>>>>>>>>>>>>>>> 11:55:00 UTC+2:
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Hello,
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> as there is actual coding of new software going on, maybe 
>>>>>>>>>>>>>>>>> one can iron out a deficiency in the Hugin lens model. At 
>>>>>>>>>>>>>>>>> least lay the 
>>>>>>>>>>>>>>>>> groundwork for it.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> The Brown-Conrady model parameters are sound, but the 
>>>>>>>>>>>>>>>>> intersection with the abc-Hugin parameter set contains only 
>>>>>>>>>>>>>>>>> one (1) 
>>>>>>>>>>>>>>>>> non-trivial distortion parameter.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> I suggest to add further Brown-Conrady parameters to the 
>>>>>>>>>>>>>>>>> software code you are currently writing. Now.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Best regards
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Klaus
>>>>>>>>>>>>>>>>> On 11.06.21 18:20, Florian Königstein wrote:
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Monkey, I much appreciate your software.
>>>>>>>>>>>>>>>>> I like it because I like big panoramas ... and the speedup 
>>>>>>>>>>>>>>>>> is welcome.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> For big panoramas there's another issue: Geometrical 
>>>>>>>>>>>>>>>>> optimization is slow.
>>>>>>>>>>>>>>>>> I developed a fork for the libpano library that I called 
>>>>>>>>>>>>>>>>> fastPTOptimizer.
>>>>>>>>>>>>>>>>> For large panoramas the speedup factor for optimization 
>>>>>>>>>>>>>>>>> can be 100 or more.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> I integrated both your multiblend and my fastPTOptimizer 
>>>>>>>>>>>>>>>>> into a "development version" of Hugin.
>>>>>>>>>>>>>>>>> Multiblend is now the default enblend-like program (in the 
>>>>>>>>>>>>>>>>> GUI is still written "enblend"
>>>>>>>>>>>>>>>>> but you can see that multiblend is used by choosing 
>>>>>>>>>>>>>>>>> Preferences / Programs).
>>>>>>>>>>>>>>>>> Only the CMakeLists.txt files are not updated so that 
>>>>>>>>>>>>>>>>> Multiblend is automatically integrated
>>>>>>>>>>>>>>>>> because I'm not yet so familiar with creating files for 
>>>>>>>>>>>>>>>>> CMake.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> My version of Hugin is here:
>>>>>>>>>>>>>>>>> https://sourceforge.net/projects/huginplusplus/files/
>>>>>>>>>>>>>>>>> development/
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Monkey schrieb am Samstag, 10. April 2021 um 22:00:35 
>>>>>>>>>>>>>>>>> UTC+2:
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> Has anyone out there tried either the x64 or x86 versions 
>>>>>>>>>>>>>>>>>> of Multiblend 2.0 on Windows XP or Windows Vista? Someone's 
>>>>>>>>>>>>>>>>>> reporting 
>>>>>>>>>>>>>>>>>> vcredist problems and I'm not sure if it's because I built 
>>>>>>>>>>>>>>>>>> using the latest 
>>>>>>>>>>>>>>>>>> platform toolset.
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> On Sunday, 4 April 2021 at 17:11:16 UTC+1 
>>>>>>>>>>>>>>>>>> [email protected] wrote:
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> I'll give it a shot. last time I used it it for a aerial 
>>>>>>>>>>>>>>>>>>> 360 it removed cars and other ground objects.
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> On Friday, March 5, 2021 at 3:57:30 PM UTC-8 Monkey 
>>>>>>>>>>>>>>>>>>> wrote:
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> *(* for a Gigapixel mosaic, anyway; it's complicated, 
>>>>>>>>>>>>>>>>>>>> see below)*
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> http://horman.net/multiblend/
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> It seems Groups won't let me post the quasi-essay I had 
>>>>>>>>>>>>>>>>>>>> written, complete with images, so the link above will have 
>>>>>>>>>>>>>>>>>>>> to suffice.
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> Here's Multiblend 2.0, faster, better, more... blendy. 
>>>>>>>>>>>>>>>>>>>> I'm calling it a Release Candidate because there's only so 
>>>>>>>>>>>>>>>>>>>> much testing I 
>>>>>>>>>>>>>>>>>>>> can stand to do, and I've hit a dead-end with features, so 
>>>>>>>>>>>>>>>>>>>> I thought I'd 
>>>>>>>>>>>>>>>>>>>> put it out there for people to try. I expect some bugs to 
>>>>>>>>>>>>>>>>>>>> be found pretty 
>>>>>>>>>>>>>>>>>>>> quickly, which I'll hopefully fix pretty quickly.
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> It's released under GPLv3.
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> -- 
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> A list of frequently asked questions is available at: 
>>>>>>>>>>>>>>>>> http://wiki.panotools.org/Hugin_FAQ
>>>>>>>>>>>>>>>>> --- 
>>>>>>>>>>>>>>>>> You received this message because you are subscribed to 
>>>>>>>>>>>>>>>>> the Google Groups "hugin and other free panoramic software" 
>>>>>>>>>>>>>>>>> group.
>>>>>>>>>>>>>>>>> To unsubscribe from this group and stop receiving emails 
>>>>>>>>>>>>>>>>> from it, send an email to [email protected].
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> To view this discussion on the web visit 
>>>>>>>>>>>>>>>>> https://groups.google.com/d/msgid/hugin-ptx/2699006b-895d-42f4-bd19-6ed0d3f3863bn%40googlegroups.com
>>>>>>>>>>>>>>>>>  
>>>>>>>>>>>>>>>>> <https://groups.google.com/d/msgid/hugin-ptx/2699006b-895d-42f4-bd19-6ed0d3f3863bn%40googlegroups.com?utm_medium=email&utm_source=footer>
>>>>>>>>>>>>>>>>> .
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> -- 
>>>>>>>>>>>>>>> Sent from my Android device with K-9 Mail. Please excuse my 
>>>>>>>>>>>>>>> brevity.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>

-- 
A list of frequently asked questions is available at: 
http://wiki.panotools.org/Hugin_FAQ
--- 
You received this message because you are subscribed to the Google Groups 
"hugin and other free panoramic software" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to [email protected].
To view this discussion on the web visit 
https://groups.google.com/d/msgid/hugin-ptx/4c7ed130-06eb-4a3e-bf8d-efc6253f4e17n%40googlegroups.com.
int hist_red[256];
int hist_grn[256];
int hist_blu[256];

class Channel {
public:
        Channel(size_t _bytes) : bytes(_bytes) {
                data = MapAlloc::Alloc(bytes);
        };
        ~Channel() { MapAlloc::Free(data); };
        void* data;
        size_t bytes;
        FILE* file = NULL;
};

class Image {
public:
        Image(char* _filename);
        ~Image();
        char* filename;
        ImageType type;
        int width;
        int height;
        int xpos;
        int ypos;
        int xpos_add = 0;
        int ypos_add = 0;
        std::vector<Channel*> channels;
        Pyramid* pyramid = NULL;
        GeoTIFFInfo geotiff;
        int tiff_width;
        int tiff_height;
        int tiff_u_height;
        int rows_per_strip;
        int first_strip;
        int end_strip;
        uint16_t bpp;
        uint16_t spp;
        void Open();
        void Read(void* data, bool gamma);
//      size_t untrimmed_pixels;
        size_t untrimmed_bytes;
        Flex* tiff_mask;
        float tiff_xres, tiff_yres;
        uint64_t mask_state;
        int mask_count;
        int mask_limit;
        bool seam_present;
        std::vector<Flex*> masks;
        void MaskPng(int i);

private:
        TIFF* tiff;
        FILE* file;
        struct jpeg_decompress_struct cinfo;
        struct jpeg_error_mgr jerr;
        png_structp png_ptr;
};

Image::Image(char* _filename) : filename(_filename) {
}

Image::~Image() {
        for (auto it = channels.begin(); it < channels.end(); ++it) delete 
(*it);
        for (auto it = masks.begin(); it < masks.end(); ++it) delete (*it);
        channels.clear();
        masks.clear();

        delete pyramid;
}

/***********************************************************************
************************************************************************
** Open
************************************************************************
***********************************************************************/
void Image::Open() {
        float tiff_xpos, tiff_ypos;
        uint16_t compression;

        char* ext = strrchr(filename, '.');
        if (!ext) {
                die("Could not identify file extension: %s", filename);
        }
        ++ext;

        if (!_stricmp(ext, "tif") || !_stricmp(ext, "tiff")) {
                type = ImageType::MB_TIFF;
        } else if (!_stricmp(ext, "jpg") || !_stricmp(ext, "jpeg")) {
                type = ImageType::MB_JPEG;
        } else if (!_stricmp(ext, "png")) {
                type = ImageType::MB_PNG;
        } else {
                die("Unknown file extension: %s", filename);
        }

        switch (type) {
                case ImageType::MB_TIFF: {
                        tiff = TIFFOpen(filename, "r");
                        if (!tiff) die("Could not open %s", filename);

                        if (!TIFFGetField(tiff, TIFFTAG_XPOSITION, &tiff_xpos)) 
tiff_xpos = -1;
                        if (!TIFFGetField(tiff, TIFFTAG_YPOSITION, &tiff_ypos)) 
tiff_ypos = -1;
                        if (!TIFFGetField(tiff, TIFFTAG_XRESOLUTION, 
&tiff_xres)) tiff_xres = -1;
                        if (!TIFFGetField(tiff, TIFFTAG_YRESOLUTION, 
&tiff_yres)) tiff_yres = -1;
                        TIFFGetField(tiff, TIFFTAG_IMAGEWIDTH, &tiff_width);
                        TIFFGetField(tiff, TIFFTAG_IMAGELENGTH, &tiff_height);
                        TIFFGetField(tiff, TIFFTAG_BITSPERSAMPLE, &bpp);
                        TIFFGetField(tiff, TIFFTAG_SAMPLESPERPIXEL, &spp);
                        TIFFGetField(tiff, TIFFTAG_ROWSPERSTRIP, 
&rows_per_strip);
                        TIFFGetField(tiff, TIFFTAG_COMPRESSION, &compression);

                        if (bpp != 8 && bpp != 16) {
                                printf("Invalid bpp %d (%s)", bpp, filename);
                                printf("%d, %d\n", tiff_width, tiff_height);
                                TIFFGetField(tiff, TIFFTAG_BITSPERSAMPLE, &bpp);
                                if (bpp != 8 && bpp != 16) die("Invalid bpp %d 
(%s)", bpp, filename);
                        }
//                      if (spp != 4) die("Images must be RGBA (%s)", filename);

                        geotiff.set = false;

                        if (tiff_xpos == -1 && tiff_ypos == -1) {
                                // try to read geotiff tags
                                if (geotiff_read(tiff, &geotiff)) {
                                        xpos = (int)(geotiff.XGeoRef / 
geotiff.XCellRes);
                                        ypos = (int)(geotiff.YGeoRef / 
geotiff.YCellRes);
                                } else {
                                        xpos = 0;
                                        ypos = 0;
                                }
                        } else {
                                if (tiff_xpos != -1 && tiff_xres > 0) xpos = 
(int)(tiff_xpos * tiff_xres + 0.5);
                                if (tiff_ypos != -1 && tiff_yres > 0) ypos = 
(int)(tiff_ypos * tiff_yres + 0.5);
                        }

                        first_strip = 0;
                        end_strip = TIFFNumberOfStrips(tiff);

                        int s;
                        tmsize_t temp;

                        if (tiff_xpos <= 0 && tiff_ypos <= 0 && compression != 
1 && spp == 4) {
                                tmsize_t min_stripsize = 0xffffffff;
                                int min_stripcount = 0;
                                for (s = 0; s < (int)TIFFNumberOfStrips(tiff) - 
1; ++s) {
                                        temp = TIFFRawStripSize(tiff, s);
                                        if (temp < min_stripsize) { 
min_stripsize = temp; min_stripcount = 1; } else if (temp == min_stripsize) 
min_stripcount++;
                                }

                                if (min_stripcount > 2) {
                                        first_strip = -1;
                                        for (s = 0; s < 
(int)TIFFNumberOfStrips(tiff); ++s) {
                                                temp = TIFFRawStripSize(tiff, 
s);
                                                if (temp != min_stripsize) {
                                                        if (first_strip == -1) 
first_strip = s;
                                                        end_strip = s + 1;
                                                }
                                        }
                                        if (first_strip == -1) first_strip = 0;
                                }
                        }

                        if (first_strip || end_strip != 
TIFFNumberOfStrips(tiff)) { // double check that min strips are (probably) 
transparent
                                tdata_t buf = 
_TIFFmalloc(TIFFScanlineSize(tiff));
                                if (first_strip) TIFFReadScanline(tiff, buf, 
0); else TIFFReadScanline(tiff, buf, tiff_u_height - 1);
                                bool trans;
                                switch (bpp) {
                                        case 8: trans = !(((uint32_t*)buf)[0] 
&& 0xff000000); break;
                                        case 16: trans = !(((uint64_t*)buf)[0] 
&& 0xffff000000000000); break;
                                }
                                if (!trans) {
                                        first_strip = 0;
                                        end_strip = TIFFNumberOfStrips(tiff);
                                }
                                _TIFFfree(buf);
                        }

                        ypos += first_strip * rows_per_strip;
                        int rows_missing = TIFFNumberOfStrips(tiff) * 
rows_per_strip - tiff_height;

                        tiff_u_height = (end_strip - first_strip) * 
rows_per_strip;
                        if (end_strip == TIFFNumberOfStrips(tiff)) 
tiff_u_height -= rows_missing;
                } break;
                case ImageType::MB_JPEG: {
                        fopen_s(&file, filename, "rb");
                        if (!file) die("Could not open %s", filename);

                        cinfo.err = jpeg_std_error(&jerr);
                        jpeg_create_decompress(&cinfo);
                        jpeg_stdio_src(&cinfo, file);
                        jpeg_read_header(&cinfo, TRUE);
                        jpeg_start_decompress(&cinfo);

                        if (!cinfo.output_width || !cinfo.output_height) {
                                die("Unknown JPEG format (%s)", filename);
                        }

                        if (cinfo.out_color_components != 3) {
                                die("Unknown JPEG format (%s)", filename);
                        }

                        tiff_width = cinfo.output_width;
                        tiff_height = tiff_u_height = cinfo.output_height;

                        bpp = 8;
                        spp = 3;

                        xpos = ypos = 0;
                        tiff_xpos = tiff_ypos = 0;
                        tiff_xres = tiff_yres = 90;
                } break;
                case ImageType::MB_PNG: {
                        fopen_s(&file, filename, "rb");
                        if (!file) die("Could not open %s", filename);

                        uint8_t sig[8];
                        size_t r = fread(sig, 1, 8, file); // assignment 
suppresses g++ -Ofast warning
                        if (!png_check_sig(sig, 8)) die("Bad PNG signature 
(%s)", filename);

                        png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, 
NULL, NULL, NULL);
                        if (!png_ptr) die("Error: libpng problem");
                        png_infop info_ptr = png_create_info_struct(png_ptr);
                        if (!info_ptr) die("Error: libpng problem");

                        png_init_io(png_ptr, file);
                        png_set_sig_bytes(png_ptr, 8);
                        png_read_info(png_ptr, info_ptr);

                        int png_colour;
                        uint32_t png_width, png_height;
                        int _bpp;
                        png_get_IHDR(png_ptr, info_ptr, &png_width, 
&png_height, &_bpp, &png_colour, NULL, NULL, NULL);
                        bpp = _bpp;
                        tiff_width = png_width;
                        tiff_u_height = tiff_height = png_height;

                        switch (png_colour) {
                                case PNG_COLOR_TYPE_RGB: spp = 3; break;
                                case PNG_COLOR_TYPE_RGBA: spp = 4; break;
                                default: die("Bad PNG colour type (%s)", 
filename);
                        }

                        if (bpp != 8 && bpp != 16) {
                                die("Bad bit depth (%s)", filename);
                        }

                        xpos = ypos = 0;
                        tiff_xpos = tiff_ypos = 0;
                        tiff_xres = tiff_yres = 90;
                } break;
        }

        xpos += xpos_add;
        ypos += ypos_add;

        size_t untrimmed_pixels = (size_t)tiff_u_height * tiff_width;
        untrimmed_bytes = (untrimmed_pixels * spp ) << (bpp >> 4);
}

/***********************************************************************
************************************************************************
** Read
************************************************************************
***********************************************************************/
void Image::Read(void* data, bool gamma) {
        Output(1, "Processing %s...", filename);

        switch (type) {
                case ImageType::MB_TIFF: {
                        char* pointer = (char*)data;

                        for (int s = first_strip; s < end_strip; s++) {
                                auto strip_size = TIFFReadEncodedStrip(tiff, s, 
pointer, -1);
                                pointer += strip_size;
                        }
                } break;
                case ImageType::MB_JPEG: {
                        uint8_t* pointer = (uint8_t*)data;

                        while (cinfo.output_scanline < cinfo.output_height) {
                                jpeg_read_scanlines(&cinfo, &pointer, 1);
                                pointer += (tiff_width * spp) << (bpp >> 4);
                        }
                } break;
                case ImageType::MB_PNG: {
                        uint8_t* pointer = (uint8_t*)data;

                        for (int y = 0; y < tiff_height; ++y) {
                                png_read_row(png_ptr, pointer, NULL);
                                pointer += (tiff_width * spp) << (bpp >> 4);
                        }
                } break;
        }

/***********************************************************************
* Trim
***********************************************************************/
        int x, y;
        int top, left = -1, bottom, right = -1;

        if (spp == 4) {
                switch (bpp) {
                        case 8: {
                                uint32_t* p = (uint32_t*)data;
                                p--;

                                for (y = 0; y < tiff_u_height; ++y) {
                                        for (x = 0; x < tiff_width; ++x) {
                                                if (*++p >= 0xff000000) {
                                                        left = x;
                                                        top = y;
                                                        y = tiff_u_height;
                                                        break;
                                                }
                                        }
                                }

                                p = ((uint32_t*)data) + tiff_width * 
tiff_u_height;
                                for (y = tiff_u_height - 1; y >= 0; --y) {
                                        for (x = tiff_width - 1; x >= 0; --x) {
                                                if (*--p >= 0xff000000) {
                                                        right = x;
                                                        bottom = y;
                                                        y = -1;
                                                        break;
                                                }
                                        }
                                }

                                if(-1 == left || -1 == right)
                                {
                                        width = height = 0;   // Error: All 
pixels are transparent. This image must not be used
                                        return;
                                }
                                else
                                {
                                        if (right < left) std::swap(left, 
right);

                                        p = ((uint32_t*)data) + tiff_width * 
top;
                                        for (y = top; y <= bottom; ++y) {
                                                for (x = 0; x < left; ++x) {
                                                        if (p[x] >= 0xff000000) 
{
                                                                left = x;
                                                                break;
                                                        }
                                                }

                                                for (x = tiff_width - 1; x > 
right; --x) {
                                                        if (p[x] >= 0xff000000) 
{
                                                                right = x;
                                                                break;
                                                        }
                                                }

                                                p += tiff_width;
                                        }
                                }
                        } break;
                        case 16: {
                                uint16_t* p = ((uint16_t*)data) + 3;

                                for (y = 0; y < tiff_u_height; ++y) {
                                        for (x = 0; x < tiff_width; ++x) {
                                                if (p[x << 2] == 0xffff) {
                                                        left = x;
                                                        top = y;
                                                        y = tiff_u_height;
                                                        break;
                                                }
                                        }

                                        p += tiff_width << 2;
                                }

                                p = ((uint16_t*)data) + (tiff_width << 2) * 
(tiff_u_height - 1) + 3;
                                for (y = tiff_u_height - 1; y >= 0; --y) {
                                        for (x = tiff_width - 1; x >= 0; --x) {
                                                if (p[x << 2] == 0xffff) {
                                                        right = x;
                                                        bottom = y;
                                                        y = -1;
                                                        break;
                                                }
                                        }

                                        p -= tiff_width << 2;
                                }

                                if (-1 == left || -1 == right)
                                {
                                        width = height = 0;   // Error: All 
pixels are transparent. This image must not be used
                                        return;
                                }
                                else
                                {
                                        if (right < left) std::swap(left, 
right);

                                        p = ((uint16_t*)data) + (tiff_width << 
2) * top + 3;
                                        for (y = top; y <= bottom; ++y) {
                                                for (x = 0; x < left; ++x) {
                                                        if (p[x << 2] == 
0xffff) {
                                                                left = x;
                                                                break;
                                                        }
                                                }

                                                for (x = tiff_width - 1; x > 
right; --x) {
                                                        if (p[x << 2] == 
0xffff) {
                                                                right = x;
                                                                break;
                                                        }
                                                }

                                                p += tiff_width << 2;
                                        }
                                }
                        } break;
                }

                width = right + 1 - left;
                height = bottom + 1 - top;

                xpos += left;
                ypos += top;

/***********************************************************************
* Inpaint
***********************************************************************/
                int temp_copy;
                uint32_t a, b, c, d;
                uint32_t* this_line = NULL;
                uint32_t* prev_line = NULL;
                Threadpool* threadpool = Threadpool::GetInstance();

                tiff_mask = new Flex(width, height);
                Flex* dt = new Flex(width, height);
                int mc;

                int n_threads = (std::max)(2, threadpool->GetNThreads());

                uint32_t** thread_lines = new uint32_t * [n_threads];
                uint32_t** thread_comp_lines = new uint32_t * [n_threads];

                std::mutex* flex_mutex_p = new std::mutex;
                std::condition_variable* flex_cond_p = new 
std::condition_variable;

                for (int i = 0; i < n_threads; ++i) {
                        thread_lines[i] = new uint32_t[width];
                        thread_comp_lines[i] = new uint32_t[width];
                }

                uint32_t* bitmap32 = NULL;
                uint64_t* bitmap64 = NULL;
                if (bpp == 8) bitmap32 = ((uint32_t*)data) + top * tiff_width + 
left; else bitmap64 = ((uint64_t*)data) + top * tiff_width + left;

                for (y = 0; y < height; ++y) {
                        int t = y % n_threads;
                        this_line = thread_lines[t];
                        uint32_t* comp = thread_comp_lines[t];
                        bool first;

                        x = 0;

                        { // make sure the last compression thread to use this 
chunk of memory is finished
                                std::unique_lock<std::mutex> 
mlock(*flex_mutex_p);
                                flex_cond_p->wait(mlock, [=] { return dt->y > y 
- n_threads; });
                        }

                        while (x < width) {
                                mc = 0;
                                first = true;

                                while (x < width && (bpp == 8 ? (bitmap32[x] < 
0xff000000) : (bitmap64[x] < 0xffff000000000000))) {
                                        if (y == 0) {
                                                if (x == 0) {
                                                        this_line[0] = 
0x80000000;
                                                } else {
                                                        this_line[x] = 
this_line[x - 1] + 3;
                                                        if (bpp == 8) {
                                                                bitmap32[x] = 
bitmap32[x - 1];
                                                        } else {
                                                                bitmap64[x] = 
bitmap64[x - 1];
                                                        }
                                                }
                                        } else {
                                                if (x == 0) {
                                                        b = prev_line[x] + 3;
                                                        c = prev_line[x + 1] + 
4;

                                                        if (b < c) {
                                                                d = b;
                                                                temp_copy = x - 
tiff_width;
                                                        } else {
                                                                d = c;
                                                                temp_copy = x - 
tiff_width + 1;
                                                        }
                                                        d = b < c ? b : c;

                                                        this_line[x] = d;
                                                        if (bpp == 8) {
                                                                bitmap32[x] = 
bitmap32[temp_copy];
                                                        } else {
                                                                bitmap64[x] = 
bitmap64[temp_copy];
                                                        }

                                                        a = b + 1;
                                                        b = c - 1;
                                                        d += 3;
                                                } else {
                                                        if (first) {
                                                                a = prev_line[x 
- 1] + 4;
                                                                b = 
prev_line[x] + 3;
                                                                d = this_line[x 
- 1] + 3;
                                                                first = false;
                                                        }

                                                        if (x < width - 1) c = 
prev_line[x + 1] + 4; else c = 0xffffffff;
                                                        temp_copy = x - 1;

                                                        if (a < d) { d = a; 
temp_copy = x - tiff_width - 1; }
                                                        if (b < d) { d = b; 
temp_copy = x - tiff_width; }
                                                        if (c < d) { d = c; 
temp_copy = x - tiff_width + 1; }

                                                        this_line[x] = d;
                                                        if (bpp == 8) {
                                                                bitmap32[x] = 
bitmap32[temp_copy];
                                                        } else {
                                                                bitmap64[x] = 
bitmap64[temp_copy];
                                                        }

                                                        a = b + 1;
                                                        b = c - 1;
                                                        d += 3;
                                                }
                                        }
                                        ++x;
                                        ++mc;
                                }

                                if (mc) {
                                        tiff_mask->Write32(mc);
                                        mc = 0;
                                }

                                switch (bpp) {
                                        case 8: {
                                                while (x < width && 
(bitmap32[x] >= 0xff000000)) {
                                                        this_line[x++] = 0;
                                                        ++mc;
                                                }
                                        } break;
                                        case 16: {
                                                while (x < width && 
(bitmap64[x] >= 0xffff000000000000)) {
                                                        this_line[x++] = 0;
                                                        ++mc;
                                                }
//                                              }
                                        } break;
                                }

                                if (mc) {
                                        tiff_mask->Write32(0x80000000 | mc);
                                        mc = 0;
                                }
                        }

                        if (bpp == 8) bitmap32 += tiff_width; else bitmap64 += 
tiff_width;

                        if (y < height - 1) {
                                threadpool->Queue([=] {
                                        int p = CompressDTLine(this_line, 
(uint8_t*)comp, width);
                                        {
                                                std::unique_lock<std::mutex> 
mlock(*flex_mutex_p);
                                                flex_cond_p->wait(mlock, [=] { 
return dt->y == y; });
                                                dt->Copy((uint8_t*)comp, p);
                                                dt->NextLine();
                                        }
                                        flex_cond_p->notify_all();
                                });
                        }

                        tiff_mask->NextLine();
                        prev_line = this_line;
                }

                threadpool->Wait();

                // backward
                int current_count = 0;
                int current_step;
                uint32_t dt_val;

                uint32_t mask;

                prev_line = thread_lines[(y - 2) % n_threads];

                // first line
                x = width - 1;
                if (bpp == 8) bitmap32 -= tiff_width; else bitmap64 -= 
tiff_width;

                while (x >= 0) {
                        mask = tiff_mask->ReadBackwards32();
                        if (mask & 0x80000000) { // solid
                                x -= mask & 0x7fffffff;
                                d = 3;
                        } else {
                                if (x == width - 1) {
                                        d = this_line[x] + 3;
                                        --mask;
                                        --x;
                                }
                                while (mask) {
                                        uint32_t best = this_line[x];
                                        if (d < best) {
                                                this_line[x] = d;
                                                if (bpp == 8) bitmap32[x] = 
bitmap32[x + 1]; else bitmap64[x] = bitmap64[x + 1];
                                                d += 3;
                                        } else {
                                                d = best + 3;
                                        }
                                        --x;
                                        --mask;
                                }
                        }
                }

                std::swap(this_line, prev_line);

                // other lines
                for (y = height - 2; y >= 0; --y) {
                        x = width - 1;
                        c = d = 0x80000000;
                        if (bpp == 8) bitmap32 -= tiff_width; else bitmap64 -= 
tiff_width;

                        while (x >= 0) {
                                mask = tiff_mask->ReadBackwards32();

                                if (mask & 0x80000000) { // solid
                                        mc = mask & 0x7fffffff;
                                        x -= mc;
                                        memset(&this_line[x + 1], 0, mc << 2);
                                        d = 3;
                                } else {
                                        b = prev_line[x] + 3;
                                        c = (x < width - 1) ? prev_line[x + 1] 
+ 4 : 0x80000000;
                                        while (mask) {
                                                a = (x > 0) ? prev_line[x - 1] 
+ 4 : 0x80000000;
                                                INPAINT_DT;
                                                int copy = 0;
                                                uint32_t best = dt_val;
                                                if (a < best) { best = a; copy 
= tiff_width + x - 1; }
                                                if (b < best) { best = b; copy 
= tiff_width + x; }
                                                if (c < best) { best = c; copy 
= tiff_width + x + 1; }
                                                if (d < best) { best = d; copy 
= x + 1; }

                                                if (copy) {
                                                        if (bpp == 8) 
bitmap32[x] = bitmap32[copy]; else bitmap64[x] = bitmap64[copy];
                                                }

                                                this_line[x--] = best;
                                                c = b + 1;
                                                d = best + 3;
                                                b = a - 1;
                                                mask--;
                                        }
                                }
                        }

                        std::swap(this_line, prev_line);
                }

                delete flex_mutex_p;
                delete flex_cond_p;

                for (int i = 0; i < n_threads; ++i) {
                        delete[] thread_lines[i];
                        delete[] thread_comp_lines[i];
                }

                delete[] thread_lines;
                delete[] thread_comp_lines;
        } else {
                width = tiff_width;
                height = tiff_height;

                tiff_mask = new Flex(width, height);

                for (y = 0; y < height; ++y) {
                        tiff_mask->Write32(0x80000000 | width);
                        tiff_mask->NextLine();
                }
        }

/***********************************************************************
* Extract channels
***********************************************************************/
        size_t channel_bytes = ((size_t)width * height) << (bpp >> 4);

        for (int c = 0; c < 3; ++c) {
                channels.push_back(new Channel(channel_bytes));
        }

        if (spp == 4) {
                switch (bpp) {
                        case 8: {
                                uint32_t* line = ((uint32_t*)data) + top * 
tiff_width + left;
                                int p = 0;
                                for (y = 0; y < height; ++y) {
                                        for (x = 0; x < width; ++x) {
                                                uint32_t pixel = line[x];
                                                ((uint8_t*)channels[0]->data)[p 
+ x] = pixel & 0xff;
                                                ((uint8_t*)channels[1]->data)[p 
+ x] = (pixel >> 8) & 0xff;
                                                ((uint8_t*)channels[2]->data)[p 
+ x] = (pixel >> 16) & 0xff;
                                        }
                                        p += width;
                                        line += tiff_width;
                                }
                        } break;
                        case 16: {
                                uint64_t* line = ((uint64_t*)data) + top * 
tiff_width + left;
                                int p = 0;
                                for (y = 0; y < height; ++y) {
                                        for (x = 0; x < width; ++x) {
                                                uint64_t pixel = line[x];
                                                
((uint16_t*)channels[0]->data)[p + x] = pixel & 0xffff;
                                                
((uint16_t*)channels[1]->data)[p + x] = (pixel >> 16) & 0xffff;
                                                
((uint16_t*)channels[2]->data)[p + x] = (pixel >> 32) & 0xffff;
                                        }
                                        p += width;
                                        line += tiff_width;
                                }
                        } break;
                }
        } else {
                switch (bpp) {
                        case 8: {
                                uint8_t byte;
                                uint8_t* bytes = (uint8_t*)data;
                                size_t p = 0;
                                for (y = 0; y < height; ++y) {
                                        for (x = 0; x < width; ++x) {
                                                byte = *bytes++;
                                                
((uint8_t*)channels[0]->data)[p] = byte;
//                                              channel_totals[0] += gamma ? 
byte*byte : byte;

                                                byte = *bytes++;
                                                
((uint8_t*)channels[1]->data)[p] = byte;
//                                              channel_totals[1] += gamma ? 
byte * byte : byte;

                                                byte = *bytes++;
                                                
((uint8_t*)channels[2]->data)[p] = byte;
//                                              channel_totals[2] += gamma ? 
byte * byte : byte;

                                                ++p;
                                        }
                                }
                        } break;
                        case 16: {
                                uint16_t word;
                                uint16_t* words = (uint16_t*)data;
                                size_t p = 0;
                                for (y = 0; y < height; ++y) {
                                        for (x = 0; x < width; ++x) {
                                                word = *words++;
                                                
((uint16_t*)channels[0]->data)[p] = word;
//                                              channel_totals[0] += gamma ? 
word * word : word;

                                                word = *words++;
                                                
((uint16_t*)channels[1]->data)[p] = word;
//                                              channel_totals[1] += gamma ? 
word * word : word;

                                                word = *words++;
                                                
((uint16_t*)channels[2]->data)[p] = word;
//                                              channel_totals[2] += gamma ? 
word * word : word;

                                                ++p;
                                        }
                                }
                        } break;
                }

//              total_pixels += width * height;
        }

        Output(1, "\n");
}

/***********************************************************************
* Debugging
***********************************************************************/
void Image::MaskPng(int i) {
        char filename[256];
        sprintf_s(filename, "masks-%d.png", i);

        int width = masks[0]->width;
        int height = masks[0]->height; // +1 + masks[1]->height;

        size_t size = (size_t)width * height;
        uint8_t* temp = (uint8_t*)malloc(size);
        memset(temp, 32, size);

        int px = 0, py = 0;
        uint32_t cur;

        for (int l = 0; l < (int)masks.size(); ++l) {
                uint32_t* data = (uint32_t*)masks[l]->data;
                uint8_t* line = temp + py * masks[0]->width + px;
                for (int y = 0; y < masks[l]->height; ++y) {
                        int x = 0;
                        while (x < masks[l]->width) {
                                float val;
                                int count;

                                cur = *data++;

                                if (cur & 0x80000000) {
                                        count = cur & 0x00ffffff;
                                        if (cur & 0x20000000) val = 
*(float*)data++; else val = (float)((cur >> 30) & 1);
                                } else {
                                        val = *((float*)&cur);
                                        count = 1;
                                }

                                int t = x + count;
                                while (x < t) line[x++] = (uint8_t)(val * 255);
                        }

                        line += masks[0]->width;
                }
                if (l & 1) px += masks[l]->width + 1; else py += 
masks[l]->height + 1;
                break;
        }

        Pnger::Quick(filename, temp, width, height, width, PNG_COLOR_TYPE_GRAY);
}
/*
        Multiblend 2.0 (c) 2021 David Horman

        This program is free software: you can redistribute it and/or modify
        it under the terms of the GNU General Public License as published by
        the Free Software Foundation, either version 3 of the License, or
        (at your option) any later version.

        This program is distributed in the hope that it will be useful,
        but WITHOUT ANY WARRANTY; without even the implied warranty of
        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
        GNU General Public License for more details.

        You should have received a copy of the GNU General Public License
        along with this program. If not, see <https://www.gnu.org/licenses/>.

        The author can be contacted at [email protected]
*/

#define NOMINMAX
#include <stdio.h>
#include <stdint.h>
#include <vector>
#include <algorithm>
#ifdef __APPLE__
#define memalign(a,b) malloc((b))
#else
#include <malloc.h>
#endif

#include "tiffio.h"
#include "jpeglib.h"

#ifndef _WIN32
#include <strings.h>
int _stricmp(const char* a, const char* b) { return strcasecmp(a, b); }
#define ZeroMemory(a,b) memset(a,0,b)
#define sprintf_s sprintf
#define sscanf_s sscanf
void* _aligned_malloc(size_t size, int boundary) { return memalign(boundary, 
size); }
void _aligned_free(void* a) { free(a); }
void fopen_s(FILE** f, const char* filename, const char* mode) { *f = 
fopen(filename, mode); }
#endif

int verbosity = 1;

#include "pnger.cpp"
#include "pyramid.cpp"
#include "functions.cpp"

#include "mapalloc.cpp"
#include "threadpool.cpp"
#include "geotiff.cpp"

class PyramidWithMasks : public Pyramid {
public:
        using Pyramid::Pyramid;
        std::vector<Flex*> masks;
};

enum class ImageType { MB_NONE, MB_TIFF, MB_JPEG, MB_PNG };

#include "image.cpp"

#ifdef _WIN32
FILE _iob[] = { *stdin, *stdout, *stderr };

extern "C" FILE * __cdecl __iob_func(void) {
        return _iob;
}
#pragma comment(lib, "legacy_stdio_definitions.lib")
// the above are required to support VS 2010 build of libjpeg-turbo 2.0.6
#pragma comment(lib, "tiff.lib")
#pragma comment(lib, "turbojpeg.lib")
#pragma comment(lib, "libpng16.lib")
#pragma comment(lib, "zlib.lib")
#pragma comment(lib, "lzma.lib")
#endif

#define MASKVAL(X) (((X) & 0x7fffffffffffffff) | images[(X) & 
0xffffffff]->mask_state)

int main(int argc, char* argv[]) {
// This is here because of a weird problem encountered during development with 
Visual Studio. It should never be triggered.
        if (verbosity != 1) {
                printf("bad compile?\n");
                exit(EXIT_FAILURE);
        }

        int i, j;
        Timer timer_all, timer;
        timer_all.Start();

        TIFFSetWarningHandler(NULL);

/***********************************************************************
* Variables
***********************************************************************/
        std::vector<Image*> images;
        int fixed_levels = 0;
        int add_levels = 0;

        int width = 0;
        int height = 0;

        bool no_mask = false;
        bool big_tiff = false;
        bool bgr = false;
        bool wideblend = false;
        bool reverse = false;
        bool timing = false;
        bool dither = true;
        bool gamma = false;
        bool all_threads = true;
        int wrap = 0;

        TIFF* tiff_file = NULL;
        FILE* jpeg_file = NULL;
        Pnger* png_file = NULL;
        ImageType output_type = ImageType::MB_NONE;
        int jpeg_quality = -1;
        int compression = -1;
        char* seamsave_filename = NULL;
        char* seamload_filename = NULL;
        char* xor_filename = NULL;
        char* output_filename = NULL;
        int output_bpp = 0;

        double images_time = 0;
        double copy_time = 0;
        double seam_time = 0;
        double shrink_mask_time = 0;
        double shrink_time = 0;
        double laplace_time = 0;
        double blend_time = 0;
        double collapse_time = 0;
        double wrap_time = 0;
        double out_time = 0;
        double write_time = 0;

/***********************************************************************
* Help
***********************************************************************/
        if (argc == 1 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") 
|| !strcmp(argv[1], "/?")) {
                Output(1, "\n");
                Output(1, "Multiblend v2.0.0 (c) 2021 David Horman        
http://horman.net/multiblend/\n";);
                Output(1, 
"----------------------------------------------------------------------------\n");

                printf("Usage: multiblend [options] [-o OUTPUT] INPUT [X,Y] 
[INPUT] [X,Y] [INPUT]...\n");
                printf("\n");
                printf("Options:\n");
                printf("  --levels X / -l X      X: set number of blending 
levels to X\n");
                printf("                        -X: decrease number of blending 
levels by X\n");
                printf("                        +X: increase number of blending 
levels by X\n");
                printf("  --depth D / -d D       Override automatic output 
image depth (8 or 16)\n");
                printf("  --bgr                  Swap RGB order\n");
                printf("  --wideblend            Calculate number of levels 
based on output image size,\n");
                printf("                         rather than input image 
size\n");
                printf("  -w, --wrap=[mode]      Blend around images boundaries 
(NONE (default),\n");
                printf("                         HORIZONTAL, VERTICAL). When 
specified without a mode,\n");
                printf("                         defaults to HORIZONTAL.\n");
                printf("  --compression=X        Output file compression. For 
TIFF output, X may be:\n");
                printf("                         NONE (default), PACKBITS, or 
LZW\n");
                printf("                         For JPEG output, X is JPEG 
quality (0-100, default 75)\n");
                printf("                         For PNG output, X is PNG 
filter (0-9, default 3)\n");
                printf("  --cache-threshold=     Allocate memory beyond X 
bytes/[K]ilobytes/\n");
                printf("      X[K/M/G]           [M]egabytes/[G]igabytes to 
disk\n");
                printf("  --no-dither            Disable dithering\n");
                printf("  --tempdir <dir>        Specify temporary directory 
(default: system temp)\n");
                printf("  --save-seams <file>    Save seams to PNG file for 
external editing\n");
                printf("  --load-seams <file>    Load seams from PNG file\n");
                printf("  --no-output            Do not blend (for use with 
--save-seams)\n");
                printf("                         Must be specified as last 
option before input images\n");
                printf("  --bigtiff              BigTIFF output\n");
                printf("  --reverse              Reverse image priority 
(last=highest) for resolving\n");
                printf("                         indeterminate pixels\n");
                printf("  --quiet                Suppress output (except 
warnings)\n");
                printf("  --all-threads          Use all available CPU 
threads\n");
                printf("  [X,Y]                  Optional position adjustment 
for previous input image\n");
                printf("  @file                  Read following command line 
arguments from lines in file.\n");
                exit(EXIT_SUCCESS);
        }

/***********************************************************************
************************************************************************
* Parse arguments
************************************************************************
***********************************************************************/
        std::vector<char*> my_argv;

        bool skip = false;
        long argbufferlength;
        char *argbuffermem = 0, *argbuffer, *curarg;
        FILE* argfile = 0;
        char *argfilename = 0;

        for (i = 1; i < argc; ++i) {

                if ('@' == argv[i][0]) {
                        if (0 == argv[i][1]) {
                                die("@ requires the filename after the @ 
character");
                        }
                        argfilename = &argv[i][1];
                        fopen_s(&argfile, argfilename, "r");
                        if (0 == argfile) {
                                die("cannot open argument file %s\n", 
argfilename);
                        }
                        if (++i < argc) {
                                die("@ requires the filename after the @ 
character, but no other arguments may follow");
                        }
                        break;
                }

                my_argv.push_back(argv[i]);

                if (!skip) {
                        int c = 0;

                        while (argv[i][c]) {
                                if (argv[i][c] == '=') {
                                        argv[i][c++] = 0;

                                        if (argv[i][c]) {
                                                my_argv.push_back(&argv[i][c]);
                                        }
                                        break;
                                }
                                ++c;
                        }

                        if (!strcmp(argv[i], "-o") || !strcmp(argv[i], 
"--output")) {
                                skip = true;
                        }
                }
        }

        if(argfile)
        {
                long prevpos, filesize;
                if(-1 == (prevpos = ftell(argfile)) || 0 != fseek(argfile, 0L, 
SEEK_END) ||
                        -1 == (filesize = ftell(argfile)) || 0 != 
fseek(argfile, prevpos, SEEK_SET)) {
                        fclose(argfile);
                        die("error reading argument file");
                }

                argbufferlength = filesize + 1;
                if(0 == (argbuffermem = (char*)malloc(argbufferlength * 
sizeof(char)))) {
                        fclose(argfile);
                        die("cannot allocate memory for argument file");
                }

                argbuffermem[argbufferlength - 1] = 0;
                argbuffer = argbuffermem;

                while (fgets(argbuffer, argbufferlength, argfile))
                {
                        curarg = argbuffer;
                        bool newline_found = false;
                        if(0 != argbuffer[0])
                        {
                                for(curarg = argbuffer + strnlen_s(argbuffer, 
argbufferlength) ; curarg-- != argbuffer ; )
                                {
                                        if ('\n' == *curarg)
                                                newline_found = true;
                                        else if ('\r' != *curarg)
                                                break;
                                        *curarg = 0;
                                }
                                for (curarg = argbuffer; '\r' == *curarg || 
'\n' == *curarg; curarg++)
                                        { }
                                if (!newline_found)
                                {
                                        printf("Warning: ignoring last line of 
file %s because it does not end with a newline character.\n", argfilename);
                                }
                        }

                        argbufferlength -= (long)(strlen(curarg) + 1);
                        argbuffer = curarg + (strlen(curarg) + 1);

                        if(newline_found)
                                my_argv.push_back(curarg);

                        if (!skip) {
                                int c = 0;

                                while (curarg[c]) {
                                        if (curarg[c] == '=') {
                                                curarg[c++] = 0;
                                                if (curarg[c]) {
                                                        
my_argv.push_back(&curarg[c]);
                                                }
                                                break;
                                        }
                                        ++c;
                                }

                                if (!strcmp(curarg, "-o") || !strcmp(curarg, 
"--output")) {
                                        skip = true;
                                }
                        }
                }

                fclose(argfile);
        }

        if ((int)my_argv.size() < 3) die("Error: Not enough arguments (try -h 
for help)");

        for (i = 0; i < (int)my_argv.size(); ++i) {
                if (!strcmp(my_argv[i], "-d") || !strcmp(my_argv[i], "--d") || 
!strcmp(my_argv[i], "--depth") || !strcmp(my_argv[i], "--bpp")) {
                        if (++i < (int)my_argv.size()) {
                                output_bpp = atoi(my_argv[i]);
                                if (output_bpp != 8 && output_bpp != 16) {
                                        die("Error: Invalid output depth 
specified");
                                }
                        } else {
                                die("Error: Missing parameter value");
                        }
                }               else if (!strcmp(my_argv[i], "-l") || 
!strcmp(my_argv[i], "--levels")) {
                        if (++i < (int)my_argv.size()) {
                                int n;
                                if (my_argv[i][0] == '+' || my_argv[i][0] == 
'-') {
                                        sscanf_s(my_argv[i], "%d%n", 
&add_levels, &n);
                                } else {
                                        sscanf_s(my_argv[i], "%d%n", 
&fixed_levels, &n);
                                        if (fixed_levels == 0) fixed_levels = 1;
                                }
                                if (my_argv[i][n]) die("Error: Bad --levels 
parameter");
                        } else {
                                die("Error: Missing parameter value");
                        }
                }               else if (!strcmp(my_argv[i], "--wrap") || 
!strcmp(my_argv[i], "-w")) {
                        if (i + 1 >= (int)my_argv.size()) {
                                die("Error: Missing parameters");
                        }
                        if (!strcmp(my_argv[i + 1], "none") || 
!strcmp(my_argv[i + 1], "open")) ++i;
                        else if (!strcmp(my_argv[i + 1], "horizontal") || 
!strcmp(my_argv[i + 1], "h")) { wrap = 1; ++i; } else if (!strcmp(my_argv[i + 
1], "vertical") || !strcmp(my_argv[i + 1], "v")) { wrap = 2; ++i; } else if 
(!strcmp(my_argv[i + 1], "both") || !strcmp(my_argv[i + 1], "hv")) { wrap = 3; 
++i; } else wrap = 1;
                }               else if (!strcmp(my_argv[i], 
"--cache-threshold")) {
                        if (i + 1 >= (int)my_argv.size()) {
                                die("Error: Missing parameters");
                        }
                        ++i;
                        int shift = 0;
                        int n = 0;
                        size_t len = strlen(my_argv[i]);
                        size_t threshold;
                        sscanf_s(my_argv[i], "%zu%n", &threshold, &n);
                        if (n != len) {
                                if (n == len - 1) {
                                        switch (my_argv[i][len - 1]) {
                                                case 'k':
                                                case 'K': shift = 10; break;
                                                case 'm':
                                                case 'M': shift = 20; break;
                                                case 'g':
                                                case 'G': shift = 30; break;
                                                default: die("Error: Bad 
--cache-threshold parameter");
                                        }
                                        threshold <<= shift;
                                } else {
                                        die("Error: Bad --cache-threshold 
parameter");
                                }
                        }
                        MapAlloc::CacheThreshold(threshold);
                }               else if (!strcmp(my_argv[i], "--nomask") || 
!strcmp(my_argv[i], "--no-mask")) no_mask = true;
                else if (!strcmp(my_argv[i], "--timing") || !strcmp(my_argv[i], 
"--timings")) timing = true;
                else if (!strcmp(my_argv[i], "--bigtiff"))   big_tiff = true;
                else if (!strcmp(my_argv[i], "--bgr"))       bgr = true;
                else if (!strcmp(my_argv[i], "--wideblend")) wideblend = true;
                else if (!strcmp(my_argv[i], "--reverse"))   reverse = true;
                else if (!strcmp(my_argv[i], "--gamma"))     gamma = true;
                else if (!strcmp(my_argv[i], "--no-dither") || 
!strcmp(my_argv[i], "--nodither")) dither = false;
                //              else if (!strcmp(my_argv[i], "--force"))     
force_coverage = true;
                else if (!strncmp(my_argv[i], "-f", 2)) Output(0, "ignoring 
Enblend option -f\n");
                else if (!strcmp(my_argv[i], "-a")) Output(0, "ignoring Enblend 
option -a\n");
                else if (!strcmp(my_argv[i], "--no-ciecam")) Output(0, 
"ignoring Enblend option --no-ciecam\n");
                else if (!strcmp(my_argv[i], "--primary-seam-generator")) {
                        Output(0, "ignoring Enblend option 
--primary-seam-generator\n");
                        ++i;
                }

                else if (!strcmp(my_argv[i], "--compression")) {
                        if (++i < (int)my_argv.size()) {
                                if (strcmp(my_argv[i], "0") == 0) jpeg_quality 
= 0;
                                else if (atoi(my_argv[i]) > 0) jpeg_quality = 
atoi(my_argv[i]);
                                else if (_stricmp(my_argv[i], "lzw") == 0) 
compression = COMPRESSION_LZW;
                                else if (_stricmp(my_argv[i], "packbits") == 0) 
compression = COMPRESSION_PACKBITS;
                                //                              else if 
(_stricmp(my_argv[i], "deflate") == 0) compression = COMPRESSION_DEFLATE;
                                else if (_stricmp(my_argv[i], "none") == 0) 
compression = COMPRESSION_NONE;
                                else die("Error: Unknown compression codec %s", 
my_argv[i]);
                        } else {
                                die("Error: Missing parameter value");
                        }
                } else if (!strcmp(my_argv[i], "-v") || !strcmp(my_argv[i], 
"--verbose")) ++verbosity;
                else if (!strcmp(my_argv[i], "-q") || !strcmp(my_argv[i], 
"--quiet")) --verbosity;
                else if ((!strcmp(my_argv[i], "--saveseams") || 
!strcmp(my_argv[i], "--save-seams")) && i < (int)my_argv.size() - 1) 
seamsave_filename = my_argv[++i];
                else if ((!strcmp(my_argv[i], "--loadseams") || 
!strcmp(my_argv[i], "--load-seams")) && i < (int)my_argv.size() - 1) 
seamload_filename = my_argv[++i];
                else if ((!strcmp(my_argv[i], "--savexor") || 
!strcmp(my_argv[i], "--save-xor")) && i < (int)my_argv.size() - 1) xor_filename 
= my_argv[++i];
                else if (!strcmp(my_argv[i], "--tempdir") || 
!strcmp(my_argv[i], "--tmpdir") && i < (int)my_argv.size() - 1) 
MapAlloc::SetTmpdir(my_argv[++i]);
                else if (!strcmp(my_argv[i], "--all-threads")) all_threads = 
true;
                else if (!strcmp(my_argv[i], "-o") || !strcmp(my_argv[i], 
"--output")) {
                        if (++i < (int)my_argv.size()) {
                                output_filename = my_argv[i];
                                char* ext = strrchr(output_filename, '.');

                                if (!ext) {
                                        die("Error: Unknown output filetype");
                                }

                                ++ext;
                                if (!(_stricmp(ext, "jpg") && _stricmp(ext, 
"jpeg"))) {
                                        output_type = ImageType::MB_JPEG;
                                        if (jpeg_quality == -1) jpeg_quality = 
75;
                                } else if (!(_stricmp(ext, "tif") && 
_stricmp(ext, "tiff"))) {
                                        output_type = ImageType::MB_TIFF;
                                } else if (!_stricmp(ext, "png")) {
                                        output_type = ImageType::MB_PNG;
                                } else {
                                        die("Error: Unknown file extension");
                                }

                                ++i;
                                break;
                        }
                } else if (!strcmp(my_argv[i], "--no-output")) {
                        ++i;
                        break;
                } else {
                        die("Error: Unknown argument \"%s\"", my_argv[i]);
                }
        }

        if (compression != -1) {
                if (output_type != ImageType::MB_TIFF) {
                        Output(0, "Warning: non-TIFF output; ignoring TIFF 
compression setting\n");
                }
        } else if (output_type == ImageType::MB_TIFF) {
                compression = COMPRESSION_LZW;
        }

        if (jpeg_quality != -1 && output_type != ImageType::MB_JPEG && 
output_type != ImageType::MB_PNG) {
                Output(0, "Warning: non-JPEG/PNG output; ignoring compression 
quality setting\n");
        }

        if ((jpeg_quality < -1 || jpeg_quality > 9) && output_type == 
ImageType::MB_PNG) {
                die("Error: Bad PNG compression quality setting\n");
        }

        if (output_type == ImageType::MB_NONE && !seamsave_filename) 
die("Error: No output file specified");
        if (seamload_filename && seamsave_filename) die("Error: Cannot load and 
save seams at the same time");
        if (wrap == 3) die("Error: Wrapping in both directions is not currently 
supported");

        if (!strcmp(my_argv[i], "--")) ++i;

/***********************************************************************
* Push remaining arguments to images vector
***********************************************************************/
        int x, y, n;

        while (i < (int)my_argv.size()) {
                if (images.size()) {
                        n = 0;
                        sscanf_s(my_argv[i], "%d,%d%n", &x, &y, &n);
                        if (!my_argv[i][n]) {
                                images.back()->xpos_add = x;
                                images.back()->ypos_add = y;
                                i++;
                                continue;
                        }
                }
                images.push_back(new Image(my_argv[i++]));
        }

        int n_images = (int)images.size();

        if (n_images == 0) die("Error: No input files specified");
        if (seamsave_filename && n_images > 256) { seamsave_filename = NULL; 
Output(0, "Warning: seam saving not possible with more than 256 images"); }
        if (seamload_filename && n_images > 256) { seamload_filename = NULL; 
Output(0, "Warning: seam loading not possible with more than 256 images"); }
        if (xor_filename && n_images > 255) { xor_filename = NULL; Output(0, 
"Warning: XOR map saving not possible with more than 255 images"); }

/***********************************************************************
* Print banner
***********************************************************************/
        Output(1, "\n");
        Output(1, "Multiblend v2.0.0 (c) 2021 David Horman        
http://horman.net/multiblend/\n";);
        Output(1, 
"----------------------------------------------------------------------------\n");

        Threadpool* threadpool = Threadpool::GetInstance(all_threads ? 2 : 0);

/***********************************************************************
************************************************************************
* Open output
************************************************************************
***********************************************************************/
        switch (output_type) {
                case ImageType::MB_TIFF: {
                        if (!big_tiff) tiff_file = TIFFOpen(output_filename, 
"w"); else tiff_file = TIFFOpen(output_filename, "w8");
                        if (!tiff_file) die("Error: Could not open output 
file");
                } break;
                case ImageType::MB_JPEG: {
                        if (output_bpp == 16) die("Error: 16bpp output is 
incompatible with JPEG output");
                        fopen_s(&jpeg_file, output_filename, "wb");
                        if (!jpeg_file) die("Error: Could not open output 
file");
                } break;
                case ImageType::MB_PNG: {
                        fopen_s(&jpeg_file, output_filename, "wb");
                        if (!jpeg_file) die("Error: Could not open output 
file");
                } break;
        }

/***********************************************************************
************************************************************************
* Process images
************************************************************************
***********************************************************************/
        timer.Start();

/***********************************************************************
* Open images to get prelimary info
***********************************************************************/
        size_t untrimmed_bytes = 0;

        for (i = 0; i < n_images; ++i) {
                images[i]->Open();
                untrimmed_bytes = std::max(untrimmed_bytes, 
images[i]->untrimmed_bytes);
        }

/***********************************************************************
* Check paramters, display warnings
***********************************************************************/
        for (i = 1; i < n_images; ++i) {
                if (images[i]->tiff_xres != images[0]->tiff_xres || 
images[i]->tiff_yres != images[0]->tiff_yres) {
                        Output(0, "Warning: TIFF resolution mismatch (%f %f/%f 
%f)\n", images[0]->tiff_xres, images[0]->tiff_yres, images[i]->tiff_xres, 
images[i]->tiff_yres);
                }
        }

        for (i = 0; i < n_images; ++i) {
                if (output_bpp == 0 && images[i]->bpp == 16) output_bpp = 16;
                if (images[i]->bpp != images[0]->bpp) {
                        die("Error: mixture of 8bpp and 16bpp images detected 
(not currently handled)\n");
                }
        }

        if (output_bpp == 0) output_bpp = 8;
        else if (output_bpp == 16 && output_type == ImageType::MB_JPEG) {
                Output(0, "Warning: 8bpp output forced by JPEG output\n");
                output_bpp = 8;
        }

/***********************************************************************
* Allocate working space for reading/trimming/extraction
***********************************************************************/
        void* untrimmed_data = MapAlloc::Alloc(untrimmed_bytes);

/***********************************************************************
* Read/trim/extract
***********************************************************************/
        for (i = j = 0; i < n_images; ++i) {
                try {
                        images[i]->Read(untrimmed_data, gamma);
                } catch (char* e) {
                        printf("\n\n");
                        printf("%s\n", e);
                        exit(EXIT_FAILURE);
                }
                if(0 == images[i]->width || 0 == images[i]->height)
                        delete images[i];
                else
                        images[j++] = images[i];
        }
        if(j < n_images)
        {
                if(0 == j)
                        die("\nError: All input images have only transparent 
pixels");
                n_images = j;
                images.resize(j);
        }

/***********************************************************************
* Clean up
***********************************************************************/
        MapAlloc::Free(untrimmed_data);

/***********************************************************************
* Tighten
***********************************************************************/
        int min_xpos = 0x7fffffff;
        int min_ypos = 0x7fffffff;
        width = 0;
        height = 0;

        for (i = 0; i < n_images; ++i) {
                min_xpos = std::min(min_xpos, images[i]->xpos);
                min_ypos = std::min(min_ypos, images[i]->ypos);
        }

        for (i = 0; i < n_images; ++i) {
                images[i]->xpos -= min_xpos;
                images[i]->ypos -= min_ypos;
                width = std::max(width, images[i]->xpos + images[i]->width);
                height = std::max(height, images[i]->ypos + images[i]->height);
        }

        images_time = timer.Read();

/***********************************************************************
* Determine number of levels
***********************************************************************/
        int blend_wh;
        int blend_levels;

        if (!fixed_levels) {
                if (!wideblend) {
                        std::vector<int> widths;
                        std::vector<int> heights;

                        for (auto image : images) {
                                widths.push_back(image->width);
                                heights.push_back(image->height);
                        }

                        std::sort(widths.begin(), widths.end());
                        std::sort(heights.begin(), heights.end());

                        size_t halfway = (widths.size() - 1) >> 1;

                        blend_wh = std::max(
                                widths.size() & 1 ? widths[halfway] : 
(widths[halfway] + widths[halfway + 1] + 1) >> 1,
                                heights.size() & 1 ? heights[halfway] : 
(heights[halfway] + heights[halfway + 1] + 1) >> 1
                        );
                } else {
                        blend_wh = (std::max)(width, height);
                }

                blend_levels = (int)floor(log2(blend_wh + 4.0f) - 1);
                if (wideblend) blend_levels++;
        } else {
                blend_levels = fixed_levels;
        }

        blend_levels += add_levels;

        if (n_images == 1) {
                blend_levels = 0;
                Output(1, "\n%d x %d, %d bpp\n\n", width, height, output_bpp);
        } else {
                Output(1, "\n%d x %d, %d levels, %d bpp\n\n", width, height, 
blend_levels, output_bpp);
        }
        
/***********************************************************************
************************************************************************
* Seaming
************************************************************************
***********************************************************************/
        timer.Start();

        Output(1, "Seaming");
        switch (((!!seamsave_filename) << 1) | !!xor_filename) {
                case 1: Output(1, " (saving XOR map)"); break;
                case 2: Output(1, " (saving seam map)"); break;
                case 3: Output(1, " (saving XOR and seam maps)"); break;
        }
        Output(1, "...\n");

        int min_count;
        int xor_count;
        int xor_image;
        uint64_t utemp;
        int stop;

        uint64_t best;
        uint64_t a, b, c, d;

#define DT_MAX 0x9000000000000000
        uint64_t* prev_line = NULL;
        uint64_t* this_line = NULL;
        bool last_pixel = false;
        bool arbitrary_seam = false;

        Flex* seam_flex = new Flex(width, height);
        int max_queue = 0;

/***********************************************************************
* Backward distance transform
***********************************************************************/
        int n_threads = std::max(2, threadpool->GetNThreads());
        uint64_t** thread_lines = new uint64_t*[n_threads];

        if (!seamload_filename) {
                std::mutex* flex_mutex_p = new std::mutex;
                std::condition_variable* flex_cond_p = new 
std::condition_variable;

                uint8_t** thread_comp_lines = new uint8_t*[n_threads];

                for (i = 0; i < n_threads; ++i) {
                        thread_lines[i] = new uint64_t[width];
                        thread_comp_lines[i] = new uint8_t[width];
                }

// set all image masks to bottom right
                for (i = 0; i < n_images; ++i) {
                        images[i]->tiff_mask->End();
                }

                for (y = height - 1; y >= 0; --y) {
                        int t = y % n_threads;
                        this_line = thread_lines[t];
                        uint8_t* comp = thread_comp_lines[t];

// set initial image mask states
                        for (i = 0; i < n_images; ++i) {
                                images[i]->mask_state = 0x8000000000000000;
                                if (y >= images[i]->ypos && y < images[i]->ypos 
+ images[i]->height) {
                                        images[i]->mask_count = width - 
(images[i]->xpos + images[i]->width);
                                        images[i]->mask_limit = images[i]->xpos;
                                } else {
                                        images[i]->mask_count = width;
                                        images[i]->mask_limit = width;
                                }
                        }

                        x = width - 1;

                        { // make sure the last compression thread to use this 
chunk of memory is finished
                                std::unique_lock<std::mutex> 
mlock(*flex_mutex_p);
                                flex_cond_p->wait(mlock, [=] { return 
seam_flex->y > (height - 1) - y - n_threads; });
                        }

                        while (x >= 0) {
                                min_count = x + 1;
                                xor_count = 0;

// update image mask states
                                for (i = 0; i < n_images; ++i) {
                                        if (!images[i]->mask_count) {
                                                if (x >= images[i]->mask_limit) 
{
                                                        utemp = 
images[i]->tiff_mask->ReadBackwards32();
                                                        images[i]->mask_state = 
((~utemp) << 32) & 0x8000000000000000;
                                                        images[i]->mask_count = 
utemp & 0x7fffffff;
                                                } else {
                                                        images[i]->mask_state = 
0x8000000000000000;
                                                        images[i]->mask_count = 
min_count;
                                                }
                                        }

                                        if (images[i]->mask_count < min_count) 
min_count = images[i]->mask_count;
                                        if (!images[i]->mask_state) { // 
mask_state is inverted
                                                ++xor_count;
                                                xor_image = i;
                                        }
                                }

                                stop = x - min_count;

                                if (xor_count == 1) {
                                        images[xor_image]->seam_present = true;
                                        while (x > stop) this_line[x--] = 
xor_image;
                                } else {
                                        if (y == height - 1) { // bottom row
                                                if (x == width - 1) { // first 
pixel(s)
                                                        while (x > stop) 
this_line[x--] = DT_MAX; // max
                                                } else {
                                                        utemp = this_line[x + 
1];
                                                        utemp = MASKVAL(utemp);
                                                        while (x > stop) {
                                                                utemp += 
0x300000000;
                                                                this_line[x--] 
= utemp; // was min(temp, DT_MAX) but this is unlikely to happen
                                                        }
                                                }
                                        } else { // other rows
                                                if (x == width - 1) { // first 
pixel(s)
                                                        utemp = prev_line[x - 
1] + 0x400000000;
                                                        a = MASKVAL(utemp);

                                                        utemp = prev_line[x] + 
0x300000000;
                                                        b = MASKVAL(utemp);

                                                        d = a < b ? a : b;

                                                        this_line[x--] = d;

                                                        if (x == stop) {
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        
images[i]->mask_count -= min_count;
                                                                }
                                                                continue;
                                                        }

                                                        c = b + 0x100000000;
                                                        b = a - 0x100000000;
                                                        d += 0x300000000;
                                                } else {
                                                        utemp = prev_line[x] + 
0x300000000;
                                                        b = MASKVAL(utemp);

                                                        utemp = prev_line[x + 
1] + 0x400000000;
                                                        c = MASKVAL(utemp);

                                                        utemp = this_line[x + 
1] + 0x300000000;
                                                        d = MASKVAL(utemp);
                                                }

                                                if (stop == -1) {
                                                        stop = 0;
                                                        last_pixel = true;
                                                }

                                                while (x > stop) {
                                                        utemp = prev_line[x - 
1] + 0x400000000;
                                                        a = MASKVAL(utemp);

                                                        if (a < d) d = a;
                                                        if (b < d) d = b;
                                                        if (c < d) d = c;
                                                        
                                                        this_line[x--] = d;

                                                        c = b + 0x100000000;
                                                        b = a - 0x100000000;
                                                        d += 0x300000000;
                                                }

                                                if (last_pixel) {
                                                        // d is the new "best" 
to compare against
                                                        if (b < d) d = b;
                                                        if (c < d) d = c;

                                                        this_line[x--] = d;

                                                        last_pixel = false;
                                                }
                                        }
                                }

                                for (i = 0; i < n_images; ++i) {
                                        images[i]->mask_count -= min_count;
                                }
                        }

                        if (y) {
                                threadpool->Queue([=] {
                                        int p = CompressSeamLine(this_line, 
comp, width);
                                        if (p > width) {
                                                printf("bad p: %d at line %d", 
p, y);
                                                exit(0);
                                        }

                                        {
                                                std::unique_lock<std::mutex> 
mlock(*flex_mutex_p);
                                                flex_cond_p->wait(mlock, [=] { 
return seam_flex->y == (height - 1) - y; });
                                                seam_flex->Copy(comp, p);
                                                seam_flex->NextLine();
                                        }
                                        flex_cond_p->notify_all();
                                });
                        }

                        prev_line = this_line;
                } // end of row loop

                threadpool->Wait();

                for (i = 0; i < n_images; ++i) {
                        if (!images[i]->seam_present) {
                                Output(1, "Warning: %s is fully obscured by 
other images\n", images[i]->filename);
                        }
                }

                for (i = 0; i < n_threads; ++i) {
                        if (i >= 2) delete[] thread_lines[i];
                        delete[] thread_comp_lines[i];
                }

                delete[] thread_comp_lines;
                delete flex_mutex_p;
                delete flex_cond_p;
        } else { // if seamload_filename:
                for (i = 0; i < n_images; ++i) {
                        images[i]->tiff_mask->Start();
                }
        }

// create top level masks
        for (i = 0; i < n_images; ++i) {
                images[i]->masks.push_back(new Flex(width, height));
        }

        Pnger* xor_map = xor_filename ? new Pnger(xor_filename, "XOR map", 
width, height, PNG_COLOR_TYPE_PALETTE) : NULL;
        Pnger* seam_map = seamsave_filename ? new Pnger(seamsave_filename, 
"Seam map", width, height, PNG_COLOR_TYPE_PALETTE) : NULL;

/***********************************************************************
* Forward distance transform
***********************************************************************/
        int current_count = 0;
        int64 current_step;
        uint64_t dt_val;

        prev_line = thread_lines[1];

        uint64_t total_pixels = 0;
        uint64_t channel_totals[3] = { 0 };

        Flex full_mask(width, height);
        Flex xor_mask(width, height);

        bool alpha = false;

        for (y = 0; y < height; ++y) {
                for (i = 0; i < n_images; ++i) {
                        images[i]->mask_state = 0x8000000000000000;
                        if (y >= images[i]->ypos && y < images[i]->ypos + 
images[i]->height) {
                                images[i]->mask_count = images[i]->xpos;
                                images[i]->mask_limit = images[i]->xpos + 
images[i]->width;
                        } else {
                                images[i]->mask_count = width;
                                images[i]->mask_limit = width;
                        }
                }

                x = 0;
                int mc = 0;
                int prev_i = -1;
                int current_i = -1;
                int best_temp;

                while (x < width) {
                        min_count = width - x;
                        xor_count = 0;

                        for (i = 0; i < n_images; ++i) {
                                if (!images[i]->mask_count) {
                                        if (x < images[i]->mask_limit) {
                                                utemp = 
images[i]->tiff_mask->ReadForwards32();
                                                images[i]->mask_state = 
((~utemp) << 32) & 0x8000000000000000;
                                                images[i]->mask_count = utemp & 
0x7fffffff;
                                        } else {
                                                images[i]->mask_state = 
0x8000000000000000;
                                                images[i]->mask_count = 
min_count;
                                        }
                                }

                                if (images[i]->mask_count < min_count) 
min_count = images[i]->mask_count;
                                if (!images[i]->mask_state) {
                                        ++xor_count;
                                        xor_image = i;
                                }
                        }

                        stop = x + min_count;

                        if (!xor_count) {
                                alpha = true;
                        }
                        full_mask.MaskWrite(min_count, xor_count);
                        xor_mask.MaskWrite(min_count, xor_count == 1);

                        if (xor_count == 1) {
                                if (xor_map) memset(&xor_map->line[x], 
xor_image, min_count);

                                size_t p = (y - images[xor_image]->ypos) * 
images[xor_image]->width + (x - images[xor_image]->xpos);

                                int total_count = min_count;
                                total_pixels += total_count;
                                if (gamma) {
                                        switch (images[xor_image]->bpp) {
                                                case 8: {
                                                        uint16_t v;
                                                        while (total_count--) {
                                                                v = 
((uint8_t*)images[xor_image]->channels[0]->data)[p];
                                                                
channel_totals[0] += v * v;
                                                                v = 
((uint8_t*)images[xor_image]->channels[1]->data)[p];
                                                                
channel_totals[1] += v * v;
                                                                v = 
((uint8_t*)images[xor_image]->channels[2]->data)[p];
                                                                
channel_totals[2] += v * v;
                                                                ++p;
                                                        }
                                                } break;
                                                case 16: {
                                                        uint32_t v;
                                                        while (total_count--) {
                                                                v = 
((uint16_t*)images[xor_image]->channels[0]->data)[p];
                                                                
channel_totals[0] += v * v;
                                                                v = 
((uint16_t*)images[xor_image]->channels[1]->data)[p];
                                                                
channel_totals[1] += v * v;
                                                                v = 
((uint16_t*)images[xor_image]->channels[2]->data)[p];
                                                                
channel_totals[2] += v * v;
                                                                ++p;
                                                        }
                                                } break;
                                        }
                                } else {
                                        switch (images[xor_image]->bpp) {
                                                case 8: {
                                                        while (total_count--) {
                                                                
channel_totals[0] += ((uint8_t*)images[xor_image]->channels[0]->data)[p];
                                                                
channel_totals[1] += ((uint8_t*)images[xor_image]->channels[1]->data)[p];
                                                                
channel_totals[2] += ((uint8_t*)images[xor_image]->channels[2]->data)[p];
                                                                ++p;
                                                        }
                                                } break;
                                                case 16: {
                                                        while (total_count--) {
                                                                
channel_totals[0] += ((uint16_t*)images[xor_image]->channels[0]->data)[p];
                                                                
channel_totals[1] += ((uint16_t*)images[xor_image]->channels[1]->data)[p];
                                                                
channel_totals[2] += ((uint16_t*)images[xor_image]->channels[2]->data)[p];
                                                                ++p;
                                                        }
                                                } break;
                                        }
                                }

                                if (!seamload_filename) {
                                        RECORD(xor_image, min_count);
                                        while (x < stop) {
                                                this_line[x++] = xor_image;
                                        }
                                } else {
                                        x = stop;
                                }

                                best = xor_image;
                        } else {
                                if (xor_map) memset(&xor_map->line[x], 0xff, 
min_count);

                                if (!seamload_filename) {
                                        if (y == 0) {
                                                // top row
                                                while (x < stop) {
                                                        best = this_line[x];

                                                        if (x > 0) {
                                                                utemp = 
this_line[x - 1] + 0x300000000;
                                                                d = 
MASKVAL(utemp);

                                                                if (d < best) 
best = d;
                                                        }

                                                        if (best & 
0x8000000000000000 && xor_count) {
                                                                arbitrary_seam 
= true;
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        if 
(!images[i]->mask_state) {
                                                                                
best = 0x8000000000000000 | i;
                                                                                
if (!reverse) break;
                                                                        }
                                                                }
                                                        }

                                                        best_temp = best & 
0xffffffff;
                                                        RECORD(best_temp, 1);
                                                        this_line[x++] = best;
                                                }
                                        } else {
                                                // other rows
                                                if (x == 0) {
                                                        SEAM_DT;
                                                        best = dt_val;

                                                        utemp = *prev_line + 
0x300000000;
                                                        b = MASKVAL(utemp);
                                                        if (b < best) best = b;

                                                        utemp = prev_line[1] + 
0x400000000;
                                                        c = MASKVAL(utemp);
                                                        if (c < best) best = c;

                                                        if (best & 
0x8000000000000000 && xor_count) {
                                                                arbitrary_seam 
= true;
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        if 
(!images[i]->mask_state) {
                                                                                
best = 0x8000000000000000 | i;
                                                                                
if (!reverse) break;
                                                                        }
                                                                }
                                                        }

                                                        best_temp = best & 
0xffffffff;
                                                        RECORD(best_temp, 1);
                                                        this_line[x++] = best;

                                                        if (x == stop) {
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        
images[i]->mask_count -= min_count;
                                                                }
                                                                continue;
                                                        }

                                                        a = b + 0x100000000;
                                                        b = c - 0x100000000;
                                                } else {
                                                        utemp = prev_line[x - 
1] + 0x400000000;
                                                        a = MASKVAL(utemp);

                                                        utemp = prev_line[x] + 
0x300000000;
                                                        b = MASKVAL(utemp);
                                                }

                                                utemp = best + 0x300000000;
                                                d = MASKVAL(utemp);

                                                if (stop == width) {
                                                        stop--;
                                                        last_pixel = true;
                                                }

                                                while (x < stop) {
                                                        utemp = prev_line[x + 
1] + 0x400000000;
                                                        c = MASKVAL(utemp);

                                                        SEAM_DT;
                                                        best = dt_val;

                                                        if (a < best) best = a;
                                                        if (b < best) best = b;
                                                        if (c < best) best = c;
                                                        if (d < best) best = d;

                                                        if (best & 
0x8000000000000000 && xor_count) {
                                                                arbitrary_seam 
= true;
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        if 
(!images[i]->mask_state) {
                                                                                
best = 0x8000000000000000 | i;
                                                                                
if (!reverse) break;
                                                                        }
                                                                }
                                                        }

                                                        best_temp = best & 
0xffffffff;
                                                        RECORD(best_temp, 1);
                                                        this_line[x++] = best; 
// best;

                                                        a = b + 0x100000000;
                                                        b = c - 0x100000000;
                                                        d = best + 0x300000000;
                                                }

                                                if (last_pixel) {
                                                        SEAM_DT;
                                                        best = dt_val;

                                                        if (a < best) best = a;
                                                        if (b < best) best = b;
                                                        if (d < best) best = d;

                                                        if (best & 
0x8000000000000000 && xor_count) {
                                                                arbitrary_seam 
= true;
                                                                for (i = 0; i < 
n_images; ++i) {
                                                                        if 
(!images[i]->mask_state) {
                                                                                
best = 0x8000000000000000 | i;
                                                                                
if (!reverse) break;
                                                                        }
                                                                }
                                                        }

                                                        best_temp = best & 
0xffffffff;
                                                        RECORD(best_temp, 1);
                                                        this_line[x++] = best; 
// best;

                                                        last_pixel = false;
                                                }
                                        }
                                } else { // if (seamload_filename)...
                                        x = stop;
                                }
                        }

                        for (i = 0; i < n_images; ++i) {
                                images[i]->mask_count -= min_count;
                        }
                }
                        
                if (!seamload_filename) {
                        RECORD(-1, 0);

                        for (i = 0; i < n_images; ++i) {
                                images[i]->masks[0]->NextLine();
                        }
                }

                full_mask.NextLine();
                xor_mask.NextLine();

                if (xor_map) xor_map->Write();
                if (seam_map) seam_map->Write();

                std::swap(this_line, prev_line);
        }

        if (!seamload_filename) {
                delete[] thread_lines[0];
                delete[] thread_lines[1];
                delete[] thread_lines;
        }

        delete xor_map;
        delete seam_map;

        if (!alpha || output_type == ImageType::MB_JPEG) no_mask = true;

/***********************************************************************
* Seam load
***********************************************************************/
        if (seamload_filename) {
                int png_depth, png_colour;
                png_uint_32 png_width, png_height;
                uint8_t sig[8];
                png_structp png_ptr;
                png_infop info_ptr;
                FILE* f;

                fopen_s(&f, seamload_filename, "rb");
                if (!f) die("Error: Couldn't open seam file");

                size_t r = fread(sig, 1, 8, f); // assignment suppresses g++ 
-Ofast warning
                if (!png_check_sig(sig, 8)) die("Error: Bad PNG signature");

                png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, 
NULL, NULL);
                if (!png_ptr) die("Error: Seam PNG problem");
                info_ptr = png_create_info_struct(png_ptr);
                if (!info_ptr) die("Error: Seam PNG problem");

                png_init_io(png_ptr, f);
                png_set_sig_bytes(png_ptr, 8);
                png_read_info(png_ptr, info_ptr);
                png_get_IHDR(png_ptr, info_ptr, &png_width, &png_height, 
&png_depth, &png_colour, NULL, NULL, NULL);

                if (png_width != width || png_height != png_height) die("Error: 
Seam PNG dimensions don't match workspace");
                if (png_depth != 8 || png_colour != PNG_COLOR_TYPE_PALETTE) 
die("Error: Incorrect seam PNG format");

                png_bytep png_line = (png_bytep)malloc(width);

                for (y = 0; y < height; ++y) {
                        png_read_row(png_ptr, png_line, NULL);

                        int ms = 0;
                        int mc = 0;
                        int prev_i = -1;
                        int current_i = -1;

                        for (x = 0; x < width; ++x) {
                                if (png_line[x] > n_images) die("Error: Bad 
pixel found in seam file: %d,%d", x, y);
                                RECORD(png_line[x], 1);
                        }

                        RECORD(-1, 0);

                        for (i = 0; i < n_images; ++i) {
                                images[i]->masks[0]->NextLine();
                        }

                }

                free(png_line);
        }

        seam_time = timer.Read();

/***********************************************************************
* No output?
***********************************************************************/
        void* output_channels[3] = { NULL, NULL, NULL };

        if (output_type != ImageType::MB_NONE) {

/***********************************************************************
* Shrink masks
***********************************************************************/
                Output(1, "Shrinking masks...\n");

                timer.Start();

                for (i = 0; i < n_images; ++i) {
                        threadpool->Queue([=] {
                                ShrinkMasks(images[i]->masks, blend_levels);
                                });
                }
                threadpool->Wait();

                shrink_mask_time = timer.Read();

/***********************************************************************
* Create shared input pyramids
***********************************************************************/
// wrapping
                std::vector<PyramidWithMasks*> wrap_pyramids;
                int wrap_levels_h = 0;
                int wrap_levels_v = 0;

                if (wrap & 1) {
                        wrap_levels_h = (int)floor(log2((width >> 1) + 4.0f) - 
1);
                        wrap_pyramids.push_back(new PyramidWithMasks(width >> 
1, height, wrap_levels_h, 0, 0, true));
                        wrap_pyramids.push_back(new PyramidWithMasks((width + 
1) >> 1, height, wrap_levels_h, width >> 1, 0, true));
                }

                if (wrap & 2) {
                        wrap_levels_v = (int)floor(log2((height >> 1) + 4.0f) - 
1);
                        wrap_pyramids.push_back(new PyramidWithMasks(width, 
height >> 1, wrap_levels_v, 0, 0, true));
                        wrap_pyramids.push_back(new PyramidWithMasks(width, 
(height + 1) >> 1, wrap_levels_v, 0, height >> 1, true));
                }

// masks
                for (auto& py : wrap_pyramids) {
                        threadpool->Queue([=] {
                                py->masks.push_back(new Flex(width, height));
                                for (int y = 0; y < height; ++y) {
                                        if (y < py->GetY() || y >= py->GetY() + 
py->GetHeight()) {
                                                
py->masks[0]->Write32(0x80000000 | width);
                                        } else {
                                                if (py->GetX()) {
                                                        
py->masks[0]->Write32(0x80000000 | py->GetX());
                                                        
py->masks[0]->Write32(0xc0000000 | py->GetWidth());
                                                } else {
                                                        
py->masks[0]->Write32(0xc0000000 | py->GetWidth());
                                                        if (py->GetWidth() != 
width) py->masks[0]->Write32(0x80000000 | (width - py->GetWidth()));
                                                }
                                        }
                                        py->masks[0]->NextLine();
                                }

                                ShrinkMasks(py->masks, py->GetWidth() == width 
? wrap_levels_v : wrap_levels_h);
                                });
                }

                threadpool->Wait();
// end wrapping

                int total_levels = std::max({ blend_levels, wrap_levels_h, 
wrap_levels_v, 1 });

                for (int i = 0; i < n_images; ++i) {
                        images[i]->pyramid = new Pyramid(images[i]->width, 
images[i]->height, blend_levels, images[i]->xpos, images[i]->ypos, true);
                }

                for (int l = total_levels - 1; l >= 0; --l) {
                        size_t max_bytes = 0;

                        if (l < blend_levels) {
                                for (auto& image : images) {
                                        max_bytes = std::max(max_bytes, 
image->pyramid->GetLevel(l).bytes);
                                }
                        }

                        for (auto& py : wrap_pyramids) {
                                if (l < py->GetNLevels()) max_bytes = 
std::max(max_bytes, py->GetLevel(l).bytes);
                        }

                        float* temp;

                        try {
                                temp = (float*)MapAlloc::Alloc(max_bytes);
                        } catch (char* e) {
                                printf("%s\n", e);
                                exit(EXIT_FAILURE);
                        }

                        if (l < blend_levels) {
                                for (auto& image : images) {
                                        image->pyramid->GetLevel(l).data = temp;
                                }
                        }

                        for (auto& py : wrap_pyramids) {
                                if (l < py->GetNLevels()) py->GetLevel(l).data 
= temp;
                        }
                }

/***********************************************************************
* Create output pyramid
***********************************************************************/
                Pyramid* output_pyramid = NULL;

                output_pyramid = new Pyramid(width, height, total_levels, 0, 0, 
true);

                for (int l = total_levels - 1; l >= 0; --l) {
                        float* temp;

                        try {
                                temp = 
(float*)MapAlloc::Alloc(output_pyramid->GetLevel(l).bytes);
                        } catch (char* e) {
                                printf("%s\n", e);
                                exit(EXIT_FAILURE);
                        }

                        output_pyramid->GetLevel(l).data = temp;
                }

/***********************************************************************
* Blend
***********************************************************************/
                if (n_images == 1) {
                        if (wrap) Output(1, "Wrapping...\n"); else Output(1, 
"Processing...\n");
                } else {
                        if (wrap) Output(1, "Blending/wrapping...\n"); else 
Output(1, "Blending...\n");
                }

                for (int c = 0; c < 3; ++c) {
                        if (n_images > 1) {
                                for (i = 0; i < n_images; ++i) {
                                        timer.Start();

                                        
images[i]->pyramid->Copy((uint8_t*)images[i]->channels[c]->data, 1, 
images[i]->width, gamma, images[i]->bpp);
                                        if (output_bpp != images[i]->bpp) 
images[i]->pyramid->Multiply(0, gamma ? (output_bpp == 8 ? 1.0f / 66049 : 
66049) : (output_bpp == 8 ? 1.0f / 257 : 257));

                                        delete images[i]->channels[c];
                                        images[i]->channels[c] = NULL;

                                        copy_time += timer.Read();

                                        timer.Start();
                                        images[i]->pyramid->Shrink();
                                        shrink_time += timer.Read();

                                        timer.Start();
                                        images[i]->pyramid->Laplace();
                                        laplace_time += timer.Read();

                // blend into output pyramid...

                                        timer.Start();

                                        for (int l = 0; l < blend_levels; ++l) {
                                                auto in_level = 
images[i]->pyramid->GetLevel(l);
                                                auto out_level = 
output_pyramid->GetLevel(l);

                                                int x_offset = (in_level.x - 
out_level.x) >> l;
                                                int y_offset = (in_level.y - 
out_level.y) >> l;

                                                for (int b = 0; b < 
(int)out_level.bands.size() - 1; ++b) {
                                                        int sy = 
out_level.bands[b];
                                                        int ey = 
out_level.bands[b + 1];

                                                        threadpool->Queue([=] {
                                                                for (int y = 
sy; y < ey; ++y) {
                                                                        int 
in_line = y - y_offset;
                                                                        if 
(in_line < 0) in_line = 0; else if (in_line > in_level.height - 1) in_line = 
in_level.height - 1;
                                                                        float* 
input_p = in_level.data + (size_t)in_line * in_level.pitch;
                                                                        float* 
output_p = out_level.data + (size_t)y * out_level.pitch;

                                                                        
CompositeLine(input_p, output_p, i, x_offset, in_level.width, out_level.width, 
out_level.pitch, images[i]->masks[l]->data, images[i]->masks[l]->rows[y]);
                                                                }
                                                                });
                                                }

                                                threadpool->Wait();
                                        }

                                        blend_time += timer.Read();
                                }

                                timer.Start();
                                output_pyramid->Collapse(blend_levels);
                                collapse_time += timer.Read();
                        } else {
                                timer.Start();

                                
output_pyramid->Copy((uint8_t*)images[0]->channels[c]->data, 1, 
images[0]->width, gamma, images[0]->bpp);
                                if (output_bpp != images[0]->bpp) 
output_pyramid->Multiply(0, gamma ? (output_bpp == 8 ? 1.0f / 66049 : 66049) : 
(output_bpp == 8 ? 1.0f / 257 : 257));

                                delete images[0]->channels[c];
                                images[0]->channels[c] = NULL;

                                copy_time += timer.Read();
                        }

/***********************************************************************
* Wrapping
***********************************************************************/
                        if (wrap) {
                                timer.Start();

                                int p = 0;

                                for (int w = 1; w <= 2; ++w) {
                                        if (wrap & w) {

                                                if (w == 1) {
                                                        SwapH(output_pyramid);
                                                } else {
                                                        SwapV(output_pyramid);
                                                }

                                                int wrap_levels = (w == 1) ? 
wrap_levels_h : wrap_levels_v;
                                                for (int wp = 0; wp < 2; ++wp) {
                                                        
wrap_pyramids[p]->Copy((uint8_t*)(output_pyramid->GetData() + 
wrap_pyramids[p]->GetX() + wrap_pyramids[p]->GetY() * 
(int64)output_pyramid->GetPitch()), 1, output_pyramid->GetPitch(), false, 32);
                                                        
wrap_pyramids[p]->Shrink();
                                                        
wrap_pyramids[p]->Laplace();

                                                        for (int l = 0; l < 
wrap_levels; ++l) {
                                                                auto in_level = 
wrap_pyramids[p]->GetLevel(l);
                                                                auto out_level 
= output_pyramid->GetLevel(l);

                                                                int x_offset = 
(in_level.x - out_level.x) >> l;
                                                                int y_offset = 
(in_level.y - out_level.y) >> l;

                                                                for (int b = 0; 
b < (int)out_level.bands.size() - 1; ++b) {
                                                                        int sy 
= out_level.bands[b];
                                                                        int ey 
= out_level.bands[b + 1];

                                                                        
threadpool->Queue([=] {
                                                                                
for (int y = sy; y < ey; ++y) {
                                                                                
        int in_line = y - y_offset;
                                                                                
        if (in_line < 0) in_line = 0; else if (in_line > in_level.height - 1) 
in_line = in_level.height - 1;
                                                                                
        float* input_p = in_level.data + (size_t)in_line * in_level.pitch;
                                                                                
        float* output_p = out_level.data + (size_t)y * out_level.pitch;

                                                                                
        CompositeLine(input_p, output_p, wp + (l == 0), x_offset, 
in_level.width, out_level.width, out_level.pitch, 
wrap_pyramids[p]->masks[l]->data, wrap_pyramids[p]->masks[l]->rows[y]);
                                                                                
}
                                                                                
});
                                                                }

                                                                
threadpool->Wait();
                                                        }
                                                        ++p;
                                                }

                                                
output_pyramid->Collapse(wrap_levels);

                                                if (w == 1) {
                                                        UnswapH(output_pyramid);
                                                } else {
                                                        UnswapV(output_pyramid);
                                                }
                                        } // if (wrap & w)
                                } // w loop

                                wrap_time += timer.Read();
                        } // if (wrap)
// end wrapping

/***********************************************************************
* Offset correction
***********************************************************************/
                        if (total_pixels) {
                                double channel_total = 0; // must be a double
                                float* data = output_pyramid->GetData();
                                xor_mask.Start();

                                for (y = 0; y < height; ++y) {
                                        x = 0;
                                        while (x < width) {
                                                uint32_t v = 
xor_mask.ReadForwards32();
                                                if (v & 0x80000000) {
                                                        v = x + v & 0x7fffffff;
                                                        while (x < (int)v) {
                                                                channel_total 
+= data[x++];
                                                        }
                                                } else {
                                                        x += v;
                                                }
                                        }

                                        data += output_pyramid->GetPitch();
                                }

                                float avg = (float)channel_totals[c] / 
total_pixels;
                                if (output_bpp != images[0]->bpp) {
                                        switch (output_bpp) {
                                                case 8: avg /= 256; break;
                                                case 16: avg *= 256; break;
                                        }
                                }
                                float output_avg = (float)channel_total / 
total_pixels;
                                output_pyramid->Add(avg - output_avg, 1);
                        }

/***********************************************************************
* Output
***********************************************************************/
                        timer.Start();

                        try {
                                output_channels[c] = 
MapAlloc::Alloc(((size_t)width * height) << (output_bpp >> 4));
                        } catch (char* e) {
                                printf("%s\n", e);
                                exit(EXIT_FAILURE);
                        }

                        switch (output_bpp) {
                                case 8: 
output_pyramid->Out((uint8_t*)output_channels[c], width, gamma, dither, true); 
break;
                                case 16: 
output_pyramid->Out((uint16_t*)output_channels[c], width, gamma, dither, true); 
break;
                        }

                        out_time += timer.Read();
                }

/***********************************************************************
* Write
***********************************************************************/
#define ROWS_PER_STRIP 64

                Output(1, "Writing %s...\n", output_filename);

                timer.Start();

                struct jpeg_compress_struct cinfo;
                struct jpeg_error_mgr jerr;

                JSAMPARRAY scanlines = NULL;

                int spp = no_mask ? 3 : 4;

                int bytes_per_pixel = spp << (output_bpp >> 4);
                int bytes_per_row = bytes_per_pixel * width;

                int n_strips = (int)((height + ROWS_PER_STRIP - 1) / 
ROWS_PER_STRIP);
                int remaining = height;
                void* strip = malloc((ROWS_PER_STRIP * (int64)width) * 
bytes_per_pixel);
                void* oc_p[3] = { output_channels[0], output_channels[1], 
output_channels[2] };
                if (bgr) std::swap(oc_p[0], oc_p[2]);

                switch (output_type) {
                        case ImageType::MB_TIFF: {
                                TIFFSetField(tiff_file, TIFFTAG_IMAGEWIDTH, 
width);
                                TIFFSetField(tiff_file, TIFFTAG_IMAGELENGTH, 
height);
                                TIFFSetField(tiff_file, TIFFTAG_COMPRESSION, 
compression);
                                TIFFSetField(tiff_file, TIFFTAG_PLANARCONFIG, 
PLANARCONFIG_CONTIG);
                                TIFFSetField(tiff_file, TIFFTAG_ROWSPERSTRIP, 
ROWS_PER_STRIP);
                                TIFFSetField(tiff_file, TIFFTAG_BITSPERSAMPLE, 
output_bpp);
                                if (no_mask) {
                                        TIFFSetField(tiff_file, 
TIFFTAG_SAMPLESPERPIXEL, 3);
                                } else {
                                        TIFFSetField(tiff_file, 
TIFFTAG_SAMPLESPERPIXEL, 4);
                                        uint16_t out[1] = { 
EXTRASAMPLE_UNASSALPHA };
                                        TIFFSetField(tiff_file, 
TIFFTAG_EXTRASAMPLES, 1, &out);
                                }

                                TIFFSetField(tiff_file, TIFFTAG_PHOTOMETRIC, 
PHOTOMETRIC_RGB);
                                if (images[0]->tiff_xres != -1) { 
TIFFSetField(tiff_file, TIFFTAG_XRESOLUTION, images[0]->tiff_xres); 
TIFFSetField(tiff_file, TIFFTAG_XPOSITION, (float)(min_xpos / 
images[0]->tiff_xres)); }
                                if (images[0]->tiff_yres != -1) { 
TIFFSetField(tiff_file, TIFFTAG_YRESOLUTION, images[0]->tiff_yres); 
TIFFSetField(tiff_file, TIFFTAG_YPOSITION, (float)(min_ypos / 
images[0]->tiff_yres)); }

                                if (images[0]->geotiff.set) {
                                        // if we got a georeferenced input, 
store the geotags in the output
                                        GeoTIFFInfo info(images[0]->geotiff);
                                        info.XGeoRef = min_xpos * 
images[0]->geotiff.XCellRes;
                                        info.YGeoRef = -min_ypos * 
images[0]->geotiff.YCellRes;
                                        Output(1, "Output georef: UL: %f %f, 
pixel size: %f %f\n", info.XGeoRef, info.YGeoRef, info.XCellRes, info.YCellRes);
                                        geotiff_write(tiff_file, &info);
                                }
                        } break;
                        case ImageType::MB_JPEG: {
                                cinfo.err = jpeg_std_error(&jerr);
                                jpeg_create_compress(&cinfo);
                                jpeg_stdio_dest(&cinfo, jpeg_file);

                                cinfo.image_width = width;
                                cinfo.image_height = height;
                                cinfo.input_components = 3;
                                cinfo.in_color_space = JCS_RGB;

                                jpeg_set_defaults(&cinfo);
                                jpeg_set_quality(&cinfo, jpeg_quality, true);
                                jpeg_start_compress(&cinfo, true);
                        } break;
                        case ImageType::MB_PNG: {
                                png_file = new Pnger(output_filename, NULL, 
width, height, no_mask ? PNG_COLOR_TYPE_RGB : PNG_COLOR_TYPE_RGB_ALPHA, 
output_bpp, jpeg_file, jpeg_quality);
                        } break;
                }

                if (output_type == ImageType::MB_PNG || output_type == 
ImageType::MB_JPEG) {
                        scanlines = new JSAMPROW[ROWS_PER_STRIP];
                        for (i = 0; i < ROWS_PER_STRIP; ++i) {
                                scanlines[i] = (JSAMPROW) & ((uint8_t*)strip)[i 
* bytes_per_row];
                        }
                }

                full_mask.Start();

                for (int s = 0; s < n_strips; ++s) {
                        int strip_p = 0;
                        int rows = std::min(remaining, ROWS_PER_STRIP);

                        for (int strip_y = 0; strip_y < rows; ++strip_y) {
                                x = 0;
                                while (x < width) {
                                        uint32_t cur = 
full_mask.ReadForwards32();
                                        if (cur & 0x80000000) {
                                                int lim = x + (cur & 
0x7fffffff);
                                                switch (output_bpp) {
                                                        case 8: {
                                                                while (x < lim) 
{
                                                                        
((uint8_t*)strip)[strip_p++] = ((uint8_t*)(oc_p[0]))[x];
                                                                        
((uint8_t*)strip)[strip_p++] = ((uint8_t*)(oc_p[1]))[x];
                                                                        
((uint8_t*)strip)[strip_p++] = ((uint8_t*)(oc_p[2]))[x];
                                                                        if 
(!no_mask) ((uint8_t*)strip)[strip_p++] = 0xff;
                                                                        ++x;
                                                                }
                                                        } break;
                                                        case 16: {
                                                                while (x < lim) 
{
                                                                        
((uint16_t*)strip)[strip_p++] = ((uint16_t*)(oc_p[0]))[x];
                                                                        
((uint16_t*)strip)[strip_p++] = ((uint16_t*)(oc_p[1]))[x];
                                                                        
((uint16_t*)strip)[strip_p++] = ((uint16_t*)(oc_p[2]))[x];
                                                                        if 
(!no_mask) ((uint16_t*)strip)[strip_p++] = 0xffff;
                                                                        ++x;
                                                                }
                                                        } break;
                                                }
                                        } else {
                                                size_t t = (size_t)cur * 
bytes_per_pixel;
                                                switch (output_bpp) {
                                                        case 8: {
                                                                
ZeroMemory(&((uint8_t*)strip)[strip_p], t);
                                                        } break;
                                                        case 16: {
                                                                
ZeroMemory(&((uint16_t*)strip)[strip_p], t);
                                                        } break;
                                                }
                                                strip_p += cur * spp;
                                                x += cur;
                                        }
                                }

                                switch (output_bpp) {
                                        case 8: {
                                                oc_p[0] = 
&((uint8_t*)(oc_p[0]))[width];
                                                oc_p[1] = 
&((uint8_t*)(oc_p[1]))[width];
                                                oc_p[2] = 
&((uint8_t*)(oc_p[2]))[width];
                                        } break;
                                        case 16: {
                                                oc_p[0] = 
&((uint16_t*)(oc_p[0]))[width];
                                                oc_p[1] = 
&((uint16_t*)(oc_p[1]))[width];
                                                oc_p[2] = 
&((uint16_t*)(oc_p[2]))[width];
                                        } break;
                                }
                        }

                        switch (output_type) {
                                case ImageType::MB_TIFF: {
                                        TIFFWriteEncodedStrip(tiff_file, s, 
strip, rows * (int64)bytes_per_row);
                                } break;
                                case ImageType::MB_JPEG: {
                                        jpeg_write_scanlines(&cinfo, scanlines, 
rows);
                                } break;
                                case ImageType::MB_PNG: {
                                        png_file->WriteRows(scanlines, rows);
                                } break;
                        }

                        remaining -= ROWS_PER_STRIP;
                }

                switch (output_type) {
                        case ImageType::MB_TIFF: {
                                TIFFClose(tiff_file);
                        } break;
                        case ImageType::MB_JPEG: {
                                jpeg_finish_compress(&cinfo);
                                jpeg_destroy_compress(&cinfo);
                                fclose(jpeg_file);
                        } break;
                }

                write_time = timer.Read();
        }

/***********************************************************************
* Timing
***********************************************************************/
        if (timing) {
                printf("\n");
                printf("Images:   %.3fs\n", images_time);
                printf("Seaming:  %.3fs\n", seam_time);
                if (output_type != ImageType::MB_NONE) {
                        printf("Masks:    %.3fs\n", shrink_mask_time);
                        printf("Copy:     %.3fs\n", copy_time);
                        printf("Shrink:   %.3fs\n", shrink_time);
                        printf("Laplace:  %.3fs\n", laplace_time);
                        printf("Blend:    %.3fs\n", blend_time);
                        printf("Collapse: %.3fs\n", collapse_time);
                        if (wrap) printf("Wrapping: %.3fs\n", wrap_time);
                        printf("Output:   %.3fs\n", out_time);
                        printf("Write:    %.3fs\n", write_time);
                }
        }

/***********************************************************************
* Clean up
***********************************************************************/
        if (timing) {
                if (output_type == ImageType::MB_NONE) {
                        timer_all.Report("\nExecution complete. Total execution 
time");
                } else {
                        timer_all.Report("\nBlend complete. Total execution 
time");
                }
        }

        if (argbuffermem) {
                free(argbuffermem);
        }

        for (i = 0; i < n_images; ++i) {
                delete images[i];
        }

        return EXIT_SUCCESS;
}

Reply via email to