covid-sim
Bitmap.cpp
1 #include <errno.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <math.h>
6 
7 #include "BinIO.h"
8 #include "Bitmap.h"
9 #include "Error.h"
10 #include "Param.h"
11 #include "Model.h"
12 
15 
16 #ifdef _WIN32
17 //HAVI avi;
18 static ULONG_PTR m_gdiplusToken;
19 static HBITMAP bmpdib;
20 static CLSID encoderClsid;
21 #endif
22 
23 #ifndef _WIN32
24 #include <sys/stat.h> // for mkdir
25 #endif
26 
27 static unsigned char* bmf, *bmPixels, *bmp;
28 // externs from CovidSim.cpp
29 // TODO: move these to a header files
30 extern char OutFile[1024], OutFileBase[1024];
31 
32 void CaptureBitmap()
33 {
34  int x, y, f, mi;
35  unsigned j;
36  static double logMaxPop;
37  static int fst = 1;
38  double prev;
39 
40  mi = (int)(P.bwidth * P.bheight);
41  if (fst)
42  {
43  fst = 0;
44  int32_t maxPop = 0;
45  for (int i = 0; i < mi; i++) bmPopulation[i] = 0;
46  for (int i = 0; i < P.PopSize; i++)
47  {
48  x = ((int)(Households[Hosts[i].hh].loc_x * P.scalex)) - P.bminx;
49  y = ((int)(Households[Hosts[i].hh].loc_y * P.scaley)) - P.bminy;
50  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
51  {
52  j = y * bmh->width + x;
53  if ((j < bmh->imagesize) && (j >= 0))
54  {
55  bmPopulation[j]++;
56  if (bmPopulation[j] > maxPop) maxPop = bmPopulation[j];
57  }
58  }
59  }
60  logMaxPop = log(1.001 * (double)maxPop);
61  for (int i = 0; i < P.NMC; i++)
62  if (Mcells[i].n > 0)
63  {
64  f = 0;
65  if ((i < P.NMC - 1) && (i / P.get_number_of_micro_cells_high() == (i + 1) / P.get_number_of_micro_cells_high()) && (Mcells[i + 1].n > 0) && ((Mcells[i].country != Mcells[i + 1].country)
66  || ((P.DoAdunitBoundaryOutput) && ((AdUnits[Mcells[i].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor != (AdUnits[Mcells[i + 1].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor)))) f = 1;
67  if ((i > 0) && (i / P.get_number_of_micro_cells_high() == (i - 1) / P.get_number_of_micro_cells_high()) && (Mcells[i - 1].n > 0) && (Mcells[i].country != Mcells[i - 1].country)) f = 1;
68  if ((i < P.NMC - P.get_number_of_micro_cells_high()) && (Mcells[i + P.get_number_of_micro_cells_high()].n > 0) && ((Mcells[i].country != Mcells[i + P.get_number_of_micro_cells_high()].country)
69  || ((P.DoAdunitBoundaryOutput) && ((AdUnits[Mcells[i].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor != (AdUnits[Mcells[i + P.get_number_of_micro_cells_high()].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor)))) f = 1;
70  if ((i >= P.get_number_of_micro_cells_high()) && (Mcells[i - P.get_number_of_micro_cells_high()].n > 0) && (Mcells[i].country != Mcells[i - P.get_number_of_micro_cells_high()].country)) f = 1;
71  if (f)
72  {
73  x = (int)(P.in_microcells_.width_ * (((double)(i / P.get_number_of_micro_cells_high())) + 0.5) * P.scalex) - P.bminx;
74  y = (int)(P.in_microcells_.height_ * (((double)(i % P.get_number_of_micro_cells_high())) + 0.5) * P.scaley) - P.bminy;
75  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
76  {
77  j = y * bmh->width + x;
78  if ((j < bmh->imagesize) && (j >= 0)) bmPopulation[j] = -1;
79  }
80  }
81  }
82  for (int i = 0; i < P.bwidth / 2; i++)
83  {
84  prev = floor(3.99999 * ((double)i) * BWCOLS / ((double)P.bwidth) * 2);
85  f = ((int)prev);
86  for (j = 0; j < 10; j++)
87  {
88  bmPixels[(j + P.bheight + 5) * bmh->width + P.bwidth / 4 + i] = f;
89  }
90  }
91  }
92 #pragma omp parallel for schedule(static,5000) default(none) \
93  shared(mi, bmPixels, bmPopulation, bmInfected, bmTreated, bmRecovered, logMaxPop)
94  for (int i = 0; i < mi; i++)
95  {
96  if (bmPopulation[i] == -1)
97  bmPixels[i] = BWCOLS - 1; /* black for country boundary */
98  else if (bmInfected[i] > 0)
99  bmPixels[i] = (unsigned char)(BWCOLS + BWCOLS * log((double)bmInfected[i]) / logMaxPop); /* red for infected */
100  else if (bmTreated[i] > 0)
101  bmPixels[i] = (unsigned char)(2 * BWCOLS + BWCOLS * log((double)bmTreated[i]) / logMaxPop); /* blue for treated */
102  else if (bmRecovered[i] > 0)
103  bmPixels[i] = (unsigned char)(3 * BWCOLS + BWCOLS * log((double)bmRecovered[i]) / logMaxPop); /* green for recovered */
104  else if (bmPopulation[i] > 0)
105  bmPixels[i] = (unsigned char)(BWCOLS * log((double)bmPopulation[i]) / logMaxPop); /* grey for just people */
106  else
107  bmPixels[i] = 0;
108  }
109 }
110 
111 void OutputBitmap(int tp)
112 {
113  char buf[3000], OutF[3000];
114  int j = 0;
115  static int cn1 = 0, cn2 = 0, cn3 = 0, cn4 = 0;
116 
117  char *OutBaseName = strrchr(OutFile, '/');
118  char *OutBaseName2 = strrchr(OutFile, '\\');
119  if (OutBaseName2 != nullptr && (OutBaseName == nullptr || OutBaseName2 > OutBaseName)) {
120  OutBaseName = OutBaseName2;
121  }
122  if (OutBaseName == nullptr) {
123  OutBaseName = OutFile;
124  }
125 
126  if (tp == 0)
127  {
128  j = cn1;
129  cn1++;
130  sprintf(OutF, "%s.ge" DIRECTORY_SEPARATOR "%s", OutFile, OutBaseName);
131  }
132  else if (tp == 1)
133  {
134  j = cn2;
135  cn2++;
136  sprintf(OutF, "%s.ge" DIRECTORY_SEPARATOR "Mean.%s", OutFile, OutBaseName);
137  }
138  else if (tp == 2)
139  {
140  j = cn3;
141  cn3++;
142  sprintf(OutF, "%s.ge" DIRECTORY_SEPARATOR "Min.%s", OutFile, OutBaseName);
143  }
144  else if (tp == 3)
145  {
146  j = cn4;
147  cn4++;
148  sprintf(OutF, "%s.ge" DIRECTORY_SEPARATOR "Max.%s", OutFile, OutBaseName);
149  }
150 
151  if (P.BitmapFormat == BF_PNG)
152  {
153 #ifdef IMAGE_MAGICK
154  FILE* dat;
155  using namespace Magick;
156  fprintf(stderr, "\noutputing ImageMagick stuff");
157  sprintf(buf, "%s.bmp", OutF);
158  if (!(dat = fopen(buf, "wb"))) ERR_CRITICAL("Unable to open bitmap file\n");
159  fprintf(dat, "BM");
160  //fwrite_big((void *) &bmf,sizeof(unsigned char),(sizeof(bitmap_header)/sizeof(unsigned char))+bmh->imagesize,dat);
161  fwrite_big((void*)bmf, sizeof(bitmap_header), 1, dat);
162  for (int i = 0; i < bmh->imagesize; i++) fputc(bmPixels[i], dat);
163  fclose(dat);
164  Image bmap(buf);
165  sprintf(buf, "%s.%d.png", OutF, j);
166  ColorRGB white(1.0, 1.0, 1.0);
167  bmap.transparent(white);
168  bmap.write(buf);
169 #elif defined(_WIN32)
170  //Windows specific bitmap manipulation code - could be recoded using LIBGD or another unix graphics library
171  using namespace Gdiplus;
172 
173  wchar_t wbuf[1024];
174  size_t a;
175 
176  //Add new bitmap to AVI
177  //if ((P.OutputBitmap == 1) && (tp == 0)) AddAviFrame(avi, bmpdib, (unsigned char*)(&bmh->palette[0][0]));
178 
179  //This transfers HBITMAP to GDI+ Bitmap object
180  Bitmap* gdip_bmp = Bitmap::FromHBITMAP(bmpdib, NULL);
181  //Now change White in palette (first entry) to be transparent
182  if ((cn1 == 1) && (tp == 0))
183  {
184  static UINT palsize;
185  static ColorPalette* palette;
186  palsize = gdip_bmp->GetPaletteSize();
187  palette = (ColorPalette*)malloc(palsize);
188  if (!palette) ERR_CRITICAL("Unable to allocate palette memory\n");
189  (void)gdip_bmp->GetPalette(palette, palsize);
190  palette->Flags = PaletteFlagsHasAlpha;
191  palette->Entries[0] = 0x00ffffff; // Transparent white
192  gdip_bmp->SetPalette(palette);
193  }
194  //Now save as png
195  sprintf(buf, "%s.%05i.png", OutF, j + 1); //sprintf(buf,"%s.ge" DIRECTORY_SEPARATOR "%s.%05i.png",OutFileBase,OutF,j+1);
196  mbstowcs_s(&a, wbuf, strlen(buf) + 1, buf, _TRUNCATE);
197  gdip_bmp->Save(wbuf, &encoderClsid, NULL);
198  delete gdip_bmp;
199 #else
200  fprintf(stderr, "Do not know how to output PNG\n");
201 #endif
202  }
203  else if (P.BitmapFormat == BF_BMP) {
204  sprintf(buf, "%s.%05i.bmp", OutF, j);
205  FILE* dat;
206  if (!(dat = fopen(buf, "wb"))) {
207  char* errMsg = strerror(errno);
208  if (errMsg == nullptr) {
209  ERR_CRITICAL("strerror failed.\n");
210  }
211  ERR_CRITICAL_FMT("Unable to open bitmap file %s (%d): %s\n", buf, errno, errMsg);
212  }
213  fprintf(dat, "BM");
214  fwrite_big((void*)bmf, sizeof(unsigned char), sizeof(BitmapHeader) / sizeof(unsigned char) + bmh->imagesize, dat);
215  fclose(dat);
216  }
217  else
218  {
219  fprintf(stderr, "Unknown Bitmap format: %d\n", (int)P.BitmapFormat);
220  }
221 }
222 void InitBMHead()
223 {
224  int i, j, k, k2, value;
225 
226  fprintf(stderr, "Initialising bitmap\n");
227  k = P.bwidth * P.bheight2;
228  k2 = sizeof(BitmapHeader) / sizeof(unsigned char);
229 
230  if (!(bmf = (unsigned char*)calloc((size_t)k + k2, sizeof(unsigned char))))
231  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
232  bmPixels = &(bmf[k2]);
233  bmp = &(bmf[12]);
234  bmh = (BitmapHeader*)bmf;
235  bmh->spare = 0;
236  bmh->boffset = 2 + sizeof(BitmapHeader);
237  bmh->headersize = 40; // BITMAPINFOHEADER
238  bmh->width = P.bwidth;
239  bmh->height = P.bheight2;
240  bmh->PlanesAndBitspp = 1 // Number of colour planes; must be 1
241  + (8 << 16); // Colour depth: 8 bits per pixel
242  bmh->compr = 0; // No compression (BI_RGB)
243  bmh->imagesize = bmh->width * bmh->height;
244  bmh->filesize = 2 // "BM"
245  + ((unsigned int) sizeof(BitmapHeader)) // BITMAP_HEADER
246  + bmh->imagesize; // Image data
247  bmh->hres = bmh->vres = (int)(bmh->width * 10); // Resolution, in pixels per metre
248  bmh->colours = BWCOLS * 4; // Number of colours in the palette
249  bmh->impcol = 0; // Every colour is important
250  for (i = 0; i < BWCOLS * 4; i++)
251  bmh->palette[i][3] = 0;
252  for (j = 0; j < BWCOLS; j++)
253  {
254  value = 255 - 255 * j / (BWCOLS - 1);
255  // Shades of gray:
256  bmh->palette[j][0] = bmh->palette[j][1] = bmh->palette[j][2] = (unsigned char)value;
257  // Shades of red:
258  bmh->palette[BWCOLS + j][0] = 0;
259  bmh->palette[BWCOLS + j][1] = 0;
260  bmh->palette[BWCOLS + j][2] = (unsigned char)value;
261  // Shades of blue:
262  bmh->palette[2 * BWCOLS + j][0] = (unsigned char)value;
263  bmh->palette[2 * BWCOLS + j][1] = 0;
264  bmh->palette[2 * BWCOLS + j][2] = 0;
265  // Shades of green:
266  bmh->palette[3 * BWCOLS + j][0] = 0;
267  bmh->palette[3 * BWCOLS + j][1] = (unsigned char)value;
268  bmh->palette[3 * BWCOLS + j][2] = 0;
269  }
270  if (!(bmPopulation = (int32_t*)malloc(bmh->imagesize * sizeof(int32_t))))
271  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
272  if (!(bmInfected = (int32_t*)malloc(bmh->imagesize * sizeof(int32_t))))
273  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
274  if (!(bmRecovered = (int32_t*)malloc(bmh->imagesize * sizeof(int32_t))))
275  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
276  if (!(bmTreated = (int32_t*)malloc(bmh->imagesize * sizeof(int32_t))))
277  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
278 
279  if (P.BitmapFormat == BF_PNG)
280  {
281 #ifdef _WIN32
282  bmpdib = CreateDIBSection(GetDC(NULL), (BITMAPINFO*)bmp, DIB_RGB_COLORS, (void**)&bmPixels, NULL, NULL);
283  Gdiplus::GdiplusStartupInput gdiplusStartupInput;
284  Gdiplus::GdiplusStartup(&m_gdiplusToken, &gdiplusStartupInput, NULL);
285 
286  UINT num = 0; // number of image encoders
287  UINT size = 0; // size of the image encoder array in bytes
288 
289  Gdiplus::ImageCodecInfo* pImageCodecInfo = NULL;
290  Gdiplus::GetImageEncodersSize(&num, &size);
291  if (!(pImageCodecInfo = (Gdiplus::ImageCodecInfo*)(malloc(size))))
292  ERR_CRITICAL("Unable to allocate storage for bitmap\n");
293  Gdiplus::GetImageEncoders(num, size, pImageCodecInfo);
294  for (UINT j = 0; j < num; ++j) {
295  // Visual Studio Analyze incorrectly reports this because it doesn't understand Gdiplus::GetImageEncodersSize()
296  // warning C6385: Reading invalid data from 'pImageCodecInfo': the readable size is 'size' bytes, but '208' bytes may be read.
297 #pragma warning( suppress: 6385 )
298  const WCHAR* type = pImageCodecInfo[j].MimeType;
299  if (wcscmp(type, L"image/png") == 0) {
300  encoderClsid = pImageCodecInfo[j].Clsid;
301  j = num;
302  }
303  }
304  free(pImageCodecInfo);
305 #endif
306  }
307 
308  char buf[1024+3];
309  sprintf(buf, "%s.ge", OutFileBase);
310 #ifdef _WIN32
311  if (!(CreateDirectory(buf, NULL))) fprintf(stderr, "Unable to create directory %s\n", buf);
312 #else
313  if (!(mkdir(buf, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH))) fprintf(stderr, "Unable to create directory %s\n", buf);
314 #endif
315 }
316 
317 void Bitmap_Finalise()
318 {
319  if (P.BitmapFormat == BF_PNG)
320  {
321 #ifdef _WIN32
322  Gdiplus::GdiplusShutdown(m_gdiplusToken);
323 #endif
324  }
325 }
DomainSize in_microcells_
Size of spatial domain in microcells.
Definition: Param.h:108
double height_
The height.
Definition: Param.h:24
double width_
The width.
Definition: Param.h:21
int PopSize
Definition: Param.h:33