gdal_read is not filling out the pixels correctly (?)

[PaulWessel]

In debugging new module grdmix in branch mix-imgs-module I read in a 50x50 all red PNG. GDAL correctly reports 50x50 with 3 bands, and we allocate 50x50x3 = 7500 bytes for it. However, the I->data array returned by GMT_Read_Data, if we loop over the 2500 pixels and print the R/G/B values, yields this result:

Pixel    0 is 255/000/000
Pixel    1 is 255/000/000
Pixel    2 is 255/000/000
Pixel    3 is 255/000/000
...
Pixel  832 is 255/000/000
Pixel  833 is 255/000/000
Pixel  834 is 000/000/000
...
Pixel 2499 is 000/000/000

All pixels after 833 are 0/0/0. So instead of 2500 red pixels, we get 2500/3 ~ 834 red pixels and the rest are all zero. Somewhere in gmt_gdalread there is a division of 3 that causes problems. I note that if I plot the red.png in psimage it comes out correctly, but this may be because postscriptlight detects this as an indexed image with 1 color so it makes a 1-bit image. Hard to trace. However, clearly the 3*size array returned by GMT_Read_Data is flawed. Need @joa-quim to have a look. Attached is red.png.

[joa-quim]

Haven't dive into the code but when I read the red.png in Matlab the output is correct

I = gmtmex('gmtread -Ti red.png')

and it reports it to be layout: 'TCBa'
Doing the same in Julia also plots a correct red image but layout is "BRPa". Maybe you are falling into this mem layout shit?

[PaulWessel]

Can you point me to a place where (GMT or elsewhere) where BRP, TRP etc are explained? It should not matter what the source image is, the I->data returned should be consistently the same, and then output can be set to whatever, but internally it makes no sense to have different layouts. That would mean every program needing to work on images must have all those cases there.
BTW, using ppm on input, grdmix produces the correct output, weather written as ppm or png etc. However, when GDAL writes the result I get correct image but

ERROR 1: mem, band 2: Failed to compute statistics, no valid pixels found in sampling.

[joa-quim]

We invented that:

1st char T(op) or B(ott) to indicate if the array is topdown ot bottomup
2nd char R(ow) or C(olumn) indicates if ordering is row or column major
3rd char B(and), P(ixel) or L(ine) indicated the interleaving

Depending on final usage (PS (pixel interleaving) is different from exporting to Matlab (band interleaving)) the image is read stored directly with the final memory layout. Having an internal fixed storage would mean doing intermediate copies and memory shuflings on potentially big images.

[PaulWessel]

Thanks, we should copy/paste that definition into the code somewhere.

So what do you recommend? That grdmix or other tools that read 1-3 images must for each one have a case-laden switch to rearrange those images to a common layout, do the compositing, then rearrange to another layout that GDAL can write? Seems like a lot of work that has to be repeated everywhere.

Just so I understand, when gmt_gdalread returns it has set the mem_layout string to reflect the actual layout of I->data? Judging from grdimage, I guess gmt_gdalread never strips off the alpha and places it in I->alpha on input since in grdimage we are dealing with possible RGBA quadruplets.

[PaulWessel]

What about the 4th character? Seems to be 'a' or ' '? Why lower case a and not A, or are the other options?

[joa-quim]

Thanks, we should copy/paste that definition into the code somewhere.

It is but always hard to find.

In gdalread we have this

	if (!metadata_only && !prhs->O.mem_layout[0]) {	/* If caller did not ask the layout, assign it to the true one used here. */
		prhs->O.mem_layout[0] = topdown  ? 'T' : 'B';
		prhs->O.mem_layout[1] = rowmajor ? 'R' : 'C';
		prhs->O.mem_layout[2] = do_BIP   ? 'P' : 'B';
		prhs->O.mem_layout[3] = 'a';
	}

It means by default the mem layout is TRB, (I think). The Band interleaving is the easiest to reason about. So my guess is that the best is to leave it like that, and yes the I->data should be arranged according to the layout. At the end gdalwrite should know how to write the computed image with no further layout changes.

Yes, I->alpha should contain the alpha band when it exists. the 'a' in the mem layout code is probably never used. It was added for future use but it's not impossible that we don't need it anymore. Don't even remember if it was supposed to become a 'A' in case we have an alpha channel.

[PaulWessel]

So grdmix calls this up front:

GMT_Set_Default (API, "API_IMAGE_LAYOUT", "TRPa"); /* Set GMT image memory layout */

I then read a PNG. It returns and says layout is TRPa. But the I->data does not match that, given only the first 1/3 of the pixels are set. So what can I do?

[joa-quim]

I would need to debug and to night I don't know when I'll have time but I think you want "TRBa". Band not Pixel interleaving.

[PaulWessel]

Change the initial call to

GMT_Set_Default (API, "API_IMAGE_LAYOUT", "TRB "); /* Set GMT image memory layout */

My dump shows that I->data has the separate bands and now all the data are there. Good.

However, upon calling GMT_Write_Data I get these

grdmix [ERROR]: TRB memory layout is not supported, for now only: T(op)C(ol)B(and) or TRP & BRP

This happens whether I try to reset API_IMAGE_LAYOUT to TRPa or not. Do I mess with the I->header->mem_layout string?

[PaulWessel]

Or do I have to shuffle I->data myself to TRPa, update the I->header->mem_layout to say that, and then write?

[PaulWessel]

Seems you are supporting TRB in input but not output?

[PaulWessel]

FYI, I am working on adding TRB to gmt_gdalwrite.c

[joa-quim]

Sorry, had a long zoom. The idea is that you select the layout, do the processing and save, without any need to the shuffle the image. But it might happen that the writing layout is not implemented. But I think we have functions to convert some of the layouts.

[PaulWessel]

Dont the TRB out would seem simple. I follow your code but build TPR from TRB. First attempt not work of course. Try, try again.
If you know where we have functions let me know - buried in gmt_api.c?

[joa-quim]

I remember that we have stuff to deal with images and that you wrote some of those cryptic codes with bitwise operations but don't remember more.

[joa-quim]

Found the thread GMT_Change_Layout()

[PaulWessel]

Hah. OK, yes, not used anywhere but it is there. So I guess I should call this before writing and request TRP (whatever mode that might be)?

[PaulWessel]

OK, so the required TRB to TRP is not one of the supported ones... So I guess I will add that and try this.

[joa-quim]

It is used by gmtmex_parser

[PaulWessel]

OK, now I get it. I found the bit-stuff you mentioned and all. It is like a goldmine - long forgotten. IT is actually pretty nice (!) if we can add more cases. I will focus on the ability to go from TRB to TRP for now and then do it via this function.

[PaulWessel]

I got GMT_Change_Layout to work. I added a few more cases but have not tested them. I did a few of the TR? cases. I guess TC? are for Matlab. Not sure what does B??.

[stale]

This issue has been automatically marked as stale because it has not had activity in the last 90 days. It will be closed if no further activity occurs within 7 days. Thank you for your contributions.