Add realtime decoding

Change calibration to better suit constant gain
Add -Wextra to Makefile
Add link to model of the PLL
Max number of rows for the image is now increasible

Makefile

@@ -1,7 +1,7 @@
 CC = gcc
 BIN = /usr/bin
 INCLUDES = -I.
-CFLAGS = -O3 -DNDEBUG -Wall $(INCLUDES)
+CFLAGS = -O3 -DNDEBUG -Wall -Wextra $(INCLUDES)
 OBJS = main.o image.o dsp.o filter.o reg.o fcolor.o pngio.o median.o color.o
 
 aptdec:	$(OBJS)

README.md

@@ -51,7 +51,7 @@ To uninstall
 ## Options
 
 ```
--i [r|a|b|c|t|l]
+-i [r|a|b|c|t|m]
 Output image type
 Raw (r), Channel A (a), Channel B (b), False Color (c), Temperature (t) or MCIR (m)
 Default: "ab"
@@ -65,7 +65,7 @@ Satellite number
 For temperature calibration
 Default: "19"
 
--e [t|h]
+-e [t|h|d|p]
 Effects
 Histogram equalise (h), Crop Telemetry (t), Denoise (d) or Precipitation (p)
 Defaults: off
@@ -76,6 +76,9 @@ Map file generated by wxmap
 -c <file>
 Use configuration file for false color generation
 Default: Internal defaults
+
+-r
+Realtime decode. When decoding in realtime it is highly recommended to choose a plain raw image.
 ```
 
 ## Output
@@ -97,12 +100,28 @@ Currently there are 4 available effects:
  - `d` for a median denoise filter
  - `p` for a precipitation overlay
 
-## Example
+## Examples
 
 `aptdec -d images -i ab *.wav`
 
 This will process all `.wav` files in the current directory, generate calibrated channel A and B images and put them in the `images` directory.
 
+`aptdec -e dh -i b audio.wav`
+
+Decode `audio.wav` with denoise and histogram equalisation and save it into the current directory.
+
+## Realtime decoding
+
+As of recently a realtime output was added allowing realtime decoding of images.
+
+```
+mkfifo /tmp/aptaudio
+aptdec /tmp/aptaudio
+sox -t pulseaudio alsa_output.pci-0000_00_1b.0.analog-stereo.monitor -c 1 -t wav /tmp/aptaudio
+```
+
+Perform a realtime decode with the audio being played out of`alsa_output.pci-0000_00_1b.0.analog`.
+
 ## Further reading
 
 [https://noaasis.noaa.gov/NOAASIS/pubs/Users_Guide-Building_Receive_Stations_March_2009.pdf](https://noaasis.noaa.gov/NOAASIS/pubs/Users_Guide-Building_Receive_Stations_March_2009.pdf)  

dsp.c

@@ -97,9 +97,8 @@ static inline double Phase(double I, double Q) {
 }
 
 /* Phase locked loop
- * https://en.wikipedia.org/wiki/Phase-locked_loop
  * https://arachnoid.com/phase_locked_loop/
- * https://simple.wikipedia.org/wiki/Phase-locked_loop
+ * Model of this filter here https://www.desmos.com/calculator/m0uadgkoee
  */
 static double pll(double I, double Q) {
 	// PLL coefficient

image.c

@@ -32,7 +32,7 @@ typedef struct {
 } rgparam_t;
 
 typedef struct {
-	float *prow[3000]; // Row buffers
+	float *prow[MAX_HEIGHT]; // Row buffers
 	int nrow; // Number of rows
 	int chA, chB; // ID of each channel
 	char name[256]; // Stripped filename
@@ -44,6 +44,7 @@ typedef struct {
 	int   satnum; // The satellite number
 	char *map; // Path to a map file
 	char *path; // Output directory
+	int   realtime;
 } options_t;
 
 extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]);
@@ -70,7 +71,6 @@ static double rgcal(float x, rgparam_t *rgpr) {
 
 static double tele[16];
 static double Cs;
-static int nbtele;
 
 void histogramEqualise(float **prow, int nrow, int offset, int width){
 	// Plot histogram
@@ -101,77 +101,59 @@ void histogramEqualise(float **prow, int nrow, int offset, int width){
 }
 
 // Brightness calibrate, including telemetry
-void calibrateBrightness(float **prow, int nrow, int offset, int width, int telestart, rgparam_t regr[30]){
+void calibrateImage(float **prow, int nrow, int offset, int width, rgparam_t regr){
 	offset -= SYNC_WIDTH+SPC_WIDTH;
 
 	for (int n = 0; n < nrow; n++) {
 		float *pixelv = prow[n];
 
 		for (int i = 0; i < width+SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH; i++) {
-			float pv = pixelv[i + offset];
+			float pv = rgcal(pixelv[i + offset], &regr);
 
-			// Blend between the calculated regression curves
-			/* FIXME: this can actually make the image look *worse*
-			 * if the signal has a constant input gain.
-			 */
-			/*int k, kof;
-			k = (n - telestart) / FRAME_LEN;
-			if (k >= nbtele)
-				k = nbtele - 1;
-			kof = (n - telestart) % FRAME_LEN;
-
-			if (kof < 64) {
-				if (k < 1) {*/
-					pv = rgcal(pv, &(regr[4]));
-				/*} else {
-					pv = rgcal(pv, &(regr[k])) * (64 + kof) / FRAME_LEN +
-						 rgcal(pv, &(regr[k - 1])) * (64 - kof) / FRAME_LEN;
-				}
-			} else {
-				if ((k + 1) >= nbtele) {
-					pv = rgcal(pv, &(regr[k]));
-				} else {
-					pv = rgcal(pv, &(regr[k])) * (192 - kof) / FRAME_LEN +
-						 rgcal(pv, &(regr[k + 1])) * (kof - 64) / FRAME_LEN;
-				}
-			}
-*/
-			pv = CLIP(pv, 0, 255);
-			pixelv[i + offset] = pv;
+			pixelv[i + offset] = CLIP(pv, 0, 255);
 		}
 	}
 }
 
+double teleNoise(double wedges[16]){
+	int pattern[9] = { 31, 63, 95, 127, 159, 191, 223, 255, 0 };
+	double noise = 0;
+	for(int i = 0; i < 9; i++)
+		noise += fabs(wedges[i] - (double)pattern[i]);
+	
+	return noise;
+}
+
 // Get telemetry data for thermal calibration/equalization
 int calibrate(float **prow, int nrow, int offset, int width) {
-    double teleline[3000] = { 0.0 };
+    double teleline[MAX_HEIGHT] = { 0.0 };
     double wedge[16];
     rgparam_t regr[30];
-    int n, k;
-    int mtelestart = 0, telestart;
-    int channel = -1;
+    int telestart, mtelestart = 0;
+	int channel = -1;
+
+	// The minimum rows required to decode a full frame
+    if (nrow < 192) {
+		fprintf(stderr, ERR_TELE_ROW);
+		return 0;
+    }
 
 	// Calculate average of a row of telemetry
-    for (n = 0; n < nrow; n++) {
+    for (int n = 0; n < nrow; n++) {
 		float *pixelv = prow[n];
 
 		// Average the center 40px
-		for (int i = 3; i < 43; i++) teleline[n] += pixelv[i + offset + width];
+		for (int i = 3; i < 43; i++)
+			teleline[n] += pixelv[i + offset + width];
 		teleline[n] /= 40.0;
     }
 
-	// The minimum rows required to decode a full frame
-    if (nrow < 192) {
-		fprintf(stderr, ERR_TELE_ROW);
-		return(0);
-    }
-
-	/* Wedge 7 is white and 8 is black, which will have the largest
+	/* Wedge 7 is white and 8 is black, this will have the largest
 	 * difference in brightness, this will always be in the center of
 	 * the frame and can thus be used to find the start of the frame
 	 */
 	double max = 0.0;
-    for (n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) {
+    for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) {
 		float df;
 
 		// (sum 4px below) / (sum 4px above)
@@ -194,68 +176,50 @@ int calibrate(float **prow, int nrow, int offset, int width) {
 		return(0);
     }
 
-	// For each frame
-    for (n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) {
-		float *pixelv = prow[n];
-		int j;
+	// Find the least noisy frame
+	double minNoise = -1;
+	int bestFrame = telestart;
+    for (int n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) {
+		// Turn pixels into wedge values
+		for (int j = 0; j < 16; j++) {
+			wedge[j] = 0.0;
 
-		// Turn each wedge into a value
-		for (j = 0; j < 16; j++) {
 			// Average the middle 6px
-			wedge[j] = 0.0;
-			for (int i = 1; i < 7; i++) wedge[j] += teleline[(j * 8) + n + i];
+			for (int i = 1; i < 7; i++)
+				wedge[j] += teleline[(j * 8) + i + n];
 			wedge[j] /= 6;
 		}
 
-		// Compute regression on the wedges
-		rgcomp(wedge, &(regr[k]));
+		double noise = teleNoise(wedge);
+		if(noise < minNoise || minNoise == -1){
+			minNoise = noise;
+			bestFrame = k;
 
-		// Read the telemetry values from the middle of the image
-		if (k == nrow / (2*FRAME_LEN)) {
-			int l;
-
-			// Equalise
-			for (j = 0; j < 16; j++) tele[j] = rgcal(wedge[j], &(regr[k]));
+			// Compute & apply regression on the wedges
+			rgcomp(wedge, &regr[k]);
+			for (int j = 0; j < 16; j++)
+				tele[j] = rgcal(wedge[j], &regr[k]);
 
 			/* Compare the channel ID wedge to the reference
 			 * wedges, the wedge with the closest match will
 			 * be the channel ID
 			 */
 			float min = -1;
-			for (j = 0; j < 6; j++) {
-				float df;
-
-				df = tele[15] - tele[j];
+			for (int j = 0; j < 6; j++) {
+				float df = tele[15] - tele[j];
 				df *= df;
+
 				if (df < min || min == -1) {
 					channel = j;
 					min = df;
 				}
 			}
-
-			// Cs computation, still have no idea what this does
-			int i;
-			for (Cs = 0.0, i = 0, j = n; j < n + FRAME_LEN; j++) {
-				double csline;
-
-				for (csline = 0.0, l = 3; l < 43; l++)
-					csline += pixelv[l + offset - SPC_WIDTH];
-				
-				csline /= 40.0;
-				if (csline > 50.0) {
-					Cs += csline;
-					i++;
-				}
-			}
-			Cs /= i;
-			Cs = rgcal(Cs, &(regr[k]));
 		}
     }
-    nbtele = k;
 
-	calibrateBrightness(prow, nrow, offset, width, telestart, regr);
+	calibrateImage(prow, nrow, offset, width, regr[bestFrame]);
 
-    return(channel + 1);
+    return channel + 1;
 }
 
 // --- Temperature Calibration --- //

main.c

@@ -38,10 +38,11 @@ typedef struct {
 	int   satnum; // The satellite number
 	char *map; // Path to a map file
 	char *path; // Output directory
+	int   realtime;
 } options_t;
 
 typedef struct {
-	float *prow[3000]; // Row buffers
+	float *prow[MAX_HEIGHT]; // Row buffers
 	int nrow; // Number of rows
 	int chA, chB; // ID of each channel
 	char name[256]; // Stripped filename
@@ -55,6 +56,9 @@ extern int init_dsp(double F);
 extern int readfcconf(char *file);
 extern int readRawImage(char *filename, float **prow, int *nrow);
 extern int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char *chid, char *palette);
+extern void closeWriter();
+extern void pushRow(float *row, int width);
+extern int initWriter(options_t *opts, image_t *img, int width, int height, char *desc, char *chid);
 
 // Image functions
 extern int calibrate(float **prow, int nrow, int offset, int width);
@@ -90,11 +94,11 @@ int main(int argc, char **argv) {
 		usage();
 	}
 
-	options_t opts = { "r", "", 19, "", "." };
+	options_t opts = { "r", "", 19, "", ".", 0 };
 
 	// Parse arguments
 	int opt;
-    while ((opt = getopt(argc, argv, "c:m:d:i:s:e:")) != EOF) {
+    while ((opt = getopt(argc, argv, "c:m:d:i:s:e:r")) != EOF) {
 		switch (opt) {
 			case 'd':
 				opts.path = optarg;
@@ -118,6 +122,9 @@ int main(int argc, char **argv) {
 			case 'e':
 				opts.effects = optarg;
 				break;
+			case 'r':
+				opts.realtime = 1;
+				break;
 			default:
 				usage();
 		}
@@ -153,6 +160,8 @@ static int processAudio(char *filename, options_t *opts){
 	strcpy(path, dirname(path));
 	sscanf(basename(filename), "%[^.].%s", img.name, extension);
 
+	if(opts->realtime) initWriter(opts, &img, IMG_WIDTH, MAX_HEIGHT, "Unprocessed realtime image", "r");
+
 	if(strcmp(extension, "png") == 0){
 		// Read PNG into image buffer
 		printf("Reading %s", filename);
@@ -166,7 +175,7 @@ static int processAudio(char *filename, options_t *opts){
 			exit(EPERM);
 
 		// Build image
-		for (img.nrow = 0; img.nrow < 3000; img.nrow++) {
+		for (img.nrow = 0; img.nrow < MAX_HEIGHT; img.nrow++) {
 			// Allocate memory for this row
 			img.prow[img.nrow] = (float *) malloc(sizeof(float) * 2150);
 				
@@ -174,6 +183,8 @@ static int processAudio(char *filename, options_t *opts){
 			if (getpixelrow(img.prow[img.nrow], img.nrow, &zenith) == 0)
 				break;
 
+			if(opts->realtime) pushRow(img.prow[img.nrow], IMG_WIDTH);
+
 			fprintf(stderr, "Row: %d\r", img.nrow);
 			fflush(stderr);
 		}
@@ -182,6 +193,8 @@ static int processAudio(char *filename, options_t *opts){
 		sf_close(audioFile);
 	}
 
+	if(opts->realtime) closeWriter();
+
 	printf("\nTotal rows: %d\n", img.nrow);
 
 	// Fallback for detecting the zenith
@@ -248,7 +261,7 @@ static int processAudio(char *filename, options_t *opts){
 	// Channel B
 	if (CONTAINS(opts->type, 'b')) {
 		sprintf(desc, "%s (%s)", ch.id[img.chB], ch.name[img.chB]);
-		ImageOut(opts, &img, CHB_OFFSET, CH_WIDTH, desc, ch.id[img.chA], NULL);
+		ImageOut(opts, &img, CHB_OFFSET, CH_WIDTH, desc, ch.id[img.chB], NULL);
 	}
 
 	// Distribution image
@@ -292,7 +305,7 @@ static int initsnd(char *filename) {
 
 // Read samples from the wave file
 int getsample(float *sample, int nb) {
-    return sf_read_float(audioFile, sample, nb);
+	return sf_read_float(audioFile, sample, nb);
 }
 
 static void usage(void) {
@@ -304,7 +317,7 @@ static void usage(void) {
 	"     h: Histogram equalise\n"
 	"     d: Denoise\n"
 	"     p: Precipitation\n"
-	" -i [r|a|b|c|t] Output image\n"
+	" -i [r|a|b|c|t|m] Output image\n"
 	"     r: Raw\n"
 	"     a: Channel A\n"
 	"     b: Channel B\n"
@@ -314,7 +327,8 @@ static void usage(void) {
 	" -d <dir>       Image destination directory.\n"
 	" -s [15-19]     Satellite number\n"
 	" -c <file>      False color config file\n"
-	" -m <file>      Map file\n");
+	" -m <file>      Map file\n"
+	" -r             Realtime decode");
 
    	exit(EINVAL);
 }

offsets.h

@@ -27,6 +27,7 @@
 #define CHA_OFFSET (SYNC_WIDTH+SPC_WIDTH)
 #define CHB_OFFSET (SYNC_WIDTH+SPC_WIDTH+CH_WIDTH+TELE_WIDTH+SYNC_WIDTH+SPC_WIDTH)
 #define TOTAL_TELE (SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH+SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH)
+#define MAX_HEIGHT 3000
 
 #define CLIP(val, bottom, top) (val > top ? top : (val > bottom ? val : bottom))
 #define CONTAINS(str, char) (strchr(str, (int) char) != NULL)

pngio.c

@@ -30,10 +30,27 @@ typedef struct {
     float r, g, b;
 } rgb_t;
 
+typedef struct {
+	float *prow[MAX_HEIGHT]; // Row buffers
+	int nrow; // Number of rows
+	int chA, chB; // ID of each channel
+	char name[256]; // Stripped filename
+} image_t;
+
+typedef struct {
+	char *type; // Output image type
+	char *effects;
+	int   satnum; // The satellite number
+	char *map; // Path to a map file
+	char *path; // Output directory
+	int   realtime;
+} options_t;
+
 extern int zenith;
 extern char PrecipPalette[256*3];
 extern rgb_t applyPalette(char *palette, int val);
 extern rgb_t RGBcomposite(rgb_t top, float top_a, rgb_t bottom, float bottom_a);
+extern rgb_t falsecolor(float vis, float temp);
 
 int mapOverlay(char *filename, rgb_t **crow, int nrow, int zenith, int MCIR) {
 	FILE *fp = fopen(filename, "rb");
@@ -118,7 +135,6 @@ int mapOverlay(char *filename, rgb_t **crow, int nrow, int zenith, int MCIR) {
 
 			// Map overlay on channel A
 			crow[y][cha] = RGBcomposite(map, alpha, crow[y][cha], 1);
-
 			// Map overlay on channel B
 			if(!MCIR)
 				crow[y][chb] = RGBcomposite(map, alpha, crow[y][chb], 1);
@@ -188,33 +204,20 @@ int readRawImage(char *filename, float **prow, int *nrow) {
 	return 1;
 }
 
-typedef struct {
-	float *prow[3000]; // Row buffers
-	int nrow; // Number of rows
-	int chA, chB; // ID of each channel
-	char name[256]; // Stripped filename
-} image_t;
+png_text text_ptr[] = {
+	{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION},
+	{PNG_TEXT_COMPRESSION_NONE, "Channel", "Unknown", 7},
+	{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}
+};
 
-typedef struct {
-	char *type; // Output image type
-	char *effects;
-	int   satnum; // The satellite number
-	char *map; // Path to a map file
-	char *path; // Output directory
-} options_t;
-
-//int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, char *palette, char *effects, char *mapFile) {
 int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char *chid, char *palette){
 	char outName[384];
     sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid);
 	
-	FILE *pngfile;
+	text_ptr[1].text = desc;
+	text_ptr[1].text_length = sizeof(desc);
 
-	png_text text_ptr[] = {
-		{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION},
-		{PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)},
-		{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}
-	};
+	FILE *pngfile;
 
 	// Reduce the width of the image to componsate for the missing telemetry
 	int fcimage = strcmp(desc, "False Color") == 0;
@@ -229,13 +232,13 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
     if (!png_ptr) {
 		png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
 		fprintf(stderr, ERR_PNG_WRITE);
-		return(0);
+		return 0;
     }
     png_infop info_ptr = png_create_info_struct(png_ptr);
     if (!info_ptr) {
 		png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
 		fprintf(stderr, ERR_PNG_INFO);
-		return(0);
+		return 0;
     }
 
 	int greyscale = 0;
@@ -254,35 +257,34 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
 	}
 
     png_set_text(png_ptr, info_ptr, text_ptr, 3);
-
-	// Channel = 25cm wide
     png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER);
 
 	// Init I/O
     pngfile = fopen(outName, "wb");
     if (!pngfile) {
 		fprintf(stderr, ERR_FILE_WRITE, outName);
-		return(1);
+		return 1;
     }
     png_init_io(png_ptr, pngfile);
     png_write_info(png_ptr, info_ptr);
 
 	// Move prow into crow, crow ~ color rows
 	rgb_t *crow[img->nrow];
-	for(int y = 0; y < img->nrow; y++){
-		crow[y] = (rgb_t *) malloc(sizeof(rgb_t) * IMG_WIDTH);
+	if(!greyscale){
+		for(int y = 0; y < img->nrow; y++){
+			crow[y] = (rgb_t *) malloc(sizeof(rgb_t) * IMG_WIDTH);
 
-		for(int x = 0; x < IMG_WIDTH; x++){
-			if(palette == NULL){				
-				crow[y][x].r = crow[y][x].g = crow[y][x].b = CLIP(img->prow[y][x], 0, 255);
-			}else{
-				crow[y][x] = applyPalette(palette, img->prow[y][x]);
+			for(int x = 0; x < IMG_WIDTH; x++){
+				if(palette == NULL)	
+					crow[y][x].r = crow[y][x].g = crow[y][x].b = CLIP(img->prow[y][x], 0, 255);
+				else
+					crow[y][x] = applyPalette(palette, img->prow[y][x]);
 			}
 		}
 	}
 
 	// Precipitation
-	// TODO: use temperature calibration for accuracy
+	// TODO: use temperature calibration
 	if(CONTAINS(opts->effects, 'p')){
 		for(int y = 0; y < img->nrow; y++){
 			for(int x = 0; x < CH_WIDTH; x++){
@@ -292,6 +294,7 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
 		}
 	}
 
+	// Map stuff
 	if(opts->map != NULL && opts->map[0] != '\0'){
 		if(mapOverlay(opts->map, crow, img->nrow, zenith, strcmp(chid, "MCIR") == 0) == 0){
 			fprintf(stderr, "Skipping MCIR generation; see above.\n");
@@ -304,11 +307,10 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
 
 	printf("Writing %s", outName);
 
-	extern rgb_t falsecolor(float vis, float temp);
-
 	// Build image
     for (int y = 0; y < img->nrow; y++) {
 		png_color pix[width];
+		png_byte gpix[width];
 		
 		int skip = 0;
 		for (int x = 0; x < width; x++) {
@@ -320,31 +322,24 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
 					case CH_WIDTH:
 						skip += TELE_WIDTH + SYNC_WIDTH + SPC_WIDTH;
 						break;
-					case CH_WIDTH*2:
-						skip += TELE_WIDTH;
-						break;
 				}
 			}
+
 			if(greyscale){
-				// Horrific but works
-				if(x % 3 == 0){
-					pix[x/3].red   = img->prow[y][x + skip + offset  ];
-					pix[x/3].green = img->prow[y][x + skip + offset+1];
-					pix[x/3].blue  = img->prow[y][x + skip + offset+2];
-				}
+				gpix[x] = img->prow[y][x + skip + offset];
 			}else if(fcimage){
 				rgb_t pixel  = falsecolor(img->prow[y][x + CHA_OFFSET], img->prow[y][x + CHB_OFFSET]);
-				pix[x].red   = pixel.r;
-				pix[x].green = pixel.g;
-				pix[x].blue  = pixel.b;
+				pix[x] = (png_color){pixel.r, pixel.g, pixel.b};
 			}else{
-				pix[x].red   = crow[y][x + skip + offset].r;
-				pix[x].green = crow[y][x + skip + offset].g;
-				pix[x].blue  = crow[y][x + skip + offset].b;
+				pix[x] = (png_color){crow[y][x + skip + offset].r, crow[y][x + skip + offset].g, crow[y][x + skip + offset].b};
 			}
 		}
 		
-		png_write_row(png_ptr, (png_bytep) pix);
+		if(greyscale){
+			png_write_row(png_ptr, (png_bytep) gpix);
+		}else{
+			png_write_row(png_ptr, (png_bytep) pix);
+		}
     }
 
 	// Tidy up
@@ -353,6 +348,67 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c
     printf("\nDone\n");
     png_destroy_write_struct(&png_ptr, &info_ptr);
 
-    return(1);
+    return 1;
+}
+
+// TODO: remove these from the global scope
+png_structp rt_png_ptr;
+png_infop rt_info_ptr;
+FILE *rt_pngfile;
+
+int initWriter(options_t *opts, image_t *img, int width, int height, char *desc, char *chid){
+	char outName[384];
+    sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid);
+
+	text_ptr[1].text = desc;
+	text_ptr[1].text_length = sizeof(desc);
+
+	// Create writer
+    rt_png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+    if (!rt_png_ptr) {
+		png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL);
+		fprintf(stderr, ERR_PNG_WRITE);
+		return 0;
+    }
+    rt_info_ptr = png_create_info_struct(rt_png_ptr);
+    if (!rt_info_ptr) {
+		png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL);
+		fprintf(stderr, ERR_PNG_INFO);
+		return 0;
+    }
+
+	// Greyscale image
+    png_set_IHDR(rt_png_ptr, rt_info_ptr, width, height,
+				 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE,
+			 	 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+	
+	png_set_text(rt_png_ptr, rt_info_ptr, text_ptr, 3);
+
+	// Channel = 25cm wide
+    png_set_pHYs(rt_png_ptr, rt_info_ptr, 3636, 3636, PNG_RESOLUTION_METER);
+
+	// Init I/O
+    rt_pngfile = fopen(outName, "wb");
+    if (!rt_pngfile) {
+		fprintf(stderr, ERR_FILE_WRITE, outName);
+		return 0;
+    }
+    png_init_io(rt_png_ptr, rt_pngfile);
+    png_write_info(rt_png_ptr, rt_info_ptr);
+	
+	return 1;
+}
+
+void pushRow(float *row, int width){
+	png_byte pix[width];
+	for(int i = 0; i < width; i++)
+		pix[i] = row[i];
+	
+	png_write_row(rt_png_ptr, (png_bytep) pix);
 }
 
+void closeWriter(){
+	png_write_end(rt_png_ptr, rt_info_ptr);
+    fclose(rt_pngfile);
+    png_destroy_write_struct(&rt_png_ptr, &rt_info_ptr);
+}